1 /*=========================================================================
2 Program: vv http://www.creatis.insa-lyon.fr/rio/vv
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
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.
13 It is distributed under dual licence
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"
35 template <class TInputImage, class TOutputImage>
36 SpatioTemporalMultiResolutionPyramidImageFilter<TInputImage, TOutputImage>
37 ::SpatioTemporalMultiResolutionPyramidImageFilter()
40 this->SetNumberOfLevels( 2 );
42 m_UseShrinkImageFilter = false;
47 * Set the number of computation levels
49 template <class TInputImage, class TOutputImage>
51 SpatioTemporalMultiResolutionPyramidImageFilter<TInputImage, TOutputImage>
55 if( m_NumberOfLevels == num )
62 // clamp value to be at least one
63 m_NumberOfLevels = num;
64 if( m_NumberOfLevels < 1 ) m_NumberOfLevels = 1;
66 // resize the schedules
67 ScheduleType temp( m_NumberOfLevels, ImageDimension );
71 // determine initial shrink factor
72 unsigned int startfactor = 1;
73 startfactor = startfactor << ( m_NumberOfLevels - 1 );
75 // set the starting shrink factors
76 this->SetStartingShrinkFactors( startfactor );
78 // set the required number of outputs
79 this->SetNumberOfRequiredOutputs( m_NumberOfLevels );
81 unsigned int numOutputs = static_cast<unsigned int>( this->GetNumberOfOutputs() );
83 if( numOutputs < m_NumberOfLevels )
86 for( idx = numOutputs; idx < m_NumberOfLevels; idx++ )
88 typename itk::DataObject::Pointer output =
89 this->MakeOutput( idx );
90 this->SetNthOutput( idx, output.GetPointer() );
94 else if( numOutputs > m_NumberOfLevels )
96 // remove extra outputs
97 for( idx = m_NumberOfLevels; idx < numOutputs; idx++ )
99 typename itk::DataObject::Pointer output =
100 this->GetOutputs()[idx];
101 this->RemoveOutput( output );
109 * Set the starting shrink factors
111 template <class TInputImage, class TOutputImage>
113 SpatioTemporalMultiResolutionPyramidImageFilter<TInputImage, TOutputImage>
114 ::SetStartingShrinkFactors(
115 unsigned int factor )
118 unsigned int array[ImageDimension];
119 //JV temporal dimension always 1
120 for( unsigned int dim = 0; dim < ImageDimension-1; ++dim )
124 array[ImageDimension-1]=1;
126 this->SetStartingShrinkFactors( array );
132 * Set the starting shrink factors
134 template <class TInputImage, class TOutputImage>
136 SpatioTemporalMultiResolutionPyramidImageFilter<TInputImage, TOutputImage>
137 ::SetStartingShrinkFactors(
138 unsigned int * factors )
141 for( unsigned int dim = 0; dim < ImageDimension-1; ++dim )
143 m_Schedule[0][dim] = factors[dim];
144 if( m_Schedule[0][dim] == 0 )
146 m_Schedule[0][dim] = 1;
149 //JV temporal dimension always 1
150 m_Schedule[0][ImageDimension-1]=1;
152 for( unsigned int level = 1; level < m_NumberOfLevels; ++level )
154 //JV temporal dimension always 1
155 for( unsigned int dim = 0; dim < ImageDimension-1; ++dim )
157 m_Schedule[level][dim] = m_Schedule[level-1][dim] / 2;
158 if( m_Schedule[level][dim] == 0 )
160 m_Schedule[level][dim] = 1;
163 m_Schedule[level][ImageDimension-1]=1;
172 * Get the starting shrink factors
174 template <class TInputImage, class TOutputImage>
176 SpatioTemporalMultiResolutionPyramidImageFilter<TInputImage, TOutputImage>
177 ::GetStartingShrinkFactors() const
179 return ( m_Schedule.data_block() );
184 * Set the multi-resolution schedule
186 template <class TInputImage, class TOutputImage>
188 SpatioTemporalMultiResolutionPyramidImageFilter<TInputImage, TOutputImage>
190 const ScheduleType& schedule )
193 if( schedule.rows() != m_NumberOfLevels ||
194 schedule.columns() != ImageDimension )
196 itkDebugMacro(<< "Schedule has wrong dimensions" );
200 if( schedule == m_Schedule )
206 unsigned int level, dim;
207 for( level = 0; level < m_NumberOfLevels; level++ )
209 //JV temporal dimension always 1
210 for( dim = 0; dim < ImageDimension-1; dim++ )
213 m_Schedule[level][dim] = schedule[level][dim];
215 // set schedule to max( 1, min(schedule[level],
216 // schedule[level-1] );
219 m_Schedule[level][dim] = std::min( m_Schedule[level][dim], m_Schedule[level-1][dim] );
222 if( m_Schedule[level][dim] < 1 )
224 m_Schedule[level][dim] = 1;
227 m_Schedule[level][ImageDimension-1]=1;
233 * Is the schedule downward divisible ?
235 template <class TInputImage, class TOutputImage>
237 SpatioTemporalMultiResolutionPyramidImageFilter<TInputImage, TOutputImage>
238 ::IsScheduleDownwardDivisible( const ScheduleType& schedule )
241 unsigned int ilevel, idim;
242 for( ilevel = 0; ilevel < schedule.rows() - 1; ilevel++ )
244 for( idim = 0; idim < schedule.columns(); idim++ )
246 if( schedule[ilevel][idim] == 0 )
250 if( ( schedule[ilevel][idim] % schedule[ilevel+1][idim] ) > 0 )
261 * GenerateData for non downward divisible schedules
263 template <class TInputImage, class TOutputImage>
265 SpatioTemporalMultiResolutionPyramidImageFilter<TInputImage, TOutputImage>
268 // Get the input and output pointers
269 InputImageConstPointer inputPtr = this->GetInput();
271 // Create caster, smoother and resampleShrinker filters
272 typedef itk::CastImageFilter<TInputImage, TOutputImage> CasterType;
273 typedef itk::DiscreteGaussianImageFilter<TOutputImage, TOutputImage> SmootherType;
275 typedef itk::ImageToImageFilter<TOutputImage,TOutputImage> ImageToImageType;
276 typedef itk::ResampleImageFilter<TOutputImage,TOutputImage> ResampleShrinkerType;
277 typedef itk::ShrinkImageFilter<TOutputImage,TOutputImage> ShrinkerType;
279 typename CasterType::Pointer caster = CasterType::New();
280 typename SmootherType::Pointer smoother = SmootherType::New();
282 typename ImageToImageType::Pointer shrinkerFilter;
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;
289 if(this->GetUseShrinkImageFilter())
291 shrinker = ShrinkerType::New();
292 shrinkerFilter = shrinker.GetPointer();
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();
306 caster->SetInput( inputPtr );
308 smoother->SetUseImageSpacing( false );
309 smoother->SetInput( caster->GetOutput() );
310 smoother->SetMaximumError( m_MaximumError );
312 shrinkerFilter->SetInput( smoother->GetOutput() );
314 unsigned int ilevel, idim;
315 unsigned int factors[ImageDimension];
316 double variance[ImageDimension];
318 for( ilevel = 0; ilevel < m_NumberOfLevels; ilevel++ )
320 this->UpdateProgress( static_cast<float>( ilevel ) /
321 static_cast<float>( m_NumberOfLevels ) );
323 // Allocate memory for each output
324 OutputImagePointer outputPtr = this->GetOutput( ilevel );
325 outputPtr->SetBufferedRegion( outputPtr->GetRequestedRegion() );
326 outputPtr->Allocate();
328 // compute shrink factors and variances
329 for( idim = 0; idim < ImageDimension; idim++ )
331 factors[idim] = m_Schedule[ilevel][idim];
332 variance[idim] = std::sqr( 0.5 *
333 static_cast<float>( factors[idim] ) );
336 if(!this->GetUseShrinkImageFilter())
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);
347 shrinker->SetShrinkFactors(factors);
349 // use mini-pipeline to compute output
350 smoother->SetVariance( variance );
352 shrinkerFilter->GraftOutput( outputPtr );
354 // force to always update in case shrink factors are the same
355 shrinkerFilter->Modified();
356 shrinkerFilter->UpdateLargestPossibleRegion();
357 this->GraftNthOutput( ilevel, shrinkerFilter->GetOutput() );
364 template <class TInputImage, class TOutputImage>
366 SpatioTemporalMultiResolutionPyramidImageFilter<TInputImage, TOutputImage>
367 ::PrintSelf(std::ostream& os, itk::Indent indent) const
369 Superclass::PrintSelf(os,indent);
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;
380 * GenerateOutputInformation
382 template <class TInputImage, class TOutputImage>
384 SpatioTemporalMultiResolutionPyramidImageFilter<TInputImage, TOutputImage>
385 ::GenerateOutputInformation()
388 // call the superclass's implementation of this method
389 Superclass::GenerateOutputInformation();
391 // get pointers to the input and output
392 InputImageConstPointer inputPtr = this->GetInput();
396 itkExceptionMacro( << "Input has not been set" );
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();
410 typedef typename OutputImageType::SizeType SizeType;
411 typedef typename SizeType::SizeValueType SizeValueType;
412 typedef typename OutputImageType::IndexType IndexType;
413 typedef typename IndexType::IndexValueType IndexValueType;
415 OutputImagePointer outputPtr;
416 typename OutputImageType::PointType outputOrigin;
417 typename OutputImageType::SpacingType outputSpacing;
419 IndexType outputStartIndex;
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++ )
425 outputPtr = this->GetOutput( ilevel );
426 if( !outputPtr ) { continue; }
428 for(unsigned int idim = 0; idim < OutputImageType::ImageDimension; idim++ )
430 const double shrinkFactor = static_cast<double>( m_Schedule[ilevel][idim] );
431 outputSpacing[idim] = inputSpacing[idim] * shrinkFactor;
433 outputSize[idim] = static_cast<SizeValueType>(
434 std::floor(static_cast<double>(inputSize[idim]) / shrinkFactor ) );
435 if( outputSize[idim] < 1 ) { outputSize[idim] = 1; }
437 outputStartIndex[idim] = static_cast<IndexValueType>(
438 std::ceil(static_cast<double>(inputStartIndex[idim]) / shrinkFactor ) );
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++ )
445 outputOrigin[idim]=inputOrigin[idim]+outputOriginOffset[idim];
448 typename OutputImageType::RegionType outputLargestPossibleRegion;
449 outputLargestPossibleRegion.SetSize( outputSize );
450 outputLargestPossibleRegion.SetIndex( outputStartIndex );
452 outputPtr->SetLargestPossibleRegion( outputLargestPossibleRegion );
453 outputPtr->SetOrigin ( outputOrigin );
454 outputPtr->SetSpacing( outputSpacing );
455 outputPtr->SetDirection( inputDirection );//Output Direction should be same as input.
461 * GenerateOutputRequestedRegion
463 template <class TInputImage, class TOutputImage>
465 SpatioTemporalMultiResolutionPyramidImageFilter<TInputImage, TOutputImage>
466 ::GenerateOutputRequestedRegion(itk::DataObject * refOutput )
468 // call the superclass's implementation of this method
469 Superclass::GenerateOutputRequestedRegion( refOutput );
471 // find the index for this output
472 unsigned int refLevel = refOutput->GetSourceOutputIndex();
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;
481 TOutputImage * ptr = static_cast<TOutputImage*>( refOutput );
484 itkExceptionMacro( << "Could not cast refOutput to TOutputImage*." );
487 unsigned int ilevel, idim;
489 if ( ptr->GetRequestedRegion() == ptr->GetLargestPossibleRegion() )
491 // set the requested regions for the other outputs to their
494 for( ilevel = 0; ilevel < m_NumberOfLevels; ilevel++ )
496 if( ilevel == refLevel ) { continue; }
497 if( !this->GetOutput(ilevel) ) { continue; }
498 this->GetOutput(ilevel)->SetRequestedRegionToLargestPossibleRegion();
503 // compute requested regions for the other outputs based on
504 // the requested region of the reference output
505 IndexType outputIndex;
507 RegionType outputRegion;
508 IndexType baseIndex = ptr->GetRequestedRegion().GetIndex();
509 SizeType baseSize = ptr->GetRequestedRegion().GetSize();
511 for( idim = 0; idim < TOutputImage::ImageDimension; idim++ )
513 unsigned int factor = m_Schedule[refLevel][idim];
514 baseIndex[idim] *= static_cast<IndexValueType>( factor );
515 baseSize[idim] *= static_cast<SizeValueType>( factor );
518 for( ilevel = 0; ilevel < m_NumberOfLevels; ilevel++ )
520 if( ilevel == refLevel ) { continue; }
521 if( !this->GetOutput(ilevel) ) { continue; }
523 for( idim = 0; idim < TOutputImage::ImageDimension; idim++ )
526 double factor = static_cast<double>( m_Schedule[ilevel][idim] );
528 outputSize[idim] = static_cast<SizeValueType>(
529 std::floor(static_cast<double>(baseSize[idim]) / factor ) );
530 if( outputSize[idim] < 1 ) { outputSize[idim] = 1; }
532 outputIndex[idim] = static_cast<IndexValueType>(
533 std::ceil(static_cast<double>(baseIndex[idim]) / factor ) );
537 outputRegion.SetIndex( outputIndex );
538 outputRegion.SetSize( outputSize );
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 );
552 * GenerateInputRequestedRegion
554 template <class TInputImage, class TOutputImage>
556 SpatioTemporalMultiResolutionPyramidImageFilter<TInputImage, TOutputImage>
557 ::GenerateInputRequestedRegion()
559 // call the superclass' implementation of this method
560 Superclass::GenerateInputRequestedRegion();
562 // get pointers to the input and output
563 InputImagePointer inputPtr =
564 const_cast< InputImageType * >( this->GetInput() );
567 itkExceptionMacro( << "Input has not been set." );
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;
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;
583 for( idim = 0; idim < ImageDimension; idim++ )
585 unsigned int factor = m_Schedule[refLevel][idim];
586 baseIndex[idim] *= static_cast<IndexValueType>( factor );
587 baseSize[idim] *= static_cast<SizeValueType>( factor );
589 baseRegion.SetIndex( baseIndex );
590 baseRegion.SetSize( baseSize );
592 // compute requirements for the smoothing part
593 typedef typename TOutputImage::PixelType OutputPixelType;
594 typedef typename itk::GaussianOperator<OutputPixelType,ImageDimension> OperatorType;
596 OperatorType *oper = new OperatorType;
598 typename TInputImage::SizeType radius;
600 RegionType inputRequestedRegion = baseRegion;
603 for( idim = 0; idim < TInputImage::ImageDimension; idim++ )
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];
614 inputRequestedRegion.PadByRadius( radius );
616 // make sure the requested region is within the largest possible
617 inputRequestedRegion.Crop( inputPtr->GetLargestPossibleRegion() );
619 // set the input requested region
620 inputPtr->SetRequestedRegion( inputRequestedRegion );