tcuMatrix.hpp revision cb82ed72dcbbfd8a6d07736c3259605227bc984f 
1#ifndef _TCUMATRIX_HPP
2#define _TCUMATRIX_HPP
3/*-------------------------------------------------------------------------
4 * drawElements Quality Program Tester Core
5 * ----------------------------------------
6 *
7 * Copyright 2014 The Android Open Source Project
8 *
9 * Licensed under the Apache License, Version 2.0 (the "License");
10 * you may not use this file except in compliance with the License.
11 * You may obtain a copy of the License at
12 *
13 *      http://www.apache.org/licenses/LICENSE-2.0
14 *
15 * Unless required by applicable law or agreed to in writing, software
16 * distributed under the License is distributed on an "AS IS" BASIS,
17 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
18 * See the License for the specific language governing permissions and
19 * limitations under the License.
20 *
21 *//*!
22 * \file
23 * \brief 	Templatized matrix class.
24 *//*--------------------------------------------------------------------*/
25
26#include "tcuDefs.hpp"
27#include "tcuVector.hpp"
28#include "tcuArray.hpp"
29
30namespace tcu
31{
32
33// Templated matrix class.
34template <typename T, int Rows, int Cols>
35class Matrix
36{
37public:
38	typedef	Vector<T, Rows>			Element;
39	typedef	T						Scalar;
40
41	enum
42	{
43		SIZE = Cols,
44		ROWS = Rows,
45		COLS = Cols,
46	};
47
48									Matrix				(void);
49	explicit						Matrix				(const T& src);
50	explicit						Matrix				(const T src[Rows*Cols]);
51									Matrix				(const Vector<T, Rows>& src);
52									Matrix				(const Matrix<T, Rows, Cols>& src);
53									~Matrix				(void);
54
55	Matrix<T, Rows, Cols>&			operator=			(const Matrix<T, Rows, Cols>& src);
56	Matrix<T, Rows, Cols>&			operator*=			(const Matrix<T, Rows, Cols>& src);
57
58	void							setRow				(int rowNdx, const Vector<T, Cols>& vec);
59	void							setColumn			(int colNdx, const Vector<T, Rows>& vec);
60
61	Vector<T, Cols>					getRow				(int ndx) const;
62	Vector<T, Rows>&				getColumn			(int ndx);
63	const Vector<T, Rows>&			getColumn			(int ndx) const;
64
65	Vector<T, Rows>&				operator[]			(int ndx)					{ return getColumn(ndx);	}
66	const Vector<T, Rows>&			operator[]			(int ndx) const				{ return getColumn(ndx);	}
67
68	inline const T&					operator()			(int row, int col) const	{ return m_data[col][row];	}
69	inline T&						operator()			(int row, int col)			{ return m_data[col][row];	}
70
71	Array<T, Rows*Cols>				getRowMajorData		(void) const;
72	Array<T, Rows*Cols>				getColumnMajorData	(void) const;
73
74private:
75	Vector<Vector<T, Rows>, Cols>	m_data;
76};
77
78// Operators.
79
80// Mat * Mat.
81template <typename T, int Rows0, int Cols0, int Rows1, int Cols1>
82Matrix<T, Rows0, Cols1> operator* (const Matrix<T, Rows0, Cols0>& a, const Matrix<T, Rows1, Cols1>& b);
83
84// Mat * Vec (column vector).
85template <typename T, int Rows, int Cols>
86Vector<T, Rows> operator* (const Matrix<T, Rows, Cols>& mtx, const Vector<T, Cols>& vec);
87
88// Vec * Mat (row vector).
89template <typename T, int Rows, int Cols>
90Vector<T, Cols> operator* (const Vector<T, Rows>& vec, const Matrix<T, Rows, Cols>& mtx);
91
92// Further operations
93
94template <typename T, int Size>
95struct SquareMatrixOps
96{
97	static T						doDeterminant	(const Matrix<T, Size, Size>& mat);
98	static Matrix<T, Size, Size>	doInverse		(const Matrix<T, Size, Size>& mat);
99};
100
101template <typename T>
102struct SquareMatrixOps<T, 2>
103{
104	static T						doDeterminant	(const Matrix<T, 2, 2>& mat);
105	static Matrix<T, 2, 2>			doInverse		(const Matrix<T, 2, 2>& mat);
106};
107
108template <typename T>
109struct SquareMatrixOps<T, 3>
110{
111	static T						doDeterminant	(const Matrix<T, 3, 3>& mat);
112	static Matrix<T, 3, 3>			doInverse		(const Matrix<T, 3, 3>& mat);
113};
114
115template <typename T>
116struct SquareMatrixOps<T, 4>
117{
118	static T						doDeterminant	(const Matrix<T, 4, 4>& mat);
119	static Matrix<T, 4, 4>			doInverse		(const Matrix<T, 4, 4>& mat);
120};
121
122namespace matrix
123{
124
125template <typename T, int Size>
126T determinant (const Matrix<T, Size, Size>& mat)
127{
128	return SquareMatrixOps<T, Size>::doDeterminant(mat);
129}
130
131template <typename T, int Size>
132Matrix<T, Size, Size> inverse (const Matrix<T, Size, Size>& mat)
133{
134	return SquareMatrixOps<T, Size>::doInverse(mat);
135}
136
137} // matrix
138
139// Template implementations.
140
141template <typename T>
142T SquareMatrixOps<T, 2>::doDeterminant (const Matrix<T, 2, 2>& mat)
143{
144	return mat(0,0) * mat(1,1) - mat(1,0) * mat(0,1);
145}
146
147template <typename T>
148T SquareMatrixOps<T, 3>::doDeterminant (const Matrix<T, 3, 3>& mat)
149{
150	return	+ mat(0,0) * mat(1,1) * mat(2,2)
151			+ mat(0,1) * mat(1,2) * mat(2,0)
152			+ mat(0,2) * mat(1,0) * mat(2,1)
153			- mat(0,0) * mat(1,2) * mat(2,1)
154			- mat(0,1) * mat(1,0) * mat(2,2)
155			- mat(0,2) * mat(1,1) * mat(2,0);
156}
157
158template <typename T>
159T SquareMatrixOps<T, 4>::doDeterminant (const Matrix<T, 4, 4>& mat)
160{
161	using matrix::determinant;
162
163	const T minorMatrices[4][3*3] =
164	{
165		{
166			mat(1,1),	mat(2,1),	mat(3,1),
167			mat(1,2),	mat(2,2),	mat(3,2),
168			mat(1,3),	mat(2,3),	mat(3,3),
169		},
170		{
171			mat(1,0),	mat(2,0),	mat(3,0),
172			mat(1,2),	mat(2,2),	mat(3,2),
173			mat(1,3),	mat(2,3),	mat(3,3),
174		},
175		{
176			mat(1,0),	mat(2,0),	mat(3,0),
177			mat(1,1),	mat(2,1),	mat(3,1),
178			mat(1,3),	mat(2,3),	mat(3,3),
179		},
180		{
181			mat(1,0),	mat(2,0),	mat(3,0),
182			mat(1,1),	mat(2,1),	mat(3,1),
183			mat(1,2),	mat(2,2),	mat(3,2),
184		}
185	};
186
187	return	+ mat(0,0) * determinant(Matrix<T, 3, 3>(minorMatrices[0]))
188			- mat(0,1) * determinant(Matrix<T, 3, 3>(minorMatrices[1]))
189			+ mat(0,2) * determinant(Matrix<T, 3, 3>(minorMatrices[2]))
190			- mat(0,3) * determinant(Matrix<T, 3, 3>(minorMatrices[3]));
191}
192
193template <typename T>
194Matrix<T, 2, 2> SquareMatrixOps<T, 2>::doInverse (const Matrix<T, 2, 2>& mat)
195{
196	using matrix::determinant;
197
198	const T			det		= determinant(mat);
199	Matrix<T, 2, 2>	retVal;
200
201	retVal(0, 0) =  mat(1, 1) / det;
202	retVal(0, 1) = -mat(0, 1) / det;
203	retVal(1, 0) = -mat(1, 0) / det;
204	retVal(1, 1) =  mat(0, 0) / det;
205
206	return retVal;
207}
208
209template <typename T>
210Matrix<T, 3, 3> SquareMatrixOps<T, 3>::doInverse (const Matrix<T, 3, 3>& mat)
211{
212	// Blockwise inversion
213	using matrix::inverse;
214
215	const T areaA[2*2] =
216	{
217		mat(0,0),	mat(0,1),
218		mat(1,0),	mat(1,1)
219	};
220	const T areaB[2] =
221	{
222		mat(0,2),
223		mat(1,2),
224	};
225	const T areaC[2] =
226	{
227		mat(2,0),	mat(2,1),
228	};
229	const T areaD[1] =
230	{
231		mat(2,2)
232	};
233	const T nullField[4] = { T(0.0f) };
234
235	const Matrix<T, 2, 2>	invA = inverse(Matrix<T, 2, 2>(areaA));
236	const Matrix<T, 2, 1>	matB =         Matrix<T, 2, 1>(areaB);
237	const Matrix<T, 1, 2>	matC =         Matrix<T, 1, 2>(areaC);
238	const Matrix<T, 1, 1>	matD =         Matrix<T, 1, 1>(areaD);
239
240	const T					schurComplement = T(1.0f) / (matD - matC*invA*matB)(0,0);
241	const Matrix<T, 2, 2>	zeroMat         = Matrix<T, 2, 2>(nullField);
242
243	const Matrix<T, 2, 2>	blockA = invA + invA*matB*schurComplement*matC*invA;
244	const Matrix<T, 2, 1>	blockB = (zeroMat-invA)*matB*schurComplement;
245	const Matrix<T, 1, 2>	blockC = matC*invA*(-schurComplement);
246	const T					blockD = schurComplement;
247
248	const T result[3*3] =
249	{
250		blockA(0,0),	blockA(0,1),	blockB(0,0),
251		blockA(1,0),	blockA(1,1),	blockB(1,0),
252		blockC(0,0),	blockC(0,1),	blockD,
253	};
254
255	return Matrix<T, 3, 3>(result);
256}
257
258template <typename T>
259Matrix<T, 4, 4> SquareMatrixOps<T, 4>::doInverse (const Matrix<T, 4, 4>& mat)
260{
261	// Blockwise inversion
262	using matrix::inverse;
263
264	const T areaA[2*2] =
265	{
266		mat(0,0),	mat(0,1),
267		mat(1,0),	mat(1,1)
268	};
269	const T areaB[2*2] =
270	{
271		mat(0,2),	mat(0,3),
272		mat(1,2),	mat(1,3)
273	};
274	const T areaC[2*2] =
275	{
276		mat(2,0),	mat(2,1),
277		mat(3,0),	mat(3,1)
278	};
279	const T areaD[2*2] =
280	{
281		mat(2,2),	mat(2,3),
282		mat(3,2),	mat(3,3)
283	};
284	const T nullField[4] = { T(0.0f) };
285
286	const Matrix<T, 2, 2> invA = inverse(Matrix<T, 2, 2>(areaA));
287	const Matrix<T, 2, 2> matB =         Matrix<T, 2, 2>(areaB);
288	const Matrix<T, 2, 2> matC =         Matrix<T, 2, 2>(areaC);
289	const Matrix<T, 2, 2> matD =         Matrix<T, 2, 2>(areaD);
290
291	const Matrix<T, 2, 2> schurComplement = inverse(matD - matC*invA*matB);
292	const Matrix<T, 2, 2> zeroMat         = Matrix<T, 2, 2>(nullField);
293
294	const Matrix<T, 2, 2> blockA = invA + invA*matB*schurComplement*matC*invA;
295	const Matrix<T, 2, 2> blockB = (zeroMat-invA)*matB*schurComplement;
296	const Matrix<T, 2, 2> blockC = (zeroMat-schurComplement)*matC*invA;
297	const Matrix<T, 2, 2> blockD = schurComplement;
298
299	const T result[4*4] =
300	{
301		blockA(0,0),	blockA(0,1),	blockB(0,0),	blockB(0,1),
302		blockA(1,0),	blockA(1,1),	blockB(1,0),	blockB(1,1),
303		blockC(0,0),	blockC(0,1),	blockD(0,0),	blockD(0,1),
304		blockC(1,0),	blockC(1,1),	blockD(1,0),	blockD(1,1),
305	};
306
307	return Matrix<T, 4, 4>(result);
308}
309
310// Initialize to identity.
311template <typename T, int Rows, int Cols>
312Matrix<T, Rows, Cols>::Matrix (void)
313{
314	for (int row = 0; row < Rows; row++)
315		for (int col = 0; col < Cols; col++)
316			(*this)(row, col) = (row == col) ? T(1) : T(0);
317}
318
319// Initialize to diagonal matrix.
320template <typename T, int Rows, int Cols>
321Matrix<T, Rows, Cols>::Matrix (const T& src)
322{
323	for (int row = 0; row < Rows; row++)
324		for (int col = 0; col < Cols; col++)
325			(*this)(row, col) = (row == col) ? src : T(0);
326}
327
328// Initialize from data array.
329template <typename T, int Rows, int Cols>
330Matrix<T, Rows, Cols>::Matrix (const T src[Rows*Cols])
331{
332	for (int row = 0; row < Rows; row++)
333		for (int col = 0; col < Cols; col++)
334			(*this)(row, col) = src[row*Cols + col];
335}
336
337// Initialize to diagonal matrix.
338template <typename T, int Rows, int Cols>
339Matrix<T, Rows, Cols>::Matrix (const Vector<T, Rows>& src)
340{
341	DE_STATIC_ASSERT(Rows == Cols);
342	for (int row = 0; row < Rows; row++)
343		for (int col = 0; col < Cols; col++)
344			(*this)(row, col) = (row == col) ? src.m_data[row] : T(0);
345}
346
347// Copy constructor.
348template <typename T, int Rows, int Cols>
349Matrix<T, Rows, Cols>::Matrix (const Matrix<T, Rows, Cols>& src)
350{
351	*this = src;
352}
353
354// Destructor.
355template <typename T, int Rows, int Cols>
356Matrix<T, Rows, Cols>::~Matrix (void)
357{
358}
359
360// Assignment operator.
361template <typename T, int Rows, int Cols>
362Matrix<T, Rows, Cols>& Matrix<T, Rows, Cols>::operator= (const Matrix<T, Rows, Cols>& src)
363{
364	for (int row = 0; row < Rows; row++)
365		for (int col = 0; col < Cols; col++)
366			(*this)(row, col) = src(row, col);
367	return *this;
368}
369
370// Multipy and assign op
371template <typename T, int Rows, int Cols>
372Matrix<T, Rows, Cols>& Matrix<T, Rows, Cols>::operator*= (const Matrix<T, Rows, Cols>& src)
373{
374	*this = *this * src;
375	return *this;
376}
377
378template <typename T, int Rows, int Cols>
379void Matrix<T, Rows, Cols>::setRow (int rowNdx, const Vector<T, Cols>& vec)
380{
381	for (int col = 0; col < Cols; col++)
382		(*this)(rowNdx, col) = vec.m_data[col];
383}
384
385template <typename T, int Rows, int Cols>
386void Matrix<T, Rows, Cols>::setColumn (int colNdx, const Vector<T, Rows>& vec)
387{
388	m_data[colNdx] = vec;
389}
390
391template <typename T, int Rows, int Cols>
392Vector<T, Cols> Matrix<T, Rows, Cols>::getRow (int rowNdx) const
393{
394	Vector<T, Cols> res;
395	for (int col = 0; col < Cols; col++)
396		res[col] = (*this)(rowNdx, col);
397	return res;
398}
399
400template <typename T, int Rows, int Cols>
401Vector<T, Rows>& Matrix<T, Rows, Cols>::getColumn (int colNdx)
402{
403	return m_data[colNdx];
404}
405
406template <typename T, int Rows, int Cols>
407const Vector<T, Rows>& Matrix<T, Rows, Cols>::getColumn (int colNdx) const
408{
409	return m_data[colNdx];
410}
411
412template <typename T, int Rows, int Cols>
413Array<T, Rows*Cols> Matrix<T, Rows, Cols>::getColumnMajorData (void) const
414{
415	Array<T, Rows*Cols> a;
416	T* dst = a.getPtr();
417	for (int col = 0; col < Cols; col++)
418		for (int row = 0; row < Rows; row++)
419			*dst++ = (*this)(row, col);
420	return a;
421}
422
423template <typename T, int Rows, int Cols>
424Array<T, Rows*Cols> Matrix<T, Rows, Cols>::getRowMajorData (void) const
425{
426	Array<T, Rows*Cols> a;
427	T* dst = a.getPtr();
428	for (int row = 0; row < Rows; row++)
429		for (int col = 0; col < Cols; col++)
430			*dst++ = (*this)(row, col);
431	return a;
432}
433
434// Multiplication of two matrices.
435template <typename T, int Rows0, int Cols0, int Rows1, int Cols1>
436Matrix<T, Rows0, Cols1> operator* (const Matrix<T, Rows0, Cols0>& a, const Matrix<T, Rows1, Cols1>& b)
437{
438	DE_STATIC_ASSERT(Cols0 == Rows1);
439	Matrix<T, Rows0, Cols1> res;
440	for (int row = 0; row < Rows0; row++)
441	{
442		for (int col = 0; col < Cols1; col++)
443		{
444			T v = T(0);
445			for (int ndx = 0; ndx < Cols0; ndx++)
446				v += a(row,ndx) * b(ndx,col);
447			res(row,col) = v;
448		}
449	}
450	return res;
451}
452
453// Multiply of matrix with column vector.
454template <typename T, int Rows, int Cols>
455Vector<T, Rows> operator* (const Matrix<T, Rows, Cols>& mtx, const Vector<T, Cols>& vec)
456{
457	Vector<T, Rows> res;
458	for (int row = 0; row < Rows; row++)
459	{
460		T v = T(0);
461		for (int col = 0; col < Cols; col++)
462			v += mtx(row,col) * vec.m_data[col];
463		res.m_data[row] = v;
464	}
465	return res;
466}
467
468// Multiply of matrix with row vector.
469template <typename T, int Rows, int Cols>
470Vector<T, Cols> operator* (const Vector<T, Rows>& vec, const Matrix<T, Rows, Cols>& mtx)
471{
472	Vector<T, Cols> res;
473	for (int col = 0; col < Cols; col++)
474	{
475		T v = T(0);
476		for (int row = 0; row < Rows; row++)
477			v += mtx(row,col) * vec.m_data[row];
478		res.m_data[col] = v;
479	}
480	return res;
481}
482
483// Common typedefs.
484typedef Matrix<float, 2, 2>		Matrix2f;
485typedef Matrix<float, 3, 3>		Matrix3f;
486typedef Matrix<float, 4, 4>		Matrix4f;
487typedef Matrix<double, 2, 2>	Matrix2d;
488typedef Matrix<double, 3, 3>	Matrix3d;
489typedef Matrix<double, 4, 4>	Matrix4d;
490
491// GLSL-style naming \note CxR.
492typedef Matrix2f			Mat2;
493typedef Matrix<float, 3, 2>	Mat2x3;
494typedef Matrix<float, 4, 2>	Mat2x4;
495typedef Matrix<float, 2, 3>	Mat3x2;
496typedef Matrix3f			Mat3;
497typedef Matrix<float, 4, 3>	Mat3x4;
498typedef Matrix<float, 2, 4>	Mat4x2;
499typedef Matrix<float, 3, 4>	Mat4x3;
500typedef Matrix4f			Mat4;
501
502// Matrix-scalar operators.
503
504template <typename T, int Rows, int Cols>
505Matrix<T, Rows, Cols> operator+ (const Matrix<T, Rows, Cols>& mtx, T scalar)
506{
507	Matrix<T, Rows, Cols> res;
508	for (int col = 0; col < Cols; col++)
509		for (int row = 0; row < Rows; row++)
510			res(row, col) = mtx(row, col) + scalar;
511	return res;
512}
513
514template <typename T, int Rows, int Cols>
515Matrix<T, Rows, Cols> operator- (const Matrix<T, Rows, Cols>& mtx, T scalar)
516{
517	Matrix<T, Rows, Cols> res;
518	for (int col = 0; col < Cols; col++)
519		for (int row = 0; row < Rows; row++)
520			res(row, col) = mtx(row, col) - scalar;
521	return res;
522}
523
524template <typename T, int Rows, int Cols>
525Matrix<T, Rows, Cols> operator* (const Matrix<T, Rows, Cols>& mtx, T scalar)
526{
527	Matrix<T, Rows, Cols> res;
528	for (int col = 0; col < Cols; col++)
529		for (int row = 0; row < Rows; row++)
530			res(row, col) = mtx(row, col) * scalar;
531	return res;
532}
533
534template <typename T, int Rows, int Cols>
535Matrix<T, Rows, Cols> operator/ (const Matrix<T, Rows, Cols>& mtx, T scalar)
536{
537	Matrix<T, Rows, Cols> res;
538	for (int col = 0; col < Cols; col++)
539		for (int row = 0; row < Rows; row++)
540			res(row, col) = mtx(row, col) / scalar;
541	return res;
542}
543
544// Matrix-matrix component-wise operators.
545
546template <typename T, int Rows, int Cols>
547Matrix<T, Rows, Cols> operator+ (const Matrix<T, Rows, Cols>& a, const Matrix<T, Rows, Cols>& b)
548{
549	Matrix<T, Rows, Cols> res;
550	for (int col = 0; col < Cols; col++)
551		for (int row = 0; row < Rows; row++)
552			res(row, col) = a(row, col) + b(row, col);
553	return res;
554}
555
556template <typename T, int Rows, int Cols>
557Matrix<T, Rows, Cols> operator- (const Matrix<T, Rows, Cols>& a, const Matrix<T, Rows, Cols>& b)
558{
559	Matrix<T, Rows, Cols> res;
560	for (int col = 0; col < Cols; col++)
561		for (int row = 0; row < Rows; row++)
562			res(row, col) = a(row, col) - b(row, col);
563	return res;
564}
565
566template <typename T, int Rows, int Cols>
567Matrix<T, Rows, Cols> operator/ (const Matrix<T, Rows, Cols>& a, const Matrix<T, Rows, Cols>& b)
568{
569	Matrix<T, Rows, Cols> res;
570	for (int col = 0; col < Cols; col++)
571		for (int row = 0; row < Rows; row++)
572			res(row, col) = a(row, col) / b(row, col);
573	return res;
574}
575
576} // tcu
577
578#endif // _TCUMATRIX_HPP
579