[Insight-users] question on itkWarpImageFilter

Tue May 24 09:05:41 EDT 2011

Thank you for clarifying your point,
Yes you are right, the itk::WarpImageFilter uses the deformation field to map a
point from the output space to the input space. This is basically resides in the difference
between a transformation and a warping, as in transformation a point from the input space
is mapped to the output space with no guarantees for back transformations or size 
correspondence which is only natural since there are numerous types of transformations
with high levels of complexity and that would make it a fairly hard task to obtain the 
inverse transformation (as you pointed out).
On the opposite side, the whole point of warping is to map points from the output space to
the input space, hence the itk::WarpImageFilter does its work with inverse mapping, this
also avoids creation of holes and overlaps in the output image.
Best regards,
Dawood

>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<

--- On Tue, 5/24/11, Gao, Yi <gaoyi.cn at gmail.com> wrote:

From: Gao, Yi <gaoyi.cn at gmail.com>
Subject: Re: question on itkWarpImageFilter
To: "Dawood Masslawi" <masslawi at gmail.com>
Cc: insight-users at itk.org
Date: Tuesday, May 24, 2011, 8:10 AM

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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/pipermail/insight-users/attachments/20110524/0d98a388/attachment.htm>