[Insight-users] Bug in PyramidFIlter?

[sgerber]

Hi

I am using the RecursiveMultiResolutionPyramidImageFilter and it
introduces an offset of 1/2 pixel. Using the 
DiscreteGaussianImageFilter
and ResampleImageFilter individually recursively does not introduce 
this
offset.

I attached the "hand-made" pyramid code below.

Kind regards
Sam




   std::string outFile = "pyramid-";

   ImageType::Pointer input = //your image;


   unsigned int numberOfLevels = 5;

   for(int i = numberOfLevels-1; i >= 0; i--){
     std::stringstream ss;
     ss << outFile << i << ".nrrd";
     std::cout << "Writing " << ss.str() << std::endl;
     typedef itk::ImageFileWriter<ImageType> FileWriterType;
     FileWriterType::Pointer writer = FileWriterType::New();
     writer->SetFileName(ss.str());
     writer->SetInput( input );
     writer->Update();

     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 outputSpacing;
     outputSpacing[0] = input->GetSpacing()[0] *
(static_cast<double>(inputSize[0]) / 
static_cast<double>(outputSize[0]));
     outputSpacing[1] = input->GetSpacing()[1] *
(static_cast<double>(inputSize[1]) / 
static_cast<double>(outputSize[1]));

     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(input->GetOrigin() );
     resample->SetOutputSpacing(outputSpacing);
     resample->SetTransform(TransformType::New());
     resample->UpdateLargestPossibleRegion();
     resample->Update();

     input = resample->GetOutput();

   }

   return EXIT_SUCCESS;