Error under 1 cm, also added more statistics, file writing and bug fixes.

SeekTransformModule/SeekTransformModule.py

@@ -37,7 +37,7 @@ class SeekTransformModuleWidget(ScriptedLoadableModuleWidget):
 
         # Dropdown menu za izbiro metode
         self.rotationMethodComboBox = qt.QComboBox()
-        self.rotationMethodComboBox.addItems(["Kabsch", "Horn", "Iterative Closest Point (Kabsch)"])
+        self.rotationMethodComboBox.addItems(["Kabsch", "Horn", "Iterative Closest Point (Horn)"])
         self.layout.addWidget(self.rotationMethodComboBox)
 
         # Checkboxi za transformacije
@@ -308,7 +308,7 @@ class MyTransformModuleLogic(ScriptedLoadableModuleLogic):
                 nearest_points = fixed[nearest_indices]
 
                 # Step 2: Compute the optimal rotation and translation
-                R_new = compute_Kabsch_rotation(moving, nearest_points)
+                R_new = compute_Horn_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)
@@ -332,16 +332,14 @@ class MyTransformModuleLogic(ScriptedLoadableModuleLogic):
 
             return moving, R, t
 
-        def match_points(cbct_points, ct_points, auto_weights=True, fallback_if_worse=True, normalize_lengths=True):
+        def match_points(cbct_points, ct_points, auto_weights=True, fallback_if_worse=True, normalize_lengths=True, normalize_angles=False, min_distance=5):
+
             def side_lengths(points):
                 lengths = [
                     np.linalg.norm(points[0] - points[1]),
                     np.linalg.norm(points[1] - points[2]),
                     np.linalg.norm(points[2] - points[0])
                 ]
-                if normalize_lengths:
-                    total = sum(lengths) + 1e-6  # da ne delimo z 0
-                    lengths = [l / total for l in lengths]
                 return lengths
 
             def triangle_angles(points):
@@ -353,22 +351,27 @@ class MyTransformModuleLogic(ScriptedLoadableModuleLogic):
                 angle_C = np.pi - angle_A - angle_B
                 return [angle_A, angle_B, angle_C]
 
+            def normalize(vec):
+                norm = np.linalg.norm(vec)
+                return [v / norm for v in vec] if norm > 0 else vec
+
             def permutation_score(perm, ct_lengths, ct_angles, w_len, w_ang):
                 perm_lengths = side_lengths(perm)
                 perm_angles = triangle_angles(perm)
 
-                # normaliziraj
-                def normalize(vec):
-                    norm = np.linalg.norm(vec)
-                    return vec / norm if norm > 0 else vec
+                # Filter za minimum razdalje
+                if min(perm_lengths) < min_distance:
+                    return float('inf')
+
+                lengths_1 = normalize(perm_lengths) if normalize_lengths else perm_lengths
+                lengths_2 = normalize(ct_lengths)   if normalize_lengths else ct_lengths
 
-                perm_lengths = normalize(perm_lengths)
-                perm_angles = normalize(perm_angles)
-                ct_lengths_n = normalize(ct_lengths)
-                ct_angles_n = normalize(ct_angles)
+                angles_1 = normalize(perm_angles) if normalize_angles else perm_angles
+                angles_2 = normalize(ct_angles)   if normalize_angles else ct_angles
+
+                score_len = sum(abs(a - b) for a, b in zip(lengths_1, lengths_2))
+                score_ang = sum(abs(a - b) for a, b in zip(angles_1, angles_2))
 
-                score_len = sum(abs(a - b) for a, b in zip(perm_lengths, ct_lengths_n))
-                score_ang = sum(abs(a - b) for a, b in zip(perm_angles, ct_angles_n))
                 return w_len * score_len + w_ang * score_ang
 
             cbct_points = list(cbct_points)
@@ -387,22 +390,30 @@ class MyTransformModuleLogic(ScriptedLoadableModuleLogic):
 
             best_perm = None
             best_score = float('inf')
-
+            #print(f"CT points: {ct_points}")
             for perm in itertools.permutations(cbct_points):
                 perm = np.array(perm)
                 score = permutation_score(perm, ct_lengths, ct_angles, weight_length, weight_angle)
+
+                #print(f"Perm: {perm.tolist()} -> Score: {score:.4f}")
+
                 if score < best_score:
                     best_score = score
                     best_perm = perm
 
+            #print("CT centroid:", np.mean(ct_points, axis=0))
+            #print("CBCT centroid (best perm):", np.mean(best_perm, axis=0))
+            
             if fallback_if_worse:
                 original_score = permutation_score(np.array(cbct_points), ct_lengths, ct_angles, weight_length, weight_angle)
-                if original_score <= best_score or np.allclose(cbct_points, best_perm):
-                    print("Fallback to original points due to worse score of the permutation.")
+                #print("Original score: ", original_score)
+                if original_score <= best_score:
+                    #print("Fallback to original points due to worse score of the permutation.")
                     return list(cbct_points)
 
             return list(best_perm)
 
+
         def compute_translation(moving_points, fixed_points, rotation_matrix):
             """
             Computes the translation vector to align moving_points to fixed_points given a rotation matrix.
@@ -428,7 +439,7 @@ class MyTransformModuleLogic(ScriptedLoadableModuleLogic):
 
             return translation
 
-        def create_vtk_transform(rotation_matrix, translation_vector): 
+        def create_vtk_transform(rotation_matrix, translation_vector, study_name=None, cbct_volume_name=None):
             """
             Creates a vtkTransform from a rotation matrix and a translation vector.
             """
@@ -437,6 +448,13 @@ class MyTransformModuleLogic(ScriptedLoadableModuleLogic):
             transform_matrix[:3, :3] = rotation_matrix  # Set rotation part
             transform_matrix[:3, 3] = translation_vector  # Set translation part
 
+            
+            global latest_transform_matrix
+            # Save to global variable
+            latest_transform_matrix = transform_matrix.copy()
+            if study_name and cbct_volume_name:
+                save_transform_matrix(transform_matrix, study_name, cbct_volume_name)
+            
             # Convert to vtkMatrix4x4
             vtk_matrix = vtk.vtkMatrix4x4()
             for i in range(4):
@@ -446,12 +464,28 @@ class MyTransformModuleLogic(ScriptedLoadableModuleLogic):
             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, 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):
+        def save_transform_matrix(matrix, study_name, cbct_volume_name):
+            """
+            Appends the given 4x4 matrix to a text file under the given study folder.
+            """
+            base_folder = "Transformacijske matrike"
+            study_folder = os.path.join(base_folder, study_name)
+            os.makedirs(study_folder, exist_ok=True)  # Create folders if they don't exist
+
+            filename = os.path.join(study_folder, f"{cbct_volume_name}.txt")
+            with open(filename, "a") as f:
+                f.write("Transformacija:\n")
+                for row in matrix:
+                    f.write(" ".join(f"{elem:.6f}" for elem in row) + "\n")
+                f.write("\n")  # Dodaj prazen vrstico med transformacijami
+                
+        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=50, 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.")
@@ -522,107 +556,104 @@ class MyTransformModuleLogic(ScriptedLoadableModuleLogic):
 
             return unique_centroids
 
-        def find_table_top_z(ct_volume_name, writefilecheck, yesCbct):
+        def find_table_top_z(ct_volume_name, writefilecheck, makemarkerscheck, yesCbct):
             """
-            Najde višino zgornjega roba mize v CT/CBCT volumnu in doda markerje.
-            
-            :param ct_volume_name: Ime volumna v slicerju
-            :param writefilecheck: Če je True, zapiše rezultat v CSV
-            :param yesCbct: Če je True, uporabi CBCT thresholde
-            :return: Višina zgornjega roba mize v mm
+            Najde višino zgornjega roba mize v CT/CBCT volumnu in (opcija) doda marker ter shrani podatke.
             """
-            # Pridobi volumen
+            # Pridobi volumen in osnovne lastnosti
             ct_volume_node = slicer.util.getNode(ct_volume_name)
             image_data = ct_volume_node.GetImageData()
             spacing = ct_volume_node.GetSpacing()
             dims = image_data.GetDimensions()
-            #origin = ct_volume_node.GetOrigin()
-
-            # Pretvori volumen v numpy array
             np_array = slicer.util.arrayFromVolume(ct_volume_node)
 
-            # Izračunaj sredinske IJK koordinate
-            mid_ijk = [dims[0] // 2, dims[1] // 2, dims[2] // 2]
+            # Izračun sredinskih koordinat
+            mid_ijk = [max(0, min(dims[i] - 1, dims[i] // 2)) for i in range(3)]
+            mid_z_voxel = mid_ijk[2]
+            mid_x_voxel = max(0, mid_ijk[0] - 15)  # levo od središča
 
-            # Preveri, da so IJK indeksi v mejah volumna
-            mid_ijk = [max(0, min(dims[i] - 1, mid_ijk[i])) for i in range(3)]
+            # Izvoz slice in stolpca intenzitet
+            slice_data = np_array[mid_z_voxel, :, :]
+            column_values = slice_data[:, mid_x_voxel]
 
-            # Pretvorba IJK → RAS
+            # Pridobi transformacijo IJK → RAS
             ijkToRasMatrix = vtk.vtkMatrix4x4()
             ct_volume_node.GetIJKToRASMatrix(ijkToRasMatrix)
 
-            
-
-            # Sredinski Z slice
-            mid_z_voxel = mid_ijk[2]
-            slice_data = np_array[mid_z_voxel, :, :]  # (Y, X)
-
-            # Sredinski stolpec
-            mid_x_voxel = mid_ijk[0] - 15  # 15 pikslov levo od sredine da ne merimo pri hrbtenici
-            column_values = slice_data[:, mid_x_voxel]  # Y smer
-
-            # Doda marker v RAS koordinatah
-            #mid_node = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLMarkupsFiducialNode", f"Sredina_{ct_volume_name}")
-            #mid_ras = np.array(ijkToRasMatrix.MultiplyPoint([*mid_ijk, 1]))[:3]
-            #mid_node.AddControlPoint(mid_ras)
-
-            # Določi threshold glede na CBCT ali CT
-            threshold = -300 if yesCbct else -100 # če je cbct iščemo vrednost -300, sicer pri CT 0
-
-            # Poišči rob mize (drugi dvig)
+            # Iskanje roba mize
+            threshold = -300 if yesCbct else -100
+            min_jump = 100 if yesCbct else 50
             previous_value = -1000
             edge_count = 0
             table_top_y = None
-            min_jump = 100 if yesCbct else 50  # Minimalni skok za CBCT in CT
 
-            for y in range(len(column_values) - 1, -1, -1):  # Od spodaj navzgor
+            for y in range(len(column_values) - 1, -1, -1):
                 intensity = column_values[y]
-                #if column_values[y] > threshold and previous_value <= threshold:
-                if (intensity - previous_value) > min_jump and intensity > threshold:  # Namesto primerjave s threshold
+                if (intensity - previous_value) > min_jump and intensity > threshold:
                     if yesCbct:
-                        table_top_y = y + 1 #Da nismo že v trupu
-                        #print(f"Zgornji rob mize najden pri Y = {table_top_y}") # CBCT
-                        break
-                    if edge_count == 0 or (edge_count == 1 and previous_value < -200):  # Check if intensity is back lower than -400
+                        if column_values[y-1] > -300:
+                            table_top_y = y
+                            break
+                    if edge_count == 0 or (edge_count == 1 and previous_value < -200):
                         edge_count += 1
-                        #print(f"Zaznan rob mize pri X, Y, Z = {mid_x_voxel},{y}, {mid_z_voxel}")
-                        if edge_count == 2:  # Drugi dvig = zgornji rob mize
-                            table_top_y = y + 1
-                            #print(f"Zgornji rob mize najden pri Y = {table_top_y}")
+                        if edge_count == 2:
+                            table_top_y = y
                             break
-                previous_value = column_values[y]
+                previous_value = intensity
 
             if table_top_y is None:
                 print("❌ Zgornji rob mize ni bil najden!")
                 return None
 
-            # Pretvorba Y IJK → RAS
+            # Pretvorba v RAS prostor
             table_ijk = [mid_x_voxel, table_top_y, mid_z_voxel]
             table_ras = np.array(ijkToRasMatrix.MultiplyPoint([*table_ijk, 1]))[:3]
+            z_value = table_ras[2]  # shranimo Z komponento
 
-            # 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)
-
-            # Izračun višine v mm
-            #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)} 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
+            # Dodaj marker, če zahtevano
+            if makemarkerscheck:
+                table_node = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLMarkupsFiducialNode", f"VišinaMize_{ct_volume_name}")
+                table_node.AddControlPoint(table_ras)
+                table_node.SetDisplayVisibility(False)
+
+            # 📝 Shrani osnovni rezultat v heightdata.csv
             if writefilecheck:
-                file_path = os.path.join(os.path.dirname(__file__), "heightdata.csv")
-                with open(file_path, mode='a', newline='') as file:
+                height_file = os.path.join(os.path.dirname(__file__), "heightdata.csv")
+                with open(height_file, 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 = {table_ras[1]:.2f} mm"])
 
-            return table_ras[1], pixel_offset
+            # 📁 Shrani razširjene informacije v table_z_debug.csv
+            if writefilecheck:
+                debug_file = os.path.join(os.path.dirname(__file__), "table_z_debug.csv")
+                neighborhood = 30
+                start = max(0, table_top_y - neighborhood)
+                end = min(len(column_values), table_top_y + 2*neighborhood + 1)
+                values_around = column_values[start:end]
+
+                modality = "CBCT" if yesCbct else "CT"
+                table_z_values[modality] = z_value  # shrani trenutno vrednost
+
+                dz_string = ""
+                if "CT" in table_z_values and "CBCT" in table_z_values:
+                    dz = abs(table_z_values["CBCT"] - table_z_values["CT"])
+                    dz_string = f"{dz:.2f} mm"
+                else:
+                    dz_string = "—"
+
+                with open(debug_file, mode='a', newline='') as file:
+                    writer = csv.writer(file)
+                    writer.writerow([
+                        modality,
+                        ct_volume_name,
+                        f"{z_value:.2f}",        # RAS-Z vrednost
+                        f"slice={table_top_y}",
+                        f"dz={dz_string}"
+                    ] + list(values_around))
+
+            return table_ras[1], table_top_y  # vrni samo Z komponento
+
         
         def align_cbct_to_ct(volumeNode, scan_type, offset, CT_offset=None, CT_spacing=None):
             """
@@ -666,6 +697,28 @@ class MyTransformModuleLogic(ScriptedLoadableModuleLogic):
                 volumeNode.SetAndObserveTransformNodeID(transformNode.GetID())
                 slicer.vtkSlicerTransformLogic().hardenTransform(volumeNode)
                 slicer.mrmlScene.RemoveNode(transformNode)
+                
+                # Poskusi najti ustrezen marker in ga premakniti
+                marker_name = f"VišinaMize_{volumeNode.GetName()}"
+                
+                # Robustno iskanje markerja po imenu
+                table_node = None
+                for node in slicer.util.getNodesByClass("vtkMRMLMarkupsFiducialNode"):
+                    if node.GetName() == marker_name:
+                        table_node = node
+                        break
+
+                if table_node is not None:
+                    current_point = [0, 0, 0]
+                    table_node.GetNthControlPointPosition(0, current_point)
+
+                    moved_point = [
+                        current_point[0],
+                        current_point[1] + alignment_offset_mm,
+                        current_point[2]
+                    ]
+
+                    table_node.SetNthControlPointPosition(0, *moved_point)
 
                 return alignment_offset_mm
     
@@ -734,6 +787,10 @@ class MyTransformModuleLogic(ScriptedLoadableModuleLogic):
                 print(f"✅ Using centroid-based translation (error: {error_centroid:.2f} mm)")
                 return translation_centroid, transformed_centroid, "centroid"
 
+        def rescale_points_to_match_spacing(points, source_spacing, target_spacing):
+            scale_factors = np.array(target_spacing) / np.array(source_spacing)
+            return np.array(points) * scale_factors
+
         def visualize_point_matches_in_slicer(cbct_points, ct_points, study_name="MatchVisualization"):
             assert len(cbct_points) == len(ct_points), "Mora biti enako število točk!"
 
@@ -781,9 +838,24 @@ class MyTransformModuleLogic(ScriptedLoadableModuleLogic):
 
             print(f"✅ Vizualizacija dodana za {study_name} (točke + povezave)")
         
+        def remove_lowest_marker(points, axis=1):
+            """
+            Odstrani točko, ki je najnižja po dani osi (default Y-os v RAS prostoru).
+            """
+            if len(points) <= 3:
+                return points  # Ni potrebe po brisanju
+
+            arr = np.array(points)
+            lowest_index = np.argmin(arr[:, axis])
+            removed = points.pop(lowest_index)
+            print(f"⚠️ Odstranjena najnižja točka (os {axis}): {removed}")
+            return points
+       
+       
         # Globalni seznami za končno statistiko
         prostate_size_est = []
         ctcbct_distance = []
+        table_z_values = {}
 
         # Pridobimo SubjectHierarchyNode
         shNode = slicer.vtkMRMLSubjectHierarchyNode.GetSubjectHierarchyNode(slicer.mrmlScene)
@@ -815,8 +887,12 @@ class MyTransformModuleLogic(ScriptedLoadableModuleLogic):
                 for k in range(finalVolumeItems.GetNumberOfIds()):
                     volumeItem = finalVolumeItems.GetId(k)
                     volumeNode = shNode.GetItemDataNode(volumeItem)
-                    
-                    dicomUIDs = volumeNode.GetAttribute("DICOM.instanceUIDs")
+                    try:
+                        dicomUIDs = volumeNode.GetAttribute("DICOM.instanceUIDs")
+                    except AttributeError:
+                        print(f"⚠️ Volume node '{volumeNode}' ni bil najden ali nima atributa 'DICOM.instanceUIDs'. Preskakujem.")
+                        dicomUIDs = None
+                        continue  # Preskoči, če ni veljaven volume
                     if not dicomUIDs:
                         print("❌ This is an NRRD volume!")
                         continue  # Preskoči, če ni DICOM volume             
@@ -873,16 +949,17 @@ class MyTransformModuleLogic(ScriptedLoadableModuleLogic):
             # Če imamo oba tipa volumna (CBCT in CT) **znotraj istega studyja**
             if cbct_list and ct_list:
                 ct_volume_name = ct_list[0]  # Uporabi prvi CT kot referenco
+                ct_volume_Node = slicer.util.getNode(ct_volume_name)
                 print(f"\nProcessing CT: {ct_volume_name}")
                 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)                                    
+                makemarkerscheck = False                               
+                mm_offset, pixel_offset = find_table_top_z(ct_volume_name, writefilecheck, makemarkerscheck, yesCbct)
+                CT_offset, CT_spacing = align_cbct_to_ct(ct_volume_Node, "CT", mm_offset)                                    
                      
                 
                 ct_points = [centroid for centroid, _ in volume_points_dict[("CT", ct_volume_name)]]
-                
-                print(f"CT points: {ct_points}")
+                ct_points = remove_lowest_marker(ct_points) #odstrani marker v riti, če obstaja
+                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. Points: {len(ct_points)}")
                     continue
@@ -890,26 +967,23 @@ class MyTransformModuleLogic(ScriptedLoadableModuleLogic):
                     for cbct_volume_name in cbct_list:
                         print(f"\nProcessing CBCT Volume: {cbct_volume_name}")                        
                         yesCbct = True
+                        latest_transform_matrix = None
                         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)                   
+                        makemarkerscheck = False  # Ustvari markerje pred poravnavo
+                        mm_offset, pixel_offset = find_table_top_z(cbct_volume_name, writefilecheck, makemarkerscheck, 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)
+                        # 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)   
@@ -920,46 +994,47 @@ class MyTransformModuleLogic(ScriptedLoadableModuleLogic):
                             print(f"CBCT Volume '{cbct_volume_name}' doesn't have enough points for registration. Points: {len(cbct_points)}")
                             continue
                         
-                        
-                        
+                        cbct_spacing = cbct_volume_node.GetSpacing()
+                        ct_spacing = ct_volume_Node.GetSpacing()
+                        cbct_points = rescale_points_to_match_spacing(cbct_points, cbct_spacing, ct_spacing)
                         
                         #Sortiramo točke po X/Y/Z da se izognemo težavam pri poravnavi
                         cbct_points = match_points(cbct_points, ct_points)
                         
                         #visualize_point_matches_in_slicer(cbct_points, ct_points, studyName) #poveže pare markerjev
         
-                        
-                        # Shranjevanje razdalj
-                        distances_ct_cbct = []
-                        distances_internal = {"A-B": [], "B-C": [], "C-A": []}
+                        if writefilecheck:
+                            # Shranjevanje razdalj
+                            distances_ct_cbct = []
+                            distances_internal = {"A-B": [], "B-C": [], "C-A": []}
 
-                        cbct_volume_node = slicer.util.getNode(cbct_volume_name)
-                        
-                        # 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])]
+                            cbct_volume_node = slicer.util.getNode(cbct_volume_name)
+                            
+                            # 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])]
 
-                        # Razdalje med CT in CBCT (SORTIRANE točke!)
-                        d_ct_cbct = np.linalg.norm(cbct_points_sorted - ct_points, axis=1)
-                        distances_ct_cbct.append(d_ct_cbct)
+                            # Razdalje med CT in CBCT (SORTIRANE točke!)
+                            d_ct_cbct = np.linalg.norm(cbct_points_sorted - ct_points, axis=1)
+                            distances_ct_cbct.append(d_ct_cbct)
 
-                        # Razdalje med točkami znotraj SORTIRANIH cbct_points
-                        d_ab = np.linalg.norm(cbct_points_sorted[0] - cbct_points_sorted[1])
-                        d_bc = np.linalg.norm(cbct_points_sorted[1] - cbct_points_sorted[2])
-                        d_ca = np.linalg.norm(cbct_points_sorted[2] - cbct_points_sorted[0])
+                            # Razdalje med točkami znotraj SORTIRANIH cbct_points
+                            d_ab = np.linalg.norm(cbct_points_sorted[0] - cbct_points_sorted[1])
+                            d_bc = np.linalg.norm(cbct_points_sorted[1] - cbct_points_sorted[2])
+                            d_ca = np.linalg.norm(cbct_points_sorted[2] - cbct_points_sorted[0])
 
-                        # Sortiramo razdalje po velikosti, da so vedno v enakem vrstnem redu
-                        sorted_distances = sorted([d_ab, d_bc, d_ca])
+                            # Sortiramo razdalje po velikosti, da so vedno v enakem vrstnem redu
+                            sorted_distances = sorted([d_ab, d_bc, d_ca])
 
-                        distances_internal["A-B"].append(sorted_distances[0])
-                        distances_internal["B-C"].append(sorted_distances[1])
-                        distances_internal["C-A"].append(sorted_distances[2])
-                        
-                        # Dodamo ime študije za v statistiko
-                        studyName = shNode.GetItemName(studyItem)
-                        
-                        # **Shrani razdalje v globalne sezname**
-                        prostate_size_est.append({"Study": studyName, "Distances": sorted_distances})
-                        ctcbct_distance.append({"Study": studyName, "Distances": list(distances_ct_cbct[-1])})  # Pretvorimo v seznam
+                            distances_internal["A-B"].append(sorted_distances[0])
+                            distances_internal["B-C"].append(sorted_distances[1])
+                            distances_internal["C-A"].append(sorted_distances[2])
+                            
+                            # Dodamo ime študije za v statistiko
+                            studyName = shNode.GetItemName(studyItem)
+                            
+                            # **Shrani razdalje v globalne sezname**
+                            prostate_size_est.append({"Study": studyName, "Distances": sorted_distances})
+                            ctcbct_distance.append({"Study": studyName, "Distances": list(distances_ct_cbct[-1])})  # Pretvorimo v seznam
 
                         # Izberi metodo glede na uporabnikov izbor
                         chosen_rotation_matrix = np.eye(3)
@@ -983,47 +1058,95 @@ class MyTransformModuleLogic(ScriptedLoadableModuleLogic):
                             elif selectedMethod == "Iterative Closest Point (Kabsch)":
                                 _, chosen_rotation_matrix, _ = icp_algorithm(cbct_points, ct_points)
                             #print("Rotation Matrix:\n", chosen_rotation_matrix)
-
-                        if applyTranslation:
-                            chosen_translation_vector, cbct_points_transformed, method_used = choose_best_translation(cbct_points, ct_points, chosen_rotation_matrix) #Izbere optimalno translacijo
                             
-                            per_axis_error = np.abs(np.array(cbct_points_transformed) - np.array(ct_points))
-                            dz_mean = np.mean(per_axis_error[:, 2])
-                            print(f"per-axis error: {per_axis_error}, dz_mean err: {dz_mean:.2f} mm")
+                        fine_shift = np.zeros(3)  # Inicializiraj fine premike
+                        if applyTranslation:
+                            chosen_translation_vector, cbct_points_transformed, method_used = choose_best_translation(
+                                cbct_points, ct_points, chosen_rotation_matrix)
                             
-                            #chosen_translation_vector = compute_translation(cbct_points, ct_points, chosen_rotation_matrix)
-                            #print("Translation Vector:\n", chosen_translation_vector)
+                            # Sistematična razlika (signed shift)
+                            rotated_cbct = np.dot(cbct_points, chosen_rotation_matrix.T)
+                            translated_cbct = rotated_cbct + chosen_translation_vector
+                            delta_y_list = [ct[1] - cbct[1] for ct, cbct in zip(ct_points, translated_cbct)]
+                            mean_delta_y = np.mean(delta_y_list)
 
-                        
-                        
-                        
-                        # Ustvari vtkTransformNode in ga poveži z CBCT volumenom
+
+                            # Uporabi sistematični shift za dodatno poravnavo v y-osi
+                            fine_shift = np.array([0.0, mean_delta_y, 0.0])  # samo Y-os
+                            cbct_points_transformed += fine_shift
+
+                            print(f"🔧 Fine correction shifts: ΔX={fine_shift[0]:.2f} mm, ΔY={fine_shift[1]:.2f} mm, ΔZ={fine_shift[2]:.2f} mm")
+
+
+
+                        # ✅ Kombinirana transformacija
+                        total_translation = chosen_translation_vector + fine_shift
+                        chosen_translation_vector = total_translation
+                        vtk_transform = create_vtk_transform(chosen_rotation_matrix, chosen_translation_vector, studyName, cbct_volume_name) #Tukaj se tudi izpiše transformacijska matrika
+
+                        # 🔄 Pripni transformacijo
                         imeTransformNoda = cbct_volume_name + " Transform"
                         transform_node = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLTransformNode", imeTransformNoda)
-
-                        # Kreiraj transformacijo in jo uporabi
-                        vtk_transform = create_vtk_transform(chosen_rotation_matrix, chosen_translation_vector)
                         transform_node.SetAndObserveTransformToParent(vtk_transform)
 
-                        # Pridobi CBCT volumen in aplikacijo transformacije
                         cbct_volume_node = slicer.util.getNode(cbct_volume_name)
                         cbct_volume_node.SetAndObserveTransformNodeID(transform_node.GetID())
 
-                        # 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)
-                        ustvari_marker = True  # Ustvari markerje
+                        # 🔨 Uporabi (ali shrani transformacijo kasneje)
+                        slicer.vtkSlicerTransformLogic().hardenTransform(cbct_volume_node)
+                        slicer.mrmlScene.RemoveNode(transform_node)
+
+                        # 📍 Detekcija markerjev po transformaciji
+                        ustvari_marker = False
                         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
+                        #popravek v x osi
+                        delta_x_list = [ct[0] - cbct[0] for ct, cbct in zip(ct_points, cbct_points)]
+                        mean_delta_x = np.mean(delta_x_list)
+                        #popravek v y osi
+                        delta_y_list = [ct[1] - cbct[1] for ct, cbct in zip(ct_points, cbct_points)]
+                        mean_delta_y = np.mean(delta_y_list)
+                        #popravek v z osi
+                        delta_z_list = [ct[2] - cbct[2] for ct, cbct in zip(ct_points, cbct_points)]
+                        mean_delta_z = np.mean(delta_z_list)
+
+                        # Uporabi sistematični shift za dodatno poravnavo
+                        fine_shift = np.array([mean_delta_x, mean_delta_y, mean_delta_z])
+                        #cbct_points_transformed += fine_shift
+                                                
+                    
+                        chosen_rotation_matrix = np.eye(3) #tokrat brez rotacije
+                        vtk_transform = create_vtk_transform(chosen_rotation_matrix, fine_shift, studyName, cbct_volume_name) #Tukaj se tudi izpiše transformacijska matrika
+
+                        # 🔄 Pripni transformacijo
+                        imeTransformNoda = cbct_volume_name + " Transform"
+                        transform_node = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLTransformNode", imeTransformNoda)
+                        transform_node.SetAndObserveTransformToParent(vtk_transform)
+
+                        cbct_volume_node = slicer.util.getNode(cbct_volume_name)
+                        cbct_volume_node.SetAndObserveTransformNodeID(transform_node.GetID())
+
+                        # 🔨 Uporabi (ali shrani transformacijo kasneje)
+                        slicer.vtkSlicerTransformLogic().hardenTransform(cbct_volume_node)
+                        slicer.mrmlScene.RemoveNode(transform_node)
+
+                        
+                        ustvari_marker = True
+                        cbct_points = [centroid for centroid, _ in detect_points_region_growing(cbct_volume_name, yesCbct, ustvari_marker, intensity_threshold=3000)]
+                        print(f"🔧 Fine correction shifts: ΔX={fine_shift[0]:.2f} mm, ΔY={fine_shift[1]:.2f} mm, ΔZ={fine_shift[2]:.2f} mm")
+
+                        # 📏 Izračun napake
                         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))
                         
+                        for i, (cbct, ct) in enumerate(zip(cbct_points, ct_points)):
+                            diff = np.array(cbct) - np.array(ct)
+                            print(f"Točka {i+1}: ΔX={diff[0]:.2f} mm, ΔY={diff[1]:.2f} mm, ΔZ={diff[2]:.2f} mm")
+                        
                         
 
             else:
@@ -1033,8 +1156,8 @@ class MyTransformModuleLogic(ScriptedLoadableModuleLogic):
         
         
         if writefilecheck:
-            print("Distances between CT & CBCT markers: ", ctcbct_distance)
-            print("Distances between pairs of markers for each volume: ", prostate_size_est)
+            #print("Distances between CT & CBCT markers: ", ctcbct_distance)
+            #print("Distances between pairs of markers for each volume: ", prostate_size_est)
             
             # Define file paths
             prostate_size_file = os.path.join(os.path.dirname(__file__), "prostate_size.csv")
@@ -1046,7 +1169,7 @@ class MyTransformModuleLogic(ScriptedLoadableModuleLogic):
                 writer.writerow(["Prostate Size"])
                 for size in prostate_size_est:
                     writer.writerow([size])
-            print("Prostate size file written at ", prostate_size_file)
+            #print("Prostate size file written at ", prostate_size_file)
 
             # Write CT-CBCT distance data
             with open(ctcbct_distance_file, mode='w', newline='') as file:
@@ -1054,7 +1177,7 @@ class MyTransformModuleLogic(ScriptedLoadableModuleLogic):
                 writer.writerow(["CT-CBCT Distance"])
                 for distance in ctcbct_distance:
                     writer.writerow([distance])
-            print("CT-CBCT distance file written at ", ctcbct_distance_file)
+            #print("CT-CBCT distance file written at ", ctcbct_distance_file)
 
         end_time = time.time()
         # Calculate and print elapsed time