[{"id":29746,"web_url":"https://patchwork.libcamera.org/comment/29746/","msgid":"<171741311095.205609.4614422647110576903@ping.linuxembedded.co.uk>","date":"2024-06-03T11:11:50","subject":"Re: [PATCH v4 2/4] ipa: libipa: Copy pwl from rpi","submitter":{"id":4,"url":"https://patchwork.libcamera.org/api/people/4/","name":"Kieran Bingham","email":"kieran.bingham@ideasonboard.com"},"content":"Quoting Paul Elder (2024-05-31 15:42:59)\n> Copy the piecewise linear function code from Raspberry Pi.\n> \n> Signed-off-by: Paul Elder <paul.elder@ideasonboard.com>\n> Reviewed-by: Stefan Klug <stefan.klug@ideasonboard.com>\n> Acked-by: David Plowman <david.plowman@raspberrypi.com>\n> \n> ---\n> Changes in v4:\n> - update the copy\n\nFiles are identical in my check so:\n\n\nReviewed-by: Kieran Bingham <kieran.bingham@ideasonboard.com>\n\n> \n> No change in v3\n> \n> No change in v2\n> ---\n>  src/ipa/libipa/meson.build |   2 +\n>  src/ipa/libipa/pwl.cpp     | 269 +++++++++++++++++++++++++++++++++++++\n>  src/ipa/libipa/pwl.h       | 127 +++++++++++++++++\n>  3 files changed, 398 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/meson.build b/src/ipa/libipa/meson.build\n> index 7ce885da7..8ec9c7847 100644\n> --- a/src/ipa/libipa/meson.build\n> +++ b/src/ipa/libipa/meson.build\n> @@ -8,6 +8,7 @@ libipa_headers = files([\n>      'fc_queue.h',\n>      'histogram.h',\n>      'module.h',\n> +    'pwl.h',\n>  ])\n>  \n>  libipa_sources = files([\n> @@ -18,6 +19,7 @@ libipa_sources = files([\n>      'fc_queue.cpp',\n>      'histogram.cpp',\n>      'module.cpp',\n> +    'pwl.cpp',\n>  ])\n>  \n>  libipa_includes = include_directories('..')\n> diff --git a/src/ipa/libipa/pwl.cpp b/src/ipa/libipa/pwl.cpp\n> new file mode 100644\n> index 000000000..e39123767\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 <cassert>\n> +#include <cmath>\n> +#include <stdexcept>\n> +\n> +#include \"pwl.h\"\n> +\n> +using namespace RPiController;\n> +\n> +int Pwl::read(const libcamera::YamlObject &params)\n> +{\n> +       if (!params.size() || params.size() % 2)\n> +               return -EINVAL;\n> +\n> +       const auto &list = params.asList();\n> +\n> +       for (auto it = list.begin(); it != list.end(); it++) {\n> +               auto x = it->get<double>();\n> +               if (!x)\n> +                       return -EINVAL;\n> +               if (it != list.begin() && *x <= points_.back().x)\n> +                       return -EINVAL;\n> +\n> +               auto y = (++it)->get<double>();\n> +               if (!y)\n> +                       return -EINVAL;\n> +\n> +               points_.push_back(Point(*x, *y));\n> +       }\n> +\n> +       return 0;\n> +}\n> +\n> +void Pwl::append(double x, double y, const double eps)\n> +{\n> +       if (points_.empty() || points_.back().x + eps < x)\n> +               points_.push_back(Point(x, y));\n> +}\n> +\n> +void Pwl::prepend(double x, double y, const double eps)\n> +{\n> +       if (points_.empty() || points_.front().x - eps > x)\n> +               points_.insert(points_.begin(), Point(x, y));\n> +}\n> +\n> +Pwl::Interval Pwl::domain() const\n> +{\n> +       return Interval(points_[0].x, points_[points_.size() - 1].x);\n> +}\n> +\n> +Pwl::Interval Pwl::range() const\n> +{\n> +       double lo = points_[0].y, hi = lo;\n> +       for (auto &p : points_)\n> +               lo = std::min(lo, p.y), hi = std::max(hi, p.y);\n> +       return Interval(lo, hi);\n> +}\n> +\n> +bool Pwl::empty() const\n> +{\n> +       return points_.empty();\n> +}\n> +\n> +double Pwl::eval(double x, int *spanPtr, bool updateSpan) const\n> +{\n> +       int span = findSpan(x, spanPtr && *spanPtr != -1 ? *spanPtr : points_.size() / 2 - 1);\n> +       if (spanPtr && updateSpan)\n> +               *spanPtr = span;\n> +       return points_[span].y +\n> +              (x - points_[span].x) * (points_[span + 1].y - points_[span].y) /\n> +                      (points_[span + 1].x - points_[span].x);\n> +}\n> +\n> +int Pwl::findSpan(double x, int span) const\n> +{\n> +       /*\n> +        * Pwls are generally small, so linear search may well be faster than\n> +        * binary, though could review this if large PWls start turning up.\n> +        */\n> +       int lastSpan = points_.size() - 2;\n> +       /*\n> +        * some algorithms may call us with span pointing directly at the last\n> +        * control point\n> +        */\n> +       span = std::max(0, std::min(lastSpan, span));\n> +       while (span < lastSpan && x >= points_[span + 1].x)\n> +               span++;\n> +       while (span && x < points_[span].x)\n> +               span--;\n> +       return span;\n> +}\n> +\n> +Pwl::PerpType Pwl::invert(Point const &xy, Point &perp, int &span,\n> +                         const double eps) const\n> +{\n> +       assert(span >= -1);\n> +       bool prevOffEnd = false;\n> +       for (span = span + 1; span < (int)points_.size() - 1; span++) {\n> +               Point spanVec = points_[span + 1] - points_[span];\n> +               double t = ((xy - points_[span]) % spanVec) / spanVec.len2();\n> +               if (t < -eps) /* off the start of this span */\n> +               {\n> +                       if (span == 0) {\n> +                               perp = points_[span];\n> +                               return PerpType::Start;\n> +                       } else if (prevOffEnd) {\n> +                               perp = points_[span];\n> +                               return PerpType::Vertex;\n> +                       }\n> +               } else if (t > 1 + eps) /* off the end of this span */\n> +               {\n> +                       if (span == (int)points_.size() - 2) {\n> +                               perp = points_[span + 1];\n> +                               return PerpType::End;\n> +                       }\n> +                       prevOffEnd = true;\n> +               } else /* a true perpendicular */\n> +               {\n> +                       perp = points_[span] + spanVec * t;\n> +                       return PerpType::Perpendicular;\n> +               }\n> +       }\n> +       return PerpType::None;\n> +}\n> +\n> +Pwl Pwl::inverse(bool *trueInverse, const double eps) const\n> +{\n> +       bool appended = false, prepended = false, neither = false;\n> +       Pwl inverse;\n> +\n> +       for (Point const &p : points_) {\n> +               if (inverse.empty())\n> +                       inverse.append(p.y, p.x, eps);\n> +               else if (std::abs(inverse.points_.back().x - p.y) <= eps ||\n> +                        std::abs(inverse.points_.front().x - p.y) <= eps)\n> +                       /* do nothing */;\n> +               else if (p.y > inverse.points_.back().x) {\n> +                       inverse.append(p.y, p.x, eps);\n> +                       appended = true;\n> +               } else if (p.y < inverse.points_.front().x) {\n> +                       inverse.prepend(p.y, p.x, eps);\n> +                       prepended = true;\n> +               } else\n> +                       neither = true;\n> +       }\n> +\n> +       /*\n> +        * This is not a proper inverse if we found ourselves putting points\n> +        * onto both ends of the inverse, or if there were points that couldn't\n> +        * go on either.\n> +        */\n> +       if (trueInverse)\n> +               *trueInverse = !(neither || (appended && prepended));\n> +\n> +       return inverse;\n> +}\n> +\n> +Pwl Pwl::compose(Pwl const &other, const double eps) const\n> +{\n> +       double thisX = points_[0].x, thisY = points_[0].y;\n> +       int thisSpan = 0, otherSpan = other.findSpan(thisY, 0);\n> +       Pwl result({ { thisX, other.eval(thisY, &otherSpan, false) } });\n> +       while (thisSpan != (int)points_.size() - 1) {\n> +               double dx = points_[thisSpan + 1].x - points_[thisSpan].x,\n> +                      dy = points_[thisSpan + 1].y - points_[thisSpan].y;\n> +               if (std::abs(dy) > eps &&\n> +                   otherSpan + 1 < (int)other.points_.size() &&\n> +                   points_[thisSpan + 1].y >=\n> +                           other.points_[otherSpan + 1].x + eps) {\n> +                       /*\n> +                        * next control point in result will be where this\n> +                        * function's y reaches the next span in other\n> +                        */\n> +                       thisX = points_[thisSpan].x +\n> +                               (other.points_[otherSpan + 1].x -\n> +                                points_[thisSpan].y) *\n> +                                       dx / dy;\n> +                       thisY = other.points_[++otherSpan].x;\n> +               } else if (std::abs(dy) > eps && otherSpan > 0 &&\n> +                          points_[thisSpan + 1].y <=\n> +                                  other.points_[otherSpan - 1].x - eps) {\n> +                       /*\n> +                        * next control point in result will be where this\n> +                        * function's y reaches the previous span in other\n> +                        */\n> +                       thisX = points_[thisSpan].x +\n> +                               (other.points_[otherSpan + 1].x -\n> +                                points_[thisSpan].y) *\n> +                                       dx / dy;\n> +                       thisY = other.points_[--otherSpan].x;\n> +               } else {\n> +                       /* we stay in the same span in other */\n> +                       thisSpan++;\n> +                       thisX = points_[thisSpan].x,\n> +                       thisY = points_[thisSpan].y;\n> +               }\n> +               result.append(thisX, other.eval(thisY, &otherSpan, false),\n> +                             eps);\n> +       }\n> +       return result;\n> +}\n> +\n> +void Pwl::map(std::function<void(double x, double y)> f) const\n> +{\n> +       for (auto &pt : points_)\n> +               f(pt.x, pt.y);\n> +}\n> +\n> +void Pwl::map2(Pwl const &pwl0, Pwl const &pwl1,\n> +              std::function<void(double x, double y0, double y1)> f)\n> +{\n> +       int span0 = 0, span1 = 0;\n> +       double x = std::min(pwl0.points_[0].x, pwl1.points_[0].x);\n> +       f(x, pwl0.eval(x, &span0, false), pwl1.eval(x, &span1, false));\n> +       while (span0 < (int)pwl0.points_.size() - 1 ||\n> +              span1 < (int)pwl1.points_.size() - 1) {\n> +               if (span0 == (int)pwl0.points_.size() - 1)\n> +                       x = pwl1.points_[++span1].x;\n> +               else if (span1 == (int)pwl1.points_.size() - 1)\n> +                       x = pwl0.points_[++span0].x;\n> +               else if (pwl0.points_[span0 + 1].x > pwl1.points_[span1 + 1].x)\n> +                       x = pwl1.points_[++span1].x;\n> +               else\n> +                       x = pwl0.points_[++span0].x;\n> +               f(x, pwl0.eval(x, &span0, false), pwl1.eval(x, &span1, false));\n> +       }\n> +}\n> +\n> +Pwl Pwl::combine(Pwl const &pwl0, Pwl const &pwl1,\n> +                std::function<double(double x, double y0, double y1)> f,\n> +                const double eps)\n> +{\n> +       Pwl result;\n> +       map2(pwl0, pwl1, [&](double x, double y0, double y1) {\n> +               result.append(x, f(x, y0, y1), eps);\n> +       });\n> +       return result;\n> +}\n> +\n> +void Pwl::matchDomain(Interval const &domain, bool clip, const double eps)\n> +{\n> +       int span = 0;\n> +       prepend(domain.start, eval(clip ? points_[0].x : domain.start, &span),\n> +               eps);\n> +       span = points_.size() - 2;\n> +       append(domain.end, eval(clip ? points_.back().x : domain.end, &span),\n> +              eps);\n> +}\n> +\n> +Pwl &Pwl::operator*=(double d)\n> +{\n> +       for (auto &pt : points_)\n> +               pt.y *= d;\n> +       return *this;\n> +}\n> +\n> +void Pwl::debug(FILE *fp) const\n> +{\n> +       fprintf(fp, \"Pwl {\\n\");\n> +       for (auto &p : points_)\n> +               fprintf(fp, \"\\t(%g, %g)\\n\", p.x, p.y);\n> +       fprintf(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 000000000..7d5e7e4d3\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 <functional>\n> +#include <math.h>\n> +#include <vector>\n> +\n> +#include \"libcamera/internal/yaml_parser.h\"\n> +\n> +namespace RPiController {\n> +\n> +class Pwl\n> +{\n> +public:\n> +       struct Interval {\n> +               Interval(double _start, double _end)\n> +                       : start(_start), end(_end)\n> +               {\n> +               }\n> +               double start, end;\n> +               bool contains(double value)\n> +               {\n> +                       return value >= start && value <= end;\n> +               }\n> +               double clip(double value)\n> +               {\n> +                       return value < start ? start\n> +                                            : (value > end ? end : value);\n> +               }\n> +               double len() const { return end - start; }\n> +       };\n> +       struct Point {\n> +               Point() : x(0), y(0) {}\n> +               Point(double _x, double _y)\n> +                       : x(_x), y(_y) {}\n> +               double x, y;\n> +               Point operator-(Point const &p) const\n> +               {\n> +                       return Point(x - p.x, y - p.y);\n> +               }\n> +               Point operator+(Point const &p) const\n> +               {\n> +                       return Point(x + p.x, y + p.y);\n> +               }\n> +               double operator%(Point const &p) const\n> +               {\n> +                       return x * p.x + y * p.y;\n> +               }\n> +               Point operator*(double f) const { return Point(x * f, y * f); }\n> +               Point operator/(double f) const { return Point(x / f, y / f); }\n> +               double len2() const { return x * x + y * y; }\n> +               double len() const { return sqrt(len2()); }\n> +       };\n> +       Pwl() {}\n> +       Pwl(std::vector<Point> const &points) : points_(points) {}\n> +       int read(const libcamera::YamlObject &params);\n> +       void append(double x, double y, const double eps = 1e-6);\n> +       void prepend(double x, double y, const double eps = 1e-6);\n> +       Interval domain() const;\n> +       Interval range() const;\n> +       bool empty() const;\n> +       /*\n> +        * Evaluate Pwl, optionally supplying an initial guess for the\n> +        * \"span\". The \"span\" may be optionally be updated.  If you want to know\n> +        * the \"span\" value but don't have an initial guess you can set it to\n> +        * -1.\n> +        */\n> +       double eval(double x, int *spanPtr = nullptr,\n> +                   bool updateSpan = true) const;\n> +       /*\n> +        * Find perpendicular closest to xy, starting from span+1 so you can\n> +        * call it repeatedly to check for multiple closest points (set span to\n> +        * -1 on the first call). Also returns \"pseudo\" perpendiculars; see\n> +        * PerpType enum.\n> +        */\n> +       enum class PerpType {\n> +               None, /* no perpendicular found */\n> +               Start, /* start of Pwl is closest point */\n> +               End, /* end of Pwl is closest point */\n> +               Vertex, /* vertex of Pwl is closest point */\n> +               Perpendicular /* true perpendicular found */\n> +       };\n> +       PerpType invert(Point const &xy, Point &perp, int &span,\n> +                       const double eps = 1e-6) const;\n> +       /*\n> +        * Compute the inverse function. Indicate if it is a proper (true)\n> +        * inverse, or only a best effort (e.g. input was non-monotonic).\n> +        */\n> +       Pwl inverse(bool *trueInverse = nullptr, const double eps = 1e-6) const;\n> +       /* Compose two Pwls together, doing \"this\" first and \"other\" after. */\n> +       Pwl compose(Pwl const &other, const double eps = 1e-6) const;\n> +       /* Apply function to (x,y) values at every control point. */\n> +       void map(std::function<void(double x, double y)> f) const;\n> +       /*\n> +        * Apply function to (x, y0, y1) values wherever either Pwl has a\n> +        * control point.\n> +        */\n> +       static void map2(Pwl const &pwl0, Pwl const &pwl1,\n> +                        std::function<void(double x, double y0, double y1)> f);\n> +       /*\n> +        * Combine two Pwls, meaning we create a new Pwl where the y values are\n> +        * given by running f wherever either has a knot.\n> +        */\n> +       static Pwl\n> +       combine(Pwl const &pwl0, Pwl const &pwl1,\n> +               std::function<double(double x, double y0, double y1)> f,\n> +               const double eps = 1e-6);\n> +       /*\n> +        * Make \"this\" match (at least) the given domain. Any extension my be\n> +        * clipped or linear.\n> +        */\n> +       void matchDomain(Interval const &domain, bool clip = true,\n> +                        const double eps = 1e-6);\n> +       Pwl &operator*=(double d);\n> +       void debug(FILE *fp = stdout) const;\n> +\n> +private:\n> +       int findSpan(double x, int span) const;\n> +       std::vector<Point> points_;\n> +};\n> +\n> +} /* namespace RPiController */\n> -- \n> 2.39.2\n>","headers":{"Return-Path":"<libcamera-devel-bounces@lists.libcamera.org>","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 77392BD87C\n\tfor <parsemail@patchwork.libcamera.org>;\n\tMon,  3 Jun 2024 11:11:56 +0000 (UTC)","from lancelot.ideasonboard.com (localhost [IPv6:::1])\n\tby lancelot.ideasonboard.com (Postfix) with ESMTP id 7E8FE634C5;\n\tMon,  3 Jun 2024 13:11:55 +0200 (CEST)","from perceval.ideasonboard.com (perceval.ideasonboard.com\n\t[213.167.242.64])\n\tby lancelot.ideasonboard.com (Postfix) with ESMTPS id F01F861A3B\n\tfor <libcamera-devel@lists.libcamera.org>;\n\tMon,  3 Jun 2024 13:11:53 +0200 (CEST)","from pendragon.ideasonboard.com\n\t(cpc89244-aztw30-2-0-cust6594.18-1.cable.virginm.net [86.31.185.195])\n\tby perceval.ideasonboard.com (Postfix) with ESMTPSA id CD04F908;\n\tMon,  3 Jun 2024 13:11:46 +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=\"ER7Qh+zW\"; dkim-atps=neutral","DKIM-Signature":"v=1; a=rsa-sha256; c=relaxed/simple; d=ideasonboard.com;\n\ts=mail; t=1717413106;\n\tbh=KlgO7gZvNIKiBpZRurWTtOosaTviR31R+lIbJ+b2/yc=;\n\th=In-Reply-To:References:Subject:From:Cc:To:Date:From;\n\tb=ER7Qh+zWobdulUT0BOkadwSJOAUdcd/Z2Rczt8Kz86MyBaMJBaxn6wf8rKzLPDbOF\n\tb7QQs2QLFU4QRQdsFEsl/a0YI+Hp9XOqty672TwkzJdQJpxJDYKT29NVRelJuIvLhC\n\tjPHEZVF82aSTFyYzoBiEUK9xPm6bw7BjeEdnp0nA=","Content-Type":"text/plain; charset=\"utf-8\"","MIME-Version":"1.0","Content-Transfer-Encoding":"quoted-printable","In-Reply-To":"<20240531144301.3950115-3-paul.elder@ideasonboard.com>","References":"<20240531144301.3950115-1-paul.elder@ideasonboard.com>\n\t<20240531144301.3950115-3-paul.elder@ideasonboard.com>","Subject":"Re: [PATCH v4 2/4] ipa: libipa: Copy pwl from rpi","From":"Kieran Bingham <kieran.bingham@ideasonboard.com>","Cc":"Paul Elder <paul.elder@ideasonboard.com>,\n\tStefan Klug <stefan.klug@ideasonboard.com>,\n\tDavid Plowman <david.plowman@raspberrypi.com>","To":"Paul Elder <paul.elder@ideasonboard.com>,\n\tlibcamera-devel@lists.libcamera.org","Date":"Mon, 03 Jun 2024 12:11:50 +0100","Message-ID":"<171741311095.205609.4614422647110576903@ping.linuxembedded.co.uk>","User-Agent":"alot/0.10","X-BeenThere":"libcamera-devel@lists.libcamera.org","X-Mailman-Version":"2.1.29","Precedence":"list","List-Id":"<libcamera-devel.lists.libcamera.org>","List-Unsubscribe":"<https://lists.libcamera.org/options/libcamera-devel>,\n\t<mailto:libcamera-devel-request@lists.libcamera.org?subject=unsubscribe>","List-Archive":"<https://lists.libcamera.org/pipermail/libcamera-devel/>","List-Post":"<mailto:libcamera-devel@lists.libcamera.org>","List-Help":"<mailto:libcamera-devel-request@lists.libcamera.org?subject=help>","List-Subscribe":"<https://lists.libcamera.org/listinfo/libcamera-devel>,\n\t<mailto:libcamera-devel-request@lists.libcamera.org?subject=subscribe>","Errors-To":"libcamera-devel-bounces@lists.libcamera.org","Sender":"\"libcamera-devel\" <libcamera-devel-bounces@lists.libcamera.org>"}}]