registration/clitkRecursiveSpatioTemporalMultiResolutionPyramidImageFilter.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://www.centreleonberard.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 __clitkRecursiveSpatioTemporalMultiResolutionPyramidImageFilter_txx
  19 #define __clitkRecursiveSpatioTemporalMultiResolutionPyramidImageFilter_txx
  20 #include "clitkRecursiveSpatioTemporalMultiResolutionPyramidImageFilter.h"
  21 #include "itkShrinkImageFilter.h"
  22 #include "itkGaussianOperator.h"
  23 #include "itkCastImageFilter.h"
  24 #include "itkDiscreteGaussianImageFilter.h"
  25 #include "itkExceptionObject.h"
  26 #include "itkResampleImageFilter.h"
  27 #include "itkShrinkImageFilter.h"
  28 #include "itkIdentityTransform.h"
  29
  30 namespace clitk
  31 {
  32
  33 /**
  34  * Constructor
  35  */
  36 template <class TInputImage, class TOutputImage>
  37 RecursiveSpatioTemporalMultiResolutionPyramidImageFilter<TInputImage, TOutputImage>
  38 ::RecursiveSpatioTemporalMultiResolutionPyramidImageFilter()
  39 {
  40   this->Superclass::m_UseShrinkImageFilter = true;
  41 }
  42
  43 /**
  44  * GenerateData
  45  */
  46 template <class TInputImage, class TOutputImage>
  47 void
  48 RecursiveSpatioTemporalMultiResolutionPyramidImageFilter<TInputImage, TOutputImage>
  49 ::GenerateData()
  50 {
  51
  52   if( !this->IsScheduleDownwardDivisible( this->GetSchedule() ) )
  53     {
  54     // use the Superclass implemenation
  55     this->Superclass::GenerateData();
  56     return;
  57     }
  58
  59   // Get the input and output pointers
  60   InputImageConstPointer  inputPtr = this->GetInput();
  61
  62   // Create caster, smoother and resampleShrink filters
  63   typedef itk::CastImageFilter<TInputImage, TOutputImage>              CasterType;
  64   typedef itk::CastImageFilter<TOutputImage, TOutputImage>             CopierType;
  65   typedef itk::DiscreteGaussianImageFilter<TOutputImage, TOutputImage> SmootherType;
  66
  67   typedef itk::ImageToImageFilter<TOutputImage,TOutputImage>           ImageToImageType;
  68   typedef itk::ResampleImageFilter<TOutputImage,TOutputImage>          ResampleShrinkerType;
  69   typedef itk::ShrinkImageFilter<TOutputImage,TOutputImage>            ShrinkerType;
  70
  71   typename CasterType::Pointer   caster   = CasterType::New();
  72   typename CopierType::Pointer   copier   = CopierType::New();
  73   typename SmootherType::Pointer smoother = SmootherType::New();
  74
  75
  76   typename ImageToImageType::Pointer shrinkerFilter;
  77   //
  78   // only one of these pointers is going to be valid, depending on the
  79   // value of UseShrinkImageFilter flag
  80   typename ResampleShrinkerType::Pointer resampleShrinker;
  81   typename ShrinkerType::Pointer shrinker;
  82
  83   if(this->GetUseShrinkImageFilter())
  84     {
  85     shrinker = ShrinkerType::New();
  86     shrinkerFilter = shrinker.GetPointer();
  87     }
  88   else
  89     {
  90     resampleShrinker = ResampleShrinkerType::New();
  91     typedef itk::LinearInterpolateImageFunction< OutputImageType, double >
  92       LinearInterpolatorType;
  93     typename LinearInterpolatorType::Pointer interpolator =
  94       LinearInterpolatorType::New();
  95     typedef itk::IdentityTransform<double,OutputImageType::ImageDimension>
  96       IdentityTransformType;
  97     typename IdentityTransformType::Pointer identityTransform =
  98       IdentityTransformType::New();
  99     resampleShrinker->SetInterpolator( interpolator );
 100     resampleShrinker->SetDefaultPixelValue( 0 );
 101     resampleShrinker->SetTransform(identityTransform);
 102     shrinkerFilter = resampleShrinker.GetPointer();
 103     }
 104
 105   int ilevel;
 106   unsigned int idim;
 107   unsigned int factors[ImageDimension];
 108   double       variance[ImageDimension];
 109
 110   bool         allOnes;
 111   OutputImagePointer   outputPtr;
 112   OutputImagePointer swapPtr;
 113   typename TOutputImage::RegionType  LPRegion;
 114
 115   smoother->SetUseImageSpacing( false );
 116   smoother->SetMaximumError( this->GetMaximumError() );
 117   shrinkerFilter->SetInput( smoother->GetOutput() );
 118
 119
 120   // recursively compute outputs starting from the last one
 121   for( ilevel = this->GetNumberOfLevels() - 1; ilevel > -1; ilevel--)
 122     {
 123
 124     this->UpdateProgress( 1.0 - static_cast<float>( 1 + ilevel ) /
 125                           static_cast<float>( this->GetNumberOfLevels() ) );
 126
 127     // Allocate memory for each output
 128     outputPtr = this->GetOutput( ilevel );
 129     outputPtr->SetBufferedRegion( outputPtr->GetRequestedRegion() );
 130     outputPtr->Allocate();
 131
 132     // cached a copy of the largest possible region
 133     LPRegion = outputPtr->GetLargestPossibleRegion();
 134
 135     // Check shrink factors and compute variances
 136     allOnes = true;
 137     for( idim = 0; idim < ImageDimension; idim++ )
 138       {
 139       if( ilevel == static_cast<int>(this->GetNumberOfLevels()) - 1)
 140         {
 141         factors[idim] = this->GetSchedule()[ilevel][idim];
 142         }
 143       else
 144         {
 145         factors[idim] = this->GetSchedule()[ilevel][idim] /
 146           this->GetSchedule()[ilevel+1][idim];
 147         }
 148       variance[idim] = std::sqr( 0.5 *
 149                                      static_cast<float>( factors[idim] ) );
 150       if( factors[idim] != 1 )
 151         {
 152         allOnes = false;
 153         }
 154       else
 155         {
 156         variance[idim] = 0.0;
 157         }
 158       }
 159
 160
 161     if( allOnes && ilevel == static_cast<int>(this->GetNumberOfLevels()) - 1 )
 162       {
 163       // just copy the input over
 164       caster->SetInput( inputPtr );
 165       caster->GraftOutput( outputPtr );
 166       // ensure only the requested region is updated
 167       caster->UpdateOutputInformation();
 168       caster->GetOutput()->SetRequestedRegion(outputPtr->GetRequestedRegion());
 169       caster->GetOutput()->PropagateRequestedRegion();
 170       caster->GetOutput()->UpdateOutputData();
 171
 172       swapPtr = caster->GetOutput();
 173
 174
 175       }
 176     else if( allOnes )
 177       {
 178       // just copy the data over
 179       copier->SetInput( swapPtr );
 180       copier->GraftOutput( outputPtr );
 181       // ensure only the requested region is updated
 182       copier->GetOutput()->UpdateOutputInformation();
 183       copier->GetOutput()->SetRequestedRegion(outputPtr->GetRequestedRegion());
 184       copier->GetOutput()->PropagateRequestedRegion();
 185       copier->GetOutput()->UpdateOutputData();
 186
 187       swapPtr = copier->GetOutput();
 188
 189       }
 190     else
 191       {
 192       if( ilevel == static_cast<int>(this->GetNumberOfLevels()) - 1 )
 193         {
 194         // use caster -> smoother -> shrinker piepline
 195         caster->SetInput( inputPtr );
 196         smoother->SetInput( caster->GetOutput() );
 197         }
 198       else
 199         {
 200         // use smoother -> shrinker pipeline
 201         smoother->SetInput( swapPtr );
 202         }
 203
 204       smoother->SetVariance( variance );
 205
 206       //      shrinker->SetShrinkFactors( factors );
 207       //      shrinker->GraftOutput( outputPtr );
 208       if(!this->GetUseShrinkImageFilter())
 209         {
 210         resampleShrinker->SetOutputParametersFromImage(outputPtr);
 211         }
 212       else
 213         {
 214         shrinker->SetShrinkFactors(factors);
 215         }
 216       shrinkerFilter->GraftOutput(outputPtr);
 217       shrinkerFilter->Modified();
 218       // ensure only the requested region is updated
 219       shrinkerFilter->GetOutput()->UpdateOutputInformation();
 220       shrinkerFilter->GetOutput()->SetRequestedRegion(outputPtr->GetRequestedRegion());
 221       shrinkerFilter->GetOutput()->PropagateRequestedRegion();
 222       shrinkerFilter->GetOutput()->UpdateOutputData();
 223
 224       swapPtr = shrinkerFilter->GetOutput();
 225
 226       }
 227
 228     // graft pipeline output back onto this filter's output
 229     swapPtr->SetLargestPossibleRegion( LPRegion );
 230     this->GraftNthOutput( ilevel, swapPtr );
 231
 232     // disconnect from pipeline to stop cycle
 233     swapPtr->DisconnectPipeline();
 234
 235     }
 236
 237 }
 238
 239 /**
 240  * PrintSelf method
 241  */
 242 template <class TInputImage, class TOutputImage>
 243 void
 244 RecursiveSpatioTemporalMultiResolutionPyramidImageFilter<TInputImage, TOutputImage>
 245 ::PrintSelf(std::ostream& os, itk::Indent indent) const
 246 {
 247   Superclass::PrintSelf(os,indent);
 248 }
 249
 250
 251 /*
 252  * GenerateOutputRequestedRegion
 253  */
 254 template <class TInputImage, class TOutputImage>
 255 void
 256 RecursiveSpatioTemporalMultiResolutionPyramidImageFilter<TInputImage, TOutputImage>
 257 ::GenerateOutputRequestedRegion(itk::DataObject * ptr )
 258 {
 259
 260   // call the superclass's implementation of this method
 261   Superclass::GenerateOutputRequestedRegion( ptr );
 262
 263   TOutputImage * refOutputPtr = static_cast<TOutputImage*>( ptr );
 264   if( !refOutputPtr )
 265     {
 266     itkExceptionMacro( << "Could not cast ptr to TOutputImage*." );
 267     }
 268
 269   // find the index for this output
 270   unsigned int refLevel;
 271   refLevel = refOutputPtr->GetSourceOutputIndex();
 272
 273   typedef typename TOutputImage::PixelType                 OutputPixelType;
 274   typedef itk::GaussianOperator<OutputPixelType,ImageDimension> OperatorType;
 275
 276   OperatorType * oper = new OperatorType;
 277   oper->SetMaximumError( this->GetMaximumError() );
 278
 279   typedef typename OutputImageType::SizeType    SizeType;
 280   typedef typename SizeType::SizeValueType      SizeValueType;
 281   typedef typename OutputImageType::IndexType   IndexType;
 282   typedef typename IndexType::IndexValueType    IndexValueType;
 283   typedef typename OutputImageType::RegionType  RegionType;
 284
 285   int ilevel, idim;
 286   unsigned int factors[ImageDimension];
 287
 288   typename TInputImage::SizeType radius;
 289
 290   RegionType requestedRegion;
 291   SizeType   requestedSize;
 292   IndexType  requestedIndex;
 293
 294   // compute requested regions for lower levels
 295   for( ilevel = refLevel + 1; ilevel < static_cast<int>(this->GetNumberOfLevels());
 296        ilevel++ )
 297     {
 298
 299     requestedRegion = this->GetOutput( ilevel - 1 )->GetRequestedRegion();
 300     requestedSize = requestedRegion.GetSize();
 301     requestedIndex = requestedRegion.GetIndex();
 302
 303     for( idim = 0; idim < static_cast<int>(ImageDimension); idim++ )
 304       {
 305       factors[idim] = this->GetSchedule()[ilevel-1][idim] / this->GetSchedule()[ilevel][idim];
 306
 307       // take into account shrink component
 308       requestedSize[idim] *= static_cast<SizeValueType>(factors[idim]);
 309       requestedIndex[idim] *= static_cast<IndexValueType>(factors[idim]);
 310
 311       // take into account smoothing component
 312       if( factors[idim] > 1 )
 313         {
 314         oper->SetDirection( idim );
 315         oper->SetVariance( std::sqr( 0.5 *
 316                                          static_cast<float>( factors[idim] ) ) );
 317         oper->CreateDirectional();
 318         radius[idim] = oper->GetRadius()[idim];
 319         }
 320       else
 321         {
 322         radius[idim] = 0;
 323         }
 324       }
 325
 326     requestedRegion.SetSize( requestedSize );
 327     requestedRegion.SetIndex( requestedIndex );
 328     requestedRegion.PadByRadius( radius );
 329     requestedRegion.Crop( this->GetOutput(ilevel)->
 330                           GetLargestPossibleRegion() );
 331
 332     this->GetOutput(ilevel)->SetRequestedRegion( requestedRegion );
 333
 334     }
 335
 336
 337   // compute requested regions for higher levels
 338   for( ilevel = refLevel - 1; ilevel > -1; ilevel-- )
 339     {
 340     requestedRegion = this->GetOutput( ilevel + 1 )->GetRequestedRegion();
 341     requestedSize = requestedRegion.GetSize();
 342     requestedIndex = requestedRegion.GetIndex();
 343
 344     for( idim = 0; idim < static_cast<int>(ImageDimension); idim++ )
 345       {
 346
 347       factors[idim] = this->GetSchedule()[ilevel][idim] / this->GetSchedule()[ilevel+1][idim];
 348
 349       // take into account smoothing component
 350       if( factors[idim] > 1 )
 351         {
 352         oper->SetDirection( idim );
 353         oper->SetVariance( std::sqr( 0.5 *
 354                                          static_cast<float>( factors[idim] ) ) );
 355         oper->CreateDirectional();
 356         radius[idim] = oper->GetRadius()[idim];
 357         }
 358       else
 359         {
 360         radius[idim] = 0;
 361         }
 362
 363       requestedSize[idim] -= static_cast<SizeValueType>(
 364         2 * radius[idim] );
 365       requestedIndex[idim] += radius[idim];
 366
 367       // take into account shrink component
 368       requestedSize[idim] = static_cast<SizeValueType>( std::floor(
 369                                                           static_cast<double>(requestedSize[idim]) /
 370                                                           static_cast<double>(factors[idim]) ) );
 371       if( requestedSize[idim] < 1 )
 372         {
 373         requestedSize[idim] = 1;
 374         }
 375       requestedIndex[idim] = static_cast<IndexValueType>( std::ceil(
 376                                                             static_cast<double>(requestedIndex[idim]) /
 377                                                             static_cast<double>(factors[idim]) ) );
 378
 379       }
 380
 381     requestedRegion.SetSize( requestedSize );
 382     requestedRegion.SetIndex( requestedIndex );
 383     requestedRegion.Crop( this->GetOutput(ilevel)->
 384                           GetLargestPossibleRegion() );
 385
 386     this->GetOutput(ilevel)->SetRequestedRegion( requestedRegion );
 387
 388     }
 389
 390   // clean up
 391   delete oper;
 392
 393 }
 394
 395 /**
 396  * GenerateInputRequestedRegion
 397  */
 398 template <class TInputImage, class TOutputImage>
 399 void
 400 RecursiveSpatioTemporalMultiResolutionPyramidImageFilter<TInputImage, TOutputImage>
 401 ::GenerateInputRequestedRegion()
 402 {
 403
 404   // call the superclass' implementation of this method
 405   Superclass::GenerateInputRequestedRegion();
 406
 407   // get pointers to the input and output
 408   InputImagePointer  inputPtr =
 409     const_cast< InputImageType *>( this->GetInput() );
 410   if ( !inputPtr )
 411     {
 412     itkExceptionMacro( << "Input has not been set." );
 413     }
 414
 415   // compute baseIndex and baseSize
 416   typedef typename OutputImageType::SizeType    SizeType;
 417   typedef typename SizeType::SizeValueType      SizeValueType;
 418   typedef typename OutputImageType::IndexType   IndexType;
 419   typedef typename IndexType::IndexValueType    IndexValueType;
 420   typedef typename OutputImageType::RegionType  RegionType;
 421
 422   unsigned int refLevel = this->GetNumberOfLevels() - 1;
 423   SizeType baseSize = this->GetOutput(refLevel)->GetRequestedRegion().GetSize();
 424   IndexType baseIndex = this->GetOutput(refLevel)->GetRequestedRegion().GetIndex();
 425   RegionType baseRegion;
 426
 427   unsigned int idim;
 428   for( idim = 0; idim < ImageDimension; idim++ )
 429     {
 430     unsigned int factor = this->GetSchedule()[refLevel][idim];
 431     baseIndex[idim] *= static_cast<IndexValueType>( factor );
 432     baseSize[idim] *= static_cast<SizeValueType>( factor );
 433     }
 434   baseRegion.SetIndex( baseIndex );
 435   baseRegion.SetSize( baseSize );
 436
 437   // compute requirements for the smoothing part
 438   typedef typename TOutputImage::PixelType                 OutputPixelType;
 439   typedef itk::GaussianOperator<OutputPixelType,ImageDimension> OperatorType;
 440
 441   OperatorType *oper = new OperatorType;
 442
 443   typename TInputImage::SizeType radius;
 444
 445   RegionType inputRequestedRegion = baseRegion;
 446   refLevel = 0;
 447
 448   for( idim = 0; idim < TInputImage::ImageDimension; idim++ )
 449     {
 450     oper->SetDirection(idim);
 451     oper->SetVariance( std::sqr( 0.5 * static_cast<float>(
 452                                        this->GetSchedule()[refLevel][idim] ) ) );
 453     oper->SetMaximumError( this->GetMaximumError() );
 454     oper->CreateDirectional();
 455     radius[idim] = oper->GetRadius()[idim];
 456     if( this->GetSchedule()[refLevel][idim] <= 1 )
 457       {
 458       radius[idim] = 0;
 459       }
 460
 461     }
 462   delete oper;
 463
 464   inputRequestedRegion.PadByRadius( radius );
 465
 466   // make sure the requested region is within the largest possible
 467   inputRequestedRegion.Crop( inputPtr->GetLargestPossibleRegion() );
 468
 469   // set the input requested region
 470   inputPtr->SetRequestedRegion( inputRequestedRegion );
 471
 472 }
 473
 474
 475 } // namespace clitk
 476
 477 #endif