aimsalgo  5.0.5
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 >
52 class MomentInvariant
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 
82  void doit( AimsData< T >&, T, int dir=(int)Moment< T >::mAdd );
83  void doit( AimsSurfaceTriangle& );
84  void doit( Moment< T > * );
85 
86  private:
87 
88  void init();
89  void toComplex();
90  void toInvariant();
91 
92  Moment< T > *_mom;
93 
94  double _inv[ 12 ];
95 
96  std::complex< double > _n0;
97  std::complex< double > _n1[ 3 ];
98  std::complex< double > _n2[ 5 ];
99  std::complex< double > _n3[ 7 ];
100 
101  double _c00;
102  double _c10;
103  double _c11;
104  double _c20;
105  double _c22;
106  double _c30;
107  double _c31;
108  double _c32;
109  double _c33;
110 
111  double _cc22;
112  double _cc222;
113  double _cc11;
114  double _cc33;
115  double _ce;
116  std::complex< double > _p1;
117  std::complex< double > _p2;
118  std::complex< double > _p3;
119  std::complex< double > _p4;
120  std::complex< double > _p5;
121  std::complex< double > _p6;
122  std::complex< double > _p7;
123  std::complex< double > _p8;
124  std::complex< double > _q1;
125  std::complex< double > _q2;
126  std::complex< double > _q4;
127  std::complex< double > _q5;
128  std::complex< double > _r1;
129  std::complex< double > _r2;
130 
131  std::complex< double > _de;
132 };
133 
134 
135 template< class T > inline
137 {
138  _mom = mom;
139  init();
140 }
141 
142 
143 template< class T > inline
145 {
146  _mom = &other.moment();
147 
148  init();
149 
150  for ( int i=0; i<12; i++ ) _inv[ i ] = other.invariant()[ i ];
151 }
152 
153 
154 template< class T > inline
156 {
157  // Constants for complex moment computation
158  _c00 = 2.0 * sqrt( M_PI ) / 3.0;
159 
160  _c10 = 2.0 * sqrt( 3.0 * M_PI ) / 5.0;
161  _c11 = sqrt( 6.0 * M_PI ) / 5.0;
162 
163  _c20 = 2.0 * sqrt( M_PI / 5.0 ) / 3.0;
164  _c22 = sqrt( 2.0 * M_PI / 15.0 );
165 
166  _c30 = 2.0 * sqrt( M_PI / 7.0 ) / 5.0;
167  _c31 = sqrt( 3.0 * M_PI / 7.0 ) / 5.0;
168  _c32 = sqrt( 6.0 * M_PI / 35.0 );
169  _c33 = sqrt( M_PI / 35.0 );
170 
171  // Constants for moment invariant computation
172  _cc22 = 16.0 * M_PI / 45.0 / sqrt( 5.0 );
173  _cc222 = 32.0 * M_PI * sqrt( 2.0 * M_PI / 7.0 ) / 675.0;
174  _cc11 = -4.0 * M_PI * sqrt( 3.0 ) / 25.0;
175  _cc33 = -16.0 * M_PI / 175.0 / sqrt( 7.0 );
176 
177  _ce = 1.0 / sqrt( 5.0 );
178 
179  _p1 = std::complex< double >( sqrt( 10.0 / 21.0 ), 0.0 );
180  _p2 = std::complex< double >( -2.0 * sqrt( 5.0 / 21.0 ), 0.0 );
181  _p3 = std::complex< double >( sqrt( 2.0 / 7.0 ), 0.0 );
182  _p4 = std::complex< double >( 5.0 / sqrt( 21.0 ), 0.0 );
183  _p5 = std::complex< double >( -sqrt( 5.0 / 7.0 ), 0.0 );
184  _p6 = std::complex< double >( sqrt( 2.0 / 21.0 ), 0.0 );
185  _p7 = std::complex< double >( -sqrt( 3.0 / 7.0 ), 0.0 );
186  _p8 = std::complex< double >( 2.0 / sqrt( 21.0 ), 0.0 );
187 
188  _q1 = std::complex< double >( -sqrt( 5.0 / 21.0 ), 0.0 );
189  _q2 = std::complex< double >( 1.0 / sqrt( 21.0 ), 0.0 );
190  _q4 = std::complex< double >( 1.0 / sqrt( 7.0 ), 0.0 );
191  _q5 = std::complex< double >( -2.0 * sqrt( 2.0 / 21.0 ), 0.0 );
192 
193  _r1 = std::complex< double >( sqrt( 2.0 ), 0.0 );
194  _r2 = std::complex< double >( sqrt( 2.0 / 3.0 ), 0.0 );
195 
196  _de = std::complex< double >( 2.0, 0.0 );
197 }
198 
199 
200 template< class T > inline
201 void MomentInvariant< T >::update( double x, double y, double z, int dir )
202 {
203  if ( _mom )
204  {
205  _mom->update( x, y, z, dir );
206  toComplex();
207  toInvariant();
208  }
209 }
210 
211 
212 template< class T > inline
213 void MomentInvariant< T >::doit( AimsData< T >& d, T label, int dir )
214 {
215  if (_mom )
216  {
217  _mom->doit( d, label, dir );
218  toComplex();
219  toInvariant();
220  }
221 }
222 
223 
224 template< class T > inline
226 {
227  if ( _mom )
228  {
229  _mom->doit( surfTri );
230  toComplex();
231  toInvariant();
232  }
233 }
234 
235 
236 template< class T > inline
238 {
239  if ( m )
240  {
241  _mom = m;
242  toComplex();
243  toInvariant();
244  }
245 }
246 
247 
248 template< class T > inline
250 {
251  _n0 = std::complex< double >( _c00 * ( _mom->m2()[ 0 ] + _mom->m2()[ 1 ] +
252  _mom->m2()[ 2 ] ) , 0.0 );
253 
254  _n1[ 0 ] = std::complex< double >( _c11 * ( _mom->m3()[ 0 ] + _mom->m3()[ 5 ] +
255  _mom->m3()[ 7 ] ),
256  -_c11 * ( _mom->m3()[ 1 ] + _mom->m3()[ 3 ] +
257  _mom->m3()[ 8 ] ) );
258  _n1[ 1 ] = std::complex< double >( _c10 * ( _mom->m3()[ 2 ] + _mom->m3()[ 4 ] +
259  _mom->m3()[ 6 ] ), 0.0 );
260  _n1[ 2 ] = std::complex< double >( -_n1[ 0 ].real(), _n1[ 0 ].imag() );
261 
262  _n2[ 0 ] = std::complex< double >( _c22 * ( _mom->m2()[ 0 ] - _mom->m2()[ 1 ] ),
263  -2.0 * _c22 * _mom->m2()[ 3 ] );
264  _n2[ 1 ] = std::complex< double >( 2.0 * _c22 * _mom->m2()[ 4 ],
265  -2.0 * _c22 * _mom->m2()[ 5 ] );
266  _n2[ 2 ] = std::complex< double >( _c20 * ( 2.0 * _mom->m2()[ 2 ] -
267  _mom->m2()[ 0 ] - _mom->m2()[ 1 ] ),
268  0.0 );
269  _n2[ 3 ] = std::complex< double >( -_n2[ 1 ].real(), _n2[ 1 ].imag() );
270  _n2[ 4 ] = conj( _n2[ 0 ] );
271 
272  _n3[ 0 ] = std::complex< double >( _c33 * ( _mom->m3()[ 0 ] -
273  3.0 * _mom->m3()[ 5 ] ),
274  _c33 * ( _mom->m3()[ 1 ] -
275  3.0 * _mom->m3()[ 3 ] ) );
276  _n3[ 1 ] = std::complex< double >( _c32 * ( _mom->m3()[ 4 ] - _mom->m3()[ 6 ] ),
277  -2.0 * _c32 * _mom->m3()[ 9 ] );
278  _n3[ 2 ] = std::complex< double >( -_c31 * ( _mom->m3()[ 0 ] + _mom->m3()[ 5 ] -
279  4.0 * _mom->m3()[ 7 ] ),
280  _c31 * ( _mom->m3()[ 1 ] + _mom->m3()[ 3 ] -
281  4.0 * _mom->m3()[ 8 ] ) );
282  _n3[ 3 ] = std::complex< double >( _c30 * ( 2.0 * _mom->m3()[ 2 ] -
283  3.0 * ( _mom->m3()[ 4 ] +
284  _mom->m3()[ 6 ] ) ), 0.0 );
285  _n3[ 4 ] = std::complex< double >( -_n3[ 2 ].real(), _n3[ 2 ].imag() );
286  _n3[ 5 ] = conj( _n3[ 1 ] );
287  _n3[ 6 ] = std::complex< double >( -_n3[ 0 ].real(), _n3[ 0 ].imag() );
288 }
289 
290 
291 template< class T > inline
293 {
294  double mo200 = _mom->m2()[ 0 ];
295  double mo020 = _mom->m2()[ 1 ];
296  double mo002 = _mom->m2()[ 2 ];
297  double mo110 = _mom->m2()[ 3 ];
298  double mo101 = _mom->m2()[ 4 ];
299  double mo011 = _mom->m2()[ 5 ];
300 
301  double mo300 = _mom->m3()[ 0 ];
302  double mo030 = _mom->m3()[ 1 ];
303  double mo003 = _mom->m3()[ 2 ];
304  double mo210 = _mom->m3()[ 3 ];
305  double mo201 = _mom->m3()[ 4 ];
306  double mo120 = _mom->m3()[ 5 ];
307  double mo021 = _mom->m3()[ 6 ];
308  double mo102 = _mom->m3()[ 7 ];
309  double mo012 = _mom->m3()[ 8 ];
310  double mo111 = _mom->m3()[ 9 ];
311 
312  double m200_2 = mo200 * mo200;
313  double m020_2 = mo020 * mo020;
314  double m002_2 = mo002 * mo002;
315  double m110_2 = mo110 * mo110;
316  double m101_2 = mo101 * mo101;
317  double m011_2 = mo011 * mo011;
318 
319  double m200m020 = mo200 * mo020;
320  double m200m002 = mo200 * mo002;
321  double m020m002 = mo020 * mo002;
322 
323  double m200_3 = m200_2 * mo200;
324  double m020_3 = m020_2 * mo020;
325  double m002_3 = m002_2 * mo002;
326  double m200_2m020 = m200_2 * mo020;
327  double m200_2m002 = m200_2 * mo002;
328  double m020_2m200 = m020_2 * mo200;
329  double m020_2m002 = m020_2 * mo002;
330  double m002_2m200 = m002_2 * mo200;
331  double m002_2m020 = m002_2 * mo020;
332  double m200m110_2 = mo200 * m110_2;
333  double m200m101_2 = mo200 * m101_2;
334  double m200m011_2 = mo200 * m011_2;
335  double m020m110_2 = mo020 * m110_2;
336  double m020m101_2 = mo020 * m101_2;
337  double m020m011_2 = mo020 * m011_2;
338  double m002m110_2 = mo002 * m110_2;
339  double m002m101_2 = mo002 * m101_2;
340  double m002m011_2 = mo002 * m011_2;
341  double m200m020m002 = m200m020 * mo002;
342  double m110m101m011 = mo110 * mo101 * mo011;
343 
344  double m300_2 = mo300 * mo300;
345  double m030_2 = mo030 * mo030;
346  double m003_2 = mo003 * mo003;
347  double m210_2 = mo210 * mo210;
348  double m201_2 = mo201 * mo201;
349  double m120_2 = mo120 * mo120;
350  double m021_2 = mo021 * mo021;
351  double m102_2 = mo102 * mo102;
352  double m012_2 = mo012 * mo012;
353  double m111_2 = mo111 * mo111;
354 
355  double m300m120 = mo300 * mo120;
356  double m300m102 = mo300 * mo102;
357  double m030m210 = mo030 * mo210;
358  double m030m012 = mo030 * mo012;
359  double m003m201 = mo003 * mo201;
360  double m003m021 = mo003 * mo021;
361  double m120m102 = mo120 * mo102;
362  double m210m012 = mo210 * mo012;
363  double m201m021 = mo201 * mo021;
364 
365  double m0 = _mom->m0();
366  double d1 = pow( m0, 5.0 / 3.0 );
367  double d2 = pow( m0, 10.0 / 3.0 );
368  double d3 = pow( m0, 5.0 );
369  double d4 = pow( m0, 4.0 );
370  double d5 = pow( m0, 17.0 / 3.0 );
371  double d6 = pow( m0, 8.0 );
372 
373  // Moment invariant of order 2 derived from algebraic calculus
374  double I0 = _c00 * ( mo200 + mo020 + mo002 ) / d1;
375  double aI0 = fabs( I0 );
376  _inv[ 0 ] = I0 * sqrt( aI0 ) / aI0;
377 
378  double I1 = _cc22 * ( m200_2 + m020_2 + m002_2 +
379  3.0 * ( m110_2 + m101_2 + m011_2 ) -
380  m200m020 - m200m002 - m020m002 ) / d2;
381  double aI1 = fabs( I1 );
382  _inv[ 1 ] = I1 * pow( aI1, 0.25 ) / aI1;
383 
384  double I2 = _cc222 * ( -2.0 * ( m200_3 + m020_3 + m002_3 ) +
385  3 * ( m200_2m002 + m200_2m020 + m020_2m200 +
386  m020_2m002 + m002_2m200 + m002_2m020 ) -
387  9.0 * ( m200m110_2 + m200m101_2 + m020m110_2 +
388  m020m011_2 + m002m101_2 + m002m011_2 ) +
389  18.0 * ( m200m011_2 + m020m101_2 + m002m110_2 ) -
390  12.0 * m200m020m002 - 54.0 * m110m101m011 ) / d3;
391  double aI2 = fabs( I2 );
392  _inv[ 2 ] = I2 * pow( aI2, 1.0 / 6.0 ) / aI2;
393 
394  // Moment invariant of order 3 derived from quantum mechanic calculus
395  double I3 = _cc11 * ( m300_2 + m030_2 + m003_2 + m210_2 + m201_2 + m120_2 +
396  m021_2 + m102_2 + m012_2 +
397  2.0 * ( m300m120 + m300m102 + m030m210 + m030m012 +
398  m003m201 + m003m021 + m210m012 + m201m021 +
399  m120m102 ) ) / d4;
400  double aI3 = fabs( I3 );
401  _inv[ 3 ] = I3 * pow( aI3, 1.0 / 6.0 ) / aI3;
402 
403  double I4 = _cc33 * ( m300_2 + m030_2 + m003_2 +
404  6.0 * ( m210_2 + m201_2 + m120_2 + m021_2 + m102_2 +
405  m012_2 ) + 15.0 * m111_2 -
406  3.0 * ( m300m120 + m300m102 + m030m210 + m030m012 +
407  m003m201 + m003m021 + m210m012 + m201m021 +
408  m120m102 ) ) / d4;
409  double aI4 = fabs( I4 );
410  _inv[ 4 ] = I4 * pow( aI4, 1.0 / 6.0 ) / aI4;
411 
412  std::complex< double > zeta[ 5 ];
413  zeta[ 0 ] = _p1 * _n3[ 0 ] * _n3[ 4 ] + _p2 * _n3[ 1 ] * _n3[ 3 ] +
414  _p3 * _n3[ 2 ] * _n3[ 2 ];
415  zeta[ 1 ] = _p4 * _n3[ 0 ] * _n3[ 5 ] + _p5 * _n3[ 1 ] * _n3[ 4 ] +
416  _p6 * _n3[ 2 ] * _n3[ 3 ];
417  zeta[ 2 ] = _p4 * _n3[ 0 ] * _n3[ 6 ] + _p7 * _n3[ 2 ] * _n3[ 4 ] +
418  _p8 * _n3[ 3 ] * _n3[ 3 ];
419  zeta[ 3 ] = _p4 * _n3[ 1 ] * _n3[ 6 ] + _p5 * _n3[ 2 ] * _n3[ 5 ] +
420  _p6 * _n3[ 3 ] * _n3[ 4 ];
421  zeta[ 4 ] = _p1 * _n3[ 2 ] * _n3[ 6 ] + _p2 * _n3[ 3 ] * _n3[ 5 ] +
422  _p3 * _n3[ 4 ] * _n3[ 4 ];
423 
424  std::complex< double > xi[ 5 ];
425  xi[ 0 ] = _q2 * _n3[ 2 ] * _n1[ 0 ] + _q1 * _n3[ 1 ] * _n1[ 1 ] -
426  _p5 * _n3[ 0 ] * _n1[ 2 ];
427  xi[ 1 ] = _q4 * _n3[ 3 ] * _n1[ 0 ] + _q5 * _n3[ 2 ] * _n1[ 1 ] +
428  _p1 * _n3[ 1 ] * _n1[ 2 ];
429  xi[ 2 ] = _p3 * _n3[ 4 ] * _n1[ 0 ] + _p7 * _n3[ 3 ] * _n1[ 1 ] +
430  _p3 * _n3[ 2 ] * _n1[ 2 ];
431  xi[ 3 ] = _p1 * _n3[ 5 ] * _n1[ 0 ] + _q5 * _n3[ 4 ] * _n1[ 1 ] +
432  _q4 * _n3[ 3 ] * _n1[ 2 ];
433  xi[ 4 ] = -_p5 * _n3[ 6 ] * _n1[ 0 ] + _q1 * _n3[ 5 ] * _n1[ 1 ] +
434  _q2 * _n3[ 4 ] * _n1[ 2 ];
435 
436  std::complex< double > nu[ 5 ];
437  nu[ 0 ] = _n1[ 0 ] * _n1[ 0 ];
438  nu[ 1 ] = _r1 * _n1[ 0 ] * _n1[ 1 ];
439  nu[ 2 ] = _r2 * ( _n1[ 0 ] * _n1[ 2 ] + _n1[ 1 ] * _n1[ 1 ] );
440  nu[ 3 ] = _r1 * _n1[ 1 ] * _n1[ 2 ];
441  nu[ 4 ] = _n1[ 2 ] * _n1[ 2 ];
442 
443  std::complex< double > dI5 = xi[ 4 ] * nu[ 0 ] - xi[ 3 ] * nu[ 1 ] +
444  xi[ 2 ] * nu[ 2 ] - xi[ 1 ] * nu[ 3 ] +
445  xi[ 0 ] * nu[ 4 ];
446  double I5 = _ce * dI5.real() / d6;
447  double aI5 = fabs( I5 );
448  _inv[ 5 ] = I5 * pow( aI5, 1.0 / 12.0 ) / aI5;
449 
450  std::complex< double > dI6 = _de * xi[ 0 ] * xi[ 4 ] - _de * xi[ 1 ] * xi[ 3 ] +
451  xi[ 2 ] * xi[ 2 ];
452  double I6 = _ce * dI6.real() / d6;
453  double aI6 = fabs( I6 );
454  _inv[ 6 ] = I6 * pow( aI6, 1.0 / 12.0 ) / aI6;
455 
456  std::complex< double > dI7 = zeta[ 4 ] * xi[ 0 ] - zeta[ 3 ] * xi[ 1 ] +
457  zeta[ 2 ] * xi[ 2 ] - zeta[ 1 ] * xi[ 3 ] +
458  zeta[ 0 ] * xi[ 4 ];
459  double I7 = _ce * dI7.real() / d6;
460  double aI7 = fabs( I7 );
461  _inv[ 7 ] = I7 * pow( aI7, 1.0 / 12.0 ) / aI7;
462 
463  std::complex< double > dI8 = _de * zeta[ 0 ] * zeta[ 4 ] -
464  _de * zeta[ 1 ] * zeta[ 3 ] + zeta[ 2 ] * zeta[ 2 ];
465  double I8 = _ce * dI8.real() / d6;
466  double aI8 = fabs( I8 );
467  _inv[ 8 ] = I8 * pow( aI8, 1.0 / 12.0 ) / aI8;
468 
469  std::complex< double > dI9 = nu[ 4 ] * _n2[ 0 ] - nu[ 3 ] * _n2[ 1 ] +
470  nu[ 2 ] * _n2[ 2 ] - nu[ 1 ] * _n2[ 3 ] +
471  nu[ 0 ] * _n2[ 4 ];
472  double I9 = _ce * dI9.real() / d5;
473  double aI9 = fabs( I9 );
474  _inv[ 9 ] = I9 * pow( aI9, 0.125 ) / aI9;
475 
476  std::complex< double > dI10 = xi[ 4 ] * _n2[ 0 ] - xi[ 3 ] * _n2[ 1 ] +
477  xi[ 2 ] * _n2[ 2 ] - xi[ 1 ] * _n2[ 3 ] +
478  xi[ 0 ] * _n2[ 4 ];
479  double I10 = _ce * dI10.real() / d5;
480  double aI10 = fabs( I10 );
481  _inv[ 10 ] = I10 * pow( aI10, 0.125 ) / aI10;
482 
483  std::complex< double > dI11 = zeta[ 4 ] * _n2[ 0 ] - zeta[ 3 ] * _n2[ 1 ] +
484  zeta[ 2 ] * _n2[ 2 ] - zeta[ 1 ] * _n2[ 3 ] +
485  zeta[ 0 ] * _n2[ 4 ];
486  double I11 = _ce * dI11.real() / d5;
487  double aI11 = fabs( I11 );
488  _inv[ 11 ] = I11 * pow( aI11, 0.125 ) / aI11;
489 }
490 
491 
492 template< class T > inline
493 std::ostream& operator << ( std::ostream& out,
494  const MomentInvariant< T >& mom )
495 {
496  for ( int i=0; i<12; i++ ) out << mom.invariant()[ i ] << " ";
497 
498  return out;
499 }
500 
501 #endif
virtual ~MomentInvariant()
Definition: momInvariant.h:74
AIMSDATA_API AimsTimeSurface< 3, Void > AimsSurfaceTriangle
void doit(AimsData< T > &, T, int dir=(int) Moment< T >::mAdd)
Definition: momInvariant.h:213
Definition: moment.h:47
const double * invariant() const
Definition: momInvariant.h:78
Moment< T > & moment() const
Definition: momInvariant.h:76
MomentInvariant(Moment< T > *m=0)
Definition: momInvariant.h:136
std::ostream & operator<<(std::ostream &, const MomentInvariant< T > &)
Definition: momInvariant.h:493
void update(double, double, double, int dir=(int) Moment< T >::mAdd)
Definition: momInvariant.h:201