[PATCH 03/10] libcamera: matrix: Add inverse() function

[Stefan Klug]

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 at 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