I did not understand @dgobbi’s comment about either. We will store image origin, spacing, and axis directions as 4x4 IJK to World matrix, the same way as it is done in ITK and many other places in VTK.
Thanks, @estan. Your understanding is correct. I do not want the Origin (or the Spacing) to be incorporated into the 4x4 matrix. I want an additional shift to be incorporated into the 4x4 matrix.
In order to keep this from being a breaking change, we need to have three coordinate systems associated with the images:
- the IJK coords, i.e. what we refer to as Structured Coords in VTK
- the image data coords as we currently understand them, incorporating the Origin and Spacing
- new oriented data coords, which are the current coords multiplied by a new 4x4 matrix
I’ve already given reasons for this approach in my earlier comments, but I can try to give it another go.
I’m thinking we would support methods such as
IndexToModel(index, pos)
ModelToIndex(pos, index)
GetModelToIndexMatrix(double 16)
GetIndexToModelMatrix(double 16)
SetIndexToModelMatrix(double 16)
along with
GetScale
SetScale
GetRotation
SetRotation
etc
So that we can work with the full matrix, or the individual components of it (such as scale, rotation, etc).
Then we have origin and position. David I want to make sure I understand what you are saying here. Normally in a vtkProp3D Origin is a center of rotation/scaling and Position is a final translation. Currently in ImageData though we use Origin as a Position. What you are suggesting is I think
What we currently have
ijk --> * spacing --> + origin
but moving to
ijk --> * spacing --> + origin --> rotation --> + position
in contrast to vtkProp3D which does
data -> - origin --> *scale --> rotation --> + origin --> +position
is that correct?
The API could remain mostly the same, we just need to add a few extra Get/Set methods:
- Origin(double[3])
- Spacing(double[3])
- Directions(vtkMatrix3x3*)
- IndexToVolumeMatrix(vtkMatrix4x4*)
We might add a few more for convenience/performance, something like these:
- VolumeToIndexMatrix(vtkMatrix4x4*)
- VolumeToIndexPosition(double[3], double[3])
- IndexToVolumePosition(double[3], double[3])
“Rotation” does not sound relevant to me. We don’t rotate anything but we directly define the a coordinate system axis directions.
Calling the volume’s physical coordinate system “Model” does not seem very intuitive to me. “Volume” may be a bit better, but there could be better ones.
a) body/space
b) material/spatial
c) Lagrangian/Eulerian
frames are the terms I’m familiar with.
With all these methods and attributes, I don’t understand the resistance to introducing a vtkCoordinateSystem class.
Yes, that is what I am proposing. Take our current image coords (which use Spacing and Origin) and then add on a rotation and position shift (and I’m proposing that the rotation and position shift be stored in a 4x4 matrix). Therefore, it would be a simple matrix multiplication to go from our current implementation to the new implementation. Also, I imagine that eventually “Origin” will fall into disuse for images, i.e. people will start setting it to zero and preferentially use the “Position” instead.
The situation I want to avoid is what Andras is proposing (from my understanding):
ijk --> * spacing --> rotation --> + origin
See, since the ‘rotation’ is inserted between the ‘spacing’ and ‘origin’ it is not straightforward to go from our current image coords to the new image coords via a simple matrix multiplication. In my proposal, I very specifically chose the new transformation parameters to be applied after any existing transformation parameters. So for e.g. volume mapping, marching cubes, etc. we can simply take the current results and apply a 4x4 matrix to get the new results.
Andras, my main contention is that your proposal changes the definition of “Origin” by moving it from one frame to another. I’d rather introduce a new “Position” with a clear definition and let “Origin” gradually go out of style.
Just wanted to mention that I am greatly in support of the idea of adding orientation to vtkImageData. We encounter this situation frequently in radiation therapy. My preference would be to have the solution integrated into vtkImageData.
I don’t see how this would help, since:
- Filters that accept oriented images would need to be updated anyway and they would probably simply use the IndexToSpace matrix (they would not rely on origin or position vectors).
- Filters that do not accept oriented images would work the same way for non-oriented images (since direction matrix is identity, it does not matter if applied before or after translation); and they would throw an error if they receive oriented images.
Having separate “origin” and “position” vectors would introduce some ambiguity. For example, how would you set origin and position vectors when vtkImageData::SetIndexToSpace(vtkMatrix4x4*) is called?
From my understanding, ITK and VTK share the same definition of origin.
VTK uses extent, and ITK uses indices but they both refer to the integer type used as an index in the data container.
the image data coords as we currently understand them, incorporating the Origin and Spacing
The world coordinates from image data are right now incomplete, why is it needed an extra unoriented coordinate? We should only care about two systems: indices space, associated with the data container itself, and knowing nothing about metadata, and world coordinates.
The transformation from index to world might depend on many many parameters, i.e. metadata of the image, in VTK, right now, it uses origin and spacing/scaling. ITK adds a direction/rotation. We can imagine as many transforms we want, but in the end, a 3x3 matrix + translation, or a 4x4 matrix in homogeneous coordinates will encode the information of any concatenation of transforms. As said by others in this thread, algorithms make assumptions when working in the data/indices space (like orthogonal axis), so support for every transformation out there is tricky.
The ITK approach is exactly the same as the VTK approach, but adding a matrix to take into account the rotation of axis. Not sure about the value of having a non-world, non-data coordinate system, and having to maintain an extra variable (that extra shift) to reconcile between systems.
2 Likes
(1) data -> *scale -> rotation -> +origin
versus
(2) data -> *scale -> +origin -> rotation -> +position
My gut is that (1) will make sense to people, and as David suggested with (2) the goal would be over time for origin to fall into disuse resulting basically in (1) but with position instead of origin. So we all want to eventually land at (1) give or take naming.
I’m not sure I understand how (2) would help in terms of transition. Any algorithm dealing with world/model coords would need to be updated and I think the algorithm changes to accommodate (1) would be just as easy as the changes to accommodate (2) unless I’m missing something (which I may be). My temptation would be to try (1) and see if something hits the fan so to speak.
5 Likes
The true solution is to use ITK where interpolation and any other tranformation is in world space ultimately, even for regular image data! So trully use ITK for image-processing and VTK for visualisation 
VTK is simply ill suited 
I usually agree 100% with any comment from you @lassoan, this one is no exception, but let me add something here which caused quite grief for me 
The IJK-2-XYZ transformation should keep native scanner space as a definition for world space instead of LPS or RAS. Those should be optained further down with extra transforms, often just flips!
Because, currently both ITK and VTK IO classes are useless since both define further tranformation to the one coming from the image file format header info. So in my cause where I want to stay in native scanner space is not really possible
Add extra poly method that can supply those ‘stadnard’ LPS/RAS space optionally only !
ITK’s transformation infrastructure is very limited (there is no support for dynamic non-linear transformation chains with arbitrary combination and real-time inverse computation), so it definitely cannot replace VTK’s resampling filters.
I don’t know what you mean by native scanner space. DICOM uses LPS for defining image position and orientation. Probably each device has some internal coordinate systems but normally you don’t have access to that from outside.
The main goal of the activity described above is that you can we should be able to remain in patient space. Eventually, all VTK filters will be updated to work well with arbitrarily oriented images (or they will report that they require axis-aligned image input).
The native header info 1-to-1 (the voxel to world matrix) without any further transformation, which usually correspond to a space defined by manufacturer. It is always patient space, but what you refer is a standart patient space LPS/RAS
To be more precise, looking at vtkNIFTIImageReader.cxx line#745 you can see that the raw SFormMatrix of a nifti file is what i want but it is that with further tranformation what I get when I call the method: GetSFormMatrix
I see now it is due to some further misconception from the nifti standard and not trying to put it in LPS/RAS space , so maybe my bad, but still my point is valid that I cannot get the raw voxel to world-scanner space
Really not sure, what you mean here, ITK support pretty complex non-liniar transforms, also upd/down sampling is done on a lattice mesh in world space which has many benfits compared to image-space, interpolators etc.
We have been trying various registration methods and shear for some reason takes the dimensionality away from the image. We applied affine transformation and would move an image from say origin to Orgin+10, then the afire transformation would apply shear such that the dimensional accuracy of the volume also got changed.
ITK does not come close to what VTK can do in some aspects that are essential interactive visualization.
For example, ITK cuts off displacement to 0 outside the boundary, which prevents inverse transform computation near the boundary. VTK provides meaningful displacements in the entire spatial domain. ITK can only invert a displacement field non-linear transform and it may take minutes. If you want to transform an image and mesh using a transform that you adjust interactively then it’s just not feasible with ITK. In contrast, VTK can compute inverse of a transform in real-time (and any combination of any transforms, even if they are in a hierarchical pipeline).
Overall, if you just do bulk computation then ITK’s transform infrastructure is sufficient but for real-time interactive visualization involving non-linear transforms, then you need VTK.
ohhh like that, indeed you are right again, I was only talking about data processing, image data processing!
Correct me if I am wrong but no reason to stick to VTK for that !?
For non-interactive, performance-insensitive applications ITK’s transform infrastructure is mostly sufficient. Lack of transformation domain boundary handling can cause problems in inverse computation (which may prevent using a transform for warping both images and points and meshes) but I guess it is always possible to find workarounds.
What is the status on this proposal @ken-martin? Will it ever be implemented or can I assume that an image data is always axis aligned? I’m currently writing a filter for image data that would need to be written differently if axis alignment is not warranted.
If we decide to go forward, would it make sense to start to implement the getters you proposed @ken-martin so, even if we cannot currently change the image alignment, we can write filters that will take it into account the day it would be implemented?