[Insight-users] ResampleImageFilter Producing Garbled Output When Image Size Changed

[Geoff Topping]

Hi.

I'm trying to use ResampleImageFilter to resize, crop or expand the volume of an image using ITK.  My input image is 128x128x95 voxels with spacings of 0.865759x0.865759x0.796, 16-bit signed integer Analyze75 format.

When I use the following code (see full listing below for context) to set the output image size and scaling, things work fine, and I get back an image very similar to my input image (although stored as floats instead of signed ints, and with the calibration factor ignored by the casting between formats):

    resampleFilter->SetSize(                inputImage->GetLargestPossibleRegion().GetSize());
    resampleFilter->SetOutputOrigin(        inputImage->GetOrigin());
    resampleFilter->SetOutputSpacing(       inputImage->GetSpacing());

When I switch to using hard coded image size and spacings (128x128x95 and 0.865759x0.865759x0.796), I also get an effectively identical output image:

    double spacing[] = {0.865759, 0.865759, 0.796};
    resampleFilter->SetOutputSpacing(spacing);
    
    double origin[] = {0.0, 0.0, 0.0};
    resampleFilter->SetOutputOrigin(origin);
    
    ImageType::SizeType size = {128, 128, 95};
    resampleFilter->SetSize(size);

However, if I change the hard-coded image size, I get garbled output.  Setting hard-coded size to 256x256x190, with or without dividing all the spacings by 2.0 results in most, but not all, of the slices of the resulting image being quite garbled, with apparent tilings of a few slices over and over, occasional mirroring of a series of slices, or what looks like the early frames of a high-speed camera observing an exploding object, instead of the original image contents.  Examining the resulting images, I see four copies of the (distorted) image in some of the transverse slices, tiled ul, ur, ll, lr.  In other slice orientations I see two or four tiled copies in a row.

Any ideas what might be going wrong or how I can fix or debug this?

Thanks,
Geoff Topping


My current code:


The program takes two filenames from the command line: input image and output image file name.  It reads the input image, applies an identity transform while resampling, and outputs the resulting image.

***********************************************************************

#include "itkVersorRigid3DTransform.h"
#include "itkIdentityTransform.h"

#include "itkNearestNeighborInterpolateImageFunction.h"
#include "itkLinearInterpolateImageFunction.h"

#include "itkImage.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"

#include "itkResampleImageFilter.h"
#include "itkShiftScaleImageFilter.h"


#include "vcl_cmath.h"

namespace {
    const unsigned int                                              Dimension = 3;
    typedef float                                                   PixelType;
    typedef itk::Image<PixelType, Dimension>                        ImageType;
    typedef itk::ImageFileReader<ImageType>                         ImageReaderType;
    typedef itk::ImageFileWriter<ImageType>                         ImageWriterType;

    typedef itk::VersorRigid3DTransform<double>                     RigidTransformType;
    typedef itk::IdentityTransform<double, Dimension>               IdentityTransformType;
    typedef IdentityTransformType                                   TransformType;

    typedef itk::LinearInterpolateImageFunction<ImageType, double>  LinearInterpolatorType;
    typedef itk::NearestNeighborInterpolateImageFunction<
                ImageType, double>                                  NearestInterpolatorType;
    typedef NearestInterpolatorType                                 InterpolatorType;    

    typedef itk::ResampleImageFilter<ImageType, ImageType>          ResampleFilterType;
}


int main( int argc, char *argv[] ) {
    if( argc < 2 ) {
        std::cerr << "Missing Parameters " << std::endl;
        std::cerr << "Usage: " << argv[0] << " fixedImageFile movingImageFile outputImagefile" << std::endl;
        return EXIT_FAILURE;
    }


    ImageReaderType::Pointer imageReader =          ImageReaderType::New();
    imageReader->SetFileName(argv[1]);
    imageReader->Update();

    
    TransformType::Pointer transform =              TransformType::New();

    InterpolatorType::Pointer interpolator =        InterpolatorType::New();


    ImageType::Pointer inputImage = imageReader->GetOutput();


    ResampleFilterType::Pointer resampleFilter =    ResampleFilterType::New();
    resampleFilter->SetInput(               inputImage);
    resampleFilter->SetTransform(           transform);
    resampleFilter->SetInterpolator(        interpolator);
    resampleFilter->SetDefaultPixelValue(   0);
    resampleFilter->SetOutputDirection(     inputImage->GetDirection());


#ifdef FROM_INPUT
    resampleFilter->SetSize(                inputImage->GetLargestPossibleRegion().GetSize());
    resampleFilter->SetOutputOrigin(        inputImage->GetOrigin());
    resampleFilter->SetOutputSpacing(       inputImage->GetSpacing());
#endif

#ifdef PRINT_INPUT_STUFF
    const ImageType::SpacingType& inputImageSpacing = inputImage->GetSpacing();
    std::cout << "input image spacings: " << inputImageSpacing[0] << ", "
                                          << inputImageSpacing[1] << ", "
                                          << inputImageSpacing[2] << std::endl;
    
    const ImageType::PointType& inputImageOrigin = inputImage->GetOrigin();
    std::cout << "input image origin: " << inputImageOrigin[0] << ", "
                                        << inputImageOrigin[1] << ", "
                                        << inputImageOrigin[2] << std::endl;

    const ImageType::SizeType& inputImageSize = inputImage->GetLargestPossibleRegion().GetSize();
    std::cout << "input image size: " << inputImageSize[0] << ", "
                                      << inputImageSize[1] << ", "
                                      << inputImageSize[2] << std::endl;
#endif

//#ifdef HARD_CODED
    double spacing[] = {0.865759, 0.865759, 0.796};
    resampleFilter->SetOutputSpacing(spacing);
    
    double origin[] = {0.0, 0.0, 0.0};
    resampleFilter->SetOutputOrigin(origin);
    
    ImageType::SizeType size = {128, 128, 95};
    resampleFilter->SetSize(size);
//#endif


    ImageWriterType::Pointer writer =               ImageWriterType::New();
    writer->SetFileName(argv[2]);
    writer->SetInput(resampleFilter->GetOutput());

    try { 
        writer->Update(); 
    } catch (itk::ExceptionObject& err) { 
        std::cerr << "ExceptionObject caught !" << std::endl; 
        std::cerr << err << std::endl; 
        return EXIT_FAILURE;
    } 

    return EXIT_SUCCESS;
}


_________________________________________________________________
Reinvent how you stay in touch with the new Windows Live Messenger.
http://go.microsoft.com/?linkid=9650731
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/pipermail/insight-users/attachments/20090408/b8f41c60/attachment-0001.htm>