| Message ID | 20250217100203.297894-4-stefan.klug@ideasonboard.com |
|---|---|
| State | New |
| Headers | show |
| Series |
|
| Related | show |
Hi Stefan, Thank you for the patch. On Mon, Feb 17, 2025 at 11:01:44AM +0100, Stefan Klug wrote: > For calculations in upcoming algorithm patches, the inverse of a matrix > is required. Add an implementation of the inverse() function for 3x3 > matrices. > > Signed-off-by: Stefan Klug <stefan.klug@ideasonboard.com> > --- > include/libcamera/internal/matrix.h | 24 ++++++++++++++++++++++++ > src/libcamera/matrix.cpp | 10 ++++++++++ > 2 files changed, 34 insertions(+) > > diff --git a/include/libcamera/internal/matrix.h b/include/libcamera/internal/matrix.h > index ca81dcda93c4..053336e9cfa4 100644 > --- a/include/libcamera/internal/matrix.h > +++ b/include/libcamera/internal/matrix.h > @@ -90,6 +90,30 @@ public: > return *this; > } > > + Matrix<T, Rows, Cols> inverse() const > + { > + static_assert(Rows == 3 && Cols == 3, "Matrix must be 3x3"); You know what I'll ask, right ? :-) Inverting the matrix only makes sense when T is a floating point type. I'd add a static assertion or enable_if for this. > + > + const auto &m = *this; > + double det = m[0][0] * (m[1][1] * m[2][2] - m[2][1] * m[1][2]) - > + m[0][1] * (m[1][0] * m[2][2] - m[1][2] * m[2][0]) + > + m[0][2] * (m[1][0] * m[2][1] - m[1][1] * m[2][0]); Calculating the determinant should be moved to its own function, it's useful by itself. > + > + double invdet = 1 / det; What if det is zero (or very close) ? > + > + Matrix<T, Rows, Cols> ret; > + ret[0][0] = (m[1][1] * m[2][2] - m[2][1] * m[1][2]) * invdet; > + ret[0][1] = (m[0][2] * m[2][1] - m[0][1] * m[2][2]) * invdet; > + ret[0][2] = (m[0][1] * m[1][2] - m[0][2] * m[1][1]) * invdet; > + ret[1][0] = (m[1][2] * m[2][0] - m[1][0] * m[2][2]) * invdet; > + ret[1][1] = (m[0][0] * m[2][2] - m[0][2] * m[2][0]) * invdet; > + ret[1][2] = (m[1][0] * m[0][2] - m[0][0] * m[1][2]) * invdet; > + ret[2][0] = (m[1][0] * m[2][1] - m[2][0] * m[1][1]) * invdet; > + ret[2][1] = (m[2][0] * m[0][1] - m[0][0] * m[2][1]) * invdet; > + ret[2][2] = (m[0][0] * m[1][1] - m[1][0] * m[0][1]) * invdet; That's a lot of calculation for an inline function, please move it to the .cpp file. > + return ret; > + } > + > template<typename T2> > Matrix<T2, Rows, Cols> cast() const > { > diff --git a/src/libcamera/matrix.cpp b/src/libcamera/matrix.cpp > index 91a3f16405a3..b17a1938f1ba 100644 > --- a/src/libcamera/matrix.cpp > +++ b/src/libcamera/matrix.cpp > @@ -77,6 +77,16 @@ LOG_DEFINE_CATEGORY(Matrix) > * \return Row \a i from the matrix, as a Span > */ > > +/** > + * \fn Matrix::inverse() const > + * \brief Compute the inverse of the matrix > + * > + * This function computes the inverse of the matrix. It is only implemented for > + * 3x3 matrices. > + * > + * \return The inverse of the matrix > + */ > + > /** > * \fn template<typename T2> Matrix<T2, Rows, Cols> Matrix::cast() const > * \brief Cast the matrix to a different type
Hi Laurent, Thank you for the review. On Mon, Feb 17, 2025 at 12:50:24PM +0200, Laurent Pinchart wrote: > Hi Stefan, > > Thank you for the patch. > > On Mon, Feb 17, 2025 at 11:01:44AM +0100, Stefan Klug wrote: > > For calculations in upcoming algorithm patches, the inverse of a matrix > > is required. Add an implementation of the inverse() function for 3x3 > > matrices. > > > > Signed-off-by: Stefan Klug <stefan.klug@ideasonboard.com> > > --- > > include/libcamera/internal/matrix.h | 24 ++++++++++++++++++++++++ > > src/libcamera/matrix.cpp | 10 ++++++++++ > > 2 files changed, 34 insertions(+) > > > > diff --git a/include/libcamera/internal/matrix.h b/include/libcamera/internal/matrix.h > > index ca81dcda93c4..053336e9cfa4 100644 > > --- a/include/libcamera/internal/matrix.h > > +++ b/include/libcamera/internal/matrix.h > > @@ -90,6 +90,30 @@ public: > > return *this; > > } > > > > + Matrix<T, Rows, Cols> inverse() const > > + { > > + static_assert(Rows == 3 && Cols == 3, "Matrix must be 3x3"); > > You know what I'll ask, right ? :-) > > Inverting the matrix only makes sense when T is a floating point type. > I'd add a static assertion or enable_if for this. Yes, I'll add a static assert for that. Even though we then remove the theoretical ability to use that function with a custom defined type. But as that is most likely never going to happen, I can restrict it to double/float... > > > + > > + const auto &m = *this; > > + double det = m[0][0] * (m[1][1] * m[2][2] - m[2][1] * m[1][2]) - > > + m[0][1] * (m[1][0] * m[2][2] - m[1][2] * m[2][0]) + > > + m[0][2] * (m[1][0] * m[2][1] - m[1][1] * m[2][0]); > > Calculating the determinant should be moved to its own function, it's > useful by itself. Yes, I can add that. > > > + > > + double invdet = 1 / det; > > What if det is zero (or very close) ? Then the values get very large :-) Jokes aside, what do you expect to happen? Should we add an optional output parameter for the determinant, so that the caller could check it? Or limit it to something close above 0? > > > + > > + Matrix<T, Rows, Cols> ret; > > + ret[0][0] = (m[1][1] * m[2][2] - m[2][1] * m[1][2]) * invdet; > > + ret[0][1] = (m[0][2] * m[2][1] - m[0][1] * m[2][2]) * invdet; > > + ret[0][2] = (m[0][1] * m[1][2] - m[0][2] * m[1][1]) * invdet; > > + ret[1][0] = (m[1][2] * m[2][0] - m[1][0] * m[2][2]) * invdet; > > + ret[1][1] = (m[0][0] * m[2][2] - m[0][2] * m[2][0]) * invdet; > > + ret[1][2] = (m[1][0] * m[0][2] - m[0][0] * m[1][2]) * invdet; > > + ret[2][0] = (m[1][0] * m[2][1] - m[2][0] * m[1][1]) * invdet; > > + ret[2][1] = (m[2][0] * m[0][1] - m[0][0] * m[2][1]) * invdet; > > + ret[2][2] = (m[0][0] * m[1][1] - m[1][0] * m[0][1]) * invdet; > > That's a lot of calculation for an inline function, please move it to > the .cpp file. Somehow I don't get it. How would you do that? One could possibly move the function definition into the cpp and explicitly instantiate Matrix<float, 3, 3> and Matrix<double, 3, 3> there. Something along the lines: template<typename T, unsigned int Rows, unsigned int Cols> Matrix<T, Rows, Cols> Matrix<T, Rows, Cols>::inverse() { ... } template class Matrix<float, 3, 3>; template class Matrix<double, 3, 3>; I gave that a quick try and failed with "invalid use of incomplete type class". And if it worked, would that be worth it? Best regards, Stefan > > > + return ret; > > + } > > + > > template<typename T2> > > Matrix<T2, Rows, Cols> cast() const > > { > > diff --git a/src/libcamera/matrix.cpp b/src/libcamera/matrix.cpp > > index 91a3f16405a3..b17a1938f1ba 100644 > > --- a/src/libcamera/matrix.cpp > > +++ b/src/libcamera/matrix.cpp > > @@ -77,6 +77,16 @@ LOG_DEFINE_CATEGORY(Matrix) > > * \return Row \a i from the matrix, as a Span > > */ > > > > +/** > > + * \fn Matrix::inverse() const > > + * \brief Compute the inverse of the matrix > > + * > > + * This function computes the inverse of the matrix. It is only implemented for > > + * 3x3 matrices. > > + * > > + * \return The inverse of the matrix > > + */ > > + > > /** > > * \fn template<typename T2> Matrix<T2, Rows, Cols> Matrix::cast() const > > * \brief Cast the matrix to a different type > > -- > Regards, > > Laurent Pinchart
On Mon, Feb 17, 2025 at 01:04:42PM +0100, Stefan Klug wrote: > Hi Laurent, > > Thank you for the review. > > On Mon, Feb 17, 2025 at 12:50:24PM +0200, Laurent Pinchart wrote: > > Hi Stefan, > > > > Thank you for the patch. > > > > On Mon, Feb 17, 2025 at 11:01:44AM +0100, Stefan Klug wrote: > > > For calculations in upcoming algorithm patches, the inverse of a matrix > > > is required. Add an implementation of the inverse() function for 3x3 > > > matrices. > > > > > > Signed-off-by: Stefan Klug <stefan.klug@ideasonboard.com> > > > --- > > > include/libcamera/internal/matrix.h | 24 ++++++++++++++++++++++++ > > > src/libcamera/matrix.cpp | 10 ++++++++++ > > > 2 files changed, 34 insertions(+) > > > > > > diff --git a/include/libcamera/internal/matrix.h b/include/libcamera/internal/matrix.h > > > index ca81dcda93c4..053336e9cfa4 100644 > > > --- a/include/libcamera/internal/matrix.h > > > +++ b/include/libcamera/internal/matrix.h > > > @@ -90,6 +90,30 @@ public: > > > return *this; > > > } > > > > > > + Matrix<T, Rows, Cols> inverse() const > > > + { > > > + static_assert(Rows == 3 && Cols == 3, "Matrix must be 3x3"); > > > > You know what I'll ask, right ? :-) By this comment I meant can the function be generalized to other sizes ? > > Inverting the matrix only makes sense when T is a floating point type. > > I'd add a static assertion or enable_if for this. > > Yes, I'll add a static assert for that. Even though we then remove the > theoretical ability to use that function with a custom defined type. But > as that is most likely never going to happen, I can restrict it to > double/float... > > > > + > > > + const auto &m = *this; > > > + double det = m[0][0] * (m[1][1] * m[2][2] - m[2][1] * m[1][2]) - > > > + m[0][1] * (m[1][0] * m[2][2] - m[1][2] * m[2][0]) + > > > + m[0][2] * (m[1][0] * m[2][1] - m[1][1] * m[2][0]); > > > > Calculating the determinant should be moved to its own function, it's > > useful by itself. > > Yes, I can add that. > > > > + > > > + double invdet = 1 / det; > > > > What if det is zero (or very close) ? > > Then the values get very large :-) Jokes aside, what do you expect to > happen? Should we add an optional output parameter for the determinant, > so that the caller could check it? Or limit it to something close above 0? I'm concerned about divisions by zero. As long as it doesn't crash or cause undefined behaviour, that's fine. For use cases such as inverting CCM matrices I don't think it will be a practical issue, but it would still be good to handle it properly in the API. An output parameter for the determinant can be useful to inform the caller. > > > + > > > + Matrix<T, Rows, Cols> ret; > > > + ret[0][0] = (m[1][1] * m[2][2] - m[2][1] * m[1][2]) * invdet; > > > + ret[0][1] = (m[0][2] * m[2][1] - m[0][1] * m[2][2]) * invdet; > > > + ret[0][2] = (m[0][1] * m[1][2] - m[0][2] * m[1][1]) * invdet; > > > + ret[1][0] = (m[1][2] * m[2][0] - m[1][0] * m[2][2]) * invdet; > > > + ret[1][1] = (m[0][0] * m[2][2] - m[0][2] * m[2][0]) * invdet; > > > + ret[1][2] = (m[1][0] * m[0][2] - m[0][0] * m[1][2]) * invdet; > > > + ret[2][0] = (m[1][0] * m[2][1] - m[2][0] * m[1][1]) * invdet; > > > + ret[2][1] = (m[2][0] * m[0][1] - m[0][0] * m[2][1]) * invdet; > > > + ret[2][2] = (m[0][0] * m[1][1] - m[1][0] * m[0][1]) * invdet; > > > > That's a lot of calculation for an inline function, please move it to > > the .cpp file. > > Somehow I don't get it. How would you do that? One could possibly move > the function definition into the cpp and explicitly instantiate > Matrix<float, 3, 3> and Matrix<double, 3, 3> there. > > Something along the lines: > > template<typename T, unsigned int Rows, unsigned int Cols> > Matrix<T, Rows, Cols> Matrix<T, Rows, Cols>::inverse() > { > ... > } > > template class Matrix<float, 3, 3>; > template class Matrix<double, 3, 3>; > > I gave that a quick try and failed with "invalid use of incomplete type > class". And if it worked, would that be worth it? I'd make a function that operates on data_, and takes the size as an argument. > > > + return ret; > > > + } > > > + > > > template<typename T2> > > > Matrix<T2, Rows, Cols> cast() const > > > { > > > diff --git a/src/libcamera/matrix.cpp b/src/libcamera/matrix.cpp > > > index 91a3f16405a3..b17a1938f1ba 100644 > > > --- a/src/libcamera/matrix.cpp > > > +++ b/src/libcamera/matrix.cpp > > > @@ -77,6 +77,16 @@ LOG_DEFINE_CATEGORY(Matrix) > > > * \return Row \a i from the matrix, as a Span > > > */ > > > > > > +/** > > > + * \fn Matrix::inverse() const > > > + * \brief Compute the inverse of the matrix > > > + * > > > + * This function computes the inverse of the matrix. It is only implemented for > > > + * 3x3 matrices. > > > + * > > > + * \return The inverse of the matrix > > > + */ > > > + > > > /** > > > * \fn template<typename T2> Matrix<T2, Rows, Cols> Matrix::cast() const > > > * \brief Cast the matrix to a different type
diff --git a/include/libcamera/internal/matrix.h b/include/libcamera/internal/matrix.h index ca81dcda93c4..053336e9cfa4 100644 --- a/include/libcamera/internal/matrix.h +++ b/include/libcamera/internal/matrix.h @@ -90,6 +90,30 @@ public: return *this; } + Matrix<T, Rows, Cols> inverse() const + { + static_assert(Rows == 3 && Cols == 3, "Matrix must be 3x3"); + + const auto &m = *this; + double det = m[0][0] * (m[1][1] * m[2][2] - m[2][1] * m[1][2]) - + m[0][1] * (m[1][0] * m[2][2] - m[1][2] * m[2][0]) + + m[0][2] * (m[1][0] * m[2][1] - m[1][1] * m[2][0]); + + double invdet = 1 / det; + + Matrix<T, Rows, Cols> ret; + ret[0][0] = (m[1][1] * m[2][2] - m[2][1] * m[1][2]) * invdet; + ret[0][1] = (m[0][2] * m[2][1] - m[0][1] * m[2][2]) * invdet; + ret[0][2] = (m[0][1] * m[1][2] - m[0][2] * m[1][1]) * invdet; + ret[1][0] = (m[1][2] * m[2][0] - m[1][0] * m[2][2]) * invdet; + ret[1][1] = (m[0][0] * m[2][2] - m[0][2] * m[2][0]) * invdet; + ret[1][2] = (m[1][0] * m[0][2] - m[0][0] * m[1][2]) * invdet; + ret[2][0] = (m[1][0] * m[2][1] - m[2][0] * m[1][1]) * invdet; + ret[2][1] = (m[2][0] * m[0][1] - m[0][0] * m[2][1]) * invdet; + ret[2][2] = (m[0][0] * m[1][1] - m[1][0] * m[0][1]) * invdet; + return ret; + } + template<typename T2> Matrix<T2, Rows, Cols> cast() const { diff --git a/src/libcamera/matrix.cpp b/src/libcamera/matrix.cpp index 91a3f16405a3..b17a1938f1ba 100644 --- a/src/libcamera/matrix.cpp +++ b/src/libcamera/matrix.cpp @@ -77,6 +77,16 @@ LOG_DEFINE_CATEGORY(Matrix) * \return Row \a i from the matrix, as a Span */ +/** + * \fn Matrix::inverse() const + * \brief Compute the inverse of the matrix + * + * This function computes the inverse of the matrix. It is only implemented for + * 3x3 matrices. + * + * \return The inverse of the matrix + */ + /** * \fn template<typename T2> Matrix<T2, Rows, Cols> Matrix::cast() const * \brief Cast the matrix to a different type
For calculations in upcoming algorithm patches, the inverse of a matrix is required. Add an implementation of the inverse() function for 3x3 matrices. Signed-off-by: Stefan Klug <stefan.klug@ideasonboard.com> --- include/libcamera/internal/matrix.h | 24 ++++++++++++++++++++++++ src/libcamera/matrix.cpp | 10 ++++++++++ 2 files changed, 34 insertions(+)