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"
29 #include "vnl/vnl_math.h"
37 template <class TInputImage, class TOutputImage>
38 SpatioTemporalMultiResolutionPyramidImageFilter<TInputImage, TOutputImage>
39 ::SpatioTemporalMultiResolutionPyramidImageFilter()
42 this->SetNumberOfLevels( 2 );
44 m_UseShrinkImageFilter = false;
49 * Set the number of computation levels
51 template <class TInputImage, class TOutputImage>
53 SpatioTemporalMultiResolutionPyramidImageFilter<TInputImage, TOutputImage>
57 if( m_NumberOfLevels == num )
64 // clamp value to be at least one
65 m_NumberOfLevels = num;
66 if( m_NumberOfLevels < 1 ) m_NumberOfLevels = 1;
68 // resize the schedules
69 ScheduleType temp( m_NumberOfLevels, ImageDimension );
73 // determine initial shrink factor
74 unsigned int startfactor = 1;
75 startfactor = startfactor << ( m_NumberOfLevels - 1 );
77 // set the starting shrink factors
78 this->SetStartingShrinkFactors( startfactor );
80 // set the required number of outputs
81 this->SetNumberOfRequiredOutputs( m_NumberOfLevels );
83 unsigned int numOutputs = static_cast<unsigned int>( this->GetNumberOfOutputs() );
85 if( numOutputs < m_NumberOfLevels )
88 for( idx = numOutputs; idx < m_NumberOfLevels; idx++ )
90 typename itk::DataObject::Pointer output =
91 this->MakeOutput( idx );
92 this->SetNthOutput( idx, output.GetPointer() );
96 else if( numOutputs > m_NumberOfLevels )
98 // remove extra outputs
99 for( idx = m_NumberOfLevels; idx < numOutputs; idx++ )
101 typename itk::DataObject::Pointer output =
102 this->GetOutputs()[idx];
103 this->RemoveOutput( output );
111 * Set the starting shrink factors
113 template <class TInputImage, class TOutputImage>
115 SpatioTemporalMultiResolutionPyramidImageFilter<TInputImage, TOutputImage>
116 ::SetStartingShrinkFactors(
117 unsigned int factor )
120 unsigned int array[ImageDimension];
121 //JV temporal dimension always 1
122 for( unsigned int dim = 0; dim < ImageDimension-1; ++dim )
126 array[ImageDimension-1]=1;
128 this->SetStartingShrinkFactors( array );
134 * Set the starting shrink factors
136 template <class TInputImage, class TOutputImage>
138 SpatioTemporalMultiResolutionPyramidImageFilter<TInputImage, TOutputImage>
139 ::SetStartingShrinkFactors(
140 unsigned int * factors )
143 for( unsigned int dim = 0; dim < ImageDimension-1; ++dim )
145 m_Schedule[0][dim] = factors[dim];
146 if( m_Schedule[0][dim] == 0 )
148 m_Schedule[0][dim] = 1;
151 //JV temporal dimension always 1
152 m_Schedule[0][ImageDimension-1]=1;
154 for( unsigned int level = 1; level < m_NumberOfLevels; ++level )
156 //JV temporal dimension always 1
157 for( unsigned int dim = 0; dim < ImageDimension-1; ++dim )
159 m_Schedule[level][dim] = m_Schedule[level-1][dim] / 2;
160 if( m_Schedule[level][dim] == 0 )
162 m_Schedule[level][dim] = 1;
165 m_Schedule[level][ImageDimension-1]=1;
174 * Get the starting shrink factors
176 template <class TInputImage, class TOutputImage>
178 SpatioTemporalMultiResolutionPyramidImageFilter<TInputImage, TOutputImage>
179 ::GetStartingShrinkFactors() const
181 return ( m_Schedule.data_block() );
186 * Set the multi-resolution schedule
188 template <class TInputImage, class TOutputImage>
190 SpatioTemporalMultiResolutionPyramidImageFilter<TInputImage, TOutputImage>
192 const ScheduleType& schedule )
195 if( schedule.rows() != m_NumberOfLevels ||
196 schedule.columns() != ImageDimension )
198 itkDebugMacro(<< "Schedule has wrong dimensions" );
202 if( schedule == m_Schedule )
208 unsigned int level, dim;
209 for( level = 0; level < m_NumberOfLevels; level++ )
211 //JV temporal dimension always 1
212 for( dim = 0; dim < ImageDimension-1; dim++ )
215 m_Schedule[level][dim] = schedule[level][dim];
217 // set schedule to max( 1, min(schedule[level],
218 // schedule[level-1] );
221 m_Schedule[level][dim] = vnl_math_min( m_Schedule[level][dim], m_Schedule[level-1][dim] );
224 if( m_Schedule[level][dim] < 1 )
226 m_Schedule[level][dim] = 1;
229 m_Schedule[level][ImageDimension-1]=1;
235 * Is the schedule downward divisible ?
237 template <class TInputImage, class TOutputImage>
239 SpatioTemporalMultiResolutionPyramidImageFilter<TInputImage, TOutputImage>
240 ::IsScheduleDownwardDivisible( const ScheduleType& schedule )
243 unsigned int ilevel, idim;
244 for( ilevel = 0; ilevel < schedule.rows() - 1; ilevel++ )
246 for( idim = 0; idim < schedule.columns(); idim++ )
248 if( schedule[ilevel][idim] == 0 )
252 if( ( schedule[ilevel][idim] % schedule[ilevel+1][idim] ) > 0 )
263 * GenerateData for non downward divisible schedules
265 template <class TInputImage, class TOutputImage>
267 SpatioTemporalMultiResolutionPyramidImageFilter<TInputImage, TOutputImage>
270 // Get the input and output pointers
271 InputImageConstPointer inputPtr = this->GetInput();
273 // Create caster, smoother and resampleShrinker filters
274 typedef itk::CastImageFilter<TInputImage, TOutputImage> CasterType;
275 typedef itk::DiscreteGaussianImageFilter<TOutputImage, TOutputImage> SmootherType;
277 typedef itk::ImageToImageFilter<TOutputImage,TOutputImage> ImageToImageType;
278 typedef itk::ResampleImageFilter<TOutputImage,TOutputImage> ResampleShrinkerType;
279 typedef itk::ShrinkImageFilter<TOutputImage,TOutputImage> ShrinkerType;
281 typename CasterType::Pointer caster = CasterType::New();
282 typename SmootherType::Pointer smoother = SmootherType::New();
284 typename ImageToImageType::Pointer shrinkerFilter;
286 // only one of these pointers is going to be valid, depending on the
287 // value of UseShrinkImageFilter flag
288 typename ResampleShrinkerType::Pointer resampleShrinker;
289 typename ShrinkerType::Pointer shrinker;
291 if(this->GetUseShrinkImageFilter())
293 shrinker = ShrinkerType::New();
294 shrinkerFilter = shrinker.GetPointer();
298 resampleShrinker = ResampleShrinkerType::New();
299 typedef itk::LinearInterpolateImageFunction< OutputImageType, double >
300 LinearInterpolatorType;
301 typename LinearInterpolatorType::Pointer interpolator =
302 LinearInterpolatorType::New();
303 resampleShrinker->SetInterpolator( interpolator );
304 resampleShrinker->SetDefaultPixelValue( 0 );
305 shrinkerFilter = resampleShrinker.GetPointer();
308 caster->SetInput( inputPtr );
310 smoother->SetUseImageSpacing( false );
311 smoother->SetInput( caster->GetOutput() );
312 smoother->SetMaximumError( m_MaximumError );
314 shrinkerFilter->SetInput( smoother->GetOutput() );
316 unsigned int ilevel, idim;
317 unsigned int factors[ImageDimension];
318 double variance[ImageDimension];
320 for( ilevel = 0; ilevel < m_NumberOfLevels; ilevel++ )
322 this->UpdateProgress( static_cast<float>( ilevel ) /
323 static_cast<float>( m_NumberOfLevels ) );
325 // Allocate memory for each output
326 OutputImagePointer outputPtr = this->GetOutput( ilevel );
327 outputPtr->SetBufferedRegion( outputPtr->GetRequestedRegion() );
328 outputPtr->Allocate();
330 // compute shrink factors and variances
331 for( idim = 0; idim < ImageDimension; idim++ )
333 factors[idim] = m_Schedule[ilevel][idim];
334 variance[idim] = vnl_math_sqr( 0.5 *
335 static_cast<float>( factors[idim] ) );
338 if(!this->GetUseShrinkImageFilter())
340 typedef itk::IdentityTransform<double,OutputImageType::ImageDimension>
341 IdentityTransformType;
342 typename IdentityTransformType::Pointer identityTransform =
343 IdentityTransformType::New();
344 resampleShrinker->SetOutputParametersFromImage( outputPtr );
345 resampleShrinker->SetTransform(identityTransform);
349 shrinker->SetShrinkFactors(factors);
351 // use mini-pipeline to compute output
352 smoother->SetVariance( variance );
354 shrinkerFilter->GraftOutput( outputPtr );
356 // force to always update in case shrink factors are the same
357 shrinkerFilter->Modified();
358 shrinkerFilter->UpdateLargestPossibleRegion();
359 this->GraftNthOutput( ilevel, shrinkerFilter->GetOutput() );
366 template <class TInputImage, class TOutputImage>
368 SpatioTemporalMultiResolutionPyramidImageFilter<TInputImage, TOutputImage>
369 ::PrintSelf(std::ostream& os, itk::Indent indent) const
371 Superclass::PrintSelf(os,indent);
373 os << indent << "MaximumError: " << m_MaximumError << std::endl;
374 os << indent << "No. levels: " << m_NumberOfLevels << std::endl;
375 os << indent << "Schedule: " << std::endl;
376 os << m_Schedule << std::endl;
377 os << "Use ShrinkImageFilter= " << m_UseShrinkImageFilter << std::endl;
382 * GenerateOutputInformation
384 template <class TInputImage, class TOutputImage>
386 SpatioTemporalMultiResolutionPyramidImageFilter<TInputImage, TOutputImage>
387 ::GenerateOutputInformation()
390 // call the superclass's implementation of this method
391 Superclass::GenerateOutputInformation();
393 // get pointers to the input and output
394 InputImageConstPointer inputPtr = this->GetInput();
398 itkExceptionMacro( << "Input has not been set" );
401 const typename InputImageType::PointType&
402 inputOrigin = inputPtr->GetOrigin();
403 const typename InputImageType::SpacingType&
404 inputSpacing = inputPtr->GetSpacing();
405 const typename InputImageType::DirectionType&
406 inputDirection = inputPtr->GetDirection();
407 const typename InputImageType::SizeType& inputSize =
408 inputPtr->GetLargestPossibleRegion().GetSize();
409 const typename InputImageType::IndexType& inputStartIndex =
410 inputPtr->GetLargestPossibleRegion().GetIndex();
412 typedef typename OutputImageType::SizeType SizeType;
413 typedef typename SizeType::SizeValueType SizeValueType;
414 typedef typename OutputImageType::IndexType IndexType;
415 typedef typename IndexType::IndexValueType IndexValueType;
417 OutputImagePointer outputPtr;
418 typename OutputImageType::PointType outputOrigin;
419 typename OutputImageType::SpacingType outputSpacing;
421 IndexType outputStartIndex;
423 // we need to compute the output spacing, the output image size,
424 // and the output image start index
425 for(unsigned int ilevel = 0; ilevel < m_NumberOfLevels; ilevel++ )
427 outputPtr = this->GetOutput( ilevel );
428 if( !outputPtr ) { continue; }
430 for(unsigned int idim = 0; idim < OutputImageType::ImageDimension; idim++ )
432 const double shrinkFactor = static_cast<double>( m_Schedule[ilevel][idim] );
433 outputSpacing[idim] = inputSpacing[idim] * shrinkFactor;
435 outputSize[idim] = static_cast<SizeValueType>(
436 vcl_floor(static_cast<double>(inputSize[idim]) / shrinkFactor ) );
437 if( outputSize[idim] < 1 ) { outputSize[idim] = 1; }
439 outputStartIndex[idim] = static_cast<IndexValueType>(
440 vcl_ceil(static_cast<double>(inputStartIndex[idim]) / shrinkFactor ) );
442 //Now compute the new shifted origin for the updated levels;
443 const typename OutputImageType::PointType::VectorType outputOriginOffset
444 =(inputDirection*(outputSpacing-inputSpacing))*0.5;
445 for(unsigned int idim = 0; idim < OutputImageType::ImageDimension; idim++ )
447 outputOrigin[idim]=inputOrigin[idim]+outputOriginOffset[idim];
450 typename OutputImageType::RegionType outputLargestPossibleRegion;
451 outputLargestPossibleRegion.SetSize( outputSize );
452 outputLargestPossibleRegion.SetIndex( outputStartIndex );
454 outputPtr->SetLargestPossibleRegion( outputLargestPossibleRegion );
455 outputPtr->SetOrigin ( outputOrigin );
456 outputPtr->SetSpacing( outputSpacing );
457 outputPtr->SetDirection( inputDirection );//Output Direction should be same as input.
463 * GenerateOutputRequestedRegion
465 template <class TInputImage, class TOutputImage>
467 SpatioTemporalMultiResolutionPyramidImageFilter<TInputImage, TOutputImage>
468 ::GenerateOutputRequestedRegion(itk::DataObject * refOutput )
470 // call the superclass's implementation of this method
471 Superclass::GenerateOutputRequestedRegion( refOutput );
473 // find the index for this output
474 unsigned int refLevel = refOutput->GetSourceOutputIndex();
476 // compute baseIndex and baseSize
477 typedef typename OutputImageType::SizeType SizeType;
478 typedef typename SizeType::SizeValueType SizeValueType;
479 typedef typename OutputImageType::IndexType IndexType;
480 typedef typename IndexType::IndexValueType IndexValueType;
481 typedef typename OutputImageType::RegionType RegionType;
483 TOutputImage * ptr = static_cast<TOutputImage*>( refOutput );
486 itkExceptionMacro( << "Could not cast refOutput to TOutputImage*." );
489 unsigned int ilevel, idim;
491 if ( ptr->GetRequestedRegion() == ptr->GetLargestPossibleRegion() )
493 // set the requested regions for the other outputs to their
496 for( ilevel = 0; ilevel < m_NumberOfLevels; ilevel++ )
498 if( ilevel == refLevel ) { continue; }
499 if( !this->GetOutput(ilevel) ) { continue; }
500 this->GetOutput(ilevel)->SetRequestedRegionToLargestPossibleRegion();
505 // compute requested regions for the other outputs based on
506 // the requested region of the reference output
507 IndexType outputIndex;
509 RegionType outputRegion;
510 IndexType baseIndex = ptr->GetRequestedRegion().GetIndex();
511 SizeType baseSize = ptr->GetRequestedRegion().GetSize();
513 for( idim = 0; idim < TOutputImage::ImageDimension; idim++ )
515 unsigned int factor = m_Schedule[refLevel][idim];
516 baseIndex[idim] *= static_cast<IndexValueType>( factor );
517 baseSize[idim] *= static_cast<SizeValueType>( factor );
520 for( ilevel = 0; ilevel < m_NumberOfLevels; ilevel++ )
522 if( ilevel == refLevel ) { continue; }
523 if( !this->GetOutput(ilevel) ) { continue; }
525 for( idim = 0; idim < TOutputImage::ImageDimension; idim++ )
528 double factor = static_cast<double>( m_Schedule[ilevel][idim] );
530 outputSize[idim] = static_cast<SizeValueType>(
531 vcl_floor(static_cast<double>(baseSize[idim]) / factor ) );
532 if( outputSize[idim] < 1 ) { outputSize[idim] = 1; }
534 outputIndex[idim] = static_cast<IndexValueType>(
535 vcl_ceil(static_cast<double>(baseIndex[idim]) / factor ) );
539 outputRegion.SetIndex( outputIndex );
540 outputRegion.SetSize( outputSize );
542 // make sure the region is within the largest possible region
543 outputRegion.Crop( this->GetOutput( ilevel )->
544 GetLargestPossibleRegion() );
545 // set the requested region
546 this->GetOutput( ilevel )->SetRequestedRegion( outputRegion );
554 * GenerateInputRequestedRegion
556 template <class TInputImage, class TOutputImage>
558 SpatioTemporalMultiResolutionPyramidImageFilter<TInputImage, TOutputImage>
559 ::GenerateInputRequestedRegion()
561 // call the superclass' implementation of this method
562 Superclass::GenerateInputRequestedRegion();
564 // get pointers to the input and output
565 InputImagePointer inputPtr =
566 const_cast< InputImageType * >( this->GetInput() );
569 itkExceptionMacro( << "Input has not been set." );
572 // compute baseIndex and baseSize
573 typedef typename OutputImageType::SizeType SizeType;
574 typedef typename SizeType::SizeValueType SizeValueType;
575 typedef typename OutputImageType::IndexType IndexType;
576 typedef typename IndexType::IndexValueType IndexValueType;
577 typedef typename OutputImageType::RegionType RegionType;
579 unsigned int refLevel = m_NumberOfLevels - 1;
580 SizeType baseSize = this->GetOutput(refLevel)->GetRequestedRegion().GetSize();
581 IndexType baseIndex = this->GetOutput(refLevel)->GetRequestedRegion().GetIndex();
582 RegionType baseRegion;
585 for( idim = 0; idim < ImageDimension; idim++ )
587 unsigned int factor = m_Schedule[refLevel][idim];
588 baseIndex[idim] *= static_cast<IndexValueType>( factor );
589 baseSize[idim] *= static_cast<SizeValueType>( factor );
591 baseRegion.SetIndex( baseIndex );
592 baseRegion.SetSize( baseSize );
594 // compute requirements for the smoothing part
595 typedef typename TOutputImage::PixelType OutputPixelType;
596 typedef typename itk::GaussianOperator<OutputPixelType,ImageDimension> OperatorType;
598 OperatorType *oper = new OperatorType;
600 typename TInputImage::SizeType radius;
602 RegionType inputRequestedRegion = baseRegion;
605 for( idim = 0; idim < TInputImage::ImageDimension; idim++ )
607 oper->SetDirection(idim);
608 oper->SetVariance( vnl_math_sqr( 0.5 * static_cast<float>(
609 m_Schedule[refLevel][idim] ) ) );
610 oper->SetMaximumError( m_MaximumError );
611 oper->CreateDirectional();
612 radius[idim] = oper->GetRadius()[idim];
616 inputRequestedRegion.PadByRadius( radius );
618 // make sure the requested region is within the largest possible
619 inputRequestedRegion.Crop( inputPtr->GetLargestPossibleRegion() );
621 // set the input requested region
622 inputPtr->SetRequestedRegion( inputRequestedRegion );