[{"id":33082,"web_url":"https://patchwork.libcamera.org/comment/33082/","msgid":"","date":"2025-01-14T22:40:24","subject":"Re: [PATCH v1 09/11] libipa: Add bayesian AWB algorithm","submitter":{"id":17,"url":"https://patchwork.libcamera.org/api/people/17/","name":"Paul Elder","email":"paul.elder@ideasonboard.com"},"content":"On Thu, Jan 09, 2025 at 12:54:00PM +0100, Stefan Klug wrote:\n> The bayesian AWB algorithm is an AWB algorithm that takes prior\n> probabilities for a given light source dependent on the current lux\n> level into account.\n> \n> The biggest improvement compared to the grey world model comes from the\n> search of the ideal white point on the CT curve. The algorithm walks the\n> CT curve to minimize the colour error for a given statistics. After the\n> minimium is found it additionally tries to search the area around that\n> spot and also off the curve. So even without defined prior probabilities\n> this algorithm provides much better results than the grey world\n> algorithm.\n> \n> The logic for this code was taken from the RaspberryPi implementation.\n> The logic was only minimally adjusted for usage with the rkisp1 and a\n> few things were left out (see doxygen doc for the AwbBayes class). The\n> code is refactored to better fit the libcamera code style and to make\n> use of the syntactic sugar provided by the Interpolator and Vector\n> classes.\n> \n> Signed-off-by: Stefan Klug \n> ---\n> src/ipa/libipa/awb_bayes.cpp | 457 +++++++++++++++++++++++++++++++++++\n> src/ipa/libipa/awb_bayes.h | 67 +++++\n> src/ipa/libipa/meson.build | 2 +\n> 3 files changed, 526 insertions(+)\n> create mode 100644 src/ipa/libipa/awb_bayes.cpp\n> create mode 100644 src/ipa/libipa/awb_bayes.h\n> \n> diff --git a/src/ipa/libipa/awb_bayes.cpp b/src/ipa/libipa/awb_bayes.cpp\n> new file mode 100644\n> index 000000000000..1e69ecd3e3f3\n> --- /dev/null\n> +++ b/src/ipa/libipa/awb_bayes.cpp\n> @@ -0,0 +1,457 @@\n> +/* SPDX-License-Identifier: LGPL-2.1-or-later */\n> +/*\n> + * Copyright (C) 2019, Raspberry Pi Ltd\n> + * Copyright (C) 2024 Ideas on Board Oy\n> + *\n> + * Implementation of a bayesian AWB algorithm\n> + */\n> +\n> +#include \"awb_bayes.h\"\n> +\n> +#include \n> +\n> +#include \n> +#include \n> +\n> +#include \"colours.h\"\n> +\n> +/**\n> + * \\file awb_bayes.h\n> + * \\brief Implementation of bayesian auto white balance algorithm\n> + *\n> + * This implementation is based on the initial implementation done by\n> + * RaspberryPi.\n> + * \\todo: Documentation\n> + *\n> + * \\todo As the statistics module of the rkisp1 provides less data than the one\n> + * from the RaspberryPi (vc4). The RaspberryPi statistics measure a grid of\n> + * zones while the rkisp1 ony measures a single area. Therefore this algorithm\n> + * doesn't contain all the features implemented by RaspberryPi.\n> + * The following parts are not implemented:\n> + *\n> + * - min_pixels: minimum proportion of pixels counted within AWB region for it\n> + * to be \"useful\"\n> + * - min_g: minimum G value of those pixels, to be regarded a \"useful\"\n> + * - min_regions: number of AWB regions that must be \"useful\" in order to do the\n> + * AWB calculation\n> + * - deltaLimit: clamp on colour error term (so as not to penalize non-grey\n> + * excessively)\n> + * - bias_proportion: The biasProportion parameter adds a small proportion of\n> + * the counted pixels to a region biased to the biasCT colour temperature.\n> + * A typical value for biasProportion would be between 0.05 to 0.1.\n> + * - bias_ct: CT target for the search bias\n> + * - sensitivityR: red sensitivity ratio (set to canonical sensor's R/G divided\n> + * by this sensor's R/G)\n> + * - sensitivityB: blue sensitivity ratio (set to canonical sensor's B/G divided\n> + * by this sensor's B/G)\n> + */\n> +\n> +namespace libcamera {\n> +\n> +LOG_DECLARE_CATEGORY(Awb)\n> +\n> +namespace ipa {\n> +\n> +/**\n> + * \\brief Step size control for CT search\n> + */\n> +constexpr double kSearchStep = 0.2;\n> +\n> +/**\n> + * \\copydoc Interpolator::interpolate()\n> + */\n> +template<>\n> +void Interpolator::interpolate(const Pwl &a, const Pwl &b, Pwl &dest, double lambda)\n> +{\n> +\tdest = Pwl::combine(a, b,\n> +\t\t\t [=](double /*x*/, double y0, double y1) -> double {\n> +\t\t\t\t return y0 * (1.0 - lambda) + y1 * lambda;\n> +\t\t\t });\n> +}\n> +\n> +/**\n> + * \\class AwbBayes\n> + * \\brief Implementation of a bayesian auto white balance algorithm\n> + *\n> + * In a bayesian AWB algorithm the auto white balance estimation is improved by\n> + * taking the likelihood of a given lightsource based on the estimated lux level\n> + * into account. E.g. If it is very bright we can assume that we are outside and\n> + * that colour temperatures around 6500 are preferred.\n> + *\n> + * The second part of this algorithm is the search for the most likely colour\n> + * temperature. It is implemented in AwbBayes::coarseSearch() and in\n> + * AwbBayes::fineSearch(). The search works very well without prior likelihoods\n> + * and therefore the algorithm itself provides very good results even without\n> + * prior likelihoods.\n> + */\n> +\n> +/**\n> + * \\var AwbBayes::transversePos_\n> + * \\brief How far to wander off CT curve towards \"more purple\"\n> + */\n> +\n> +/**\n> + * \\var AwbBayes::transverseNeg_\n> + * \\brief How far to wander off CT curve towards \"more green\"\n> + */\n> +\n> +/**\n> + * \\var AwbBayes::currentMode_\n> + * \\brief The currently selected mode\n> + */\n> +\n> +int AwbBayes::init(const YamlObject &tuningData)\n> +{\n> +\tint ret = colourGainCurve_.readYaml(tuningData[\"colourGains\"], \"ct\", \"gains\");\n> +\tif (ret) {\n> +\t\tLOG(Awb, Error)\n> +\t\t\t<< \"Failed to parse 'colourGains' \"\n> +\t\t\t<< \"parameter from tuning file\";\n> +\t\treturn ret;\n> +\t}\n> +\n> +\tctR_.clear();\n> +\tctB_.clear();\n> +\tfor (const auto &[ct, g] : colourGainCurve_.data()) {\n> +\t\tctR_.append(ct, 1.0 / g[0]);\n> +\t\tctB_.append(ct, 1.0 / g[1]);\n> +\t}\n> +\n> +\t/* We will want the inverse functions of these too. */\n> +\tctRInverse_ = ctR_.inverse().first;\n> +\tctBInverse_ = ctB_.inverse().first;\n> +\n> +\tret = readPriors(tuningData);\n> +\tif (ret) {\n> +\t\tLOG(Awb, Error) << \"Failed to read priors\";\n> +\t\treturn ret;\n> +\t}\n> +\n> +\tret = parseModeConfigs(tuningData);\n> +\tif (ret) {\n> +\t\tLOG(Awb, Error)\n> +\t\t\t<< \"Failed to parse mode parameter from tuning file\";\n> +\t\treturn ret;\n> +\t}\n> +\tcurrentMode_ = &modes_[controls::AwbAuto];\n> +\n> +\ttransversePos_ = tuningData[\"transversePos\"].get(0.01);\n> +\ttransverseNeg_ = tuningData[\"transverseNeg\"].get(0.01);\n> +\tif (transversePos_ <= 0 || transverseNeg_ <= 0) {\n> +\t\tLOG(Awb, Error) << \"AwbConfig: transversePos/Neg must be > 0\";\n> +\t\treturn -EINVAL;\n> +\t}\n> +\n> +\treturn 0;\n> +}\n> +\n> +int AwbBayes::readPriors(const YamlObject &tuningData)\n> +{\n> +\tconst auto &priorsList = tuningData[\"priors\"];\n> +\tstd::map priors;\n> +\n> +\tif (!priorsList) {\n> +\t\tLOG(Awb, Error) << \"No priors specified\";\n> +\t\treturn -EINVAL;\n> +\t}\n> +\n> +\tfor (const auto &p : priorsList.asList()) {\n> +\t\tif (!p.contains(\"lux\")) {\n> +\t\t\tLOG(Awb, Error) << \"Missing lux value\";\n> +\t\t\treturn -EINVAL;\n> +\t\t}\n> +\n> +\t\tuint32_t lux = p[\"lux\"].get(0);\n> +\t\tif (priors.count(lux)) {\n> +\t\t\tLOG(Awb, Error) << \"Duplicate prior for lux value \" << lux;\n> +\t\t\treturn -EINVAL;\n> +\t\t}\n> +\n> +\t\tstd::vector temperatures =\n> +\t\t\tp[\"ct\"].getList().value_or(std::vector{});\n> +\t\tstd::vector probabilites =\n> +\t\t\tp[\"probability\"].getList().value_or(std::vector{});\n> +\n> +\t\tif (temperatures.size() != probabilites.size()) {\n> +\t\t\tLOG(Awb, Error)\n> +\t\t\t\t<< \"Ct and probability array sizes are unequal\";\n> +\t\t\treturn -EINVAL;\n> +\t\t}\n> +\n> +\t\tif (temperatures.empty()) {\n> +\t\t\tLOG(Awb, Error)\n> +\t\t\t\t<< \"Ct and probability arrays are empty\";\n> +\t\t\treturn -EINVAL;\n> +\t\t}\n> +\n> +\t\tstd::map ctToProbability;\n> +\t\tfor (unsigned int i = 0; i < temperatures.size(); i++) {\n> +\t\t\tint t = temperatures[i];\n> +\t\t\tif (ctToProbability.count(t)) {\n> +\t\t\t\tLOG(Awb, Error) << \"Duplicate ct value\";\n> +\t\t\t\treturn -EINVAL;\n> +\t\t\t}\n> +\n> +\t\t\tctToProbability[t] = probabilites[i];\n> +\t\t}\n> +\n> +\t\tauto &pwl = priors[lux];\n> +\t\tfor (const auto &[ct, prob] : ctToProbability) {\n> +\t\t\tpwl.append(ct, prob);\n> +\t\t}\n> +\t}\n> +\n> +\tif (priors.empty()) {\n> +\t\tLOG(Awb, Error) << \"No priors specified\";\n> +\t\t;\n\nI think we don't need this line.\n\n\nOther than that, looks good to me.\n\nReviewed-by: Paul Elder \n\n> +\t\treturn -EINVAL;\n> +\t}\n> +\n> +\tpriors_.setData(std::move(priors));\n> +\n> +\treturn 0;\n> +}\n> +\n> +void AwbBayes::handleControls(const ControlList &controls)\n> +{\n> +\tauto mode = controls.get(controls::AwbMode);\n> +\tif (mode) {\n> +\t\tauto it = modes_.find(static_cast(*mode));\n> +\t\tif (it != modes_.end())\n> +\t\t\tcurrentMode_ = &it->second;\n> +\t\telse\n> +\t\t\tLOG(Awb, Error) << \"Unknown AWB mode\";\n> +\t}\n> +}\n> +\n> +RGB AwbBayes::gainsFromColourTemperature(double colourTemperature)\n> +{\n> +\t/*\n> +\t * \\todo: In the RaspberryPi code, the ct curve was interpolated in\n> +\t * the white point space (1/x) not in gains space. This feels counter\n> +\t * intuitive, as the gains are in linear space. But I can't prove it.\n> +\t */\n> +\tconst auto &gains = colourGainCurve_.getInterpolated(colourTemperature);\n> +\treturn { { gains[0], 1.0, gains[1] } };\n> +}\n> +\n> +AwbResult AwbBayes::calculateAwb(const AwbStats &stats, int lux)\n> +{\n> +\tipa::Pwl prior;\n> +\tif (lux > 0) {\n> +\t\tprior = priors_.getInterpolated(lux);\n> +\t\tprior.map([](double x, double y) {\n> +\t\t\tLOG(Awb, Debug) << \"(\" << x << \",\" << y << \")\";\n> +\t\t});\n> +\t} else {\n> +\t\tprior.append(0, 1.0);\n> +\t}\n> +\n> +\tdouble t = coarseSearch(prior, stats);\n> +\tdouble r = ctR_.eval(t);\n> +\tdouble b = ctB_.eval(t);\n> +\tLOG(Awb, Debug)\n> +\t\t<< \"After coarse search: r \" << r << \" b \" << b << \" (gains r \"\n> +\t\t<< 1 / r << \" b \" << 1 / b << \")\";\n> +\n> +\t/*\n> +\t * Original comment from RaspberryPi:\n> +\t * Not entirely sure how to handle the fine search yet. Mostly the\n> +\t * estimated CT is already good enough, but the fine search allows us to\n> +\t * wander transversely off the CT curve. Under some illuminants, where\n> +\t * there may be more or less green light, this may prove beneficial,\n> +\t * though I probably need more real datasets before deciding exactly how\n> +\t * this should be controlled and tuned.\n> +\t */\n> +\tfineSearch(t, r, b, prior, stats);\n> +\tLOG(Awb, Debug)\n> +\t\t<< \"After fine search: r \" << r << \" b \" << b << \" (gains r \"\n> +\t\t<< 1 / r << \" b \" << 1 / b << \")\";\n> +\n> +\treturn { { { 1.0 / r, 1.0, 1.0 / b } }, t };\n> +}\n> +\n> +double AwbBayes::coarseSearch(const ipa::Pwl &prior, const AwbStats &stats) const\n> +{\n> +\tstd::vector points;\n> +\tsize_t bestPoint = 0;\n> +\tdouble t = currentMode_->ctLo;\n> +\tint spanR = -1;\n> +\tint spanB = -1;\n> +\n> +\t/* Step down the CT curve evaluating log likelihood. */\n> +\twhile (true) {\n> +\t\tdouble r = ctR_.eval(t, &spanR);\n> +\t\tdouble b = ctB_.eval(t, &spanB);\n> +\t\tRGB gains({ 1 / r, 1.0, 1 / b });\n> +\t\tdouble delta2Sum = stats.computeColourError(gains);\n> +\t\tdouble priorLogLikelihood = prior.eval(prior.domain().clamp(t));\n> +\t\tdouble finalLogLikelihood = delta2Sum - priorLogLikelihood;\n> +\n> +\t\tLOG(Awb, Debug) << \"Coarse search t: \" << t\n> +\t\t\t\t<< \" gains: \" << gains\n> +\t\t\t\t<< \" error: \" << delta2Sum\n> +\t\t\t\t<< \" prior: \" << priorLogLikelihood\n> +\t\t\t\t<< \" likelihood: \" << finalLogLikelihood;\n> +\n> +\t\tpoints.push_back({ { t, finalLogLikelihood } });\n> +\t\tif (points.back().y() < points[bestPoint].y())\n> +\t\t\tbestPoint = points.size() - 1;\n> +\n> +\t\tif (t == currentMode_->ctHi)\n> +\t\t\tbreak;\n> +\n> +\t\t/*\n> +\t\t * Ensure even steps along the r/b curve by scaling them by the\n> +\t\t * current t.\n> +\t\t */\n> +\t\tt = std::min(t + t / 10 * kSearchStep, currentMode_->ctHi);\n> +\t}\n> +\n> +\tt = points[bestPoint].x();\n> +\tLOG(Awb, Debug) << \"Coarse search found CT \" << t;\n> +\n> +\t/*\n> +\t * We have the best point of the search, but refine it with a quadratic\n> +\t * interpolation around its neighbors.\n> +\t */\n> +\tif (points.size() > 2) {\n> +\t\tbestPoint = std::clamp(bestPoint, 1ul, points.size() - 2);\n> +\t\tt = interpolateQuadratic(points[bestPoint - 1],\n> +\t\t\t\t\t points[bestPoint],\n> +\t\t\t\t\t points[bestPoint + 1]);\n> +\t\tLOG(Awb, Debug)\n> +\t\t\t<< \"After quadratic refinement, coarse search has CT \"\n> +\t\t\t<< t;\n> +\t}\n> +\n> +\treturn t;\n> +}\n> +\n> +void AwbBayes::fineSearch(double &t, double &r, double &b, ipa::Pwl const &prior, const AwbStats &stats) const\n> +{\n> +\tint spanR = -1;\n> +\tint spanB = -1;\n> +\tdouble step = t / 10 * kSearchStep * 0.1;\n> +\tint nsteps = 5;\n> +\tctR_.eval(t, &spanR);\n> +\tctB_.eval(t, &spanB);\n> +\tdouble rDiff = ctR_.eval(t + nsteps * step, &spanR) -\n> +\t\t ctR_.eval(t - nsteps * step, &spanR);\n> +\tdouble bDiff = ctB_.eval(t + nsteps * step, &spanB) -\n> +\t\t ctB_.eval(t - nsteps * step, &spanB);\n> +\tPwl::Point transverse({ bDiff, -rDiff });\n> +\tif (transverse.length2() < 1e-6)\n> +\t\treturn;\n> +\t/*\n> +\t * transverse is a unit vector orthogonal to the b vs. r function\n> +\t * (pointing outwards with r and b increasing)\n> +\t */\n> +\ttransverse = transverse / transverse.length();\n> +\tdouble bestLogLikelihood = 0;\n> +\tdouble bestT = 0;\n> +\tPwl::Point bestRB;\n> +\tdouble transverseRange = transverseNeg_ + transversePos_;\n> +\tconst int maxNumDeltas = 12;\n> +\n> +\t/* a transverse step approximately every 0.01 r/b units */\n> +\tint numDeltas = floor(transverseRange * 100 + 0.5) + 1;\n> +\tnumDeltas = std::clamp(numDeltas, 3, maxNumDeltas);\n> +\n> +\t/*\n> +\t * Step down CT curve. March a bit further if the transverse range is\n> +\t * large.\n> +\t */\n> +\tnsteps += numDeltas;\n> +\tfor (int i = -nsteps; i <= nsteps; i++) {\n> +\t\tdouble tTest = t + i * step;\n> +\t\tdouble priorLogLikelihood =\n> +\t\t\tprior.eval(prior.domain().clamp(tTest));\n> +\t\tPwl::Point rbStart{ { ctR_.eval(tTest, &spanR),\n> +\t\t\t\t ctB_.eval(tTest, &spanB) } };\n> +\t\tPwl::Point samples[maxNumDeltas];\n> +\t\tint bestPoint = 0;\n> +\n> +\t\t/*\n> +\t\t * Sample numDeltas points transversely *off* the CT curve\n> +\t\t * in the range [-transverseNeg, transversePos].\n> +\t\t * The x values of a sample contains the distance and the y\n> +\t\t * value contains the corresponding log likelihood.\n> +\t\t */\n> +\t\tdouble transverseStep = transverseRange / (numDeltas - 1);\n> +\t\tfor (int j = 0; j < numDeltas; j++) {\n> +\t\t\tauto &p = samples[j];\n> +\t\t\tp.x() = -transverseNeg_ + transverseStep * j;\n> +\t\t\tPwl::Point rbTest = rbStart + transverse * p.x();\n> +\t\t\tRGB gains({ 1 / rbTest[0], 1.0, 1 / rbTest[1] });\n> +\t\t\tdouble delta2Sum = stats.computeColourError(gains);\n> +\t\t\tp.y() = delta2Sum - priorLogLikelihood;\n> +\n> +\t\t\tif (p.y() < samples[bestPoint].y())\n> +\t\t\t\tbestPoint = j;\n> +\t\t}\n> +\n> +\t\t/*\n> +\t\t * We have all samples transversely across the CT curve,\n> +\t\t * now let's do a quadratic interpolation for the best result.\n> +\t\t */\n> +\t\tbestPoint = std::clamp(bestPoint, 1, numDeltas - 2);\n> +\t\tdouble bestOffset = interpolateQuadratic(samples[bestPoint - 1],\n> +\t\t\t\t\t\t\t samples[bestPoint],\n> +\t\t\t\t\t\t\t samples[bestPoint + 1]);\n> +\t\tPwl::Point rbTest = rbStart + transverse * bestOffset;\n> +\t\tRGB gains({ 1 / rbTest[0], 1.0, 1 / rbTest[1] });\n> +\t\tdouble delta2Sum = stats.computeColourError(gains);\n> +\t\tdouble finalLogLikelihood = delta2Sum - priorLogLikelihood;\n> +\t\tLOG(Awb, Debug)\n> +\t\t\t<< \"Fine search t: \" << tTest\n> +\t\t\t<< \" r: \" << rbTest[0]\n> +\t\t\t<< \" b: \" << rbTest[1]\n> +\t\t\t<< \" offset: \" << bestOffset\n> +\t\t\t<< \" likelihood: \" << finalLogLikelihood\n> +\t\t\t<< (finalLogLikelihood < bestLogLikelihood ? \" NEW BEST\" : \"\");\n> +\n> +\t\tif (bestT == 0 || finalLogLikelihood < bestLogLikelihood) {\n> +\t\t\tbestLogLikelihood = finalLogLikelihood;\n> +\t\t\tbestT = tTest;\n> +\t\t\tbestRB = rbTest;\n> +\t\t}\n> +\t}\n> +\n> +\tt = bestT;\n> +\tr = bestRB[0];\n> +\tb = bestRB[1];\n> +\tLOG(Awb, Debug)\n> +\t\t<< \"Fine search found t \" << t << \" r \" << r << \" b \" << b;\n> +}\n> +\n> +/**\n> + * \\brief Find extremum of function\n> + * \\param[in] a First point\n> + * \\param[in] b Second point\n> + * \\param[in] c Third point\n> + *\n> + * Given 3 points on a curve, find the extremum of the function in that interval\n> + * by fitting a quadratic.\n> + *\n> + * \\return The x value of the extremum clamped to the interval [a.x(), c.x()]\n> + */\n> +double AwbBayes::interpolateQuadratic(Pwl::Point const &a, Pwl::Point const &b,\n> +\t\t\t\t Pwl::Point const &c) const\n> +{\n> +\tconst double eps = 1e-3;\n> +\tPwl::Point ca = c - a;\n> +\tPwl::Point ba = b - a;\n> +\tdouble denominator = 2 * (ba.y() * ca.x() - ca.y() * ba.x());\n> +\tif (std::abs(denominator) > eps) {\n> +\t\tdouble numerator = ba.y() * ca.x() * ca.x() - ca.y() * ba.x() * ba.x();\n> +\t\tdouble result = numerator / denominator + a.x();\n> +\t\treturn std::max(a.x(), std::min(c.x(), result));\n> +\t}\n> +\t/* has degenerated to straight line segment */\n> +\treturn a.y() < c.y() - eps ? a.x() : (c.y() < a.y() - eps ? c.x() : b.x());\n> +}\n> +\n> +} /* namespace ipa */\n> +\n> +} /* namespace libcamera */\n> diff --git a/src/ipa/libipa/awb_bayes.h b/src/ipa/libipa/awb_bayes.h\n> new file mode 100644\n> index 000000000000..5284f0ca4cc0\n> --- /dev/null\n> +++ b/src/ipa/libipa/awb_bayes.h\n> @@ -0,0 +1,67 @@\n> +/* SPDX-License-Identifier: LGPL-2.1-or-later */\n> +/*\n> + * Copyright (C) 2024 Ideas on Board Oy\n> + *\n> + * Base class for bayes AWB algorithms\n> + */\n> +\n> +#pragma once\n> +\n> +#include \n> +#include \n> +#include \n> +#include \n> +\n> +#include \n> +\n> +#include \n> +#include \n> +\n> +#include \"libcamera/internal/yaml_parser.h\"\n> +\n> +#include \"awb.h\"\n> +#include \"interpolator.h\"\n> +#include \"pwl.h\"\n> +#include \"vector.h\"\n> +\n> +namespace libcamera {\n> +\n> +namespace ipa {\n> +\n> +class AwbBayes : public AwbAlgorithm\n> +{\n> +public:\n> +\tAwbBayes() = default;\n> +\n> +\tint init(const YamlObject &tuningData) override;\n> +\tAwbResult calculateAwb(const AwbStats &stats, int lux) override;\n> +\tRGB gainsFromColourTemperature(double temperatureK) override;\n> +\tvoid handleControls(const ControlList &controls) override;\n> +\n> +private:\n> +\tint readPriors(const YamlObject &tuningData);\n> +\n> +\tvoid fineSearch(double &t, double &r, double &b, ipa::Pwl const &prior,\n> +\t\t\tconst AwbStats &stats) const;\n> +\tdouble coarseSearch(const ipa::Pwl &prior, const AwbStats &stats) const;\n> +\tdouble interpolateQuadratic(ipa::Pwl::Point const &a,\n> +\t\t\t\t ipa::Pwl::Point const &b,\n> +\t\t\t\t ipa::Pwl::Point const &c) const;\n> +\n> +\tInterpolator priors_;\n> +\tInterpolator> colourGainCurve_;\n> +\n> +\tipa::Pwl ctR_;\n> +\tipa::Pwl ctB_;\n> +\tipa::Pwl ctRInverse_;\n> +\tipa::Pwl ctBInverse_;\n> +\n> +\tdouble transversePos_;\n> +\tdouble transverseNeg_;\n> +\n> +\tModeConfig *currentMode_ = nullptr;\n> +};\n> +\n> +} /* namespace ipa */\n> +\n> +} /* namespace libcamera */\n> diff --git a/src/ipa/libipa/meson.build b/src/ipa/libipa/meson.build\n> index c550a6eb45b6..b519c8146d7c 100644\n> --- a/src/ipa/libipa/meson.build\n> +++ b/src/ipa/libipa/meson.build\n> @@ -3,6 +3,7 @@\n> libipa_headers = files([\n> 'agc_mean_luminance.h',\n> 'algorithm.h',\n> + 'awb_bayes.h',\n> 'awb_grey.h',\n> 'awb.h',\n> 'camera_sensor_helper.h',\n> @@ -22,6 +23,7 @@ libipa_headers = files([\n> libipa_sources = files([\n> 'agc_mean_luminance.cpp',\n> 'algorithm.cpp',\n> + 'awb_bayes.cpp',\n> 'awb_grey.cpp',\n> 'awb.cpp',\n> 'camera_sensor_helper.cpp',\n> -- \n> 2.43.0\n>","headers":{"Return-Path":"","X-Original-To":"parsemail@patchwork.libcamera.org","Delivered-To":"parsemail@patchwork.libcamera.org","Received":["from lancelot.ideasonboard.com (lancelot.ideasonboard.com\n\t[92.243.16.209])\n\tby patchwork.libcamera.org (Postfix) with ESMTPS id 95807C3301\n\tfor ;\n\tTue, 14 Jan 2025 22:41:03 +0000 (UTC)","from lancelot.ideasonboard.com (localhost [IPv6:::1])\n\tby lancelot.ideasonboard.com (Postfix) with ESMTP id 8CE9F68521;\n\tTue, 14 Jan 2025 23:41:02 +0100 (CET)","from perceval.ideasonboard.com (perceval.ideasonboard.com\n\t[213.167.242.64])\n\tby lancelot.ideasonboard.com (Postfix) with ESMTPS id BB72C607D6\n\tfor ;\n\tTue, 14 Jan 2025 23:41:00 +0100 (CET)","from pyrite.rasen.tech (unknown [173.16.167.215])\n\tby perceval.ideasonboard.com (Postfix) with ESMTPSA id 590AABEB;\n\tTue, 14 Jan 2025 23:40:02 +0100 (CET)"],"Authentication-Results":"lancelot.ideasonboard.com; dkim=pass (1024-bit key;\n\tunprotected) header.d=ideasonboard.com header.i=@ideasonboard.com\n\theader.b=\"m3fagTFM\"; dkim-atps=neutral","DKIM-Signature":"v=1; a=rsa-sha256; c=relaxed/simple; d=ideasonboard.com;\n\ts=mail; t=1736894403;\n\tbh=JX/RfvtIytwpecqXWKHheg6g/cCm50rMR4IHioXnZFc=;\n\th=Date:From:To:Cc:Subject:References:In-Reply-To:From;\n\tb=m3fagTFMSDRJzOrdDL0Y/d/qsgNBDbKBQ0h//+BV+OZNyom1EX1maQ0yV0YrTlibV\n\tJw2YQWf5Pgo3w3PezIqvlvaKETFhGFm/MGjJ3Xvj9CX9h761oQ0Qwekw4pGHOW59wq\n\tiN+93jofmU3ID6UD6LQ2X3xc+b0L9OwnHrNePRPU=","Date":"Tue, 14 Jan 2025 16:40:24 -0600","From":"Paul Elder ","To":"Stefan Klug ","Cc":"libcamera-devel@lists.libcamera.org","Subject":"Re: [PATCH v1 09/11] libipa: Add bayesian AWB algorithm","Message-ID":"","References":"<20250109115412.356768-1-stefan.klug@ideasonboard.com>\n\t<20250109115412.356768-10-stefan.klug@ideasonboard.com>","MIME-Version":"1.0","Content-Type":"text/plain; charset=us-ascii","Content-Disposition":"inline","In-Reply-To":"<20250109115412.356768-10-stefan.klug@ideasonboard.com>","X-BeenThere":"libcamera-devel@lists.libcamera.org","X-Mailman-Version":"2.1.29","Precedence":"list","List-Id":"","List-Unsubscribe":",\n\t","List-Archive":"","List-Post":"","List-Help":"","List-Subscribe":",\n\t","Errors-To":"libcamera-devel-bounces@lists.libcamera.org","Sender":"\"libcamera-devel\" "}},{"id":33113,"web_url":"https://patchwork.libcamera.org/comment/33113/","msgid":"<6efffa60-5f9f-4b99-bda4-0922f19985e8@ideasonboard.com>","date":"2025-01-21T10:57:13","subject":"Re: [PATCH v1 09/11] libipa: Add bayesian AWB algorithm","submitter":{"id":156,"url":"https://patchwork.libcamera.org/api/people/156/","name":"Dan Scally","email":"dan.scally@ideasonboard.com"},"content":"On 09/01/2025 11:54, Stefan Klug wrote:\n> The bayesian AWB algorithm is an AWB algorithm that takes prior\n> probabilities for a given light source dependent on the current lux\n> level into account.\n>\n> The biggest improvement compared to the grey world model comes from the\n> search of the ideal white point on the CT curve. The algorithm walks the\n> CT curve to minimize the colour error for a given statistics. After the\n> minimium is found it additionally tries to search the area around that\n> spot and also off the curve. So even without defined prior probabilities\n> this algorithm provides much better results than the grey world\n> algorithm.\n>\n> The logic for this code was taken from the RaspberryPi implementation.\n> The logic was only minimally adjusted for usage with the rkisp1 and a\n> few things were left out (see doxygen doc for the AwbBayes class). The\n> code is refactored to better fit the libcamera code style and to make\n> use of the syntactic sugar provided by the Interpolator and Vector\n> classes.\n>\n> Signed-off-by: Stefan Klug \n> ---\n> src/ipa/libipa/awb_bayes.cpp | 457 +++++++++++++++++++++++++++++++++++\n> src/ipa/libipa/awb_bayes.h | 67 +++++\n> src/ipa/libipa/meson.build | 2 +\n> 3 files changed, 526 insertions(+)\n> create mode 100644 src/ipa/libipa/awb_bayes.cpp\n> create mode 100644 src/ipa/libipa/awb_bayes.h\n>\n> diff --git a/src/ipa/libipa/awb_bayes.cpp b/src/ipa/libipa/awb_bayes.cpp\n> new file mode 100644\n> index 000000000000..1e69ecd3e3f3\n> --- /dev/null\n> +++ b/src/ipa/libipa/awb_bayes.cpp\n> @@ -0,0 +1,457 @@\n> +/* SPDX-License-Identifier: LGPL-2.1-or-later */\n> +/*\n> + * Copyright (C) 2019, Raspberry Pi Ltd\n> + * Copyright (C) 2024 Ideas on Board Oy\n> + *\n> + * Implementation of a bayesian AWB algorithm\n> + */\n> +\n> +#include \"awb_bayes.h\"\n> +\n> +#include \n> +\n> +#include \n> +#include \n> +\n> +#include \"colours.h\"\n> +\n> +/**\n> + * \\file awb_bayes.h\n> + * \\brief Implementation of bayesian auto white balance algorithm\n> + *\n> + * This implementation is based on the initial implementation done by\n> + * RaspberryPi.\n> + * \\todo: Documentation\n> + *\n> + * \\todo As the statistics module of the rkisp1 provides less data than the one\n> + * from the RaspberryPi (vc4). The RaspberryPi statistics measure a grid of\n> + * zones while the rkisp1 ony measures a single area.\n\n\nI don't think I'd mention the rkisp1 here...perhaps just that not all the features of the \nRaspberryPi implementation are ported over because they rely on features of its hardware that aren't \nuniversally applicable.\n\n> Therefore this algorithm\n> + * doesn't contain all the features implemented by RaspberryPi.\n> + * The following parts are not implemented:\n> + *\n> + * - min_pixels: minimum proportion of pixels counted within AWB region for it\n> + * to be \"useful\"\n> + * - min_g: minimum G value of those pixels, to be regarded a \"useful\"\n> + * - min_regions: number of AWB regions that must be \"useful\" in order to do the\n> + * AWB calculation\n> + * - deltaLimit: clamp on colour error term (so as not to penalize non-grey\n> + * excessively)\n> + * - bias_proportion: The biasProportion parameter adds a small proportion of\n> + * the counted pixels to a region biased to the biasCT colour temperature.\n> + * A typical value for biasProportion would be between 0.05 to 0.1.\n> + * - bias_ct: CT target for the search bias\n> + * - sensitivityR: red sensitivity ratio (set to canonical sensor's R/G divided\n> + * by this sensor's R/G)\n> + * - sensitivityB: blue sensitivity ratio (set to canonical sensor's B/G divided\n> + * by this sensor's B/G)\n> + */\n> +\n> +namespace libcamera {\n> +\n> +LOG_DECLARE_CATEGORY(Awb)\n> +\n> +namespace ipa {\n> +\n> +/**\n> + * \\brief Step size control for CT search\n> + */\n> +constexpr double kSearchStep = 0.2;\n> +\n> +/**\n> + * \\copydoc Interpolator::interpolate()\n> + */\n> +template<>\n> +void Interpolator::interpolate(const Pwl &a, const Pwl &b, Pwl &dest, double lambda)\n> +{\n> +\tdest = Pwl::combine(a, b,\n> +\t\t\t [=](double /*x*/, double y0, double y1) -> double {\nWhy is x commented here?\n> +\t\t\t\t return y0 * (1.0 - lambda) + y1 * lambda;\n> +\t\t\t });\n> +}\n> +\n> +/**\n> + * \\class AwbBayes\n> + * \\brief Implementation of a bayesian auto white balance algorithm\n> + *\n> + * In a bayesian AWB algorithm the auto white balance estimation is improved by\n> + * taking the likelihood of a given lightsource based on the estimated lux level\n> + * into account. E.g. If it is very bright we can assume that we are outside and\n> + * that colour temperatures around 6500 are preferred.\n> + *\n> + * The second part of this algorithm is the search for the most likely colour\n> + * temperature. It is implemented in AwbBayes::coarseSearch() and in\n> + * AwbBayes::fineSearch(). The search works very well without prior likelihoods\n> + * and therefore the algorithm itself provides very good results even without\n> + * prior likelihoods.\n> + */\n> +\n> +/**\n> + * \\var AwbBayes::transversePos_\n> + * \\brief How far to wander off CT curve towards \"more purple\"\n> + */\n> +\n> +/**\n> + * \\var AwbBayes::transverseNeg_\n> + * \\brief How far to wander off CT curve towards \"more green\"\n> + */\n> +\n> +/**\n> + * \\var AwbBayes::currentMode_\n> + * \\brief The currently selected mode\n> + */\n> +\nPersonally I think that that some example in the documentation of the tuning data format would be \nnice...if we're expecting people to stick to using libtuning to generate them then that might seem \nlike noise that's unnecessary but I think it would make it much easier to follow the function with \nthe format available. Up to you though.\n> +int AwbBayes::init(const YamlObject &tuningData)\n> +{\n> +\tint ret = colourGainCurve_.readYaml(tuningData[\"colourGains\"], \"ct\", \"gains\");\n> +\tif (ret) {\n> +\t\tLOG(Awb, Error)\n> +\t\t\t<< \"Failed to parse 'colourGains' \"\n> +\t\t\t<< \"parameter from tuning file\";\n> +\t\treturn ret;\n> +\t}\n> +\n> +\tctR_.clear();\n> +\tctB_.clear();\nThis shouldn't be necessary should it?\n> +\tfor (const auto &[ct, g] : colourGainCurve_.data()) {\n> +\t\tctR_.append(ct, 1.0 / g[0]);\n> +\t\tctB_.append(ct, 1.0 / g[1]);\n> +\t}\n> +\n> +\t/* We will want the inverse functions of these too. */\n> +\tctRInverse_ = ctR_.inverse().first;\n> +\tctBInverse_ = ctB_.inverse().first;\n> +\n> +\tret = readPriors(tuningData);\n> +\tif (ret) {\n> +\t\tLOG(Awb, Error) << \"Failed to read priors\";\n> +\t\treturn ret;\n> +\t}\n> +\n> +\tret = parseModeConfigs(tuningData);\n> +\tif (ret) {\n> +\t\tLOG(Awb, Error)\n> +\t\t\t<< \"Failed to parse mode parameter from tuning file\";\n> +\t\treturn ret;\n> +\t}\n> +\tcurrentMode_ = &modes_[controls::AwbAuto];\n> +\n> +\ttransversePos_ = tuningData[\"transversePos\"].get(0.01);\n> +\ttransverseNeg_ = tuningData[\"transverseNeg\"].get(0.01);\n> +\tif (transversePos_ <= 0 || transverseNeg_ <= 0) {\n> +\t\tLOG(Awb, Error) << \"AwbConfig: transversePos/Neg must be > 0\";\nShouldn't need \"AwbConfig\" here I think\n> +\t\treturn -EINVAL;\n> +\t}\n> +\n> +\treturn 0;\n> +}\n> +\n> +int AwbBayes::readPriors(const YamlObject &tuningData)\n> +{\n> +\tconst auto &priorsList = tuningData[\"priors\"];\n> +\tstd::map priors;\n> +\n> +\tif (!priorsList) {\n> +\t\tLOG(Awb, Error) << \"No priors specified\";\n> +\t\treturn -EINVAL;\n> +\t}\n> +\n> +\tfor (const auto &p : priorsList.asList()) {\n> +\t\tif (!p.contains(\"lux\")) {\n> +\t\t\tLOG(Awb, Error) << \"Missing lux value\";\n> +\t\t\treturn -EINVAL;\n> +\t\t}\n> +\n> +\t\tuint32_t lux = p[\"lux\"].get(0);\n> +\t\tif (priors.count(lux)) {\n> +\t\t\tLOG(Awb, Error) << \"Duplicate prior for lux value \" << lux;\n> +\t\t\treturn -EINVAL;\n> +\t\t}\n> +\n> +\t\tstd::vector temperatures =\n> +\t\t\tp[\"ct\"].getList().value_or(std::vector{});\n> +\t\tstd::vector probabilites =\n> +\t\t\tp[\"probability\"].getList().value_or(std::vector{});\n> +\n> +\t\tif (temperatures.size() != probabilites.size()) {\n> +\t\t\tLOG(Awb, Error)\n> +\t\t\t\t<< \"Ct and probability array sizes are unequal\";\n> +\t\t\treturn -EINVAL;\n> +\t\t}\n> +\n> +\t\tif (temperatures.empty()) {\n> +\t\t\tLOG(Awb, Error)\n> +\t\t\t\t<< \"Ct and probability arrays are empty\";\n> +\t\t\treturn -EINVAL;\n> +\t\t}\n> +\n> +\t\tstd::map ctToProbability;\n> +\t\tfor (unsigned int i = 0; i < temperatures.size(); i++) {\n> +\t\t\tint t = temperatures[i];\n> +\t\t\tif (ctToProbability.count(t)) {\n> +\t\t\t\tLOG(Awb, Error) << \"Duplicate ct value\";\n> +\t\t\t\treturn -EINVAL;\n> +\t\t\t}\n> +\n> +\t\t\tctToProbability[t] = probabilites[i];\n> +\t\t}\n> +\n> +\t\tauto &pwl = priors[lux];\n> +\t\tfor (const auto &[ct, prob] : ctToProbability) {\n> +\t\t\tpwl.append(ct, prob);\n> +\t\t}\n> +\t}\n> +\n> +\tif (priors.empty()) {\n> +\t\tLOG(Awb, Error) << \"No priors specified\";\n> +\t\t;\n> +\t\treturn -EINVAL;\n> +\t}\n> +\n> +\tpriors_.setData(std::move(priors));\n> +\n> +\treturn 0;\n> +}\n> +\n> +void AwbBayes::handleControls(const ControlList &controls)\n> +{\n> +\tauto mode = controls.get(controls::AwbMode);\n> +\tif (mode) {\n> +\t\tauto it = modes_.find(static_cast(*mode));\n> +\t\tif (it != modes_.end())\n> +\t\t\tcurrentMode_ = &it->second;\n> +\t\telse\n> +\t\t\tLOG(Awb, Error) << \"Unknown AWB mode\";\n> +\t}\n> +}\n\n\nThis could be part of the base class I think, and I think I would lean towards \"Unsupported\" rather \nthan \"Unknown\", as it seems more likely we'll be asked to set a mode the tuning data doesn't \ninstantiate rather than one we don't recognise at all.\n\n> +\n> +RGB AwbBayes::gainsFromColourTemperature(double colourTemperature)\n> +{\n> +\t/*\n> +\t * \\todo: In the RaspberryPi code, the ct curve was interpolated in\n> +\t * the white point space (1/x) not in gains space. This feels counter\n> +\t * intuitive, as the gains are in linear space. But I can't prove it.\n> +\t */\n> +\tconst auto &gains = colourGainCurve_.getInterpolated(colourTemperature);\n> +\treturn { { gains[0], 1.0, gains[1] } };\n> +}\n> +\n> +AwbResult AwbBayes::calculateAwb(const AwbStats &stats, int lux)\n> +{\n> +\tipa::Pwl prior;\n> +\tif (lux > 0) {\n> +\t\tprior = priors_.getInterpolated(lux);\n> +\t\tprior.map([](double x, double y) {\n> +\t\t\tLOG(Awb, Debug) << \"(\" << x << \",\" << y << \")\";\n> +\t\t});\n> +\t} else {\n> +\t\tprior.append(0, 1.0);\n> +\t}\n> +\n> +\tdouble t = coarseSearch(prior, stats);\n> +\tdouble r = ctR_.eval(t);\n> +\tdouble b = ctB_.eval(t);\n> +\tLOG(Awb, Debug)\n> +\t\t<< \"After coarse search: r \" << r << \" b \" << b << \" (gains r \"\n> +\t\t<< 1 / r << \" b \" << 1 / b << \")\";\n> +\n> +\t/*\n> +\t * Original comment from RaspberryPi:\n> +\t * Not entirely sure how to handle the fine search yet. Mostly the\n> +\t * estimated CT is already good enough, but the fine search allows us to\n> +\t * wander transversely off the CT curve. Under some illuminants, where\n> +\t * there may be more or less green light, this may prove beneficial,\n> +\t * though I probably need more real datasets before deciding exactly how\n> +\t * this should be controlled and tuned.\n> +\t */\n> +\tfineSearch(t, r, b, prior, stats);\n> +\tLOG(Awb, Debug)\n> +\t\t<< \"After fine search: r \" << r << \" b \" << b << \" (gains r \"\n> +\t\t<< 1 / r << \" b \" << 1 / b << \")\";\n> +\n> +\treturn { { { 1.0 / r, 1.0, 1.0 / b } }, t };\n> +}\n> +\n> +double AwbBayes::coarseSearch(const ipa::Pwl &prior, const AwbStats &stats) const\n> +{\n> +\tstd::vector points;\n> +\tsize_t bestPoint = 0;\n> +\tdouble t = currentMode_->ctLo;\n> +\tint spanR = -1;\n> +\tint spanB = -1;\n> +\n> +\t/* Step down the CT curve evaluating log likelihood. */\n> +\twhile (true) {\n> +\t\tdouble r = ctR_.eval(t, &spanR);\n> +\t\tdouble b = ctB_.eval(t, &spanB);\n> +\t\tRGB gains({ 1 / r, 1.0, 1 / b });\n> +\t\tdouble delta2Sum = stats.computeColourError(gains);\n> +\t\tdouble priorLogLikelihood = prior.eval(prior.domain().clamp(t));\n> +\t\tdouble finalLogLikelihood = delta2Sum - priorLogLikelihood;\n> +\n> +\t\tLOG(Awb, Debug) << \"Coarse search t: \" << t\n> +\t\t\t\t<< \" gains: \" << gains\n> +\t\t\t\t<< \" error: \" << delta2Sum\n> +\t\t\t\t<< \" prior: \" << priorLogLikelihood\n> +\t\t\t\t<< \" likelihood: \" << finalLogLikelihood;\n> +\n> +\t\tpoints.push_back({ { t, finalLogLikelihood } });\n> +\t\tif (points.back().y() < points[bestPoint].y())\n> +\t\t\tbestPoint = points.size() - 1;\n> +\n> +\t\tif (t == currentMode_->ctHi)\n> +\t\t\tbreak;\n> +\n> +\t\t/*\n> +\t\t * Ensure even steps along the r/b curve by scaling them by the\n> +\t\t * current t.\n> +\t\t */\n> +\t\tt = std::min(t + t / 10 * kSearchStep, currentMode_->ctHi);\n> +\t}\n> +\n> +\tt = points[bestPoint].x();\n> +\tLOG(Awb, Debug) << \"Coarse search found CT \" << t;\n> +\n> +\t/*\n> +\t * We have the best point of the search, but refine it with a quadratic\n> +\t * interpolation around its neighbors.\n> +\t */\n> +\tif (points.size() > 2) {\n> +\t\tbestPoint = std::clamp(bestPoint, 1ul, points.size() - 2);\n> +\t\tt = interpolateQuadratic(points[bestPoint - 1],\n> +\t\t\t\t\t points[bestPoint],\n> +\t\t\t\t\t points[bestPoint + 1]);\n> +\t\tLOG(Awb, Debug)\n> +\t\t\t<< \"After quadratic refinement, coarse search has CT \"\n> +\t\t\t<< t;\n> +\t}\n> +\n> +\treturn t;\n> +}\n> +\n> +void AwbBayes::fineSearch(double &t, double &r, double &b, ipa::Pwl const &prior, const AwbStats &stats) const\n> +{\n> +\tint spanR = -1;\n> +\tint spanB = -1;\n> +\tdouble step = t / 10 * kSearchStep * 0.1;\n> +\tint nsteps = 5;\n> +\tctR_.eval(t, &spanR);\n> +\tctB_.eval(t, &spanB);\n> +\tdouble rDiff = ctR_.eval(t + nsteps * step, &spanR) -\n> +\t\t ctR_.eval(t - nsteps * step, &spanR);\n> +\tdouble bDiff = ctB_.eval(t + nsteps * step, &spanB) -\n> +\t\t ctB_.eval(t - nsteps * step, &spanB);\n> +\tPwl::Point transverse({ bDiff, -rDiff });\n> +\tif (transverse.length2() < 1e-6)\n> +\t\treturn;\n> +\t/*\n> +\t * transverse is a unit vector orthogonal to the b vs. r function\n> +\t * (pointing outwards with r and b increasing)\n> +\t */\n> +\ttransverse = transverse / transverse.length();\n> +\tdouble bestLogLikelihood = 0;\n> +\tdouble bestT = 0;\n> +\tPwl::Point bestRB;\n> +\tdouble transverseRange = transverseNeg_ + transversePos_;\n> +\tconst int maxNumDeltas = 12;\n> +\n> +\t/* a transverse step approximately every 0.01 r/b units */\n> +\tint numDeltas = floor(transverseRange * 100 + 0.5) + 1;\n> +\tnumDeltas = std::clamp(numDeltas, 3, maxNumDeltas);\n> +\n> +\t/*\n> +\t * Step down CT curve. March a bit further if the transverse range is\n> +\t * large.\n> +\t */\n> +\tnsteps += numDeltas;\n> +\tfor (int i = -nsteps; i <= nsteps; i++) {\n> +\t\tdouble tTest = t + i * step;\n> +\t\tdouble priorLogLikelihood =\n> +\t\t\tprior.eval(prior.domain().clamp(tTest));\n> +\t\tPwl::Point rbStart{ { ctR_.eval(tTest, &spanR),\n> +\t\t\t\t ctB_.eval(tTest, &spanB) } };\n> +\t\tPwl::Point samples[maxNumDeltas];\n> +\t\tint bestPoint = 0;\n> +\n> +\t\t/*\n> +\t\t * Sample numDeltas points transversely *off* the CT curve\n> +\t\t * in the range [-transverseNeg, transversePos].\n> +\t\t * The x values of a sample contains the distance and the y\n> +\t\t * value contains the corresponding log likelihood.\n> +\t\t */\n> +\t\tdouble transverseStep = transverseRange / (numDeltas - 1);\n> +\t\tfor (int j = 0; j < numDeltas; j++) {\n> +\t\t\tauto &p = samples[j];\n> +\t\t\tp.x() = -transverseNeg_ + transverseStep * j;\n> +\t\t\tPwl::Point rbTest = rbStart + transverse * p.x();\n> +\t\t\tRGB gains({ 1 / rbTest[0], 1.0, 1 / rbTest[1] });\n> +\t\t\tdouble delta2Sum = stats.computeColourError(gains);\n> +\t\t\tp.y() = delta2Sum - priorLogLikelihood;\n> +\n> +\t\t\tif (p.y() < samples[bestPoint].y())\n> +\t\t\t\tbestPoint = j;\n> +\t\t}\n> +\n> +\t\t/*\n> +\t\t * We have all samples transversely across the CT curve,\n> +\t\t * now let's do a quadratic interpolation for the best result.\n> +\t\t */\n> +\t\tbestPoint = std::clamp(bestPoint, 1, numDeltas - 2);\n> +\t\tdouble bestOffset = interpolateQuadratic(samples[bestPoint - 1],\n> +\t\t\t\t\t\t\t samples[bestPoint],\n> +\t\t\t\t\t\t\t samples[bestPoint + 1]);\n> +\t\tPwl::Point rbTest = rbStart + transverse * bestOffset;\n> +\t\tRGB gains({ 1 / rbTest[0], 1.0, 1 / rbTest[1] });\n> +\t\tdouble delta2Sum = stats.computeColourError(gains);\n> +\t\tdouble finalLogLikelihood = delta2Sum - priorLogLikelihood;\n> +\t\tLOG(Awb, Debug)\n> +\t\t\t<< \"Fine search t: \" << tTest\n> +\t\t\t<< \" r: \" << rbTest[0]\n> +\t\t\t<< \" b: \" << rbTest[1]\n> +\t\t\t<< \" offset: \" << bestOffset\n> +\t\t\t<< \" likelihood: \" << finalLogLikelihood\n> +\t\t\t<< (finalLogLikelihood < bestLogLikelihood ? \" NEW BEST\" : \"\");\n> +\n> +\t\tif (bestT == 0 || finalLogLikelihood < bestLogLikelihood) {\n> +\t\t\tbestLogLikelihood = finalLogLikelihood;\n> +\t\t\tbestT = tTest;\n> +\t\t\tbestRB = rbTest;\n> +\t\t}\n> +\t}\n> +\n> +\tt = bestT;\n> +\tr = bestRB[0];\n> +\tb = bestRB[1];\n> +\tLOG(Awb, Debug)\n> +\t\t<< \"Fine search found t \" << t << \" r \" << r << \" b \" << b;\n> +}\n> +\n> +/**\n> + * \\brief Find extremum of function\n> + * \\param[in] a First point\n> + * \\param[in] b Second point\n> + * \\param[in] c Third point\n> + *\n> + * Given 3 points on a curve, find the extremum of the function in that interval\n> + * by fitting a quadratic.\n> + *\n> + * \\return The x value of the extremum clamped to the interval [a.x(), c.x()]\n> + */\n> +double AwbBayes::interpolateQuadratic(Pwl::Point const &a, Pwl::Point const &b,\n> +\t\t\t\t Pwl::Point const &c) const\n> +{\n> +\tconst double eps = 1e-3;\n> +\tPwl::Point ca = c - a;\n> +\tPwl::Point ba = b - a;\n> +\tdouble denominator = 2 * (ba.y() * ca.x() - ca.y() * ba.x());\n> +\tif (std::abs(denominator) > eps) {\n> +\t\tdouble numerator = ba.y() * ca.x() * ca.x() - ca.y() * ba.x() * ba.x();\n> +\t\tdouble result = numerator / denominator + a.x();\n> +\t\treturn std::max(a.x(), std::min(c.x(), result));\n> +\t}\n> +\t/* has degenerated to straight line segment */\n> +\treturn a.y() < c.y() - eps ? a.x() : (c.y() < a.y() - eps ? c.x() : b.x());\n> +}\n\nA candidate for a helper perhaps, though doesn't have to be done now.\n\n\nAll really minor points:\n\n\nReviewed-by: Daniel Scally \n\n> +\n> +} /* namespace ipa */\n> +\n> +} /* namespace libcamera */\n> diff --git a/src/ipa/libipa/awb_bayes.h b/src/ipa/libipa/awb_bayes.h\n> new file mode 100644\n> index 000000000000..5284f0ca4cc0\n> --- /dev/null\n> +++ b/src/ipa/libipa/awb_bayes.h\n> @@ -0,0 +1,67 @@\n> +/* SPDX-License-Identifier: LGPL-2.1-or-later */\n> +/*\n> + * Copyright (C) 2024 Ideas on Board Oy\n> + *\n> + * Base class for bayes AWB algorithms\n> + */\n> +\n> +#pragma once\n> +\n> +#include \n> +#include \n> +#include \n> +#include \n> +\n> +#include \n> +\n> +#include \n> +#include \n> +\n> +#include \"libcamera/internal/yaml_parser.h\"\n> +\n> +#include \"awb.h\"\n> +#include \"interpolator.h\"\n> +#include \"pwl.h\"\n> +#include \"vector.h\"\n> +\n> +namespace libcamera {\n> +\n> +namespace ipa {\n> +\n> +class AwbBayes : public AwbAlgorithm\n> +{\n> +public:\n> +\tAwbBayes() = default;\n> +\n> +\tint init(const YamlObject &tuningData) override;\n> +\tAwbResult calculateAwb(const AwbStats &stats, int lux) override;\n> +\tRGB gainsFromColourTemperature(double temperatureK) override;\n> +\tvoid handleControls(const ControlList &controls) override;\n> +\n> +private:\n> +\tint readPriors(const YamlObject &tuningData);\n> +\n> +\tvoid fineSearch(double &t, double &r, double &b, ipa::Pwl const &prior,\n> +\t\t\tconst AwbStats &stats) const;\n> +\tdouble coarseSearch(const ipa::Pwl &prior, const AwbStats &stats) const;\n> +\tdouble interpolateQuadratic(ipa::Pwl::Point const &a,\n> +\t\t\t\t ipa::Pwl::Point const &b,\n> +\t\t\t\t ipa::Pwl::Point const &c) const;\n> +\n> +\tInterpolator priors_;\n> +\tInterpolator> colourGainCurve_;\n> +\n> +\tipa::Pwl ctR_;\n> +\tipa::Pwl ctB_;\n> +\tipa::Pwl ctRInverse_;\n> +\tipa::Pwl ctBInverse_;\n> +\n> +\tdouble transversePos_;\n> +\tdouble transverseNeg_;\n> +\n> +\tModeConfig *currentMode_ = nullptr;\n> +};\n> +\n> +} /* namespace ipa */\n> +\n> +} /* namespace libcamera */\n> diff --git a/src/ipa/libipa/meson.build b/src/ipa/libipa/meson.build\n> index c550a6eb45b6..b519c8146d7c 100644\n> --- a/src/ipa/libipa/meson.build\n> +++ b/src/ipa/libipa/meson.build\n> @@ -3,6 +3,7 @@\n> libipa_headers = files([\n> 'agc_mean_luminance.h',\n> 'algorithm.h',\n> + 'awb_bayes.h',\n> 'awb_grey.h',\n> 'awb.h',\n> 'camera_sensor_helper.h',\n> @@ -22,6 +23,7 @@ libipa_headers = files([\n> libipa_sources = files([\n> 'agc_mean_luminance.cpp',\n> 'algorithm.cpp',\n> + 'awb_bayes.cpp',\n> 'awb_grey.cpp',\n> 'awb.cpp',\n> 'camera_sensor_helper.cpp',","headers":{"Return-Path":"","X-Original-To":"parsemail@patchwork.libcamera.org","Delivered-To":"parsemail@patchwork.libcamera.org","Received":["from lancelot.ideasonboard.com (lancelot.ideasonboard.com\n\t[92.243.16.209])\n\tby patchwork.libcamera.org (Postfix) with ESMTPS id C2F9CBD160\n\tfor ;\n\tTue, 21 Jan 2025 10:57:19 +0000 (UTC)","from lancelot.ideasonboard.com (localhost [IPv6:::1])\n\tby lancelot.ideasonboard.com (Postfix) with ESMTP id E4D1668543;\n\tTue, 21 Jan 2025 11:57:18 +0100 (CET)","from perceval.ideasonboard.com (perceval.ideasonboard.com\n\t[213.167.242.64])\n\tby lancelot.ideasonboard.com (Postfix) with ESMTPS id 8A8DB68516\n\tfor ;\n\tTue, 21 Jan 2025 11:57:16 +0100 (CET)","from [192.168.0.43]\n\t(cpc141996-chfd3-2-0-cust928.12-3.cable.virginm.net [86.13.91.161])\n\tby perceval.ideasonboard.com (Postfix) with ESMTPSA id 607DF788\n\tfor ;\n\tTue, 21 Jan 2025 11:56:14 +0100 (CET)"],"Authentication-Results":"lancelot.ideasonboard.com; dkim=pass (1024-bit key;\n\tunprotected) header.d=ideasonboard.com header.i=@ideasonboard.com\n\theader.b=\"iNoTrD+O\"; dkim-atps=neutral","DKIM-Signature":"v=1; a=rsa-sha256; c=relaxed/simple; d=ideasonboard.com;\n\ts=mail; t=1737456974;\n\tbh=LuDfnkJn40O7H0kUIi8JykqOXK5PrpH8HgHDAbomRu4=;\n\th=Date:Subject:To:References:From:In-Reply-To:From;\n\tb=iNoTrD+O/ja7sBE1e6QLPOee0LeYOMkPEBx5nJjuOdV5QDW1zUu2xO16JF1S/i+nw\n\tooweVBR8+oAXdnM4h41dPR7AAZpGoxMDQ3OoG2yDoETelPLZD3JRnKfXu6YWd4HoU5\n\t8wiLngiP6QnvqmxeZt3mfYahRF/xuGbi1QfIDOGg=","Message-ID":"<6efffa60-5f9f-4b99-bda4-0922f19985e8@ideasonboard.com>","Date":"Tue, 21 Jan 2025 10:57:13 +0000","MIME-Version":"1.0","User-Agent":"Mozilla Thunderbird","Subject":"Re: [PATCH v1 09/11] libipa: Add bayesian AWB algorithm","To":"libcamera-devel@lists.libcamera.org","References":"<20250109115412.356768-1-stefan.klug@ideasonboard.com>\n\t<20250109115412.356768-10-stefan.klug@ideasonboard.com>","Content-Language":"en-US","From":"Dan Scally ","Autocrypt":"addr=dan.scally@ideasonboard.com; keydata=\n\txsFNBGLydlEBEADa5O2s0AbUguprfvXOQun/0a8y2Vk6BqkQALgeD6KnXSWwaoCULp18etYW\n\tB31bfgrdphXQ5kUQibB0ADK8DERB4wrzrUb5CMxLBFE7mQty+v5NsP0OFNK9XTaAOcmD+Ove\n\teIjYvqurAaro91jrRVrS1gBRxIFqyPgNvwwL+alMZhn3/2jU2uvBmuRrgnc/e9cHKiuT3Dtq\n\tMHGPKL2m+plk+7tjMoQFfexoQ1JKugHAjxAhJfrkXh6uS6rc01bYCyo7ybzg53m1HLFJdNGX\n\tsUKR+dQpBs3SY4s66tc1sREJqdYyTsSZf80HjIeJjU/hRunRo4NjRIJwhvnK1GyjOvvuCKVU\n\tRWpY8dNjNu5OeAfdrlvFJOxIE9M8JuYCQTMULqd1NuzbpFMjc9524U3Cngs589T7qUMPb1H1\n\tNTA81LmtJ6Y+IV5/kiTUANflpzBwhu18Ok7kGyCq2a2jsOcVmk8gZNs04gyjuj8JziYwwLbf\n\tvzABwpFVcS8aR+nHIZV1HtOzyw8CsL8OySc3K9y+Y0NRpziMRvutrppzgyMb9V+N31mK9Mxl\n\t1YkgaTl4ciNWpdfUe0yxH03OCuHi3922qhPLF4XX5LN+NaVw5Xz2o3eeWklXdouxwV7QlN33\n\tu4+u2FWzKxDqO6WLQGjxPE0mVB4Gh5Pa1Vb0ct9Ctg0qElvtGQARAQABzShEYW4gU2NhbGx5\n\tIDxkYW4uc2NhbGx5QGlkZWFzb25ib2FyZC5jb20+wsGNBBMBCAA3FiEEsdtt8OWP7+8SNfQe\n\tkiQuh/L+GMQFAmLydlIFCQWjmoACGwMECwkIBwUVCAkKCwUWAgMBAAAKCRCSJC6H8v4YxDI2\n\tEAC2Gz0iyaXJkPInyshrREEWbo0CA6v5KKf3I/HlMPqkZ48bmGoYm4mEQGFWZJAT3K4ir8bg\n\tcEfs9V54gpbrZvdwS4abXbUK4WjKwEs8HK3XJv1WXUN2bsz5oEJWZUImh9gD3naiLLI9QMMm\n\tw/aZkT+NbN5/2KvChRWhdcha7+2Te4foOY66nIM+pw2FZM6zIkInLLUik2zXOhaZtqdeJZQi\n\tHSPU9xu7TRYN4cvdZAnSpG7gQqmLm5/uGZN1/sB3kHTustQtSXKMaIcD/DMNI3JN/t+RJVS7\n\tc0Jh/ThzTmhHyhxx3DRnDIy7kwMI4CFvmhkVC2uNs9kWsj1DuX5kt8513mvfw2OcX9UnNKmZ\n\tnhNCuF6DxVrL8wjOPuIpiEj3V+K7DFF1Cxw1/yrLs8dYdYh8T8vCY2CHBMsqpESROnTazboh\n\tAiQ2xMN1cyXtX11Qwqm5U3sykpLbx2BcmUUUEAKNsM//Zn81QXKG8vOx0ZdMfnzsCaCzt8f6\n\t9dcDBBI3tJ0BI9ByiocqUoL6759LM8qm18x3FYlxvuOs4wSGPfRVaA4yh0pgI+ModVC2Pu3y\n\tejE/IxeatGqJHh6Y+iJzskdi27uFkRixl7YJZvPJAbEn7kzSi98u/5ReEA8Qhc8KO/B7wprj\n\txjNMZNYd0Eth8+WkixHYj752NT5qshKJXcyUU87BTQRi8nZSARAAx0BJayh1Fhwbf4zoY56x\n\txHEpT6DwdTAYAetd3yiKClLVJadYxOpuqyWa1bdfQWPb+h4MeXbWw/53PBgn7gI2EA7ebIRC\n\tPJJhAIkeym7hHZoxqDQTGDJjxFEL11qF+U3rhWiL2Zt0Pl+zFq0eWYYVNiXjsIS4FI2+4m16\n\ttPbDWZFJnSZ828VGtRDQdhXfx3zyVX21lVx1bX4/OZvIET7sVUufkE4hrbqrrufre7wsjD1t\n\t8MQKSapVrr1RltpzPpScdoxknOSBRwOvpp57pJJe5A0L7+WxJ+vQoQXj0j+5tmIWOAV1qBQp\n\thyoyUk9JpPfntk2EKnZHWaApFp5TcL6c5LhUvV7F6XwOjGPuGlZQCWXee9dr7zym8iR3irWT\n\t+49bIh5PMlqSLXJDYbuyFQHFxoiNdVvvf7etvGfqFYVMPVjipqfEQ38ST2nkzx+KBICz7uwj\n\tJwLBdTXzGFKHQNckGMl7F5QdO/35An/QcxBnHVMXqaSd12tkJmoRVWduwuuoFfkTY5mUV3uX\n\txGj3iVCK4V+ezOYA7c2YolfRCNMTza6vcK/P4tDjjsyBBZrCCzhBvd4VVsnnlZhVaIxoky4K\n\taL+AP+zcQrUZmXmgZjXOLryGnsaeoVrIFyrU6ly90s1y3KLoPsDaTBMtnOdwxPmo1xisH8oL\n\ta/VRgpFBfojLPxMAEQEAAcLBfAQYAQgAJhYhBLHbbfDlj+/vEjX0HpIkLofy/hjEBQJi8nZT\n\tBQkFo5qAAhsMAAoJEJIkLofy/hjEXPcQAMIPNqiWiz/HKu9W4QIf1OMUpKn3YkVIj3p3gvfM\n\tRes4fGX94Ji599uLNrPoxKyaytC4R6BTxVriTJjWK8mbo9jZIRM4vkwkZZ2bu98EweSucxbp\n\tvjESsvMXGgxniqV/RQ/3T7LABYRoIUutARYq58p5HwSP0frF0fdFHYdTa2g7MYZl1ur2JzOC\n\tFHRpGadlNzKDE3fEdoMobxHB3Lm6FDml5GyBAA8+dQYVI0oDwJ3gpZPZ0J5Vx9RbqXe8RDuR\n\tdu90hvCJkq7/tzSQ0GeD3BwXb9/R/A4dVXhaDd91Q1qQXidI+2jwhx8iqiYxbT+DoAUkQRQy\n\txBtoCM1CxH7u45URUgD//fxYr3D4B1SlonA6vdaEdHZOGwECnDpTxecENMbz/Bx7qfrmd901\n\tD+N9SjIwrbVhhSyUXYnSUb8F+9g2RDY42Sk7GcYxIeON4VzKqWM7hpkXZ47pkK0YodO+dRKM\n\tyMcoUWrTK0Uz6UzUGKoJVbxmSW/EJLEGoI5p3NWxWtScEVv8mO49gqQdrRIOheZycDmHnItt\n\t9Qjv00uFhEwv2YfiyGk6iGF2W40s2pH2t6oeuGgmiZ7g6d0MEK8Ql/4zPItvr1c1rpwpXUC1\n\tu1kQWgtnNjFHX3KiYdqjcZeRBiry1X0zY+4Y24wUU0KsEewJwjhmCKAsju1RpdlPg2kC","In-Reply-To":"<20250109115412.356768-10-stefan.klug@ideasonboard.com>","Content-Type":"text/plain; charset=UTF-8; format=flowed","Content-Transfer-Encoding":"7bit","X-BeenThere":"libcamera-devel@lists.libcamera.org","X-Mailman-Version":"2.1.29","Precedence":"list","List-Id":"","List-Unsubscribe":",\n\t","List-Archive":"","List-Post":"","List-Help":"","List-Subscribe":",\n\t","Errors-To":"libcamera-devel-bounces@lists.libcamera.org","Sender":"\"libcamera-devel\" "}}]