aimsalgo  5.1.2
Neuroimaging image processing
momInvariant.h
Go to the documentation of this file.
1 /* This software and supporting documentation are distributed by
2  * Institut Federatif de Recherche 49
3  * CEA/NeuroSpin, Batiment 145,
4  * 91191 Gif-sur-Yvette cedex
5  * France
6  *
7  * This software is governed by the CeCILL-B license under
8  * French law and abiding by the rules of distribution of free software.
9  * You can use, modify and/or redistribute the software under the
10  * terms of the CeCILL-B license as circulated by CEA, CNRS
11  * and INRIA at the following URL "http://www.cecill.info".
12  *
13  * As a counterpart to the access to the source code and rights to copy,
14  * modify and redistribute granted by the license, users are provided only
15  * with a limited warranty and the software's author, the holder of the
16  * economic rights, and the successive licensors have only limited
17  * liability.
18  *
19  * In this respect, the user's attention is drawn to the risks associated
20  * with loading, using, modifying and/or developing or reproducing the
21  * software by the user in light of its specific status of free software,
22  * that may mean that it is complicated to manipulate, and that also
23  * therefore means that it is reserved for developers and experienced
24  * professionals having in-depth computer knowledge. Users are therefore
25  * encouraged to load and test the software's suitability as regards their
26  * requirements in conditions enabling the security of their systems and/or
27  * data to be ensured and, more generally, to use and operate it in the
28  * same conditions as regards security.
29  *
30  * The fact that you are presently reading this means that you have had
31  * knowledge of the CeCILL-B license and that you accept its terms.
32  */
33 
34 
35 #ifndef AIMS_MOMENT_MOMINVARIANT_H
36 #define AIMS_MOMENT_MOMINVARIANT_H
37 
38 #include <cstdlib>
39 #include <aims/moment/moment.h>
40 #include <math.h>
41 #include <complex>
42 
43 
44 template< class T > class MomentInvariant;
45 
46 
47 template< class T >
48 std::ostream& operator << ( std::ostream&, const MomentInvariant< T >& );
49 
50 
51 template< class T >
53 {
54  public:
55 
57  {
58  I00_2 = 0,
59  I22_2 = 1,
60  I222_3 = 2,
61  I11_3 = 3,
62  I33_3 = 4,
63  I3111_3 = 5,
64  I3131_3 = 6,
65  I3331_3 = 7,
66  I3333_3 = 8,
67  I112_23 = 9,
68  I123_23 = 10,
69  I233_23 = 11
70  };
71 
74  virtual ~MomentInvariant() { }
75 
76  Moment< T >& moment() const { return *_mom; }
77 
78  const double *invariant() const { return &_inv[ 0 ]; }
79 
80  void update( double, double, double, int dir=(int)Moment< T >::mAdd );
81 
83  int dir=(int)Moment< T >::mAdd );
84  void doit( AimsSurfaceTriangle& );
85  void doit( Moment< T > * );
86 
87  private:
88 
89  void init();
90  void toComplex();
91  void toInvariant();
92 
93  Moment< T > *_mom;
94 
95  double _inv[ 12 ];
96 
97  std::complex< double > _n0;
98  std::complex< double > _n1[ 3 ];
99  std::complex< double > _n2[ 5 ];
100  std::complex< double > _n3[ 7 ];
101 
102  double _c00;
103  double _c10;
104  double _c11;
105  double _c20;
106  double _c22;
107  double _c30;
108  double _c31;
109  double _c32;
110  double _c33;
111 
112  double _cc22;
113  double _cc222;
114  double _cc11;
115  double _cc33;
116  double _ce;
117  std::complex< double > _p1;
118  std::complex< double > _p2;
119  std::complex< double > _p3;
120  std::complex< double > _p4;
121  std::complex< double > _p5;
122  std::complex< double > _p6;
123  std::complex< double > _p7;
124  std::complex< double > _p8;
125  std::complex< double > _q1;
126  std::complex< double > _q2;
127  std::complex< double > _q4;
128  std::complex< double > _q5;
129  std::complex< double > _r1;
130  std::complex< double > _r2;
131 
132  std::complex< double > _de;
133 };
134 
135 
136 template< class T > inline
138 {
139  _mom = mom;
140  init();
141 }
142 
143 
144 template< class T > inline
146 {
147  _mom = &other.moment();
148 
149  init();
150 
151  for ( int i=0; i<12; i++ ) _inv[ i ] = other.invariant()[ i ];
152 }
153 
154 
155 template< class T > inline
157 {
158  // Constants for complex moment computation
159  _c00 = 2.0 * sqrt( M_PI ) / 3.0;
160 
161  _c10 = 2.0 * sqrt( 3.0 * M_PI ) / 5.0;
162  _c11 = sqrt( 6.0 * M_PI ) / 5.0;
163 
164  _c20 = 2.0 * sqrt( M_PI / 5.0 ) / 3.0;
165  _c22 = sqrt( 2.0 * M_PI / 15.0 );
166 
167  _c30 = 2.0 * sqrt( M_PI / 7.0 ) / 5.0;
168  _c31 = sqrt( 3.0 * M_PI / 7.0 ) / 5.0;
169  _c32 = sqrt( 6.0 * M_PI / 35.0 );
170  _c33 = sqrt( M_PI / 35.0 );
171 
172  // Constants for moment invariant computation
173  _cc22 = 16.0 * M_PI / 45.0 / sqrt( 5.0 );
174  _cc222 = 32.0 * M_PI * sqrt( 2.0 * M_PI / 7.0 ) / 675.0;
175  _cc11 = -4.0 * M_PI * sqrt( 3.0 ) / 25.0;
176  _cc33 = -16.0 * M_PI / 175.0 / sqrt( 7.0 );
177 
178  _ce = 1.0 / sqrt( 5.0 );
179 
180  _p1 = std::complex< double >( sqrt( 10.0 / 21.0 ), 0.0 );
181  _p2 = std::complex< double >( -2.0 * sqrt( 5.0 / 21.0 ), 0.0 );
182  _p3 = std::complex< double >( sqrt( 2.0 / 7.0 ), 0.0 );
183  _p4 = std::complex< double >( 5.0 / sqrt( 21.0 ), 0.0 );
184  _p5 = std::complex< double >( -sqrt( 5.0 / 7.0 ), 0.0 );
185  _p6 = std::complex< double >( sqrt( 2.0 / 21.0 ), 0.0 );
186  _p7 = std::complex< double >( -sqrt( 3.0 / 7.0 ), 0.0 );
187  _p8 = std::complex< double >( 2.0 / sqrt( 21.0 ), 0.0 );
188 
189  _q1 = std::complex< double >( -sqrt( 5.0 / 21.0 ), 0.0 );
190  _q2 = std::complex< double >( 1.0 / sqrt( 21.0 ), 0.0 );
191  _q4 = std::complex< double >( 1.0 / sqrt( 7.0 ), 0.0 );
192  _q5 = std::complex< double >( -2.0 * sqrt( 2.0 / 21.0 ), 0.0 );
193 
194  _r1 = std::complex< double >( sqrt( 2.0 ), 0.0 );
195  _r2 = std::complex< double >( sqrt( 2.0 / 3.0 ), 0.0 );
196 
197  _de = std::complex< double >( 2.0, 0.0 );
198 }
199 
200 
201 template< class T > inline
202 void MomentInvariant< T >::update( double x, double y, double z, int dir )
203 {
204  if ( _mom )
205  {
206  _mom->update( x, y, z, dir );
207  toComplex();
208  toInvariant();
209  }
210 }
211 
212 
213 template< class T > inline
215  T label, int dir )
216 {
217  if (_mom )
218  {
219  _mom->doit( d, label, dir );
220  toComplex();
221  toInvariant();
222  }
223 }
224 
225 
226 template< class T > inline
228 {
229  if ( _mom )
230  {
231  _mom->doit( surfTri );
232  toComplex();
233  toInvariant();
234  }
235 }
236 
237 
238 template< class T > inline
240 {
241  if ( m )
242  {
243  _mom = m;
244  toComplex();
245  toInvariant();
246  }
247 }
248 
249 
250 template< class T > inline
252 {
253  _n0 = std::complex< double >( _c00 * ( _mom->m2()[ 0 ] + _mom->m2()[ 1 ] +
254  _mom->m2()[ 2 ] ) , 0.0 );
255 
256  _n1[ 0 ] = std::complex< double >( _c11 * ( _mom->m3()[ 0 ] + _mom->m3()[ 5 ] +
257  _mom->m3()[ 7 ] ),
258  -_c11 * ( _mom->m3()[ 1 ] + _mom->m3()[ 3 ] +
259  _mom->m3()[ 8 ] ) );
260  _n1[ 1 ] = std::complex< double >( _c10 * ( _mom->m3()[ 2 ] + _mom->m3()[ 4 ] +
261  _mom->m3()[ 6 ] ), 0.0 );
262  _n1[ 2 ] = std::complex< double >( -_n1[ 0 ].real(), _n1[ 0 ].imag() );
263 
264  _n2[ 0 ] = std::complex< double >( _c22 * ( _mom->m2()[ 0 ] - _mom->m2()[ 1 ] ),
265  -2.0 * _c22 * _mom->m2()[ 3 ] );
266  _n2[ 1 ] = std::complex< double >( 2.0 * _c22 * _mom->m2()[ 4 ],
267  -2.0 * _c22 * _mom->m2()[ 5 ] );
268  _n2[ 2 ] = std::complex< double >( _c20 * ( 2.0 * _mom->m2()[ 2 ] -
269  _mom->m2()[ 0 ] - _mom->m2()[ 1 ] ),
270  0.0 );
271  _n2[ 3 ] = std::complex< double >( -_n2[ 1 ].real(), _n2[ 1 ].imag() );
272  _n2[ 4 ] = conj( _n2[ 0 ] );
273 
274  _n3[ 0 ] = std::complex< double >( _c33 * ( _mom->m3()[ 0 ] -
275  3.0 * _mom->m3()[ 5 ] ),
276  _c33 * ( _mom->m3()[ 1 ] -
277  3.0 * _mom->m3()[ 3 ] ) );
278  _n3[ 1 ] = std::complex< double >( _c32 * ( _mom->m3()[ 4 ] - _mom->m3()[ 6 ] ),
279  -2.0 * _c32 * _mom->m3()[ 9 ] );
280  _n3[ 2 ] = std::complex< double >( -_c31 * ( _mom->m3()[ 0 ] + _mom->m3()[ 5 ] -
281  4.0 * _mom->m3()[ 7 ] ),
282  _c31 * ( _mom->m3()[ 1 ] + _mom->m3()[ 3 ] -
283  4.0 * _mom->m3()[ 8 ] ) );
284  _n3[ 3 ] = std::complex< double >( _c30 * ( 2.0 * _mom->m3()[ 2 ] -
285  3.0 * ( _mom->m3()[ 4 ] +
286  _mom->m3()[ 6 ] ) ), 0.0 );
287  _n3[ 4 ] = std::complex< double >( -_n3[ 2 ].real(), _n3[ 2 ].imag() );
288  _n3[ 5 ] = conj( _n3[ 1 ] );
289  _n3[ 6 ] = std::complex< double >( -_n3[ 0 ].real(), _n3[ 0 ].imag() );
290 }
291 
292 
293 template< class T > inline
295 {
296  double mo200 = _mom->m2()[ 0 ];
297  double mo020 = _mom->m2()[ 1 ];
298  double mo002 = _mom->m2()[ 2 ];
299  double mo110 = _mom->m2()[ 3 ];
300  double mo101 = _mom->m2()[ 4 ];
301  double mo011 = _mom->m2()[ 5 ];
302 
303  double mo300 = _mom->m3()[ 0 ];
304  double mo030 = _mom->m3()[ 1 ];
305  double mo003 = _mom->m3()[ 2 ];
306  double mo210 = _mom->m3()[ 3 ];
307  double mo201 = _mom->m3()[ 4 ];
308  double mo120 = _mom->m3()[ 5 ];
309  double mo021 = _mom->m3()[ 6 ];
310  double mo102 = _mom->m3()[ 7 ];
311  double mo012 = _mom->m3()[ 8 ];
312  double mo111 = _mom->m3()[ 9 ];
313 
314  double m200_2 = mo200 * mo200;
315  double m020_2 = mo020 * mo020;
316  double m002_2 = mo002 * mo002;
317  double m110_2 = mo110 * mo110;
318  double m101_2 = mo101 * mo101;
319  double m011_2 = mo011 * mo011;
320 
321  double m200m020 = mo200 * mo020;
322  double m200m002 = mo200 * mo002;
323  double m020m002 = mo020 * mo002;
324 
325  double m200_3 = m200_2 * mo200;
326  double m020_3 = m020_2 * mo020;
327  double m002_3 = m002_2 * mo002;
328  double m200_2m020 = m200_2 * mo020;
329  double m200_2m002 = m200_2 * mo002;
330  double m020_2m200 = m020_2 * mo200;
331  double m020_2m002 = m020_2 * mo002;
332  double m002_2m200 = m002_2 * mo200;
333  double m002_2m020 = m002_2 * mo020;
334  double m200m110_2 = mo200 * m110_2;
335  double m200m101_2 = mo200 * m101_2;
336  double m200m011_2 = mo200 * m011_2;
337  double m020m110_2 = mo020 * m110_2;
338  double m020m101_2 = mo020 * m101_2;
339  double m020m011_2 = mo020 * m011_2;
340  double m002m110_2 = mo002 * m110_2;
341  double m002m101_2 = mo002 * m101_2;
342  double m002m011_2 = mo002 * m011_2;
343  double m200m020m002 = m200m020 * mo002;
344  double m110m101m011 = mo110 * mo101 * mo011;
345 
346  double m300_2 = mo300 * mo300;
347  double m030_2 = mo030 * mo030;
348  double m003_2 = mo003 * mo003;
349  double m210_2 = mo210 * mo210;
350  double m201_2 = mo201 * mo201;
351  double m120_2 = mo120 * mo120;
352  double m021_2 = mo021 * mo021;
353  double m102_2 = mo102 * mo102;
354  double m012_2 = mo012 * mo012;
355  double m111_2 = mo111 * mo111;
356 
357  double m300m120 = mo300 * mo120;
358  double m300m102 = mo300 * mo102;
359  double m030m210 = mo030 * mo210;
360  double m030m012 = mo030 * mo012;
361  double m003m201 = mo003 * mo201;
362  double m003m021 = mo003 * mo021;
363  double m120m102 = mo120 * mo102;
364  double m210m012 = mo210 * mo012;
365  double m201m021 = mo201 * mo021;
366 
367  double m0 = _mom->m0();
368  double d1 = pow( m0, 5.0 / 3.0 );
369  double d2 = pow( m0, 10.0 / 3.0 );
370  double d3 = pow( m0, 5.0 );
371  double d4 = pow( m0, 4.0 );
372  double d5 = pow( m0, 17.0 / 3.0 );
373  double d6 = pow( m0, 8.0 );
374 
375  // Moment invariant of order 2 derived from algebraic calculus
376  double I0 = _c00 * ( mo200 + mo020 + mo002 ) / d1;
377  double aI0 = fabs( I0 );
378  _inv[ 0 ] = I0 * sqrt( aI0 ) / aI0;
379 
380  double I1 = _cc22 * ( m200_2 + m020_2 + m002_2 +
381  3.0 * ( m110_2 + m101_2 + m011_2 ) -
382  m200m020 - m200m002 - m020m002 ) / d2;
383  double aI1 = fabs( I1 );
384  _inv[ 1 ] = I1 * pow( aI1, 0.25 ) / aI1;
385 
386  double I2 = _cc222 * ( -2.0 * ( m200_3 + m020_3 + m002_3 ) +
387  3 * ( m200_2m002 + m200_2m020 + m020_2m200 +
388  m020_2m002 + m002_2m200 + m002_2m020 ) -
389  9.0 * ( m200m110_2 + m200m101_2 + m020m110_2 +
390  m020m011_2 + m002m101_2 + m002m011_2 ) +
391  18.0 * ( m200m011_2 + m020m101_2 + m002m110_2 ) -
392  12.0 * m200m020m002 - 54.0 * m110m101m011 ) / d3;
393  double aI2 = fabs( I2 );
394  _inv[ 2 ] = I2 * pow( aI2, 1.0 / 6.0 ) / aI2;
395 
396  // Moment invariant of order 3 derived from quantum mechanic calculus
397  double I3 = _cc11 * ( m300_2 + m030_2 + m003_2 + m210_2 + m201_2 + m120_2 +
398  m021_2 + m102_2 + m012_2 +
399  2.0 * ( m300m120 + m300m102 + m030m210 + m030m012 +
400  m003m201 + m003m021 + m210m012 + m201m021 +
401  m120m102 ) ) / d4;
402  double aI3 = fabs( I3 );
403  _inv[ 3 ] = I3 * pow( aI3, 1.0 / 6.0 ) / aI3;
404 
405  double I4 = _cc33 * ( m300_2 + m030_2 + m003_2 +
406  6.0 * ( m210_2 + m201_2 + m120_2 + m021_2 + m102_2 +
407  m012_2 ) + 15.0 * m111_2 -
408  3.0 * ( m300m120 + m300m102 + m030m210 + m030m012 +
409  m003m201 + m003m021 + m210m012 + m201m021 +
410  m120m102 ) ) / d4;
411  double aI4 = fabs( I4 );
412  _inv[ 4 ] = I4 * pow( aI4, 1.0 / 6.0 ) / aI4;
413 
414  std::complex< double > zeta[ 5 ];
415  zeta[ 0 ] = _p1 * _n3[ 0 ] * _n3[ 4 ] + _p2 * _n3[ 1 ] * _n3[ 3 ] +
416  _p3 * _n3[ 2 ] * _n3[ 2 ];
417  zeta[ 1 ] = _p4 * _n3[ 0 ] * _n3[ 5 ] + _p5 * _n3[ 1 ] * _n3[ 4 ] +
418  _p6 * _n3[ 2 ] * _n3[ 3 ];
419  zeta[ 2 ] = _p4 * _n3[ 0 ] * _n3[ 6 ] + _p7 * _n3[ 2 ] * _n3[ 4 ] +
420  _p8 * _n3[ 3 ] * _n3[ 3 ];
421  zeta[ 3 ] = _p4 * _n3[ 1 ] * _n3[ 6 ] + _p5 * _n3[ 2 ] * _n3[ 5 ] +
422  _p6 * _n3[ 3 ] * _n3[ 4 ];
423  zeta[ 4 ] = _p1 * _n3[ 2 ] * _n3[ 6 ] + _p2 * _n3[ 3 ] * _n3[ 5 ] +
424  _p3 * _n3[ 4 ] * _n3[ 4 ];
425 
426  std::complex< double > xi[ 5 ];
427  xi[ 0 ] = _q2 * _n3[ 2 ] * _n1[ 0 ] + _q1 * _n3[ 1 ] * _n1[ 1 ] -
428  _p5 * _n3[ 0 ] * _n1[ 2 ];
429  xi[ 1 ] = _q4 * _n3[ 3 ] * _n1[ 0 ] + _q5 * _n3[ 2 ] * _n1[ 1 ] +
430  _p1 * _n3[ 1 ] * _n1[ 2 ];
431  xi[ 2 ] = _p3 * _n3[ 4 ] * _n1[ 0 ] + _p7 * _n3[ 3 ] * _n1[ 1 ] +
432  _p3 * _n3[ 2 ] * _n1[ 2 ];
433  xi[ 3 ] = _p1 * _n3[ 5 ] * _n1[ 0 ] + _q5 * _n3[ 4 ] * _n1[ 1 ] +
434  _q4 * _n3[ 3 ] * _n1[ 2 ];
435  xi[ 4 ] = -_p5 * _n3[ 6 ] * _n1[ 0 ] + _q1 * _n3[ 5 ] * _n1[ 1 ] +
436  _q2 * _n3[ 4 ] * _n1[ 2 ];
437 
438  std::complex< double > nu[ 5 ];
439  nu[ 0 ] = _n1[ 0 ] * _n1[ 0 ];
440  nu[ 1 ] = _r1 * _n1[ 0 ] * _n1[ 1 ];
441  nu[ 2 ] = _r2 * ( _n1[ 0 ] * _n1[ 2 ] + _n1[ 1 ] * _n1[ 1 ] );
442  nu[ 3 ] = _r1 * _n1[ 1 ] * _n1[ 2 ];
443  nu[ 4 ] = _n1[ 2 ] * _n1[ 2 ];
444 
445  std::complex< double > dI5 = xi[ 4 ] * nu[ 0 ] - xi[ 3 ] * nu[ 1 ] +
446  xi[ 2 ] * nu[ 2 ] - xi[ 1 ] * nu[ 3 ] +
447  xi[ 0 ] * nu[ 4 ];
448  double I5 = _ce * dI5.real() / d6;
449  double aI5 = fabs( I5 );
450  _inv[ 5 ] = I5 * pow( aI5, 1.0 / 12.0 ) / aI5;
451 
452  std::complex< double > dI6 = _de * xi[ 0 ] * xi[ 4 ] - _de * xi[ 1 ] * xi[ 3 ] +
453  xi[ 2 ] * xi[ 2 ];
454  double I6 = _ce * dI6.real() / d6;
455  double aI6 = fabs( I6 );
456  _inv[ 6 ] = I6 * pow( aI6, 1.0 / 12.0 ) / aI6;
457 
458  std::complex< double > dI7 = zeta[ 4 ] * xi[ 0 ] - zeta[ 3 ] * xi[ 1 ] +
459  zeta[ 2 ] * xi[ 2 ] - zeta[ 1 ] * xi[ 3 ] +
460  zeta[ 0 ] * xi[ 4 ];
461  double I7 = _ce * dI7.real() / d6;
462  double aI7 = fabs( I7 );
463  _inv[ 7 ] = I7 * pow( aI7, 1.0 / 12.0 ) / aI7;
464 
465  std::complex< double > dI8 = _de * zeta[ 0 ] * zeta[ 4 ] -
466  _de * zeta[ 1 ] * zeta[ 3 ] + zeta[ 2 ] * zeta[ 2 ];
467  double I8 = _ce * dI8.real() / d6;
468  double aI8 = fabs( I8 );
469  _inv[ 8 ] = I8 * pow( aI8, 1.0 / 12.0 ) / aI8;
470 
471  std::complex< double > dI9 = nu[ 4 ] * _n2[ 0 ] - nu[ 3 ] * _n2[ 1 ] +
472  nu[ 2 ] * _n2[ 2 ] - nu[ 1 ] * _n2[ 3 ] +
473  nu[ 0 ] * _n2[ 4 ];
474  double I9 = _ce * dI9.real() / d5;
475  double aI9 = fabs( I9 );
476  _inv[ 9 ] = I9 * pow( aI9, 0.125 ) / aI9;
477 
478  std::complex< double > dI10 = xi[ 4 ] * _n2[ 0 ] - xi[ 3 ] * _n2[ 1 ] +
479  xi[ 2 ] * _n2[ 2 ] - xi[ 1 ] * _n2[ 3 ] +
480  xi[ 0 ] * _n2[ 4 ];
481  double I10 = _ce * dI10.real() / d5;
482  double aI10 = fabs( I10 );
483  _inv[ 10 ] = I10 * pow( aI10, 0.125 ) / aI10;
484 
485  std::complex< double > dI11 = zeta[ 4 ] * _n2[ 0 ] - zeta[ 3 ] * _n2[ 1 ] +
486  zeta[ 2 ] * _n2[ 2 ] - zeta[ 1 ] * _n2[ 3 ] +
487  zeta[ 0 ] * _n2[ 4 ];
488  double I11 = _ce * dI11.real() / d5;
489  double aI11 = fabs( I11 );
490  _inv[ 11 ] = I11 * pow( aI11, 0.125 ) / aI11;
491 }
492 
493 
494 template< class T > inline
495 std::ostream& operator << ( std::ostream& out,
496  const MomentInvariant< T >& mom )
497 {
498  for ( int i=0; i<12; i++ ) out << mom.invariant()[ i ] << " ";
499 
500  return out;
501 }
502 
503 #endif
void update(double, double, double, int dir=(int) Moment< T >::mAdd)
Definition: momInvariant.h:202
void doit(carto::rc_ptr< carto::Volume< T > > &, T, int dir=(int) Moment< T >::mAdd)
Definition: momInvariant.h:214
MomentInvariant(Moment< T > *m=0)
Definition: momInvariant.h:137
virtual ~MomentInvariant()
Definition: momInvariant.h:74
const double * invariant() const
Definition: momInvariant.h:78
Moment< T > & moment() const
Definition: momInvariant.h:76
Definition: moment.h:56
std::ostream & operator<<(std::ostream &, const MomentInvariant< T > &)
Definition: momInvariant.h:495
AIMSDATA_API AimsTimeSurface< 3, Void > AimsSurfaceTriangle