Geogram Version 1.8.5
A programming library of geometric algorithms
Loading...
Searching...
No Matches
matrix.h
Go to the documentation of this file.
1/*
2 * Copyright (c) 2000-2022 Inria
3 * All rights reserved.
4 *
5 * Redistribution and use in source and binary forms, with or without
6 * modification, are permitted provided that the following conditions are met:
7 *
8 * * Redistributions of source code must retain the above copyright notice,
9 * this list of conditions and the following disclaimer.
10 * * Redistributions in binary form must reproduce the above copyright notice,
11 * this list of conditions and the following disclaimer in the documentation
12 * and/or other materials provided with the distribution.
13 * * Neither the name of the ALICE Project-Team nor the names of its
14 * contributors may be used to endorse or promote products derived from this
15 * software without specific prior written permission.
16 *
17 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
18 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
19 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
20 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
21 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
22 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
23 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
24 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
25 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
26 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
27 * POSSIBILITY OF SUCH DAMAGE.
28 *
29 * Contact: Bruno Levy
30 *
31 * https://www.inria.fr/fr/bruno-levy
32 *
33 * Inria,
34 * Domaine de Voluceau,
35 * 78150 Le Chesnay - Rocquencourt
36 * FRANCE
37 *
38 */
39
40#ifndef GEOGRAM_BASIC_MATRIX
41#define GEOGRAM_BASIC_MATRIX
42
44#include <geogram/basic/vecg.h>
45
51namespace GEO {
52
53 /************************************************************************/
54
55
64 template <index_t DIM, class FT>
65 class Matrix {
66 public:
69
71 typedef FT value_type;
72
74 static const index_t dim = DIM;
75
81 inline Matrix() {
83 }
84
91 explicit Matrix(const FT* vals) {
92 for(index_t i = 0; i < DIM; i++) {
93 for(index_t j = 0; j < DIM; j++) {
94 coeff_[i][j] = *vals;
95 ++vals;
96 }
97 }
98 }
99
104 inline index_t dimension() const {
105 return DIM;
106 }
107
112 inline void load_zero() {
113 for(index_t i = 0; i < DIM; i++) {
114 for(index_t j = 0; j < DIM; j++) {
115 coeff_[i][j] = FT(0);
116 }
117 }
118 }
119
125 inline void load_identity() {
126 for(index_t i = 0; i < DIM; i++) {
127 for(index_t j = 0; j < DIM; j++) {
128 coeff_[i][j] = (i == j) ? FT(1) : FT(0);
129 }
130 }
131 }
132
138 inline bool is_identity() const {
139 for(index_t i = 0; i < DIM; i++) {
140 for(index_t j = 0; j < DIM; j++) {
141 FT rhs = ((i == j) ? FT(1) : FT(0));
142 if(coeff_[i][j] != rhs) {
143 return false;
144 }
145 }
146 }
147 return true;
148 }
149
158 inline FT& operator() (index_t i, index_t j) {
159 geo_debug_assert(i < DIM);
160 geo_debug_assert(j < DIM);
161 return coeff_[i][j];
162 }
163
172 inline const FT& operator() (index_t i, index_t j) const {
173 geo_debug_assert(i < DIM);
174 geo_debug_assert(j < DIM);
175 return coeff_[i][j];
176 }
177
185 for(index_t i = 0; i < DIM; i++) {
186 for(index_t j = 0; j < DIM; j++) {
187 coeff_[i][j] += m.coeff_[i][j];
188 }
189 }
190 return *this;
191 }
192
200 for(index_t i = 0; i < DIM; i++) {
201 for(index_t j = 0; j < DIM; j++) {
202 coeff_[i][j] -= m.coeff_[i][j];
203 }
204 }
205 return *this;
206 }
207
215 inline matrix_type& operator*= (FT val) {
216 for(index_t i = 0; i < DIM; i++) {
217 for(index_t j = 0; j < DIM; j++) {
218 coeff_[i][j] *= val;
219 }
220 }
221 return *this;
222 }
223
231 inline matrix_type& operator/= (FT val) {
232 for(index_t i = 0; i < DIM; i++) {
233 for(index_t j = 0; j < DIM; j++) {
234 coeff_[i][j] /= val;
235 }
236 }
237 return *this;
238 }
239
246 inline matrix_type operator+ (const matrix_type& m) const {
247 matrix_type result = *this;
248 result += m;
249 return result;
250 }
251
258 inline matrix_type operator- (const matrix_type& m) const {
259 matrix_type result = *this;
260 result -= m;
261 return result;
262 }
263
271 inline matrix_type operator* (FT val) const {
272 matrix_type result = *this;
273 result *= val;
274 return result;
275 }
276
284 inline matrix_type operator/ (FT val) const {
285 matrix_type result = *this;
286 result /= val;
287 return result;
288 }
289
298 matrix_type result;
299 for(index_t i = 0; i < DIM; i++) {
300 for(index_t j = 0; j < DIM; j++) {
301 result.coeff_[i][j] = FT(0);
302 for(index_t k = 0; k < DIM; k++) {
303 result.coeff_[i][j] += coeff_[i][k] * m.coeff_[k][j];
304 }
305 }
306 }
307 return result;
308 }
309
316 matrix_type result;
317 bool invertible = compute_inverse(result);
318 geo_assert(invertible);
319 return result;
320 }
321
322
330 bool compute_inverse(matrix_type& result) const {
331 FT val=FT(0.0), val2=FT(0.0);
332 matrix_type tmp = (*this);
333
334 result.load_identity();
335
336 for(index_t i = 0; i != DIM; i++) {
337 val = tmp(i, i); /* find pivot */
338 index_t ind = i;
339 for(index_t j = i + 1; j != DIM; j++) {
340 if(fabs(tmp(j, i)) > fabs(val)) {
341 ind = j;
342 val = tmp(j, i);
343 }
344 }
345
346 if(ind != i) {
347 for(index_t j = 0; j != DIM; j++) {
348 val2 = result(i, j);
349 result(i, j) = result(ind, j);
350 result(ind, j) = val2; /* swap columns */
351 val2 = tmp(i, j);
352 tmp(i, j) = tmp(ind, j);
353 tmp(ind, j) = val2;
354 }
355 }
356
357 if(val == 0.0) {
358 return false;
359 }
360
361 for(index_t j = 0; j != DIM; j++) {
362 tmp(i, j) /= val;
363 result(i, j) /= val;
364 }
365
366 for(index_t j = 0; j != DIM; j++) {
367 if(j == i) {
368 continue; /* eliminate column */
369 }
370 val = tmp(j, i);
371 for(index_t k = 0; k != DIM; k++) {
372 tmp(j, k) -= tmp(i, k) * val;
373 result(j, k) -= result(i, k) * val;
374 }
375 }
376 }
377
378 return true;
379 }
380
386 matrix_type result;
387 for(index_t i = 0; i < DIM; i++) {
388 for(index_t j = 0; j < DIM; j++) {
389 result(i, j) = (* this)(j, i);
390 }
391 }
392 return result;
393 }
394
401 inline const FT* data() const {
402 return &(coeff_[0][0]);
403 }
404
411 inline FT* data() {
412 return &(coeff_[0][0]);
413 }
414
422 void get_lower_triangle(FT* store) const {
423 for(index_t i = 0; i < DIM; i++) {
424 for(index_t j = 0; j <= i; j++) {
425 *store++ = coeff_[i][j];
426 }
427 }
428 }
429
430 private:
431 FT coeff_[DIM][DIM];
432 };
433
434 /************************************************************************/
435
445 template <index_t DIM, class FT>
446 inline std::ostream& operator<< (
447 std::ostream& output, const Matrix<DIM, FT>& m
448 ) {
449 const char* sep = "";
450 for(index_t i = 0; i < DIM; i++) {
451 for(index_t j = 0; j < DIM; j++) {
452 output << sep << m(i, j);
453 sep = " ";
454 }
455 }
456 return output;
457 }
458
468 template <index_t DIM, class FT>
469 inline std::istream& operator>> (
470 std::istream& input, Matrix<DIM, FT>& m
471 ) {
472 for(index_t i = 0; i < DIM; i++) {
473 for(index_t j = 0; j < DIM; j++) {
474 input >> m(i, j);
475 }
476 }
477 return input;
478 }
479
480 /************************************************************************/
481
495 template <index_t DIM, class FT> inline
496 void mult(const Matrix<DIM, FT>& M, const FT* x, FT* y) {
497 for(index_t i = 0; i < DIM; i++) {
498 y[i] = 0;
499 for(index_t j = 0; j < DIM; j++) {
500 y[i] += M(i, j) * x[j];
501 }
502 }
503 }
504
505 /************************************************************************/
506
515 template <index_t DIM, class FT> inline
517 const Matrix<DIM, FT>& M, const vecng<DIM,FT>& x
518 ) {
520 for(index_t i = 0; i < DIM; i++) {
521 y[i] = 0;
522 for(index_t j = 0; j < DIM; j++) {
523 y[i] += M(i, j) * x[j];
524 }
525 }
526 return y;
527 }
528
537 template <index_t DIM, class FT> inline
539 const Matrix<DIM, FT>& M, const vecng<DIM,FT>& x
540 ) {
542 for(index_t i = 0; i < DIM; i++) {
543 y[i] = 0;
544 for(index_t j = 0; j < DIM; j++) {
545 y[i] += M(i, j) * x[j];
546 }
547 }
548 return y;
549 }
550
551 /************************************************************************/
552
553}
554
555#endif
556
#define geo_assert(x)
Verifies that a condition is met.
Definition assert.h:149
#define geo_debug_assert(x)
Verifies that a condition is met.
Definition assert.h:195
Common include file, providing basic definitions. Should be included before anything else by all head...
A matrix type.
Definition matrix.h:65
matrix_type transpose() const
Computes the transposed matrix.
Definition matrix.h:385
matrix_type & operator*=(FT val)
Multiplies by a scalar in place.
Definition matrix.h:215
void get_lower_triangle(FT *store) const
Gets the lower triangle of the matrix.
Definition matrix.h:422
matrix_type & operator-=(const matrix_type &m)
Subtracts a matrix in place.
Definition matrix.h:199
Matrix(const FT *vals)
Constructs a matrix from an array of values.
Definition matrix.h:91
FT value_type
Definition matrix.h:71
index_t dimension() const
Gets the matrix dimension.
Definition matrix.h:104
Matrix< DIM, FT > matrix_type
Definition matrix.h:68
matrix_type operator+(const matrix_type &m) const
Adds 2 matrices.
Definition matrix.h:246
const FT * data() const
Gets non-modifiable matrix data.
Definition matrix.h:401
static const index_t dim
Definition matrix.h:74
FT & operator()(index_t i, index_t j)
Gets a modifiable element.
Definition matrix.h:158
matrix_type operator*(FT val) const
Multiplies a matrix by a scalar.
Definition matrix.h:271
void load_zero()
Clears the matrix.
Definition matrix.h:112
matrix_type & operator/=(FT val)
Divides by a scalar in place.
Definition matrix.h:231
bool is_identity() const
Tests whether a matrix is the identity matrix.
Definition matrix.h:138
matrix_type operator-(const matrix_type &m) const
Subtracts 2 matrices.
Definition matrix.h:258
void mult(const Matrix< DIM, FT > &M, const FT *x, FT *y)
Multiplies a matrix by a vector.
Definition matrix.h:496
FT * data()
Gets modifiable matrix data.
Definition matrix.h:411
matrix_type inverse() const
Computes the inverse matrix.
Definition matrix.h:315
void load_identity()
Sets the matrix to identity.
Definition matrix.h:125
Matrix()
Default constructor.
Definition matrix.h:81
bool compute_inverse(matrix_type &result) const
Computes the inverse matrix.
Definition matrix.h:330
matrix_type & operator+=(const matrix_type &m)
Adds a matrix in place.
Definition matrix.h:184
matrix_type operator/(FT val) const
Divides a matrix by a scalar.
Definition matrix.h:284
Generic maths vector.
Definition vecg.h:69
Global Vorpaline namespace.
Definition algorithm.h:64
std::istream & operator>>(std::istream &in, Quaternion &q)
Reads a Quaternion from a stream.
Definition quaternion.h:225
std::ostream & operator<<(std::ostream &out, const Quaternion &q)
Writes a Quaternion to a stream.
Definition quaternion.h:213
geo_index_t index_t
The type for storing and manipulating indices.
Definition numeric.h:287
vecng< DIM, FT > mult(const Matrix< DIM, FT > &M, const vecng< DIM, FT > &x)
Computes a matrix vector product.
Definition matrix.h:538
vecng< DIM, FT > operator*(const Matrix< DIM, FT > &M, const vecng< DIM, FT > &x)
Computes a matrix vector product.
Definition matrix.h:516
Generic implementation of geometric vectors.