A.I.M.S algorithms


masklinresampler.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_MASKLINRESAMPLER_H
36 #define AIMS_RESAMPLING_MASKLINRESAMPLER_H
37 
38 #include <aims/resampling/motion.h>
40 #include <aims/data/data.h>
41 #include <aims/vector/vector.h>
42 
43 
44 template <class T>
45 class MaskLinearResampler : public Resampler< T >
46 {
47  public:
48 
51 
52  void doit( const Motion& motion, AimsData<T>& thing );
53  AimsData<T> doit( const Motion& motion,int dimX, int dimY, int dimZ,
54  const Point3df& resolution );
55 
56  protected:
57 
58  void _sliceResamp( AimsData<T>& resamp, T* out,
59  const Point3df& start,
60  int t, const AimsData<float>& Rinv );
61 
62  void
63  doResample( const AimsData< T > &, const Motion &, const T &,
64  const Point3df &, T &, int ) {}
65 };
66 
67 
68 template <class T> inline
71  int dimX, int dimY, int dimZ,
72  const Point3df& resolution )
73 {
74  ASSERT( this->_ref );
75 
76  AimsData<T> thing( dimX, dimY, dimZ, this->_ref->dimT() );
77  thing.setSizeXYZT( resolution[ 0 ], resolution[ 1 ], resolution[ 2 ],
78  this->_ref->sizeT() );
79 
80  doit( motion, thing );
81 
82  return thing;
83 }
84 
85 
86 template <class T> inline
88  AimsData<T>& thing )
89 {
90  ASSERT( this->_ref );
91  ASSERT( thing.dimT() == this->_ref->dimT() && thing.borderWidth() == 0 );
92 
93  Motion dirMotion = motion;
94  Motion invMotion = motion.inverse();
95 
96  // scale motion
97  Point3df sizeFrom( this->_ref->sizeX(), this->_ref->sizeY(), this->_ref->sizeZ() );
98  Point3df sizeTo( thing.sizeX(), thing.sizeY(), thing.sizeZ() );
99  dirMotion.scale( sizeFrom, sizeTo );
100  invMotion.scale( sizeTo, sizeFrom );
101 
102  //
103  Point3df start;
104  typename AimsData<T>::iterator it = thing.begin() + thing.oFirstPoint();
105 
106 #ifdef DEBUG
107  std::cout << "Resampling volume [ 1] / slice [ 1]" << std::flush;
108 #endif
109  for (int t = 1; t <= thing.dimT(); t++ )
110  {
111  start = invMotion.translation();
112  for ( int s = 1; s <= thing.dimZ(); s++ )
113  {
114 #ifdef DEBUG
115  std::cout << "\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b"
116  << std::setw( 3 ) << t << "] / slice ["
117  << std::setw( 3 ) << s << "]" << std::flush;
118 #endif
119  _sliceResamp( thing, it, start, t - 1, invMotion.rotation() );
120  it += thing.oSlice();
121  start[ 0 ] += invMotion.rotation()( 0, 2 );
122  start[ 1 ] += invMotion.rotation()( 1, 2 );
123  start[ 2 ] += invMotion.rotation()( 2, 2 );
124  }
125  }
126  std::cout << std::endl;
127 }
128 
129 
130 template <class T> inline
132  T* out,
133  const Point3df& start, int t,
134  const AimsData<float>& Rinv )
135 {
137  pOrig = this->_ref->begin() + this->_ref->oFirstPoint() +
138  t * this->_ref->oVolume();
139 
140  int dimX1 = this->_ref->dimX();
141  int dimX2 = resamp.dimX();
142  int dimY2 = resamp.dimY();
143 
144  const int SIXTEEN = 16;
145  const int TWO_THEN_SIXTEEN = 65536;
146 
147  int maxX = TWO_THEN_SIXTEEN * ( dimX1 - 1 );
148  int maxY = TWO_THEN_SIXTEEN * ( this->_ref->dimY() - 1 );
149  int maxZ = TWO_THEN_SIXTEEN * ( this->_ref->dimZ() - 1 );
150 
151  int xLinCurrent = ( int )( 65536.0 * start[ 0 ] );
152  int yLinCurrent = ( int )( 65536.0 * start[ 1 ] );
153  int zLinCurrent = ( int )( 65536.0 * start[ 2 ] );
154 
155  int xu = ( int )( 65536.0 * Rinv( 0, 0 ) );
156  int yu = ( int )( 65536.0 * Rinv( 1, 0 ) );
157  int zu = ( int )( 65536.0 * Rinv( 2, 0 ) );
158  int xv = ( int )( 65536.0 * Rinv( 0, 1 ) );
159  int yv = ( int )( 65536.0 * Rinv( 1, 1 ) );
160  int zv = ( int )( 65536.0 * Rinv( 2, 1 ) );
161 
162  int xCurrent, yCurrent, zCurrent;
163  int dx, dy, dz;
164  typename AimsData<T>::const_iterator it;
165  int incr_vois[8] = {0, 1,
166  dimX1+1,dimX1,
167  this->_ref->oSlice(),this->_ref->oSlice()+1,
168  this->_ref->oSlice()+ dimX1+1,this->_ref->oSlice() + dimX1};
169  float stock1, stock2, stock3;
170  for ( int v = dimY2; v--; )
171  {
172  xCurrent = xLinCurrent;
173  yCurrent = yLinCurrent;
174  zCurrent = zLinCurrent;
175 
176  for ( int u = dimX2; u--; )
177  {
178  if ( xCurrent >= 0 && xCurrent < maxX &&
179  yCurrent >= 0 && yCurrent < maxY &&
180  zCurrent >= 0 && zCurrent < maxZ )
181  {
182  dx = xCurrent & 0xffff;
183  dy = yCurrent & 0xffff;
184  dz = zCurrent & 0xffff;
185 
186  it = pOrig + ( zCurrent >> SIXTEEN ) * this->_ref->oSlice()
187  + ( yCurrent >> SIXTEEN ) * dimX1
188  + ( xCurrent >> SIXTEEN );
189 
190  if (*(it + incr_vois[0]) == -32768 ||
191  *(it + incr_vois[1]) == -32768 ||
192  *(it + incr_vois[2]) == -32768 ||
193  *(it + incr_vois[3]) == -32768 ||
194  *(it + incr_vois[4]) == -32768 ||
195  *(it + incr_vois[5]) == -32768 ||
196  *(it + incr_vois[6]) == -32768 ||
197  *(it + incr_vois[7]) == -32768 )
198  {
199  *out++ = (T) -32768;
200  }
201  else
202  {
203 
204 
205  stock1 = *it * ( TWO_THEN_SIXTEEN - dx );
206  stock1 += *(++it) * dx;
207  stock1 /= TWO_THEN_SIXTEEN;
208 
209  it += dimX1;
210  stock2 = *it * dx;
211  stock2 += *(--it) * ( TWO_THEN_SIXTEEN - dx );
212  stock2 /= TWO_THEN_SIXTEEN;
213 
214  stock1 *= ( TWO_THEN_SIXTEEN - dy );
215  stock1 += stock2 * dy;
216  stock1 /= TWO_THEN_SIXTEEN;
217 
218  it += this->_ref->oSlice() - dimX1;
219  stock2 = *it * ( TWO_THEN_SIXTEEN - dx );
220  stock2 += *(++it) * dx;
221  stock2 /= TWO_THEN_SIXTEEN;
222 
223  it += dimX1;
224  stock3 = *it * dx;
225  stock3 += *(--it) * ( TWO_THEN_SIXTEEN - dx );
226  stock3 /= TWO_THEN_SIXTEEN;
227 
228  stock2 *= TWO_THEN_SIXTEEN - dy;
229  stock2 += stock3 * dy;
230  stock2 /= TWO_THEN_SIXTEEN;
231 
232  stock1 *= ( TWO_THEN_SIXTEEN - dz );
233  stock1 += stock2 * dz;
234 
235  *out++ = ( T )( stock1 / TWO_THEN_SIXTEEN );
236  }
237  }
238  else
239  *out++ = ( T ) -32768;
240 
241  xCurrent += xu;
242  yCurrent += yu;
243  zCurrent += zu;
244  }
245  xLinCurrent += xv;
246  yLinCurrent += yv;
247  zLinCurrent += zv;
248  }
249 }
250 
251 #endif
virtual void scale(const Point3df &sizeFrom, const Point3df &sizeTo)
void doResample(const AimsData< T > &, const Motion &, const T &, const Point3df &, T &, int)
const T * const_iterator
int dimX() const
float sizeY() const
iterator begin()
T * iterator
Resampler resamples an input data to build or fill an output data, using an affine transformation...
Definition: resampler.h:52
int dimT() const
void _sliceResamp(AimsData< T > &resamp, T *out, const Point3df &start, int t, const AimsData< float > &Rinv)
void doit(const Motion &motion, AimsData< T > &thing)
Resample the reference input data (set via setRef()) into an existing output data.
int borderWidth() const
AimsData< float > & rotation()
float sizeZ() const
float sizeX() const
int dimZ() const
void setSizeXYZT(float sizex=1.0f, float sizey=1.0f, float sizez=1.0f, float sizet=1.0f)
int dimY() const
#define ASSERT(EX)
int oSlice() const
int oFirstPoint() const
AffineTransformation3d inverse() const