[Insight-users] Re: FMM spacing question : Fast Marching : Resampling

[Luis Ibanez]

This is a multi-part message in MIME format.
--------------010406080702040505040803
Content-Type: text/plain; charset=us-ascii; format=flowed
Content-Transfer-Encoding: 7bit

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.
>>
>>

--------------010406080702040505040803
Content-Type: text/plain;
 name="ResampleVolumes.cxx"
Content-Transfer-Encoding: 7bit
Content-Disposition: inline;
 filename="ResampleVolumes.cxx"

/*=========================================================================

  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;
}


--------------010406080702040505040803--