aimsalgo 6.0.0
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
44template< class T > class MomentInvariant;
45
46
47template< class T >
48std::ostream& operator << ( std::ostream&, const MomentInvariant< T >& );
49
50
51template< 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,
68 I123_23 = 10,
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 );
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
136template< class T > inline
138{
139 _mom = mom;
140 init();
141}
142
143
144template< 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
155template< class T > inline
156void MomentInvariant< T >::init()
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
201template< class T > inline
202void 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
213template< 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
226template< class T > inline
228{
229 if ( _mom )
230 {
231 _mom->doit( surfTri );
232 toComplex();
233 toInvariant();
234 }
235}
236
237
238template< class T > inline
240{
241 if ( m )
242 {
243 _mom = m;
244 toComplex();
245 toInvariant();
246 }
247}
248
249
250template< class T > inline
251void MomentInvariant< T >::toComplex()
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
293template< class T > inline
294void MomentInvariant< T >::toInvariant()
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
494template< class T > inline
495std::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)
void doit(carto::rc_ptr< carto::Volume< T > > &, T, int dir=(int) Moment< T >::mAdd)
MomentInvariant(Moment< T > *m=0)
virtual ~MomentInvariant()
Moment< T > & moment() const
const double * invariant() const
@ mAdd
Definition moment.h:62
std::ostream & operator<<(std::ostream &, const MomentInvariant< T > &)
AIMSDATA_API AimsTimeSurface< 3, Void > AimsSurfaceTriangle