[Insight-users] question on itkWarpImageFilter

[Gao, Yi]

Hi Dawood,

Thanks for the answer.

The reason I was asking is because that when previously dealing with
image transformation, I thought we need to know "where it is from"
from the transformation information, which tells "where to go".

For example, if image A is to be translated by 10 pixel to the RIGHT.
Then, we go through the new image, and for each pixel we find out
"where it is from": 10 pixel from the LEFT. Note in such cases, the
output image do NOT have to be the same size as the input one.
Moreover, if A is a simple binary image with an object in the center,
then in the transformed image, the object would appear to be a little
on the right (on the screen, and assuming the output region is the
same as the input one.)

Therefore I thought for transformation with general vector deformation
field, this would be the same. However getting the inverse transform
from a vector field is in general difficult. So I asked the previous
question.

However, for transformation with vector deformation field, it does not
seems to be working that way. The transformation information (vector
field) directly tells "where it is from". I did the following
experiment:

1. I generate a 2D image A of the width 100 and height 50. In the
middle (middle of the left-right) 1/3 region it has some content.

2. I generate a vector field D, the 1st component of each vector is
10, and the 2nd component is 0. Therefore I hope this vector field to
do the same as the translation above. And I expect the transformed
image, B, has the object moved to the RIGHT.

3. However, the result is opposites to that of the 2D Translation
Transformation. The object appears to the LEFT.

4. No problem, because this is consistent with the Software Guide:
p_in = p_out + deformationVector.

5. I think this confirms the behavior that the vector field tells
"where is from".

6. So seems that the we should say the deformation is defined on the
domain of the *output image*?



Here is my testing code:

#include "itkImage.h"
#include "itkImageFileWriter.h"
#include "itkWarpImageFilter.h"
#include "itkLinearInterpolateImageFunction.h"


template<typename itkImage_t>
typename itkImage_t::Pointer generate_image(int nx, int ny);

template<typename deformation_t>
typename deformation_t::Pointer generate_deformation_field(int nx, int ny);

template<typename itkImage_t>
void write_image(typename itkImage_t::Pointer im, std::string file_name);


int main( int argc, char * argv[] )
{
  const     unsigned int   Dimension = 2;

  typedef   float VectorComponentType;
  typedef   itk::Vector< VectorComponentType, Dimension > VectorPixelType;
  typedef   itk::Image< VectorPixelType,  Dimension >   DeformationFieldType;

  typedef   unsigned char  PixelType;
  typedef   itk::Image< PixelType,  Dimension >   ImageType;


  int nx = 100;
  int ny = 50;
  ImageType::Pointer img = generate_image<ImageType>(nx, ny);

  DeformationFieldType::Pointer deformationField =
generate_deformation_field<DeformationFieldType>(nx, ny);

  typedef itk::LinearInterpolateImageFunction< ImageType, double >
InterpolatorType;
  InterpolatorType::Pointer interpolator = InterpolatorType::New();

  typedef itk::WarpImageFilter< ImageType, ImageType,
DeformationFieldType  >  FilterType;
  FilterType::Pointer filter = FilterType::New();
  filter->SetInterpolator( interpolator );
  filter->SetOutputSpacing( deformationField->GetSpacing() );
  filter->SetOutputOrigin(  deformationField->GetOrigin() );
  filter->SetOutputDirection(  deformationField->GetDirection() );
  filter->SetDeformationField( deformationField );
  filter->SetInput( img );
  filter->Update();

  write_image<ImageType>(img, "a.nrrd");
  write_image<ImageType>(filter->GetOutput(), "b.nrrd");
  write_image<DeformationFieldType>(deformationField, "c.nrrd");

  return EXIT_SUCCESS;
}

template<typename itkImage_t>
typename itkImage_t::Pointer generate_image(int nx, int ny)
{
  typename itkImage_t::IndexType start = {{0, 0}};
  typename itkImage_t::SizeType size = {{nx, ny}};

  typename itkImage_t::RegionType region;
  region.SetSize(size);
  region.SetIndex(start);

  typename itkImage_t::Pointer img = itkImage_t::New();
  img->SetRegions(region);
  img->Allocate();
  img->FillBuffer(0);

  for (int ix = nx/3; ix < 2*nx/3; ++ix)
    {
      for (int iy = 0; iy < ny; ++iy)
        {
          typename itkImage_t::IndexType idx = {{ix, iy}};

          img->SetPixel(idx, static_cast<typename
itkImage_t::PixelType>(127*(sin(ix/2.0) + 1)));
        }
    }

  return img;
}



template<typename deformation_t>
typename deformation_t::Pointer generate_deformation_field(int nx, int ny)
{
  typename deformation_t::IndexType start = {{0, 0}};
  typename deformation_t::SizeType size = {{nx, ny}};

  typename deformation_t::RegionType region;
  region.SetSize(size);
  region.SetIndex(start);

  typename deformation_t::PixelType v;
  v[0] = 0.0;
  v[1] = 0.0;

  typename deformation_t::Pointer deformationField = deformation_t::New();
  deformationField->SetRegions(region);
  deformationField->Allocate();
  deformationField->FillBuffer(v);

  for (int ix = 0; ix < nx; ++ix)
    {
      for (int iy = 0; iy < ny; ++iy)
        {
          typename deformation_t::IndexType idx = {{ix, iy}};

          //v[0] = ix/5.0;
          v[0] = 10;
          v[1] = 0;

          deformationField->SetPixel(idx, v);
        }
    }

  return deformationField;
}


template<typename itkImage_t>
void write_image(typename itkImage_t::Pointer im, std::string file_name)
{
  typedef itk::ImageFileWriter< itkImage_t >  WriterType;
  typename WriterType::Pointer writer = WriterType::New();
  writer->SetFileName( file_name.c_str() );
  writer->SetInput( im );

  try
    {
      writer->Update();
    }
  catch( itk::ExceptionObject & excp )
    {
    std::cerr << "Exception thrown " << std::endl;
    std::cerr << excp << std::endl;
    }

  return;
}







On Tue, May 24, 2011 at 4:38 AM, Dawood Masslawi <masslawi at gmail.com> wrote:
>
> Hi Yi,
> Each vector in the deformation filed image has the same dimension as the input image,
> meaning that it has equal number of dimensions for each pixel corresponding to movement
> of the input pixel. Naturally the deformation field image would have the same size as
> the input image so it wold cover all the pixels in the input image.
> As for your third question, note that the input and output image would also have the same
> size, so using the deformation field it is clear that which pixels correspond to each other.
> Hope this answers your questions,
> Dawood
>
> >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
> >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
>
> Dear group,
>
> I'm a little confused about the itkWarpImageFilter....
>
> 1. I understand that in computing the deformed (output) image, we
> should go through the domain of the *output* image, and inversely
> trace back (according to the deformation vector field) to find the
> intensity in the input image.
>
> 2. However, on page 242 of the software guide, above equation 6.21, we
> see that the deformation vector field (as a vector image) should have
> the same size as the *input* image. So the deformation field is
> defined on the domain of the *input* image.
>
> 3. Then, go back to item 1, how do we achieve the inversion of the
> vector field so as it defines on the domain of the *output* image?
>
> Is there anything I misunderstand? Thanks for any hint!
>
> Best,
> yi