rotation of both image and landmarks issue

[Jakub_Mitura]

Hello I try to rotate an image and landmarks , around image center. I tried resampling landmarks boh by resample physical point and manually using rotation matrix (I would prefer to stic with second solution as I need to make it differentiable in the end). However it do not work Image get rotated but fiducial points are not staying in the same anatomical positions (they get generally get outside of the body after rotation).
code below, I attach also image and landmark points

import SimpleITK as sitk
import numpy as np
from jax.scipy.spatial.transform import Rotation

def save_landmarks_for_slicer(points, out_file_path,original_json_path):
    with open(original_json_path) as f:
        data = json.load(f)
    for i in range(len(cps)):
        data['markups'][0]['controlPoints'][i]['position']=points[i]

    with open(out_file_path, 'w') as f:
        json.dump(data, f)

landmarks=np.load('/workspaces/pilot_lymphoma/data/pat_3/lin_transf/From.npy')
image=sitk.ReadImage('/workspaces/pilot_lymphoma/data/pat_3/lin_transf/study_0_SUVS.nii.gz')
# Assume 'image' is your image and 'landmarks' is your list of landmark points
# weights=jnp.ones(6)*30
weights=np.ones(6)*30
Rotationn=Rotation.from_euler('xyz', weights[0:3], degrees=True)
quaternion=np.array(Rotationn.as_quat())

# Get the center of the image
center = image.TransformContinuousIndexToPhysicalPoint([(index - 1) / 2.0 for index in image.GetSize()])

# Create a VersorRigid3DTransform
transform = sitk.VersorRigid3DTransform()
transform.SetCenter(center)
transform.SetRotation([float(quaternion[0]),float(quaternion[1]),float(quaternion[2]),float(quaternion[3])])

# Resample the image
resampler = sitk.ResampleImageFilter()
resampler.SetOutputDirection(image.GetDirection())
resampler.SetOutputOrigin(image.GetOrigin())
resampler.SetOutputSpacing(image.GetSpacing())
resampler.SetSize(image.GetSize())
resampler.SetTransform(transform)
resampled_image = resampler.Execute(image)

sitk.WriteImage(resampled_image, "/workspaces/pilot_lymphoma/data/transformed_image.nii.gz")

rotation_matrix =np.array(Rotationn.as_matrix())
resampled_landmarks= (landmarks - center) @ rotation_matrix.T + center 


save_landmarks_for_slicer(resampled_landmarks, "/workspaces/pilot_lymphoma/data/From_transformed.mrk.json",'/workspaces/pilot_lymphoma/data/pat_3/lin_transf/From.mrk.json')

From.npy (272 Bytes)
From.mrk.json (5.1 KB)

study_0_SUVS.nii.gz

[dzenanz (Dženan Zukić)]

Resampling transforms coordinates from moving image to the fixed image (the opposite of intuitive), so you need to use the inverse of that matrix to transform the landmarks.

[Jakub_Mitura]

Thanks @dzenanz for response Hovewer It did not resolved the issue I modified line

rotation_matrix =np.array(Rotationn.as_matrix())

into

rotation_matrix =np.linalg.inv(np.array(Rotationn.as_matrix()))

[dzenanz (Dženan Zukić)]

@Jakub_Mitura

resampled_landmarks= (landmarks - center) @ rotation_matrix.T + center

You probably also need to subtract/add the center the opposite way, maybe:
resampled_landmarks= (landmarks + center) @ rotation_matrix.T - center

[dzenanz (Dženan Zukić)]

And why fiddle with this manually? Why not:

iTransform = transform.GetInverseTransform()
resampled_landmarks = [iTransform.TransformPoint(p) for p in landmarks]

[Jakub_Mitura]

Thanks for response! Hovewer VersorRigid3DTransform has no inverse but this is not main issue i need the landmarks change to be differentiable (resampling of image do not need to be). I tried changing + to - but still it do not show correctly. I think weather it may be related to non isovolumetric voxels (I suppose it should have not as we are working on word coordinates), or maybe one need to save it diffrently for slicer to display it correctly?

[blowekamp (Bradley Lowekamp)]

@Jakub_Mitura

Hovewer VersorRigid3DTransform has no inverse but this is not main issue i need the landmarks change to be differentiable (resampling of image do not need to be).

You can convert the transform to an AffineTransform before taking the inverse:

affine = sitk.AffineTransform(vr3d.GetMatrix(), vr3d.GetTranslation(), vr3d.GetCenter())

[Jakub_Mitura]

@dzenanz @blowekamp thanks for discussion I still did not found a solution maybe this will give you some idea?
I am adding some visualizations.
I checked scipy rotation matrix and sitk rotation matrix is the same. Images suggests that simple itk transform and manual perform similarly; I modified first post where I write the landmarks into json instead of fcsv as it is theorethically better

1. original landmarks and 3d model of image
   

   

2. inverse sitk affine

affine = sitk.AffineTransform(transform.GetMatrix(), transform.GetTranslation(), transform.GetCenter())
iTransform = affine.GetInverse()
resampled_landmarks = [iTransform.TransformPoint(p) for p in landmarks]

3. Inverse numpy first add center then subtract

rotation_matrix =np.linalg.inv(np.array(Rotationn.as_matrix()))
# rotation_matrix =np.array(Rotationn.as_matrix())
resampled_landmarks= (landmarks + center) @ rotation_matrix.T - center 
resampled_landmarks=np.array(resampled_landmarks)
resampled_landmarks=list(map(tuple,resampled_landmarks))

4. Inverse numpy first subtract center then add

rotation_matrix =np.linalg.inv(np.array(Rotationn.as_matrix()))
# rotation_matrix =np.array(Rotationn.as_matrix())
resampled_landmarks= (landmarks - center) @ rotation_matrix.T + center 
resampled_landmarks=np.array(resampled_landmarks)
resampled_landmarks=list(map(tuple,resampled_landmarks))

5) manually without inverse

rotation_matrix =np.array(Rotationn.as_matrix())
resampled_landmarks= (landmarks - center) @ rotation_matrix.T + center 
resampled_landmarks=np.array(resampled_landmarks)
resampled_landmarks=list(map(tuple,resampled_landmarks))

6. affine no inverse

affine = sitk.AffineTransform(transform.GetMatrix(), transform.GetTranslation(), transform.GetCenter())
# iTransform = affine.GetInverse()
resampled_landmarks = [affine.TransformPoint(p) for p in landmarks]
print(resampled_landmarks)

Conclusion:
Inverse seem to work better as you suggested as one can see in the image below in some axis it rotates correctly

rotation_matrix =np.linalg.inv(np.array(Rotationn.as_matrix()))
resampled_landmarks= (landmarks - center) @ rotation_matrix.T + center 
resampled_landmarks=np.array(resampled_landmarks)
resampled_landmarks=list(map(tuple,resampled_landmarks))

yellow and grey - not rotated; blue and red rotated (3 images of the same case)

but in other axes it do not work well

[blowekamp (Bradley Lowekamp)]

How are you importing and exporting coordinates between SimpleITK and 3D Slicer? There is a transformation that needs to occur there, that may or may not happen depending on how you are doing it. I don’t have a good reference off hand. Maybe @lassoan could provide the reference or the 3DSlicer discourse.

[lassoan (Andras Lasso)]

It seems that you are not using LPS and RAS coordinate systems consistently. It can be confusing because Slicer internally uses RAS coordinate system (for historical reasons), while it uses LPS coordinate system in all files (images, points, transforms, etc. for compatibility with all DICOM-based software). SimpleITK image reader will provide images and transforms in LPS coordinate system. Conversion between RAS and LPS is trivial, you just need to keep track of what coordinate system you are using.

Note that you can very easily apply a transform to the image and landmarks in Slicer: Slicer automatically uses the modeling transform for landmark points and the resampling transform to images. It can invert any transform on-the-fly, even non-linear warping transforns. It also takes care of all necessary LPS/RAS conversions.

[Jakub_Mitura]

Thanks for response @lassoan ! I
There is something in what you say as From2 and From 3 are indicating acetabulum of hip joints in original 2 is on the left and in rotated on the right as one can see below
(blue imae and red landmarks are those rotated ; black image and yellow landmarks are original)

I had ensured now that image is LPS oriented by adding

resampled_image=sitk.DICOMOrient(resampled_image,"LPS")

in Landmarks json it is kept as was

"coordinateSystem": "LPS"

Hovewer the problem persist.