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