[v7,1/3] ipa: libipa: Add Matrix class

20261 diff mbox series

Add a class to represent a Matrix object and operations for adding
matrices, multipling a matrix by a scalar, and multiplying two matrices.

Signed-off-by: Paul Elder <paul.elder@ideasonboard.com>
Reviewed-by: Stefan Klug <stefan.klug@ideasonboard.com>
Reviewed-by: Kieran Bingham <kieran.bingham@ideasonboard.com>

---
Changes in v7:
- fix copyright and license

Changes in v6:
- fix doxygen

Changes in v5:
- add documentation

Changes in v4:
- remove stray semicolons
- add operator<<
- clean up/optimize constructor
- replace get() and set() with operator[] (and a second [] can be used
  as operator[] returns a Span)

Changes in v3:
- fix template parameters of operator* to allow different types for the
  scalar multiplier and the matrix's number type
- clear data in constructors
- fix assert in constructor

Changes v2:
- make rows and columns into template arguments
- initialize to identity matrix on construction
- add getter and setter
- change from struct to class
- fix matrix multiplication
- clean up unused includes
- avoid dereferencing an absent std::optional
---
 src/ipa/libipa/matrix.cpp  | 123 ++++++++++++++++++++++++++
 src/ipa/libipa/matrix.h    | 172 +++++++++++++++++++++++++++++++++++++
 src/ipa/libipa/meson.build |   2 +
 3 files changed, 297 insertions(+)
 create mode 100644 src/ipa/libipa/matrix.cpp
 create mode 100644 src/ipa/libipa/matrix.h

Hi Paul,

Thank you for the patch.

On Tue, Jun 11, 2024 at 11:02:05PM +0900, Paul Elder wrote:
> Add a class to represent a Matrix object and operations for adding
> matrices, multipling a matrix by a scalar, and multiplying two matrices.
> 
> Signed-off-by: Paul Elder <paul.elder@ideasonboard.com>
> Reviewed-by: Stefan Klug <stefan.klug@ideasonboard.com>
> Reviewed-by: Kieran Bingham <kieran.bingham@ideasonboard.com>
> 
> ---
> Changes in v7:
> - fix copyright and license
> 
> Changes in v6:
> - fix doxygen
> 
> Changes in v5:
> - add documentation
> 
> Changes in v4:
> - remove stray semicolons
> - add operator<<
> - clean up/optimize constructor
> - replace get() and set() with operator[] (and a second [] can be used
>   as operator[] returns a Span)
> 
> Changes in v3:
> - fix template parameters of operator* to allow different types for the
>   scalar multiplier and the matrix's number type
> - clear data in constructors
> - fix assert in constructor
> 
> Changes v2:
> - make rows and columns into template arguments
> - initialize to identity matrix on construction
> - add getter and setter
> - change from struct to class
> - fix matrix multiplication
> - clean up unused includes
> - avoid dereferencing an absent std::optional
> ---
>  src/ipa/libipa/matrix.cpp  | 123 ++++++++++++++++++++++++++
>  src/ipa/libipa/matrix.h    | 172 +++++++++++++++++++++++++++++++++++++
>  src/ipa/libipa/meson.build |   2 +
>  3 files changed, 297 insertions(+)
>  create mode 100644 src/ipa/libipa/matrix.cpp
>  create mode 100644 src/ipa/libipa/matrix.h
> 
> diff --git a/src/ipa/libipa/matrix.cpp b/src/ipa/libipa/matrix.cpp
> new file mode 100644
> index 000000000000..02091a0b4ce1
> --- /dev/null
> +++ b/src/ipa/libipa/matrix.cpp
> @@ -0,0 +1,123 @@
> +/* SPDX-License-Identifier: LGPL-2.1-or-later */
> +/*
> + * Copyright (C) 2024, Paul Elder <paul.elder@ideasonboard.com>
> + *
> + * Matrix and related operations
> + */
> +
> +#include "matrix.h"
> +
> +#include <libcamera/base/log.h>
> +
> +/**
> + * \file matrix.h
> + * \brief Matrix class
> + */
> +
> +namespace libcamera {
> +
> +LOG_DEFINE_CATEGORY(Matrix)
> +
> +namespace ipa {
> +
> +/**
> + * \class Matrix
> + * \brief Matrix class
> + * \tparam T Type of numerical values to be stored in the matrix
> + * \tparam R Number of rows in the matrix
> + * \tparam C Number of columns in the matrix
> + */
> +
> +/**
> + * \fn Matrix::Matrix()
> + * \brief Construct an identity matrix

I would have expected the default constructor to construct a zero
matrix. Is there a reason to go for identity instead ?

As identity matrices are useful too, I would add a static
Matrix::identity() function that constructs and returns an identity
matrix.

> + */
> +
> +/**
> + * \fn Matrix::Matrix(const std::vector<T> &data)

Turning the argument into a Span would allow passing any type of
contiguous data (vector, array, C array, ...) without requiring the
caller to construct a std::vector.

> + * \brief Construct matrix from supplied data

s/matrix/a matrix/

> + * \param data Data from which to construct a matrix

\param[in]

Same below.

> + *
> + * \a data is a one-dimensional vector and will be turned into a matrix in
> + * row-major order. The size of \a data must be equal to the product of the
> + * number of rows and columns of the matrix (RxC).
> + */
> +
> +/**
> + * \fn Matrix::readYaml

s/$/()/

Same below.

> + * \brief Populate the matrix with yaml data
> + * \param yaml Yaml data to populate the matrix with
> + *
> + * Any existing data in the matrix will be overwritten. The size of the data
> + * read from \a yaml must be equal to the product of the number of rows and
> + * columns of the matrix (RxC).
> + *
> + * The yaml data is expected to be a list with elements of type T.

in row-major order.

> + *
> + * \return 0 on success, negative error code otherwise
> + */
> +
> +/**
> + * \fn Matrix::toString
> + * \brief Assemble and return a string describing the matrix
> + * \return A string describing the matrix
> + */
> +
> +/**
> + * \fn Span<const T, C> Matrix::operator[](size_t i) const
> + * \brief Index to a row in the matrix
> + * \param i Index of row to retrieve
> + *
> + * This operator[] returns a Span, which can then be indexed into again with
> + * another operator[], allowing a convenient m[i][j] to access elements of the
> + * matrix. Note that the lifetime of the Span returned by this first-level
> + * operator[] is bound to that of the Matrix itself, so it is not recommended
> + * to save the Span that is the result of this operator[].

That's a clever idea, I like it.

Longer term we may need a different (and
unfortunately much more complex) implementation to allow using Matrix[i]
in APIs that expect a Vector (which itself may need to be reimplemented
as a matrix with a single column). And the next thing we know is that
we'll reimplement uBLAS or Eigen3 :-) I'll cross my fingers and wish we
won't reach that point.

https://stackoverflow.com/questions/1380371/what-are-the-most-widely-used-c-vector-matrix-math-linear-algebra-libraries-a

"It seems that many projects slowly come upon a need to do matrix math,
and fall into the trap of first building some vector classes and slowly
adding in functionality until they get caught building a half-assed
custom linear algebra library, and depending on it."

> + *
> + * \return Row \a i from the matrix, as a Span
> + */
> +
> +/**
> + * \fn Matrix::operator[](size_t i)
> + * \copydoc Matrix::operator[](size_t i) const
> + */
> +
> +/**
> + * \fn Matrix::Matrix<U, R, C> operator*(T d, const Matrix<U, R, C> &m)
> + * \brief Scalar product
> + * \tparam T Type of the numerical scalar value
> + * \tparam U Type of numerical values in the matrix
> + * \tparam R Number of rows in the matrix
> + * \tparam C Number of columns in the matrix
> + * \param d Scalar
> + * \param m Matrix
> + * \return Product of scalar \a d and matrix \a m
> + */
> +
> +/**
> + * \fn Matrix<T, R1, C2> operator*(const Matrix<T, R1, C1> &m1, const Matrix<T, R2, C2> &m2)
> + * \brief Matrix multiplication
> + * \tparam T Type of numerical values in the matrices
> + * \tparam R1 Number of rows in the first matrix
> + * \tparam C1 Number of columns in the first matrix
> + * \tparam R2 Number of rows in the second matrix
> + * \tparam C2 Number of columns in the second matrix
> + * \param m1 Multiplicand matrix
> + * \param m2 Multiplier matrix
> + * \return Matrix product of matrices \a m1 and \a m2
> + */
> +
> +/**
> + * \fn Matrix<T, R, C> operator+(const Matrix<T, R, C> &m1, const Matrix<T, R, C> &m2)
> + * \brief Matrix addition
> + * \tparam T Type of numerical values in the matrices
> + * \tparam R Number of rows in the matrices
> + * \tparam C Number of columns in the matrices
> + * \param m1 Summand matrix
> + * \param m2 Summand matrix
> + * \return Matrix sum of matrices \a m1 and \a m2
> + */
> +
> +} /* namespace ipa */
> +
> +} /* namespace libcamera */
> diff --git a/src/ipa/libipa/matrix.h b/src/ipa/libipa/matrix.h
> new file mode 100644
> index 000000000000..90eaea03bd14
> --- /dev/null
> +++ b/src/ipa/libipa/matrix.h
> @@ -0,0 +1,172 @@
> +/* SPDX-License-Identifier: LGPL-2.1-or-later */
> +/*
> + * Copyright (C) 2024, Paul Elder <paul.elder@ideasonboard.com>
> + *
> + * Matrix and related operations
> + */
> +#pragma once
> +
> +#include <algorithm>
> +#include <cmath>
> +#include <sstream>
> +#include <vector>
> +
> +#include <libcamera/base/log.h>
> +#include <libcamera/base/span.h>
> +
> +#include "libcamera/internal/yaml_parser.h"
> +
> +namespace libcamera {
> +
> +LOG_DECLARE_CATEGORY(Matrix)
> +
> +namespace ipa {
> +
> +#ifndef __DOXYGEN__
> +template<typename T, unsigned int R, unsigned int C,

Following the Vector class,

template<typename T, unsigned int Rows, unsigned int Cols,

> +	 std::enable_if_t<std::is_arithmetic_v<T>> * = nullptr>
> +#else
> +template<typename T, unsigned int R, unsigned int C>
> +#endif /* __DOXYGEN__ */
> +class Matrix
> +{
> +public:
> +	Matrix()
> +		: data_(R * C, static_cast<T>(false))

false ? That's a strange one, anything wrong with 0 ?

> +	{
> +		for (size_t i = 0; i < std::min(R, C); i++)
> +			(*this)[i][i] = static_cast<T>(true);

And 1 here.

> +	}
> +
> +	Matrix(const std::vector<T> &data)
> +	{
> +		ASSERT(data.size() == R * C);
> +
> +		data_.clear();
> +		for (const T &x : data)
> +			data_.push_back(x);

That will perform quite a few reallocations. I think you can simply
write

		std::copy(data.begin(), data.end(), data_.begin());

> +	}
> +
> +	~Matrix() = default;
> +
> +	int readYaml(const libcamera::YamlObject &yaml)
> +	{
> +		if (yaml.size() != R * C) {
> +			LOG(Matrix, Error)
> +				<< "Wrong number of values in matrix: expected "
> +				<< R * C << ", got " << yaml.size();
> +			return -EINVAL;
> +		}
> +
> +		unsigned int i = 0;
> +		for (const auto &x : yaml.asList()) {
> +			auto value = x.get<T>();
> +			if (!value) {
> +				LOG(Matrix, Error) << "Failed to read matrix value";
> +				return -EINVAL;
> +			}
> +
> +			data_[i++] = *value;
> +		}
> +
> +		return 0;
> +	}
> +
> +	const std::string toString() const
> +	{
> +		std::stringstream out;
> +
> +		out << "Matrix { ";
> +		for (unsigned int i = 0; i < R; i++) {
> +			out << "[ ";
> +			for (unsigned int j = 0; j < C; j++) {
> +				out << (*this)[i][j];
> +				out << ((j + 1 < C) ? ", " : " ");
> +			}
> +			out << ((i + 1 < R) ? "], " : "]");
> +		}
> +		out << " }";
> +
> +		return out.str();
> +	}
> +
> +	Span<const T, C> operator[](size_t i) const
> +	{
> +		return Span<const T, C>{ &data_.data()[i * C], C };
> +	}
> +
> +	Span<T, C> operator[](size_t i)
> +	{
> +		return Span<T, C>{ &data_.data()[i * C], C };
> +	}
> +
> +private:
> +	std::vector<T> data_;

Should this be a std::array<T, Rows * Cols>, like for the Vector class ?

> +};
> +
> +#ifndef __DOXYGEN__
> +template<typename T, typename U, unsigned int R, unsigned int C,
> +	 std::enable_if_t<std::is_arithmetic_v<T> && std::is_arithmetic_v<U>> * = nullptr>

The latter is ensured by the Matrix class.

You need an

#else
template<typename T, typename U, unsigned int R, unsigned int C>

Same below.

> +#endif /* __DOXYGEN__ */
> +Matrix<U, R, C> operator*(T d, const Matrix<U, R, C> &m)

I'm surprised by the order of the operands, as I would have expected
users to write

	Matrix m;
	Matrix r = m * factor;

and not

	Matrix m;
	Matrix r = factor * m;

I recommend supporting both, you can easily implement one based on the
other.

It would also be useful to support the following with a member operator:

	Matrix m;
	Matrix r;
	r *= factor;

> +{
> +	Matrix<U, R, C> result;
> +
> +	for (unsigned int i = 0; i < R; i++)
> +		for (unsigned int j = 0; j < C; j++)
> +			result[i][j] = d * m[i][j];

I think

	for (unsigned int i = 0; i < R * C; i++)
		result.data_[i] = d * m.data_[i];

would be more efficient.

> +
> +	return result;
> +}
> +
> +#ifndef __DOXYGEN__
> +template<typename T,
> +	 unsigned int R1, unsigned int C1,
> +	 unsigned int R2, unsigned int C2,
> +	 std::enable_if_t<std::is_arithmetic_v<T> && C1 == R2> * = nullptr>
> +#endif /* __DOXYGEN__ */
> +Matrix<T, R1, C2> operator*(const Matrix<T, R1, C1> &m1, const Matrix<T, R2, C2> &m2)

This can also be simplified in a similar way as above.

> +{
> +	Matrix<T, R1, C2> result;
> +
> +	for (unsigned int i = 0; i < R1; i++) {
> +		for (unsigned int j = 0; j < C2; j++) {
> +			T sum = 0;
> +
> +			for (unsigned int k = 0; k < C1; k++)
> +				sum += m1[i][k] * m2[k][j];
> +
> +			result[i][j] = sum;
> +		}
> +	}
> +
> +	return result;
> +}
> +
> +#ifndef __DOXYGEN__
> +template<typename T, unsigned int R, unsigned int C,
> +	 std::enable_if_t<std::is_arithmetic_v<T>> * = nullptr>

Ensured by the Matrix type.

> +#endif /* __DOXYGEN__ */
> +Matrix<T, R, C> operator+(const Matrix<T, R, C> &m1, const Matrix<T, R, C> &m2)
> +{
> +	Matrix<T, R, C> result;
> +
> +	for (unsigned int i = 0; i < R; i++)
> +		for (unsigned int j = 0; j < C; j++)
> +			result[i][j] = m1[i][j] + m2[i][j];

I think

	for (unsigned int i = 0; i < R * C; i++)
		result.data_[i] = m1.data_[i] + m2.data_[i];

would be more efficient.

> +
> +	return result;
> +}
> +
> +} /* namespace ipa */
> +
> +#ifndef __DOXYGEN__
> +template<typename T, unsigned int R, unsigned int C>
> +std::ostream &operator<<(std::ostream &out, const ipa::Matrix<T, R, C> &m)
> +{
> +	out << m.toString();
> +	return out;
> +}
> +#endif /* __DOXYGEN__ */
> +
> +} /* namespace libcamera */
> diff --git a/src/ipa/libipa/meson.build b/src/ipa/libipa/meson.build
> index b663afc7d9fe..067d0d273e0a 100644
> --- a/src/ipa/libipa/meson.build
> +++ b/src/ipa/libipa/meson.build
> @@ -7,6 +7,7 @@ libipa_headers = files([
>      'exposure_mode_helper.h',
>      'fc_queue.h',
>      'histogram.h',
> +    'matrix.h',
>      'module.h',
>  ])
>  
> @@ -17,6 +18,7 @@ libipa_sources = files([
>      'exposure_mode_helper.cpp',
>      'fc_queue.cpp',
>      'histogram.cpp',
> +    'matrix.cpp',
>      'module.cpp',
>  ])
>

On Wed, Jun 12, 2024 at 02:33:33AM +0300, Laurent Pinchart wrote:
> Hi Paul,
> 
> Thank you for the patch.
> 
> On Tue, Jun 11, 2024 at 11:02:05PM +0900, Paul Elder wrote:
> > Add a class to represent a Matrix object and operations for adding
> > matrices, multipling a matrix by a scalar, and multiplying two matrices.
> > 
> > Signed-off-by: Paul Elder <paul.elder@ideasonboard.com>
> > Reviewed-by: Stefan Klug <stefan.klug@ideasonboard.com>
> > Reviewed-by: Kieran Bingham <kieran.bingham@ideasonboard.com>
> > 
> > ---
> > Changes in v7:
> > - fix copyright and license
> > 
> > Changes in v6:
> > - fix doxygen
> > 
> > Changes in v5:
> > - add documentation
> > 
> > Changes in v4:
> > - remove stray semicolons
> > - add operator<<
> > - clean up/optimize constructor
> > - replace get() and set() with operator[] (and a second [] can be used
> >   as operator[] returns a Span)
> > 
> > Changes in v3:
> > - fix template parameters of operator* to allow different types for the
> >   scalar multiplier and the matrix's number type
> > - clear data in constructors
> > - fix assert in constructor
> > 
> > Changes v2:
> > - make rows and columns into template arguments
> > - initialize to identity matrix on construction
> > - add getter and setter
> > - change from struct to class
> > - fix matrix multiplication
> > - clean up unused includes
> > - avoid dereferencing an absent std::optional
> > ---
> >  src/ipa/libipa/matrix.cpp  | 123 ++++++++++++++++++++++++++
> >  src/ipa/libipa/matrix.h    | 172 +++++++++++++++++++++++++++++++++++++
> >  src/ipa/libipa/meson.build |   2 +
> >  3 files changed, 297 insertions(+)
> >  create mode 100644 src/ipa/libipa/matrix.cpp
> >  create mode 100644 src/ipa/libipa/matrix.h
> > 
> > diff --git a/src/ipa/libipa/matrix.cpp b/src/ipa/libipa/matrix.cpp
> > new file mode 100644
> > index 000000000000..02091a0b4ce1
> > --- /dev/null
> > +++ b/src/ipa/libipa/matrix.cpp
> > @@ -0,0 +1,123 @@
> > +/* SPDX-License-Identifier: LGPL-2.1-or-later */
> > +/*
> > + * Copyright (C) 2024, Paul Elder <paul.elder@ideasonboard.com>
> > + *
> > + * Matrix and related operations
> > + */
> > +
> > +#include "matrix.h"
> > +
> > +#include <libcamera/base/log.h>
> > +
> > +/**
> > + * \file matrix.h
> > + * \brief Matrix class
> > + */
> > +
> > +namespace libcamera {
> > +
> > +LOG_DEFINE_CATEGORY(Matrix)
> > +
> > +namespace ipa {
> > +
> > +/**
> > + * \class Matrix
> > + * \brief Matrix class
> > + * \tparam T Type of numerical values to be stored in the matrix
> > + * \tparam R Number of rows in the matrix
> > + * \tparam C Number of columns in the matrix
> > + */
> > +
> > +/**
> > + * \fn Matrix::Matrix()
> > + * \brief Construct an identity matrix
> 
> I would have expected the default constructor to construct a zero
> matrix. Is there a reason to go for identity instead ?
> 
> As identity matrices are useful too, I would add a static
> Matrix::identity() function that constructs and returns an identity
> matrix.
> 
> > + */
> > +
> > +/**
> > + * \fn Matrix::Matrix(const std::vector<T> &data)
> 
> Turning the argument into a Span would allow passing any type of
> contiguous data (vector, array, C array, ...) without requiring the
> caller to construct a std::vector.
> 

That's what I wanted to do but it was refusing to work.

> > + * \brief Construct matrix from supplied data
> 
> s/matrix/a matrix/
> 
> > + * \param data Data from which to construct a matrix
> 
> \param[in]
> 
> Same below.
> 
> > + *
> > + * \a data is a one-dimensional vector and will be turned into a matrix in
> > + * row-major order. The size of \a data must be equal to the product of the
> > + * number of rows and columns of the matrix (RxC).
> > + */
> > +
> > +/**
> > + * \fn Matrix::readYaml
> 
> s/$/()/
> 
> Same below.
> 
> > + * \brief Populate the matrix with yaml data
> > + * \param yaml Yaml data to populate the matrix with
> > + *
> > + * Any existing data in the matrix will be overwritten. The size of the data
> > + * read from \a yaml must be equal to the product of the number of rows and
> > + * columns of the matrix (RxC).
> > + *
> > + * The yaml data is expected to be a list with elements of type T.
> 
> in row-major order.
> 
> > + *
> > + * \return 0 on success, negative error code otherwise
> > + */
> > +
> > +/**
> > + * \fn Matrix::toString
> > + * \brief Assemble and return a string describing the matrix
> > + * \return A string describing the matrix
> > + */
> > +
> > +/**
> > + * \fn Span<const T, C> Matrix::operator[](size_t i) const
> > + * \brief Index to a row in the matrix
> > + * \param i Index of row to retrieve
> > + *
> > + * This operator[] returns a Span, which can then be indexed into again with
> > + * another operator[], allowing a convenient m[i][j] to access elements of the
> > + * matrix. Note that the lifetime of the Span returned by this first-level
> > + * operator[] is bound to that of the Matrix itself, so it is not recommended
> > + * to save the Span that is the result of this operator[].
> 
> That's a clever idea, I like it.
> 
> Longer term we may need a different (and
> unfortunately much more complex) implementation to allow using Matrix[i]
> in APIs that expect a Vector (which itself may need to be reimplemented
> as a matrix with a single column). And the next thing we know is that
> we'll reimplement uBLAS or Eigen3 :-) I'll cross my fingers and wish we
> won't reach that point.
> 
> https://stackoverflow.com/questions/1380371/what-are-the-most-widely-used-c-vector-matrix-math-linear-algebra-libraries-a
> 
> "It seems that many projects slowly come upon a need to do matrix math,
> and fall into the trap of first building some vector classes and slowly
> adding in functionality until they get caught building a half-assed
> custom linear algebra library, and depending on it."
> 

Oops... :)

> > + *
> > + * \return Row \a i from the matrix, as a Span
> > + */
> > +
> > +/**
> > + * \fn Matrix::operator[](size_t i)
> > + * \copydoc Matrix::operator[](size_t i) const
> > + */
> > +
> > +/**
> > + * \fn Matrix::Matrix<U, R, C> operator*(T d, const Matrix<U, R, C> &m)
> > + * \brief Scalar product
> > + * \tparam T Type of the numerical scalar value
> > + * \tparam U Type of numerical values in the matrix
> > + * \tparam R Number of rows in the matrix
> > + * \tparam C Number of columns in the matrix
> > + * \param d Scalar
> > + * \param m Matrix
> > + * \return Product of scalar \a d and matrix \a m
> > + */
> > +
> > +/**
> > + * \fn Matrix<T, R1, C2> operator*(const Matrix<T, R1, C1> &m1, const Matrix<T, R2, C2> &m2)
> > + * \brief Matrix multiplication
> > + * \tparam T Type of numerical values in the matrices
> > + * \tparam R1 Number of rows in the first matrix
> > + * \tparam C1 Number of columns in the first matrix
> > + * \tparam R2 Number of rows in the second matrix
> > + * \tparam C2 Number of columns in the second matrix
> > + * \param m1 Multiplicand matrix
> > + * \param m2 Multiplier matrix
> > + * \return Matrix product of matrices \a m1 and \a m2
> > + */
> > +
> > +/**
> > + * \fn Matrix<T, R, C> operator+(const Matrix<T, R, C> &m1, const Matrix<T, R, C> &m2)
> > + * \brief Matrix addition
> > + * \tparam T Type of numerical values in the matrices
> > + * \tparam R Number of rows in the matrices
> > + * \tparam C Number of columns in the matrices
> > + * \param m1 Summand matrix
> > + * \param m2 Summand matrix
> > + * \return Matrix sum of matrices \a m1 and \a m2
> > + */
> > +
> > +} /* namespace ipa */
> > +
> > +} /* namespace libcamera */
> > diff --git a/src/ipa/libipa/matrix.h b/src/ipa/libipa/matrix.h
> > new file mode 100644
> > index 000000000000..90eaea03bd14
> > --- /dev/null
> > +++ b/src/ipa/libipa/matrix.h
> > @@ -0,0 +1,172 @@
> > +/* SPDX-License-Identifier: LGPL-2.1-or-later */
> > +/*
> > + * Copyright (C) 2024, Paul Elder <paul.elder@ideasonboard.com>
> > + *
> > + * Matrix and related operations
> > + */
> > +#pragma once
> > +
> > +#include <algorithm>
> > +#include <cmath>
> > +#include <sstream>
> > +#include <vector>
> > +
> > +#include <libcamera/base/log.h>
> > +#include <libcamera/base/span.h>
> > +
> > +#include "libcamera/internal/yaml_parser.h"
> > +
> > +namespace libcamera {
> > +
> > +LOG_DECLARE_CATEGORY(Matrix)
> > +
> > +namespace ipa {
> > +
> > +#ifndef __DOXYGEN__
> > +template<typename T, unsigned int R, unsigned int C,
> 
> Following the Vector class,
> 
> template<typename T, unsigned int Rows, unsigned int Cols,
> 
> > +	 std::enable_if_t<std::is_arithmetic_v<T>> * = nullptr>
> > +#else
> > +template<typename T, unsigned int R, unsigned int C>
> > +#endif /* __DOXYGEN__ */
> > +class Matrix
> > +{
> > +public:
> > +	Matrix()
> > +		: data_(R * C, static_cast<T>(false))
> 
> false ? That's a strange one, anything wrong with 0 ?
> 
> > +	{
> > +		for (size_t i = 0; i < std::min(R, C); i++)
> > +			(*this)[i][i] = static_cast<T>(true);
> 
> And 1 here.
> 
> > +	}
> > +
> > +	Matrix(const std::vector<T> &data)
> > +	{
> > +		ASSERT(data.size() == R * C);
> > +
> > +		data_.clear();
> > +		for (const T &x : data)
> > +			data_.push_back(x);
> 
> That will perform quite a few reallocations. I think you can simply
> write
> 
> 		std::copy(data.begin(), data.end(), data_.begin());
> 
> > +	}
> > +
> > +	~Matrix() = default;
> > +
> > +	int readYaml(const libcamera::YamlObject &yaml)
> > +	{
> > +		if (yaml.size() != R * C) {
> > +			LOG(Matrix, Error)
> > +				<< "Wrong number of values in matrix: expected "
> > +				<< R * C << ", got " << yaml.size();
> > +			return -EINVAL;
> > +		}
> > +
> > +		unsigned int i = 0;
> > +		for (const auto &x : yaml.asList()) {
> > +			auto value = x.get<T>();
> > +			if (!value) {
> > +				LOG(Matrix, Error) << "Failed to read matrix value";
> > +				return -EINVAL;
> > +			}
> > +
> > +			data_[i++] = *value;
> > +		}
> > +
> > +		return 0;
> > +	}
> > +
> > +	const std::string toString() const
> > +	{
> > +		std::stringstream out;
> > +
> > +		out << "Matrix { ";
> > +		for (unsigned int i = 0; i < R; i++) {
> > +			out << "[ ";
> > +			for (unsigned int j = 0; j < C; j++) {
> > +				out << (*this)[i][j];
> > +				out << ((j + 1 < C) ? ", " : " ");
> > +			}
> > +			out << ((i + 1 < R) ? "], " : "]");
> > +		}
> > +		out << " }";
> > +
> > +		return out.str();
> > +	}
> > +
> > +	Span<const T, C> operator[](size_t i) const
> > +	{
> > +		return Span<const T, C>{ &data_.data()[i * C], C };
> > +	}
> > +
> > +	Span<T, C> operator[](size_t i)
> > +	{
> > +		return Span<T, C>{ &data_.data()[i * C], C };
> > +	}
> > +
> > +private:
> > +	std::vector<T> data_;
> 
> Should this be a std::array<T, Rows * Cols>, like for the Vector class ?
> 
> > +};
> > +
> > +#ifndef __DOXYGEN__
> > +template<typename T, typename U, unsigned int R, unsigned int C,
> > +	 std::enable_if_t<std::is_arithmetic_v<T> && std::is_arithmetic_v<U>> * = nullptr>
> 
> The latter is ensured by the Matrix class.
> 
> You need an
> 
> #else
> template<typename T, typename U, unsigned int R, unsigned int C>
> 
> Same below.
> 
> > +#endif /* __DOXYGEN__ */
> > +Matrix<U, R, C> operator*(T d, const Matrix<U, R, C> &m)
> 
> I'm surprised by the order of the operands, as I would have expected
> users to write
> 
> 	Matrix m;
> 	Matrix r = m * factor;
> 
> and not
> 
> 	Matrix m;
> 	Matrix r = factor * m;
> 
> I recommend supporting both, you can easily implement one based on the
> other.
> 
> It would also be useful to support the following with a member operator:
> 
> 	Matrix m;
> 	Matrix r;
> 	r *= factor;
> 
> > +{
> > +	Matrix<U, R, C> result;
> > +
> > +	for (unsigned int i = 0; i < R; i++)
> > +		for (unsigned int j = 0; j < C; j++)
> > +			result[i][j] = d * m[i][j];
> 
> I think
> 
> 	for (unsigned int i = 0; i < R * C; i++)
> 		result.data_[i] = d * m.data_[i];
> 
> would be more efficient.

Yeah but data_ is private...


Paul

> 
> > +
> > +	return result;
> > +}
> > +
> > +#ifndef __DOXYGEN__
> > +template<typename T,
> > +	 unsigned int R1, unsigned int C1,
> > +	 unsigned int R2, unsigned int C2,
> > +	 std::enable_if_t<std::is_arithmetic_v<T> && C1 == R2> * = nullptr>
> > +#endif /* __DOXYGEN__ */
> > +Matrix<T, R1, C2> operator*(const Matrix<T, R1, C1> &m1, const Matrix<T, R2, C2> &m2)
> 
> This can also be simplified in a similar way as above.
> 
> > +{
> > +	Matrix<T, R1, C2> result;
> > +
> > +	for (unsigned int i = 0; i < R1; i++) {
> > +		for (unsigned int j = 0; j < C2; j++) {
> > +			T sum = 0;
> > +
> > +			for (unsigned int k = 0; k < C1; k++)
> > +				sum += m1[i][k] * m2[k][j];
> > +
> > +			result[i][j] = sum;
> > +		}
> > +	}
> > +
> > +	return result;
> > +}
> > +
> > +#ifndef __DOXYGEN__
> > +template<typename T, unsigned int R, unsigned int C,
> > +	 std::enable_if_t<std::is_arithmetic_v<T>> * = nullptr>
> 
> Ensured by the Matrix type.
> 
> > +#endif /* __DOXYGEN__ */
> > +Matrix<T, R, C> operator+(const Matrix<T, R, C> &m1, const Matrix<T, R, C> &m2)
> > +{
> > +	Matrix<T, R, C> result;
> > +
> > +	for (unsigned int i = 0; i < R; i++)
> > +		for (unsigned int j = 0; j < C; j++)
> > +			result[i][j] = m1[i][j] + m2[i][j];
> 
> I think
> 
> 	for (unsigned int i = 0; i < R * C; i++)
> 		result.data_[i] = m1.data_[i] + m2.data_[i];
> 
> would be more efficient.
> 
> > +
> > +	return result;
> > +}
> > +
> > +} /* namespace ipa */
> > +
> > +#ifndef __DOXYGEN__
> > +template<typename T, unsigned int R, unsigned int C>
> > +std::ostream &operator<<(std::ostream &out, const ipa::Matrix<T, R, C> &m)
> > +{
> > +	out << m.toString();
> > +	return out;
> > +}
> > +#endif /* __DOXYGEN__ */
> > +
> > +} /* namespace libcamera */
> > diff --git a/src/ipa/libipa/meson.build b/src/ipa/libipa/meson.build
> > index b663afc7d9fe..067d0d273e0a 100644
> > --- a/src/ipa/libipa/meson.build
> > +++ b/src/ipa/libipa/meson.build
> > @@ -7,6 +7,7 @@ libipa_headers = files([
> >      'exposure_mode_helper.h',
> >      'fc_queue.h',
> >      'histogram.h',
> > +    'matrix.h',
> >      'module.h',
> >  ])
> >  
> > @@ -17,6 +18,7 @@ libipa_sources = files([
> >      'exposure_mode_helper.cpp',
> >      'fc_queue.cpp',
> >      'histogram.cpp',
> > +    'matrix.cpp',
> >      'module.cpp',
> >  ])
> >

On Wed, Jun 12, 2024 at 04:59:39PM +0900, Paul Elder wrote:
> On Wed, Jun 12, 2024 at 02:33:33AM +0300, Laurent Pinchart wrote:
> > On Tue, Jun 11, 2024 at 11:02:05PM +0900, Paul Elder wrote:
> > > Add a class to represent a Matrix object and operations for adding
> > > matrices, multipling a matrix by a scalar, and multiplying two matrices.
> > > 
> > > Signed-off-by: Paul Elder <paul.elder@ideasonboard.com>
> > > Reviewed-by: Stefan Klug <stefan.klug@ideasonboard.com>
> > > Reviewed-by: Kieran Bingham <kieran.bingham@ideasonboard.com>
> > > 
> > > ---
> > > Changes in v7:
> > > - fix copyright and license
> > > 
> > > Changes in v6:
> > > - fix doxygen
> > > 
> > > Changes in v5:
> > > - add documentation
> > > 
> > > Changes in v4:
> > > - remove stray semicolons
> > > - add operator<<
> > > - clean up/optimize constructor
> > > - replace get() and set() with operator[] (and a second [] can be used
> > >   as operator[] returns a Span)
> > > 
> > > Changes in v3:
> > > - fix template parameters of operator* to allow different types for the
> > >   scalar multiplier and the matrix's number type
> > > - clear data in constructors
> > > - fix assert in constructor
> > > 
> > > Changes v2:
> > > - make rows and columns into template arguments
> > > - initialize to identity matrix on construction
> > > - add getter and setter
> > > - change from struct to class
> > > - fix matrix multiplication
> > > - clean up unused includes
> > > - avoid dereferencing an absent std::optional
> > > ---
> > >  src/ipa/libipa/matrix.cpp  | 123 ++++++++++++++++++++++++++
> > >  src/ipa/libipa/matrix.h    | 172 +++++++++++++++++++++++++++++++++++++
> > >  src/ipa/libipa/meson.build |   2 +
> > >  3 files changed, 297 insertions(+)
> > >  create mode 100644 src/ipa/libipa/matrix.cpp
> > >  create mode 100644 src/ipa/libipa/matrix.h
> > > 
> > > diff --git a/src/ipa/libipa/matrix.cpp b/src/ipa/libipa/matrix.cpp
> > > new file mode 100644
> > > index 000000000000..02091a0b4ce1
> > > --- /dev/null
> > > +++ b/src/ipa/libipa/matrix.cpp
> > > @@ -0,0 +1,123 @@
> > > +/* SPDX-License-Identifier: LGPL-2.1-or-later */
> > > +/*
> > > + * Copyright (C) 2024, Paul Elder <paul.elder@ideasonboard.com>
> > > + *
> > > + * Matrix and related operations
> > > + */
> > > +
> > > +#include "matrix.h"
> > > +
> > > +#include <libcamera/base/log.h>
> > > +
> > > +/**
> > > + * \file matrix.h
> > > + * \brief Matrix class
> > > + */
> > > +
> > > +namespace libcamera {
> > > +
> > > +LOG_DEFINE_CATEGORY(Matrix)
> > > +
> > > +namespace ipa {
> > > +
> > > +/**
> > > + * \class Matrix
> > > + * \brief Matrix class
> > > + * \tparam T Type of numerical values to be stored in the matrix
> > > + * \tparam R Number of rows in the matrix
> > > + * \tparam C Number of columns in the matrix
> > > + */
> > > +
> > > +/**
> > > + * \fn Matrix::Matrix()
> > > + * \brief Construct an identity matrix
> > 
> > I would have expected the default constructor to construct a zero
> > matrix. Is there a reason to go for identity instead ?
> > 
> > As identity matrices are useful too, I would add a static
> > Matrix::identity() function that constructs and returns an identity
> > matrix.
> > 
> > > + */
> > > +
> > > +/**
> > > + * \fn Matrix::Matrix(const std::vector<T> &data)
> > 
> > Turning the argument into a Span would allow passing any type of
> > contiguous data (vector, array, C array, ...) without requiring the
> > caller to construct a std::vector.
> 
> That's what I wanted to do but it was refusing to work.

Show me the failure and I'll help :-)

> > > + * \brief Construct matrix from supplied data
> > 
> > s/matrix/a matrix/
> > 
> > > + * \param data Data from which to construct a matrix
> > 
> > \param[in]
> > 
> > Same below.
> > 
> > > + *
> > > + * \a data is a one-dimensional vector and will be turned into a matrix in
> > > + * row-major order. The size of \a data must be equal to the product of the
> > > + * number of rows and columns of the matrix (RxC).
> > > + */
> > > +
> > > +/**
> > > + * \fn Matrix::readYaml
> > 
> > s/$/()/
> > 
> > Same below.
> > 
> > > + * \brief Populate the matrix with yaml data
> > > + * \param yaml Yaml data to populate the matrix with
> > > + *
> > > + * Any existing data in the matrix will be overwritten. The size of the data
> > > + * read from \a yaml must be equal to the product of the number of rows and
> > > + * columns of the matrix (RxC).
> > > + *
> > > + * The yaml data is expected to be a list with elements of type T.
> > 
> > in row-major order.
> > 
> > > + *
> > > + * \return 0 on success, negative error code otherwise
> > > + */
> > > +
> > > +/**
> > > + * \fn Matrix::toString
> > > + * \brief Assemble and return a string describing the matrix
> > > + * \return A string describing the matrix
> > > + */
> > > +
> > > +/**
> > > + * \fn Span<const T, C> Matrix::operator[](size_t i) const
> > > + * \brief Index to a row in the matrix
> > > + * \param i Index of row to retrieve
> > > + *
> > > + * This operator[] returns a Span, which can then be indexed into again with
> > > + * another operator[], allowing a convenient m[i][j] to access elements of the
> > > + * matrix. Note that the lifetime of the Span returned by this first-level
> > > + * operator[] is bound to that of the Matrix itself, so it is not recommended
> > > + * to save the Span that is the result of this operator[].
> > 
> > That's a clever idea, I like it.
> > 
> > Longer term we may need a different (and
> > unfortunately much more complex) implementation to allow using Matrix[i]
> > in APIs that expect a Vector (which itself may need to be reimplemented
> > as a matrix with a single column). And the next thing we know is that
> > we'll reimplement uBLAS or Eigen3 :-) I'll cross my fingers and wish we
> > won't reach that point.
> > 
> > https://stackoverflow.com/questions/1380371/what-are-the-most-widely-used-c-vector-matrix-math-linear-algebra-libraries-a
> > 
> > "It seems that many projects slowly come upon a need to do matrix math,
> > and fall into the trap of first building some vector classes and slowly
> > adding in functionality until they get caught building a half-assed
> > custom linear algebra library, and depending on it."
> 
> Oops... :)
> 
> > > + *
> > > + * \return Row \a i from the matrix, as a Span
> > > + */
> > > +
> > > +/**
> > > + * \fn Matrix::operator[](size_t i)
> > > + * \copydoc Matrix::operator[](size_t i) const
> > > + */
> > > +
> > > +/**
> > > + * \fn Matrix::Matrix<U, R, C> operator*(T d, const Matrix<U, R, C> &m)
> > > + * \brief Scalar product
> > > + * \tparam T Type of the numerical scalar value
> > > + * \tparam U Type of numerical values in the matrix
> > > + * \tparam R Number of rows in the matrix
> > > + * \tparam C Number of columns in the matrix
> > > + * \param d Scalar
> > > + * \param m Matrix
> > > + * \return Product of scalar \a d and matrix \a m
> > > + */
> > > +
> > > +/**
> > > + * \fn Matrix<T, R1, C2> operator*(const Matrix<T, R1, C1> &m1, const Matrix<T, R2, C2> &m2)
> > > + * \brief Matrix multiplication
> > > + * \tparam T Type of numerical values in the matrices
> > > + * \tparam R1 Number of rows in the first matrix
> > > + * \tparam C1 Number of columns in the first matrix
> > > + * \tparam R2 Number of rows in the second matrix
> > > + * \tparam C2 Number of columns in the second matrix
> > > + * \param m1 Multiplicand matrix
> > > + * \param m2 Multiplier matrix
> > > + * \return Matrix product of matrices \a m1 and \a m2
> > > + */
> > > +
> > > +/**
> > > + * \fn Matrix<T, R, C> operator+(const Matrix<T, R, C> &m1, const Matrix<T, R, C> &m2)
> > > + * \brief Matrix addition
> > > + * \tparam T Type of numerical values in the matrices
> > > + * \tparam R Number of rows in the matrices
> > > + * \tparam C Number of columns in the matrices
> > > + * \param m1 Summand matrix
> > > + * \param m2 Summand matrix
> > > + * \return Matrix sum of matrices \a m1 and \a m2
> > > + */
> > > +
> > > +} /* namespace ipa */
> > > +
> > > +} /* namespace libcamera */
> > > diff --git a/src/ipa/libipa/matrix.h b/src/ipa/libipa/matrix.h
> > > new file mode 100644
> > > index 000000000000..90eaea03bd14
> > > --- /dev/null
> > > +++ b/src/ipa/libipa/matrix.h
> > > @@ -0,0 +1,172 @@
> > > +/* SPDX-License-Identifier: LGPL-2.1-or-later */
> > > +/*
> > > + * Copyright (C) 2024, Paul Elder <paul.elder@ideasonboard.com>
> > > + *
> > > + * Matrix and related operations
> > > + */
> > > +#pragma once
> > > +
> > > +#include <algorithm>
> > > +#include <cmath>
> > > +#include <sstream>
> > > +#include <vector>
> > > +
> > > +#include <libcamera/base/log.h>
> > > +#include <libcamera/base/span.h>
> > > +
> > > +#include "libcamera/internal/yaml_parser.h"
> > > +
> > > +namespace libcamera {
> > > +
> > > +LOG_DECLARE_CATEGORY(Matrix)
> > > +
> > > +namespace ipa {
> > > +
> > > +#ifndef __DOXYGEN__
> > > +template<typename T, unsigned int R, unsigned int C,
> > 
> > Following the Vector class,
> > 
> > template<typename T, unsigned int Rows, unsigned int Cols,
> > 
> > > +	 std::enable_if_t<std::is_arithmetic_v<T>> * = nullptr>
> > > +#else
> > > +template<typename T, unsigned int R, unsigned int C>
> > > +#endif /* __DOXYGEN__ */
> > > +class Matrix
> > > +{
> > > +public:
> > > +	Matrix()
> > > +		: data_(R * C, static_cast<T>(false))
> > 
> > false ? That's a strange one, anything wrong with 0 ?
> > 
> > > +	{
> > > +		for (size_t i = 0; i < std::min(R, C); i++)
> > > +			(*this)[i][i] = static_cast<T>(true);
> > 
> > And 1 here.
> > 
> > > +	}
> > > +
> > > +	Matrix(const std::vector<T> &data)
> > > +	{
> > > +		ASSERT(data.size() == R * C);
> > > +
> > > +		data_.clear();
> > > +		for (const T &x : data)
> > > +			data_.push_back(x);
> > 
> > That will perform quite a few reallocations. I think you can simply
> > write
> > 
> > 		std::copy(data.begin(), data.end(), data_.begin());
> > 
> > > +	}
> > > +
> > > +	~Matrix() = default;
> > > +
> > > +	int readYaml(const libcamera::YamlObject &yaml)
> > > +	{
> > > +		if (yaml.size() != R * C) {
> > > +			LOG(Matrix, Error)
> > > +				<< "Wrong number of values in matrix: expected "
> > > +				<< R * C << ", got " << yaml.size();
> > > +			return -EINVAL;
> > > +		}
> > > +
> > > +		unsigned int i = 0;
> > > +		for (const auto &x : yaml.asList()) {
> > > +			auto value = x.get<T>();
> > > +			if (!value) {
> > > +				LOG(Matrix, Error) << "Failed to read matrix value";
> > > +				return -EINVAL;
> > > +			}
> > > +
> > > +			data_[i++] = *value;
> > > +		}
> > > +
> > > +		return 0;
> > > +	}
> > > +
> > > +	const std::string toString() const
> > > +	{
> > > +		std::stringstream out;
> > > +
> > > +		out << "Matrix { ";
> > > +		for (unsigned int i = 0; i < R; i++) {
> > > +			out << "[ ";
> > > +			for (unsigned int j = 0; j < C; j++) {
> > > +				out << (*this)[i][j];
> > > +				out << ((j + 1 < C) ? ", " : " ");
> > > +			}
> > > +			out << ((i + 1 < R) ? "], " : "]");
> > > +		}
> > > +		out << " }";
> > > +
> > > +		return out.str();
> > > +	}
> > > +
> > > +	Span<const T, C> operator[](size_t i) const
> > > +	{
> > > +		return Span<const T, C>{ &data_.data()[i * C], C };
> > > +	}
> > > +
> > > +	Span<T, C> operator[](size_t i)
> > > +	{
> > > +		return Span<T, C>{ &data_.data()[i * C], C };
> > > +	}
> > > +
> > > +private:
> > > +	std::vector<T> data_;
> > 
> > Should this be a std::array<T, Rows * Cols>, like for the Vector class ?
> > 
> > > +};
> > > +
> > > +#ifndef __DOXYGEN__
> > > +template<typename T, typename U, unsigned int R, unsigned int C,
> > > +	 std::enable_if_t<std::is_arithmetic_v<T> && std::is_arithmetic_v<U>> * = nullptr>
> > 
> > The latter is ensured by the Matrix class.
> > 
> > You need an
> > 
> > #else
> > template<typename T, typename U, unsigned int R, unsigned int C>
> > 
> > Same below.
> > 
> > > +#endif /* __DOXYGEN__ */
> > > +Matrix<U, R, C> operator*(T d, const Matrix<U, R, C> &m)
> > 
> > I'm surprised by the order of the operands, as I would have expected
> > users to write
> > 
> > 	Matrix m;
> > 	Matrix r = m * factor;
> > 
> > and not
> > 
> > 	Matrix m;
> > 	Matrix r = factor * m;
> > 
> > I recommend supporting both, you can easily implement one based on the
> > other.
> > 
> > It would also be useful to support the following with a member operator:
> > 
> > 	Matrix m;
> > 	Matrix r;
> > 	r *= factor;
> > 
> > > +{
> > > +	Matrix<U, R, C> result;
> > > +
> > > +	for (unsigned int i = 0; i < R; i++)
> > > +		for (unsigned int j = 0; j < C; j++)
> > > +			result[i][j] = d * m[i][j];
> > 
> > I think
> > 
> > 	for (unsigned int i = 0; i < R * C; i++)
> > 		result.data_[i] = d * m.data_[i];
> > 
> > would be more efficient.
> 
> Yeah but data_ is private...

You can add a friend statement to make it accessible to the out-of-class
operators. See the other Vector class I shared with you recently.

> > > +
> > > +	return result;
> > > +}
> > > +
> > > +#ifndef __DOXYGEN__
> > > +template<typename T,
> > > +	 unsigned int R1, unsigned int C1,
> > > +	 unsigned int R2, unsigned int C2,
> > > +	 std::enable_if_t<std::is_arithmetic_v<T> && C1 == R2> * = nullptr>
> > > +#endif /* __DOXYGEN__ */
> > > +Matrix<T, R1, C2> operator*(const Matrix<T, R1, C1> &m1, const Matrix<T, R2, C2> &m2)
> > 
> > This can also be simplified in a similar way as above.
> > 
> > > +{
> > > +	Matrix<T, R1, C2> result;
> > > +
> > > +	for (unsigned int i = 0; i < R1; i++) {
> > > +		for (unsigned int j = 0; j < C2; j++) {
> > > +			T sum = 0;
> > > +
> > > +			for (unsigned int k = 0; k < C1; k++)
> > > +				sum += m1[i][k] * m2[k][j];
> > > +
> > > +			result[i][j] = sum;
> > > +		}
> > > +	}
> > > +
> > > +	return result;
> > > +}
> > > +
> > > +#ifndef __DOXYGEN__
> > > +template<typename T, unsigned int R, unsigned int C,
> > > +	 std::enable_if_t<std::is_arithmetic_v<T>> * = nullptr>
> > 
> > Ensured by the Matrix type.
> > 
> > > +#endif /* __DOXYGEN__ */
> > > +Matrix<T, R, C> operator+(const Matrix<T, R, C> &m1, const Matrix<T, R, C> &m2)
> > > +{
> > > +	Matrix<T, R, C> result;
> > > +
> > > +	for (unsigned int i = 0; i < R; i++)
> > > +		for (unsigned int j = 0; j < C; j++)
> > > +			result[i][j] = m1[i][j] + m2[i][j];
> > 
> > I think
> > 
> > 	for (unsigned int i = 0; i < R * C; i++)
> > 		result.data_[i] = m1.data_[i] + m2.data_[i];
> > 
> > would be more efficient.
> > 
> > > +
> > > +	return result;
> > > +}
> > > +
> > > +} /* namespace ipa */
> > > +
> > > +#ifndef __DOXYGEN__
> > > +template<typename T, unsigned int R, unsigned int C>
> > > +std::ostream &operator<<(std::ostream &out, const ipa::Matrix<T, R, C> &m)
> > > +{
> > > +	out << m.toString();
> > > +	return out;
> > > +}
> > > +#endif /* __DOXYGEN__ */
> > > +
> > > +} /* namespace libcamera */
> > > diff --git a/src/ipa/libipa/meson.build b/src/ipa/libipa/meson.build
> > > index b663afc7d9fe..067d0d273e0a 100644
> > > --- a/src/ipa/libipa/meson.build
> > > +++ b/src/ipa/libipa/meson.build
> > > @@ -7,6 +7,7 @@ libipa_headers = files([
> > >      'exposure_mode_helper.h',
> > >      'fc_queue.h',
> > >      'histogram.h',
> > > +    'matrix.h',
> > >      'module.h',
> > >  ])
> > >  
> > > @@ -17,6 +18,7 @@ libipa_sources = files([
> > >      'exposure_mode_helper.cpp',
> > >      'fc_queue.cpp',
> > >      'histogram.cpp',
> > > +    'matrix.cpp',
> > >      'module.cpp',
> > >  ])
> > >