AvtomatskaPoravnavaMarkerjev/SeekTransformModule/SeekTransformModule.py

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 slicer.util
import itertools
import DICOMLib
from DICOMLib import DICOMUtils
import DicomRtImportExportPlugin
from collections import deque, Counter
import vtk
from slicer.ScriptedLoadableModule import *
import qt
import matplotlib.pyplot as plt
import csv
import time
import logging
import matplotlib
matplotlib.use('Agg')  # << to dodaš ZGORAJ, da omogoči PNG zapis brez GUI
from mpl_toolkits.mplot3d import Axes3D  


#exec(open("C:/Users/lkomar/Documents/Prostata/FirstTryRegister.py").read())

cumulative_matrices = {}

class SeekTransformModule(ScriptedLoadableModule):
    """
    Module description shown in the module panel.
    """
    def __init__(self, parent):
        ScriptedLoadableModule.__init__(self, parent)
        self.parent.title = "Seek Transform module"
        self.parent.categories = ["Image Processing"]
        self.parent.contributors = ["Luka Komar (Onkološki Inštitut Ljubljana, Fakulteta za Matematiko in Fiziko Ljubljana)"]
        self.parent.helpText = "This module applies rigid transformations to CBCT volumes based on reference CT volumes."
        self.parent.acknowledgementText = "Supported by doc. Primož Peterlin & prof. Andrej Studen"

class SeekTransformModuleWidget(ScriptedLoadableModuleWidget):
    """
    GUI of the module.
    """
    def setup(self):
        ScriptedLoadableModuleWidget.setup(self)

        # Dropdown menu za izbiro metode
        self.rotationMethodComboBox = qt.QComboBox()
        self.rotationMethodComboBox.addItems(["Kabsch", "Horn", "Iterative Closest Point (Horn)"])
        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)

        self.translationCheckBox = qt.QCheckBox("Translation")
        self.translationCheckBox.setChecked(True)
        self.layout.addWidget(self.translationCheckBox)

        self.markersCheckBox = qt.QCheckBox("Place control points for detected markers")
        self.markersCheckBox.setChecked(False)
        self.layout.addWidget(self.markersCheckBox)
        
        self.writefileCheckBox = qt.QCheckBox("Write data to csv file")
        self.writefileCheckBox.setChecked(True)
        self.layout.addWidget(self.writefileCheckBox)
        
        self.tableCheckBox = qt.QCheckBox("Find top of the table and match height")
        self.tableCheckBox.setChecked(True)
        self.layout.addWidget(self.tableCheckBox)
        
        self.DicomCheckBox = qt.QCheckBox("Save transformed DICOM")
        self.DicomCheckBox.setChecked(True)
        self.layout.addWidget(self.DicomCheckBox)


        # Load button
        self.applyButton = qt.QPushButton("Find markers and transform")
        self.applyButton.toolTip = "Finds markers, computes optimal rigid transform and applies it to CBCT volumes."
        self.applyButton.enabled = True
        self.layout.addWidget(self.applyButton)

        # Connect button to logic
        self.applyButton.connect('clicked(bool)', self.onApplyButton)

        self.layout.addStretch(1)

    def onApplyButton(self):
                # Nastavi globalni logger
        log_file_path = os.path.join("C:/Users/lkomar/Documents/Prostata", "seektransform_log.txt")
        logging.basicConfig(filename=log_file_path, level=logging.INFO, format='%(asctime)s - %(message)s')
        
        try:
            logging.info("▶️ onApplyButton pressed.")
        except Exception as e:
            print("❌ Logging setup failed:", e)
        
        logic = MyTransformModuleLogic()
        selectedMethod = self.rotationMethodComboBox.currentText  # izberi metodo izračuna rotacije

        # Preberi stanje checkboxov
        applyRotation = self.rotationCheckBox.isChecked()
        applyTranslation = self.translationCheckBox.isChecked()
        applyScaling = self.scalingCheckBox.isChecked()
        applyMarkers = self.markersCheckBox.isChecked()
        writefilecheck = self.writefileCheckBox.isChecked()
        tablefind = self.tableCheckBox.isChecked()
        saveasdicom = self.DicomCheckBox.isChecked()

        # Pokliči logiko z izbranimi nastavitvami
        logic.run(selectedMethod, applyRotation, applyTranslation, applyScaling, applyMarkers, writefilecheck, tablefind, saveasdicom)


class MyTransformModuleLogic(ScriptedLoadableModuleLogic):
    """
    Core logic of the module.
    """
    
    
    def run(self, selectedMethod, applyRotation, applyTranslation, applyScaling, applymarkers, writefilecheck, tablefind, save_as_dicom):
        start_time = time.time()
        print("Calculating...")
        #slicer.util.delayDisplay(f"Starting", 1000)

        
        
        def group_points(points, threshold):
            # Function to group points that are close to each other
            grouped_points = []
            while points:
                point = points.pop()  # Take one point from the list
                group = [point]  # Start a new group
                
                # Find all points close to this one
                distances = cdist([point], points)  # Calculate distances from this point to others
                close_points = [i for i, dist in enumerate(distances[0]) if dist < threshold]
                
                # Add the close points to the group
                group.extend([points[i] for i in close_points])
                
                # Remove the grouped points from the list
                points = [point for i, point in enumerate(points) if i not in close_points]
                
                # Add the group to the result
                grouped_points.append(group)
            
            return grouped_points

        def region_growing(image_data, seed, intensity_threshold, max_distance):
            dimensions = image_data.GetDimensions()
            visited = set()
            region = []
            queue = deque([seed])

            while queue:
                x, y, z = queue.popleft()
                if (x, y, z) in visited:
                    continue

                visited.add((x, y, z))
                voxel_value = image_data.GetScalarComponentAsDouble(x, y, z, 0)
                
                if voxel_value >= intensity_threshold:
                    region.append((x, y, z))
                    # Add neighbors within bounds
                    for dx, dy, dz in [(1, 0, 0), (-1, 0, 0), (0, 1, 0), (0, -1, 0), (0, 0, 1), (0, 0, -1)]:
                        nx, ny, nz = x + dx, y + dy, z + dz
                        if 0 <= nx < dimensions[0] and 0 <= ny < dimensions[1] and 0 <= nz < dimensions[2]:
                            if (nx, ny, nz) not in visited:
                                queue.append((nx, ny, nz))

            return region

        def compute_scaling_stddev(moving_points, fixed_points):
            moving = np.array(moving_points)
            fixed = np.array(fixed_points)

            # 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.
                scaling_factors (tuple): Scaling factors (sx, sy, sz) for each axis.

            Returns:
                np.ndarray: Scaled CBCT points.
            """
            sx, sy, sz = scaling_factors  # Extract scaling factors
            scaling_matrix = np.diag([sx, sy, sz])  # Create diagonal scaling matrix

            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

        def compute_Kabsch_rotation(moving_points, fixed_points):
            """
            Computes the optimal rotation matrix to align moving_points to fixed_points.
            
            Parameters:
            moving_points (list or ndarray): List of points to be rotated CBCT
            fixed_points (list or ndarray): List of reference points CT

            Returns:
            ndarray: Optimal rotation matrix.
            """
            assert len(moving_points) == len(fixed_points), "Point lists must be the same length."

            # Convert to numpy arrays
            moving = np.array(moving_points)
            fixed = np.array(fixed_points)

            # Compute centroids
            centroid_moving = np.mean(moving, axis=0)
            centroid_fixed = np.mean(fixed, axis=0)

            # Center the points
            moving_centered = moving - centroid_moving
            fixed_centered = fixed - centroid_fixed

            # Compute covariance matrix
            H = np.dot(moving_centered.T, fixed_centered)

            # SVD decomposition
            U, _, Vt = np.linalg.svd(H)
            Rotate_optimal = np.dot(Vt.T, U.T)

            # Correct improper rotation (reflection)
            if np.linalg.det(Rotate_optimal) < 0:
                Vt[-1, :] *= -1
                Rotate_optimal = np.dot(Vt.T, U.T)

            return Rotate_optimal

        def compute_Horn_rotation(moving_points, fixed_points):
            """
            Computes the optimal rotation matrix using quaternions.

            Parameters:
            moving_points (list or ndarray): List of points to be rotated.
            fixed_points (list or ndarray): List of reference points.

            Returns:
            ndarray: Optimal rotation matrix.
            """
            assert len(moving_points) == len(fixed_points), "Point lists must be the same length."
            
            moving = np.array(moving_points)
            fixed = np.array(fixed_points)
            
            # Compute centroids
            centroid_moving = np.mean(moving, axis=0)
            centroid_fixed = np.mean(fixed, axis=0)
            
            # Center the points
            moving_centered = moving - centroid_moving
            fixed_centered = fixed - centroid_fixed
            
            # Construct the cross-dispersion matrix
            M = np.dot(moving_centered.T, fixed_centered)
            
            # Construct the N matrix for quaternion solution
            A = M - M.T
            delta = np.array([A[1, 2], A[2, 0], A[0, 1]])
            trace = np.trace(M)
            
            N = np.zeros((4, 4))
            N[0, 0] = trace
            N[1:, 0] = delta
            N[0, 1:] = delta
            N[1:, 1:] = M + M.T - np.eye(3) * trace
            
            # Compute the eigenvector corresponding to the maximum eigenvalue
            eigvals, eigvecs = np.linalg.eigh(N)
            q_optimal = eigvecs[:, np.argmax(eigvals)]  # Optimal quaternion
            
            # Convert quaternion to rotation matrix
            w, x, y, z = q_optimal
            R = np.array([
                [1 - 2*(y**2 + z**2), 2*(x*y - z*w), 2*(x*z + y*w)],
                [2*(x*y + z*w), 1 - 2*(x**2 + z**2), 2*(y*z - x*w)],
                [2*(x*z - y*w), 2*(y*z + x*w), 1 - 2*(x**2 + y**2)]
            ])
            
            return R

        def icp_algorithm(moving_points, fixed_points, max_iterations=100, tolerance=1e-5):
            """
            Iterative Closest Point (ICP) algorithm to align moving_points to fixed_points.
            
            Parameters:
            moving_points (list or ndarray): List of points to be aligned.
            fixed_points (list or ndarray): List of reference points.
            max_iterations (int): Maximum number of iterations.
            tolerance (float): Convergence tolerance.

            Returns:
            ndarray: Transformed moving points.
            ndarray: Optimal rotation matrix.
            ndarray: Optimal translation vector.
            """
            # Convert to numpy arrays
            moving = np.array(moving_points)
            fixed = np.array(fixed_points)

            # Initialize transformation
            R = np.eye(3)  # Identity matrix for rotation
            t = np.zeros(3)  # Zero vector for translation

            prev_error = np.inf  # Initialize previous error to a large value

            for iteration in range(max_iterations):
                # Step 1: Find the nearest neighbors (correspondences)
                distances = np.linalg.norm(moving[:, np.newaxis] - fixed, axis=2)
                nearest_indices = np.argmin(distances, axis=1)
                nearest_points = fixed[nearest_indices]

                # Step 2: Compute the optimal rotation and translation
                R_new = compute_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)

                # Step 3: Apply the transformation
                moving = np.dot(moving, R_new.T) + t_new

                # Update the cumulative transformation
                R = np.dot(R_new, R)
                t = np.dot(R_new, t) + t_new

                # Step 4: Check for convergence
                mean_error = np.mean(np.linalg.norm(moving - nearest_points, axis=1))
                if np.abs(prev_error - mean_error) < tolerance:
                    print(f"ICP converged after {iteration + 1} iterations.")
                    break
                prev_error = mean_error

            else:
                print(f"ICP reached maximum iterations ({max_iterations}).")

            return moving, R, t

        def match_points(cbct_points, ct_points, auto_weights=False, fallback_if_worse=False, normalize_lengths=True, normalize_angles=False, min_distance=5, w_order=1.0):

            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])
                ]
                return lengths

            def triangle_angles(points):
                a = np.linalg.norm(points[1] - points[2])
                b = np.linalg.norm(points[0] - points[2])
                c = np.linalg.norm(points[0] - points[1])
                angle_A = np.arccos(np.clip((b**2 + c**2 - a**2) / (2 * b * c), -1.0, 1.0))
                angle_B = np.arccos(np.clip((a**2 + c**2 - b**2) / (2 * a * c), -1.0, 1.0))
                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, penalty_angle_thresh=np.deg2rad(10)):
                perm_lengths = side_lengths(perm)
                perm_angles = triangle_angles(perm)

                # 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

                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))
                
                order_penalty = order_mismatch_penalty(ct_points, perm, axis='z')
                return w_len * score_len + w_ang * score_ang + order_penalty + w_order * order_penalty
                
            def smart_sort_cbct_points(cbct_points, z_threshold=5.0):
                """
                Sortira točke tako, da poskusi najprej po Z. Če so razlike po Z manjše
                od praga, sortira po (Y, X), sicer sortira po Z.
                """
                z_values = [pt[2] for pt in cbct_points]
                z_range = max(z_values) - min(z_values)

                if z_range < z_threshold:
                    # Sortiraj po Y, nato X (če so točke v isti ravnini po Z)
                    return sorted(cbct_points, key=lambda pt: (pt[1], pt[0]))
                else:
                    # Sortiraj po Z, nato Y, nato X
                    return sorted(cbct_points, key=lambda pt: (pt[2], pt[1], pt[0]))
            
            def order_mismatch_penalty(ct_points, perm, axis='z'):
                axis_idx = {'x': 0, 'y': 1, 'z': 2}[axis]
                ct_sorted = np.argsort([pt[axis_idx] for pt in ct_points])  
                perm_sorted = np.argsort([pt[axis_idx] for pt in perm])
                return sum(1 for a, b in zip(ct_sorted, perm_sorted) if a != b)
            
            
            cbct_points = list(cbct_points)
            print("CBCT points:", cbct_points)
            ct_lengths = side_lengths(np.array(ct_points))
            ct_angles = triangle_angles(np.array(ct_points))

            if auto_weights:
                var_len = np.var(ct_lengths)
                var_ang = np.var(ct_angles)
                total_var = var_len + var_ang + 1e-6
                weight_length = (1 - var_len / total_var)
                weight_angle = (1 - var_ang / total_var)
            else:
                weight_length = 0.8
                weight_angle = 0.2

            cbct_sorted = smart_sort_cbct_points(cbct_points)
            original_score = permutation_score(np.array(cbct_sorted), ct_lengths, ct_angles, weight_length, weight_angle)

            # Če je ta rezultat dovolj dober, uporabi
            best_score = float('inf')
            best_perm = None

            if original_score < float('inf'):  # lahko dodaš prag če želiš
                best_score = original_score
                best_perm = np.array(cbct_sorted)

            # Nato preveri vse permutacije (vključno s prvotnim vrstnim redom, če fallback_if_worse=True)
            for perm in itertools.permutations(cbct_points):
                perm = np.array(perm)
                score = permutation_score(perm, ct_lengths, ct_angles, weight_length, weight_angle)
                if score < best_score:
                    best_score = score
                    best_perm = perm
                    print(f"New best permutation found with 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)
                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.
            
            Parameters:
            moving_points (list or ndarray): List of points to be translated.
            fixed_points (list or ndarray): List of reference points.
            rotation_matrix (ndarray): Rotation matrix.

            Returns:
            ndarray: Translation vector.
            """
            # Convert to numpy arrays
            moving = np.array(moving_points)
            fixed = np.array(fixed_points)

            # Compute centroids
            centroid_moving = np.mean(moving, axis=0)
            centroid_fixed = np.mean(fixed, axis=0)

            # Compute translation
            translation = centroid_fixed - np.dot(centroid_moving, rotation_matrix)

            return translation

        def create_vtk_transform(rotation_matrix, translation_vector, tablefound, study_name=None, cbct_volume_name=None, scaling_factors=None):
            """
            Creates a vtkTransform from scaling, rotation, and translation.
            Shrani tudi kumulativno matriko v globalni slovar cumulative_matrices.
            """
            # ----- Inicializacija -----
            global cumulative_matrices
            transform = vtk.vtkTransform()

            # ----- 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, 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)))

            # ----- 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):
            """
            Appends the given 4x4 matrix to a text file under the given study folder.
            """
            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
            safe_cbct_name = re.sub(r'[<>:"/\\|?*]', '_', cbct_volume_name)
            # Preveri ali je CT miza najdena
            filename = os.path.join(study_folder, f"{safe_cbct_name}.txt")
            with open(filename, "w") 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
            #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 = find_volume_node_by_partial_name(volume_name)
            if not volume_node:
                raise RuntimeError(f"Volume {volume_name} not found.")
            
            image_data = volume_node.GetImageData()
            matrix = vtk.vtkMatrix4x4()
            volume_node.GetIJKToRASMatrix(matrix)

            dimensions = image_data.GetDimensions()
            #detected_regions = []

            if yesCbct: #je cbct ali ct?
                valid_x_min, valid_x_max = 0, dimensions[0] - 1
                valid_y_min, valid_y_max = 0, dimensions[1] - 1
                valid_z_min, valid_z_max = 0, dimensions[2] - 1
            else:
                valid_x_min, valid_x_max = max(x_min, 0), min(x_max, dimensions[0] - 1)
                valid_y_min, valid_y_max = max(y_min, 0), min(y_max, dimensions[1] - 1)
                valid_z_min, valid_z_max = max(z_min, 0), min(z_max, dimensions[2] - 1)

            visited = set()

            def grow_region(x, y, z):
                if (x, y, z) in visited:
                    return None

                voxel_value = image_data.GetScalarComponentAsDouble(x, y, z, 0)
                if voxel_value < intensity_threshold:
                    return None

                region = region_growing(image_data, (x, y, z), intensity_threshold, max_distance=max_distance)
                if region:
                    for point in region:
                        visited.add(tuple(point))
                    return region
                return None

            regions = []
            for z in range(valid_z_min, valid_z_max + 1):
                for y in range(valid_y_min, valid_y_max + 1):
                    for x in range(valid_x_min, valid_x_max + 1):
                        region = grow_region(x, y, z)
                        if region:
                            regions.append(region)

            # Collect centroids using intensity-weighted average
            centroids = []
            for region in regions:
                points = np.array([matrix.MultiplyPoint([*point, 1])[:3] for point in region])
                intensities = np.array([image_data.GetScalarComponentAsDouble(*point, 0) for point in region])
                
                if intensities.sum() > 0:
                    weighted_centroid = np.average(points, axis=0, weights=intensities)
                    max_intensity = intensities.max()
                    centroids.append((np.round(weighted_centroid, 2), max_intensity))

            unique_centroids = []
            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))
            
            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

        def find_table_top_z(ct_volume_name, writefilecheck, makemarkerscheck, yesCbct):
            """
            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 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)

            # --- 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

            # --- 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)}/11")

            
            """ # --- 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(f"⚠️ Rob mize ni bil najden (X = {x_index})")
                print("Column values (partial):", column_values.tolist())
                return None

            # --- Pretvorba v RAS koordinato ---
            table_ijk = [x_index, table_top_y, z_index]
            table_ras = np.array(ijkToRasMatrix.MultiplyPoint([*table_ijk, 1]))[:3]

            # --- 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 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 table edge at Z = {table_ras[1]:.2f} mm"])

            return table_ras[1], table_top_y

        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

                #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)
                
                # 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
    
        def print_orientation(volume_name):
            node = find_volume_node_by_partial_name(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)])

        def prealign_by_centroid(cbct_points, ct_points):
            """
            Predporavna CBCT markerje na CT markerje glede na centrične točke.

            Args:
                cbct_points: List ali ndarray točk iz CBCT.
                ct_points: List ali ndarray točk iz CT.

            Returns:
                List: CBCT točke premaknjene tako, da so centrične točke usklajene.
            """
            cbct_points = np.array(cbct_points)
            ct_points = np.array(ct_points)
            cbct_centroid = np.mean(cbct_points, axis=0)
            ct_centroid = np.mean(ct_points, axis=0)
            translation_vector = ct_centroid - cbct_centroid
            aligned_cbct = cbct_points + translation_vector
            return aligned_cbct, translation_vector
        
        def choose_best_translation(cbct_points, ct_points, rotation_matrix):
            """
            Izbere boljšo translacijo: centroidno ali povprečno po rotaciji (retranslation).
            
            Args:
                cbct_points (array-like): Točke iz CBCT (še ne rotirane).
                ct_points (array-like): Ciljne CT točke.
                rotation_matrix (ndarray): Rotacijska matrika.

            Returns:
                tuple: (best_translation_vector, transformed_cbct_points, used_method)
            """
            cbct_points = np.array(cbct_points)
            ct_points = np.array(ct_points)
            
            # 1. Rotiraj CBCT točke
            rotated_cbct = np.dot(cbct_points, rotation_matrix.T)
            
            # 2. Centroid translacija
            centroid_moving = np.mean(cbct_points, axis=0)
            centroid_fixed = np.mean(ct_points, axis=0)
            translation_centroid = centroid_fixed - np.dot(centroid_moving, rotation_matrix)
            transformed_centroid = rotated_cbct + translation_centroid
            error_centroid = np.mean(np.linalg.norm(transformed_centroid - ct_points, axis=1))

            # 3. Retranslacija (srednja razlika)
            translation_recomputed = np.mean(ct_points - rotated_cbct, axis=0)
            transformed_recomputed = rotated_cbct + translation_recomputed
            error_recomputed = np.mean(np.linalg.norm(transformed_recomputed - ct_points, axis=1))

            # 4. Izberi boljšo
            if error_recomputed < error_centroid:
                #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)")
                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!"

            # Ustvari markups za CBCT
            cbct_node = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLMarkupsFiducialNode", f"{study_name}_CBCT")
            cbct_node.GetDisplayNode().SetSelectedColor(0, 0, 1)  # modra

            # Ustvari markups za CT
            ct_node = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLMarkupsFiducialNode", f"{study_name}_CT")
            ct_node.GetDisplayNode().SetSelectedColor(1, 0, 0)  # rdeča

            # Dodaj točke
            for i, (cbct, ct) in enumerate(zip(cbct_points, ct_points)):
                cbct_node.AddControlPoint(*cbct, f"CBCT_{i}")
                ct_node.AddControlPoint(*ct, f"CT_{i}")

            # Ustvari model z linijami med pari
            points = vtk.vtkPoints()
            lines = vtk.vtkCellArray()

            for i, (p1, p2) in enumerate(zip(cbct_points, ct_points)):
                id1 = points.InsertNextPoint(p1)
                id2 = points.InsertNextPoint(p2)

                line = vtk.vtkLine()
                line.GetPointIds().SetId(0, id1)
                line.GetPointIds().SetId(1, id2)
                lines.InsertNextCell(line)

            polyData = vtk.vtkPolyData()
            polyData.SetPoints(points)
            polyData.SetLines(lines)

            # Model node
            modelNode = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLModelNode", f"{study_name}_Connections")
            modelNode.SetAndObservePolyData(polyData)

            modelDisplay = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLModelDisplayNode")
            modelDisplay.SetColor(0, 0, 0)  # črna
            modelDisplay.SetLineWidth(2)
            modelDisplay.SetVisibility(True)

            modelNode.SetAndObserveDisplayNodeID(modelDisplay.GetID())
            modelNode.SetAndObservePolyData(polyData)

            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).
            """
            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
       
        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", "Table", "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.")
                
        def convert_rtstruct_to_segmentation_nodes(ct_volume_node):
            shNode = slicer.vtkMRMLSubjectHierarchyNode.GetSubjectHierarchyNode(slicer.mrmlScene)
            segmentation_nodes = []

            # ✅ Najdi vse SegmentationNode-e, ki imajo segmente
            for seg_node in slicer.util.getNodesByClass("vtkMRMLSegmentationNode"):
                num_segments = seg_node.GetSegmentation().GetNumberOfSegments()
                if num_segments == 0:
                    continue

                # Če še nima referenceVolume, jo nastavimo
                if not seg_node.GetNodeReferenceID("referenceVolume"):
                    seg_node.SetReferenceImageGeometryParameterFromVolumeNode(ct_volume_node)
                    seg_node.SetNodeReferenceID("referenceVolume", ct_volume_node.GetID())
                    print(f"🔗 Nastavljena referenca na CT za: {seg_node.GetName()}")

                segmentation_nodes.append(seg_node)
                print(f"📎 Najdena segmentacija z vsebino: {seg_node.GetName()}, segmentov: {num_segments}")

            if segmentation_nodes:
                return segmentation_nodes

            for item_id in range(shNode.GetNumberOfItems()):
                modality = shNode.GetItemAttribute(item_id, "DICOM.Modality")
                if modality != "RTSTRUCT":
                    continue

                rtstruct_node = shNode.GetItemDataNode(item_id)
                if not rtstruct_node:
                    continue

                print(f"📎 Najden RTSTRUCT: {rtstruct_node.GetName()}")

                # Ustvari nov SegmentationNode
                seg_node = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLSegmentationNode", f"Seg_{rtstruct_node.GetName()}")
                seg_node.SetReferenceImageGeometryParameterFromVolumeNode(ct_volume_node)
                seg_node.SetNodeReferenceID("referenceVolume", ct_volume_node.GetID())



                # Najdi child elemente v SubjectHierarchy (strukture)
                rtstruct_item_id = shNode.GetItemByDataNode(rtstruct_node)
                structure_ids = vtk.vtkIdList()
                shNode.GetItemChildren(rtstruct_item_id, structure_ids)

                segment_count = 0

                for i in range(structure_ids.GetNumberOfIds()):
                    structure_item_id = structure_ids.GetId(i)
                    name = shNode.GetItemName(structure_item_id)
                    associated_node = shNode.GetItemDataNode(structure_item_id)
                    if associated_node and associated_node.IsA("vtkMRMLModelNode"):
                        print(f" ➕ Dodajam strukturo: {name}")
                        slicer.modules.segmentations.logic().ImportModelToSegmentationNode(associated_node, seg_node)
                        seg_node.GetSegmentation().GetSegment(seg_node.GetSegmentation().GetNumberOfSegments() - 1).SetName(name)
                        segment_count += 1

                # Če ni bilo nič dodano iz SH: poskusi uvoziti modele iz scene
                if segment_count == 0:
                    print("⚠️ RTSTRUCT nima struktur v SubjectHierarchy – poskus uvoza vseh modelov iz scene.")
                    for model_node in slicer.util.getNodesByClass("vtkMRMLModelNode"):
                        if "RTSTRUCT" in model_node.GetName().upper() or model_node.GetName().startswith("Model"):
                            print(f" ➕ [fallback] Uvoz modela: {model_node.GetName()}")
                            slicer.modules.segmentations.logic().ImportModelToSegmentationNode(model_node, seg_node)
                            seg_node.GetSegmentation().GetSegment(seg_node.GetSegmentation().GetNumberOfSegments() - 1).SetName(model_node.GetName())
                
                print(f"📊 Segmentov v {seg_node.GetName()}: {seg_node.GetSegmentation().GetNumberOfSegments()}")
                segmentation_nodes.append(seg_node)

            return segmentation_nodes

        def apply_cumulative_transform_to_segmentation(segmentation_node, matrix):
            transform = vtk.vtkTransform()
            vtk_matrix = vtk.vtkMatrix4x4()
            for i in range(4):
                for j in range(4):
                    vtk_matrix.SetElement(i, j, matrix[i, j])
            transform.SetMatrix(vtk_matrix)

            transform_node = slicer.vtkMRMLTransformNode()
            slicer.mrmlScene.AddNode(transform_node)
            transform_node.SetAndObserveTransformToParent(transform)
            segmentation_node.SetAndObserveTransformNodeID(transform_node.GetID())

            slicer.vtkSlicerTransformLogic().hardenTransform(segmentation_node)
            slicer.mrmlScene.RemoveNode(transform_node)

        def convert_all_models_to_segmentation(reference_volume_name: str, prefix: str = "Imported_"):
            """
            Pretvori vse modele (ModelNode) v sceni v enoten vtkMRMLSegmentationNode.

            📥 VHODI:
            ----------
            reference_volume_name : str
                Ime obstoječega CT volumna (npr. "CT_1"), ki določa geometrijo za segmentacijo.
                Ta volumen mora biti že naložen v sceni.

            prefix : str
                Predpona za ime novega segmentation noda (npr. "Imported_").
                Ime novega noda bo nekaj kot: "Imported_Segmentation".

            📤 IZHOD:
            ----------
            segmentation_node : vtkMRMLSegmentationNode
                Nov nod, ki vsebuje en segment za vsak najden model v sceni.
                Ta segmentacijski nod je pripravljen za transformacijo in DICOM export (vsebuje BinaryLabelmap).
            """
            import slicer

            # Pridobi referenčni volumen (CT)
            reference_volume = slicer.util.getNode(reference_volume_name)

            # Ustvari nov segmentacijski nod
            segmentation_node = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLSegmentationNode", prefix + "Segmentation")
            segmentation_node.SetReferenceImageGeometryParameterFromVolumeNode(reference_volume)
            segmentation_node.SetNodeReferenceID("referenceVolume", reference_volume.GetID())

            # Najdi vse modele
            model_nodes = slicer.util.getNodesByClass("vtkMRMLModelNode")
            #print(f"📦 Najdenih modelov: {len(model_nodes)}")
            
            for model_node in model_nodes:
                name = model_node.GetName()
                #print(f"🔍 Model: {name}")
                skip = (
                    name.lower().startswith("segmentation")
                    or name.lower().startswith("surface")
                    or name.lower() in ["red volume slice", "green volume slice", "yellow volume slice"]
                    or "rtstruct" in name.lower()
                )
                if skip:
                    continue


                # Uvozi model kot segment
                success = slicer.modules.segmentations.logic().ImportModelToSegmentationNode(model_node, segmentation_node)
                if not success:
                    #print(f"✅ Model '{name}' uvožen kot segment.")
                    print(f"❌ Napaka pri uvozu modela: {name}")
                    

            # Ustvari BinaryLabelmap reprezentacijo (nujno za DICOM export)
            created = segmentation_node.GetSegmentation().CreateRepresentation("BinaryLabelmap")
            if created:
                print("✅ BinaryLabelmap reprezentacija uspešno ustvarjena.")
                segmentation_node.GetSegmentation().SetMasterRepresentationName("BinaryLabelmap")
            #else:
                #print("❌ Pretvorba v BinaryLabelmap ni uspela.")

            return segmentation_node

        def apply_cumulative_transform_to_volume(volume_node, matrix):
            transform = vtk.vtkTransform()
            vtk_matrix = vtk.vtkMatrix4x4()
            for i in range(4):
                for j in range(4):
                    vtk_matrix.SetElement(i, j, matrix[i, j])
            transform.SetMatrix(vtk_matrix)

            transform_node = slicer.vtkMRMLTransformNode()
            slicer.mrmlScene.AddNode(transform_node)
            transform_node.SetAndObserveTransformToParent(transform)
            volume_node.SetAndObserveTransformNodeID(transform_node.GetID())

            slicer.vtkSlicerTransformLogic().hardenTransform(volume_node)
            slicer.mrmlScene.RemoveNode(transform_node)

        def triangle_similarity(p1, p2):
            """
            Primerja dva trikotnika (vsak definiran s 3 točkami v 3D) na podlagi:
            - razmerij dolžin
            - razlik med koti

            :param p1: seznam treh točk (npr. ct_points)
            :param p2: seznam treh točk (npr. cbct_points)
            :return: dict z razliko dolžin, razliko kotov, povprečno napako
            """
            def side_lengths(pts):
                a = np.linalg.norm(pts[1] - pts[0])
                b = np.linalg.norm(pts[2] - pts[1])
                c = np.linalg.norm(pts[0] - pts[2])
                return np.array([a, b, c])

            def angles(pts):
                # uporabimo kosinusni izrek
                a, b, c = side_lengths(pts)
                angle_A = np.arccos((b**2 + c**2 - a**2) / (2*b*c))
                angle_B = np.arccos((a**2 + c**2 - b**2) / (2*a*c))
                angle_C = np.pi - angle_A - angle_B
                return np.degrees([angle_A, angle_B, angle_C])

            l1 = side_lengths(p1)
            l2 = side_lengths(p2)
            a1 = angles(p1)
            a2 = angles(p2)

            length_diff = np.abs(l1 - l2)
            angle_diff = np.abs(a1 - a2)

            return {
                "side_length_diff_mm": length_diff,
                "angle_diff_deg": angle_diff,
                "mean_length_error_mm": np.mean(length_diff),
                "mean_angle_error_deg": np.mean(angle_diff),
                "TSR": 1 / (1 + np.mean(length_diff) + np.mean(angle_diff) / 10)
            }

        def save_triangle_visualization(ct_pts, cbct_pts, outpath):
            fig = plt.figure(figsize=(10, 8))
            ax = fig.add_subplot(111, projection='3d')

            ct_closed = np.vstack([ct_pts, ct_pts[0]])
            cbct_closed = np.vstack([cbct_pts, cbct_pts[0]])

            ax.plot(ct_closed[:, 0], ct_closed[:, 1], ct_closed[:, 2], 'b-', label='CT trikotnik')
            ax.scatter(ct_pts[:, 0], ct_pts[:, 1], ct_pts[:, 2], c='blue')

            ax.plot(cbct_closed[:, 0], cbct_closed[:, 1], cbct_closed[:, 2], 'r--', label='CBCT trikotnik')
            ax.scatter(cbct_pts[:, 0], cbct_pts[:, 1], cbct_pts[:, 2], c='red')

            ct_centroid = np.mean(ct_pts, axis=0)
            cbct_centroid = np.mean(cbct_pts, axis=0)
            ax.scatter(*ct_centroid, c='blue', marker='x', s=60)
            ax.scatter(*cbct_centroid, c='red', marker='x', s=60)

            
            ax.set_title("Trikotnik markerjev: CT (modro) vs CBCT (rdeče)")
            ax.set_xlabel('X (mm)')
            ax.set_ylabel('Y (mm)')
            ax.set_zlabel('Z (mm)')
            ax.legend()
            ax.view_init(elev=20, azim=30)
            plt.tight_layout()

            try:
                fig.savefig(outpath)
                print(f"🖼 Triangle visualization saved to: {outpath}")
            except Exception as e:
                print(f"❌ Failed to save triangle visualization: {e}")
            finally:
                plt.close(fig)
            
        # Globalni seznami za končno statistiko
        prostate_size_est = []
        ctcbct_distance = []
        table_z_values = {}

        # Pridobimo SubjectHierarchyNode
        shNode = slicer.vtkMRMLSubjectHierarchyNode.GetSubjectHierarchyNode(slicer.mrmlScene)
        
        studyItems = vtk.vtkIdList()
        shNode.GetItemChildren(shNode.GetSceneItemID(), studyItems)
        
        #slicer.util.delayDisplay(f"[DEBUG] Starting the loops", 1000)

        
        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}")
            # **LOKALNI** seznami, resetirajo se pri vsakem study-ju
            cbct_list = []
            ct_list = []
            volume_points_dict = {}
            CT_offset = 0

            # Get child items of the study item
            volumeItems = vtk.vtkIdList()
            shNode.GetItemChildren(studyItem, volumeItems)
            
            # Iteracija čez vse volumne v posameznem studyju
            for j in range(volumeItems.GetNumberOfIds()):
                intermediateItem = volumeItems.GetId(j)

                finalVolumeItems = vtk.vtkIdList()
                shNode.GetItemChildren(intermediateItem, finalVolumeItems)  # Išči globlje!

                for k in range(finalVolumeItems.GetNumberOfIds()):
                    volumeItem = finalVolumeItems.GetId(k)
                    volumeNode = shNode.GetItemDataNode(volumeItem)
                    try:
                        dicomUIDs = volumeNode.GetAttribute("DICOM.instanceUIDs")
                    except AttributeError:
                        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:
                        print("❌ This is an NRRD volume!")
                        continue  # Preskoči, če ni DICOM volume             
                    
                                            
                    volumeName = volumeNode.GetName()
                    imageItem = shNode.GetItemByDataNode(volumeNode)
                    modality = shNode.GetItemAttribute(imageItem, "DICOM.Modality")                                 #deluje!

                    #dimensions = volumeNode.GetImageData().GetDimensions()
                    #spacing = volumeNode.GetSpacing()
                    #print(f"Volume {volumeNode.GetName()} - Dimenzije: {dimensions}, Spacing: {spacing}")

                    if modality != "CT":
                        print("Not a CT")
                        continue  # Preskoči, če ni CT

                    # Preveri, če volume obstaja v sceni
                    if not slicer.mrmlScene.IsNodePresent(volumeNode):
                        print(f"Volume {volumeName} not present in the scene.")
                        continue

                    # Preverimo proizvajalca (DICOM metapodatki)
                    manufacturer = shNode.GetItemAttribute(imageItem, 'DICOM.Manufacturer')
                    #manufacturer = volumeNode.GetAttribute("DICOM.Manufacturer")
                    #manufacturer = slicer.dicomDatabase.fileValue(uid, "0008,0070")
                    #print(manufacturer)
                    
                    # Določimo, ali gre za CBCT ali CT
                    if "varian" in manufacturer.lower() or "elekta" in manufacturer.lower():
                        cbct_list.append(volumeName)
                        scan_type = "CBCT"
                        yesCbct = True
                    else:  # Siemens ali Philips
                        ct_list.append(volumeName)
                        scan_type = "CT"
                        yesCbct = False

                    if volumeNode and volumeNode.IsA("vtkMRMLScalarVolumeNode"):
                        print(f"✔️ {scan_type} {volumeNode.GetName()} (ID: {volumeItem})")
                
                    if not volumeNode or not volumeNode.IsA("vtkMRMLScalarVolumeNode"):
                        print("Can't find volumeNode")
                        #continue  # Preskoči, če ni veljaven volume
                        
                        
                    # Detekcija točk v volumnu
                    ustvari_marker = not yesCbct  # Ustvari markerje pred poravnavo na mizo
                    grouped_points = detect_points_region_growing(volumeName, yesCbct, ustvari_marker, intensity_threshold=3000)
                    if grouped_points is None or len(grouped_points) < 3:
                        print(f"⚠️ Volume {volumeName} doesn't have enough points for registration. Points: {len(grouped_points)}")
                        continue
                    if not yesCbct:
                        
                        # loči koordinate in intenzitete
                        coords_only = [pt for pt, _ in grouped_points]
                        intensities = [intensity for _, intensity in grouped_points]

                        # permutiraj koordinate (npr. zaradi boljšega ujemanja)
                        coords_sorted = match_points(coords_only, coords_only)

                        # ponovno sestavi pare (točka, intenziteta)
                        grouped_points = list(zip(coords_sorted, intensities))
                        
                    #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

            # Če imamo oba tipa volumna (CBCT in CT) **znotraj istega studyja**
            end_io = time.time()
            if cbct_list and ct_list:
                fixing = fixing_end = 0
                table1_time = table1end_time = 0
                table2_time = table2end_time = 0
                start_scaling = end_scaling = 0
                start_align = end_align = 0
                start_rotation = end_rotation = 0
                start_translation = end_translation = 0
                start_transform = end_transform = 0
                study_start_time = study_end_time = 0
                
                ct_volume_name = ct_list[0]  # Uporabi prvi CT kot referenco
                ct_volume_Node = find_volume_node_by_partial_name(ct_volume_name)
                print(f"\nProcessing CT: {ct_volume_name}")
                yesCbct = False
                if(tablefind):
                    table1_time = time.time()
                    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)
                    table1end_time = time.time()                                    
                     
                
                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. Points: {len(ct_points)}")
                    continue
                else:
                    # if len(ct_points) == 4:
                    #     ct_points = remove_lowest_marker(ct_points) #odstrani marker v riti, če obstaja
                    
                    
                    for cbct_volume_name in cbct_list:
                        print(f"\nProcessing CBCT Volume: {cbct_volume_name}")                        
                        yesCbct = True
                        scan_type = "CBCT" #redundant but here we are
                        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()
                        if(tablefind):
                            table2_time = 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"CT offset: {CT_offset}, cbct offset: {mm_offset}")
                                    
                                    #ne rabimo več CT_offset, le še cbct offset. 
                                    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
                            table2end_time = time.time()            
                                            
                        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
                        
                        if len(cbct_points) != len(ct_points):
                            if(len(cbct_points) + 1 == len(ct_points)):
                                ct_points = remove_lowest_marker(ct_points) #odstrani marker v riti, če obstaja
                            else:
                                print(f"Neujemajoče število točk! CBCT: {len(cbct_points)}, CT: {len(ct_points)}")
                                continue  
                        
                        # 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}")
                        
                        
                        
                        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
                                                
                                                
                        if len(cbct_points) < 3:
                            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)
                        fixing_end = time.time()
                        #visualize_point_matches_in_slicer(cbct_points, ct_points, studyName) #poveže pare markerjev
        
                        if writefilecheck:
                            # Shranjevanje razdalj
                            distances_ct_cbct = []
                            distances_internal = {"A-B": [], "B-C": [], "C-A": []}

                            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])]
                            
                            # Razdalje med CT in CBCT (SORTIRANE točke!) 
                            if cbct_points_sorted.shape[0] != len(ct_points):
                                print(f"⚠️ Število točk CBCT ({cbct_points_sorted.shape[0]}) != CT ({len(ct_points)}), preskakujem izračun.")
                            else:
                                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])

                            # 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

                        # 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)
                        start_scaling = time.time()
                        scaling_factors = None
                        if applyScaling:
                            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, 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)
                            elif selectedMethod == "Horn":
                                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)
                        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(
                                cbct_points, ct_points, chosen_rotation_matrix)
                            
                            # 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)


                            # 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

                        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, scaling_factors)
                        
                        
                        
                        
                        
                        combined_matrix = np.eye(4)

                        # 1. Rotacija
                        if chosen_rotation_matrix is not None:
                            combined_matrix[:3, :3] = chosen_rotation_matrix

                        # 2. Skaliranje
                        if scaling_factors is not None:
                            # Pomnoži rotacijo s skalirnim faktorjem po vsaki osi
                            scaled_rotation = combined_matrix[:3, :3] * scaling_factors  # broadcasting vsake vrstice
                            combined_matrix[:3, :3] = scaled_rotation

                        # 3. Translacija
                        if chosen_translation_vector is not None:
                            combined_matrix[:3, 3] = chosen_translation_vector + transvector  # združena translacija

                        
                        #preverjanje determinante
                        rot_part = combined_matrix[:3, :3]
                        det = np.linalg.det(rot_part)
                        if not np.isclose(det, 1.0, atol=0.01):
                            print(f"⚠️ Neortogonalna rotacija! Determinanta: {det}")
                        
                        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 = 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)
                        cbct_points = match_points(cbct_points, ct_points)
                        #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
                        
                        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
                        
                        #### TEST ROTACIJA ########
                        angle_deg = 0
                        angle_rad = np.deg2rad(angle_deg)

                        chosen_rotation_matrix = np.array([
                            [np.cos(angle_rad), -np.sin(angle_rad), 0],
                            [np.sin(angle_rad),  np.cos(angle_rad), 0],
                            [0,                 0,                 1]
                        ])
                        ###KONEC TESTA###
                        
                        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 = 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)

                        
                        table_shift_matrix_ct = np.eye(4)
                        table_shift_matrix_ct[1, 3] = CT_offset  # ali skupni_offset, če je potreben simetričen premik
                        
                        cumulative_matrices[(studyName, cbct_volume_name)] = np.dot(
                            table_shift_matrix,
                            cumulative_matrices[(studyName, cbct_volume_name)]
                        )
                        
                        #ustvari_marker = True
                        cbct_points = [centroid for centroid, _ in detect_points_region_growing(cbct_volume_name, yesCbct, applymarkers, 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")

                        #shrani transformacijsko matriko v datoteko
                        save_transform_matrix(cumulative_matrices[(studyName, cbct_volume_name)], studyName, cbct_volume_name)
                        
                        if save_as_dicom:
                            # Apply transform to CT
                            logging.info(f"[{studyName}] Start applying transform to CT volume.")
                            slicer.util.delayDisplay(f"[DEBUG] Start transform on CT", 1000)
                            apply_cumulative_transform_to_volume(ct_volume_Node, cumulative_matrices[(studyName, cbct_volume_name)])
                            
                            # Convert RTSTRUCT to SegmentationNode if needed
                            logging.info(f"[{studyName}] Getting segmentation nodes from RTSTRUCT.")
                            #slicer.util.delayDisplay(f"[DEBUG] Start segmentation conversion", 1000)
                            #seg_nodes = convert_rtstruct_to_segmentation_nodes(ct_volume_Node)
                            seg_node = convert_all_models_to_segmentation(ct_volume_name, prefix="Imported_")
                            #new_seg_node = slicer.util.getNode(new_seg_name)
                            apply_cumulative_transform_to_segmentation(seg_node, cumulative_matrices[(studyName, cbct_volume_name)])

                            plugin = DicomRtImportExportPlugin.DicomRtImportExportPluginClass()
                            shNode = slicer.vtkMRMLSubjectHierarchyNode.GetSubjectHierarchyNode(slicer.mrmlScene)
                            ct_itemID = shNode.GetItemByDataNode(ct_volume_Node)
                            seg_itemID = shNode.GetItemByDataNode(seg_node)

                            
                            
                            cbct_item = shNode.GetItemByDataNode(cbct_volume_node)
                            cbct_date = shNode.GetItemAttribute(cbct_item, "DICOM.SeriesDate") or "unknownDate"

                            #cbct_date = cbct_volume_node.GetAttribute("DICOM.AcquisitionDate") or "unknownDate"
                            results_base = os.path.join(os.path.dirname(__file__), "Rezultati")
                            study_folder = studyName.replace('^', '_')
                            study_dir = os.path.join(results_base, study_folder)
                            os.makedirs(study_dir, exist_ok=True)  

                            export_dir = os.path.join(study_dir, f"{cbct_date}_DICOM")
                            os.makedirs(export_dir, exist_ok=True)
                            
                            exportables = []
                            exportables += plugin.examineForExport(ct_itemID)
                            exportables += plugin.examineForExport(seg_itemID)
                            for exp in exportables:
                                 # Kopiraj podatke iz volumetričnega noda v exportable
                                if ct_volume_Node.GetAttribute("DICOM.PatientName"):
                                    exp.setTag("0010,0010", ct_volume_Node.GetAttribute("DICOM.PatientName"))  # PatientName
                                if ct_volume_Node.GetAttribute("DICOM.PatientID"):
                                    exp.setTag("0010,0020", ct_volume_Node.GetAttribute("DICOM.PatientID"))    # PatientID
                                if ct_volume_Node.GetAttribute("DICOM.StudyInstanceUID"):
                                    exp.setTag("0020,000D", ct_volume_Node.GetAttribute("DICOM.StudyInstanceUID"))  # StudyInstanceUID
                                exp.directory = export_dir
                                print(f"[{studyName}] ✅ Export successful")
                            plugin.export(exportables)

                                            
                        # 📏 Izračun napake
                        errors = [np.linalg.norm(cbct - ct) for cbct, ct in zip(cbct_points, ct_points)]
                        mean_error = np.mean(errors)

                        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"Specific marker errors {i+1}: ΔX={diff[0]:.2f} mm, ΔY={diff[1]:.2f} mm, ΔZ={diff[2]:.2f} mm")
                        
                        ct_pts = np.array(ct_points)
                        cbct_pts = np.array(cbct_points)

                        if len(ct_pts) == 3 and len(cbct_pts) == 3:
                            sim = triangle_similarity(ct_pts, cbct_pts)
                            print("\n📐 Triangle similarity analysis:")
                            print(f"Side length differences (mm): {sim['side_length_diff_mm']}")
                            print(f"Angle differences (deg): {sim['angle_diff_deg']}")
                            print(f"Mean length error: {sim['mean_length_error_mm']:.2f} mm")
                            print(f"Mean angle error: {sim['mean_angle_error_deg']:.2f}°")
                            print(f"Triangle Similarity Ratio (TSR): {sim['TSR']:.3f}")
                        else:
                            print("⚠️ Za primerjavo trikotnikov potrebne točno 3 točke.")
                        
                        if writefilecheck:

 
                            errorsfile = os.path.join(study_dir, f"{cbct_date}_errors.csv")

                            with open(errorsfile, mode='a', newline='') as file:
                                writer = csv.writer(file)

                                # Glava za študijo
                                writer.writerow(["Study", studyName])
                                writer.writerow(["Total Individual Errors (mm)"])
                                
                                for error in errors:
                                    writer.writerow(["", f"{error:.2f}"])

                                writer.writerow(["Average Error (mm)", f"{mean_error:.2f}"])

                                # Specifične napake markerjev
                                writer.writerow(["Specific Marker Errors"])
                                writer.writerow(["Marker", "ΔX (mm)", "ΔY (mm)", "ΔZ (mm)"])

                                for i, (cbct, ct) in enumerate(zip(cbct_points, ct_points)):
                                    diff = np.array(cbct) - np.array(ct)
                                    writer.writerow([f"Marker {i+1}", f"{diff[0]:.2f}", f"{diff[1]:.2f}", f"{diff[2]:.2f}"])

                                if len(cbct_points) == 3 and len(ct_points) == 3:
                                    writer.writerow([])
                                    writer.writerow(["Triangle Similarity"])
                                    writer.writerow(["Side Length Diff (mm)", *[f"{v:.2f}" for v in sim["side_length_diff_mm"]]])
                                    writer.writerow(["Angle Diff (deg)", *[f"{v:.2f}" for v in sim["angle_diff_deg"]]])
                                    writer.writerow(["Mean Length Error (mm)", f"{sim['mean_length_error_mm']:.2f}"])
                                    writer.writerow(["Mean Angle Error (deg)", f"{sim['mean_angle_error_deg']:.2f}"])
                                    writer.writerow(["Triangle Similarity Ratio (TSR)", f"{sim['TSR']:.3f}"])
                                    writer.writerow([])
                                 
                            
                            graphs_dir = os.path.join(results_base, study_folder)
                            os.makedirs(graphs_dir, exist_ok=True)

                            pngfile = os.path.join(graphs_dir, f"{cbct_date}_trikotnika.png")
                            save_triangle_visualization(np.array(ct_points), np.array(cbct_points), pngfile)
                                

                            #print("Prostate size file written at ", prostate_size_file)
                        
                        
                        
                
            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,
                "Table": ((table1end_time - table1_time)+(table2end_time - table2_time)) if tablefind else 0,
                "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)
            
            # 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 prostate size data
            with open(prostate_size_file, mode='w', newline='') as file:
                writer = csv.writer(file)
                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_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
        # Convert to minutes and seconds
        minutes = int(elapsed_time // 60)
        seconds = elapsed_time % 60

        print(f"Execution time: {minutes} minutes and {seconds:.6f} seconds")