[Insight-users] Bug in PyramidFIlter?

[sgerber]

Below is a self contained example that demonstrates the undesired 
behavior. Run the example and open the resulting images in Seg3D (or any 
other image viewer that renders the image spacing units / origins 
correctly and allows you to overlay the images), the edge in the image 
downsampled by RecursiveMultiResolutionPyramidImageFilter is not 
centered on the edge of the input image.

The culprit appears to be the ShrinkImageFilter. If the pyramid filter 
is set UseShrinkImageFilterOff() the outputs match. Maybe that is the 
intended behavior but for a multiresolution image pyramid I would claim 
that this is typically not the desired output.

Kind regards
Sam

#include "itkImage.h"
#include "itkImageFileWriter.h"
#include "itkRescaleIntensityImageFilter.h"
#include "itkRecursiveMultiResolutionPyramidImageFilter.h"
#include "itkResampleImageFilter.h"
#include "itkDiscreteGaussianImageFilter.h"




typedef itk::Image<double> ImageType;

ImageType::Pointer CreateImage(){
   // Create a black image with 2 white regions
   ImageType::Pointer image = ImageType::New();
   ImageType::IndexType start;
   start.Fill(0);

   ImageType::SizeType size;
   size.Fill(200);

   ImageType::RegionType region(start,size);
   image->SetRegions(region);
   image->Allocate();
   image->FillBuffer(0);

   // Make a square
   for(unsigned int r = 20; r < 80; r++)
   {
     for(unsigned int c = 30; c < 100; c++)
     {
       ImageType::IndexType pixelIndex;
       pixelIndex[0] = r;
       pixelIndex[1] = c;

       image->SetPixel(pixelIndex, 255);

     }
   }
   return image;
};


int main( int argc, char ** argv ){


   ImageType::Pointer image = CreateImage();
   {
     typedef itk::ImageFileWriter<ImageType> FileWriterType;
     FileWriterType::Pointer writer = FileWriterType::New();
     writer->SetFileName("test-input.nrrd");
     writer->SetInput( image );
     writer->Update();
   }

   unsigned int numberOfLevels = 3;

   {
     typedef itk::RecursiveMultiResolutionPyramidImageFilter<
       ImageType, ImageType>  
RecursiveMultiResolutionPyramidImageFilterType;
     RecursiveMultiResolutionPyramidImageFilterType::Pointer 
pyramidFilter =
       RecursiveMultiResolutionPyramidImageFilterType::New();
     pyramidFilter->SetInput(image);
     pyramidFilter->SetNumberOfLevels(numberOfLevels);
     pyramidFilter->Update();


     typedef itk::ImageFileWriter<ImageType> FileWriterType;
     FileWriterType::Pointer writer = FileWriterType::New();
     writer->SetFileName("test1.nrrd");
     writer->SetInput( pyramidFilter->GetOutput() );
     writer->Update();


   }


   {
     ImageType::Pointer input = image;
     for(int i = 1; i < numberOfLevels; i++){

       typedef itk::DiscreteGaussianImageFilter<ImageType, ImageType> 
Blur;
       Blur::Pointer blur = Blur::New();
       blur->SetUseImageSpacingOff();
       blur->SetVariance(1);
       blur->SetInput(input);

       blur->Update();
       input = blur->GetOutput();

       // Resize
       ImageType::SizeType inputSize = 
input->GetLargestPossibleRegion().GetSize();

       ImageType::SizeType outputSize;
       outputSize[0] = inputSize[0] / 2;
       outputSize[1] = inputSize[1] / 2;

       ImageType::SpacingType inputSpacing = input->GetSpacing();
       ImageType::SpacingType outputSpacing;
       outputSpacing[0] = inputSpacing[0] * 
(static_cast<double>(inputSize[0]) / 
static_cast<double>(outputSize[0]));
       outputSpacing[1] = inputSpacing[1] * 
(static_cast<double>(inputSize[1]) / 
static_cast<double>(outputSize[1]));

       ImageType::PointType origin = input->GetOrigin();
       origin[0] -= inputSpacing[0]/2 - outputSpacing[0]/2;
       origin[1] -= inputSpacing[1]/2 - outputSpacing[1]/2;

       typedef itk::IdentityTransform<double, 2> TransformType;
       typedef itk::ResampleImageFilter<ImageType, ImageType> 
ResampleImageFilterType;
       ResampleImageFilterType::Pointer resample = 
ResampleImageFilterType::New();
       resample->SetInput(input);
       resample->SetSize(outputSize);
       resample->SetOutputOrigin(origin);
       resample->SetOutputSpacing(outputSpacing);
       resample->SetTransform(TransformType::New());
       resample->UpdateLargestPossibleRegion();
       resample->Update();

       input = resample->GetOutput();

     }

     typedef itk::ImageFileWriter<ImageType> FileWriterType;
     FileWriterType::Pointer writer = FileWriterType::New();
     writer->SetFileName("test2.nrrd");
     writer->SetInput( input );
     writer->Update();
   }
   return EXIT_SUCCESS;
}