Out of beta, working version.

find_table_top_z finetuned,
scaling improved,
added writing transform matrix to file
elapse time keeping
more statistic

SeekTransformModule/SeekTransformModule.py

@@ -1,12 +1,13 @@
 import os
 import numpy as np
 import scipy
+import re
 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
+from collections import deque, Counter
 import vtk
 from slicer.ScriptedLoadableModule import *
 import qt
@@ -16,6 +17,8 @@ import time
 
 #exec(open("C:/Users/lkomar/Documents/Prostata/FirstTryRegister.py").read())
 
+cumulative_matrices = {}
+
 class SeekTransformModule(ScriptedLoadableModule):
     """
     Module description shown in the module panel.
@@ -41,6 +44,10 @@ class SeekTransformModuleWidget(ScriptedLoadableModuleWidget):
         self.layout.addWidget(self.rotationMethodComboBox)
 
         # Checkboxi za transformacije
+        self.scalingCheckBox = qt.QCheckBox("Scaling")
+        self.scalingCheckBox.setChecked(True)
+        self.layout.addWidget(self.scalingCheckBox)
+        
         self.rotationCheckBox = qt.QCheckBox("Rotation")
         self.rotationCheckBox.setChecked(True)
         self.layout.addWidget(self.rotationCheckBox)
@@ -49,9 +56,9 @@ class SeekTransformModuleWidget(ScriptedLoadableModuleWidget):
         self.translationCheckBox.setChecked(True)
         self.layout.addWidget(self.translationCheckBox)
 
-        self.scalingCheckBox = qt.QCheckBox("Scaling")
-        self.scalingCheckBox.setChecked(True)
-        self.layout.addWidget(self.scalingCheckBox)
+        self.markersCheckBox = qt.QCheckBox("Place control points for detected markers")
+        self.markersCheckBox.setChecked(True)
+        self.layout.addWidget(self.markersCheckBox)
         
         self.writefileCheckBox = qt.QCheckBox("Write distances to csv file")
         self.writefileCheckBox.setChecked(True)
@@ -76,10 +83,11 @@ class SeekTransformModuleWidget(ScriptedLoadableModuleWidget):
         applyRotation = self.rotationCheckBox.isChecked()
         applyTranslation = self.translationCheckBox.isChecked()
         applyScaling = self.scalingCheckBox.isChecked()
+        applyMarkers = self.markersCheckBox.isChecked()
         writefilecheck = self.writefileCheckBox.isChecked()
 
         # Pokliči logiko z izbranimi nastavitvami
-        logic.run(selectedMethod, applyRotation, applyTranslation, applyScaling, writefilecheck)
+        logic.run(selectedMethod, applyRotation, applyTranslation, applyScaling, applyMarkers, writefilecheck)
 
 
 class MyTransformModuleLogic(ScriptedLoadableModuleLogic):
@@ -88,7 +96,7 @@ class MyTransformModuleLogic(ScriptedLoadableModuleLogic):
     """
     
     
-    def run(self, selectedMethod, applyRotation, applyTranslation, applyScaling, writefilecheck):
+    def run(self, selectedMethod, applyRotation, applyTranslation, applyScaling, applymarkers, writefilecheck):
         start_time = time.time()
         print("Calculating...")
         
@@ -139,35 +147,17 @@ class MyTransformModuleLogic(ScriptedLoadableModuleLogic):
 
             return region
 
-        
-        def compute_optimal_scaling_per_axis(moving_points, fixed_points):
-            """Computes optimal scaling factors for each axis (X, Y, Z) to align moving points (CBCT) to fixed points (CT).
-
-            Args:
-                moving_points (list of lists): List of (x, y, z) moving points (CBCT).
-                fixed_points (list of lists): List of (x, y, z) fixed points (CT).
-
-            Returns:
-                tuple: Scaling factors (sx, sy, sz).
-            """
-            moving_points_np = np.array(moving_points)
-            fixed_points_np = np.array(fixed_points)
-
-            # Compute centroids
-            centroid_moving = np.mean(moving_points_np, axis=0)
-            centroid_fixed = np.mean(fixed_points_np, axis=0)
-
-            # Compute absolute distances of each point from its centroid along each axis
-            distances_moving = np.abs(moving_points_np - centroid_moving)
-            distances_fixed = np.abs(fixed_points_np - centroid_fixed)
-
-            # Compute scaling factors as the ratio of mean absolute distances per axis
-            scale_factors = np.mean(distances_fixed, axis=0) / np.mean(distances_moving, axis=0)
+        def compute_scaling_stddev(moving_points, fixed_points):
+            moving = np.array(moving_points)
+            fixed = np.array(fixed_points)
 
-            return tuple(scale_factors)
+            # Standard deviation around centroid, po osi
+            scaling_factors = np.std(fixed, axis=0) / np.std(moving, axis=0)
+            return tuple(scaling_factors)
 
         def compute_scaling(cbct_points, scaling_factors):
             """Applies non-uniform scaling to CBCT points.
+            
 
             Args:
                 cbct_points (list of lists): List of (x, y, z) points.
@@ -181,6 +171,11 @@ class MyTransformModuleLogic(ScriptedLoadableModuleLogic):
 
             cbct_points_np = np.array(cbct_points)  # Convert to numpy array
             scaled_points = cbct_points_np @ scaling_matrix.T  # Apply scaling
+            
+            scaling_4x4 = np.eye(4)
+            scaling_4x4[0, 0] = sx
+            scaling_4x4[1, 1] = sy
+            scaling_4x4[2, 2] = sz
 
             return scaled_points.tolist()  # Convert back to list
 
@@ -223,7 +218,6 @@ class MyTransformModuleLogic(ScriptedLoadableModuleLogic):
 
             return Rotate_optimal
 
-
         def compute_Horn_rotation(moving_points, fixed_points):
             """
             Computes the optimal rotation matrix using quaternions.
@@ -413,7 +407,6 @@ class MyTransformModuleLogic(ScriptedLoadableModuleLogic):
 
             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.
@@ -439,54 +432,85 @@ class MyTransformModuleLogic(ScriptedLoadableModuleLogic):
 
             return translation
 
-        def create_vtk_transform(rotation_matrix, translation_vector, study_name=None, cbct_volume_name=None):
+        def create_vtk_transform(rotation_matrix, translation_vector, tablefound, study_name=None, cbct_volume_name=None, scaling_factors=None):
             """
-            Creates a vtkTransform from a rotation matrix and a translation vector.
+            Creates a vtkTransform from scaling, rotation, and translation.
+            Shrani tudi kumulativno matriko v globalni slovar cumulative_matrices.
             """
-            # Create a 4x4 transformation matrix
-            transform_matrix = np.eye(4)  # Start with an identity matrix
-            transform_matrix[:3, :3] = rotation_matrix  # Set rotation part
-            transform_matrix[:3, 3] = translation_vector  # Set translation part
+            # ----- Inicializacija -----
+            global cumulative_matrices
+            transform = vtk.vtkTransform()
 
-            
-            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
+            # ----- 1. Skaliranje -----
+            if scaling_factors is not None:
+                sx, sy, sz = scaling_factors
+                transform.Scale(sx, sy, sz)
+
+            # ----- 2. Rotacija -----
+            # Rotacijsko matriko in translacijo pretvori v homogeno matriko
+            affine_matrix = np.eye(4)
+            affine_matrix[:3, :3] = rotation_matrix
+            affine_matrix[:3, 3] = translation_vector
+
+            # Vstavi v vtkMatrix4x4
             vtk_matrix = vtk.vtkMatrix4x4()
             for i in range(4):
                 for j in range(4):
-                    vtk_matrix.SetElement(i, j, transform_matrix[i, j])
+                    vtk_matrix.SetElement(i, j, affine_matrix[i, j])
+
+            transform.Concatenate(vtk_matrix)
+
+            # ----- 3. Debug izpis -----
             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)
+                print(" ".join(f"{vtk_matrix.GetElement(i, j):.6f}" for j in range(4)))
+
+            # ----- 4. Shrani v kumulativni matriki -----
+            if study_name and cbct_volume_name:
+                key = (study_name, cbct_volume_name)
+
+                if key not in cumulative_matrices:
+                    cumulative_matrices[key] = np.eye(4)
+
+                cumulative_matrices[key] = np.dot(cumulative_matrices[key], affine_matrix)
+
             return transform
         
-        def save_transform_matrix(matrix, study_name, cbct_volume_name):
+        def matrix_to_array(vtk_transform):
+            """
+            Converts a vtkTransform to a 4x4 numpy array.
+            """
+            vtk_matrix = vtk.vtkMatrix4x4()
+            vtk_transform.GetMatrix(vtk_matrix)
+            np_matrix = np.zeros((4, 4))
+            for i in range(4):
+                for j in range(4):
+                    np_matrix[i, j] = vtk_matrix.GetElement(i, j)
+            return np_matrix
+        
+        def save_transform_matrix(matrix, study_name, cbct_volume_name, ct_table_found):
             """
             Appends the given 4x4 matrix to a text file under the given study folder.
             """
-            base_folder = "Transformacijske matrike"
+            base_folder = os.path.join(os.path.dirname(__file__), "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")
+            safe_cbct_name = re.sub(r'[<>:"/\\|?*]', '_', cbct_volume_name)
+            # Preveri ali je CT miza najdena
+            if ct_table_found:
+                safe_cbct_name += "_MizaPoravnava"
+            else:
+                safe_cbct_name += "_NIMizaPoravnana"
+            filename = os.path.join(study_folder, f"{safe_cbct_name}.txt")
             with open(filename, "a") as f:
-                f.write("Transformacija:\n")
+                #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
+            print(f"Transform matrix saved to {filename}")
                 
         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)
+            volume_node = find_volume_node_by_partial_name(volume_name)
             if not volume_node:
                 raise RuntimeError(f"Volume {volume_name} not found.")
             
@@ -558,101 +582,111 @@ class MyTransformModuleLogic(ScriptedLoadableModuleLogic):
 
         def find_table_top_z(ct_volume_name, writefilecheck, makemarkerscheck, yesCbct):
             """
-            Najde višino zgornjega roba mize v CT/CBCT volumnu in (opcija) doda marker ter shrani podatke.
+            Najde višino zgornjega roba mize v CT/CBCT volumnu in po želji doda marker v sceno.
+
+            Args:
+                ct_volume_name (str): Ime volumna.
+                writefilecheck (bool): Ali naj se rezultat shrani v .csv.
+                makemarkerscheck (bool): Ali naj se doda marker v 3D Slicer.
+                yesCbct (bool): True, če je CBCT; False, če je CT.
+
+            Returns:
+                (float, int): Z komponenta v RAS prostoru, in Y indeks v slicerjevem volumnu.
             """
-            # 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()
-            np_array = slicer.util.arrayFromVolume(ct_volume_node)
-
-            # 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
-
-            # Izvoz slice in stolpca intenzitet
-            slice_data = np_array[mid_z_voxel, :, :]
-            column_values = slice_data[:, mid_x_voxel]
-
-            # Pridobi transformacijo IJK → RAS
+            # --- Pridobi volume node ---
+            ct_volume_node = find_volume_node_by_partial_name(ct_volume_name)
+            np_array = slicer.util.arrayFromVolume(ct_volume_node)  # (Z, Y, X)
             ijkToRasMatrix = vtk.vtkMatrix4x4()
             ct_volume_node.GetIJKToRASMatrix(ijkToRasMatrix)
 
-            # Iskanje roba mize
-            threshold = -300 if yesCbct else -100
-            min_jump = 100 if yesCbct else 50
-            previous_value = -1000
-            edge_count = 0
+            # --- Določimo lokacijo stolpca ---
+            z_index = np_array.shape[0] // 2  # srednji slice
+            y_size = np_array.shape[1]
+            # x_index = int(np_array.shape[2] * 0.15)
+            # x_index = max(0, min(x_index, np_array.shape[2] - 1))
+            
+
+            # --- Izračun spodnje tretjine (spodnji del slike) ---
+            y_start = int(y_size * 2 / 3)
+            slice_data = np_array[z_index, :, :]  # (Y, X)
+            y_end = y_size  # Do dna slike
+            #column_values = slice_data[y_start:y_end, x_index]  # (Y)
+
+            # --- Parametri za rob ---
+            threshold_high = -300 if yesCbct else -100
+            threshold_low = -700 if yesCbct else -350
+            min_jump = 100 if yesCbct else 100
+            window_size = 4  # število voxelov nad/pod
+            #previous_value = column_values[-1]
             table_top_y = None
 
-            for y in range(len(column_values) - 1, -1, -1):
-                intensity = column_values[y]
-                if (intensity - previous_value) > min_jump and intensity > threshold:
-                    if yesCbct:
-                        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
-                        if edge_count == 2:
-                            table_top_y = y
-                            break
-                previous_value = intensity
+            # --- Več stolpcev okoli x_index ---
+            x_center = np_array.shape[2] // 2
+            x_offset = 30  # 30 levo od sredine
+            x_index_base = max(0, x_center - x_offset)
+            candidate_y_values = []
+
+            search_range = range(-5, 6)  # od -5 do +5 stolpcev
+            
+            for dx in search_range:
+                x_index = x_index_base + dx
+                if x_index < 0 or x_index >= np_array.shape[2]:
+                    continue
+
+                column_values = slice_data[y_start:y_end, x_index]
+
+                for i in range(window_size, len(column_values) - window_size):
+                    curr = column_values[i]
+                    above_avg = np.mean(column_values[i - window_size:i])
+                    below_avg = np.mean(column_values[i + 1:i + 1 + window_size])
+
+                    if (threshold_low < curr < threshold_high
+                        and (above_avg - below_avg) > min_jump
+                        and below_avg < -400
+                        and above_avg > -300):
+                        y_found = y_start + i
+                        candidate_y_values.append(y_found)
+                        break  # samo prvi zadetek v stolpcu
+            if candidate_y_values:
+                most_common_y, _ = Counter(candidate_y_values).most_common(1)[0]
+                table_top_y = most_common_y
+                print(f"✅ Rob mize (najpogostejši Y): {table_top_y}, pojavitev: {candidate_y_values.count(table_top_y)}×")
+
+            
+            """ # --- Poišči skok navzdol pod prag (od spodaj navzgor) ---
+            for i in range(len(column_values) - 2, -1, -1):  # od spodaj proti vrhu
+                intensity = column_values[i]
+                if (intensity - previous_value) > min_jump and intensity < thresholdhigh and intensity > thresholdlow:
+                    table_top_y = y_start + i - 1
+                    print(f"✅ Rob mize najden pri Y = {table_top_y}, intenziteta = {intensity}")
+                    print("Column values (partial):", column_values.tolist())
+                    break
+                previous_value = intensity """
 
             if table_top_y is None:
-                print("❌ Zgornji rob mize ni bil najden!")
+                print(f"⚠️ Rob mize ni bil najden (X = {x_index})")
+                print("Column values (partial):", column_values.tolist())
                 return None
 
-            # Pretvorba v RAS prostor
-            table_ijk = [mid_x_voxel, table_top_y, mid_z_voxel]
+            # --- Pretvorba v RAS koordinato ---
+            table_ijk = [x_index, table_top_y, z_index]
             table_ras = np.array(ijkToRasMatrix.MultiplyPoint([*table_ijk, 1]))[:3]
-            z_value = table_ras[2]  # shranimo Z komponento
 
-            # Dodaj marker, če zahtevano
+            # --- Marker ---
             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
+            # --- Shrani v CSV ---
             if writefilecheck:
                 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"])
-
-            # 📁 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))
+                    modality = "CBCT" if yesCbct else "CT"
+                    writer.writerow([modality, ct_volume_name, f"Upper table edge at Z = {table_ras[1]:.2f} mm"])
 
-            return table_ras[1], table_top_y  # vrni samo Z komponento
+            return table_ras[1], table_top_y
 
         
         def align_cbct_to_ct(volumeNode, scan_type, offset, CT_offset=None, CT_spacing=None):
@@ -723,7 +757,7 @@ class MyTransformModuleLogic(ScriptedLoadableModuleLogic):
                 return alignment_offset_mm
     
         def print_orientation(volume_name):
-            node = slicer.util.getNode(volume_name)
+            node = find_volume_node_by_partial_name(volume_name)
             matrix = vtk.vtkMatrix4x4()
             node.GetIJKToRASMatrix(matrix)
             print(f"{volume_name} IJK→RAS:")
@@ -747,7 +781,7 @@ class MyTransformModuleLogic(ScriptedLoadableModuleLogic):
             ct_centroid = np.mean(ct_points, axis=0)
             translation_vector = ct_centroid - cbct_centroid
             aligned_cbct = cbct_points + translation_vector
-            return aligned_cbct
+            return aligned_cbct, translation_vector
         
         def choose_best_translation(cbct_points, ct_points, rotation_matrix):
             """
@@ -781,10 +815,10 @@ class MyTransformModuleLogic(ScriptedLoadableModuleLogic):
 
             # 4. Izberi boljšo
             if error_recomputed < error_centroid:
-                print(f"✅ Using retranslation (error: {error_recomputed:.2f} mm)")
+                #print(f"✅ Using retranslation (error: {error_recomputed:.2f} mm)")
                 return translation_recomputed, transformed_recomputed, "retranslation"
             else:
-                print(f"✅ Using centroid-based translation (error: {error_centroid:.2f} mm)")
+                #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):
@@ -851,7 +885,27 @@ class MyTransformModuleLogic(ScriptedLoadableModuleLogic):
             print(f"⚠️ Odstranjena najnižja točka (os {axis}): {removed}")
             return points
        
-       
+        def update_timing_csv(timing_data, study_name):
+            file_path = os.path.join(os.path.dirname(__file__), "timing_summary.csv")
+            file_exists = os.path.isfile(file_path)
+
+            with open(file_path, mode='a', newline='') as csvfile:
+                fieldnames = ["Study", "IO", "Fixing", "Scaling", "CentroidAlign", "Rotation", "Translation", "Transform", "FileSave", "Total"]
+                writer = csv.DictWriter(csvfile, fieldnames=fieldnames)
+
+                if not file_exists:
+                    writer.writeheader()
+
+                row = {"Study": study_name}
+                row.update(timing_data)
+                writer.writerow(row)
+                
+        def find_volume_node_by_partial_name(partial_name):
+            for node in slicer.util.getNodesByClass("vtkMRMLScalarVolumeNode"):
+                if partial_name in node.GetName():
+                    return node
+            raise RuntimeError(f"❌ Volume with name containing '{partial_name}' not found.")
+                
         # Globalni seznami za končno statistiko
         prostate_size_est = []
         ctcbct_distance = []
@@ -862,8 +916,10 @@ class MyTransformModuleLogic(ScriptedLoadableModuleLogic):
         
         studyItems = vtk.vtkIdList()
         shNode.GetItemChildren(shNode.GetSceneItemID(), studyItems)
-
+        
         for i in range(studyItems.GetNumberOfIds()):
+            study_start_time = time.time()
+            start_io = time.time()
             studyItem = studyItems.GetId(i)
             studyName = shNode.GetItemName(studyItem)
             print(f"\nProcessing study: {studyName}")
@@ -890,7 +946,7 @@ class MyTransformModuleLogic(ScriptedLoadableModuleLogic):
                     try:
                         dicomUIDs = volumeNode.GetAttribute("DICOM.instanceUIDs")
                     except AttributeError:
-                        print(f"⚠️ Volume node '{volumeNode}' ni bil najden ali nima atributa 'DICOM.instanceUIDs'. Preskakujem.")
+                        print(f"⚠️ Volume node '{volumeNode}' not found or no attribute 'DICOM.instanceUIDs'. Skip.")
                         dicomUIDs = None
                         continue  # Preskoči, če ni veljaven volume
                     if not dicomUIDs:
@@ -947,13 +1003,23 @@ class MyTransformModuleLogic(ScriptedLoadableModuleLogic):
                     #print(volume_points_dict)  # Check if the key is correctly added
 
             # Če imamo oba tipa volumna (CBCT in CT) **znotraj istega studyja**
+            end_io = time.time()
             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)
+                ct_volume_Node = find_volume_node_by_partial_name(ct_volume_name)
                 print(f"\nProcessing CT: {ct_volume_name}")
                 yesCbct = False
-                makemarkerscheck = False                               
-                mm_offset, pixel_offset = find_table_top_z(ct_volume_name, writefilecheck, makemarkerscheck, yesCbct)
+                makemarkerscheck = True                               
+                result = find_table_top_z(ct_volume_name, writefilecheck, makemarkerscheck, yesCbct)
+                if result is not None:
+                    mm_offset, pixel_offset = result
+                    ct_table_found = True
+                    #print(f"✔️ Poravnava z višino mize: ΔY = {mm_offset:.2f} mm")
+                    # Dodaj ΔY k translaciji ali transformaciji po potrebi
+                else:
+                    print("⚠️ Table top not found – continue without correction on Y axis.")
+                    mm_offset = 0.0  # ali None, če želiš eksplicitno ignorirati
+                    ct_table_found = False
                 CT_offset, CT_spacing = align_cbct_to_ct(ct_volume_Node, "CT", mm_offset)                                    
                      
                 
@@ -967,11 +1033,35 @@ 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)
-                        makemarkerscheck = False  # Ustvari markerje pred poravnavo
-                        mm_offset, pixel_offset = find_table_top_z(cbct_volume_name, writefilecheck, makemarkerscheck, yesCbct)                   
+                        cbct_volume_node = find_volume_node_by_partial_name(cbct_volume_name)
+                        
+                        key = (studyName, cbct_volume_name)
+                        if key not in cumulative_matrices:
+                            cumulative_matrices[key] = np.eye(4)
+                        
+                        fixing = time.time()
+                        makemarkerscheck = False  # Ustvari CBCT miza markerje pred poravnavo
+                        if(ct_table_found):
+                            result = find_table_top_z(cbct_volume_name, writefilecheck, makemarkerscheck, yesCbct) #!!!!!!!!!!!!!???????????? ct_volume_name
+                            if result is not None:
+                                mm_offset, pixel_offset = result
+                                #print(f"✔️ Poravnava z višino mize: ΔY = {mm_offset:.2f} mm")
+                                skupni_offset = align_cbct_to_ct(cbct_volume_node, scan_type, mm_offset, CT_offset, CT_spacing) #poravna CBCT in sporoči skupni offset                                
+                                
+                                table_shift_matrix = np.eye(4)
+                                table_shift_matrix[1, 3] = skupni_offset  # Premik po Y
+
+                                # Pomnoži obstoječo kumulativno matriko
+                                key = (studyName, cbct_volume_name)
+                                if key not in cumulative_matrices:
+                                    cumulative_matrices[key] = np.eye(4)
+
+                                cumulative_matrices[key] = np.dot(cumulative_matrices[key], table_shift_matrix)
+                                
+                            else:
+                                print("⚠️ Table top not found – continue without correction on Y axis.")
+                                mm_offset = 0.0             
                                             
                         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
@@ -982,7 +1072,7 @@ class MyTransformModuleLogic(ScriptedLoadableModuleLogic):
                         #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)
+                        
                         
                         ustvari_marker = False  # Ustvari markerje
                         cbct_points = [centroid for centroid, _ in detect_points_region_growing(cbct_volume_name, yesCbct, ustvari_marker, intensity_threshold=3000)]                     
@@ -1000,7 +1090,7 @@ class MyTransformModuleLogic(ScriptedLoadableModuleLogic):
                         
                         #Sortiramo točke po X/Y/Z da se izognemo težavam pri poravnavi
                         cbct_points = match_points(cbct_points, ct_points)
-                        
+                        fixing_end = time.time()
                         #visualize_point_matches_in_slicer(cbct_points, ct_points, studyName) #poveže pare markerjev
         
                         if writefilecheck:
@@ -1008,7 +1098,7 @@ class MyTransformModuleLogic(ScriptedLoadableModuleLogic):
                             distances_ct_cbct = []
                             distances_internal = {"A-B": [], "B-C": [], "C-A": []}
 
-                            cbct_volume_node = slicer.util.getNode(cbct_volume_name)
+                            cbct_volume_node = find_volume_node_by_partial_name(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])]
@@ -1039,17 +1129,29 @@ 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)
+                        #print("Markerji pred transformacijo:", cbct_points, ct_points)
+                        start_scaling = time.time()
+                        scaling_factors = None
                         if applyScaling:
-                            scaling_factors = compute_optimal_scaling_per_axis(cbct_points, ct_points)
+                            scaling_factors = compute_scaling_stddev(cbct_points, ct_points)
                             #print("Scaling factors: ", scaling_factors)
                             cbct_points = compute_scaling(cbct_points, scaling_factors)
+                            
                         
+                        
+                            
+                            
+                        end_scaling = time.time()
+                        start_align = time.time()
                         initial_error = np.mean(np.linalg.norm(np.array(cbct_points) - np.array(ct_points), axis=1))
                         if initial_error > 30:
-                            print("⚠️ Initial distance too large, applying centroid prealignment.")
-                            cbct_points = prealign_by_centroid(cbct_points, ct_points)
+                            #print("⚠️ Initial distance too large, applying centroid prealignment.")
+                            cbct_points, transvector = prealign_by_centroid(cbct_points, ct_points)
+                        else:  
+                            transvector = np.zeros(3)   
+                        end_align = time.time()
                         
+                        start_rotation = time.time()
                         if applyRotation:
                             if selectedMethod == "Kabsch":
                                 chosen_rotation_matrix = compute_Kabsch_rotation(cbct_points, ct_points)
@@ -1058,7 +1160,10 @@ 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)
-                            
+                        end_rotation = time.time()
+
+                        
+                        start_translation = time.time()
                         fine_shift = np.zeros(3)  # Inicializiraj fine premike
                         if applyTranslation:
                             chosen_translation_vector, cbct_points_transformed, method_used = choose_best_translation(
@@ -1075,31 +1180,52 @@ class MyTransformModuleLogic(ScriptedLoadableModuleLogic):
                             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")
-
-
+                        end_translation = time.time()
 
+                        start_transform = time.time()
                         # ✅ 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
+                        vtk_transform = create_vtk_transform(chosen_rotation_matrix, chosen_translation_vector, studyName, cbct_volume_name, scaling_factors)
+                        
+                        
+                        
+                        
+                        
+                        combined_matrix = np.eye(4)
+                        
+                        if scaling_factors is not None:
+                            scaling_matrix = np.diag([scaling_factors[0], scaling_factors[1], scaling_factors[2], 1.0])
+                            combined_matrix = np.dot(scaling_matrix, combined_matrix)
+                        
+                        if chosen_translation_vector is not None:
+                            combined_matrix[:3, 3] = chosen_translation_vector + transvector # glavna poravnava, sistematični y shift in groba poravnava
+                        if chosen_rotation_matrix is not None:
+                            combined_matrix[:3, :3] = chosen_rotation_matrix
+                        
+                        
+                        cumulative_matrices[(studyName, cbct_volume_name)] = np.dot(cumulative_matrices[(studyName, cbct_volume_name)], combined_matrix)
+
+
 
                         # 🔄 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 = find_volume_node_by_partial_name(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)
+                        
+                        end_transform = time.time()
 
                         # 📍 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)
+                        #print("Markerji po transformaciji:\n", cbct_points, ct_points)
                         
                         #popravek v x osi
                         delta_x_list = [ct[0] - cbct[0] for ct, cbct in zip(ct_points, cbct_points)]
@@ -1114,7 +1240,12 @@ class MyTransformModuleLogic(ScriptedLoadableModuleLogic):
                         # 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
-                                                
+                        
+                        if fine_shift is not None:
+                            shift_matrix = np.eye(4)
+                            shift_matrix[:3, 3] = fine_shift
+                            cumulative_matrices[(studyName, cbct_volume_name)] = np.dot(cumulative_matrices[(studyName, cbct_volume_name)], shift_matrix)
+                 
                     
                         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
@@ -1124,7 +1255,7 @@ class MyTransformModuleLogic(ScriptedLoadableModuleLogic):
                         transform_node = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLTransformNode", imeTransformNoda)
                         transform_node.SetAndObserveTransformToParent(vtk_transform)
 
-                        cbct_volume_node = slicer.util.getNode(cbct_volume_name)
+                        cbct_volume_node = find_volume_node_by_partial_name(cbct_volume_name)
                         cbct_volume_node.SetAndObserveTransformNodeID(transform_node.GetID())
 
                         # 🔨 Uporabi (ali shrani transformacijo kasneje)
@@ -1134,27 +1265,45 @@ class MyTransformModuleLogic(ScriptedLoadableModuleLogic):
                         
                         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")
+                        print(f"Fine correction shifts: ΔX={fine_shift[0]:.2f} mm, ΔY={fine_shift[1]:.2f} mm, ΔZ={fine_shift[2]:.2f} mm")
 
+                        #shrani transformacijsko matriko v datoteko
+                        save_transform_matrix(cumulative_matrices[(studyName, cbct_volume_name)], studyName, cbct_volume_name, ct_table_found)
+                        
+                        
                         # 📏 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))
+                        print("Total Individual errors:", errors)
+                        print("Average error: {:.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")
+                            print(f"Specific marker errors {i+1}: ΔX={diff[0]:.2f} mm, ΔY={diff[1]:.2f} mm, ΔZ={diff[2]:.2f} mm")
                         
                         
-
+                
             else:
                 print(f"Study {studyItem} doesn't have any appropriate CT or CBCT volumes.")
+                continue
+            
+            study_end_time = time.time()    
+            timing_data = {
+                "IO": end_io - start_io,
+                "Fixing": fixing_end - fixing,
+                "Scaling": end_scaling - start_scaling,
+                "CentroidAlign": end_align - start_align,
+                "Rotation": end_rotation - start_rotation,
+                "Translation": end_translation - start_translation,
+                "Transform": end_transform - start_transform,
+                "Total": study_end_time - study_start_time}
+            update_timing_csv(timing_data, studyName)
+            print(f"Timing data for {studyName}: {timing_data}")
 
         # Izpis globalne statistike
         
-        
+        start_save = time.time()
         if writefilecheck:
             #print("Distances between CT & CBCT markers: ", ctcbct_distance)
             #print("Distances between pairs of markers for each volume: ", prostate_size_est)
@@ -1172,13 +1321,14 @@ class MyTransformModuleLogic(ScriptedLoadableModuleLogic):
             #print("Prostate size file written at ", prostate_size_file)
 
             # Write CT-CBCT distance data
-            with open(ctcbct_distance_file, mode='w', newline='') as file:
+            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_save = time.time()
+        print(f"Saving time: {end_save - start_save:.2f} seconds")
         end_time = time.time()
         # Calculate and print elapsed time
         elapsed_time = end_time - start_time
@@ -1186,4 +1336,5 @@ class MyTransformModuleLogic(ScriptedLoadableModuleLogic):
         minutes = int(elapsed_time // 60)
         seconds = elapsed_time % 60
 
-        print(f"Execution time: {minutes} minutes and {seconds:.6f} seconds")
+        print(f"Execution time: {minutes} minutes and {seconds:.6f} seconds")
+