A.I.M.S algorithms


sampler.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_RESAMPLING_SAMPLER_H
36 #define AIMS_RESAMPLING_SAMPLER_H
37 
38 #include <aims/def/assert.h>
39 #include <aims/data/data.h>
40 #include <aims/vector/vector.h>
41 #include <aims/resampling/motion.h>
42 #include <aims/bucket/bucketMap.h>
43 #include <vector>
44 
45 
46 typedef struct {
47  unsigned short x, y, z;
48  int offset;
49 } PVItem;
50 
51 
52 namespace carto
53 {
54  extern template class Volume<PVItem>;
55 }
56 
57 
58 template <class T>
59 class Sampler
60 {
61  public:
62 
63  Sampler( ) : _ref( 0 ) { };
64  virtual ~Sampler() { }
65 
66  virtual void doit(const Motion& motion, AimsData<PVItem>& thing ) const;
67 
68 
69  void setRef( const AimsData<T>& ref ) { _ref = &ref; }
70 
71  protected:
72 
73  void _sliceResamp( AimsData<PVItem>& resamp,
74  PVItem* out,
75  const Point3df& start, int t,
76  const AimsData<float>& Rinv ) const;
77  const AimsData<T>* _ref;
78 };
79 
80 
81 typedef struct {
82  unsigned short x, y, z ;
83  int offset ;
84 } PVVectorItem;
85 
86 template <class T>
88 {
89  public:
90 
91  BucketSampler( ) : _ref(0) { }
92  virtual ~BucketSampler() { }
93 
94  virtual void doit( const Motion& motion,
95  std::vector< PVVectorItem >& thing,
96  const aims::BucketMap<T>& model ) const;
97 
98  void setRef( const AimsData<T>& ref ) { _ref = &ref; }
99 
100  protected:
102 };
103 
104 
105 
106 
107 template <class T> inline
108 void
109 Sampler<T>::doit( const Motion& motion, AimsData<PVItem>& thing ) const
110 {
111 
112 
113  ASSERT( _ref );
114  ASSERT( thing.dimT() == _ref->dimT() && thing.borderWidth() == 0 );
115 
116  Motion dirMotion = motion;
117  Motion invMotion = motion.inverse();
118 
119  // scale motion
120  Point3df sizeFrom( _ref->sizeX(), _ref->sizeY(), _ref->sizeZ() );
121  Point3df sizeTo( thing.sizeX(), thing.sizeY(), thing.sizeZ() );
122  dirMotion.scale( sizeFrom, sizeTo );
123  invMotion.scale( sizeTo, sizeFrom );
124 
125 
126  AimsData<PVItem>::iterator it = thing.begin() + thing.oFirstPoint();
127 
128 #ifdef DEBUG
129  cout << "Sampling volume [ 1] / slice [ 1]" << flush;
130 #endif
131  for (int t = 1; t <= thing.dimT(); t++ )
132  {
133  Point3df start = invMotion.translation();
134  for ( int s = 1; s <= thing.dimZ(); s++ )
135  {
136 #ifdef DEBUG
137  cout << "\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b"
138  << setw( 3 ) << t << "] / slice ["
139  << setw( 3 ) << s << "]" << flush;
140 #endif
141  _sliceResamp( thing, it, start, t - 1, invMotion.rotation() );
142  it += thing.oSlice();
143  start[ 0 ] += invMotion.rotation()( 0, 2 );
144  start[ 1 ] += invMotion.rotation()( 1, 2 );
145  start[ 2 ] += invMotion.rotation()( 2, 2 );
146  }
147  }
148  //cout << endl;
149 
150 
151 }
152 
153 
154 template <class T> inline
156  PVItem* out,
157  const Point3df& start, int t,
158  const AimsData<float>& Rinv )
159  const
160 {
161  int dimX1 = _ref->dimX();
162  int dimX2 = resamp.dimX();
163  int dimY2 = resamp.dimY();
164 
165  const int SIXTEEN = 16;
166  const int TWO_THEN_SIXTEEN = 65536;
167 
168  int maxX = TWO_THEN_SIXTEEN * ( dimX1 - 1 );
169  int maxY = TWO_THEN_SIXTEEN * ( _ref->dimY() - 1 );
170  int maxZ = TWO_THEN_SIXTEEN * ( _ref->dimZ() - 1 );
171 
172  int xLinCurrent = ( int )( 65536.0 * start[ 0 ] );
173  int yLinCurrent = ( int )( 65536.0 * start[ 1 ] );
174  int zLinCurrent = ( int )( 65536.0 * start[ 2 ] );
175 
176  int xu = ( int )( 65536.0 * Rinv( 0, 0 ) );
177  int yu = ( int )( 65536.0 * Rinv( 1, 0 ) );
178  int zu = ( int )( 65536.0 * Rinv( 2, 0 ) );
179  int xv = ( int )( 65536.0 * Rinv( 0, 1 ) );
180  int yv = ( int )( 65536.0 * Rinv( 1, 1 ) );
181  int zv = ( int )( 65536.0 * Rinv( 2, 1 ) );
182 
183  int xCurrent, yCurrent, zCurrent;
184  typename AimsData<T>::const_iterator it;
185  int incr_vois[8] = {0, 1,
186  dimX1+1,dimX1,
187  _ref->oSlice(),_ref->oSlice()+1,
188  _ref->oSlice()+ dimX1+1,_ref->oSlice() + dimX1};
189  for ( int v = dimY2; v--; )
190  {
191  xCurrent = xLinCurrent;
192  yCurrent = yLinCurrent;
193  zCurrent = zLinCurrent;
194 
195  for ( int u = dimX2; u--; )
196  {
197  if ( xCurrent >= 0 && xCurrent < maxX &&
198  yCurrent >= 0 && yCurrent < maxY &&
199  zCurrent >= 0 && zCurrent < maxZ )
200  {
201  it = _ref->begin() + _ref->oFirstPoint() + t * _ref->oVolume()
202  + ( zCurrent >> SIXTEEN ) * _ref->oSlice()
203  + ( yCurrent >> SIXTEEN ) * dimX1
204  + ( xCurrent >> SIXTEEN );
205 
206  if (*(it + incr_vois[0]) == -32768 ||
207  *(it + incr_vois[1]) == -32768 ||
208  *(it + incr_vois[2]) == -32768 ||
209  *(it + incr_vois[3]) == -32768 ||
210  *(it + incr_vois[4]) == -32768 ||
211  *(it + incr_vois[5]) == -32768 ||
212  *(it + incr_vois[6]) == -32768 ||
213  *(it + incr_vois[7]) == -32768 )
214  {
215  out->offset = -1;
216  }
217  else
218  {
219  out->x = xCurrent & 0xffff;
220  out->y = yCurrent & 0xffff;
221  out->z = zCurrent & 0xffff;
222  out->offset = ( zCurrent >> SIXTEEN ) * _ref->oSlice()
223  + ( yCurrent >> SIXTEEN ) * dimX1
224  + ( xCurrent >> SIXTEEN );
225  }
226  }
227  else
228  {
229  out->offset = -1;
230  }
231  out++;
232  xCurrent += xu;
233  yCurrent += yu;
234  zCurrent += zu;
235  }
236  xLinCurrent += xv;
237  yLinCurrent += yv;
238  zLinCurrent += zv;
239  }
240 }
241 
242 template <class T> inline
243 void
245  std::vector< PVVectorItem >& thing,
246  const aims::BucketMap<T>& model ) const
247 {
248  const int TWO_THEN_SIXTEEN = 65536;
249 
250 
251  thing.clear() ;
252 
253  ASSERT( _ref );
254 
255  Motion dirMotion; dirMotion = motion;
256  Motion invMotion; invMotion = motion.inverse();
257 
258 
259  // scale motion
260  Point3df sizeFrom( _ref->sizeX(), _ref->sizeY(), _ref->sizeZ() );
261  Point3df sizeTo( model.sizeX(), model.sizeY(), model.sizeZ() );
262 
263  dirMotion.scale( sizeFrom, sizeTo );
264  invMotion.scale( sizeTo, sizeFrom );
265 
267  it = model.begin()->second.begin() ;
268 
269  PVVectorItem item ;
270  int nbOfVoxels = model.begin()->second.size() ;
271 
272  thing.reserve(nbOfVoxels) ;
273  thing.resize(nbOfVoxels) ;
274 
275  Point3df x, y, z;
276 
277  float transl0 = invMotion.translation()[0];
278  float transl1 = invMotion.translation()[1];
279  float transl2 = invMotion.translation()[2];
280 
281  invMotion.translation() = 0. ;
282 
283  x = invMotion.transform( 1.F, 0.F, 0.F ) ;
284  y = invMotion.transform( 0.F, 1.F, 0.F ) ;
285  z = invMotion.transform( 0.F, 0.F, 1.F ) ;
286 
287  int maxX = _ref->dimX() - 1;
288  int maxY = _ref->dimY() - 1;
289  int maxZ = _ref->dimZ() - 1;
290  int oS= _ref->oSlice(), dX = _ref->dimX();
291 
292  for(int i = 0 ; i < nbOfVoxels ; ++i , ++it )
293  {
294  float current0, current1, current2;
295  current0 = static_cast<float>( (it->first)[0] ) ;
296  current1 = static_cast<float>( (it->first)[1] ) ;
297  current2 = static_cast<float>( (it->first)[2] ) ;
298 
299  float tmp0 = current0 * x[0] + current1 * y[0] + current2 * z[0] + transl0 ;
300  float tmp1 = current0 * x[1] + current1 * y[1] + current2 * z[1] + transl1 ;
301  float tmp2 = current0 * x[2] + current1 * y[2] + current2 * z[2] + transl2 ;
302 
303 
304  if ( tmp0 >= 0 && tmp0 < maxX &&
305  tmp1 >= 0 && tmp1 < maxY &&
306  tmp2 >= 0 && tmp2 < maxZ)
307  {
308  item.x = (unsigned short) (( tmp0 - int(tmp0) ) * TWO_THEN_SIXTEEN);
309  item.y = (unsigned short) (( tmp1 - int(tmp1) ) * TWO_THEN_SIXTEEN);
310  item.z = (unsigned short) (( tmp2 - int(tmp2) ) * TWO_THEN_SIXTEEN);
311  item.offset = int(tmp2) * oS +
312  int(tmp1) * dX +
313  int(tmp0);
314  }
315  else
316  {
317  item.offset = -1;
318  }
319 
320  thing[i] = item ;
321  }
322 }
323 
324 #endif
virtual void scale(const Point3df &sizeFrom, const Point3df &sizeTo)
const T * const_iterator
Sampler()
Definition: sampler.h:63
unsigned short x
Definition: sampler.h:82
void setRef(const AimsData< T > &ref)
Definition: sampler.h:69
virtual ~BucketSampler()
Definition: sampler.h:92
int dimX() const
float sizeY() const
iterator begin()
Definition: sampler.h:46
int dimT() const
float sizeX() const
const AimsData< T > * _ref
Definition: sampler.h:77
BucketSampler()
Definition: sampler.h:91
int borderWidth() const
float sizeZ() const
Point3dd transform(double x, double y, double z) const
int offset
Definition: sampler.h:48
void _sliceResamp(AimsData< PVItem > &resamp, PVItem *out, const Point3df &start, int t, const AimsData< float > &Rinv) const
Definition: sampler.h:155
virtual void doit(const Motion &motion, std::vector< PVVectorItem > &thing, const aims::BucketMap< T > &model) const
Definition: sampler.h:244
AimsData< float > & rotation()
float sizeZ() const
unsigned short y
Definition: sampler.h:47
unsigned short z
Definition: sampler.h:47
float sizeX() const
float sizeY() const
void setRef(const AimsData< T > &ref)
Definition: sampler.h:98
int dimZ() const
unsigned short z
Definition: sampler.h:82
virtual void doit(const Motion &motion, AimsData< PVItem > &thing) const
Definition: sampler.h:109
int dimY() const
int offset
Definition: sampler.h:83
virtual ~Sampler()
Definition: sampler.h:64
unsigned short y
Definition: sampler.h:82
unsigned short x
Definition: sampler.h:47
#define ASSERT(EX)
int oSlice() const
int oFirstPoint() const
AffineTransformation3d inverse() const
const AimsData< T > * _ref
Definition: sampler.h:101