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