Added manufacturer check: if Varian then its CBCT, else its Siemens or Phillips and its CT.

SeekTransformModule/SeekTransformModule.py

@@ -33,7 +33,7 @@ class SeekTransformModuleWidget(ScriptedLoadableModuleWidget):
 
         # Dropdown menu za izbiro metode
         self.rotationMethodComboBox = qt.QComboBox()
-        self.rotationMethodComboBox.addItems(["Kabsch", "Horn"])
+        self.rotationMethodComboBox.addItems(["Kabsch", "Horn", "Iterative Closest Point (Kabsch)"])
         self.layout.addWidget(self.rotationMethodComboBox)
         
         # Load button
@@ -268,6 +268,63 @@ class MyTransformModuleLogic(ScriptedLoadableModuleLogic):
             
             return R
 
+        def icp_algorithm(moving_points, fixed_points, max_iterations=100, tolerance=1e-5):
+            """
+            Iterative Closest Point (ICP) algorithm to align moving_points to fixed_points.
+            
+            Parameters:
+            moving_points (list or ndarray): List of points to be aligned.
+            fixed_points (list or ndarray): List of reference points.
+            max_iterations (int): Maximum number of iterations.
+            tolerance (float): Convergence tolerance.
+
+            Returns:
+            ndarray: Transformed moving points.
+            ndarray: Optimal rotation matrix.
+            ndarray: Optimal translation vector.
+            """
+            # Convert to numpy arrays
+            moving = np.array(moving_points)
+            fixed = np.array(fixed_points)
+
+            # Initialize transformation
+            R = np.eye(3)  # Identity matrix for rotation
+            t = np.zeros(3)  # Zero vector for translation
+
+            prev_error = np.inf  # Initialize previous error to a large value
+
+            for iteration in range(max_iterations):
+                # Step 1: Find the nearest neighbors (correspondences)
+                distances = np.linalg.norm(moving[:, np.newaxis] - fixed, axis=2)
+                nearest_indices = np.argmin(distances, axis=1)
+                nearest_points = fixed[nearest_indices]
+
+                # Step 2: Compute the optimal rotation and translation
+                R_new = compute_Kabsch_rotation(moving, nearest_points)
+                centroid_moving = np.mean(moving, axis=0)
+                centroid_fixed = np.mean(nearest_points, axis=0)
+                t_new = centroid_fixed - np.dot(R_new, centroid_moving)
+
+                # Step 3: Apply the transformation
+                moving = np.dot(moving, R_new.T) + t_new
+
+                # Update the cumulative transformation
+                R = np.dot(R_new, R)
+                t = np.dot(R_new, t) + t_new
+
+                # Step 4: Check for convergence
+                mean_error = np.mean(np.linalg.norm(moving - nearest_points, axis=1))
+                if np.abs(prev_error - mean_error) < tolerance:
+                    print(f"ICP converged after {iteration + 1} iterations.")
+                    break
+                prev_error = mean_error
+
+            else:
+                print(f"ICP reached maximum iterations ({max_iterations}).")
+
+            return moving, R, t
+
+
         def compute_translation(moving_points, fixed_points, rotation_matrix):
             """
             Computes the translation vector to align moving_points to fixed_points given a rotation matrix.
@@ -326,7 +383,11 @@ class MyTransformModuleLogic(ScriptedLoadableModuleLogic):
             shNode = slicer.vtkMRMLSubjectHierarchyNode.GetSubjectHierarchyNode(slicer.mrmlScene)
             imageItem = shNode.GetItemByDataNode(volumeNode)
             
+               # Pridobi vse atribute za ta element
+            attributeNames = shNode.GetItemAttributeNames(imageItem)
+            
             #modality = shNode.GetItemAttribute(imageItem, 'DICOM.Modality')
+            #manufacturer = shNode.GetItemAttribute(imageItem, 'DICOM.Manufacturer')
             #print(modality)
             
             # Check if the volume is loaded into the scene
@@ -334,11 +395,13 @@ class MyTransformModuleLogic(ScriptedLoadableModuleLogic):
                 print(f"Volume {volumeName} not present in the scene.")
                 continue
             
+            manufacturer = shNode.GetItemAttribute(imageItem, 'DICOM.Manufacturer')
+            
             # Determine scan type
-            if "cbct" in volumeName.lower():
+            if "varian" in manufacturer.lower():
                 cbct_list.append(volumeName)
                 scan_type = "CBCT"
-            else:
+            else: #Siemens or phillips imamo CTje
                 ct_list.append(volumeName)
                 scan_type = "CT"
             
@@ -388,25 +451,40 @@ class MyTransformModuleLogic(ScriptedLoadableModuleLogic):
                     # print("Rotation Matrix:\n", quaternion_rotation_matrix)
                     # print("Translation Vector:\n", quaternion_translation_vector)
 
+                    # Primer neporavnanih točk
+                    moving_points = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
+                    fixed_points = np.array([[2, 3, 4], [5, 6, 7], [8, 9, 10]])
+
+                    # Uporaba ICP algoritma
+                    aligned_points, R, t = icp_algorithm(moving_points, fixed_points)
+
+                    print("Poravnane točke:")
+                    print(aligned_points)
+                    print("Optimalna rotacijska matrika R:")
+                    print(R)
+                    print("Optimalni translacijski vektor t:")
+                    print(t)
+                    
+                    
                     # Izberi metodo glede na uporabnikov izbor
                     if selectedMethod == "Kabsch":
                         chosen_rotation_matrix = compute_Kabsch_rotation(cbct_points, ct_points)
                         chosen_translation_vector = compute_translation(cbct_points, ct_points, chosen_rotation_matrix)
-                        print("\nKabsch Method:")
-                        print("Rotation Matrix:\n", chosen_rotation_matrix)
-                        print("Translation Vector:\n", chosen_translation_vector)
+                        #print("\nKabsch Method:")
+                        #print("Rotation Matrix:\n", chosen_rotation_matrix)
+                        #print("Translation Vector:\n", chosen_translation_vector)
                     elif selectedMethod == "Horn":
                         chosen_rotation_matrix = compute_Horn_rotation(cbct_points, ct_points)
                         chosen_translation_vector = compute_translation(cbct_points, ct_points, chosen_rotation_matrix)
-                        print("\nHorn Method:")
-                        print("Rotation Matrix:\n", chosen_rotation_matrix)
-                        print("Translation Vector:\n", chosen_translation_vector)
-                    # elif selectedMethod == "Quaternion":
-                    #     chosen_rotation_matrix = compute_quaternion_rotation(cbct_points, ct_points)
-                    #     chosen_translation_vector = compute_translation(cbct_points, ct_points, chosen_rotation_matrix)
-                    #     print("\nQuaternion Method:")
-                    #     print("Rotation Matrix:\n", chosen_rotation_matrix)
-                    #     print("Translation Vector:\n", chosen_translation_vector)
+                        #print("\nHorn Method:")
+                        #print("Rotation Matrix:\n", chosen_rotation_matrix)
+                        #print("Translation Vector:\n", chosen_translation_vector)
+                    elif selectedMethod == "Iterative Closest Point (Kabsch)":
+                        new_points, chosen_rotation_matrix, chosen_translation_vector = icp_algorithm(cbct_points, ct_points)
+                        #chosen_translation_vector = compute_translation(cbct_points, ct_points, chosen_rotation_matrix)
+                        #print("\Iterative Closest Point Method:")
+                        #print("Rotation Matrix:\n", chosen_rotation_matrix)
+                        #print("Translation Vector:\n", chosen_translation_vector)
 
                     imeTransformNoda = cbct_volume_name + " Transform"
                     # Ustvari vtkTransformNode in ga poveži z CBCT volumenom
@@ -428,37 +506,4 @@ class MyTransformModuleLogic(ScriptedLoadableModuleLogic):
                     #transformed_cbct_image = create_vtk_transform(cbct_image_sitk, chosen_rotation_matrix, chosen_translation_vector)
 
         else:
-            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
+            print("CBCT ali CT volumen ni bil najden.")