Transforming Coordinates

[kleingeo]

I am wondering about transforming coordinates. If I have an image and associate coordinate (x, y, z) and I apply a transform to image (say resampling and some affine augmentation), is there a way to transform the coordinate into the new space?

[fepegar]

Hi, @kleingeo. I am modifying the Affine transform so you can get the SimpleITK transform from it. Then, you can use it to transform coordinates as you wish. Note that ITK expects coordinates in LPS orientation.

In [1]: import torchio as tio

In [2]: colin = tio.datasets.Colin27()

In [3]: transform = tio.RandomAffine()

In [4]: transformed = transform(colin)

In [5]: sitk_transform = transformed.history[0].get_affine_transform(transformed.t1)

In [6]: sitk_transform
Out[6]: <SimpleITK.SimpleITK.CompositeTransform; proxy of <Swig Object of type 'itk::simple::CompositeTransform *' at 0x7fe2313bed20> >

This should be available in v0.18.55.

[romainVala]

Hi
Sorry if I misunderstood but :

I have an image and associate coordinate (x, y, z)

not sure to understand here:
usually, with nifti image, one get an image and an affine that will allow to transform index i,j,k to space coordinate x,y,z in mm. After any torchio transform, the transformed image will have an affine which is changed correctly (to reflect the new index to voxel correspondance )...

[kleingeo]

Right now when the image is modified with a random affine transform the resulting image's affine matrix is unchanged.

[fepegar]

@romainVala, I think @kleingeo just means they have some point in the physical space of the image (e.g. bounding box or mesh or landmarks) and they want it to be transformed using the same transform that was used for the image.

@kleingeo, as discussed in #662, the matrix is not expected to change after the transformation, as the image is resampled.

[kleingeo]

That is right, there are landmarks with the image that need to be transformed as the image changes. This will work fine with the Affine transform (backing out the affine transform, taking its inverse and getting the transformed coordinates) but this becomes more complicated with the Elastic transformation since the BSpline cannot be inverted.

[fepegar]

This will work fine with the Affine transform (backing out the affine transform, taking its inverse and getting the transformed coordinates)

Why do you need to invert the transform?

the BSpline cannot be inverted.

It can!

In [1]: import torchio as tio

In [2]: colin = tio.datasets.Colin27()

In [3]: elastic = tio.RandomElasticDeformation()

In [4]: transformed = elastic(colin)

In [5]: tio_transform = transformed.history[0]

In [6]: inverse = tio_transform.inverse()

In [7]: sitk_transform = inverse.get_bspline_transform(transformed.t1.as_sitk())  # note a SimpleITK image is needed here

I will commit something in a second to make step 7 possible (in the current version, the control points need to be passed).

[romainVala]

Ok, I missed the landmarks objectiv, ... although what is it exactly; coordinate in voxel space or in mm space ?

But I still do not understand why you need an inversion here, what you want is to apply to your landmak the same affine as the RandomAffine or the same deformation field as the RandomElasticDeformation (no inversion needed ?)

but I just see fernando has the same question ...

[kleingeo]

I have both the coordinates in voxel space (but it is easier for me to get them in mm space). I do not need the inversion necessarily, but it is one way to get the transformed coordinates back into voxel space. Essentially, I have a landmark that I need to track after I apply augmentations to the sample data.

[romainVala]

to get any landmark from mm to voxel space, the only thing you need, is the affine from the voxel space of the transformed image. (and it's inverse indeed)
But to get the transformed landmarks (if they are in mm), you need to apply the direct affine, or Elastic deformation
this is what @fepegar suggest with getting the affine or the bspline from the transform's history

The alternative way to do, would be to add a torchio type 'landmarks' and let the torchio transforms take care of it ... kind of related to the surface issue ... (#537)

[kleingeo]

Right now, with SimpleITK, the way it is working for me is as follows. Right now I take the inverse of the transform. They way I understand it (and correct me if I am wrong) is:

T x index = coord -> where T is the transform, index is in voxel space and coord is in physical space
index = T^-1 x coord

compos = sitk.CompositeTransform([trans, rot, scale, translate_output_center])

resample_filter = sitk.ResampleImageFilter()
resample_filter.SetSize(size)
resample_filter.SetInterpolator(sitk.sitkLinear)
resample_filter.SetOutputSpacing(output_spacing)
resample_filter.SetOutputOrigin(image.GetOrigin())
resample_filter.SetOutputDirection(image.GetDirection())
resample_filter.SetTransform(composite_transform)

output_image = resample_filter.Execute(image)

# landmark is a list of coodinates in physical space (ie, [x, y, z] in mm)
new_landmark = output_image.TransformPhysicalPointToContinuousIndex(composite_transform.GetInverse().TransformPoint(landmark))

[romainVala]

yes the index to coord is correct,
but here you transformed the image again, did you with what fernando propose #661 (comment) ?
get the sitk transform from torchio transform's history.

In [5]: sitk_transform = transformed.history[0].get_affine_transform(transformed.t1)

In [6]: new_landmark_mm = sitk_transform(landmark)

does it do the same thing ? or may be you need to transform back to voxel with something like

In [7]: new_landmark_vox = np.linalg.inv(transformed.t1.affine) * new_landmark

[kleingeo]

but here you transformed the image again

What do you mean by transformed the image again? There is one (set) of transforms being applied to the image.

In [6]: new_landmark_mm = sitk_transform(landmark)

This does not work as sitk_transform is not callable.

The method described above does not work either.

# Random affine transform
aff_trans = tio.RandomAffine(
    scales=0.15,
    degrees=0.25 * 180 / np.pi,
    translation=30,
    center='origin',
    default_pad_value='minimum')

subject_trans = aff_trans(subject)

sitk_transform = subject_trans.history[0].get_affine_transform(subject_trans['image'])

# Go from voxel space to image space. Visually confirmed this from ITK-Snap
phys_coord = np.array([-1, -1, 1]) *  np.matmul(subject['image'].affine, list(landmarks_list[7].coords) + [1, ])[:3]

# Transform physical coordinate
new_landmark_mm = sitk_transform.TransformPoint(phys_coord)
new_landmark_vox = np.matmul(np.linalg.inv(subject_trans['image'].affine), np.array([-1, -1, 1, 1]) * np.array(list(new_landmark_mm) + [1, ]))

>>> print(new_landmark_mm)
(-228.81963634350666, 125.73227822861823, -121.10387276726081)
>>> print(new_landmark_vox)
array([ 8.5418, 33.4706, 46.9879,  1.    ])

new_landmark_mm_attempt_2 = sitk_transform.GetInverse().TransformPoint(phys_coord)
new_landmark_vox_attempt_2 = np.matmul(np.linalg.inv(subject_trans['image'].affine), np.array([-1, -1, 1, 1]) * np.array(list(new_landmark_mm_attempt_2) + [1,]))

>>> print(new_landmark_mm_attempt_2)
(-109.44018148213122, 226.55606202734995, -148.57536314161302)
>>> print(new_landmark_vox_attempt_2)
array([ 66.1696, 138.9838,  19.5164,   1.    ])


phys_coord_sitk = subject['image'].as_sitk().TransformContinuousIndexToPhysicalPoint(landmarks_list[7].coords.astype('float'))
new_landmark_mm_sitk = sitk_transform.GetInverse().TransformPoint(phys_coord_sitk)
new_landmark_mm_vox = subject_trans['image'].as_sitk().TransformPhysicalPointToContinuousIndex(new_landmark_mm_sitk)

>>> print(new_landmark_mm_sitk)
(-109.44018148213122, 226.55606202734995, -148.57536314161302)
>>> print(new_landmark_mm_vox)
(66.16960758375194, 138.98379850301927, 19.516378771458776)


image_trans = subject_trans['image']
seg_trans = subject_trans['seg']

image_trans.save('image_trans.nii.gz')
seg_trans.save('seg_trans.nii.gz')

[romainVala]

Hi
@fepegar I think there is a bug in the get_affine_transform : you forgot the translation !

To get a reproducible example we can test, we need to use the same subject and the same landmark
so something like this

import torchio as tio
import numpy as np
import torch 
torch.manual_seed(1)   #so that we get the same random transform
subject = tio.datasets.Colin27() 

# Random affine transform (without translations)
aff_trans = tio.RandomAffine(
    scales=0.15,
    degrees=0.25 * 180 / np.pi,
    translation=0,
    center='origin',
    default_pad_value='minimum')

subject_trans = aff_trans(subject)

sitk_transform = subject_trans.history[0].get_affine_transform(subject_trans['t1'])

center_landmark_vox = [90,108,90,1] 
center_landmark_mm = np.matmul(subject['t1'].affine, center_landmark_vox)[:3]

transformed_landmark_mm =  sitk_transform.TransformPoint(center_landmark_mm)

In [72]: transformed_landmark_mm                                                                                                 
Out[72]: (2.1402958874783193, -13.87748607416791, 20.090379359295813)

transformed_landmark_vox = np.matmul( np.linalg.inv(subject['t1'].affine), list(transformed_landmark_mm) +[1] )

This is look almost fine. I visually test different point, and I do see in the transformed imaged (at the predicted position by transformed_landmark_vox (or mm).
but it seems there is one voxel shift ... (although it is difficult to be 100 % sure)

[fepegar]

@fepegar I think there is a bug in the get_affine_transform : you forgot the translation !

Translation is used here:

https://github.com/fepegar/torchio/blob/97722cb9acefe5d58a19a6271fbec658201b0373/torchio/transforms/augmentation/spatial/random_affine.py#L251

Is it not working?

By the way, you might like apply_affine.

[romainVala]

yes you right (I read the code too quick !) and I may have made a mistake in my first tests ..

it seems indeed to work (for me) with translation too

[romainVala]

yes apply_affine may be easier to use, but I thing it is important to understand what need to be done
here, just a matrix multiplication : so simple, but also so many way to make it wrong ... (because of convention ect ...)

@kleingeo when you write :

# Go from voxel space to image space. Visually confirmed this from ITK-Snap
phys_coord = np.array([-1, -1, 1]) *  np.matmul(subject['image'].affine, list(landmarks_list[7].coords) + [1, ])[:3]

this is kind of strange, and the multiplication with np.array([-1, -1, 1]) should not be needeed !
your landmarks_list[7].coords) are all positive ?

Not sure why you need it, but it may be related to negativ number in the diagonal of your image's affine, ...
I do not know if itk-snap handle it properly (may be try an other viewer Slicer (or mrview from mrtrix)

to investigate and make things easier, you can reoder the data in order to have only positive number on the affine
you can do it with torchio
image_ok = tio.ToCanonical(image)

[romainVala]

just to better understand : what is the affine matrix of your images ?
to change it just try
tio.ToCanonical

[fepegar]

@romainVala, the physical space of ITK is, like DICOM, LPS+. Libraries like NiBabel use RAS+. That's why one needs to be careful with the first two axes.

• https://www.slicer.org/wiki/Coordinate_systems#Anatomical_coordinate_system
• https://nipy.org/nibabel/coordinate_systems.html#other-reference-spaces

[romainVala]

I see the point if one need to mixt the two
but within a soft, one should keep a unique choice (RAS like nibabel !) and then it does not matter...
The the data can be in RAS or LPS, this should only change the affine. no ?
by the way so what the choice within torchio ?

[romainVala]

@fepegar
do I miss something ?; don't you think we need to have a single choice within torchio for either ras or lps ?

[kleingeo]

It is confusing, but there are differences between them (in fact, this is one of the most annoying aspects of SimpleITK). Here are the affines of the images. You can see how the signs differ in the affine matrices due to the different conventions.

image = sitk.ReadImage('image.nii.gz')

sitk_affine = np.eye(4)
sitk_affine[:3, :3] = np.reshape(image.GetDirection(), (3, 3)) * np.array(image.GetSpacing())
sitk_affine[:3, -1] = np.array(image.GetOrigin())

image_tio = tio.ScalarImage('image.nii.gz')

>>> print(sitk_affine)
array([[ 1.9984e+00,  3.9956e-02,  2.2820e-08, -2.4723e+02],
       [-7.9911e-02,  9.9920e-01,  2.1926e-08,  9.2971e+01],
       [-4.3851e-08, -2.2820e-08,  1.0000e+00, -1.6809e+02],
       [ 0.0000e+00,  0.0000e+00,  0.0000e+00,  1.0000e+00]])

>>> print(image_tio.affine)
array([[-1.9984e+00, -3.9956e-02, -2.2820e-08,  2.4723e+02],
       [ 7.9911e-02, -9.9920e-01, -2.1926e-08, -9.2971e+01],
       [-4.3851e-08, -2.2820e-08,  1.0000e+00, -1.6809e+02],
       [ 0.0000e+00,  0.0000e+00,  0.0000e+00,  1.0000e+00]])

[kleingeo]

I also posted in the SimpleITK forum as to you need the transform.GetInverse().TransformPoint() as opposed to transform.TransformPoint().

Basically, the transform points from the output to the input (not sure why...) which is why you need the inverse of the transform to go from input to out.

https://discourse.itk.org/t/transform-point-from-affine-transform/4429

[romainVala]

Ok, I'll have to check again if the GetInverse make more sense

For the affine, yes there is a signs difference, As explained #661 (reply in thread). This is expected, as the internal representation is LPS for itk and RAS for nibabel (and torchio)

In torchio, you can see that the function tio.io._read_sitk is adding an extra step to convert from LPS to RAS, so at the end both torchio fonction tio.io._read_sitk and tio.io._read_nibabel retrun now the same affine

I think that is fine, (and that answer to my question #661 (reply in thread)): torchio internal data representation is RAS, (and the affine is set accordingly)

So to resume
_You should stay also with the RAS affine, and make your vox to mm (or mm to vox) transform with this affine.
(not sure, but if you stay with itksnap as viewer, and if it has the same convention as stik (LPS), then you may see some inconsistency
that is why I would test with an other viewer (like Slider, or mrview (form mrtrix))

_ Now about applying the sitk transform you get from history, it should be fine as it apply in the mm space