Major improvements and more info.

SeekTransformModule/SeekTransformModule.py

@@ -4,6 +4,7 @@ import scipy
 from scipy.spatial.distance import cdist
 from scipy.spatial.transform import Rotation as R
 import slicer
+import itertools
 from DICOMLib import DICOMUtils
 from collections import deque
 import vtk
@@ -11,6 +12,7 @@ from slicer.ScriptedLoadableModule import *
 import qt
 import matplotlib.pyplot as plt
 import csv
+import time
 
 #exec(open("C:/Users/lkomar/Documents/Prostata/FirstTryRegister.py").read())
 
@@ -87,6 +89,7 @@ class MyTransformModuleLogic(ScriptedLoadableModuleLogic):
     
     
     def run(self, selectedMethod, applyRotation, applyTranslation, applyScaling, writefilecheck):
+        start_time = time.time()
         print("Calculating...")
         
         def group_points(points, threshold):
@@ -329,6 +332,21 @@ class MyTransformModuleLogic(ScriptedLoadableModuleLogic):
 
             return moving, R, t
 
+        def match_points(cbct_points, ct_points):
+            """
+            Vrne cbct_points permutirane tako, da so najbolje poravnani s ct_points
+            glede na vsoto evklidskih razdalj.
+            """
+            best_perm = None
+            min_total_dist = float('inf')
+
+            for perm in itertools.permutations(cbct_points):
+                total_dist = sum(np.linalg.norm(np.array(p1) - np.array(p2)) for p1, p2 in zip(perm, ct_points))
+                if total_dist < min_total_dist:
+                    min_total_dist = total_dist
+                    best_perm = perm
+
+            return list(best_perm)
 
         def compute_translation(moving_points, fixed_points, rotation_matrix):
             """
@@ -369,16 +387,16 @@ class MyTransformModuleLogic(ScriptedLoadableModuleLogic):
             for i in range(4):
                 for j in range(4):
                     vtk_matrix.SetElement(i, j, transform_matrix[i, j])
-            #print("Transform matrix:")
-            #for i in range(4):
-            #    print(" ".join(f"{vtk_matrix.GetElement(i, j):.6f}" for j in range(4)))
+            print("Transform matrix:")
+            for i in range(4):
+               print(" ".join(f"{vtk_matrix.GetElement(i, j):.6f}" for j in range(4)))
             # Create vtkTransform and set the matrix
             transform = vtk.vtkTransform()
             transform.SetMatrix(vtk_matrix)
             return transform
         
         
-        def detect_points_region_growing(volume_name, yesCbct, intensity_threshold=3000, x_min=90, x_max=380, y_min=190, y_max=380, z_min=80, z_max=140, max_distance=9, centroid_merge_threshold=5):
+        def detect_points_region_growing(volume_name, yesCbct, create_marker, intensity_threshold=3000, x_min=90, x_max=380, y_min=190, y_max=380, z_min=80, z_max=140, max_distance=9, centroid_merge_threshold=5):
             volume_node = slicer.util.getNode(volume_name)
             if not volume_node:
                 raise RuntimeError(f"Volume {volume_name} not found.")
@@ -439,12 +457,13 @@ class MyTransformModuleLogic(ScriptedLoadableModuleLogic):
             for centroid, intensity in centroids:
                 if not any(np.linalg.norm(centroid - existing_centroid) < centroid_merge_threshold for existing_centroid, _ in unique_centroids):
                     unique_centroids.append((centroid, intensity))
-                    
-            markups_node = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLMarkupsFiducialNode", f"Markers_{volume_name}")
-            for centroid, intensity in unique_centroids:
-                markups_node.AddControlPoint(*centroid)
-                markups_node.SetDisplayVisibility(False)
-                #print(f"Detected Centroid (RAS): {centroid}, Max Intensity: {intensity}")
+            
+            if create_marker:        
+                markups_node = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLMarkupsFiducialNode", f"Markers_{volume_name}")
+                for centroid, intensity in unique_centroids:
+                    markups_node.AddControlPoint(*centroid)
+                    markups_node.SetDisplayVisibility(False)
+                    #print(f"Detected Centroid (RAS): {centroid}, Max Intensity: {intensity}")
 
             return unique_centroids
 
@@ -526,28 +545,83 @@ class MyTransformModuleLogic(ScriptedLoadableModuleLogic):
             table_ras = np.array(ijkToRasMatrix.MultiplyPoint([*table_ijk, 1]))[:3]
 
             # Doda marker za višino mize
-            table_node = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLMarkupsFiducialNode", f"VišinaMize_{ct_volume_name}")
-            table_node.AddControlPoint(table_ras)
-            table_node.SetDisplayVisibility(False)
+            #table_node = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLMarkupsFiducialNode", f"VišinaMize_{ct_volume_name}")
+            #table_node.AddControlPoint(table_ras)
+            #table_node.SetDisplayVisibility(False)
 
             # Izračun višine v mm
-            image_center_y = dims[1] // 2
-            pixel_offset = table_top_y - image_center_y
-            mm_offset = pixel_offset * spacing[1]
+            #image_center_y = dims[1] // 2
+            pixel_offset = table_top_y
+            #mm_offset = pixel_offset * spacing[1]
 
             #print(f"📏 Miza je {abs(mm_offset):.2f} mm {'nižja' if mm_offset > 0 else 'višja'} od središča.")
-            #print(f"📏 Miza je {abs(pixel_offset)} pixel {'nižja' if pixel_offset > 0 else 'višja'} od središča.")
+            #print(f"📏 Miza je {abs(pixel_offset)} pixlov {'nižja' if pixel_offset > 0 else 'višja'} od središča.")
+            #print(f"📏 CBCT table_top_y: {table_top_y}, image_center_y: {image_center_y}, pixel_offset: {pixel_offset}, mm_offset: {mm_offset}")
+            #print(f"🛏 Absolutna višina mize (RAS Y): {table_ras[1]} mm")
             # Shrani v CSV
             if writefilecheck:
                 file_path = os.path.join(os.path.dirname(__file__), "heightdata.csv")
                 with open(file_path, mode='a', newline='') as file:
                     writer = csv.writer(file)
                     modality = "CBCT " if yesCbct else "CT "
-                    writer.writerow([modality, ct_volume_name, f" Upper part of table detected at Z = {mm_offset:.2f} mm, {pixel_offset} pixels"])
+                    writer.writerow([modality, ct_volume_name, f" Upper part of table detected at Z = {table_ras[1]:.2f} mm"])
 
-            return mm_offset, pixel_offset
-        
+            return table_ras[1], pixel_offset
         
+        def align_cbct_to_ct(volumeNode, scan_type, offset, CT_offset=None, CT_spacing=None):
+            """
+            Aligns CBCT volume to CT volume based on height offset.
+
+            Args:
+                volumeNode (vtkMRMLScalarVolumeNode): The volume node to be aligned.
+                scan_type (str): The type of scan ("CT" or "CBCT").
+                offset (float): The height offset of the current volume from the center in mm.
+                CT_offset (float, optional): The height offset of the CT volume from the center. Required for CBCT alignment.
+                CT_spacing (float, optional): The voxel spacing of the CT volume in mm (for scaling the offset).
+
+            Returns:
+                float: The alignment offset applied to the CBCT volume (if applicable).
+            """
+            if scan_type == "CT":
+                CT_offset = offset
+                CT_spacing = volumeNode.GetSpacing()[1]
+                #print(f"CT offset set to: {CT_offset}, CT spacing: {CT_spacing} mm/voxel")
+                return CT_offset, CT_spacing
+            else:
+                if CT_offset is None or CT_spacing is None:
+                    raise ValueError("CT_offset and CT_spacing must be provided to align CBCT to CT.")
+
+                CBCT_offset = offset
+
+                # Razlika v mm brez skaliranja na CBCT_spacing
+                alignment_offset_mm = CT_offset - CBCT_offset
+
+                #print(f"CT offset: {CT_offset}, CBCT offset: {CBCT_offset}")
+                #print(f"CT spacing: {CT_spacing} mm/voxel, CBCT spacing: {volumeNode.GetSpacing()[1]} mm/voxel")
+                #print(f"Aligning CBCT with CT. Offset in mm: {alignment_offset_mm}")
+
+                # Uporabi transformacijo
+                transform = vtk.vtkTransform()
+                transform.Translate(0, alignment_offset_mm, 0)
+
+                transformNode = slicer.vtkMRMLTransformNode()
+                slicer.mrmlScene.AddNode(transformNode)
+                transformNode.SetAndObserveTransformToParent(transform)
+                volumeNode.SetAndObserveTransformNodeID(transformNode.GetID())
+                slicer.vtkSlicerTransformLogic().hardenTransform(volumeNode)
+                slicer.mrmlScene.RemoveNode(transformNode)
+
+                return alignment_offset_mm
+    
+        def print_orientation(volume_name):
+            node = slicer.util.getNode(volume_name)
+            matrix = vtk.vtkMatrix4x4()
+            node.GetIJKToRASMatrix(matrix)
+            print(f"{volume_name} IJK→RAS:")
+            for i in range(3):
+                print([matrix.GetElement(i, j) for j in range(3)])
+
+
         # Globalni seznami za končno statistiko
         prostate_size_est = []
         ctcbct_distance = []
@@ -560,7 +634,8 @@ class MyTransformModuleLogic(ScriptedLoadableModuleLogic):
 
         for i in range(studyItems.GetNumberOfIds()):
             studyItem = studyItems.GetId(i)
-            
+            studyName = shNode.GetItemName(studyItem)
+            print(f"\nProcessing study: {studyName}")
             # **LOKALNI** seznami, resetirajo se pri vsakem study-ju
             cbct_list = []
             ct_list = []
@@ -630,31 +705,11 @@ class MyTransformModuleLogic(ScriptedLoadableModuleLogic):
                     if not volumeNode or not volumeNode.IsA("vtkMRMLScalarVolumeNode"):
                         print("Can't find volumeNode")
                         #continue  # Preskoči, če ni veljaven volume
-                    
-                    mm_offset, pixel_offset = find_table_top_z(volumeName, writefilecheck, yesCbct)
-                        
-                    if scan_type == "CT":
-                        CT_offset = pixel_offset #Določi koliko je višina CT slike od središča
-                    else:
-                        # Poravnava CBCT z višino CT
-                        CBCT_offset = pixel_offset
-                        alignment_offset = CT_offset - CBCT_offset
-                        print(f"Poravnavam CBCT z CT. Offset: {alignment_offset}")
-
-                        # Uporabi alignment_offset za premik CBCT
-                        transform = vtk.vtkTransform()
-                        transform.Translate(0, alignment_offset, 0)  # Premik samo po y-osi (višina)
-                        
-                        #transformacija volumna oziroma poravnava po višini
-                        transformNode = slicer.vtkMRMLTransformNode()
-                        slicer.mrmlScene.AddNode(transformNode)
-                        transformNode.SetAndObserveTransformToParent(transform)
-                        volumeNode.SetAndObserveTransformNodeID(transformNode.GetID())  # Prilepi transformacijo na CBCT
-                        slicer.mrmlScene.RemoveNode(transformNode)  # Odstrani transformacijo iz scene
                         
                         
                     # Detekcija točk v volumnu
-                    grouped_points = detect_points_region_growing(volumeName, yesCbct, intensity_threshold=3000)
+                    ustvari_marker = False  # Ustvari markerje
+                    grouped_points = detect_points_region_growing(volumeName, yesCbct, ustvari_marker, intensity_threshold=3000)
                     #print(f"Populating volume_points_dict with key ('{scan_type}', '{volumeName}')")
                     volume_points_dict[(scan_type, volumeName)] = grouped_points
                     #print(volume_points_dict)  # Check if the key is correctly added
@@ -663,40 +718,66 @@ class MyTransformModuleLogic(ScriptedLoadableModuleLogic):
             if cbct_list and ct_list:
                 ct_volume_name = ct_list[0]  # Uporabi prvi CT kot referenco
                 print(f"\nProcessing CT: {ct_volume_name}")
-                #yesCbct = False
-                    
+                yesCbct = False
+                                               
+                mm_offset, pixel_offset = find_table_top_z(ct_volume_name, writefilecheck, yesCbct)
+                CT_offset, CT_spacing = align_cbct_to_ct(slicer.util.getNode(ct_volume_name), "CT", mm_offset)                                    
+                     
                 
                 ct_points = [centroid for centroid, _ in volume_points_dict[("CT", ct_volume_name)]]
-
+                
+                print(f"CT points: {ct_points}")
                 if len(ct_points) < 3:
-                    print(f"CT volume {ct_volume_name} doesn't have enough points for registration.")
+                    print(f"CT volume {ct_volume_name} doesn't have enough points for registration. Points: {len(ct_points)}")
+                    continue
                 else:
                     for cbct_volume_name in cbct_list:
-                        print(f"\nProcessing CBCT Volume: {cbct_volume_name}")
-                        #yesCbct = True
-                        cbct_points = [centroid for centroid, _ in volume_points_dict[("CBCT", cbct_volume_name)]]
+                        print(f"\nProcessing CBCT Volume: {cbct_volume_name}")                        
+                        yesCbct = True
+                        scan_type = "CBCT" #redundant but here we are
+                        cbct_volume_node = slicer.util.getNode(cbct_volume_name)
+                        
+                        mm_offset, pixel_offset = find_table_top_z(cbct_volume_name, writefilecheck, yesCbct)                        
+                                            
+                        cbct_points = [centroid for centroid, _ in volume_points_dict[("CBCT", cbct_volume_name)]] #zastareli podatki
+                        cbct_points_array = np.array(cbct_points)  # Pretvorba v numpy array
+                        
+                        #print_orientation(ct_volume_name)
+                        #print_orientation(cbct_volume_name)
+                        
+                        #for i, (cb, ct) in enumerate(zip(cbct_points, ct_points)):
+                        #    print(f"Pair {i}: CBCT {cb}, CT {ct}, diff: {np.linalg.norm(cb - ct):.2f}")
+                        
+                        align_cbct_to_ct(cbct_volume_node, scan_type, mm_offset, CT_offset, CT_spacing)
+                        #cbct_points = [centroid for centroid, _ in volume_points_dict[("CBCT", cbct_volume_name)]] #zastareli podatki
+                        # for i in range(len(cbct_points)):
+                        #     x, y, z = cbct_points[i]
+                        #     y_new = y + mm_offset  # Upoštevaj premik zaradi poravnave mize
+                        #     cbct_points[i] = (x, y_new, z)  # Posodobi marker s premaknjeno višino
+                        ustvari_marker = False  # Ustvari markerje
+                        cbct_points = [centroid for centroid, _ in detect_points_region_growing(cbct_volume_name, yesCbct, ustvari_marker, intensity_threshold=3000)]                     
+                        #cbct_points = detect_points_region_growing(cbct_volume_name, yesCbct, intensity_threshold=3000)   
+                        #cbct_points = [centroid for centroid, _ in volume_points_dict[("CBCT", cbct_volume_name)]] #zastareli podatki
+                                                
                                                 
-                        #find_table_top_z(cbct_volume_name, writefilecheck, yesCbct)
-
                         if len(cbct_points) < 3:
-                            print(f"CBCT Volume '{cbct_volume_name}' doesn't have enough points for registration.")
+                            print(f"CBCT Volume '{cbct_volume_name}' doesn't have enough points for registration. Points: {len(cbct_points)}")
                             continue
-
+                        
+                        
+                        
+                        
+                        #Sortiramo točke po X/Y/Z da se izognemo težavam pri poravnavi
+                        cbct_points = match_points(cbct_points, ct_points)
+                        
+                        #for i, (cb, ct) in enumerate(zip(cbct_points, ct_points)):
+                        #    print(f"Pair {i}: CBCT {cb}, CT {ct}, diff: {np.linalg.norm(cb - ct):.2f}")
+                        
                         # Shranjevanje razdalj
                         distances_ct_cbct = []
                         distances_internal = {"A-B": [], "B-C": [], "C-A": []}
 
-                        cbct_points_array = np.array(cbct_points)  # Pretvorba v numpy array
-
-                        ct_volume_node = slicer.util.getNode(ct_volume_name)
                         cbct_volume_node = slicer.util.getNode(cbct_volume_name)
-                        #ct_spacing = ct_volume_node.GetSpacing()  # (x_spacing, y_spacing, z_spacing)
-                        #cbct_spacing = cbct_volume_node.GetSpacing()  # (x_spacing, y_spacing, z_spacing)
-                        
-                        #ct_scale_factor = np.array(ct_spacing)  # Spacing za CT (x, y, z)                        
-                        #cbct_scale_factor = np.array(cbct_spacing)  # Spacing za CBCT (x, y, z)
-                        #print(ct_scale_factor, cbct_scale_factor)
-
                         
                         # Sortiramo točke po Z-koordinati (ali X/Y, če raje uporabljaš drugo os)
                         cbct_points_sorted = cbct_points_array[np.argsort(cbct_points_array[:, 2])]
@@ -727,10 +808,10 @@ class MyTransformModuleLogic(ScriptedLoadableModuleLogic):
                         # Izberi metodo glede na uporabnikov izbor
                         chosen_rotation_matrix = np.eye(3)
                         chosen_translation_vector = np.zeros(3)
-
+                        print("Markerji pred transformacijo:", cbct_points, ct_points)
                         if applyScaling:
                             scaling_factors = compute_optimal_scaling_per_axis(cbct_points, ct_points)
-                            print("Scaling factors: ", scaling_factors)
+                            #print("Scaling factors: ", scaling_factors)
                             cbct_points = compute_scaling(cbct_points, scaling_factors)
 
                         if applyRotation:
@@ -740,12 +821,13 @@ class MyTransformModuleLogic(ScriptedLoadableModuleLogic):
                                 chosen_rotation_matrix = compute_Horn_rotation(cbct_points, ct_points)
                             elif selectedMethod == "Iterative Closest Point (Kabsch)":
                                 _, chosen_rotation_matrix, _ = icp_algorithm(cbct_points, ct_points)
-                            print("Rotation Matrix:\n", chosen_rotation_matrix)
+                            #print("Rotation Matrix:\n", chosen_rotation_matrix)
 
                         if applyTranslation:
                             chosen_translation_vector = compute_translation(cbct_points, ct_points, chosen_rotation_matrix)
-                            print("Translation Vector:\n", chosen_translation_vector)
-
+                            #print("Translation Vector:\n", chosen_translation_vector)
+                        
+                        
                         # Ustvari vtkTransformNode in ga poveži z CBCT volumenom
                         imeTransformNoda = cbct_volume_name + " Transform"
                         transform_node = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLTransformNode", imeTransformNoda)
@@ -761,7 +843,17 @@ class MyTransformModuleLogic(ScriptedLoadableModuleLogic):
                         # Uporabi transformacijo na volumnu (fizična aplikacija)
                         slicer.vtkSlicerTransformLogic().hardenTransform(cbct_volume_node) #aplicira transformacijo na volumnu
                         slicer.mrmlScene.RemoveNode(transform_node) # Odstrani transformacijo iz scene
-                        print("Transform successful on ", cbct_volume_name)
+                        #print("Transform successful on ", cbct_volume_name)
+                        ustvari_marker = True  # Ustvari markerje
+                        cbct_points = [centroid for centroid, _ in detect_points_region_growing(cbct_volume_name, yesCbct, ustvari_marker, intensity_threshold=3000)]
+                        print("Markerji po transformaciji:\n", cbct_points, ct_points)
+                        
+                        # Izračun razdalj
+                        errors = [np.linalg.norm(cbct - ct) for cbct, ct in zip(cbct_points, ct_points)]
+                        mean_error = np.mean(errors)
+
+                        print("Individualne napake:", errors)
+                        print("📏 Povprečna napaka poravnave: {:.2f} mm".format(mean_error))
 
             else:
                 print(f"Study {studyItem} doesn't have any appropriate CT or CBCT volumes.")
@@ -769,19 +861,35 @@ class MyTransformModuleLogic(ScriptedLoadableModuleLogic):
         # Izpis globalne statistike
         
         
-        if(writefilecheck):
+        if writefilecheck:
             print("Distances between CT & CBCT markers: ", ctcbct_distance)
             print("Distances between pairs of markers for each volume: ", prostate_size_est)
-            # Define the file path for the CSV file
-            file_path = os.path.join(os.path.dirname(__file__), "study_data.csv")
+            
+            # Define file paths
+            prostate_size_file = os.path.join(os.path.dirname(__file__), "prostate_size.csv")
+            ctcbct_distance_file = os.path.join(os.path.dirname(__file__), "ct_cbct_distance.csv")
 
-            # Write lists to the CSV file
-            with open(file_path, mode='w', newline='') as file: #w za write, a za append
+            # Write prostate size data
+            with open(prostate_size_file, mode='w', newline='') as file:
                 writer = csv.writer(file)
-                # Write headers
-                writer.writerow(["Prostate Size", "CT-CBCT Distance"])
-                # Write data rows
-                for i in range(len(prostate_size_est)):
-                    writer.writerow([prostate_size_est[i], ctcbct_distance[i]])
-                print("File written at ", file_path)
+                writer.writerow(["Prostate Size"])
+                for size in prostate_size_est:
+                    writer.writerow([size])
+            print("Prostate size file written at ", prostate_size_file)
 
+            # Write CT-CBCT distance data
+            with open(ctcbct_distance_file, mode='w', newline='') as file:
+                writer = csv.writer(file)
+                writer.writerow(["CT-CBCT Distance"])
+                for distance in ctcbct_distance:
+                    writer.writerow([distance])
+            print("CT-CBCT distance file written at ", ctcbct_distance_file)
+
+        end_time = time.time()
+        # Calculate and print elapsed time
+        elapsed_time = end_time - start_time
+        # Convert to minutes and seconds
+        minutes = int(elapsed_time // 60)
+        seconds = elapsed_time % 60
+
+        print(f"Execution time: {minutes} minutes and {seconds:.6f} seconds")