From c425a022c310ef3527615f114a49d5e397bfd769 Mon Sep 17 00:00:00 2001 From: hrobarts Date: Tue, 11 Nov 2025 14:52:44 +0000 Subject: [PATCH 1/2] Calculate euler angles for parallel tilted geometry --- Wrappers/Python/cil/plugins/tigre/Geometry.py | 32 +++++++++++++++++-- 1 file changed, 30 insertions(+), 2 deletions(-) diff --git a/Wrappers/Python/cil/plugins/tigre/Geometry.py b/Wrappers/Python/cil/plugins/tigre/Geometry.py index 358b5bca0e..c88769062d 100644 --- a/Wrappers/Python/cil/plugins/tigre/Geometry.py +++ b/Wrappers/Python/cil/plugins/tigre/Geometry.py @@ -18,6 +18,7 @@ from cil.framework.labels import AcquisitionType, AngleUnit import numpy as np +from scipy.spatial.transform import Rotation try: from tigre.utilities.geometry import Geometry @@ -36,6 +37,32 @@ def getTIGREGeometry(ig, ag): if ag.config.angles.angle_unit == AngleUnit.DEGREE: angles *= (np.pi/180.) + if ag.geom_type == 'parallel' and ag.system_description == 'advanced': + ag_cil = ag.copy() + ag_cil.config.system.align_reference_frame('cil') # geometry in CIL frame + untilted_rotation_axis = ag_cil.config.system.rotation_axis.direction + untilted_rotation_axis /= np.linalg.norm(untilted_rotation_axis) + + tilted_rotation_axis = ag.config.system.rotation_axis.direction + tilted_rotation_axis /= np.linalg.norm(tilted_rotation_axis) + + v = np.cross(untilted_rotation_axis, tilted_rotation_axis) + v_norm = np.linalg.norm(v) + c = float(np.dot(untilted_rotation_axis, tilted_rotation_axis)) + + # if tilted and untilted rotation axes are not parallel calculate Euler angles + if v_norm > 1e-12: + theta = np.arctan2(v_norm, c) + rotation_matrix = Rotation.from_rotvec(theta * v / v_norm) + euler_angles = [] + for angle in angles: + R1 = Rotation.from_euler("z", angle, degrees=False) + combined = rotation_matrix * R1 + euler_angles.append(combined.as_euler("ZYZ", degrees=False)) + tg.euler_angles = np.array(euler_angles, dtype=np.float32) + + return tg, tg.euler_angles + #convert CIL to TIGRE angles s angles = -(angles + np.pi/2 +tg.theta ) @@ -89,13 +116,14 @@ def __init__(self, ig, ag): self.mode = 'cone' else: - if ag_in.system_description == 'advanced': - raise NotImplementedError ("CIL cannot use TIGRE to process parallel geometries with tilted axes") self.DSO = clearance_len self.DSD = 2*clearance_len self.mode = 'parallel' + # if ag_in.system_description == 'advanced': + # raise NotImplementedError ("CIL cannot use TIGRE to process parallel geometries with tilted axes") + # number of voxels (vx) self.nVoxel = np.array( [ig.voxel_num_z, ig.voxel_num_y, ig.voxel_num_x] ) # size of each voxel (mm) From 9ba6f5dbd543c2cfeffce59ec63ca63de13a461b Mon Sep 17 00:00:00 2001 From: hrobarts Date: Wed, 21 Jan 2026 11:59:49 +0000 Subject: [PATCH 2/2] Use rotation matrices not scipy euler methods --- Wrappers/Python/cil/plugins/tigre/Geometry.py | 60 +++++++++++++++++-- 1 file changed, 55 insertions(+), 5 deletions(-) diff --git a/Wrappers/Python/cil/plugins/tigre/Geometry.py b/Wrappers/Python/cil/plugins/tigre/Geometry.py index 0fc1ccd66b..9971388549 100644 --- a/Wrappers/Python/cil/plugins/tigre/Geometry.py +++ b/Wrappers/Python/cil/plugins/tigre/Geometry.py @@ -18,7 +18,6 @@ from cil.framework.labels import AcquisitionType, AngleUnit import numpy as np -from scipy.spatial.transform import Rotation try: from tigre.utilities.geometry import Geometry @@ -53,12 +52,12 @@ def getTIGREGeometry(ig, ag): # if tilted and untilted rotation axes are not parallel calculate Euler angles if v_norm > 1e-12: theta = np.arctan2(v_norm, c) - rotation_matrix = Rotation.from_rotvec(theta * v / v_norm) + rotation_matrix = rotation_matrix_from_vec(theta * v / v_norm) euler_angles = [] for angle in angles: - R1 = Rotation.from_euler("z", angle, degrees=False) - combined = rotation_matrix * R1 - euler_angles.append(combined.as_euler("ZYZ", degrees=False)) + R1 = rot_z(angle) + combined = rotation_matrix @ R1 + euler_angles.append(euler_from_matrix_zyz(combined)) tg.euler_angles = np.array(euler_angles, dtype=np.float32) return tg, tg.euler_angles @@ -75,6 +74,57 @@ def getTIGREGeometry(ig, ag): angles[i] = a return tg, angles + +def rotation_matrix_from_vec(rotvec): + rotvec = np.asarray(rotvec, dtype=float) + theta = np.linalg.norm(rotvec) + + diag = np.array([[1., 0., 0.], + [0., 1., 0.], + [0., 0., 1.]]) + + if theta == 0.0: + return diag + + axis = rotvec / theta + x, y, z = axis + K = np.array([[ 0, -z, y], + [ z, 0, -x], + [-y, x, 0] + ]) + + R = (diag + + np.sin(theta) * K + + (1 - np.cos(theta)) * (K @ K) + ) + + return R + +def rot_z(angle): + c, s = np.cos(angle), np.sin(angle) + return np.array([[c, -s, 0], + [s, c, 0], + [0, 0, 1]]) + +def euler_from_matrix_zyz(R, eps=1e-12): + """ + Equivalent to scipy Rotation.as_euler("ZYZ", degrees=False) + """ + R = np.asarray(R, dtype=float) + + beta = np.arccos(np.clip(R[2, 2], -1.0, 1.0)) + + if abs(np.sin(beta)) > eps: + alpha = np.arctan2(R[1, 2], R[0, 2]) + gamma = np.arctan2(R[2, 1], -R[2, 0]) + else: + alpha = np.arctan2(R[0, 1], R[0, 0]) + gamma = 0.0 + + alpha = (alpha + np.pi) % (2*np.pi) - np.pi + gamma = (gamma + np.pi) % (2*np.pi) - np.pi + + return np.array([alpha, beta, gamma]) class TIGREGeometry(Geometry):