[Insight-users] Problem with bspline interpolation and VTKImageToImageFilter

[Michael Kuhn]

Hi,

I have a problem using the bspline interpolation, which I think is 
related to the VTKImageToImageFilter. If I'm trying to transform an 
image (which originates from a VTKImageToImageFilter) using a resample 
image filter and a bspline interpolator, the output is the default pixel 
value for all the pixels. Using a linear interpolator, the problem does 
not occur. Since I've already successfully transformed images with a 
bspline interpolator that did not originate from a 
VTKImageToImageFilter, I think the problem is related to that one.

I have quickly set up a  simple program that demonstrates the problem 
(see code below). It generates vtk image data (a cube) and transforms it 
by an identity transform. Using a linear interpolator, the output is the 
original cube, using a bspline interpolator it consists only of 1 values 
(which is the default pixel value) (comment and uncomment the BSPLINE 
BLOCK and LINEAR BLOCK to see the difference):

Thanks,

Michael



#include "itkImage.h"
#include "itkLinearInterpolateImageFunction.h"
#include "itkBSplineInterpolateImageFunction.h"
#include "itkVTKImageToImageFilter.h"
#include "itkResampleImageFilter.h"
#include "itkAffineTransform.h"

int main(int argc, char** argv) {
    vtkImageData* data = vtkImageData::New();
    data->SetDimensions(12, 12, 12);
    data->SetScalarTypeToShort();
    data->AllocateScalars();
    for (int x = 0; x < 12; x++) {
        for (int y = 0; y < 12; y++) {
            for (int z = 0; z < 12; z++) {
                short* tmp = (short*)data->GetScalarPointer(x, y, z);
                if (x >= 2 && x < 10 && y >= 2 && y < 10 && z >= 2 && z 
< 10)
                {
                    *tmp = 5;
                } else {
                    *tmp = 0;
                }
            }
        }
    }

    typedef itk::Image<short, 3> ImageType;
    typedef itk::VTKImageToImageFilter<ImageType> VTK2ITKType;
    VTK2ITKType::Pointer vtk2itk = VTK2ITKType::New();
    vtk2itk->GetImporter()->SetReleaseDataFlag(true);
    vtk2itk->SetInput(data);
    vtk2itk->GetImporter()->Update();
    ImageType::Pointer image1 = vtk2itk->GetImporter()->GetOutput();


    typedef itk::ResampleImageFilter<ImageType, ImageType> 
ResampleFilterType;
    ResampleFilterType::Pointer resampleFilter = ResampleFilterType::New();
   
   
    // BSPLINE BLOCK
    typedef itk::BSplineInterpolateImageFunction<ImageType, double, double>
        BSplineType;
    BSplineType::Pointer interpolator = BSplineType::New();
    interpolator->SetSplineOrder(3);
   
   
    /*
    // LINEAR BLOCK
    typedef itk::LinearInterpolateImageFunction<ImageType, double>
        LinearType;
    LinearType::Pointer interpolator = LinearType::New();
    */   

    typedef itk::AffineTransform<double, 3> AffineType;
    AffineType::Pointer affineTransform = AffineType::New();

    AffineType::ParametersType affineParams(12);
    affineParams[0] = 1;
    affineParams[1] = 0;
    affineParams[2] = 0;
    affineParams[3] = 0;
    affineParams[4] = 1;
    affineParams[5] = 0;
    affineParams[6] = 0;
    affineParams[7] = 0;
    affineParams[8] = 1;
    affineParams[9] = 0;
    affineParams[10] = 0;
    affineParams[11] = 0;

    affineTransform->SetParameters(affineParams);
    resampleFilter->SetInterpolator(interpolator);
    resampleFilter->SetTransform(affineTransform);
    resampleFilter->SetInput(image1);

    resampleFilter->SetOutputOrigin(image1->GetOrigin());
    resampleFilter->SetOutputSpacing(image1->GetSpacing());
    resampleFilter->SetSize(image1->GetBufferedRegion().GetSize());
    resampleFilter->SetDefaultPixelValue(1);
    resampleFilter->Update();
    ImageType::Pointer image2 = resampleFilter->GetOutput();

    for (int x = 0; x < 12; x++) {
        for (int y = 0; y < 12; y++) {
            for (int z = 0; z < 12; z++) {
                ImageType::IndexType idx;
                idx[0] = x;
                idx[1] = y;
                idx[2] = z;
                cout << image2->GetPixel(idx) << " ";
            }
            cout << endl;
        }
        cout << endl;
        cout << endl;
    }

    return 0;
}