[PATCH v1 6/8] ipa: libipa: Add lsc polynomial class

[Stefan Klug]

Add a basic class to represent polynomials as specified in the DNG spec for
vignetting correction.

Signed-off-by: Stefan Klug <stefan.klug at ideasonboard.com>
---
 src/ipa/libipa/meson.build    |   2 +
 src/ipa/libipa/polynomial.cpp |  23 ++++++++
 src/ipa/libipa/polynomial.h   | 107 ++++++++++++++++++++++++++++++++++
 3 files changed, 132 insertions(+)
 create mode 100644 src/ipa/libipa/polynomial.cpp
 create mode 100644 src/ipa/libipa/polynomial.h

diff --git a/src/ipa/libipa/meson.build b/src/ipa/libipa/meson.build
index 3740178b2909..8b085750807a 100644
--- a/src/ipa/libipa/meson.build
+++ b/src/ipa/libipa/meson.build
@@ -10,6 +10,7 @@ libipa_headers = files([
     'interpolator.h',
     'matrix.h',
     'module.h',
+    'polynomial.h',
     'pwl.h',
     'vector.h',
 ])
@@ -24,6 +25,7 @@ libipa_sources = files([
     'interpolator.cpp',
     'matrix.cpp',
     'module.cpp',
+    'polynomial.cpp',
     'pwl.cpp',
     'vector.cpp',
 ])
diff --git a/src/ipa/libipa/polynomial.cpp b/src/ipa/libipa/polynomial.cpp
new file mode 100644
index 000000000000..39b17b29314b
--- /dev/null
+++ b/src/ipa/libipa/polynomial.cpp
@@ -0,0 +1,23 @@
+/* SPDX-License-Identifier: LGPL-2.1-or-later */
+/*
+ * Copyright (C) 2024, Ideas On Board
+ *
+ * Polynomal class to represent lens shading correction
+ */
+
+#include "polynomial.h"
+
+#include <libcamera/base/log.h>
+
+/**
+ * \file polynomial.h
+ * \brief Polynomial class
+ */
+
+namespace libcamera {
+
+LOG_DEFINE_CATEGORY(Polynomial)
+
+namespace ipa {
+}
+} // namespace libcamera
diff --git a/src/ipa/libipa/polynomial.h b/src/ipa/libipa/polynomial.h
new file mode 100644
index 000000000000..971be062297f
--- /dev/null
+++ b/src/ipa/libipa/polynomial.h
@@ -0,0 +1,107 @@
+/* SPDX-License-Identifier: LGPL-2.1-or-later */
+/*
+ * Copyright (C) 2024, Ideas On Board
+ *
+ * Vector and related operations
+ */
+#pragma once
+
+#include <algorithm>
+#include <array>
+#include <assert.h>
+#include <cmath>
+
+#include <libcamera/base/log.h>
+#include <libcamera/base/span.h>
+
+#include "libcamera/internal/yaml_parser.h"
+
+namespace libcamera {
+
+LOG_DECLARE_CATEGORY(Polynomial)
+
+namespace ipa {
+
+constexpr int kNumCoefficients = 5;
+class Polynomial
+{
+public:
+	/**
+	 * Note: cx and cy are relative to the image size
+	 */
+	Polynomial(double cx = 0.0, double cy = 0.0, double k0 = 0.0,
+		      double k1 = 0.0, double k2 = 0.0, double k3 = 0.0,
+		      double k4 = 0.0)
+		: cx_(cx), cy_(cy),
+		  coefficients_({ k0, k1, k2, k3, k4 })
+	{
+	}
+
+	double sampleAtNormalizedPixelPos(double x, double y) const
+	{
+		double dx = x - cnx_;
+		double dy = y - cny_;
+		double r = sqrt(dx * dx + dy * dy);
+		double res = 1.0;
+		for (int i = 0; i < kNumCoefficients; i++) {
+			res += coefficients_[i] * std::pow(r, (i + 1) * 2);
+		}
+		return res;
+	}
+
+	/**
+         * Returns m according to dng spec. 𝑚 represents the Euclidean distance
+	 * (in pixels) from the optical center to the farthest pixel in the
+	 * image.
+	 */
+	double getM() const
+	{
+		double cpx = imageSize_.width * cx_;
+		double cpy = imageSize_.height * cy_;
+		double mx = std::max(cpx, std::fabs(imageSize_.width - cpx));
+		double my = std::max(cpy, std::fabs(imageSize_.height - cpy));
+
+		return sqrt(mx * mx + my * my);
+	}
+
+	void setReferenceImageSize(const Size &size)
+	{
+		imageSize_ = size;
+		/* Calculate normalized centers */
+		double m = getM();
+		cnx_ = (size.width * cx_) / m;
+		cny_ = (size.height * cy_) / m;
+	}
+
+private:
+	double cx_;
+	double cy_;
+	double cnx_;
+	double cny_;
+	std::array<double, kNumCoefficients> coefficients_;
+
+	Size imageSize_;
+};
+
+} /* namespace ipa */
+
+template<>
+struct YamlObject::Getter<ipa::Polynomial> {
+	std::optional<ipa::Polynomial> get(const YamlObject &obj) const
+	{
+		std::optional<double> cx = obj["cx"].get<double>();
+		std::optional<double> cy = obj["cy"].get<double>();
+		std::optional<double> k0 = obj["k0"].get<double>();
+		std::optional<double> k1 = obj["k1"].get<double>();
+		std::optional<double> k2 = obj["k2"].get<double>();
+		std::optional<double> k3 = obj["k3"].get<double>();
+		std::optional<double> k4 = obj["k4"].get<double>();
+
+		if (!(cx && cy && k0 && k1 && k2 && k3 && k4))
+			LOG(Polynomial, Error) << "Polynomial is missing parameter ";
+
+		return ipa::Polynomial(*cx, *cy, *k0, *k1, *k2, *k3, *k4);
+	}
+};
+
+} /* namespace libcamera */
-- 
2.43.0