3d/matrix.hh
Go to the documentation of this file.
1 /* -*- mia-c++ -*-
2  *
3  * This file is part of MIA - a toolbox for medical image analysis
4  * Copyright (c) Leipzig, Madrid 1999-2013 Gert Wollny
5  *
6  * MIA is free software; you can redistribute it and/or modify
7  * it under the terms of the GNU General Public License as published by
8  * the Free Software Foundation; either version 3 of the License, or
9  * (at your option) any later version.
10  *
11  * This program is distributed in the hope that it will be useful,
12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  * GNU General Public License for more details.
15  *
16  * You should have received a copy of the GNU General Public License
17  * along with MIA; if not, see <http://www.gnu.org/licenses/>.
18  *
19  */
20 #ifndef __mia_3d_matrix_hh
21 #define __mia_3d_matrix_hh
22 
23 #include <mia/3d/vector.hh>
24 #include <mia/core/msgstream.hh>
25 
27 
38 template <typename T>
39 struct T3DMatrix: public T3DVector< T3DVector<T> > {
40 
41  T3DMatrix() = default;
42  T3DMatrix(const T3DMatrix<T>& o) = default;
43 
44 
50  static T3DMatrix<T> diagonal(T value);
51 
57  static T3DMatrix<T> diagonal(const T3DVector<T>& values);
58 
59 
65  template <typename I>
66  T3DMatrix(const T3DMatrix<I>& o);
67 
68 
74  T3DMatrix(const T3DVector< T3DVector<T> >& other);
75 
76 
83  T3DMatrix(const T3DVector< T >& x, const T3DVector< T >& y, const T3DVector< T >& z );
84 
90  T3DMatrix<T>& operator -= (const T3DMatrix<T>& other);
91 
92 
97  void print( std::ostream& os) const;
98 
99 
103  T3DMatrix<T> transposed()const;
104 
105 
109  T get_det() const;
110 
111 
115  int get_rank()const;
116 
126  int get_eigenvalues(C3DFVector& v)const;
127 
135  int get_eigenvector(float ev, C3DFVector& v)const;
136 
137 
139  static const T3DMatrix _1;
140 
142  static const T3DMatrix _0;
143 
144 };
145 
146 
149 
150 
151 template <typename T>
153  T3DVector< T >(0,1,0),
154  T3DVector< T >(0,0,1));
155 
156 template <typename T>
158 
159 
160 template <typename T>
162 {
163  return T3DMatrix<T>(T3DVector< T >(v,0,0),
164  T3DVector< T >(0,v,0),
165  T3DVector< T >(0,0,v));
166 }
167 
168 template <typename T>
170 {
171  return T3DMatrix<T>(T3DVector< T >(v.x,0,0),
172  T3DVector< T >(0,v.y,0),
173  T3DVector< T >(0,0,v.z));
174 }
175 
176 template <typename T>
177 template <typename I>
179  T3DVector<T3DVector<T> >(o.x, o.y, o.z)
180 {
181 }
182 
183 template <typename T>
185  T3DVector<T3DVector<T> >(other.x, other.y, other.z)
186 {
187 }
188 
189 template <typename T>
191  T3DVector<T3DVector<T> >(x, y, z)
192 {
193 }
194 template <typename T>
195 void T3DMatrix<T>::print( std::ostream& os) const
196 {
197  os << "<" << this->x << ", " << this->y << ", " << this->z << " >";
198 }
199 
200 template <typename T>
201 std::ostream& operator << (std::ostream& os, const T3DMatrix<T>& m)
202 {
203  m.print(os);
204  return os;
205 }
206 
207 template <typename T>
209 {
210  this->x -= o.x;
211  this->y -= o.y;
212  this->z -= o.z;
213  return *this;
214 }
215 
216 template <typename T>
218 {
219  return T3DMatrix<T>(T3DVector<T>(this->x.x, this->y.x, this->z.x),
220  T3DVector<T>(this->x.y, this->y.y, this->z.y),
221  T3DVector<T>(this->x.z, this->y.z, this->z.z));
222 }
223 
224 template <typename T>
226 {
227  return T3DVector<T>(dot(m.x, x), dot(m.y, x), dot(m.z, x));
228 }
229 
230 template <typename T>
232 {
233  return T3DVector<T>(m.x.x * x.x + m.y.x * x.y + m.z.x * x.z,
234  m.x.y * x.x + m.y.y * x.y + m.z.y * x.z,
235  m.x.z * x.x + m.y.z * x.y + m.z.z * x.z);
236 }
237 
238 template <typename T>
240 {
241  return T3DMatrix<T>(T3DVector<T>(m.x.x * x.x.x + m.x.y * x.y.x + m.x.z * x.z.x,
242  m.x.x * x.x.y + m.x.y * x.y.y + m.x.z * x.z.y,
243  m.x.x * x.x.z + m.x.y * x.y.z + m.x.z * x.z.z),
244  T3DVector<T>(m.y.x * x.x.x + m.y.y * x.y.x + m.y.z * x.z.x,
245  m.y.x * x.x.y + m.y.y * x.y.y + m.y.z * x.z.y,
246  m.y.x * x.x.z + m.y.y * x.y.z + m.y.z * x.z.z),
247  T3DVector<T>(m.z.x * x.x.x + m.z.y * x.y.x + m.z.z * x.z.x,
248  m.z.x * x.x.y + m.z.y * x.y.y + m.z.z * x.z.y,
249  m.z.x * x.x.z + m.z.y * x.y.z + m.z.z * x.z.z));
250 }
251 
252 
253 template <typename T>
255 {
256  C3DFVector ev;
257  this->get_eigenvalues(ev);
258  cvdebug()<< "Matrix = "<< *this <<", Rank: eigenvalues: " << ev << "\n";
259  int rank = 0;
260  if (ev.x != 0.0)
261  rank++;
262  if (ev.y != 0.0)
263  rank++;
264  if (ev.z != 0.0)
265  rank++;
266  return rank;
267 }
268 
269 inline double cubrt(double a)
270 {
271  if ( a == 0.0 )
272  return 0.0;
273  return a > 0.0 ? pow(a,1.0/3.0) : - pow(-a, 1.0/ 3.0);
274 }
275 
276 
277 template <class T>
279 {
280  return dot(this->x,T3DVector<T>(this->y.y * this->z.z - this->z.y * this->y.z,
281  this->y.z * this->z.x - this->y.x * this->z.z,
282  this->y.x * this->z.y - this->y.y * this->z.x));
283 }
284 
285 
286 template <typename T>
288 {
289  int retval = 0;
290 
291  double t = - get_det();
292  double s = this->x.x * this->y.y + this->z.z * this->y.y + this->x.x * this->z.z -
293  this->x.y * this->y.x - this->x.z * this->z.x - this->y.z * this->z.y;
294 
295  double r = - this->x.x - this->y.y - this->z.z;
296 
297  // a = 1;
298 
299  double p = s - r * r / 3.0;
300  double q = ( 27.0 * t - 9.0 * r * s + 2.0 * r * r * r ) / 27.0;
301 
302  double diskr = q *q / 4.0 + p * p * p / 27.0;
303 
304  cvdebug() << "discr =" << diskr << "\n";
305  if ( diskr > 1e-6 ) {
306  // complex solution
307  double sqrt_discr = sqrt(diskr);
308  double u = cubrt( - q/2.0 + sqrt_discr );
309  double v = cubrt( - q/2.0 - sqrt_discr );
310  result.x = u + v - r / 3.0;
311  result.y = -(u + v) / 2.0 - r / 3.0; // real part
312  result.z = (u - v)/2.0 * sqrt(3.0f); // imag part
313  return 1;
314 
315  }
316 
317  std::vector<double> res(3);
318  if ( diskr < -1e-6) {
319  double rho = sqrt(-p*p*p/ 27.0);
320  double cphi = - q / ( 2.0 * rho);
321  double phi = acos(cphi)/3.0;
322  double sqrt_p = 2 * cubrt(rho);
323  res[0] = sqrt_p * cos(phi)- r/3.;
324  res[1] = sqrt_p * cos(phi + M_PI * 2.0 / 3.0) - r/3.;
325  res[2] = sqrt_p * cos(phi + M_PI * 4.0 / 3.0) - r/3.;
326  retval = 3;
327  } else { // at least two values are equal, all real
328  double u = cubrt( - q/2.0 );
329  res[0] = 2.0 * u - r / 3.0;
330  res[1] = - u - r / 3.0;
331  res[2] = res[1];
332  retval = 2;
333  }
334 
335  std::sort(res.begin(), res.end(), [](double x, double y){ return std::fabs(x) > std::fabs(y);});
336 
337  for_each(res.begin(), res.end(), [](double& x){ if (std::fabs(x) < 1e-12) x = 0.0;});
338 
339  result.x = res[0];
340  result.y = res[1];
341  result.z = res[2];
342  return retval;
343 
344 }
345 
359 template<class T>
360 bool solve_2x2(T a11, T a12,T b1,T a21, T a22,T b2,T *x1,T *x2)
361 {
362  T h1 = a11 * a22 - a12 * a21;
363  if (h1 == T())
364  return false;
365 
366  *x1 = (b1 * a22 - b2 * a12) / h1;
367  *x2 = (b2 * a11 - b1 * a21) / h1;
368  return true;
369 }
370 
373  int a,b;
374 };
375 
376 
377 template <typename T>
379 {
380  const solve_lines_t l[3] = { {0,1}, {1,2}, {2,0}};
381  // T b1,b2,a11,a12,a21,a22;
382  if (ev == 0.0) {
383  return 1;
384  }
385 
386  T3DMatrix<T> M = *this - T3DMatrix<T>::diagonal(ev);
387 
388  float x = std::abs(M.x.x)+std::abs(M.y.x)+std::abs(M.z.x);
389  float y = std::abs(M.x.y)+std::abs(M.y.y)+std::abs(M.z.y);
390  float z = std::abs(M.x.z)+std::abs(M.y.z)+std::abs(M.z.z);
391 
392  if (x+y+z == 0.0) {
393  v = T3DVector<T>(1,0,0);
394  return 0;
395  }
396 
397  T3DVector<int> col;
398  T *rx;
399  T *ry;
400 
401  // thats tricky:
402  // 1st col index is 1st column in solver
403  // 2nd col index is 2nd column in solver
404  // 3th col index is right side
405  // the respective result value is 0=x,1=y,2=z
406  // the right side presenting value of result is preset with 1.0
407  // the others get a pointer
408 
409  if (x < y) {
410  if (x < z){
411  col = T3DVector<int>(1,2,0);
412  rx = &v.y;
413  ry = &v.z;
414  v.x = 1.0;
415  }else{
416  col = T3DVector<int>(0,1,2);
417  rx = &v.x;
418  ry = &v.y;
419  v.z = 1.0;
420  }
421  }else{
422  if (y < z){
423  col = T3DVector<int>(0,2,1);
424  rx = &v.x;
425  ry = &v.z;
426  v.y = 1.0;
427 
428  }else{
429  col = T3DVector<int>(0,1,2);
430  rx = &v.x;
431  ry = &v.y;
432  v.z= 1.0;
433  }
434  }
435 
436  bool good = false;
437  for (int i = 0; i < 3 && !good; i++) {
438  good = solve_2x2(M[l[i].a][col.x],M[l[i].a][col.y],-M[l[i].a][col.z],
439  M[l[i].b][col.x],M[l[i].b][col.y],-M[l[i].b][col.z],
440  rx,ry);
441  }
442  // seems there is no solution
443  if (!good)
444  return 2;
445 
446 
447  if ((M * v).norm2() > 1e-5) {
448 
449  // a solution for only two rows is not a solution
450  // but there is no better
451  fprintf(stderr,"WARNING: rank of A-ev*I\n numerical > 2");
452  return 0;
453  }
454  v /= v.norm();
455  return 0;
456 }
457 
458 
459 
460 
462 
463 #endif