[Insight-users] Bug in PyramidFIlter?

[Bradley Lowekamp]

Hello,

Keep in mind that in ITK's image representation that the origin is at the center of 0-indexed pixel and the extent goes a half pixel before this.

This offset is intentional. The Pyramid filter uses the RecursiveGaussian filters, along with the ShrinkImage Filter. If the image is divisible by the shrink factor then the edge to edge extent is maintained when the origin is adjusted. If not the origin and samples are adjusted such that the center of the images align, which is the preferred approach for neuro-imaging registration.

In your included code your origin pixels is "sampling" more image space outside the original image extent than the current approach. This may be visible by setting the DefaultPixelValue in the resample image filter to 0 and having an image of 1s.

I myself, being an author of SimpleITK, would use SimpelITK in python to play around with the filters you are using to get a better understanding of this through interactive usage.

On Jul 13, 2013, at 9:05 PM, Bradley Lowekamp <brad at lowekamp.net> wrote:

> Hello,
> 
> Keep in mind that in ITK's image representation that the origin is at the center of 0-indexed pixel and the extent goes a half pixel before this.
> 
> This offset is intentional. The Pyramid filter uses the RecursiveGaussian filters, along with the ShrinkImage Filter. If the image is divisible by the shrink factor then the edge to edge extent is maintained when the origin is adjusted. If not the origin and samples are adjusted such that the center of the images align, which is the preferred approach for neuro-imaging registration.
> 
> In your included code your origin pixels is "sampling" more image space outside the original image extent than the current approach. This may be visible by setting the DefaultPixelValue in the resample image filter to 0 and having an image of 1s.
> 
> I myself, being an author of SimpleITK, would use SimpelITK in python to play around with the filters you are using to get a better understanding of this through interactive usage.
> 
> 
> Brad
> 
> On Jul 13, 2013, at 2:37 PM, sgerber <sgerber at cs.utah.edu> wrote:
> 
>> 
>> 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;
>> _____________________________________
>> Powered by www.kitware.com
>> 
>> Visit other Kitware open-source projects at
>> http://www.kitware.com/opensource/opensource.html
>> 
>> Kitware offers ITK Training Courses, for more information visit:
>> http://www.kitware.com/products/protraining.php
>> 
>> Please keep messages on-topic and check the ITK FAQ at:
>> http://www.itk.org/Wiki/ITK_FAQ
>> 
>> Follow this link to subscribe/unsubscribe:
>> http://www.itk.org/mailman/listinfo/insight-users
>