[Insight-users] Bug in interaction between SetDirection and SetOrigin/SetSpacing methods in ITKImage?

[Alex Taylor]

All,

I'm seeing something I wouldn't have expected after using the SetDirection method of the ITKImage object to switch the orientation of the X/Y axes. What I believe I'm seeing is that SetOrigin expects the origin to be specified in the Xworld,Yworld system established by SetDirection, while SetSpacing expects the spacing to be specified as {SpacingAlongRows,SpacingAlongColumns}. This feels like an inconsistent design. I'm curious whether this was an intention design decision or a bug that will be fixed in a future release. Here is an example set of reproduction steps.

In the reproduction steps below, I define a 3x3 Image in which X rows parallel to rows and Y runs parallel to columns. I want the origin to be at world X,Y location (1,4). I want an X spacing of 1 and Y spacing of 4. There is an ambiguity in the documentation and interface in ITK Image whether the numeric array inputs for Spacing and Origin are in Row,Column or X,Y order. Either would be fine, but I wouldn't have expected that one would be row,column and the other would be X,Y.

When I execute the code below, I see the following output:

Continuous index at world location (1,4):  [0, 0]
Continuous index at world location (2,4):  [0, 1]
Continuous index at world location (1,8):  [1, 0]
Continuous index at world location (2,8):  [1, 1]

I take this as a verification that I have defined my intended coordinate system, but that it was necessary to use SetSpacing in a way I wouldn't have expected.

Is this a bug I should report or is this the intended design?

Thanks,

Alex Taylor


#include "mex.h"
#include "matrix.h"
#include "itkImage.h"

/* Define a 3 row, 3 column image. Set Y direction to parallel to rows, and X direction to be parallel to columns. Define center of first pixel (origin) to be at world X,Y location (1,4). Define X spacing to be 1
   and Y spacing to be 4. */

    typedef itk::Image<double,2> imageType;
    imageType::Pointer image = imageType::New();

    imageType::IndexType start;
    imageType::SizeType  size;
    for (mwSize i = 0; i < 2; ++i)
    {
        start[i] = 0;
        size[i]  = 3;
   }
    imageType::RegionType region;
    region.SetSize(size);
    region.SetIndex(start);
    image->SetRegions(region);

    typedef imageType::DirectionType DirectionType;
    DirectionType axisDirection = image->GetDirection();
    axisDirection[0][0] = 0;
    axisDirection[0][1] = 1;
    axisDirection[1][0] = 1;
    axisDirection[1][1] = 0;
    image->SetDirection(axisDirection);

    double origin[2] = {1,4};
    double spacing[2] = {4,1};
    // THIS IS WHAT I BELIEVE IS A BUG. I WOULD HAVE THOUGHT SPACING WOULD BE IN X,Y ORDER, OR BE CONSISTENT WITH THE DEFINITION OF ORIGIN AT A MINIMUM.

    // Flip X and Y axis orientation so that X goes across rows and Y goes across columns
    image->SetDirection(axisDirection);
    image->SetOrigin(origin);
    image->SetSpacing(spacing);

    itk::ContinuousIndex< double, 2 > pixelIndex;
    imageType::PointType worldPointOfInterest;
    worldPointOfInterest[0] = 1;      // x position of the pixel
    worldPointOfInterest[1] = 4;      // y position of the pixel
    image->TransformPhysicalPointToContinuousIndex(worldPointOfInterest,pixelIndex);

    std::cout << "Continuous index at world location (1,4):  " << pixelIndex << std::endl;

    worldPointOfInterest[0] = 2;      // x position of the pixel
    worldPointOfInterest[1] = 4;      // y position of the pixel
    image->TransformPhysicalPointToContinuousIndex(worldPointOfInterest,pixelIndex);

    std::cout << "Continuous index at world location (2,4):  " << pixelIndex << std::endl;

    worldPointOfInterest[0] = 1;      // x position of the pixel
    worldPointOfInterest[1] = 8;      // y position of the pixel
    image->TransformPhysicalPointToContinuousIndex(worldPointOfInterest,pixelIndex);


    std::cout << "Continuous index at world location (1,8):  " << pixelIndex << std::endl;

    worldPointOfInterest[0] = 2;      // x position of the pixel
    worldPointOfInterest[1] = 8;      // y position of the pixel
    image->TransformPhysicalPointToContinuousIndex(worldPointOfInterest,pixelIndex);

    std::cout << "Continuous index at world location (2,8):  " << pixelIndex << std::endl;
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/pipermail/insight-users/attachments/20130221/6fab0bd5/attachment.htm>