A.I.M.S algorithms


ffd.h
Go to the documentation of this file.
1 /* Copyright (C) 2000-2013 CEA
2  *
3  * This software and supporting documentation were developed by
4  * bioPICSEL
5  * CEA/DSV/I²BM/MIRCen/LMN, Batiment 61,
6  * 18, route du Panorama
7  * 92265 Fontenay-aux-Roses
8  * France
9  */
10 
11 
12 #ifndef AIMS_REGISTRATION_FFD_H
13 #define AIMS_REGISTRATION_FFD_H
14 
15 #include <aims/data/data.h> // AimsData
16 #include <aims/ffd/tabulSpline.h> // tabulSpline
19 #include <aims/transformation/transformation.h> // aims::Transformatin
20 #include <aims/utility/channel.h>
21 #include <aims/vector/vector.h> // Point*
22 #include <limits>
23 #include <string>
24 
25 namespace bio {
26 
27  //==========================================================================
28  // FFD TRANSFORMATION
29  //==========================================================================
30 
31  class SplineFfd:
33  {
34  public:
35  //--- Constructor --------------------------------------------------------
36  SplineFfd( int dimX = 0, int dimY = 1, int dimZ = 1,
37  float sizeX = 1., float sizeY = 1., float sizeZ = 1. );
38  template <typename T>
39  SplineFfd( int dimX, int dimY, int dimZ,
40  const AimsData<T> & test_volume );
41 
42  SplineFfd( const SplineFfd & other );
43  SplineFfd( const AimsData<Point3df> & other );
44  SplineFfd & operator=( const SplineFfd & other );
45 
46  operator const AimsData<Point3df>&() const { return _ctrlPointDelta; }
47  operator AimsData<Point3df>&() { return _ctrlPointDelta; }
48 
49  bool isIdentity() const
50  {
51  for( int z = 0; z < dimZ(); ++z )
52  for( int y = 0; y < dimY(); ++y )
53  for( int x = 0; x < dimX(); ++x )
54  {
55  if( _ctrlPointDelta(x, y, z) != Point3df(0., 0., 0.) )
56  return false;
57  }
58  return true;
59  }
60 
61  //--- Control Knots ------------------------------------------------------
62  Point3df getCtrlKnot( int nx, int ny, int nz ) const;
63  void updateCtrlKnot( int nx, int ny, int nz, const Point3df & newCtrlKnot );
64  void updateAllCtrlKnot( const AimsData<Point3df> & newCtrlKnotGrid );
65  void updateAllCtrlKnotFromDeformation( const AimsData<Point3df> & newDeformationGrid );
66 
67  //--- Modify -------------------------------------------------------------
68  void increaseResolution( const Point3d & addKnots );
69  void inverseTransform();
70  void estimateLocalDisplacement( const Point3df & voxelSize);
71 
72  //--- Deformation --------------------------------------------------------
73  Point3dd deformation( const Point3dd& p_mm ) const;
74  Point3dd ffdCoord( const Point3dd& p_mm ) const
75  { return mmToSplineVox(p_mm); }
76  double spline3( double x ) const { return _spline.spline3(x); }
77 
78  //--- Parameters ---------------------------------------------------------
79  int dimX() const { return _ctrlPointDelta.dimX(); }
80  int dimY() const { return _ctrlPointDelta.dimY(); }
81  int dimZ() const { return _ctrlPointDelta.dimZ(); }
82  float sizeX() const { return _ctrlPointDelta.sizeX(); }
83  float sizeY() const { return _ctrlPointDelta.sizeY(); }
84  float sizeZ() const { return _ctrlPointDelta.sizeZ(); }
85  bool isFlat( int i ) const
86  {
87  switch(i)
88  {
89  case 0: return dimX() == 1;
90  case 1: return dimY() == 1;
91  case 2: return dimZ() == 1;
92  default: return false;
93  }
94  }
95 
96  //--- Debug --------------------------------------------------------------
97  void printControlPointsGrid() const;
98  void writeDebugCtrlKnots( const std::string & filename ) const;
99  void writeDebugDeformations( const std::string & filename,
100  int dimX, int dimY, int dimZ,
101  float sizeX, float sizeY, float sizeZ ) const;
102 
103  //--- Output -------------------------------------------------------------
104  void write( const std::string & filename ) const;
105 
106  private:
107  //--- Protected methods --------------------------------------------------
108  Point3dd splineVoxToMm( const Point3dd& p ) const;
109  Point3dd mmToSplineVox( const Point3dd& p ) const;
110  void updateGridResolution( const AimsData<Point3df> & newGrid );
111  Point3dd transformDouble( double x, double y, double z ) const;
112 
113  AimsData<Point3df> _ctrlPointDelta;
114  TabulSpline _spline;
115  };
116 
117  template <typename T>
118  SplineFfd::SplineFfd( int dimX, int dimY, int dimZ,
119  const AimsData<T> & test_volume ):
120  _spline("SplineFFD"),
121  _ctrlPointDelta( dimX, dimY, dimZ )
122  {
123  _ctrlPointDelta = Point3df(0., 0., 0.);
124  _ctrlPointDelta.setSizeXYZT(
125  double(test_volume.dimX() - 1) / double(dimX - 1) * test_volume.sizeX(),
126  double(test_volume.dimY() - 1) / double(dimY - 1) * test_volume.sizeY(),
127  double(test_volume.dimZ() - 1) / double(dimZ - 1) * test_volume.sizeZ()
128  );
129  }
130 
131  inline Point3dd
132  SplineFfd::splineVoxToMm( const Point3dd & p ) const
133  {
134  return Point3dd( p[0] * sizeX(),
135  p[1] * sizeY(),
136  p[2] * sizeZ() );
137  }
138 
139  inline Point3dd
140  SplineFfd::mmToSplineVox( const Point3dd & p ) const
141  {
142  return Point3dd( p[0] / sizeX(),
143  p[1] / sizeY(),
144  p[2] / sizeZ() );
145  }
146 } // namespace bio
147 
148 
149 //============================================================================
150 // FFD READER/WRITER
151 //============================================================================
152 #include <aims/io/io_g.h> // aims::Reader / Writer
153 namespace aims
154 {
155 
156  template <>
157  class Reader<bio::SplineFfd>: public Reader<AimsData<Point3df> >
158  {
160  public:
161  Reader(): base() {}
162  Reader( const std::string& filename ): base(filename) {}
163  virtual ~Reader() {}
164  virtual bool read( bio::SplineFfd & obj,
165  int border=0,
166  const std::string* format = 0,
167  int frame = -1 )
168  {
169  bool res = base::read(obj, border, format, frame);
170  return res;
171  }
172  };
173 
174  template<>
175  class Writer<bio::SplineFfd> : public Writer<AimsData<Point3df> >
176  {
178  public:
179  inline Writer(): base() {}
180  inline Writer( const std::string& filename,
182  base( filename, options ) {}
183  virtual ~Writer() {}
184  virtual bool write( const bio::SplineFfd & obj,
185  bool ascii = false,
186  const std::string* format = 0 )
187  {
188  return base::write(obj, ascii, format);
189  }
190  };
191 
192 } // namespace aims
193 
194 //============================================================================
195 // R E S A M P L I N G
196 //============================================================================
197 namespace bio {
198 
205  template <class T>
207  {
208  public:
209 
210  virtual void init() = 0;
211  virtual void setRef(const AimsData<T> & ref) = 0;
212  virtual Point3df resample( const Point3df & output_location,
213  T & output_value, int t = 0 ) = 0;
214  };
215 
216  template <class T, class C = T>
217  class SplineFfdResampler : public FfdResampler<T>,
218  public CubicResampler<C>
219  {
220  public:
221 
222  virtual void init();
223  SplineFfdResampler(const SplineFfd & spline, T background = (T)0);
224  SplineFfdResampler(const SplineFfd & spline, Motion affine, T background = (T)0);
225  virtual void setRef(const AimsData<T> & ref);
226  virtual Point3df resample( const Point3df & output_location,
227  T & output_value, int t = 0 );
228 
229  private:
230  void updateCoef( int t = 0 );
231  const Motion _affine;
232  const SplineFfd & _transformation;
233  AimsData<T> _ref;
234  std::vector<AimsData<double> > _channelcoef;
235  ChannelSelector<T,C> _channelselector;
236 
237  T _background;
238  int _samples;
239  int _last_t;
240  std::vector<C> _min;
241  std::vector<C> _max;
242  };
243 
244  template <class T, class C = T>
246  public NearestNeighborResampler<C>
247  {
248  public:
249 
250  virtual void init();
251  NearestNeighborFfdResampler(const SplineFfd & spline, T background = (T)0);
252  NearestNeighborFfdResampler(const SplineFfd & spline, Motion affine, T background = (T)0);
253  virtual void setRef(const AimsData<T> & ref);
254  virtual Point3df resample( const Point3df & output_location,
255  T & output_value, int t = 0 );
256 
257  private:
258  const Motion _affine;
259  const SplineFfd & _transformation;
260  AimsData<T> _ref;
261 
262  T _background;
263  int _samples;
264  };
265 
266 } // namespace bio
267 
268 #endif
virtual Point3df resample(const Point3df &output_location, T &output_value, int t=0)
int dimY() const
Definition: ffd.h:80
int dimX() const
Definition: ffd.h:79
Point3df getCtrlKnot(int nx, int ny, int nz) const
NearestNeighborFfdResampler(const SplineFfd &spline, T background=(T) 0)
virtual bool write(const bio::SplineFfd &obj, bool ascii=false, const std::string *format=0)
Definition: ffd.h:184
Point3dd ffdCoord(const Point3dd &p_mm) const
Definition: ffd.h:74
Resampling image using Free Form Deformation transformations.
Definition: ffd.h:206
virtual void setRef(const AimsData< T > &ref)=0
int dimX() const
float sizeY() const
float sizeY() const
Definition: ffd.h:83
void increaseResolution(const Point3d &addKnots)
void estimateLocalDisplacement(const Point3df &voxelSize)
void writeDebugCtrlKnots(const std::string &filename) const
SplineFfd & operator=(const SplineFfd &other)
Writer(const std::string &filename, carto::Object options=carto::none())
Definition: ffd.h:180
float spline3(double) const
Definition: tabulSpline.h:94
virtual Point3df resample(const Point3df &output_location, T &output_value, int t=0)
const carto::Object options() const
SplineFfdResampler(const SplineFfd &spline, T background=(T) 0)
virtual bool read(bio::SplineFfd &obj, int border=0, const std::string *format=0, int frame=-1)
Definition: ffd.h:164
void printControlPointsGrid() const
void updateAllCtrlKnotFromDeformation(const AimsData< Point3df > &newDeformationGrid)
void updateAllCtrlKnot(const AimsData< Point3df > &newCtrlKnotGrid)
void writeDebugDeformations(const std::string &filename, int dimX, int dimY, int dimZ, float sizeX, float sizeY, float sizeZ) const
double spline3(double x) const
Definition: ffd.h:76
void write(const std::string &filename) const
float sizeX() const
Definition: ffd.h:82
float sizeZ() const
Definition: ffd.h:84
virtual void setRef(const AimsData< T > &ref)
Reader(const std::string &filename)
Definition: ffd.h:162
virtual void init()=0
float sizeZ() const
SplineFfd(int dimX=0, int dimY=1, int dimZ=1, float sizeX=1., float sizeY=1., float sizeZ=1.)
virtual Point3df resample(const Point3df &output_location, T &output_value, int t=0)=0
float sizeX() const
void inverseTransform()
int dimZ() const
int dimY() const
virtual void setRef(const AimsData< T > &ref)
bool isFlat(int i) const
Definition: ffd.h:85
void updateCtrlKnot(int nx, int ny, int nz, const Point3df &newCtrlKnot)
Object none()
bool isIdentity() const
Definition: ffd.h:49
Point3dd deformation(const Point3dd &p_mm) const
int dimZ() const
Definition: ffd.h:81
reference_wrapper< T > ref(T &ref)