Patch Detail
Show a patch.
GET /api/patches/23121/?format=api
{ "id": 23121, "url": "https://patchwork.libcamera.org/api/patches/23121/?format=api", "web_url": "https://patchwork.libcamera.org/patch/23121/", "project": { "id": 1, "url": "https://patchwork.libcamera.org/api/projects/1/?format=api", "name": "libcamera", "link_name": "libcamera", "list_id": "libcamera_core", "list_email": "libcamera-devel@lists.libcamera.org", "web_url": "", "scm_url": "", "webscm_url": "" }, "msgid": "<20250403154925.382973-6-stefan.klug@ideasonboard.com>", "date": "2025-04-03T15:49:10", "name": "[v3,05/16] libcamera: matrix: Add inverse() function", "commit_ref": null, "pull_url": null, "state": "accepted", "archived": false, "hash": "6d9f659271c521495095642799a1936467fe5df8", "submitter": { "id": 184, "url": "https://patchwork.libcamera.org/api/people/184/?format=api", "name": "Stefan Klug", "email": "stefan.klug@ideasonboard.com" }, "delegate": null, "mbox": "https://patchwork.libcamera.org/patch/23121/mbox/", "series": [ { "id": 5111, "url": "https://patchwork.libcamera.org/api/series/5111/?format=api", "web_url": "https://patchwork.libcamera.org/project/libcamera/list/?series=5111", "date": "2025-04-03T15:49:05", "name": "Some rkisp1 awb improvements", "version": 3, "mbox": "https://patchwork.libcamera.org/series/5111/mbox/" } ], "comments": "https://patchwork.libcamera.org/api/patches/23121/comments/", "check": "pending", "checks": "https://patchwork.libcamera.org/api/patches/23121/checks/", "tags": {}, "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 04AC2C327D\n\tfor <parsemail@patchwork.libcamera.org>;\n\tThu, 3 Apr 2025 15:49:52 +0000 (UTC)", "from lancelot.ideasonboard.com (localhost [IPv6:::1])\n\tby lancelot.ideasonboard.com (Postfix) with ESMTP id A271B689D2;\n\tThu, 3 Apr 2025 17:49:51 +0200 (CEST)", "from perceval.ideasonboard.com (perceval.ideasonboard.com\n\t[213.167.242.64])\n\tby lancelot.ideasonboard.com (Postfix) with ESMTPS id 1345E689AD\n\tfor <libcamera-devel@lists.libcamera.org>;\n\tThu, 3 Apr 2025 17:49:50 +0200 (CEST)", "from ideasonboard.com (unknown\n\t[IPv6:2a00:6020:448c:6c00:6d9d:9854:3fc1:4bb2])\n\tby perceval.ideasonboard.com (Postfix) with ESMTPSA id 62BBD11E9;\n\tThu, 3 Apr 2025 17:47:56 +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=\"mDLO6sUM\"; dkim-atps=neutral", "DKIM-Signature": "v=1; a=rsa-sha256; c=relaxed/simple; d=ideasonboard.com;\n\ts=mail; t=1743695276;\n\tbh=4F4zExUEB+OlZERtsESklQcZSlANX3zcMQjF6KU8TpY=;\n\th=From:To:Cc:Subject:Date:In-Reply-To:References:From;\n\tb=mDLO6sUM5sWaNFQ53UE2WWy+r7V8wyyKvZjnY7+/bh98BGX/BaApieWjtc4JgQzpB\n\tVDm4u1NkJ1oVSPyCyPVJ2S7cdS6k2C41tHmDEIYpNwlUJ2hsGHgUsrIofUUc9loyJ9\n\tsdfVOoGzvrDviQA6xF3uhl9jCym2zQBQKRJO7qHY=", "From": "Stefan Klug <stefan.klug@ideasonboard.com>", "To": "libcamera-devel@lists.libcamera.org", "Cc": "Stefan Klug <stefan.klug@ideasonboard.com>,\n\tLaurent Pinchart <laurent.pinchart@ideasonboard.com>", "Subject": "[PATCH v3 05/16] libcamera: matrix: Add inverse() function", "Date": "Thu, 3 Apr 2025 17:49:10 +0200", "Message-ID": "<20250403154925.382973-6-stefan.klug@ideasonboard.com>", "X-Mailer": "git-send-email 2.43.0", "In-Reply-To": "<20250403154925.382973-1-stefan.klug@ideasonboard.com>", "References": "<20250403154925.382973-1-stefan.klug@ideasonboard.com>", "MIME-Version": "1.0", "Content-Transfer-Encoding": "8bit", "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>" }, "content": "For calculations in upcoming algorithm patches, the inverse of a matrix\nis required. Add an implementation of the inverse() function for square\nmatrices.\n\nSigned-off-by: Stefan Klug <stefan.klug@ideasonboard.com>\nSigned-off-by: Laurent Pinchart <laurent.pinchart@ideasonboard.com>\n\n---\n\nChanges in v2:\n- Replaced the implementation by a generic one provided by Laurent that\n supports arbitrary square matrices instead of 2x2 and 3x3 only.\n- Moved the implementation into the cpp file.\n\nChanges inv3:\n- Added stack allocated scratchBuffers to matrixInvert so that no heap\n allocations are necessary.\n---\n include/libcamera/internal/matrix.h | 23 ++++\n src/libcamera/matrix.cpp | 166 ++++++++++++++++++++++++++++\n 2 files changed, 189 insertions(+)", "diff": "diff --git a/include/libcamera/internal/matrix.h b/include/libcamera/internal/matrix.h\nindex 6d40567af0a0..a07a47701336 100644\n--- a/include/libcamera/internal/matrix.h\n+++ b/include/libcamera/internal/matrix.h\n@@ -19,6 +19,12 @@ namespace libcamera {\n \n LOG_DECLARE_CATEGORY(Matrix)\n \n+#ifndef __DOXYGEN__\n+template<typename T>\n+bool matrixInvert(Span<const T> dataIn, Span<T> dataOut, unsigned int dim,\n+\t\t Span<T> scratchBuffer, Span<unsigned int> swapBuffer);\n+#endif /* __DOXYGEN__ */\n+\n template<typename T, unsigned int Rows, unsigned int Cols>\n class Matrix\n {\n@@ -91,6 +97,23 @@ public:\n \t\treturn *this;\n \t}\n \n+\tMatrix<T, Rows, Cols> inverse(bool *ok = nullptr) const\n+\t{\n+\t\tstatic_assert(Rows == Cols, \"Matrix must be square\");\n+\n+\t\tMatrix<T, Rows, Cols> inverse;\n+\t\tstd::array<T, Rows * Cols * 2> scratchBuffer;\n+\t\tstd::array<unsigned int, Rows> swapBuffer;\n+\t\tbool res = matrixInvert(Span<const T>(data_),\n+\t\t\t\t\tSpan<T>(inverse.data_),\n+\t\t\t\t\tRows,\n+\t\t\t\t\tSpan<T>(scratchBuffer),\n+\t\t\t\t\tSpan<unsigned int>(swapBuffer));\n+\t\tif (ok)\n+\t\t\t*ok = res;\n+\t\treturn inverse;\n+\t}\n+\n private:\n \t/*\n \t * \\todo The initializer is only necessary for the constructor to be\ndiff --git a/src/libcamera/matrix.cpp b/src/libcamera/matrix.cpp\nindex 49e2aa3b4f2c..68fc1b7bd5ac 100644\n--- a/src/libcamera/matrix.cpp\n+++ b/src/libcamera/matrix.cpp\n@@ -7,6 +7,12 @@\n \n #include \"libcamera/internal/matrix.h\"\n \n+#include <algorithm>\n+#include <assert.h>\n+#include <cmath>\n+#include <numeric>\n+#include <vector>\n+\n #include <libcamera/base/log.h>\n \n /**\n@@ -87,6 +93,20 @@ LOG_DEFINE_CATEGORY(Matrix)\n * \\return Row \\a i from the matrix, as a Span\n */\n \n+/**\n+ * \\fn Matrix::inverse(bool *ok) const\n+ * \\param[out] ok Indicate if the matrix was successfully inverted\n+ * \\brief Compute the inverse of the matrix\n+ *\n+ * This function computes the inverse of the matrix. It is only implemented for\n+ * matrices of float and double types. If \\a ok is provided it will be set to a\n+ * boolean value to indicate of the inversion was successful. This can be used\n+ * to check if the matrix is singular, in which case the function will return\n+ * an identity matrix.\n+ *\n+ * \\return The inverse of the matrix\n+ */\n+\n /**\n * \\fn Matrix::operator[](size_t i)\n * \\copydoc Matrix::operator[](size_t i) const\n@@ -142,6 +162,152 @@ LOG_DEFINE_CATEGORY(Matrix)\n */\n \n #ifndef __DOXYGEN__\n+template<typename T>\n+bool matrixInvert(Span<const T> dataIn, Span<T> dataOut, unsigned int dim,\n+\t\t Span<T> scratchBuffer, Span<unsigned int> swapBuffer)\n+{\n+\t/*\n+\t * Convenience class to access matrix data, providing a row-major (i,j)\n+\t * element accessor through the call operator, and the ability to swap\n+\t * rows without modifying the backing storage.\n+\t */\n+\tclass MatrixAccessor\n+\t{\n+\tpublic:\n+\t\tMatrixAccessor(Span<T> data, Span<unsigned int> swapBuffer, unsigned int rows, unsigned int cols)\n+\t\t\t: data_(data), swap_(swapBuffer), rows_(rows), cols_(cols)\n+\t\t{\n+\t\t\tASSERT(swap_.size() == rows);\n+\t\t\tstd::iota(swap_.begin(), swap_.end(), T{ 0 });\n+\t\t}\n+\n+\t\tT &operator()(unsigned int row, unsigned int col)\n+\t\t{\n+\t\t\tassert(row < rows_ && col < cols_);\n+\t\t\treturn data_[index(row, col)];\n+\t\t}\n+\n+\t\tvoid swap(unsigned int a, unsigned int b)\n+\t\t{\n+\t\t\tassert(a < rows_ && a < cols_);\n+\t\t\tstd::swap(swap_[a], swap_[b]);\n+\t\t}\n+\n+\tprivate:\n+\t\tunsigned int index(unsigned int row, unsigned int col) const\n+\t\t{\n+\t\t\treturn swap_[row] * cols_ + col;\n+\t\t}\n+\n+\t\tSpan<T> data_;\n+\t\tSpan<unsigned int> swap_;\n+\t\tunsigned int rows_;\n+\t\tunsigned int cols_;\n+\t};\n+\n+\t/*\n+\t * Matrix inversion using Gaussian elimination.\n+\t *\n+\t * Start by augmenting the original matrix with an identiy matrix of\n+\t * the same size.\n+\t */\n+\tASSERT(scratchBuffer.size() == dim * dim * 2);\n+\tMatrixAccessor matrix(scratchBuffer, swapBuffer, dim, dim * 2);\n+\n+\tfor (unsigned int i = 0; i < dim; ++i) {\n+\t\tfor (unsigned int j = 0; j < dim; ++j) {\n+\t\t\tmatrix(i, j) = dataIn[i * dim + j];\n+\t\t\tmatrix(i, j + dim) = T{ 0 };\n+\t\t}\n+\t\tmatrix(i, i + dim) = T{ 1 };\n+\t}\n+\n+\t/* Start by triangularizing the input . */\n+\tfor (unsigned int pivot = 0; pivot < dim; ++pivot) {\n+\t\t/*\n+\t\t * Locate the next pivot. To improve numerical stability, use\n+\t\t * the row with the largest value in the pivot's column.\n+\t\t */\n+\t\tunsigned int row;\n+\t\tT maxValue{ 0 };\n+\n+\t\tfor (unsigned int i = pivot; i < dim; ++i) {\n+\t\t\tT value = std::abs(matrix(i, pivot));\n+\t\t\tif (maxValue < value) {\n+\t\t\t\tmaxValue = value;\n+\t\t\t\trow = i;\n+\t\t\t}\n+\t\t}\n+\n+\t\t/*\n+\t\t * If no pivot is found in the column, the matrix is not\n+\t\t * invertible. Return an identity matrix.\n+\t\t */\n+\t\tif (maxValue == 0) {\n+\t\t\tstd::fill(dataOut.begin(), dataOut.end(), T{ 0 });\n+\t\t\tfor (unsigned int i = 0; i < dim; ++i)\n+\t\t\t\tdataOut[i * dim + i] = T{ 1 };\n+\t\t\treturn false;\n+\t\t}\n+\n+\t\t/* Swap rows to bring the pivot in the right location. */\n+\t\tmatrix.swap(pivot, row);\n+\n+\t\t/* Process all rows below the pivot to zero the pivot column. */\n+\t\tconst T pivotValue = matrix(pivot, pivot);\n+\n+\t\tfor (unsigned int i = pivot + 1; i < dim; ++i) {\n+\t\t\tconst T factor = matrix(i, pivot) / pivotValue;\n+\n+\t\t\t/*\n+\t\t\t * We know the element in the pivot column will be 0,\n+\t\t\t * hardcode it instead of computing it.\n+\t\t\t */\n+\t\t\tmatrix(i, pivot) = T{ 0 };\n+\n+\t\t\tfor (unsigned int j = pivot + 1; j < dim * 2; ++j)\n+\t\t\t\tmatrix(i, j) -= matrix(pivot, j) * factor;\n+\t\t}\n+\t}\n+\n+\t/*\n+\t * Then diagonalize the input, walking the diagonal backwards. There's\n+\t * no need to update the input matrix, as all the values we would write\n+\t * in the top-right triangle aren't used in further calculations (and\n+\t * would all by definition be zero).\n+\t */\n+\tfor (unsigned int pivot = dim - 1; pivot > 0; --pivot) {\n+\t\tconst T pivotValue = matrix(pivot, pivot);\n+\n+\t\tfor (unsigned int i = 0; i < pivot; ++i) {\n+\t\t\tconst T factor = matrix(i, pivot) / pivotValue;\n+\n+\t\t\tfor (unsigned int j = dim; j < dim * 2; ++j)\n+\t\t\t\tmatrix(i, j) -= matrix(pivot, j) * factor;\n+\t\t}\n+\t}\n+\n+\t/*\n+\t * Finally, normalize the diagonal and store the result in the output\n+\t * data.\n+\t */\n+\tfor (unsigned int i = 0; i < dim; ++i) {\n+\t\tconst T factor = matrix(i, i);\n+\n+\t\tfor (unsigned int j = 0; j < dim; ++j)\n+\t\t\tdataOut[i * dim + j] = matrix(i, j + dim) / factor;\n+\t}\n+\n+\treturn true;\n+}\n+\n+template bool matrixInvert<float>(Span<const float> dataIn, Span<float> dataOut,\n+\t\t\t\t unsigned int dim, Span<float> scratchBuffer,\n+\t\t\t\t Span<unsigned int> swapBuffer);\n+template bool matrixInvert<double>(Span<const double> data, Span<double> dataOut,\n+\t\t\t\t unsigned int dim, Span<double> scratchBuffer,\n+\t\t\t\t Span<unsigned int> swapBuffer);\n+\n /*\n * The YAML data shall be a list of numerical values. Its size shall be equal\n * to the product of the number of rows and columns of the matrix (Rows x\n", "prefixes": [ "v3", "05/16" ] }