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 __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"
30 #include "vnl/vnl_math.h"
38 template <class TInputImage, class TOutputImage>
39 RecursiveSpatioTemporalMultiResolutionPyramidImageFilter<TInputImage, TOutputImage>
40 ::RecursiveSpatioTemporalMultiResolutionPyramidImageFilter()
42 this->Superclass::m_UseShrinkImageFilter = true;
48 template <class TInputImage, class TOutputImage>
50 RecursiveSpatioTemporalMultiResolutionPyramidImageFilter<TInputImage, TOutputImage>
54 if( !this->IsScheduleDownwardDivisible( this->GetSchedule() ) )
56 // use the Superclass implemenation
57 this->Superclass::GenerateData();
61 // Get the input and output pointers
62 InputImageConstPointer inputPtr = this->GetInput();
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;
69 typedef itk::ImageToImageFilter<TOutputImage,TOutputImage> ImageToImageType;
70 typedef itk::ResampleImageFilter<TOutputImage,TOutputImage> ResampleShrinkerType;
71 typedef itk::ShrinkImageFilter<TOutputImage,TOutputImage> ShrinkerType;
73 typename CasterType::Pointer caster = CasterType::New();
74 typename CopierType::Pointer copier = CopierType::New();
75 typename SmootherType::Pointer smoother = SmootherType::New();
78 typename ImageToImageType::Pointer shrinkerFilter;
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;
85 if(this->GetUseShrinkImageFilter())
87 shrinker = ShrinkerType::New();
88 shrinkerFilter = shrinker.GetPointer();
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();
109 unsigned int factors[ImageDimension];
110 double variance[ImageDimension];
113 OutputImagePointer outputPtr;
114 OutputImagePointer swapPtr;
115 typename TOutputImage::RegionType LPRegion;
117 smoother->SetUseImageSpacing( false );
118 smoother->SetMaximumError( this->GetMaximumError() );
119 shrinkerFilter->SetInput( smoother->GetOutput() );
122 // recursively compute outputs starting from the last one
123 for( ilevel = this->GetNumberOfLevels() - 1; ilevel > -1; ilevel--)
126 this->UpdateProgress( 1.0 - static_cast<float>( 1 + ilevel ) /
127 static_cast<float>( this->GetNumberOfLevels() ) );
129 // Allocate memory for each output
130 outputPtr = this->GetOutput( ilevel );
131 outputPtr->SetBufferedRegion( outputPtr->GetRequestedRegion() );
132 outputPtr->Allocate();
134 // cached a copy of the largest possible region
135 LPRegion = outputPtr->GetLargestPossibleRegion();
137 // Check shrink factors and compute variances
139 for( idim = 0; idim < ImageDimension; idim++ )
141 if( ilevel == static_cast<int>(this->GetNumberOfLevels()) - 1)
143 factors[idim] = this->GetSchedule()[ilevel][idim];
147 factors[idim] = this->GetSchedule()[ilevel][idim] /
148 this->GetSchedule()[ilevel+1][idim];
150 variance[idim] = vnl_math_sqr( 0.5 *
151 static_cast<float>( factors[idim] ) );
152 if( factors[idim] != 1 )
158 variance[idim] = 0.0;
163 if( allOnes && ilevel == static_cast<int>(this->GetNumberOfLevels()) - 1 )
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();
174 swapPtr = caster->GetOutput();
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();
189 swapPtr = copier->GetOutput();
194 if( ilevel == static_cast<int>(this->GetNumberOfLevels()) - 1 )
196 // use caster -> smoother -> shrinker piepline
197 caster->SetInput( inputPtr );
198 smoother->SetInput( caster->GetOutput() );
202 // use smoother -> shrinker pipeline
203 smoother->SetInput( swapPtr );
206 smoother->SetVariance( variance );
208 // shrinker->SetShrinkFactors( factors );
209 // shrinker->GraftOutput( outputPtr );
210 if(!this->GetUseShrinkImageFilter())
212 resampleShrinker->SetOutputParametersFromImage(outputPtr);
216 shrinker->SetShrinkFactors(factors);
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();
226 swapPtr = shrinkerFilter->GetOutput();
230 // graft pipeline output back onto this filter's output
231 swapPtr->SetLargestPossibleRegion( LPRegion );
232 this->GraftNthOutput( ilevel, swapPtr );
234 // disconnect from pipeline to stop cycle
235 swapPtr->DisconnectPipeline();
244 template <class TInputImage, class TOutputImage>
246 RecursiveSpatioTemporalMultiResolutionPyramidImageFilter<TInputImage, TOutputImage>
247 ::PrintSelf(std::ostream& os, itk::Indent indent) const
249 Superclass::PrintSelf(os,indent);
254 * GenerateOutputRequestedRegion
256 template <class TInputImage, class TOutputImage>
258 RecursiveSpatioTemporalMultiResolutionPyramidImageFilter<TInputImage, TOutputImage>
259 ::GenerateOutputRequestedRegion(itk::DataObject * ptr )
262 // call the superclass's implementation of this method
263 Superclass::GenerateOutputRequestedRegion( ptr );
265 TOutputImage * refOutputPtr = static_cast<TOutputImage*>( ptr );
268 itkExceptionMacro( << "Could not cast ptr to TOutputImage*." );
271 // find the index for this output
272 unsigned int refLevel;
273 refLevel = refOutputPtr->GetSourceOutputIndex();
275 typedef typename TOutputImage::PixelType OutputPixelType;
276 typedef itk::GaussianOperator<OutputPixelType,ImageDimension> OperatorType;
278 OperatorType * oper = new OperatorType;
279 oper->SetMaximumError( this->GetMaximumError() );
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;
288 unsigned int factors[ImageDimension];
290 typename TInputImage::SizeType radius;
292 RegionType requestedRegion;
293 SizeType requestedSize;
294 IndexType requestedIndex;
296 // compute requested regions for lower levels
297 for( ilevel = refLevel + 1; ilevel < static_cast<int>(this->GetNumberOfLevels());
301 requestedRegion = this->GetOutput( ilevel - 1 )->GetRequestedRegion();
302 requestedSize = requestedRegion.GetSize();
303 requestedIndex = requestedRegion.GetIndex();
305 for( idim = 0; idim < static_cast<int>(ImageDimension); idim++ )
307 factors[idim] = this->GetSchedule()[ilevel-1][idim] / this->GetSchedule()[ilevel][idim];
309 // take into account shrink component
310 requestedSize[idim] *= static_cast<SizeValueType>(factors[idim]);
311 requestedIndex[idim] *= static_cast<IndexValueType>(factors[idim]);
313 // take into account smoothing component
314 if( factors[idim] > 1 )
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];
328 requestedRegion.SetSize( requestedSize );
329 requestedRegion.SetIndex( requestedIndex );
330 requestedRegion.PadByRadius( radius );
331 requestedRegion.Crop( this->GetOutput(ilevel)->
332 GetLargestPossibleRegion() );
334 this->GetOutput(ilevel)->SetRequestedRegion( requestedRegion );
339 // compute requested regions for higher levels
340 for( ilevel = refLevel - 1; ilevel > -1; ilevel-- )
342 requestedRegion = this->GetOutput( ilevel + 1 )->GetRequestedRegion();
343 requestedSize = requestedRegion.GetSize();
344 requestedIndex = requestedRegion.GetIndex();
346 for( idim = 0; idim < static_cast<int>(ImageDimension); idim++ )
349 factors[idim] = this->GetSchedule()[ilevel][idim] / this->GetSchedule()[ilevel+1][idim];
351 // take into account smoothing component
352 if( factors[idim] > 1 )
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];
365 requestedSize[idim] -= static_cast<SizeValueType>(
367 requestedIndex[idim] += radius[idim];
369 // take into account shrink component
370 requestedSize[idim] = static_cast<SizeValueType>( vcl_floor(
371 static_cast<double>(requestedSize[idim]) /
372 static_cast<double>(factors[idim]) ) );
373 if( requestedSize[idim] < 1 )
375 requestedSize[idim] = 1;
377 requestedIndex[idim] = static_cast<IndexValueType>( vcl_ceil(
378 static_cast<double>(requestedIndex[idim]) /
379 static_cast<double>(factors[idim]) ) );
383 requestedRegion.SetSize( requestedSize );
384 requestedRegion.SetIndex( requestedIndex );
385 requestedRegion.Crop( this->GetOutput(ilevel)->
386 GetLargestPossibleRegion() );
388 this->GetOutput(ilevel)->SetRequestedRegion( requestedRegion );
398 * GenerateInputRequestedRegion
400 template <class TInputImage, class TOutputImage>
402 RecursiveSpatioTemporalMultiResolutionPyramidImageFilter<TInputImage, TOutputImage>
403 ::GenerateInputRequestedRegion()
406 // call the superclass' implementation of this method
407 Superclass::GenerateInputRequestedRegion();
409 // get pointers to the input and output
410 InputImagePointer inputPtr =
411 const_cast< InputImageType *>( this->GetInput() );
414 itkExceptionMacro( << "Input has not been set." );
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;
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;
430 for( idim = 0; idim < ImageDimension; idim++ )
432 unsigned int factor = this->GetSchedule()[refLevel][idim];
433 baseIndex[idim] *= static_cast<IndexValueType>( factor );
434 baseSize[idim] *= static_cast<SizeValueType>( factor );
436 baseRegion.SetIndex( baseIndex );
437 baseRegion.SetSize( baseSize );
439 // compute requirements for the smoothing part
440 typedef typename TOutputImage::PixelType OutputPixelType;
441 typedef itk::GaussianOperator<OutputPixelType,ImageDimension> OperatorType;
443 OperatorType *oper = new OperatorType;
445 typename TInputImage::SizeType radius;
447 RegionType inputRequestedRegion = baseRegion;
450 for( idim = 0; idim < TInputImage::ImageDimension; idim++ )
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 )
466 inputRequestedRegion.PadByRadius( radius );
468 // make sure the requested region is within the largest possible
469 inputRequestedRegion.Crop( inputPtr->GetLargestPossibleRegion() );
471 // set the input requested region
472 inputPtr->SetRequestedRegion( inputRequestedRegion );