[Insight-users] Re: FMM spacing question : Fast Marching : Resampling
: Errata
Luis Ibanez
luis . ibanez at kitware . com
Wed, 26 Nov 2003 11:40:26 -0500
Hi Erik,
A quick fix on the code I sent you yesterday,
The Gaussian smoothing filters in X and Y were
missing the SetDirection() methods. Therefore
they were both smoothing along X :-/
The fix is simply to add the lines
> smootherX->SetDirection( 0 );
> smootherY->SetDirection( 1 );
>
in line 101.
For convenience, this code has been commited
in the Examples under
Insight/Examples/Filtering/ResampleVolumesToBeIsotropic.cxx
Please let us know if you find any problem,
Thanks
Luis
---------------------------------
Luis Ibanez wrote:
> Hi Erik,
>
> The attached program will resample 3D data in order
> to produced isotropic volumes. It is performing first
> a Gaussian smoothing along the X and Y directions, then
> it uses the Resample image filter for resampling the
> data.
>
> You may find it useful...
>
> Note that this process produces correct spacing and
> origin. The output image may then be passed to any
> of the LevelSet methods.
>
> Regards,
>
>
> Luis
>
>
> --------------------------------
> Luis Ibanez wrote:
>
>>
>> Hi Erik,
>>
>> Thanks for pointing this out.
>>
>> You are right, the FastMarching filter doesn't take
>> the spacing of the image into account. This is also
>> true for most of the level set filters. This behavior
>> is the result of a design decision. Considering the spacing
>> will introduce a hit in performance on the process of solving
>> the PDE associated with these filters.
>>
>> Notice that there is a difference between 'using' the
>> spacing and 'carring' the spacing. In the first case,
>> the spacing should be use as part of the computations
>> performed by the filter, in the second case, the filter
>> simply copies the spacing values from the input to the
>> output. The FastMarching filter is not doing one nor
>> the other. It would be easy to do just the 'carrying'
>> part, however you must be aware that the spacing is not
>> used during the computation. This is quite important since
>> in some cases the output of the fastmarching filter is
>> used as a distance map.
>>
>>
>> You could also use the ChangeInformation filter in order to
>> inject the origin and spacing information from the input
>> image into the output image.
>> http://www . itk . org/Insight/Doxygen/html/classitk_1_1ChangeInformationImageFilter . html
>>
>>
>>
>> Unless your image is significantly big, the best solution
>> at this point is for you to use the ResampleImageFilter and
>> resample your image in order to get isotropic pixels.
>> http://www . itk . org/Insight/Doxygen/html/classitk_1_1ResampleImageFilter . html
>>
>>
>>
>> Regards,
>>
>>
>> Luis
>>
>>
>> --------------------
>> Erik Bekkers wrote:
>>
>>> Luis,
>>>
>>> This should be a really quick question: I'm using the
>>> Fast Marching Filter on a 3D image w/ spacing .8 .8
>>> 2.0, but the output from the FMM filter generates an
>>> image that is spacing 1 1 1; it also seems to be
>>> computing (values for distance etc.) with these
>>> unitary spacing values. Other pipelines of filters
>>> maintain the correct spacing.... how do I tell FMM to
>>> use the spacing in the image file it was given?
>>>
>>> THanks,
>>>
>>> Erik.
>>>
>>>
>
> ------------------------------------------------------------------------
>
> /*=========================================================================
>
> Program: Insight Segmentation & Registration Toolkit
> Module: $RCSfile: ResampleImageFilter3.cxx,v $
> Language: C++
> Date: $Date: 2003/09/10 14:29:55 $
> Version: $Revision: 1.15 $
>
> Copyright (c) Insight Software Consortium. All rights reserved.
> See ITKCopyright.txt or http://www . itk . org/HTML/Copyright . htm for details.
>
> This software is distributed WITHOUT ANY WARRANTY; without even
> the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
> PURPOSE. See the above copyright notices for more information.
>
> =========================================================================*/
>
>
> #include "itkImage.h"
> #include "itkImageFileReader.h"
> #include "itkImageFileWriter.h"
> #include "itkResampleImageFilter.h"
> #include "itkIdentityTransform.h"
> #include "itkLinearInterpolateImageFunction.h"
> #include "itkRecursiveGaussianImageFilter.h"
>
>
> int main( int argc, char * argv[] )
> {
> if( argc < 3 )
> {
> std::cerr << "Usage: " << std::endl;
> std::cerr << argv[0] << " inputImageFile outputImageFile" << std::endl;
> return 1;
> }
>
> const unsigned int Dimension = 3;
>
> typedef unsigned short InputPixelType;
> typedef float InternalPixelType;
> typedef unsigned short OutputPixelType;
>
> typedef itk::Image< InputPixelType, Dimension > InputImageType;
> typedef itk::Image< InternalPixelType, Dimension > InternalImageType;
> typedef itk::Image< OutputPixelType, Dimension > OutputImageType;
>
> typedef itk::ImageFileReader< InputImageType > ReaderType;
> typedef itk::ImageFileWriter< OutputImageType > WriterType;
>
> ReaderType::Pointer reader = ReaderType::New();
> WriterType::Pointer writer = WriterType::New();
>
> reader->SetFileName( argv[1] );
> writer->SetFileName( argv[2] );
>
> try
> {
> reader->Update();
> }
> catch( itk::ExceptionObject & excep )
> {
> std::cerr << "Exception catched !" << std::endl;
> std::cerr << excep << std::endl;
> }
>
>
> InputImageType::ConstPointer inputImage = reader->GetOutput();
>
> const double * inputSpacing = inputImage->GetSpacing();
>
> const double isoSpacing = sqrt( inputSpacing[2] * inputSpacing[0] );
>
> typedef itk::RecursiveGaussianImageFilter<
> InputImageType,
> InternalImageType > GaussianFilter1Type;
>
> typedef itk::RecursiveGaussianImageFilter<
> InternalImageType,
> InternalImageType > GaussianFilter2Type;
>
> GaussianFilter1Type::Pointer smootherX = GaussianFilter1Type::New();
> GaussianFilter2Type::Pointer smootherY = GaussianFilter2Type::New();
>
> smootherX->SetInput( reader->GetOutput() );
> smootherY->SetInput( smootherX->GetOutput() );
>
> smootherX->SetSigma( isoSpacing );
> smootherY->SetSigma( isoSpacing );
>
> try
> {
> smootherY->Update();
> }
> catch( itk::ExceptionObject & excep )
> {
> std::cerr << "Exception catched !" << std::endl;
> std::cerr << excep << std::endl;
> }
>
>
> InternalImageType::ConstPointer smoothedImage = smootherY->GetOutput();
>
> typedef itk::ResampleImageFilter<
> InternalImageType, OutputImageType > ResampleFilterType;
>
> ResampleFilterType::Pointer resampler = ResampleFilterType::New();
>
> typedef itk::IdentityTransform< double, Dimension > TransformType;
>
> typedef itk::LinearInterpolateImageFunction<
> InternalImageType, double > InterpolatorType;
>
> InterpolatorType::Pointer interpolator = InterpolatorType::New();
>
> resampler->SetInterpolator( interpolator );
>
> resampler->SetDefaultPixelValue( 1000 );
>
>
> double spacing[ Dimension ];
>
> spacing[0] = isoSpacing;
> spacing[1] = isoSpacing;
> spacing[2] = isoSpacing;
>
> resampler->SetOutputSpacing( spacing );
>
> // Use the same origin
> resampler->SetOutputOrigin( inputImage->GetOrigin() );
>
>
> InputImageType::SizeType inputSize = inputImage->GetLargestPossibleRegion().GetSize();
> InputImageType::SizeType size;
>
> size[0] = inputSize[0] * inputSpacing[0] / isoSpacing;
> size[1] = inputSize[1] * inputSpacing[1] / isoSpacing;
> size[2] = ( inputSize[2] - 1 ) * inputSpacing[2] / isoSpacing;
>
> resampler->SetSize( size );
>
> resampler->SetInput( smoothedImage );
>
> writer->SetInput( resampler->GetOutput() );
>
> TransformType::Pointer transform = TransformType::New();
>
> transform->SetIdentity();
>
> resampler->SetTransform( transform );
>
>
> try
> {
> writer->Update();
> }
> catch( itk::ExceptionObject & excep )
> {
> std::cerr << "Exception catched !" << std::endl;
> std::cerr << excep << std::endl;
> }
>
>
> return 0;
> }
>