Tulip  4.1.0
Better Visualization Through Research
 All Classes Files Functions Variables Enumerations Enumerator Properties Groups Pages
Matrix.cxx
1 /*
2  *
3  * This file is part of Tulip (www.tulip-software.org)
4  *
5  * Authors: David Auber and the Tulip development Team
6  * from LaBRI, University of Bordeaux 1 and Inria Bordeaux - Sud Ouest
7  *
8  * Tulip is free software; you can redistribute it and/or modify
9  * it under the terms of the GNU Lesser General Public License
10  * as published by the Free Software Foundation, either version 3
11  * of the License, or (at your option) any later version.
12  *
13  * Tulip is distributed in the hope that it will be useful,
14  * but WITHOUT ANY WARRANTY; without even the implied warranty of
15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
16  * See the GNU General Public License for more details.
17  *
18  */
19 #include <cmath>
20 
21 //===================================================================
22 //Specialisation
23 namespace tlp {
24 template<typename Obj>
25 class Matrix<Obj,1>:public Array<Vector<Obj,1>,1> {
26 public:
27  Obj determinant() {
28  return (*this)[0][0];
29  }
30  // Matrix<Obj,1>& fill(Obj obj) {return *this;}
31  Matrix<Obj,1>& inverse() {
32  (*this)[0][0] = 1.0 / (*this)[0][0];
33  return *this;
34  }
35  Matrix<Obj,1>& transpose() {
36  return *this;
37  }
38  Matrix<Obj,1>& operator*=(const Matrix<Obj,1> &mat) {
39  (*this)[0][0] *= mat[0][0];
40  return *this;
41  }
42  // Matrix<Obj,1>& operator/=(const Obj &obj){return *this;}
43  Matrix<Obj,1> cofactor() {
44  return *this;
45  }
46 };
47 }
48 
49 
50 //===================================================================
51 template<typename Obj,unsigned int SIZE>
52 MATRIX::Matrix(const std::vector<std::vector<Obj> > &covarianceMatrix) {
53  for(unsigned int i=0; i < SIZE; i++)
54  for(unsigned int j=0; j < SIZE; j++)
55  (*this)[i][j] = covarianceMatrix[i][j] / (sqrt(covarianceMatrix[i][i] * covarianceMatrix[j][j]));
56 }
57 
58 //===================================================================
59 template<typename Obj,unsigned int SIZE>
60 MATRIX & MATRIX::fill(Obj obj) {
61  for (unsigned int i=0; i<SIZE; ++i)
62  (*this)[i].fill(obj);
63 
64  return (*this);
65 }
66 //======================================================
67 template<typename Obj,unsigned int SIZE>
68 MATRIX& MATRIX::operator+=(const MATRIX &mat) {
69  for (unsigned int i=0; i<SIZE; ++i)
70  (*this)[i] += mat[i];
71 
72  return (*this);
73 }
74 //======================================================
75 template<typename Obj,unsigned int SIZE>
76 MATRIX& MATRIX::operator-=(const MATRIX &mat) {
77  for (unsigned int i=0; i<SIZE; ++i)
78  (*this)[i] -= mat[i];
79 
80  return (*this);
81 }
82 //======================================================
83 template<typename Obj,unsigned int SIZE>
84 bool MATRIX::operator==(const MATRIX &mat) const {
85  for (unsigned int i=0; i<SIZE; ++i) {
86  if (((*this)[i] != mat[i]))
87  return false;
88  }
89 
90  return true;
91 }
92 //======================================================
93 template<typename Obj,unsigned int SIZE>
94 bool MATRIX::operator!=(const MATRIX &mat) const {
95  for (unsigned int i=0; i<SIZE; ++i) {
96  if (((*this)[i] != mat[i]))
97  return true;
98  }
99 
100  return false;
101 }
102 //===================================================================
103 template<typename Obj,unsigned int SIZE>
104 MATRIX & MATRIX::operator*=(const MATRIX &mat) {
105  (*this) = (*this) * mat;
106  return (*this);
107 }
108 //=====================================================================================
109 template<typename Obj,unsigned int SIZE>
110 MATRIX & MATRIX::operator*=(const Obj &obj) {
111  for (unsigned int i=0; i<SIZE; ++i)
112  (*this)[i] *= obj;
113 
114  return (*this);
115 }
116 //=====================================================================================
117 template<typename Obj,unsigned int SIZE>
118 MATRIX & MATRIX::operator/=(const MATRIX &mat) {
119  MATRIX tmpMat(mat);
120  tmpMat.inverse();
121  (*this) *= tmpMat;
122  return (*this);
123 }
124 //=====================================================================================
125 template<typename Obj,unsigned int SIZE>
126 MATRIX & MATRIX::operator/=(const Obj &obj) {
127  for (unsigned int i=0; i<SIZE; ++i)
128  (*this)[i] /= obj;
129 
130  return (*this);
131 }
132 //=====================================================================================
133 template<typename Obj,unsigned int SIZE>
134 Obj MATRIX::determinant() const {
135  switch (SIZE) {
136  case 2:
137  return (*this)[0][0] * (*this)[1][1] - (*this)[1][0] * (*this)[0][1];
138  break;
139 
140  case 3:
141  return (*this)[0][0] * ((*this)[1][1]*(*this)[2][2] - (*this)[1][2] * (*this)[2][1])
142  - (*this)[0][1] * ((*this)[1][0]*(*this)[2][2] - (*this)[1][2] * (*this)[2][0])
143  + (*this)[0][2] * ((*this)[1][0]*(*this)[2][1] - (*this)[1][1] * (*this)[2][0]) ;
144  break;
145 
146  default:
147  int j2;
148  Obj det = 0;
149 
150  for (unsigned int j1=0; j1<SIZE; ++j1) {
151  tlp::Matrix<Obj, SIZE - 1> m;
152 
153  for (unsigned int i=1; i<SIZE; i++) {
154  j2 = 0;
155 
156  for (unsigned int j=0; j<SIZE; ++j) {
157  if (j == j1)
158  continue;
159 
160  m[i-1][j2] = (*this)[i][j];
161  ++j2;
162  }
163  }
164 
165  if (j1 & 1)
166  det += (*this)[0][j1] * m.determinant();
167  else
168  det -= (*this)[0][j1] * m.determinant();
169  }
170 
171  return(det);
172  }
173 }
174 //=====================================================================================
175 template<typename Obj,unsigned int SIZE>
176 MATRIX MATRIX::cofactor() const {
177  MATRIX result;
178 
179  switch (SIZE) {
180  case 2:
181  (result)[0][0] = (*this)[1][1];
182  (result)[0][1] = - (*this)[1][0];
183  (result)[1][0] = - (*this)[0][1];
184  (result)[1][1] = (*this)[0][0];
185  break;
186 
187  case 3:
188  (result)[0][0] = (*this)[1][1]*(*this)[2][2] - (*this)[1][2]*(*this)[2][1];
189  (result)[0][1] = - ((*this)[1][0]*(*this)[2][2] - (*this)[2][0]*(*this)[1][2]);
190  (result)[0][2] = (*this)[1][0]*(*this)[2][1] - (*this)[1][1]*(*this)[2][0];
191  (result)[1][0] = - ((*this)[0][1]*(*this)[2][2] - (*this)[0][2]*(*this)[2][1]);
192  (result)[1][1] = (*this)[0][0]*(*this)[2][2] - (*this)[0][2]*(*this)[2][0];
193  (result)[1][2] = - ((*this)[0][0]*(*this)[2][1] - (*this)[0][1]*(*this)[2][0]);
194  (result)[2][0] = (*this)[0][1]*(*this)[1][2] - (*this)[0][2]*(*this)[1][1];
195  (result)[2][1] = - ((*this)[0][0]*(*this)[1][2] - (*this)[0][2]*(*this)[1][0]);
196  (result)[2][2] = (*this)[0][0]*(*this)[1][1] - (*this)[0][1]*(*this)[1][0];
197  break;
198 
199  default :
200  int i1,j1;
201  tlp::Matrix<Obj,SIZE - 1> c;
202 
203  for (unsigned int j=0; j<SIZE; ++j) {
204  for (unsigned int i=0; i<SIZE; ++i) {
205  i1 = 0;
206 
207  for (unsigned int ii=0; ii<SIZE; ++ii) {
208  if (ii == i)
209  continue;
210 
211  j1 = 0;
212 
213  for (unsigned int jj=0; jj<SIZE; jj++) {
214  if (jj == j)
215  continue;
216 
217  c[i1][j1] = (*this)[ii][jj];
218  ++j1;
219  }
220 
221  ++i1;
222  }
223 
224  if ((i+j) & 1) result[i][j]=c.determinant();
225  else result[i][j]=-c.determinant();
226  }
227  }
228 
229  break;
230  }
231 
232  return result;
233 }
234 //=====================================================================================
235 template<typename Obj,unsigned int SIZE>
236 MATRIX & MATRIX::transpose() {
237  register Obj tmp;
238 
239  for (unsigned int i=1; i<SIZE; ++i) {
240  for (unsigned int j=0; j<i; ++j) {
241  tmp = (*this)[i][j];
242  (*this)[i][j] = (*this)[j][i];
243  (*this)[j][i] = tmp;
244  }
245  }
246 
247  return (*this);
248 }
249 //=====================================================================================
250 template<typename Obj,unsigned int SIZE>
251 MATRIX & MATRIX::inverse() {
252  (*this) = (*this).cofactor().transpose() /= (*this).determinant();
253  return (*this);
254 }
255 //=====================================================================================
256 template<typename Obj,unsigned int SIZE>
257 MATRIX tlp::operator+(const MATRIX &mat1 ,const MATRIX &mat2) {
258  return MATRIX(mat1)+=mat2;
259 }
260 //=====================================================================================
261 template<typename Obj,unsigned int SIZE>
262 MATRIX tlp::operator-(const MATRIX &mat1 ,const MATRIX &mat2) {
263  return MATRIX(mat1)-=mat2;
264 }
265 //=====================================================================================
266 template<typename Obj,unsigned int SIZE>
267 MATRIX tlp::operator*(const MATRIX &mat1 ,const MATRIX &mat2) {
268  MATRIX result;
269 
270  for (unsigned int i=0; i<SIZE; ++i)
271  for (unsigned int j=0; j<SIZE; ++j) {
272  Obj tmpObj = mat1[i][0] * mat2[0][j];
273 
274  for (unsigned int k=1; k<SIZE; ++k)
275  tmpObj += mat1[i][k] * mat2[k][j];
276 
277  result[i][j] = tmpObj;
278  }
279 
280  return result;
281 }
282 //=====================================================================================
283 template<typename Obj,unsigned int SIZE>
284 MATRIX MATRIX::operator/(const MATRIX &mat2) const {
285  return MATRIX(*this)/=mat2;
286 }
287 //=====================================================================================
288 template<typename Obj,unsigned int SIZE>
289 MATRIX MATRIX::operator/(const Obj &obj) const {
290  return MATRIX(*this) /= obj;
291 }
292 //=====================================================================================
293 template<typename Obj,unsigned int SIZE>
294 MATRIX tlp::operator*(const MATRIX &mat ,const Obj &obj) {
295  return MATRIX(mat) *= obj;
296 }
297 //=====================================================================================
298 template<typename Obj, unsigned int SIZE>
299 tlp::Vector<Obj, SIZE> tlp::operator*( const MATRIX &mat , const tlp::Vector<Obj, SIZE> &vec) {
300  tlp::Vector<Obj,SIZE> result;
301 
302  for (unsigned int row=0; row<SIZE; ++row) {
303  result[row] = mat[row][0] * vec[0];
304  }
305 
306  for (unsigned int col=1; col<SIZE; ++col) {
307  for (unsigned int row=0; row<SIZE; ++row) {
308  result[row] += mat[row][col] * vec[col];
309  }
310  }
311 
312  return result;
313 }
314 //=====================================================================================
315 template<typename Obj,unsigned int SIZE>
316 tlp::Vector<Obj,SIZE> tlp::operator*( const tlp::Vector<Obj,SIZE> &vec, const MATRIX &mat) {
317  tlp::Vector<Obj,SIZE> result;
318 
319  for (unsigned int row=0; row<SIZE; ++row) {
320  result[row] = mat[0][row] * vec[0];
321  }
322 
323  for (unsigned int col=1; col<SIZE; ++col) {
324  for (unsigned int row=0; row<SIZE; ++row) {
325  result[row] += mat[col][row] * vec[col];
326  }
327  }
328 
329  return result;
330 }
331 
332 //=====================================================================================
333 template<typename Obj, unsigned int SIZE>
334 tlp::Vector<Obj, SIZE> MATRIX::powerIteration(const unsigned int nIterations) const {
335  tlp::Vector<Obj, SIZE> iteration;
336 
337  for(unsigned int i=0; i < SIZE; i++)
338  iteration[i] = 1;
339 
340  for(unsigned int i=0; i < nIterations; i++) {
341  iteration = (*this) * iteration;
342 
343  iteration /= iteration.norm();
344  }
345 
346  return iteration;
347 }
348 
349 //=====================================================================================
350 
351 template<typename Obj, unsigned int SIZE>
352 bool MATRIX::simplify(tlp::Matrix<Obj, 2> &simplifiedMatrix) const {
353  if (SIZE != 3) {
354  qWarning() << "Computation allowed only for 3x3 Matrices. Yours sizes : " << SIZE << "x" << SIZE;
355 
356  return false;
357  }
358 
359  // We start with a matrix representing an equation system under the following form :
360  //
361  // ax + by + cz = 0
362  // dx + ey + fz = 0
363  // gx + hy + iz = 0
364  //
365  // So M looks like :
366  //
367  // ( ax by cz ) *(e1)*
368  // M = ( dx ey fz ) *(e2)*
369  // ( gx hy iz ) *(e3)*
370  //
371  // What we want is something like that :
372  //
373  // jx + ky = 0
374  // lx + mz = 0
375  //
376  // So to reduce the matrix, we will use the Gaussian Elimination.
377  // For the first line we will apply a Gaussian Elimination between (e1) and (e2)
378  // For the second line we will apply a Gaussian Elimination between (e1) and (e3)
379 
380  float coeff;
381 
382  // First Gaussian Elimination :
383  // The pivot is z
384 
385  coeff = (*this)[1][2] / (*this)[0][2]; // fz / cz
386 
387  // After that:
388  // jx = dx - (coeff * ax)
389  // ky = ey - (coeff * by)
390  simplifiedMatrix[0][0] = (*this)[1][0] - (coeff * (*this)[0][0]);
391  simplifiedMatrix[0][1] = (*this)[1][1] - (coeff * (*this)[0][1]);
392 
393  // Second Gaussian Elimination :
394  // The pivot is y
395 
396  coeff = (*this)[2][1] / (*this)[0][1]; // hy / by
397 
398  // Idem :
399  // lx = gx - (coeff * ax)
400  // mz = iz - (coeff * cz)
401  simplifiedMatrix[1][0] = (*this)[2][0] - (coeff * (*this)[0][0]);
402  simplifiedMatrix[1][1] = (*this)[2][2] - (coeff * (*this)[0][2]);
403 
404  return true;
405 }
406 
407 //=====================================================================================
408 
409 template<typename Obj, unsigned int SIZE>
410 bool MATRIX::computeEigenVector(const float x, tlp::Vector<Obj, 3> &eigenVector) const {
411  if (SIZE != 2) {
412  qWarning() << "Computation allowed only for 2x2 Matrices. Yours sizes : " << SIZE << "x" << SIZE;
413 
414  return false;
415  }
416 
417  eigenVector[0] = x; // Fixed by user
418 
419  // We know that the matrix we are using is under that form :
420  //
421  // ( ax by )
422  // M = ( )
423  // ( cx dz )
424  //
425  // Since we have a fixed x, we can compute y and z :
426  //
427  // y = -a / b
428  // z = -c / d
429 
430  float a, b, c, d;
431 
432  a = (*this)[0][0];
433  b = (*this)[0][1];
434  c = (*this)[1][0];
435  d = (*this)[1][1];
436 
437  eigenVector[1] = (-a * x) / b;
438  eigenVector[2] = (-c * x) / d;
439 
440  return true;
441 }