]> Creatis software - clitk.git/blob - registration/clitkTransformToDeformationFieldSource.txx
small bugs corrections
[clitk.git] / registration / clitkTransformToDeformationFieldSource.txx
1 /*=========================================================================
2   Program:   vv                     http://www.creatis.insa-lyon.fr/rio/vv
3
4   Authors belong to: 
5   - University of LYON              http://www.universite-lyon.fr/
6   - Léon Bérard cancer center       http://oncora1.lyon.fnclcc.fr
7   - CREATIS CNRS laboratory         http://www.creatis.insa-lyon.fr
8
9   This software is distributed WITHOUT ANY WARRANTY; without even
10   the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
11   PURPOSE.  See the copyright notices for more information.
12
13   It is distributed under dual licence
14
15   - BSD        See included LICENSE.txt file
16   - CeCILL-B   http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
17 ======================================================================-====*/
18 #ifndef __clitkTransformToDeformationFieldSource_txx
19 #define __clitkTransformToDeformationFieldSource_txx
20 #include "clitkTransformToDeformationFieldSource.h"
21
22 #include "itkIdentityTransform.h"
23 #include "itkProgressReporter.h"
24 #include "itkImageRegionIteratorWithIndex.h"
25 #include "itkImageLinearIteratorWithIndex.h"
26
27 namespace clitk
28 {
29 /**
30  * Constructor
31  */
32   template <class TOutputImage, class TTransformPrecisionType>
33   TransformToDeformationFieldSource<TOutputImage, TTransformPrecisionType>::TransformToDeformationFieldSource()
34 {
35   this->m_OutputSpacing.Fill(1.0);
36   this->m_OutputOrigin.Fill(0.0);
37   this->m_OutputDirection.SetIdentity();
38
39   SizeType size;
40   size.Fill( 0 );
41   this->m_OutputRegion.SetSize( size );
42
43   IndexType index;
44   index.Fill( 0 );
45   this->m_OutputRegion.SetIndex( index );
46
47   this->m_Transform
48     = itk::IdentityTransform<TTransformPrecisionType, ImageDimension>::New();
49 } // end Constructor
50
51 /**
52  * Print out a description of self
53  *
54  * \todo Add details about this class
55  */
56 template <class TOutputImage, class TTransformPrecisionType>
57 void
58 TransformToDeformationFieldSource<TOutputImage, TTransformPrecisionType>
59 ::PrintSelf( std::ostream & os, itk::Indent indent ) const
60 {
61   Superclass::PrintSelf( os, indent );
62
63   os << indent << "OutputRegion: " << this->m_OutputRegion << std::endl;
64   os << indent << "OutputSpacing: " << this->m_OutputSpacing << std::endl;
65   os << indent << "OutputOrigin: " << this->m_OutputOrigin << std::endl;
66   os << indent << "OutputDirection: " << this->m_OutputDirection << std::endl;
67   os << indent << "Transform: " << this->m_Transform.GetPointer() << std::endl;
68 } // end PrintSelf()
69
70 /**
71  * Set the output image size.
72  */
73 template <class TOutputImage, class TTransformPrecisionType>
74 void
75 TransformToDeformationFieldSource<TOutputImage, TTransformPrecisionType>
76   ::SetOutputSize( const SizeType & size )
77 {
78   this->m_OutputRegion.SetSize( size );
79 }
80
81 /**
82  * Get the output image size.
83  */
84 template <class TOutputImage, class TTransformPrecisionType>
85 const typename TransformToDeformationFieldSource<TOutputImage,
86   TTransformPrecisionType>
87   ::SizeType &
88 TransformToDeformationFieldSource<TOutputImage, TTransformPrecisionType>
89   ::GetOutputSize()
90 {
91   return this->m_OutputRegion.GetSize();
92 }
93
94 /**
95  * Set the output image index.
96  */
97 template <class TOutputImage, class TTransformPrecisionType>
98 void
99 TransformToDeformationFieldSource<TOutputImage, TTransformPrecisionType>
100   ::SetOutputIndex( const IndexType & index )
101 {
102   this->m_OutputRegion.SetIndex( index );
103 }
104
105 /**
106  * Get the output image index.
107  */
108 template <class TOutputImage, class TTransformPrecisionType>
109 const typename TransformToDeformationFieldSource<TOutputImage,
110   TTransformPrecisionType>
111   ::IndexType &
112 TransformToDeformationFieldSource<TOutputImage, TTransformPrecisionType>
113   ::GetOutputIndex()
114 {
115   return this->m_OutputRegion.GetIndex();
116 }
117
118 /**
119  * Set the output image spacing.
120  */
121 template <class TOutputImage, class TTransformPrecisionType>
122 void
123 TransformToDeformationFieldSource<TOutputImage, TTransformPrecisionType>
124   ::SetOutputSpacing( const double *spacing )
125 {
126   SpacingType s( spacing );
127
128   this->SetOutputSpacing( s );
129 } // end SetOutputSpacing()
130
131 /**
132  * Set the output image origin.
133  */
134 template <class TOutputImage, class TTransformPrecisionType>
135 void
136 TransformToDeformationFieldSource<TOutputImage, TTransformPrecisionType>
137   ::SetOutputOrigin( const double *origin )
138 {
139   OriginType p( origin );
140
141   this->SetOutputOrigin( p );
142 }
143
144 /** Helper method to set the output parameters based on this image */
145 template <class TOutputImage, class TTransformPrecisionType>
146 void
147 TransformToDeformationFieldSource<TOutputImage, TTransformPrecisionType>
148   ::SetOutputParametersFromImage ( const ImageBaseType *image )
149 {
150   if ( !image )
151     {
152     itkExceptionMacro(<< "Cannot use a null image reference");
153     }
154
155   this->SetOutputOrigin( image->GetOrigin() );
156   this->SetOutputSpacing( image->GetSpacing() );
157   this->SetOutputDirection( image->GetDirection() );
158   this->SetOutputRegion( image->GetLargestPossibleRegion() );
159 } // end SetOutputParametersFromImage()
160
161 /**
162  * Set up state of filter before multi-threading.
163  * InterpolatorType::SetInputImage is not thread-safe and hence
164  * has to be set up before ThreadedGenerateData
165  */
166 template <class TOutputImage, class TTransformPrecisionType>
167 void
168 TransformToDeformationFieldSource<TOutputImage, TTransformPrecisionType>
169   ::BeforeThreadedGenerateData( void )
170 {
171   if ( !this->m_Transform )
172     {
173     itkExceptionMacro(<< "Transform not set");
174     }
175 } // end BeforeThreadedGenerateData()
176
177 /**
178  * ThreadedGenerateData
179  */
180 template <class TOutputImage, class TTransformPrecisionType>
181 void
182 TransformToDeformationFieldSource<TOutputImage, TTransformPrecisionType>
183   ::ThreadedGenerateData(
184   const OutputImageRegionType & outputRegionForThread,
185   int threadId )
186 {
187   // Check whether we can use a fast path for resampling. Fast path
188   // can be used if the transformation is linear. Transform respond
189   // to the IsLinear() call.
190   if ( this->m_Transform->IsLinear() )
191     {
192     this->LinearThreadedGenerateData( outputRegionForThread, threadId );
193     return;
194     }
195
196   // Otherwise, we use the normal method where the transform is called
197   // for computing the transformation of every point.
198   this->NonlinearThreadedGenerateData( outputRegionForThread, threadId );
199 } // end ThreadedGenerateData()
200
201 template <class TOutputImage, class TTransformPrecisionType>
202 void
203 TransformToDeformationFieldSource<TOutputImage, TTransformPrecisionType>
204   ::NonlinearThreadedGenerateData(
205   const OutputImageRegionType & outputRegionForThread,
206   int threadId )
207 {
208   // Get the output pointer
209   OutputImagePointer outputPtr = this->GetOutput();
210
211   // Create an iterator that will walk the output region for this thread.
212   typedef itk::ImageRegionIteratorWithIndex<TOutputImage> OutputIteratorType;
213   OutputIteratorType outIt( outputPtr, outputRegionForThread );
214
215   // Define a few variables that will be used to translate from an input pixel
216   // to an output pixel
217   PointType outputPoint;         // Coordinates of output pixel
218   PointType transformedPoint;    // Coordinates of transformed pixel
219   PixelType deformation;         // the difference
220
221   // Support for progress methods/callbacks
222   itk::ProgressReporter progress( this, threadId,
223                             outputRegionForThread.GetNumberOfPixels() );
224
225   // Walk the output region
226   outIt.GoToBegin();
227   while ( !outIt.IsAtEnd() )
228     {
229     // Determine the index of the current output pixel
230     outputPtr->TransformIndexToPhysicalPoint( outIt.GetIndex(), outputPoint );
231
232     // Compute corresponding input pixel position
233     transformedPoint = this->m_Transform->TransformPoint( outputPoint );
234
235     // Compute the deformation
236     for ( unsigned int i = 0; i < SpaceDimension; ++i )
237       {
238       deformation[i] = static_cast<PixelValueType>(
239         transformedPoint[i] - outputPoint[i] );
240       }
241
242     // Set it
243     outIt.Set( deformation );
244
245     // Update progress and iterator
246     progress.CompletedPixel();
247     ++outIt;
248     }
249 } // end NonlinearThreadedGenerateData()
250
251 template <class TOutputImage, class TTransformPrecisionType>
252 void
253 TransformToDeformationFieldSource<TOutputImage, TTransformPrecisionType>
254   ::LinearThreadedGenerateData(
255   const OutputImageRegionType & outputRegionForThread,
256   int threadId )
257 {
258   // Get the output pointer
259   OutputImagePointer outputPtr = this->GetOutput();
260
261   // Create an iterator that will walk the output region for this thread.
262   typedef itk::ImageLinearIteratorWithIndex<TOutputImage> OutputIteratorType;
263   OutputIteratorType outIt( outputPtr, outputRegionForThread );
264
265   outIt.SetDirection( 0 );
266
267   // Define a few indices that will be used to translate from an input pixel
268   // to an output pixel
269   PointType outputPoint;         // Coordinates of current output pixel
270   PointType transformedPoint;    // Coordinates of transformed pixel
271   PixelType deformation;         // the difference
272
273   IndexType index;
274
275   // Support for progress methods/callbacks
276   itk::ProgressReporter progress( this, threadId,
277                             outputRegionForThread.GetNumberOfPixels() );
278
279   // Determine the position of the first pixel in the scanline
280   outIt.GoToBegin();
281   index = outIt.GetIndex();
282   outputPtr->TransformIndexToPhysicalPoint( index, outputPoint );
283
284   // Compute corresponding transformed pixel position
285   transformedPoint = this->m_Transform->TransformPoint( outputPoint );
286
287   // Compare with the ResampleImageFilter
288
289   // Compute delta
290   PointType outputPointNeighbour;
291   PointType transformedPointNeighbour;
292   typedef typename PointType::VectorType VectorType;
293   VectorType delta;
294   ++index[0];
295   outputPtr->TransformIndexToPhysicalPoint( index, outputPointNeighbour );
296   transformedPointNeighbour = this->m_Transform->TransformPoint(
297     outputPointNeighbour );
298   delta = transformedPointNeighbour - transformedPoint
299           - ( outputPointNeighbour - outputPoint );
300
301   // loop over the vector image
302   while ( !outIt.IsAtEnd() )
303     {
304     // Get current point
305     index = outIt.GetIndex();
306     outputPtr->TransformIndexToPhysicalPoint( index, outputPoint );
307
308     // Compute transformed point
309     transformedPoint = this->m_Transform->TransformPoint( outputPoint );
310
311     while ( !outIt.IsAtEndOfLine() )
312       {
313       // Compute the deformation
314       for ( unsigned int i = 0; i < SpaceDimension; ++i )
315         {
316         deformation[i] = static_cast<PixelValueType>(
317           transformedPoint[i] - outputPoint[i] );
318         }
319
320       // Set it
321       outIt.Set( deformation );
322
323       // Update stuff
324       progress.CompletedPixel();
325       ++outIt;
326       transformedPoint += delta;
327       }
328
329     outIt.NextLine();
330     }
331 } // end LinearThreadedGenerateData()
332
333 /**
334  * Inform pipeline of required output region
335  */
336 template <class TOutputImage, class TTransformPrecisionType>
337 void
338 TransformToDeformationFieldSource<TOutputImage, TTransformPrecisionType>
339   ::GenerateOutputInformation( void )
340 {
341   // call the superclass' implementation of this method
342   Superclass::GenerateOutputInformation();
343
344   // get pointer to the output
345   OutputImagePointer outputPtr = this->GetOutput();
346   if ( !outputPtr )
347     {
348     return;
349     }
350
351   outputPtr->SetLargestPossibleRegion( m_OutputRegion );
352
353   outputPtr->SetSpacing( m_OutputSpacing );
354   outputPtr->SetOrigin( m_OutputOrigin );
355   outputPtr->SetDirection( m_OutputDirection );
356 } // end GenerateOutputInformation()
357
358 /**
359  * Verify if any of the components has been modified.
360  */
361 template <class TOutputImage, class TTransformPrecisionType>
362 unsigned long
363 TransformToDeformationFieldSource<TOutputImage, TTransformPrecisionType>
364   ::GetMTime( void ) const
365 {
366   unsigned long latestTime = itk::Object::GetMTime();
367
368   if ( this->m_Transform )
369     {
370     if ( latestTime < this->m_Transform->GetMTime() )
371       {
372       latestTime = this->m_Transform->GetMTime();
373       }
374     }
375
376   return latestTime;
377 } // end GetMTime()
378 } // end namespace clitk
379
380 #endif // end #ifndef _clitkTransformToDeformationFieldSource_txx