[{"id":29849,"web_url":"https://patchwork.libcamera.org/comment/29849/","msgid":"<20240611214058.GC8290@pendragon.ideasonboard.com>","date":"2024-06-11T21:40:58","subject":"Re: [PATCH v8 2/4] ipa: libipa: Copy pwl from rpi","submitter":{"id":2,"url":"https://patchwork.libcamera.org/api/people/2/","name":"Laurent Pinchart","email":"laurent.pinchart@ideasonboard.com"},"content":"Hi Paul,\n\nThank you for the patch.\n\nOn Tue, Jun 11, 2024 at 10:24:28PM +0900, Paul Elder wrote:\n> Copy the piecewise linear function code from Raspberry Pi.\n> \n> Signed-off-by: Paul Elder \n> Reviewed-by: Stefan Klug \n> Acked-by: David Plowman \n> Reviewed-by: Kieran Bingham \n\nIf we don't squash this with 3/4 before merging,\n\nReviewed-by: Laurent Pinchart \n\n> \n> ---\n> No change in v8\n> \n> No change in v7\n> \n> No change in v6\n> \n> Changes in v5:\n> - remove meson.build to prevent compilation this early in the merge\n> \n> Changes in v4:\n> - update the copy\n> \n> No change in v3\n> \n> No change in v2\n> ---\n> src/ipa/libipa/pwl.cpp | 269 +++++++++++++++++++++++++++++++++++++++++\n> src/ipa/libipa/pwl.h | 127 +++++++++++++++++++\n> 2 files changed, 396 insertions(+)\n> create mode 100644 src/ipa/libipa/pwl.cpp\n> create mode 100644 src/ipa/libipa/pwl.h\n> \n> diff --git a/src/ipa/libipa/pwl.cpp b/src/ipa/libipa/pwl.cpp\n> new file mode 100644\n> index 000000000000..e39123767aa6\n> --- /dev/null\n> +++ b/src/ipa/libipa/pwl.cpp\n> @@ -0,0 +1,269 @@\n> +/* SPDX-License-Identifier: BSD-2-Clause */\n> +/*\n> + * Copyright (C) 2019, Raspberry Pi Ltd\n> + *\n> + * piecewise linear functions\n> + */\n> +\n> +#include \n> +#include \n> +#include \n> +\n> +#include \"pwl.h\"\n> +\n> +using namespace RPiController;\n> +\n> +int Pwl::read(const libcamera::YamlObject ¶ms)\n> +{\n> +\tif (!params.size() || params.size() % 2)\n> +\t\treturn -EINVAL;\n> +\n> +\tconst auto &list = params.asList();\n> +\n> +\tfor (auto it = list.begin(); it != list.end(); it++) {\n> +\t\tauto x = it->get();\n> +\t\tif (!x)\n> +\t\t\treturn -EINVAL;\n> +\t\tif (it != list.begin() && *x <= points_.back().x)\n> +\t\t\treturn -EINVAL;\n> +\n> +\t\tauto y = (++it)->get();\n> +\t\tif (!y)\n> +\t\t\treturn -EINVAL;\n> +\n> +\t\tpoints_.push_back(Point(*x, *y));\n> +\t}\n> +\n> +\treturn 0;\n> +}\n> +\n> +void Pwl::append(double x, double y, const double eps)\n> +{\n> +\tif (points_.empty() || points_.back().x + eps < x)\n> +\t\tpoints_.push_back(Point(x, y));\n> +}\n> +\n> +void Pwl::prepend(double x, double y, const double eps)\n> +{\n> +\tif (points_.empty() || points_.front().x - eps > x)\n> +\t\tpoints_.insert(points_.begin(), Point(x, y));\n> +}\n> +\n> +Pwl::Interval Pwl::domain() const\n> +{\n> +\treturn Interval(points_[0].x, points_[points_.size() - 1].x);\n> +}\n> +\n> +Pwl::Interval Pwl::range() const\n> +{\n> +\tdouble lo = points_[0].y, hi = lo;\n> +\tfor (auto &p : points_)\n> +\t\tlo = std::min(lo, p.y), hi = std::max(hi, p.y);\n> +\treturn Interval(lo, hi);\n> +}\n> +\n> +bool Pwl::empty() const\n> +{\n> +\treturn points_.empty();\n> +}\n> +\n> +double Pwl::eval(double x, int *spanPtr, bool updateSpan) const\n> +{\n> +\tint span = findSpan(x, spanPtr && *spanPtr != -1 ? *spanPtr : points_.size() / 2 - 1);\n> +\tif (spanPtr && updateSpan)\n> +\t\t*spanPtr = span;\n> +\treturn points_[span].y +\n> +\t (x - points_[span].x) * (points_[span + 1].y - points_[span].y) /\n> +\t\t (points_[span + 1].x - points_[span].x);\n> +}\n> +\n> +int Pwl::findSpan(double x, int span) const\n> +{\n> +\t/*\n> +\t * Pwls are generally small, so linear search may well be faster than\n> +\t * binary, though could review this if large PWls start turning up.\n> +\t */\n> +\tint lastSpan = points_.size() - 2;\n> +\t/*\n> +\t * some algorithms may call us with span pointing directly at the last\n> +\t * control point\n> +\t */\n> +\tspan = std::max(0, std::min(lastSpan, span));\n> +\twhile (span < lastSpan && x >= points_[span + 1].x)\n> +\t\tspan++;\n> +\twhile (span && x < points_[span].x)\n> +\t\tspan--;\n> +\treturn span;\n> +}\n> +\n> +Pwl::PerpType Pwl::invert(Point const &xy, Point &perp, int &span,\n> +\t\t\t const double eps) const\n> +{\n> +\tassert(span >= -1);\n> +\tbool prevOffEnd = false;\n> +\tfor (span = span + 1; span < (int)points_.size() - 1; span++) {\n> +\t\tPoint spanVec = points_[span + 1] - points_[span];\n> +\t\tdouble t = ((xy - points_[span]) % spanVec) / spanVec.len2();\n> +\t\tif (t < -eps) /* off the start of this span */\n> +\t\t{\n> +\t\t\tif (span == 0) {\n> +\t\t\t\tperp = points_[span];\n> +\t\t\t\treturn PerpType::Start;\n> +\t\t\t} else if (prevOffEnd) {\n> +\t\t\t\tperp = points_[span];\n> +\t\t\t\treturn PerpType::Vertex;\n> +\t\t\t}\n> +\t\t} else if (t > 1 + eps) /* off the end of this span */\n> +\t\t{\n> +\t\t\tif (span == (int)points_.size() - 2) {\n> +\t\t\t\tperp = points_[span + 1];\n> +\t\t\t\treturn PerpType::End;\n> +\t\t\t}\n> +\t\t\tprevOffEnd = true;\n> +\t\t} else /* a true perpendicular */\n> +\t\t{\n> +\t\t\tperp = points_[span] + spanVec * t;\n> +\t\t\treturn PerpType::Perpendicular;\n> +\t\t}\n> +\t}\n> +\treturn PerpType::None;\n> +}\n> +\n> +Pwl Pwl::inverse(bool *trueInverse, const double eps) const\n> +{\n> +\tbool appended = false, prepended = false, neither = false;\n> +\tPwl inverse;\n> +\n> +\tfor (Point const &p : points_) {\n> +\t\tif (inverse.empty())\n> +\t\t\tinverse.append(p.y, p.x, eps);\n> +\t\telse if (std::abs(inverse.points_.back().x - p.y) <= eps ||\n> +\t\t\t std::abs(inverse.points_.front().x - p.y) <= eps)\n> +\t\t\t/* do nothing */;\n> +\t\telse if (p.y > inverse.points_.back().x) {\n> +\t\t\tinverse.append(p.y, p.x, eps);\n> +\t\t\tappended = true;\n> +\t\t} else if (p.y < inverse.points_.front().x) {\n> +\t\t\tinverse.prepend(p.y, p.x, eps);\n> +\t\t\tprepended = true;\n> +\t\t} else\n> +\t\t\tneither = true;\n> +\t}\n> +\n> +\t/*\n> +\t * This is not a proper inverse if we found ourselves putting points\n> +\t * onto both ends of the inverse, or if there were points that couldn't\n> +\t * go on either.\n> +\t */\n> +\tif (trueInverse)\n> +\t\t*trueInverse = !(neither || (appended && prepended));\n> +\n> +\treturn inverse;\n> +}\n> +\n> +Pwl Pwl::compose(Pwl const &other, const double eps) const\n> +{\n> +\tdouble thisX = points_[0].x, thisY = points_[0].y;\n> +\tint thisSpan = 0, otherSpan = other.findSpan(thisY, 0);\n> +\tPwl result({ { thisX, other.eval(thisY, &otherSpan, false) } });\n> +\twhile (thisSpan != (int)points_.size() - 1) {\n> +\t\tdouble dx = points_[thisSpan + 1].x - points_[thisSpan].x,\n> +\t\t dy = points_[thisSpan + 1].y - points_[thisSpan].y;\n> +\t\tif (std::abs(dy) > eps &&\n> +\t\t otherSpan + 1 < (int)other.points_.size() &&\n> +\t\t points_[thisSpan + 1].y >=\n> +\t\t\t other.points_[otherSpan + 1].x + eps) {\n> +\t\t\t/*\n> +\t\t\t * next control point in result will be where this\n> +\t\t\t * function's y reaches the next span in other\n> +\t\t\t */\n> +\t\t\tthisX = points_[thisSpan].x +\n> +\t\t\t\t(other.points_[otherSpan + 1].x -\n> +\t\t\t\t points_[thisSpan].y) *\n> +\t\t\t\t\tdx / dy;\n> +\t\t\tthisY = other.points_[++otherSpan].x;\n> +\t\t} else if (std::abs(dy) > eps && otherSpan > 0 &&\n> +\t\t\t points_[thisSpan + 1].y <=\n> +\t\t\t\t other.points_[otherSpan - 1].x - eps) {\n> +\t\t\t/*\n> +\t\t\t * next control point in result will be where this\n> +\t\t\t * function's y reaches the previous span in other\n> +\t\t\t */\n> +\t\t\tthisX = points_[thisSpan].x +\n> +\t\t\t\t(other.points_[otherSpan + 1].x -\n> +\t\t\t\t points_[thisSpan].y) *\n> +\t\t\t\t\tdx / dy;\n> +\t\t\tthisY = other.points_[--otherSpan].x;\n> +\t\t} else {\n> +\t\t\t/* we stay in the same span in other */\n> +\t\t\tthisSpan++;\n> +\t\t\tthisX = points_[thisSpan].x,\n> +\t\t\tthisY = points_[thisSpan].y;\n> +\t\t}\n> +\t\tresult.append(thisX, other.eval(thisY, &otherSpan, false),\n> +\t\t\t eps);\n> +\t}\n> +\treturn result;\n> +}\n> +\n> +void Pwl::map(std::function f) const\n> +{\n> +\tfor (auto &pt : points_)\n> +\t\tf(pt.x, pt.y);\n> +}\n> +\n> +void Pwl::map2(Pwl const &pwl0, Pwl const &pwl1,\n> +\t std::function f)\n> +{\n> +\tint span0 = 0, span1 = 0;\n> +\tdouble x = std::min(pwl0.points_[0].x, pwl1.points_[0].x);\n> +\tf(x, pwl0.eval(x, &span0, false), pwl1.eval(x, &span1, false));\n> +\twhile (span0 < (int)pwl0.points_.size() - 1 ||\n> +\t span1 < (int)pwl1.points_.size() - 1) {\n> +\t\tif (span0 == (int)pwl0.points_.size() - 1)\n> +\t\t\tx = pwl1.points_[++span1].x;\n> +\t\telse if (span1 == (int)pwl1.points_.size() - 1)\n> +\t\t\tx = pwl0.points_[++span0].x;\n> +\t\telse if (pwl0.points_[span0 + 1].x > pwl1.points_[span1 + 1].x)\n> +\t\t\tx = pwl1.points_[++span1].x;\n> +\t\telse\n> +\t\t\tx = pwl0.points_[++span0].x;\n> +\t\tf(x, pwl0.eval(x, &span0, false), pwl1.eval(x, &span1, false));\n> +\t}\n> +}\n> +\n> +Pwl Pwl::combine(Pwl const &pwl0, Pwl const &pwl1,\n> +\t\t std::function f,\n> +\t\t const double eps)\n> +{\n> +\tPwl result;\n> +\tmap2(pwl0, pwl1, [&](double x, double y0, double y1) {\n> +\t\tresult.append(x, f(x, y0, y1), eps);\n> +\t});\n> +\treturn result;\n> +}\n> +\n> +void Pwl::matchDomain(Interval const &domain, bool clip, const double eps)\n> +{\n> +\tint span = 0;\n> +\tprepend(domain.start, eval(clip ? points_[0].x : domain.start, &span),\n> +\t\teps);\n> +\tspan = points_.size() - 2;\n> +\tappend(domain.end, eval(clip ? points_.back().x : domain.end, &span),\n> +\t eps);\n> +}\n> +\n> +Pwl &Pwl::operator*=(double d)\n> +{\n> +\tfor (auto &pt : points_)\n> +\t\tpt.y *= d;\n> +\treturn *this;\n> +}\n> +\n> +void Pwl::debug(FILE *fp) const\n> +{\n> +\tfprintf(fp, \"Pwl {\\n\");\n> +\tfor (auto &p : points_)\n> +\t\tfprintf(fp, \"\\t(%g, %g)\\n\", p.x, p.y);\n> +\tfprintf(fp, \"}\\n\");\n> +}\n> diff --git a/src/ipa/libipa/pwl.h b/src/ipa/libipa/pwl.h\n> new file mode 100644\n> index 000000000000..7d5e7e4d3fda\n> --- /dev/null\n> +++ b/src/ipa/libipa/pwl.h\n> @@ -0,0 +1,127 @@\n> +/* SPDX-License-Identifier: BSD-2-Clause */\n> +/*\n> + * Copyright (C) 2019, Raspberry Pi Ltd\n> + *\n> + * piecewise linear functions interface\n> + */\n> +#pragma once\n> +\n> +#include \n> +#include \n> +#include \n> +\n> +#include \"libcamera/internal/yaml_parser.h\"\n> +\n> +namespace RPiController {\n> +\n> +class Pwl\n> +{\n> +public:\n> +\tstruct Interval {\n> +\t\tInterval(double _start, double _end)\n> +\t\t\t: start(_start), end(_end)\n> +\t\t{\n> +\t\t}\n> +\t\tdouble start, end;\n> +\t\tbool contains(double value)\n> +\t\t{\n> +\t\t\treturn value >= start && value <= end;\n> +\t\t}\n> +\t\tdouble clip(double value)\n> +\t\t{\n> +\t\t\treturn value < start ? start\n> +\t\t\t\t\t : (value > end ? end : value);\n> +\t\t}\n> +\t\tdouble len() const { return end - start; }\n> +\t};\n> +\tstruct Point {\n> +\t\tPoint() : x(0), y(0) {}\n> +\t\tPoint(double _x, double _y)\n> +\t\t\t: x(_x), y(_y) {}\n> +\t\tdouble x, y;\n> +\t\tPoint operator-(Point const &p) const\n> +\t\t{\n> +\t\t\treturn Point(x - p.x, y - p.y);\n> +\t\t}\n> +\t\tPoint operator+(Point const &p) const\n> +\t\t{\n> +\t\t\treturn Point(x + p.x, y + p.y);\n> +\t\t}\n> +\t\tdouble operator%(Point const &p) const\n> +\t\t{\n> +\t\t\treturn x * p.x + y * p.y;\n> +\t\t}\n> +\t\tPoint operator*(double f) const { return Point(x * f, y * f); }\n> +\t\tPoint operator/(double f) const { return Point(x / f, y / f); }\n> +\t\tdouble len2() const { return x * x + y * y; }\n> +\t\tdouble len() const { return sqrt(len2()); }\n> +\t};\n> +\tPwl() {}\n> +\tPwl(std::vector const &points) : points_(points) {}\n> +\tint read(const libcamera::YamlObject ¶ms);\n> +\tvoid append(double x, double y, const double eps = 1e-6);\n> +\tvoid prepend(double x, double y, const double eps = 1e-6);\n> +\tInterval domain() const;\n> +\tInterval range() const;\n> +\tbool empty() const;\n> +\t/*\n> +\t * Evaluate Pwl, optionally supplying an initial guess for the\n> +\t * \"span\". The \"span\" may be optionally be updated. If you want to know\n> +\t * the \"span\" value but don't have an initial guess you can set it to\n> +\t * -1.\n> +\t */\n> +\tdouble eval(double x, int *spanPtr = nullptr,\n> +\t\t bool updateSpan = true) const;\n> +\t/*\n> +\t * Find perpendicular closest to xy, starting from span+1 so you can\n> +\t * call it repeatedly to check for multiple closest points (set span to\n> +\t * -1 on the first call). Also returns \"pseudo\" perpendiculars; see\n> +\t * PerpType enum.\n> +\t */\n> +\tenum class PerpType {\n> +\t\tNone, /* no perpendicular found */\n> +\t\tStart, /* start of Pwl is closest point */\n> +\t\tEnd, /* end of Pwl is closest point */\n> +\t\tVertex, /* vertex of Pwl is closest point */\n> +\t\tPerpendicular /* true perpendicular found */\n> +\t};\n> +\tPerpType invert(Point const &xy, Point &perp, int &span,\n> +\t\t\tconst double eps = 1e-6) const;\n> +\t/*\n> +\t * Compute the inverse function. Indicate if it is a proper (true)\n> +\t * inverse, or only a best effort (e.g. input was non-monotonic).\n> +\t */\n> +\tPwl inverse(bool *trueInverse = nullptr, const double eps = 1e-6) const;\n> +\t/* Compose two Pwls together, doing \"this\" first and \"other\" after. */\n> +\tPwl compose(Pwl const &other, const double eps = 1e-6) const;\n> +\t/* Apply function to (x,y) values at every control point. */\n> +\tvoid map(std::function f) const;\n> +\t/*\n> +\t * Apply function to (x, y0, y1) values wherever either Pwl has a\n> +\t * control point.\n> +\t */\n> +\tstatic void map2(Pwl const &pwl0, Pwl const &pwl1,\n> +\t\t\t std::function f);\n> +\t/*\n> +\t * Combine two Pwls, meaning we create a new Pwl where the y values are\n> +\t * given by running f wherever either has a knot.\n> +\t */\n> +\tstatic Pwl\n> +\tcombine(Pwl const &pwl0, Pwl const &pwl1,\n> +\t\tstd::function f,\n> +\t\tconst double eps = 1e-6);\n> +\t/*\n> +\t * Make \"this\" match (at least) the given domain. Any extension my be\n> +\t * clipped or linear.\n> +\t */\n> +\tvoid matchDomain(Interval const &domain, bool clip = true,\n> +\t\t\t const double eps = 1e-6);\n> +\tPwl &operator*=(double d);\n> +\tvoid debug(FILE *fp = stdout) const;\n> +\n> +private:\n> +\tint findSpan(double x, int span) const;\n> +\tstd::vector points_;\n> +};\n> +\n> +} /* namespace RPiController */","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 36184BD87C\n\tfor ;\n\tTue, 11 Jun 2024 21:41:19 +0000 (UTC)","from lancelot.ideasonboard.com (localhost [IPv6:::1])\n\tby lancelot.ideasonboard.com (Postfix) with ESMTP id D863F65466;\n\tTue, 11 Jun 2024 23:41:18 +0200 (CEST)","from perceval.ideasonboard.com (perceval.ideasonboard.com\n\t[IPv6:2001:4b98:dc2:55:216:3eff:fef7:d647])\n\tby lancelot.ideasonboard.com (Postfix) with ESMTPS id D73C065446\n\tfor ;\n\tTue, 11 Jun 2024 23:41:17 +0200 (CEST)","from pendragon.ideasonboard.com (81-175-209-231.bb.dnainternet.fi\n\t[81.175.209.231])\n\tby perceval.ideasonboard.com (Postfix) with ESMTPSA id DB89B512;\n\tTue, 11 Jun 2024 23:41:04 +0200 (CEST)"],"Authentication-Results":"lancelot.ideasonboard.com; dkim=pass (1024-bit key;\n\tunprotected) header.d=ideasonboard.com header.i=@ideasonboard.com\n\theader.b=\"UP/aXydN\"; dkim-atps=neutral","DKIM-Signature":"v=1; a=rsa-sha256; c=relaxed/simple; d=ideasonboard.com;\n\ts=mail; t=1718142065;\n\tbh=C7FnYjXxfGpcq6/VXMoYSnmbhYfF2UOA+v1edy+8auU=;\n\th=Date:From:To:Cc:Subject:References:In-Reply-To:From;\n\tb=UP/aXydN95ksJ+AQQmgaN1N/pCSZBbGYZtHo/bTczJvAvn6stMv7lh5va1DZcWqYm\n\tpqN2KnY24lUoq4a2yeQZN+xiTSPGH2fEDhKJkBdGjJNMYuX/+dzymjX0bgukx2kdHw\n\tVFjf7SJHvvZgHr4jP95dpTUcpb8hXR7RkyZO9AGY=","Date":"Wed, 12 Jun 2024 00:40:58 +0300","From":"Laurent Pinchart ","To":"Paul Elder ","Cc":"libcamera-devel@lists.libcamera.org,\n\tStefan Klug ,\n\tDavid Plowman ,\n\tKieran Bingham ","Subject":"Re: [PATCH v8 2/4] ipa: libipa: Copy pwl from rpi","Message-ID":"<20240611214058.GC8290@pendragon.ideasonboard.com>","References":"<20240611132430.404814-1-paul.elder@ideasonboard.com>\n\t<20240611132430.404814-3-paul.elder@ideasonboard.com>","MIME-Version":"1.0","Content-Type":"text/plain; charset=utf-8","Content-Disposition":"inline","In-Reply-To":"<20240611132430.404814-3-paul.elder@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\" "}}]