A.I.M.S algorithms


momNormStgy.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_MOMNORMSTGY_H
36 #define AIMS_MOMENT_MOMNORMSTGY_H
37 
38 #include <aims/moment/momStgy.h>
39 
40 
41 template< class T >
43 {
44  public:
45 
47 
48  void doit( Moment< T > *, AimsData< T >&, T, int );
49  void doit( Moment< T > *, const aims::BucketMap<Void> &, int );
50 
51  private:
52  void addVoxel2( Moment<T>*, const Point3d & p );
53 };
54 
55 
56 template< class T > inline
58 {
59  double cx = m->cx();
60  double cy = m->cy();
61  double cz = m->cz();
62  double mx = m->m1()[ 0 ] / m->m0();
63  double my = m->m1()[ 1 ] / m->m0();
64  double mz = m->m1()[ 2 ] / m->m0();
65 
66  double tx = (double)p[0] * cx - mx;
67  double ty = (double)p[1] * cy - my;
68  double tz = (double)p[2] * cz - mz;
69 
70  m->m2()[ 0 ] += tx * tx;
71  m->m2()[ 1 ] += ty * ty;
72  m->m2()[ 2 ] += tz * tz;
73  m->m2()[ 3 ] += tx * ty;
74  m->m2()[ 4 ] += tx * tz;
75  m->m2()[ 5 ] += ty * tz;
76 
77  m->m3()[ 0 ] += tx * tx * tx;
78  m->m3()[ 1 ] += ty * ty * ty;
79  m->m3()[ 2 ] += tz * tz * tz;
80  m->m3()[ 3 ] += tx * tx * ty;
81  m->m3()[ 4 ] += tx * tx * tz;
82  m->m3()[ 5 ] += tx * ty * ty;
83  m->m3()[ 6 ] += ty * ty * tz;
84  m->m3()[ 7 ] += tx * tz * tz;
85  m->m3()[ 8 ] += ty * tz * tz;
86  m->m3()[ 9 ] += tx * ty * tz;
87 }
88 
89 
90 template< class T > inline
92  T label, int )
93 {
94  int x, y, z;
95  int dx = d.dimX();
96  int dy = d.dimY();
97  int dz = d.dimZ();
98  int olbs = d.oLineBetweenSlice();
99  int opbl = d.oPointBetweenLine();
100  typename AimsData< T >::const_iterator it = d.begin() + d.oFirstPoint();
101 
102  double cx = m->cx();
103  double cy = m->cy();
104  double cz = m->cz();
105  double ct = m->ct();
106 
107  for ( z=0; z<dz; z++, it+=olbs )
108  for ( y=0; y<dy; y++, it+=opbl )
109  for ( x=0; x<dx; x++, it++ )
110  if ( *it == label )
111  {
112  m->sum()++;
113 
114  m->m0()++;
115 
116  m->m1()[ 0 ] += (double)x * cx;
117  m->m1()[ 1 ] += (double)y * cy;
118  m->m1()[ 2 ] += (double)z * cz;
119  }
120 
121  double mx = m->m1()[ 0 ] / m->m0();
122  double my = m->m1()[ 1 ] / m->m0();
123  double mz = m->m1()[ 2 ] / m->m0();
124 
125  m->gravity()[ 0 ] = mx;
126  m->gravity()[ 1 ] = my;
127  m->gravity()[ 2 ] = mz;
128 
129  it = d.begin() + d.oFirstPoint();
130  for ( z=0; z<dz; z++, it+=olbs )
131  for ( y=0; y<dy; y++, it+=opbl )
132  for ( x=0; x<dx; x++, it++ )
133  if ( *it == label )
134  addVoxel2( m, Point3d( x, y, z ) );
135 
136  m->m0() *= ct; // is is OK or should it have been done before ?
137 
138  for ( x=0; x<3; x++ ) m->m1()[ x ] *= ct;
139  for ( x=0; x<6; x++ ) m->m2()[ x ] *= ct;
140  for ( x=0; x<10; x++ ) m->m3()[ x ] *= ct;
141 }
142 
143 
144 template <typename T>
146  const aims::BucketMap<Void> & b, int )
147 {
148  const aims::BucketMap<Void>::Bucket & bk = b.begin()->second;
150 
151  if( bk.size() == 0 )
152  return;
153 
154  double cx = m->cx();
155  double cy = m->cy();
156  double cz = m->cz();
157  double ct = m->ct();
158 
159  m->sum() += bk.size();
160  m->m0() += bk.size();
161 
162  for( ib=bk.begin(); ib!=eb; ++ib )
163  {
164  m->m1()[ 0 ] += (double)ib->first[0] * cx;
165  m->m1()[ 1 ] += (double)ib->first[1] * cy;
166  m->m1()[ 2 ] += (double)ib->first[2] * cz;
167  }
168 
169  double mx = m->m1()[ 0 ] / m->m0();
170  double my = m->m1()[ 1 ] / m->m0();
171  double mz = m->m1()[ 2 ] / m->m0();
172 
173  m->gravity()[ 0 ] = mx;
174  m->gravity()[ 1 ] = my;
175  m->gravity()[ 2 ] = mz;
176 
177  for( ib=bk.begin(); ib!=eb; ++ib )
178  addVoxel2( m, ib->first );
179 
180  m->m0() *= ct;
181 
182  int x;
183  for ( x=0; x<3; x++ ) m->m1()[ x ] *= ct;
184  for ( x=0; x<6; x++ ) m->m2()[ x ] *= ct;
185  for ( x=0; x<10; x++ ) m->m3()[ x ] *= ct;
186 }
187 
188 #endif
double * m3()
Definition: moment.h:119
const T * const_iterator
double cz() const
Definition: moment.h:100
int dimX() const
iterator begin()
double cx() const
Definition: moment.h:98
void doit(Moment< T > *, AimsData< T > &, T, int)
Definition: momNormStgy.h:91
double & sum()
Definition: moment.h:110
double * m1()
Definition: moment.h:115
double & m0()
Definition: moment.h:112
Definition: moment.h:47
int oPointBetweenLine() const
double cy() const
Definition: moment.h:99
std::map< Point3d, T, BucketMapLess > Bucket
double ct() const
Definition: moment.h:101
int oLineBetweenSlice() const
double * m2()
Definition: moment.h:117
double * gravity()
Definition: moment.h:107
int dimZ() const
int dimY() const
int oFirstPoint() const