Smol fix for bigger patients

FirstTryLinear.txt

@@ -6,6 +6,10 @@ from scipy.spatial.transform import Rotation as R
 import slicer
 from DICOMLib import DICOMUtils
 from collections import deque
+import SimpleITK as sitk
+import sitkUtils
+
+#exec(open("C:/Users/lkomar/Documents/Prostata/FirstTryRegister.py").read())
 
 # Define a threshold for grouping nearby points (in voxel space)
 #distance_threshold = 4  # This can be adjusted based on your dataset
@@ -131,7 +135,19 @@ def detect_points_region_growing(volume_name, intensity_threshold=3000, x_min=90
 
     return unique_centroids
 
-def compute_rigid_transform(moving_points, fixed_points):
+
+
+def compute_Kabsch_rotation(moving_points, fixed_points):
+    """
+    Computes the optimal rotation matrix to align moving_points to fixed_points.
+    
+    Parameters:
+    moving_points (list or ndarray): List of points to be rotated CBCT
+    fixed_points (list or ndarray): List of reference points CT
+
+    Returns:
+    ndarray: Optimal rotation matrix.
+    """
     assert len(moving_points) == len(fixed_points), "Point lists must be the same length."
 
     # Convert to numpy arrays
@@ -151,23 +167,164 @@ def compute_rigid_transform(moving_points, fixed_points):
 
     # SVD decomposition
     U, _, Vt = np.linalg.svd(H)
-    R_optimal = np.dot(Vt.T, U.T)
+    Rotate_optimal = np.dot(Vt.T, U.T)
 
     # Correct improper rotation (reflection)
-    if np.linalg.det(R_optimal) < 0:
+    if np.linalg.det(Rotate_optimal) < 0:
+        Vt[-1, :] *= -1
+        Rotate_optimal = np.dot(Vt.T, U.T)
+
+    return Rotate_optimal
+
+def compute_Horn_rotation(moving_points, fixed_points):
+    """
+    Computes the optimal rotation matrix using Horn's method.
+
+    Parameters:
+    moving_points (list or ndarray): List of points to be rotated, CBCT
+    fixed_points (list or ndarray): List of reference points, CT
+
+    Returns:
+    ndarray: Optimal rotation matrix.
+    """
+    assert len(moving_points) == len(fixed_points), "Point lists must be the same length."
+    
+    moving = np.array(moving_points)
+    fixed = np.array(fixed_points)
+    
+    # Compute centroids
+    centroid_moving = np.mean(moving, axis=0)
+    centroid_fixed = np.mean(fixed, axis=0)
+    
+    # Center the points
+    moving_centered = moving - centroid_moving
+    fixed_centered = fixed - centroid_fixed
+    
+    # Compute cross-dispersion matrix
+    H = np.dot(moving_centered.T, fixed_centered)
+    
+    # Compute SVD of H
+    U, _, Vt = np.linalg.svd(H)
+    
+    # Compute rotation matrix
+    R = np.dot(Vt.T, U.T)
+    
+    # Ensure a proper rotation (avoid reflection)
+    if np.linalg.det(R) < 0:
         Vt[-1, :] *= -1
-        R_optimal = np.dot(Vt.T, U.T)
+        R = np.dot(Vt.T, U.T)
+    
+    return R
+
+def compute_quaternion_rotation(moving_points, fixed_points):
+    """
+    Computes the optimal rotation matrix using quaternions.
+
+    Parameters:
+    moving_points (list or ndarray): List of points to be rotated.
+    fixed_points (list or ndarray): List of reference points.
+
+    Returns:
+    ndarray: Optimal rotation matrix.
+    """
+    assert len(moving_points) == len(fixed_points), "Point lists must be the same length."
+    
+    moving = np.array(moving_points)
+    fixed = np.array(fixed_points)
+    
+    # Compute centroids
+    centroid_moving = np.mean(moving, axis=0)
+    centroid_fixed = np.mean(fixed, axis=0)
+    
+    # Center the points
+    moving_centered = moving - centroid_moving
+    fixed_centered = fixed - centroid_fixed
+    
+    # Construct the cross-dispersion matrix
+    M = np.dot(moving_centered.T, fixed_centered)
+    
+    # Construct the N matrix for quaternion solution
+    A = M - M.T
+    delta = np.array([A[1, 2], A[2, 0], A[0, 1]])
+    trace = np.trace(M)
+    
+    N = np.zeros((4, 4))
+    N[0, 0] = trace
+    N[1:, 0] = delta
+    N[0, 1:] = delta
+    N[1:, 1:] = M + M.T - np.eye(3) * trace
+    
+    # Compute the eigenvector corresponding to the maximum eigenvalue
+    eigvals, eigvecs = np.linalg.eigh(N)
+    q_optimal = eigvecs[:, np.argmax(eigvals)]  # Optimal quaternion
+    
+    # Convert quaternion to rotation matrix
+    w, x, y, z = q_optimal
+    R = np.array([
+        [1 - 2*(y**2 + z**2), 2*(x*y - z*w), 2*(x*z + y*w)],
+        [2*(x*y + z*w), 1 - 2*(x**2 + z**2), 2*(y*z - x*w)],
+        [2*(x*z - y*w), 2*(y*z + x*w), 1 - 2*(x**2 + y**2)]
+    ])
+    
+    return R
+
+def compute_translation(moving_points, fixed_points, rotation_matrix):
+    """
+    Computes the translation vector to align moving_points to fixed_points given a rotation matrix.
+    
+    Parameters:
+    moving_points (list or ndarray): List of points to be translated.
+    fixed_points (list or ndarray): List of reference points.
+    rotation_matrix (ndarray): Rotation matrix.
+
+    Returns:
+    ndarray: Translation vector.
+    """
+    # Convert to numpy arrays
+    moving = np.array(moving_points)
+    fixed = np.array(fixed_points)
+
+    # Compute centroids
+    centroid_moving = np.mean(moving, axis=0)
+    centroid_fixed = np.mean(fixed, axis=0)
 
     # Compute translation
-    translation = centroid_fixed - np.dot(centroid_moving, R_optimal)
+    translation = centroid_fixed - np.dot(centroid_moving, rotation_matrix)
+
+    return translation
+
+# def apply_transform(points, rotation_matrix, translation_vector):
+#     points = np.array(points)
+#     transformed_points = np.dot(points, rotation_matrix.T) + translation_vector
+#     return transformed_points
 
-    return R_optimal, translation
+def apply_transform(volume, rotation_matrix, translation_vector):
+    """
+    Transforms a 3D volume using a given rotation matrix and translation vector.
 
-def apply_transform(points, rotation_matrix, translation_vector):
-    points = np.array(points)
-    transformed_points = np.dot(points, rotation_matrix.T) + translation_vector
-    return transformed_points
+    Parameters:
+    volume (sitk.Image): Input 3D volume to be transformed.
+    rotation_matrix (ndarray): 3x3 rotation matrix.
+    translation_vector (ndarray): 1x3 translation vector.
 
+    Returns:
+    sitk.Image: Transformed 3D volume.
+    """
+    # Create an affine transform
+    transform = sitk.AffineTransform(3)
+    transform.SetMatrix(rotation_matrix.flatten())
+    translation_vector = translation_vector.tolist()
+    transform.SetTranslation(translation_vector)
+
+    # Resample the volume
+    resampler = sitk.ResampleImageFilter()
+    resampler.SetReferenceImage(volume)
+    resampler.SetInterpolator(sitk.sitkLinear)  # Linear interpolation
+    resampler.SetTransform(transform)
+    resampler.SetDefaultPixelValue(0)  # Background value for areas outside the original volume
+
+    transformed_volume = resampler.Execute(volume)
+    return transformed_volume
 
 
 
@@ -210,16 +367,96 @@ for volumeNode in slicer.util.getNodesByClass("vtkMRMLScalarVolumeNode"):
 #     print(f"{scan_type} Volume '{vol_name}': {len(points)} points detected.")
 
 
-cbct_points = [centroid for centroid, _ in volume_points_dict[("CBCT", "CBCT1")]]  # Extract only centroids (coordinates)
-ct_points = [centroid for centroid, _ in volume_points_dict[("CT", "CT")]]  # Extract only centroids (coordinates)
-print("CBCT points: ", np.array(cbct_points))
-# Ensure we have enough points for registration
-if len(cbct_points) >= 3 and len(ct_points) >= 3:
-    rotation_matrix, translation_vector = compute_rigid_transform(cbct_points, ct_points)
-    transformed_cbct_points = apply_transform(cbct_points, rotation_matrix, translation_vector)
+if cbct_list and ct_list:
+    # Izberi prvi CT volumen kot referenco
+    ct_volume_name = ct_list[0]
+    ct_points = [centroid for centroid, _ in volume_points_dict[("CT", ct_volume_name)]]
+
+    if len(ct_points) < 3:
+        print("CT volumen nima dovolj točk za registracijo.")
+    else:
+        print("CT points: ", np.array(ct_points))
+        
+        for cbct_volume_name in cbct_list:
+            # Izvleci točke za trenutni CBCT volumen
+            cbct_points = [centroid for centroid, _ in volume_points_dict[("CBCT", cbct_volume_name)]]
 
-    print("Optimal Rotation Matrix:\n", rotation_matrix)
-    print("Translation Vector:\n", translation_vector)
-    print("Transformed CBCT Points:\n", transformed_cbct_points)
+            print(f"\nProcessing CBCT Volume: {cbct_volume_name}")
+            if len(cbct_points) < 3:
+                print(f"CBCT Volume '{cbct_volume_name}' nima dovolj točk za registracijo.")
+                continue
+
+            #print("CBCT points: ", np.array(cbct_points))
+
+            # Compute rotation matrices using different methods
+            svd_rotation_matrix = compute_Kabsch_rotation(cbct_points, ct_points)
+            horn_rotation_matrix = compute_Horn_rotation(cbct_points, ct_points)
+            quaternion_rotation_matrix = compute_quaternion_rotation(cbct_points, ct_points)
+
+            # Compute translations for each method
+            svd_translation_vector = compute_translation(cbct_points, ct_points, svd_rotation_matrix)
+            horn_translation_vector = compute_translation(cbct_points, ct_points, horn_rotation_matrix)
+            quaternion_translation_vector = compute_translation(cbct_points, ct_points, quaternion_rotation_matrix)
+
+            # Display the results for the current CBCT volume
+            print("\nSVD Method:")
+            print("Rotation Matrix:\n", svd_rotation_matrix)
+            print("Translation Vector:\n", svd_translation_vector)
+
+            print("\nHorn Method:")
+            print("Rotation Matrix:\n", horn_rotation_matrix)
+            print("Translation Vector:\n", horn_translation_vector)
+
+            print("\nQuaternion Method:")
+            print("Rotation Matrix:\n", quaternion_rotation_matrix)
+            print("Translation Vector:\n", quaternion_translation_vector)
+
+            # Apply chosen transformation (select one method manually)
+            chosen_rotation_matrix = svd_rotation_matrix  # Change to horn_rotation_matrix or quaternion_rotation_matrix if needed
+            chosen_translation_vector = svd_translation_vector #also change here
+
+            # Pridobitev CBCT slike kot SimpleITK
+            cbct_volume_node = slicer.util.getNode(cbct_volume_name)
+            cbct_image_sitk = sitkUtils.PullVolumeFromSlicer(cbct_volume_node)
+
+            transformed_cbct_image = apply_transform(cbct_image_sitk, chosen_rotation_matrix, chosen_translation_vector)
+
+            print(f"\nTransformed Volume '{cbct_volume_name}' using SVD Method")
+
+            # Shranimo nazaj v Slicer
+            sitkUtils.PushVolumeToSlicer(transformed_cbct_image, name=f"Transformed_{cbct_volume_name}")
 else:
-    print("Not enough points for registration.")
+    print("CBCT ali CT volumen ni bil najden.")
+
+
+# def compute_rigid_transform(moving_points, fixed_points):
+#     assert len(moving_points) == len(fixed_points), "Point lists must be the same length."
+
+#     # Convert to numpy arrays
+#     moving = np.array(moving_points)
+#     fixed = np.array(fixed_points)
+
+#     # Compute centroids
+#     centroid_moving = np.mean(moving, axis=0)
+#     centroid_fixed = np.mean(fixed, axis=0)
+
+#     # Center the points
+#     moving_centered = moving - centroid_moving
+#     fixed_centered = fixed - centroid_fixed
+
+#     # Compute covariance matrix
+#     H = np.dot(moving_centered.T, fixed_centered)
+
+#     # SVD decomposition
+#     U, _, Vt = np.linalg.svd(H)
+#     Rotate_optimal = np.dot(Vt.T, U.T)
+
+#     # Correct improper rotation (reflection)
+#     if np.linalg.det(Rotate_optimal) < 0:
+#         Vt[-1, :] *= -1
+#         Rotate_optimal = np.dot(Vt.T, U.T)
+
+#     # Compute translation
+#     translation = centroid_fixed - np.dot(centroid_moving, Rotate_optimal)
+
+#     return Rotate_optimal, translation