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

[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 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 ? :-)

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

-- 
Regards,

Laurent Pinchart