Support for orientation
This document was created based on the Slicer wiki Labs page found here originally written by @lassoan with input from @jcfr and later updated based on comments from @ken-martin , @demarle , @berkgeveci , @lisa-avila and @jcfr.
In medical image computing, we often have to work with images that are arbitrarily oriented in space. We’ll refer to images that can have arbitrary orientation as "oriented images", and images that have voxel coordinate axes aligned with the object coordinate axes as “axis-aligned images”.
Missing image directions from vtkImageData requires implementations of many workarounds in applications, which increase complexity, probability of errors, and decreases performance:
All file formats that we use (nrrd, mha, nifti, etc.) stores the axis directions, but currently when we load these images, the direction information is lost.
Most VTK filters could be easily extended to accept oriented images as input and/or generate oriented images as output, but currently they only work with axis-aligned images. This requires complicated and often computationally expensive transformations to perform computations in the voxel coordinate system of the volume and then transform the resulting data back to world coordinate system. For example:
- Image->polydata conversion: vtkContourFilter results has to be transformed using vtkPolyDataFilter
- Polydata->image conversion: we need to pre-transform the polydata to voxel coordinate system before we can call vtkPolyDataToImageStencil)
When a filter takes multiple images as inputs then those cannot be used with oriented images. There is no workaround, other then resample the images or implement new filters. For example, vtkGridTransform or vtkBSplineTransform classes would need to work with oriented images (coefficients are stored in vtkImageData) because the DICOM standard requires transformation with arbitrarily oriented grids, but if the input is an oriented image as well, there is no common voxel coordinate system where all the input data could be transformed to.
Some filters have additional inputs to be able to pass axis directions (such as Filters that are often used with oriented images (such as vtkImageReslice), but collecting the direction information separately and then injecting it into these filters correctly makes these filters very difficult to use.
While actors can make the volumes appear in arbitrary orientation, the orientation information has to be collected from the data reader and sent to the actors manually.
Linear transforms cannot be applied on vtkImageData because transformation result cannot be stored in a vtkImageData (when the transformation included rotation).
Different software applications deal with this in various ways – none of them are good solutions:
3D Slicer does not use vtkImageData origin and spacing fields and just stores the VoxelToObject transform (origin, spacing, direction) separately in a 4x4 matrix. This makes implementing VTK-based image processing pipelines quite complicated to implement.
Paraview ignores this problem and assumes that all images are axis-aligned. This makes it challenging for users to use Paraview for their medical images that are not always axis-aligned.
ITK, ITK-Snap, etc. use ITK-based classes for image storage, although they use VTK for visualization.
ITK started off as VTK, supporting only axis-aligned images, but after a couple of years, support for arbitrary oriented images was added. It was not a one-time, breaking change, but an oriented image data class was added and gradually all filters were updated to work correctly with arbitrarily oriented images. Most of the developments were completed in a short time but then it took a few years to fix up the code a few places where it was not obvious to detect the assumption of axis alignment. Make the same change in VTK is more difficult in the sense that VTK is larger than ITK, but may be easier, as only a small part is related to image processing and visualization.
A similar gradual approach would probably work well for VTK, too.
Oriented image support should be introduced to minimal changes to existing code base.
In the first phase, oriented image support would added to vtkImageData and a limited number of readers/writers, filters, etc.
When a filter that has not declared itself to be able to handle oriented image data, has to return with error when it receives an oriented image data as input.
There are two main design options:
Create a new image data class, such as vtkOrientedImageData and during the transition time use that to distinguish between oriented and axis-aligned data. Once the transition time is over (probably in a few years), the old vtkImageData class can be replaced by vtkOrientedImageData (an alias may be kept for backward compatibility).
Keep using vtkImageData but add direction information to the vtkInformation object that already stores origin and spacing information.
Option A: ITK used this approach, but maybe because ITK objects did not have vtkInformation-type generic description that can be passed through pipelines. The disadvantage of this approach that a new class has to be added, all filters and other classes has to be changed to use this new type, and once the transition is complete (after many years) the vtkOrientedImageData would need to be retired and changed back everywhere to vtkImageData.
Option B would be less invasive, requiring less changes to existing code base, but probably it would require a special check in the pipeline for presence of image rotation (for a couple of years, until VTK fully supports oriented image data).
To be investigated:
How to manage vtkStructuredPoints (which is a specialization of vtkImageData) and vtkUniformGrid (which is an extension of vtkImageData)?
If a new class is created, should it be child or parent of vtkImageData? All ImageData could be interpreted as a vtkOrientedImageData, so it seems to make sense to make vtkOrientedImageData the base class.
Phase 1. Algorithm feasibility - DONE
Implement oriented image class and filters can operate on oriented image data.
This phase is completed. We’ve found that oriented image data support greatly simplified data application code, but it was difficult to use filters that are not orientation-aware. It was not always trivial to update classes - had to spend time with completely understanding how the filter or transform worked (update how normals and derivatives, cutting plane positions and orientations are computed) - to add oriented image data support. However, the filters did not become slower or more complex by making them orientation-aware.
Related code is available here:
Phase 2. Design feasibility - WIP
Simply implementing vtkOrientedImageData class did not help much because developers need to propagate direction information manually. It is important to make the pipeline accept vtkOrientedImageData.
Plan: Implement mechanism that allows storing and passing direction data in vtkImageData. It would be nice to implement a mechanism to report error if a filter receives oriented image data as input but it can only work with axis-aligned image data.
Implementation is work in progress in this branch:
Phase 3. Oriented image data infrastructure integration into VTK
Merging minimal amount of changes to VTK that allows reading, processing, and displaying image data. Code reviews, testing, etc.
Phase 4. Extending oriented image data support
Update more VTK filters, IO classes, etc. to support oriented image data. This phase will likely to be a continuous activity for several years.
What is orientation specifically? Are we talking euler angles (3dof) or something more general that can include shears and inversions, lefthanded coords etc?
euler angles (3dof). These can be described using vtkMatrix3x3
we should treat it the same way as in ITK.
Storing as a 3x3 matrix is the way to go (it would allow fast computation of coordinate values).
We would probably expect right-handed orthonormal basis (orthogonal and vector size of 1) in most filters but in some filters it may be just as easy to accept any 3x3 matrix, so we may decide to enforce constraints at filter level. For example, in itkImageBase.h: "The columns of the Direction matrix are expected to form an orthogonal right handed coordinate system. But this is not checked nor enforced in itk::ImageBase."
What should be done for vtkStructuredPoints and vtkUniformGrid ?
Currently, vtkImageData is used as base class for vtkStructuredPoints and vtkUniformGrid. If we add orientation to vtkImageData, then these classes also get it. It would be probably generally useful, but it has to be checked what is the preference of heavy users of vtkStructuredPoints and vtkUniformGrid classes (do they want to change their base class and not even have the opportunity to orient the grids, or just disable it in those filters that cannot handle it yet).