registration/clitkSpatioTemporalMultiResolutionPyramidImageFilter.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 __clitkSpatioTemporalMultiResolutionPyramidImageFilter_txx
  19 #define __clitkSpatioTemporalMultiResolutionPyramidImageFilter_txx
  20 #include "clitkSpatioTemporalMultiResolutionPyramidImageFilter.h"
  21 #include "itkGaussianOperator.h"
  22 #include "itkCastImageFilter.h"
  23 #include "itkDiscreteGaussianImageFilter.h"
  24 #include "itkExceptionObject.h"
  25 #include "itkResampleImageFilter.h"
  26 #include "itkShrinkImageFilter.h"
  27 #include "itkIdentityTransform.h"
  28
  29 namespace clitk
  30 {
  31
  32 /**
  33  * Constructor
  34  */
  35 template <class TInputImage, class TOutputImage>
  36 SpatioTemporalMultiResolutionPyramidImageFilter<TInputImage, TOutputImage>
  37 ::SpatioTemporalMultiResolutionPyramidImageFilter()
  38 {
  39   m_NumberOfLevels = 0;
  40   this->SetNumberOfLevels( 2 );
  41   m_MaximumError = 0.1;
  42   m_UseShrinkImageFilter = false;
  43 }
  44
  45
  46 /**
  47  * Set the number of computation levels
  48  */
  49 template <class TInputImage, class TOutputImage>
  50 void
  51 SpatioTemporalMultiResolutionPyramidImageFilter<TInputImage, TOutputImage>
  52 ::SetNumberOfLevels(
  53   unsigned int num )
  54 {
  55   if( m_NumberOfLevels == num )
  56     {
  57     return;
  58     }
  59
  60   this->Modified();
  61
  62   // clamp value to be at least one
  63   m_NumberOfLevels = num;
  64   if( m_NumberOfLevels < 1 ) m_NumberOfLevels = 1;
  65
  66   // resize the schedules
  67   ScheduleType temp( m_NumberOfLevels, ImageDimension );
  68   temp.Fill( 0 );
  69   m_Schedule = temp;
  70
  71   // determine initial shrink factor
  72   unsigned int startfactor = 1;
  73   startfactor = startfactor << ( m_NumberOfLevels - 1 );
  74
  75   // set the starting shrink factors
  76   this->SetStartingShrinkFactors( startfactor );
  77
  78   // set the required number of outputs
  79   this->SetNumberOfRequiredOutputs( m_NumberOfLevels );
  80
  81   unsigned int numOutputs = static_cast<unsigned int>( this->GetNumberOfOutputs() );
  82   unsigned int idx;
  83   if( numOutputs < m_NumberOfLevels )
  84     {
  85     // add extra outputs
  86     for( idx = numOutputs; idx < m_NumberOfLevels; idx++ )
  87       {
  88       typename itk::DataObject::Pointer output =
  89         this->MakeOutput( idx );
  90       this->SetNthOutput( idx, output.GetPointer() );
  91       }
  92
  93     }
  94   else if( numOutputs > m_NumberOfLevels )
  95     {
  96     // remove extra outputs
  97     for( idx = m_NumberOfLevels; idx < numOutputs; idx++ )
  98       {
  99       typename itk::DataObject::Pointer output =
 100         this->GetOutputs()[idx];
 101       this->RemoveOutput( output );
 102       }
 103     }
 104
 105 }
 106
 107
 108 /*
 109  * Set the starting shrink factors
 110  */
 111 template <class TInputImage, class TOutputImage>
 112 void
 113 SpatioTemporalMultiResolutionPyramidImageFilter<TInputImage, TOutputImage>
 114 ::SetStartingShrinkFactors(
 115   unsigned int factor )
 116 {
 117
 118   unsigned int array[ImageDimension];
 119   //JV temporal dimension always 1
 120   for( unsigned int dim = 0; dim < ImageDimension-1; ++dim )
 121     {
 122     array[dim] = factor;
 123     }
 124   array[ImageDimension-1]=1;
 125
 126   this->SetStartingShrinkFactors( array );
 127
 128 }
 129
 130
 131 /**
 132  * Set the starting shrink factors
 133  */
 134 template <class TInputImage, class TOutputImage>
 135 void
 136 SpatioTemporalMultiResolutionPyramidImageFilter<TInputImage, TOutputImage>
 137 ::SetStartingShrinkFactors(
 138   unsigned int * factors )
 139 {
 140
 141   for( unsigned int dim = 0; dim < ImageDimension-1; ++dim )
 142     {
 143     m_Schedule[0][dim] = factors[dim];
 144     if( m_Schedule[0][dim] == 0 )
 145       {
 146       m_Schedule[0][dim] = 1;
 147       }
 148     }
 149   //JV temporal dimension always 1
 150   m_Schedule[0][ImageDimension-1]=1;
 151
 152   for( unsigned int level = 1; level < m_NumberOfLevels; ++level )
 153     {
 154       //JV temporal dimension always 1
 155       for( unsigned int dim = 0; dim < ImageDimension-1; ++dim )
 156            {
 157           m_Schedule[level][dim] = m_Schedule[level-1][dim] / 2;
 158           if( m_Schedule[level][dim] == 0 )
 159             {
 160               m_Schedule[level][dim] = 1;
 161             }
 162         }
 163       m_Schedule[level][ImageDimension-1]=1;
 164     }
 165
 166   this->Modified();
 167
 168 }
 169
 170
 171 /*
 172  * Get the starting shrink factors
 173  */
 174 template <class TInputImage, class TOutputImage>
 175 const unsigned int *
 176 SpatioTemporalMultiResolutionPyramidImageFilter<TInputImage, TOutputImage>
 177 ::GetStartingShrinkFactors() const
 178 {
 179   return ( m_Schedule.data_block() );
 180 }
 181
 182
 183 /*
 184  * Set the multi-resolution schedule
 185  */
 186 template <class TInputImage, class TOutputImage>
 187 void
 188 SpatioTemporalMultiResolutionPyramidImageFilter<TInputImage, TOutputImage>
 189 ::SetSchedule(
 190   const ScheduleType& schedule )
 191 {
 192
 193   if( schedule.rows() != m_NumberOfLevels ||
 194       schedule.columns() != ImageDimension )
 195     {
 196     itkDebugMacro(<< "Schedule has wrong dimensions" );
 197     return;
 198     }
 199
 200   if( schedule == m_Schedule )
 201     {
 202     return;
 203     }
 204
 205   this->Modified();
 206   unsigned int level, dim;
 207   for( level = 0; level < m_NumberOfLevels; level++ )
 208     {
 209       //JV temporal dimension always 1
 210       for( dim = 0; dim < ImageDimension-1; dim++ )
 211         {
 212
 213           m_Schedule[level][dim] = schedule[level][dim];
 214
 215           // set schedule to max( 1, min(schedule[level],
 216           //  schedule[level-1] );
 217           if( level > 0 )
 218             {
 219               m_Schedule[level][dim] = std::min( m_Schedule[level][dim], m_Schedule[level-1][dim] );
 220             }
 221
 222           if( m_Schedule[level][dim] < 1 )
 223             {
 224               m_Schedule[level][dim] = 1;
 225             }
 226         }
 227       m_Schedule[level][ImageDimension-1]=1;
 228     }
 229 }
 230
 231
 232 /*
 233  * Is the schedule downward divisible ?
 234  */
 235 template <class TInputImage, class TOutputImage>
 236 bool
 237 SpatioTemporalMultiResolutionPyramidImageFilter<TInputImage, TOutputImage>
 238 ::IsScheduleDownwardDivisible( const ScheduleType& schedule )
 239 {
 240
 241   unsigned int ilevel, idim;
 242   for( ilevel = 0; ilevel < schedule.rows() - 1; ilevel++ )
 243     {
 244     for( idim = 0; idim < schedule.columns(); idim++ )
 245       {
 246       if( schedule[ilevel][idim] == 0 )
 247         {
 248         return false;
 249         }
 250       if( ( schedule[ilevel][idim] % schedule[ilevel+1][idim] ) > 0 )
 251         {
 252         return false;
 253         }
 254       }
 255     }
 256
 257   return true;
 258 }
 259
 260 /*
 261  * GenerateData for non downward divisible schedules
 262  */
 263 template <class TInputImage, class TOutputImage>
 264 void
 265 SpatioTemporalMultiResolutionPyramidImageFilter<TInputImage, TOutputImage>
 266 ::GenerateData()
 267 {
 268   // Get the input and output pointers
 269   InputImageConstPointer  inputPtr = this->GetInput();
 270
 271   // Create caster, smoother and resampleShrinker filters
 272   typedef itk::CastImageFilter<TInputImage, TOutputImage>              CasterType;
 273   typedef itk::DiscreteGaussianImageFilter<TOutputImage, TOutputImage> SmootherType;
 274
 275   typedef itk::ImageToImageFilter<TOutputImage,TOutputImage>           ImageToImageType;
 276   typedef itk::ResampleImageFilter<TOutputImage,TOutputImage>          ResampleShrinkerType;
 277   typedef itk::ShrinkImageFilter<TOutputImage,TOutputImage>            ShrinkerType;
 278
 279   typename CasterType::Pointer caster = CasterType::New();
 280   typename SmootherType::Pointer smoother = SmootherType::New();
 281
 282   typename ImageToImageType::Pointer shrinkerFilter;
 283   //
 284   // only one of these pointers is going to be valid, depending on the
 285   // value of UseShrinkImageFilter flag
 286   typename ResampleShrinkerType::Pointer resampleShrinker;
 287   typename ShrinkerType::Pointer shrinker;
 288
 289   if(this->GetUseShrinkImageFilter())
 290     {
 291     shrinker = ShrinkerType::New();
 292     shrinkerFilter = shrinker.GetPointer();
 293     }
 294   else
 295     {
 296     resampleShrinker = ResampleShrinkerType::New();
 297     typedef itk::LinearInterpolateImageFunction< OutputImageType, double >
 298       LinearInterpolatorType;
 299     typename LinearInterpolatorType::Pointer interpolator =
 300       LinearInterpolatorType::New();
 301     resampleShrinker->SetInterpolator( interpolator );
 302     resampleShrinker->SetDefaultPixelValue( 0 );
 303     shrinkerFilter = resampleShrinker.GetPointer();
 304     }
 305   // Setup the filters
 306   caster->SetInput( inputPtr );
 307
 308   smoother->SetUseImageSpacing( false );
 309   smoother->SetInput( caster->GetOutput() );
 310   smoother->SetMaximumError( m_MaximumError );
 311
 312   shrinkerFilter->SetInput( smoother->GetOutput() );
 313
 314   unsigned int ilevel, idim;
 315   unsigned int factors[ImageDimension];
 316   double       variance[ImageDimension];
 317
 318   for( ilevel = 0; ilevel < m_NumberOfLevels; ilevel++ )
 319     {
 320     this->UpdateProgress( static_cast<float>( ilevel ) /
 321                           static_cast<float>( m_NumberOfLevels ) );
 322
 323     // Allocate memory for each output
 324     OutputImagePointer outputPtr = this->GetOutput( ilevel );
 325     outputPtr->SetBufferedRegion( outputPtr->GetRequestedRegion() );
 326     outputPtr->Allocate();
 327
 328     // compute shrink factors and variances
 329     for( idim = 0; idim < ImageDimension; idim++ )
 330       {
 331       factors[idim] = m_Schedule[ilevel][idim];
 332       variance[idim] = std::sqr( 0.5 *
 333                                      static_cast<float>( factors[idim] ) );
 334       }
 335
 336     if(!this->GetUseShrinkImageFilter())
 337       {
 338       typedef itk::IdentityTransform<double,OutputImageType::ImageDimension>
 339         IdentityTransformType;
 340       typename IdentityTransformType::Pointer identityTransform =
 341         IdentityTransformType::New();
 342       resampleShrinker->SetOutputParametersFromImage( outputPtr );
 343       resampleShrinker->SetTransform(identityTransform);
 344       }
 345     else
 346       {
 347       shrinker->SetShrinkFactors(factors);
 348       }
 349     // use mini-pipeline to compute output
 350     smoother->SetVariance( variance );
 351
 352     shrinkerFilter->GraftOutput( outputPtr );
 353
 354     // force to always update in case shrink factors are the same
 355     shrinkerFilter->Modified();
 356     shrinkerFilter->UpdateLargestPossibleRegion();
 357     this->GraftNthOutput( ilevel, shrinkerFilter->GetOutput() );
 358     }
 359 }
 360
 361 /**
 362  * PrintSelf method
 363  */
 364 template <class TInputImage, class TOutputImage>
 365 void
 366 SpatioTemporalMultiResolutionPyramidImageFilter<TInputImage, TOutputImage>
 367 ::PrintSelf(std::ostream& os, itk::Indent indent) const
 368 {
 369   Superclass::PrintSelf(os,indent);
 370
 371   os << indent << "MaximumError: " << m_MaximumError << std::endl;
 372   os << indent << "No. levels: " << m_NumberOfLevels << std::endl;
 373   os << indent << "Schedule: " << std::endl;
 374   os << m_Schedule << std::endl;
 375   os << "Use ShrinkImageFilter= " << m_UseShrinkImageFilter << std::endl;
 376 }
 377
 378
 379 /*
 380  * GenerateOutputInformation
 381  */
 382 template <class TInputImage, class TOutputImage>
 383 void
 384 SpatioTemporalMultiResolutionPyramidImageFilter<TInputImage, TOutputImage>
 385 ::GenerateOutputInformation()
 386 {
 387
 388   // call the superclass's implementation of this method
 389   Superclass::GenerateOutputInformation();
 390
 391   // get pointers to the input and output
 392   InputImageConstPointer inputPtr = this->GetInput();
 393
 394   if ( !inputPtr  )
 395     {
 396     itkExceptionMacro( << "Input has not been set" );
 397     }
 398
 399   const typename InputImageType::PointType&
 400     inputOrigin = inputPtr->GetOrigin();
 401   const typename InputImageType::SpacingType&
 402     inputSpacing = inputPtr->GetSpacing();
 403   const typename InputImageType::DirectionType&
 404     inputDirection = inputPtr->GetDirection();
 405   const typename InputImageType::SizeType& inputSize =
 406     inputPtr->GetLargestPossibleRegion().GetSize();
 407   const typename InputImageType::IndexType& inputStartIndex =
 408     inputPtr->GetLargestPossibleRegion().GetIndex();
 409
 410   typedef typename OutputImageType::SizeType  SizeType;
 411   typedef typename SizeType::SizeValueType    SizeValueType;
 412   typedef typename OutputImageType::IndexType IndexType;
 413   typedef typename IndexType::IndexValueType  IndexValueType;
 414
 415   OutputImagePointer outputPtr;
 416   typename OutputImageType::PointType   outputOrigin;
 417   typename OutputImageType::SpacingType outputSpacing;
 418   SizeType    outputSize;
 419   IndexType   outputStartIndex;
 420
 421   // we need to compute the output spacing, the output image size,
 422   // and the output image start index
 423   for(unsigned int ilevel = 0; ilevel < m_NumberOfLevels; ilevel++ )
 424     {
 425     outputPtr = this->GetOutput( ilevel );
 426     if( !outputPtr ) { continue; }
 427
 428     for(unsigned int idim = 0; idim < OutputImageType::ImageDimension; idim++ )
 429       {
 430       const double shrinkFactor = static_cast<double>( m_Schedule[ilevel][idim] );
 431       outputSpacing[idim] = inputSpacing[idim] * shrinkFactor;
 432
 433       outputSize[idim] = static_cast<SizeValueType>(
 434         std::floor(static_cast<double>(inputSize[idim]) / shrinkFactor ) );
 435       if( outputSize[idim] < 1 ) { outputSize[idim] = 1; }
 436
 437       outputStartIndex[idim] = static_cast<IndexValueType>(
 438         std::ceil(static_cast<double>(inputStartIndex[idim]) / shrinkFactor ) );
 439       }
 440     //Now compute the new shifted origin for the updated levels;
 441     const typename OutputImageType::PointType::VectorType outputOriginOffset
 442          =(inputDirection*(outputSpacing-inputSpacing))*0.5;
 443     for(unsigned int idim = 0; idim < OutputImageType::ImageDimension; idim++ )
 444       {
 445         outputOrigin[idim]=inputOrigin[idim]+outputOriginOffset[idim];
 446       }
 447
 448     typename OutputImageType::RegionType outputLargestPossibleRegion;
 449     outputLargestPossibleRegion.SetSize( outputSize );
 450     outputLargestPossibleRegion.SetIndex( outputStartIndex );
 451
 452     outputPtr->SetLargestPossibleRegion( outputLargestPossibleRegion );
 453     outputPtr->SetOrigin ( outputOrigin );
 454     outputPtr->SetSpacing( outputSpacing );
 455     outputPtr->SetDirection( inputDirection );//Output Direction should be same as input.
 456     }
 457 }
 458
 459
 460 /*
 461  * GenerateOutputRequestedRegion
 462  */
 463 template <class TInputImage, class TOutputImage>
 464 void
 465 SpatioTemporalMultiResolutionPyramidImageFilter<TInputImage, TOutputImage>
 466 ::GenerateOutputRequestedRegion(itk::DataObject * refOutput )
 467 {
 468   // call the superclass's implementation of this method
 469   Superclass::GenerateOutputRequestedRegion( refOutput );
 470
 471   // find the index for this output
 472   unsigned int refLevel = refOutput->GetSourceOutputIndex();
 473
 474   // compute baseIndex and baseSize
 475   typedef typename OutputImageType::SizeType    SizeType;
 476   typedef typename SizeType::SizeValueType      SizeValueType;
 477   typedef typename OutputImageType::IndexType   IndexType;
 478   typedef typename IndexType::IndexValueType    IndexValueType;
 479   typedef typename OutputImageType::RegionType  RegionType;
 480
 481   TOutputImage * ptr = static_cast<TOutputImage*>( refOutput );
 482   if( !ptr )
 483     {
 484     itkExceptionMacro( << "Could not cast refOutput to TOutputImage*." );
 485     }
 486
 487   unsigned int ilevel, idim;
 488
 489   if ( ptr->GetRequestedRegion() == ptr->GetLargestPossibleRegion() )
 490     {
 491     // set the requested regions for the other outputs to their
 492     // requested region
 493
 494     for( ilevel = 0; ilevel < m_NumberOfLevels; ilevel++ )
 495       {
 496       if( ilevel == refLevel ) { continue; }
 497       if( !this->GetOutput(ilevel) ) { continue; }
 498       this->GetOutput(ilevel)->SetRequestedRegionToLargestPossibleRegion();
 499       }
 500     }
 501   else
 502     {
 503     // compute requested regions for the other outputs based on
 504     // the requested region of the reference output
 505     IndexType outputIndex;
 506     SizeType  outputSize;
 507     RegionType outputRegion;
 508     IndexType  baseIndex = ptr->GetRequestedRegion().GetIndex();
 509     SizeType   baseSize  = ptr->GetRequestedRegion().GetSize();
 510
 511     for( idim = 0; idim < TOutputImage::ImageDimension; idim++ )
 512       {
 513       unsigned int factor = m_Schedule[refLevel][idim];
 514       baseIndex[idim] *= static_cast<IndexValueType>( factor );
 515       baseSize[idim] *= static_cast<SizeValueType>( factor );
 516       }
 517
 518     for( ilevel = 0; ilevel < m_NumberOfLevels; ilevel++ )
 519       {
 520       if( ilevel == refLevel ) { continue; }
 521       if( !this->GetOutput(ilevel) ) { continue; }
 522
 523       for( idim = 0; idim < TOutputImage::ImageDimension; idim++ )
 524         {
 525
 526         double factor = static_cast<double>( m_Schedule[ilevel][idim] );
 527
 528         outputSize[idim] = static_cast<SizeValueType>(
 529           std::floor(static_cast<double>(baseSize[idim]) / factor ) );
 530         if( outputSize[idim] < 1 ) { outputSize[idim] = 1; }
 531
 532         outputIndex[idim] = static_cast<IndexValueType>(
 533           std::ceil(static_cast<double>(baseIndex[idim]) / factor ) );
 534
 535         }
 536
 537       outputRegion.SetIndex( outputIndex );
 538       outputRegion.SetSize( outputSize );
 539
 540       // make sure the region is within the largest possible region
 541       outputRegion.Crop( this->GetOutput( ilevel )->
 542                          GetLargestPossibleRegion() );
 543       // set the requested region
 544       this->GetOutput( ilevel )->SetRequestedRegion( outputRegion );
 545       }
 546
 547     }
 548 }
 549
 550
 551 /**
 552  * GenerateInputRequestedRegion
 553  */
 554 template <class TInputImage, class TOutputImage>
 555 void
 556 SpatioTemporalMultiResolutionPyramidImageFilter<TInputImage, TOutputImage>
 557 ::GenerateInputRequestedRegion()
 558 {
 559   // call the superclass' implementation of this method
 560   Superclass::GenerateInputRequestedRegion();
 561
 562   // get pointers to the input and output
 563   InputImagePointer  inputPtr =
 564     const_cast< InputImageType * >( this->GetInput() );
 565   if ( !inputPtr )
 566     {
 567     itkExceptionMacro( << "Input has not been set." );
 568     }
 569
 570   // compute baseIndex and baseSize
 571   typedef typename OutputImageType::SizeType    SizeType;
 572   typedef typename SizeType::SizeValueType      SizeValueType;
 573   typedef typename OutputImageType::IndexType   IndexType;
 574   typedef typename IndexType::IndexValueType    IndexValueType;
 575   typedef typename OutputImageType::RegionType  RegionType;
 576
 577   unsigned int refLevel = m_NumberOfLevels - 1;
 578   SizeType baseSize = this->GetOutput(refLevel)->GetRequestedRegion().GetSize();
 579   IndexType baseIndex = this->GetOutput(refLevel)->GetRequestedRegion().GetIndex();
 580   RegionType baseRegion;
 581
 582   unsigned int idim;
 583   for( idim = 0; idim < ImageDimension; idim++ )
 584     {
 585     unsigned int factor = m_Schedule[refLevel][idim];
 586     baseIndex[idim] *= static_cast<IndexValueType>( factor );
 587     baseSize[idim] *= static_cast<SizeValueType>( factor );
 588     }
 589   baseRegion.SetIndex( baseIndex );
 590   baseRegion.SetSize( baseSize );
 591
 592   // compute requirements for the smoothing part
 593   typedef typename TOutputImage::PixelType                 OutputPixelType;
 594   typedef typename itk::GaussianOperator<OutputPixelType,ImageDimension> OperatorType;
 595
 596   OperatorType *oper = new OperatorType;
 597
 598   typename TInputImage::SizeType radius;
 599
 600   RegionType inputRequestedRegion = baseRegion;
 601   refLevel = 0;
 602
 603   for( idim = 0; idim < TInputImage::ImageDimension; idim++ )
 604     {
 605     oper->SetDirection(idim);
 606     oper->SetVariance( std::sqr( 0.5 * static_cast<float>(
 607                                        m_Schedule[refLevel][idim] ) ) );
 608     oper->SetMaximumError( m_MaximumError );
 609     oper->CreateDirectional();
 610     radius[idim] = oper->GetRadius()[idim];
 611     }
 612   delete oper;
 613
 614   inputRequestedRegion.PadByRadius( radius );
 615
 616   // make sure the requested region is within the largest possible
 617   inputRequestedRegion.Crop( inputPtr->GetLargestPossibleRegion() );
 618
 619   // set the input requested region
 620   inputPtr->SetRequestedRegion( inputRequestedRegion );
 621
 622 }
 623
 624
 625 } // namespace clitk
 626
 627 #endif