[Insight-users] VectorImage Orientation from NRRD

[Daniel Betz]

Hello Gordon,

> This patch confuses two independent issues:
>
> 1) turning on ITK_USE_ORIENTED_IMAGE_DIRECTION.  It sounds like
> you're happy with this, specifically the resulting enhancements in
> orientation behavior to itkVectorImage.  Great.
>
> 2) itkNrrdImageIO.cxx's conversion of explicitly RAS orientation to
> the implicit LPS orientation of ITK.  You're understandably not happy
> with this.
>
> The patch you suggest confuses the issues by disabling (2) when (1)
> is off<comment>you mean on?</comment>, but the two things are
> completely unrelated, so I'm not going to commit this.
You nabbed me ;) Originally I planed a long post about the pro and
contra of sign flipping in file readers, and kind and audience of
documentation. As I saw that there is not enough interest for that,
I posted only the "patch" to show up two things:

    1. For a oriented itk::Image the sign flipping is even more terrible.
    2. Provocate some informative reactions from core developers.
       (Because I knew that it would interact in some ways with
        SpatialOrientationAdapter)

> The real solution is to enable explicit run-time control of whether
> the RAS -> LPS orientation conversion happens inside itkNrrdImageIO.
> This would replace how you used ITK_USE_ORIENTED_IMAGE_DIRECTION in
> your patch, with some new member variable (and would use "if ()"
> instead of "#ifndef").
That would be a nasty work around. I think there should be changed some
other things in ITK. I want to come back on this shortly.

> Before I try anything like that, though, I need others to comment on
>
> 1) Whether there are any existing examples of state variables in file
> input, so that at run-time users can vary the handling of data after
> its read from disk, and before its returned to ITK.
I nowhere found something like that.

> 2) Documentation of ITK's use of LPS for 3-D orientation
> specification.  Is it really the policy, or not?
ITK is LPS oriented. It is explicitly documented in the comments of some
fileReaders and implicitly with the SpacialOrientationAdapter.

But what that means can not be answered in a consistent way.
I want to show you something.

The OrientImageFilter uses PermuteAxesImageFilter and than
FlipImageFilter (FlipAboutOriginOff) to change the memory layout
of the image. (The final CastImageFilter isn't of interest for
the moment.)

Please take a short look on the following program.

<code>
#include <itkOrientImageFilter.h>

int main( int argc, char ** argv )
{
  // ITK Default Orientation is (Direction Cosines Identity):
  itk::SpatialOrientation::ValidCoordinateOrientationFlags RAI =
    itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_RAI;

  // ITK Default Orientation want to be:
  itk::SpatialOrientation::ValidCoordinateOrientationFlags LPS =
    itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_LPS;

  std::cout
    << "RAI Direction Cosines: " << std::endl
    << itk::SpatialOrientationAdapter().ToDirectionCosines(RAI) << std::endl
    << "LPS Direction Cosines: " << std::endl
    << itk::SpatialOrientationAdapter().ToDirectionCosines(LPS) << std::endl;

  typedef itk::Image< unsigned short, 3 > ImageType;
  ImageType::Pointer image = ImageType::New();
  ImageType::SizeType  size  = {{ 4, 4, 4 }};
  image->SetRegions(size);
  image->Allocate();
  image->FillBuffer(0);
  ImageType::IndexType markerIndex = {{ 0, 0, 0 }};
  ImageType::PixelType markerValue = 1;
  image->SetPixel( markerIndex, markerValue );
  std::cout
    << "Direction Cosines of Image:" << std::endl << image->GetDirection() << 
std::endl
    << "Value(IJK: 0 0 0) = " << image->GetPixel(markerIndex) << std::endl << 
std::endl;

  itk::OrientImageFilter<ImageType, ImageType>::Pointer orienter =
    itk::OrientImageFilter<ImageType, ImageType>::New();

  // direction of image is identity; so the following two lines are (mostly) 
equivalent
  orienter->UseImageDirectionOn();
  //orienter->SetGivenCoordinateOrientation(RAI);
  orienter->SetDesiredCoordinateOrientation(LPS);
  orienter->SetInput(image);
  orienter->Update();
  std::cout << orienter << std::endl;
  ImageType::Pointer newImage = orienter->GetOutput();

  ImageType::IndexType newMarkerIndex = {{ 3, 3, 3 }};
  std::cout << std::endl
    << "Direction Cosines of new Image:" << std::endl  << 
newImage->GetDirection() << std::endl
    << "Value(IJK: 3 3 3) = " << newImage->GetPixel(newMarkerIndex) << 
std::endl;
}
</code>

The output of the program.

<output>
RAI Direction Cosines:
1 0 0
0 1 0
0 0 1

LPS Direction Cosines:
-1 0 0
0 -1 0
0 0 -1

Direction Cosines of Image:
1 0 0
0 1 0
0 0 1

Value(IJK: 0 0 0) = 1

OrientImageFilter (0x812c4f0)
  RTTI typeinfo:   itk::OrientImageFilter<itk::Image<unsigned short, 3u>, 
itk::Image<unsigned short, 3u> >
  Reference Count: 2
  Modified Time: 19
  Debug: Off
  Observers:
    none
  Number Of Required Inputs: 1
  Number Of Required Outputs: 1
  Number Of Threads: 2
  ReleaseDataFlag: Off
  ReleaseDataBeforeUpdateFlag: Off
  Input 0: (0x812c268)
  Output 0: (0x812ec18)
  AbortGenerateData: Off
  Progress: 0.333333
  Multithreader:
    RTTI typeinfo:   itk::MultiThreader
    Reference Count: 1
    Modified Time: 7
    Debug: Off
    Observers:
      none
    Thread Count: 2
    Global Maximum Number Of Threads: 0
  Desired Coordinate Orientation: 590851 (LPS)
  Given Coordinate Orientation: 525570 (RAI)
  Use Image Direction: 1
  Permute Axes: [0, 1, 2]
  Flip Axes: [1, 1, 1]


Direction Cosines of new Image:
-1 0 0
0 -1 0
0 0 -1

Value(IJK: 3 3 3) = 1
</output>


First of all I want to declare two terms to separate things
which doesnt't belong together.

Implicit Orientation of Image:
I imagine a cube in front of me with the ijk origin in the
left (from my point of view) upper front corner. When I place
a head in this cube, that way that I look in the face of it,
the implicit orientation of the image is LIP.
This is the orientation needed for m_GivenCoordinateOrientation
(in form of ITK_COORDINATE_ORIENTATION_RSA) in OrientImageFilter
to change the memory layout in that way that the new implicit
orientation of it complies with m_DesiredCoordinateOrientation.

Explicit Orientation of Image:
Target space for TransformIndexToPhysicalPoint(). This is also
the space of the "space" field in nrrd files. It will be reached
after applying the signes of the "space directions" field.

We can now determine some peculiarities of ITK. First of all,
ITK encodes the orientation with the position of the space origin
in opposite to almost the rest of the world (DICOM LPS is
ITK_COORDINATE_ORIENTATION_RAI). That is well documented and not
a big issue.

When I load a image, I have two questions:
    1. How to go to world space
    2. What is my world space

ITK trys to answer both questions with the direction cosins.

The comments in the listed ImageIO classes of my previous post
speak about a conversion to LPS orientation (when we can;
Even though the meant orientation is called RAI in ITK).

Here LPS orientation means:
Calculate your world space with the direction cosins and it will be
LPS (explicit orientation).

The SpacialOrientationAdapter (used e.g. in OrientImageFilter)
maps the direction cosines to the orientation.

Here LPS orientation means:
Look on the direction cosins. If it is identity we are in
LPS (implicit orientation).

If we had implicit and explicit orientation codes in
itk::ImageBase we could compute the index-to-world transform
from this. (Don't take me to serious, I am only thinking loud.)

Btw.: Besides the design level bugs, there is a normal bug
in OrientImageFilter. If I call UseImageDirectionOn() and
than shift the object to cout (*before* calling Update()) the
displayed desired orientation is ever ITK_COORDINATE_ORIENTATION_RIP
(the default value from the OrientImageAdapter constructor).
The reason is, that m_UseImageDirection is regarded first in
GenerateOutputInformation() to set m_DesiredCoordinateOrientation
and so PrintSelf() prints the default value, that will not be
used by the filter on Update().

Corrections, comments ?

Regards
Daniel Betz
-------------- next part --------------
A non-text attachment was scrubbed...
Name: not available
Type: application/pgp-signature
Size: 189 bytes
Desc: This is a digitally signed message part.
URL: <http://www.itk.org/pipermail/insight-users/attachments/20081102/37eac565/attachment.pgp>