[Insight-developers] Flaw in ITK? : orientation of images

[Koen Van Leemput]

Hi,

In ITK, itk::Image comes with an IndexToPhysicalTransform and its inverse=
=20
PhysicalToIndexTransform that defines how coordinates in the index space =
map=20
into coordinates in the "real-world" physical space and vice versa. In=20
itk::Image, these transforms are AffineTransforms that take into account =
the=20
spacing of the voxels (scaling) and the position of the origin (translati=
on).

In many medical image processing applications, however, there is a need f=
or=20
additional positioning information, defining how the patient is oriented =
with=20
respect to the image coordinate system. Examples include flipping and=20
swapping of image axes, gantry tilt in CT images, interfacing with other=20
image processing tools (SPM, 3D Slicer, ...), and proper handling of vari=
ous=20
medical image formats (DICOM, SPM, ...).=20

Because of these reasons, I decided to go for it and derive a class from=20
itk::Image that builds m_IndexToPhysicalTransform (and=20
m_PhysicalToIndexTransform by) by composing the standard scaling and=20
translation with an additional affine transform that defines orientation.=
=20
Apart from some convenience methods to manipulate this orientation, I=20
overloaded itk::Image::RebuildTransforms() and itk::Image::CopyInformatio=
n(),=20
and voila, it works: the orientation information travels nicely through t=
he=20
ITK pipeline.

BUT... things are messed up by what seems to me a design flaw in=20
itk::ImageFunction. Indeed, rather than using=20
itk::Image::TransformPhysicalPointToContinuousIndex() to convert a point =
in=20
physical space into an index in the image grid, itk::ImageFunction for so=
me=20
reason provides its own implementation=20
itk::ImageFunction::ConvertPointToContinuousIndex that noisily ignores=20
whatever PhysicalToIndexTransform the image at hand actually uses. This=20
messes up all the registraton metrics and itk::ResampleImageFilter, where=
 the=20
following construction is typically used:

    m_FixedImage->TransformIndexToPhysicalPoint( inputIndex, inputPoint )=
;
    transformedPoint =3D m_Transform->TransformPoint( inputPoint );
    if( m_Interpolator->IsInsideBuffer( transformedPoint ) )
      {
      movingValue  =3D m_Interpolator->Evaluate( transformedPoint );
      [clip clip]     =09=09=09
      }

(Note that also itk::ImageFunction::IsInsideBuffer() doesn't handle thing=
s=20
properly here neither)

As I see it, itk::ImageFunction::IsInsideBuffer() and=20
itk::ImageFunction::ConvertPointToContinuousIndex() should use the=20
implementation provided by itk::Image. To avoid inefficiencies induced by=
=20
converting from physical space into index space twice, the following=20
construction should probably be considered:

    m_FixedImage->TransformIndexToPhysicalPoint( inputIndex, inputPoint )=
;
    transformedPoint =3D m_Transform->TransformPoint( inputPoint );
    m_MovingImage->TransformPhysicalPointToContinuousIndex(
                              transformedPoint, transformedIndex );
    if( m_Interpolator->IsInsideBuffer( transformedIndex ) )
      {
      movingValue  =3D m_Interpolator->Evaluate( transformedIndex );
      [clip clip]     =09=09=09
      }


More in general, is support in ITK for orientation planned any time in th=
e=20
future? After all, the fiddling around with itk::PermuteAxesImageFilters=20
itk::FlipImageFilters in the registration examples seems quite cumbersome=
=2E..

Cheers,

Koen