Elastix idiosyncracy?

Engineering
1

Hi,
The flat array “direction” property of sitk images are specified as a row-order concatenation of the 3x3 direction matrix.
For Elastix/SimpleElastix, the Direction property specified in Transformparamterer maps is a column order flattening.
This is a little impractical and leaves room for error. I guess it is really a clash of conventions between elastix and ITK.
I just wanted to bring this up - not sure which table (elastix or ITK) it belongs on!
Self contained example below.

Soren

import SimpleITK as sitk
import numpy as np

test_image = sitk.GetImageFromArray(np.zeros((10,10,10),dtype=np.float32))
dir1 = [0,0,-1] #voxel X maps to neg Z world
dir2 = [0,1,0] #voxel Y maps to pos Y world
dir3 = [1,0,0] #voxel Z maps to pos X world

directions_serial = dir1+dir2+dir3
test_image.SetDirection(directions_serial)
sitk.WriteImage(test_image,"test.nii")
test_image_read = sitk.ReadImage("test.nii")
np.testing.assert_equal(np.array(test_image_read.GetDirection()),np.array(test_image.GetDirection()))
np.testing.assert_equal(np.array(test_image_read.GetDirection()),np.array(directions_serial))
#ok as expected. Now lets manually define a resampling space based on test_image

#first - obtain a transformparametermap - this is a litte detaour
elastix = sitk.ElastixImageFilter()
elastix.SetParameterMap(elastix.GetDefaultParameterMap("rigid",1))
elastix.SetMovingImage(test_image)
elastix.SetFixedImage(test_image)
elastix.Execute()
transformparametermap_template = elastix.GetTransformParameterMap()
# ok, we have a parametermap. The directions of this map are identical to the fixed image by design, but the flattening of the direction matrix was by column order
directions_flat_from_transform_file = np.array([float(v) for v in transformparametermap_template[0]["Direction"]])
np.testing.assert_equal(directions_flat_from_transform_file,np.array(directions_serial))
# FAILS - ordering is not identical

#so when interfacting with transformparameter maps, the flattened arrays use column ordering as opposed to row-ordering for the simple itk images

re_ordered_directions_from_transform_file = directions_flat_from_transform_file.reshape((3,3),order="F").flatten(order="C")
np.testing.assert_equal(re_ordered_directions_from_transform_file,np.array(directions_serial))
2

@Niels_Dekker and @blowekamp might pitch in.

3

Is this a row vs column issue or a transform direction issue ( forwards vs reverse ):


In [29]: tx1.SetParameters([0,0,-1,0,1,0,1,0,0, 0,0,0])

In [30]: print(tx1)
itk::simple::AffineTransform
 AffineTransform (0x13467a5e0)
   RTTI typeinfo:   itk::AffineTransform<double, 3u>
   Reference Count: 1
   Modified Time: 1655
   Debug: Off
   Object Name: 
   Observers: 
     none
   Matrix: 
     0 0 -1 
     0 1 0 
     1 0 0 
   Offset: [0, 0, 0]
   Center: [0, 0, 0]
   Translation: [0, 0, 0]
   Inverse: 
     0 0 1 
     0 1 0 
     -1 0 0 
   Singular: 0

Note the inverse of your transform matrix is the transpose.

1 Like
4

Sorry I am not sure I understand - are you saying there is no difference in the format between Elastix and SImpleITK or are you confirming there is a format difference?
As the example shows, the direction vector (flattened direction matrix) of the fixed image is not the same as the vector returned by elastix. The conventions differ and it seems counterintuitive to have 2 different conventions.

I should emphasize I am not looking at anything to do with the actual transformation from elastix. the focus is on the fixed image directions vs the directions of the resampling space (which is the fixed image directions by definition). So not an inverse transform issue imo.

5

One potential source of confusion is that in your example the image direction was not set correctly here:

You can see the direction vectors if you save the image in .nrrd file format by calling sitk.WriteImage(test_image,"test.nrrd"). The relevant line in the output file is this: space directions: (0,0,1) (0,1,0) (-1,0,0), which shows that the vectors are not what you expected.

In all ITK and VTK related libraries a matrix is stored in memory in C order: [dir1_x, dir2_x, dir3_x, dir1_y, dir2_y, dir3_y, dir1_z, dir2_z, dir3_z], so one way to set the direction values correctly would be this:

directions = np.column_stack([dir1, dir2, dir3])
directions_serial = directions.flatten().tolist()
test_image.SetDirection(directions_serial)

It would be a bit more convenient and less error prone if SimpleITK image offered direction setter methods that take a 3x3 matrix, or 3 column direction vectors.

6

Proposals and suggestions have always been welcome. Please create a SimpleITK issue to discuss further, what would be a better unambiguous and cross language way to set direction matrices.

7

Thank you @blowekamp. Let’s see what @Soren_Christensen says, but if it turns out that the issue was that he used wrong element ordering in directions vector then I’ll submit a feature request to SimpleITK to allow setting orientation using a 3x3 matrix.

8

Thank you for pointing that out Andras - you are correct, SetDirection() expects row major ordering. I have pasted a corrected example below. However, that does not change the core issue - the flattened direction array notation between SimpleITK and Elastix differ.
Here’s a summary. Direction matrix:
[ 0, 0, 1],
[ 0, 1, 0],
[-1, 0, 0]
Putting this in an sitk image:
test_image.SetDirection(direction_matrix.flatten(order=“C”).tolist())
the direction vector is now: array([ 0, 0, 1, 0, 1, 0, -1, 0, 0]), following the row order pointed out by you Andras.
Now, register any image to this image (fixed image) using elastix, and as per default elastix behavior, the image is resampled in the fixed space. Here’s the directions in the resulting transform parameter file:
array([ 0., 0., -1., 0., 1., 0., 1., 0., 0.])
So this means elastix is putting out a column order flattening.

I believe the conventions do differ. I guess it is not a bug as such, but it would probably be less error prone if the conventions were identical.

Corrected example:

import SimpleITK as sitk
import numpy as np

test_image = sitk.GetImageFromArray(np.zeros((10,10,10),dtype=np.float32))


dir1 = np.array([[0,0,-1]]).T #voxel X maps to neg Z world
dir2 = np.array([[0,1,0]]).T #voxel Y maps to pos Y world
dir3 = np.array([[1,0,0]]).T #voxel Z maps to pos X world
mydir_mat = np.concatenate( (dir1,dir2,dir3),axis=1 )

directions_serial = mydir_mat.flatten(order="C")
test_image.SetDirection(directions_serial.tolist())
sitk.WriteImage(test_image,"test.nii")
test_image_read = sitk.ReadImage("test.nii")
np.testing.assert_equal(np.array(test_image_read.GetDirection()),np.array(test_image.GetDirection()))
np.testing.assert_equal(np.array(test_image_read.GetDirection()),np.array(directions_serial))
#ok as expected. Now lets manually define a resampling space based on test_image

#first - obtain a transformparametermap - this is a litte detaour
elastix = sitk.ElastixImageFilter()
elastix.SetParameterMap(elastix.GetDefaultParameterMap("rigid",1))
elastix.SetMovingImage(test_image)
elastix.SetFixedImage(test_image)
elastix.Execute()
transformparametermap_template = elastix.GetTransformParameterMap()
# ok, we have a parametermap. The directions of this map are identical to the fixed image by design, but the flattening of the direction matrix was by column order
directions_flat_from_transform_file = np.array([float(v) for v in transformparametermap_template[0]["Direction"]])
np.testing.assert_equal(directions_flat_from_transform_file,directions_serial)
#FAILS
#x: array([0., 0., -1., 0., 1., 0., 1., 0., 0.])
#y: array([0, 0, 1, 0, 1, 0, -1, 0, 0])

#so when interfacting with transformparameter maps, the flattened arrays use column ordering as opposed to row-ordering for the simple itk images

re_ordered_directions_from_transform_file = directions_flat_from_transform_file.reshape((3,3),order="F").flatten(order="C")
np.testing.assert_equal(re_ordered_directions_from_transform_file,np.array(directions_serial))
9

I would not expect a transformation matrix and a transform parameter map to have similar ordering of items, as they contain different data. Transform parameters may store additional data, such as center of rotation, and may store less data or store it differently than matrix values, for example direction may be stored as a quaternion. It would be error-prone to make users believe that they can simply copy items between these two different lists. It would cause extensive backward incompatibility damage if someone changed order of parameters in any existing transform class.

Nevertheless, it is undeniable that you and many others have made and will make the same mistakes with the current API; and it is great that you took the time to report this problem.

Would you be able to submit pull requests that improve the documentation, examples, maybe add convenience functions for setting parameters or transforms - in a way that you think would have helped you to avoid mistakes?

1 Like
10

Hi Andras,
Let me be more specific:
I am not talking about the spatial transform - just the Direction of the resampling space.
As I am sure you are aware, the key components of a resampling operation is the spatial transform AND the definition of the target grid. The target grid is defined just like an ITK image, by position, spacing and directions. It is this “direction” I am talking about.
So in the SimpleElastix parameter maps there is “Direction” and “TransformParameters” and I am focused on Direction in this topic. For reference, here are my values below.

>>> transformparametermap_template[0]["Direction"]
('0', '0', '-1', '0', '1', '0', '1', '0', '0')
>>> transformparametermap_template[0]["TransformParameters"]
('0', '0', '0', '0', '0', '0')

The Direction is per default the same as the Direction of the fixed image. The resampling space is defined by the fixed image. Sometimes you wish to change the resampling grid and in this situation you want to change the directions prior to resampling.

As for interface changes I think we still have to agree on the issue. As I see it, it would be advantageous so have an overloaded setter accepting a 3x3 matrix type for sitk/ITK images. For the Elastix transformparametermap, I think the situation is more tricky, because it is an Elastix convention clashing with ITK. Ideally it would be changed on the Elastix side I guess, but it would probably require introduction of a replacement variable to not incorrectly interpret older transform parameter files.

11

Trying again to loop in @Niels_Dekker :smiley:

1 Like
12

For the general array and matrix implementation there two options to consider:

  • How is the matrix stored: column or row major
  • How is the matrix indexed: (i,j,k) vs (k,j,i) or [x][y][z] vs [z][y][z]

This is simply to say that a converting from a generic matrix is still an error prone operations, if the concepts and formats are not understood.

That being said the following may be a useful and clear methods to setting the Direction:

  void sitk::Image SetDirection(const std::vector<double> & column1,
               const std::vector<double> & column2,
               const std::vector<double> & column3 = {},
               const std::vector<double> & column4 = {},
               const std::vector<double> & column5 = {});

EDIT: Not sure if the variable name should start with 1 or 0 :slight_smile:

1 Like
13

Yes I see your point and your suggestion seems to lessen the scope for error. Why do you need more than 3 columns though? Wouldn’t there be max 3 spatial dimension?

Just to keep things on track - I think your suggestion is an improvement, but personally I have not been bothered by the signature of the SetDirection function (the documentation highlights the row order convention). The bigger issue to me is still the convention mismatch between simpleitk image direction and simpleelastix image direction. Making the change above will make ITK simpler to use without consulting the docs, but it won’t solve the issue on the Elastix side.

Thanks,
Soren

14

Unfortunately, Elastix developers went against conventions of the ITK and VTK community when they specified the FlatDirectionCosinesType type. This was probably just an oversight, and should be fixed, before it causes more confusion or - even worse - spreads to more places.

However, this type is used in so many places in Elastix and even used in many transformation files, so it would need to be phased out slowly (a new type and corresponding methods should be introduced, the old type and methods should be deprecated now and then removed in 5-10 years). Also, these kind of changes (API change without any functional benefit, such as bugfixes or new features) are also quite annoying for all existing Elastix users, so preferably the change should be done when there are inevitable other changes around that code anyway. Such inevitable other changes may not come up for several years.

Overall, I feel that just carefully documenting in Elastix this inconsistent type and its every usages may be a better solution in the short term (I’m quite sure that Elastix developers would welcome a pull request from you that improves their documentation). In the long term, I hope Elastix developers will keep in mind this inconsistency and if the opportunity comes then they will improve the code (but who knows, maybe in a couple of years AI will take over image registration the same way as it did with image segmentation and then everybody will just forget about these “classic” registration toolkits).

15

I’m currently looking at this issue, thanks for bringing it up, @Soren_Christensen I certainly understand that the ordering of the Direction elements as assumed by elastix deserves some extra documentation!

Update: please check: DOC: Add a note on the order of matrix elements of `Direction` parameter by N-Dekker · Pull Request #1208 · SuperElastix/elastix · GitHub

16

Hi Andras and Lasso,
Thanks for your comments and considerations! I agree documentation seems like a sufficient remedy given the widespread use of this in Elastix.
Let me know if I can help with anything!
Soren

17

I think I agree that it would have been more intuitive if the elastix Direction parameter would have the matrix elements in row major order, but it appears that it has been column order for more than 16 years already. (So it’s hard to change it now, without breaking backward compatibility.)

Can you please confirm that my pull request is OK? DOC: Add a note on the order of matrix elements of `Direction` parameter by N-Dekker · Pull Request #1208 · SuperElastix/elastix · GitHub

Also please let me know if there are other places in elastix documentation that should be further clarified, with respect to matrix element orderings.

1 Like
18

Thanks!
Just commented there. That’s the only documentation issue I am aware of.
Cheers,
Soren

1 Like
19

For the record, MetaIO also uses column-major order: MetaIO image rotation conventions

20

For your information, I just realized that elastix also assumes column-major order for the bspline transform parameter “GridDirection”. So I extended the elastix pull request to also add notes for the documentation of “GridDirection”: