35 #ifndef AIMS_MOMENT_MOMINVARIANT_H
36 #define AIMS_MOMENT_MOMINVARIANT_H
78 const double *
invariant()
const {
return &_inv[ 0 ]; }
97 std::complex< double > _n0;
98 std::complex< double > _n1[ 3 ];
99 std::complex< double > _n2[ 5 ];
100 std::complex< double > _n3[ 7 ];
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;
132 std::complex< double > _de;
136 template<
class T >
inline
144 template<
class T >
inline
151 for (
int i=0; i<12; i++ ) _inv[ i ] = other.
invariant()[ i ];
155 template<
class T >
inline
159 _c00 = 2.0 * sqrt(
M_PI ) / 3.0;
161 _c10 = 2.0 * sqrt( 3.0 *
M_PI ) / 5.0;
162 _c11 = sqrt( 6.0 *
M_PI ) / 5.0;
164 _c20 = 2.0 * sqrt(
M_PI / 5.0 ) / 3.0;
165 _c22 = sqrt( 2.0 *
M_PI / 15.0 );
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 );
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 );
178 _ce = 1.0 / sqrt( 5.0 );
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 );
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 );
194 _r1 = std::complex< double >( sqrt( 2.0 ), 0.0 );
195 _r2 = std::complex< double >( sqrt( 2.0 / 3.0 ), 0.0 );
197 _de = std::complex< double >( 2.0, 0.0 );
201 template<
class T >
inline
206 _mom->update( x, y, z, dir );
213 template<
class T >
inline
219 _mom->doit( d, label, dir );
226 template<
class T >
inline
231 _mom->doit( surfTri );
238 template<
class T >
inline
250 template<
class T >
inline
253 _n0 = std::complex< double >( _c00 * ( _mom->m2()[ 0 ] + _mom->m2()[ 1 ] +
254 _mom->m2()[ 2 ] ) , 0.0 );
256 _n1[ 0 ] = std::complex< double >( _c11 * ( _mom->m3()[ 0 ] + _mom->m3()[ 5 ] +
258 -_c11 * ( _mom->m3()[ 1 ] + _mom->m3()[ 3 ] +
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() );
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 ] ),
271 _n2[ 3 ] = std::complex< double >( -_n2[ 1 ].real(), _n2[ 1 ].imag() );
272 _n2[ 4 ] = conj( _n2[ 0 ] );
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() );
293 template<
class T >
inline
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 ];
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 ];
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;
321 double m200m020 = mo200 * mo020;
322 double m200m002 = mo200 * mo002;
323 double m020m002 = mo020 * mo002;
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;
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;
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;
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 );
376 double I0 = _c00 * ( mo200 + mo020 + mo002 ) / d1;
377 double aI0 = fabs( I0 );
378 _inv[ 0 ] = I0 * sqrt( aI0 ) / aI0;
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;
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;
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 +
402 double aI3 = fabs( I3 );
403 _inv[ 3 ] = I3 * pow( aI3, 1.0 / 6.0 ) / aI3;
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 +
411 double aI4 = fabs( I4 );
412 _inv[ 4 ] = I4 * pow( aI4, 1.0 / 6.0 ) / aI4;
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 ];
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 ];
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 ];
445 std::complex< double > dI5 = xi[ 4 ] * nu[ 0 ] - xi[ 3 ] * nu[ 1 ] +
446 xi[ 2 ] * nu[ 2 ] - xi[ 1 ] * nu[ 3 ] +
448 double I5 = _ce * dI5.real() / d6;
449 double aI5 = fabs( I5 );
450 _inv[ 5 ] = I5 * pow( aI5, 1.0 / 12.0 ) / aI5;
452 std::complex< double > dI6 = _de * xi[ 0 ] * xi[ 4 ] - _de * xi[ 1 ] * xi[ 3 ] +
454 double I6 = _ce * dI6.real() / d6;
455 double aI6 = fabs( I6 );
456 _inv[ 6 ] = I6 * pow( aI6, 1.0 / 12.0 ) / aI6;
458 std::complex< double > dI7 = zeta[ 4 ] * xi[ 0 ] - zeta[ 3 ] * xi[ 1 ] +
459 zeta[ 2 ] * xi[ 2 ] - zeta[ 1 ] * xi[ 3 ] +
461 double I7 = _ce * dI7.real() / d6;
462 double aI7 = fabs( I7 );
463 _inv[ 7 ] = I7 * pow( aI7, 1.0 / 12.0 ) / aI7;
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;
471 std::complex< double > dI9 = nu[ 4 ] * _n2[ 0 ] - nu[ 3 ] * _n2[ 1 ] +
472 nu[ 2 ] * _n2[ 2 ] - nu[ 1 ] * _n2[ 3 ] +
474 double I9 = _ce * dI9.real() / d5;
475 double aI9 = fabs( I9 );
476 _inv[ 9 ] = I9 * pow( aI9, 0.125 ) / aI9;
478 std::complex< double > dI10 = xi[ 4 ] * _n2[ 0 ] - xi[ 3 ] * _n2[ 1 ] +
479 xi[ 2 ] * _n2[ 2 ] - xi[ 1 ] * _n2[ 3 ] +
481 double I10 = _ce * dI10.real() / d5;
482 double aI10 = fabs( I10 );
483 _inv[ 10 ] = I10 * pow( aI10, 0.125 ) / aI10;
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;
494 template<
class T >
inline
498 for (
int i=0; i<12; i++ ) out << mom.
invariant()[ i ] <<
" ";
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()
const double * invariant() const
Moment< T > & moment() const
std::ostream & operator<<(std::ostream &, const MomentInvariant< T > &)
AIMSDATA_API AimsTimeSurface< 3, Void > AimsSurfaceTriangle