[Insight-users] Bug in PyramidFIlter?

[Johnson, Hans J]

The edges of down-sampled images should not align in this case.

An image extent goes from 1/2 spacing before the first sample to 1/2
spacing beyond the last sample.

If the spacing changes, the effective physical space also changes.

Hans


-----Original Message-----
From: sgerber <sgerber at cs.utah.edu>
Date: Monday, July 15, 2013 12:58 PM
To: "insight-users at itk.org" <insight-users at itk.org>
Subject: Re: [Insight-users] Bug in PyramidFIlter?

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;
}
_____________________________________
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



________________________________
Notice: This UI Health Care e-mail (including attachments) is covered by the Electronic Communications Privacy Act, 18 U.S.C. 2510-2521, is confidential and may be legally privileged.  If you are not the intended recipient, you are hereby notified that any retention, dissemination, distribution, or copying of this communication is strictly prohibited.  Please reply to the sender that you have received the message in error, then delete it.  Thank you.
________________________________