A.I.M.S algorithms


moment.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_MOMENT_H
36 #define AIMS_MOMENT_MOMENT_H
37 
38 #include <fstream>
39 #include <iostream>
40 
41 #include <aims/data/data.h>
42 #include <aims/mesh/surface.h>
43 #include <aims/bucket/bucket.h>
44 #include <aims/math/eigen.h>
45 
46 
47 template< class T > class Moment;
48 
49 
50 template< class T >
51 std::ostream& operator << ( std::ostream&, const Moment< T >& );
52 
53 
54 template< class T >
55 class Moment
56 {
57  public:
58 
59  enum Operation
60  {
61  mSub = -1,
62  mAdd = 1
63  };
64 
65  enum MomentId
66  {
67  m000 = 0,
68 
69  m100 = 0,
70  m010 = 1,
71  m001 = 2,
72 
73  m200 = 0,
74  m020 = 1,
75  m002 = 2,
76  m110 = 3,
77  m101 = 4,
78  m011 = 5,
79 
80  m300 = 0,
81  m030 = 1,
82  m003 = 2,
83  m210 = 3,
84  m201 = 4,
85  m120 = 5,
86  m021 = 6,
87  m102 = 7,
88  m012 = 8,
89  m111 = 9
90  };
91 
93  Moment();
94  Moment( const Moment< T >& );
95  virtual ~Moment() { }
96 
97  virtual void clear();
98  double cx() const { return _cx; }
99  double cy() const { return _cy; }
100  double cz() const { return _cz; }
101  double ct() const { return _ct; }
102  void setcx( double x ) { _cx = x; }
103  void setcy( double x ) { _cy = x; }
104  void setcz( double x ) { _cz = x; }
105  void setct( double x ) { _ct = x; }
106 
107  double *gravity() { return &_g[ 0 ]; }
108  const double *gravity() const { return &_g[ 0 ]; }
109 
110  double& sum() { return _sum; }
111  const double& sum() const { return _sum; }
112  double& m0() { return _m0; }
113  const double& m0() const { return _m0; }
114 
115  double *m1() { return &_m1[ 0 ]; }
116  const double *m1() const { return &_m1[ 0 ]; }
117  double *m2() { return &_m2[ 0 ]; }
118  const double *m2() const { return &_m2[ 0 ]; }
119  double *m3() { return &_m3[ 0 ]; }
120  const double *m3() const { return &_m3[ 0 ]; }
121 
122  const AimsData< double >& eigenValue() const { return _eigenValue; }
123  const AimsData< double >& eigenVector() const { return _eigenVector; }
124 
125  virtual void update( Point3df&, int );
126  virtual void update( AimsVector<double,3>&, int );
127  virtual void update( double, double, double, int );
128 
129  // TODO: fix this API, having these methods left unimplemented is asking
130  // for trouble...
131  virtual void doit( AimsData< T >&, T, int );
132  virtual void doit( AimsSurfaceTriangle& );
133  virtual void doit( const aims::BucketMap<Void> & );
134 
135  virtual void orientation();
136 
137  Moment<T> & operator += ( const Moment<T> & );
138  Moment<T> operator + ( const Moment<T> & );
139 
140  protected:
141 
142  double _cx;
143  double _cy;
144  double _cz;
145  double _ct;
146 
147  double _g[ 3 ];
148 
149  double _sum;
150  double _m0;
151 
152  double _m1[ 3 ];
153  double _m2[ 6 ];
154  double _m3[ 10 ];
155 
158 };
159 
160 
161 template< class T > inline
162 Moment< T >::Moment() : _cx( 1.0 ), _cy( 1.0 ), _cz( 1.0 ), _ct( 1.0 ),
163  _sum( 0.0 ), _m0( 0.0 )
164 {
165 }
166 
167 
168 template< class T > inline
170  : _cx( other.cx() ), _cy( other.cy() ), _cz( other.cz() ),
171  _ct( other.ct() ), _sum( other.sum() ), _m0( other.m0() ),
172  _eigenValue( other.eigenValue() ),
173  _eigenVector( other.eigenVector() )
174 {
175  int i;
176 
177  for ( i=0; i<3; i++ )
178  {
179  _g[ i ] = other.gravity()[ i ];
180  _m1[ i ] = other.m1()[ i ];
181  }
182  for ( i=0; i<6; i++ ) _m2[ i ] = other.m2()[ i ];
183  for ( i=0; i<9; i++ ) _m3[ i ] = other.m3()[ i ];
184 }
185 
186 
187 template< class T > inline
189 {
190  _m0 = 0;
191  _m1[0] = 0;
192  _m1[1] = 0;
193  _m1[2] = 0;
194  _m2[0] = 0;
195  _m2[1] = 0;
196  _m2[2] = 0;
197  _m2[3] = 0;
198  _m2[4] = 0;
199  _m2[5] = 0;
200  _m3[0] = 0;
201  _m3[1] = 0;
202  _m3[2] = 0;
203  _m3[3] = 0;
204  _m3[4] = 0;
205  _m3[5] = 0;
206  _m3[6] = 0;
207  _m3[7] = 0;
208  _m3[8] = 0;
209  _m3[9] = 0;
210  _sum = 0;
211  _g[0] = 0;
212  _g[1] = 0;
213  _g[2] = 0;
214  // cx..ct ?
215 }
216 
217 
218 // make GCC warning non-fatal
219 #pragma GCC diagnostic push
220 #pragma GCC diagnostic ignored "-Woverloaded-virtual"
221 template< class T > inline
222 void Moment< T >::update( Point3df& pt, int dir )
223 {
224  update( (double)pt[ 0 ], (double)pt[ 1 ], (double)pt[ 2 ], dir );
225 }
226 
227 
228 template< class T > inline
230 {
231  update( pt[ 0 ], pt[ 1 ], pt[ 2 ], dir );
232 }
233 
234 
235 template< class T > inline
236 void Moment< T >::update( double, double, double, int )
237 {
238  orientation();
239 }
240 #pragma GCC diagnostic pop
241 
242 
243 template< class T > inline
245 {
246  _eigenVector = AimsData< double >( 3, 3 );
247 
248  _eigenVector( 0, 0 ) = _m2[ 0 ];
249  _eigenVector( 1, 1 ) = _m2[ 1 ];
250  _eigenVector( 2, 2 ) = _m2[ 2 ];
251  _eigenVector( 0, 1 ) = _eigenVector( 1, 0 ) = _m2[ 3 ];
252  _eigenVector( 0, 2 ) = _eigenVector( 2, 0 ) = _m2[ 4 ];
253  _eigenVector( 1, 2 ) = _eigenVector( 2, 1 ) = _m2[ 5 ];
254 
255  AimsEigen< double > eigen;
256  _eigenValue = eigen.doit( _eigenVector );
257  eigen.sort( _eigenVector, _eigenValue );
258 }
259 
260 
261 template< class T > inline
262 std::ostream& operator << ( std::ostream& out , const Moment< T >& mom )
263 {
264  int i;
265 
266  out << "Image geometry: (" << mom.cx() << "x" << mom.cy() << "x" << mom.cz();
267  out << ")" << std::endl << "Volume : " << mom.m0() << "mm3" << std::endl;
268  out << "Number of points : " << mom.sum() << std::endl;
269  out << "Volume (mm3) : " << mom.m0() << std::endl;
270  out << "Order 1 moments : ";
271  for ( i=0; i<3; i++ ) out << mom.m1()[ i ] << " ";
272  out << std::endl << "Order 2 moments : ";
273  for ( i=0; i<6; i++ ) out << mom.m2()[ i ] << " ";
274  out << std::endl << "Order 3 moments : ";
275  for ( i=0; i<10; i++ ) out << mom.m3()[ i ] << " ";
276  out << std::endl << "Center of Mass : (" << mom.gravity()[ 0 ] << ",";
277  out << mom.gravity()[ 1 ] << "," << mom.gravity()[ 2 ] << ")" << std::endl;
278  out << "Orientation :" << std::endl << " Inerty : ";
279  out << mom.eigenValue()( 0, 0 ) << " ; " << mom.eigenValue()( 1, 1 );
280  out << " ; " << mom.eigenValue()( 2, 2 );
281  out << std::endl << " V1 : (" << mom.eigenVector()( 0, 0 ) << ",";
282  out << mom.eigenVector()( 1, 0 ) << "," << mom.eigenVector()( 2, 0 ) << ")";
283  out << std::endl << " V2 : (" << mom.eigenVector()( 0, 1 ) << ",";
284  out << mom.eigenVector()( 1, 1 ) << "," << mom.eigenVector()( 2, 1 ) << ")";
285  out << std::endl << " V3 : (" << mom.eigenVector()( 0, 2 ) << ",";
286  out << mom.eigenVector()( 1, 2 ) << "," << mom.eigenVector()( 2, 2 ) << ")";
287 
288  return out;
289 }
290 
291 
292 template< class T > inline
294 {
295  // un-center m2 and m3 which are centered, contrarily to m1
296  // don't ask me why (denis 22/07/2006)
297  if( _sum == 0 )
298  {
299  *this = m;
300  return *this;
301  }
302  else if( m._sum == 0 )
303  return *this;
304 
305  double g0 = _m1[0] / _m0;
306  double g1 = _m1[1] / _m0;
307  double g2 = _m1[2] / _m0;
308 
309  _m3[0] += _m0*g0*g0*g0 + g0*_m2[0]*3;
310  _m3[1] += _m0*g1*g1*g1 + g1*_m2[1]*3;
311  _m3[2] += _m0*g2*g2*g2 + g2*_m2[2]*3;
312  _m3[3] += _m0*g0*g0*g1 + g0*_m2[3]*2 + g1*_m2[0];
313  _m3[4] += _m0*g0*g0*g2 + g0*_m2[4]*2 + g2*_m2[0];
314  _m3[5] += _m0*g0*g1*g1 + g0*_m2[1] + g1*_m2[3]*2;
315  _m3[6] += _m0*g1*g1*g2 + g1*_m2[5]*2 + g2*_m2[1];
316  _m3[7] += _m0*g0*g2*g2 + g0*_m2[2] + g2*_m2[4]*2;
317  _m3[8] += _m0*g1*g2*g2 + g1*_m2[2] + g2*_m2[5]*2;
318  _m3[9] += _m0*g0*g1*g2 + g0*_m2[5] + g1*_m2[4] +g2*_m2[3];
319 
320  _m2[0] += _m0 * g0*g0;
321  _m2[1] += _m0 * g1*g1;
322  _m2[2] += _m0 * g2*g2;
323  _m2[3] += _m0 * g0*g1;
324  _m2[4] += _m0 * g0*g2;
325  _m2[5] += _m0 * g1*g2;
326 
327  double m20 = m._m2[0];
328  double m21 = m._m2[1];
329  double m22 = m._m2[2];
330  double m23 = m._m2[3];
331  double m24 = m._m2[4];
332  double m25 = m._m2[5];
333 
334  double h0 = m._m1[0] / m._m0;
335  double h1 = m._m1[1] / m._m0;
336  double h2 = m._m1[2] / m._m0;
337 
338  _sum += m._sum;
339  _m0 += m._m0;
340  _m1[0] += m._m1[0];
341  _m1[1] += m._m1[1];
342  _m1[2] += m._m1[2];
343 
344  _m3[0] += m._m3[0] + m._m0*h0*h0*h0 + h0*m20*3;
345  _m3[1] += m._m3[1] + m._m0*h1*h1*h1 + h1*m21*3;
346  _m3[2] += m._m3[2] + m._m0*h2*h2*h2 + h2*m22*3;
347  _m3[3] += m._m3[3] + m._m0*h0*h0*h1 + h0*m23*2 + h1*m20;
348  _m3[4] += m._m3[4] + m._m0*h0*h0*h2 + h0*m24*2 + h2*m20;
349  _m3[5] += m._m3[5] + m._m0*h0*h1*h1 + h0*m21 + h1*m23*2;
350  _m3[6] += m._m3[6] + m._m0*h1*h1*h2 + h1*m25*2 + h2*m21;
351  _m3[7] += m._m3[7] + m._m0*h0*h2*h2 + h0*m22 + h2*m24*2;
352  _m3[8] += m._m3[8] + m._m0*h1*h2*h2 + h1*m22 + h2*m25*2;
353  _m3[9] += m._m3[9] + m._m0*h0*h1*h2 + h0*m25 + h1*m24 +h2*m23;
354 
355  _m2[0] += m._m2[0] + m._m0 * h0*h0;
356  _m2[1] += m._m2[1] + m._m0 * h1*h1;
357  _m2[2] += m._m2[2] + m._m0 * h2*h2;
358  _m2[3] += m._m2[3] + m._m0 * h0*h1;
359  _m2[4] += m._m2[4] + m._m0 * h0*h2;
360  _m2[5] += m._m2[5] + m._m0 * h1*h2;
361 
362  // TODO: check that _cx..._ct match
363 
364  _g[0] = _m1[0] / _m0;
365  _g[1] = _m1[1] / _m0;
366  _g[2] = _m1[2] / _m0;
367 
368  g0 = _g[0];
369  g1 = _g[1];
370  g2 = _g[2];
371 
372  // center back m2 and m3
373  _m2[0] -= _m0 * g0*g0;
374  _m2[1] -= _m0 * g1*g1;
375  _m2[2] -= _m0 * g2*g2;
376  _m2[3] -= _m0 * g0*g1;
377  _m2[4] -= _m0 * g0*g2;
378  _m2[5] -= _m0 * g1*g2;
379 
380  _m3[0] -= _m0*g0*g0*g0 + g0*_m2[0]*3;
381  _m3[1] -= _m0*g1*g1*g1 + g1*_m2[1]*3;
382  _m3[2] -= _m0*g2*g2*g2 + g2*_m2[2]*3;
383  _m3[3] -= _m0*g0*g0*g1 + g0*_m2[3]*2 + g1*_m2[0];
384  _m3[4] -= _m0*g0*g0*g2 + g0*_m2[4]*2 + g2*_m2[0];
385  _m3[5] -= _m0*g0*g1*g1 + g0*_m2[1] + g1*_m2[3]*2;
386  _m3[6] -= _m0*g1*g1*g2 + g1*_m2[5]*2 + g2*_m2[1];
387  _m3[7] -= _m0*g0*g2*g2 + g0*_m2[2] + g2*_m2[4]*2;
388  _m3[8] -= _m0*g1*g2*g2 + g1*_m2[2] + g2*_m2[5]*2;
389  _m3[9] -= _m0*g0*g1*g2 + g0*_m2[5] + g1*_m2[4] +g2*_m2[3];
390 
391  orientation();
392  return *this;
393 }
394 
395 
396 template< class T > inline
398 {
399  Moment<T> mom( *this );
400  mom += m;
401  return mom;
402 }
403 
404 // make warning non-fatal (the API is broken, see comment in Moment<T>)
405 #pragma GCC diagnostic push
406 #pragma GCC diagnostic ignored "-Woverloaded-virtual"
407 template< class T >
409 {
410  std::clog << "PROGRAMMING ERROR: called unimplemented base method"
411  " Moment<T>::doit( AimsData< T >&, T, int )"
412  << std::endl;
413 #ifndef NDEBUG
414  throw std::logic_error("called unimplemented base method"
415  " Moment<T>::doit( AimsData< T >&, T, int )");
416 #endif
417 }
418 
419 template< class T >
421 {
422  std::clog << "PROGRAMMING ERROR: called unimplemented base method"
423  " Moment<T>::doit( AimsSurfaceTriangle& )"
424  << std::endl;
425 #ifndef NDEBUG
426  throw std::logic_error("called unimplemented base method"
427  " Moment<T>::doit( AimsSurfaceTriangle& )");
428 #endif
429 }
430 
431 template< class T >
433 {
434  std::clog << "PROGRAMMING ERROR: called unimplemented base method"
435  " Moment<T>::doit( const aims::BucketMap<Void> & )"
436  << std::endl;
437 #ifndef NDEBUG
438  throw std::logic_error("called unimplemented base method"
439  " Moment<T>::doit( const aims::BucketMap<Void> & )");
440 #endif
441 }
442 #pragma GCC diagnostic pop
443 
444 #endif
AimsData< double > _eigenVector
Definition: moment.h:157
double * m3()
Definition: moment.h:119
virtual ~Moment()
Definition: moment.h:95
double cz() const
Definition: moment.h:100
AIMSDATA_API AimsTimeSurface< 3, Void > AimsSurfaceTriangle
double _ct
Definition: moment.h:145
const double & m0() const
Definition: moment.h:113
double _m3[10]
Definition: moment.h:154
double _m2[6]
Definition: moment.h:153
const double & sum() const
Definition: moment.h:111
double cx() const
Definition: moment.h:98
double _g[3]
Definition: moment.h:147
double _cx
Definition: moment.h:142
double & sum()
Definition: moment.h:110
const AimsData< double > & eigenVector() const
Definition: moment.h:123
double * m1()
Definition: moment.h:115
virtual void clear()
Definition: moment.h:188
const AimsData< double > & eigenValue() const
Definition: moment.h:122
double & m0()
Definition: moment.h:112
void setcx(double x)
Definition: moment.h:102
const double * gravity() const
Definition: moment.h:108
Moment< T > & operator+=(const Moment< T > &)
Definition: moment.h:293
Definition: moment.h:47
AimsData< double > _eigenValue
Definition: moment.h:156
double _m1[3]
Definition: moment.h:152
double cy() const
Definition: moment.h:99
virtual void update(Point3df &, int)
Definition: moment.h:222
double _sum
Definition: moment.h:149
Moment< T > operator+(const Moment< T > &)
Definition: moment.h:397
void setct(double x)
Definition: moment.h:105
double _cy
Definition: moment.h:143
Moment()
doesn't necessarily clear all moments to 0. Call clear() to be sure
Definition: moment.h:162
AimsData< T > doit(AimsData< T > &, AimsData< T > *wi=NULL)
Eigen system resolution function.
double ct() const
Definition: moment.h:101
double _m0
Definition: moment.h:150
virtual void orientation()
Definition: moment.h:244
const double * m3() const
Definition: moment.h:120
void setcy(double x)
Definition: moment.h:103
virtual void doit(AimsData< T > &, T, int)
Definition: moment.h:408
double _cz
Definition: moment.h:144
double * m2()
Definition: moment.h:117
double * gravity()
Definition: moment.h:107
void setcz(double x)
Definition: moment.h:104
void sort(AimsData< T > &eigenvectors, AimsData< T > &eigenvalues, AimsData< T > *wi=NULL)
Sort the eigenvectors and eigenvalues in decreasing order.
const double * m2() const
Definition: moment.h:118
const double * m1() const
Definition: moment.h:116