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
from datetime import datetime
import csv
import time
import logging
import matplotlib
matplotlib.use('Agg')  # << to dodaš ZGORAJ, da omogoči PNG zapis brez GUI
import matplotlib.pyplot as plt
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(False)
        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.FineShiftCheckBox = qt.QCheckBox("Use extra fine shift correction")
        self.FineShiftCheckBox.setChecked(True)
        self.layout.addWidget(self.FineShiftCheckBox)
        
        self.save_as_dicomCheckBox = qt.QCheckBox("Save transformed CT and segmentations as DICOM files")
        self.save_as_dicomCheckBox.setChecked(False)
        self.layout.addWidget(self.save_as_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()
        use_fine_shift = self.FineShiftCheckBox.isChecked()
        save_as_dicom = self.save_as_dicomCheckBox.isChecked()

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


class MyTransformModuleLogic(ScriptedLoadableModuleLogic):
    """
    Core logic of the module.
    """
    
    
    def run(self, selectedMethod, applyRotation, applyTranslation, applyScaling, applymarkers, writefilecheck, tablefind, use_fine_shift, 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")
                    print(" ".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=200, y_max=380, z_min=25, 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(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 ---
            volume_node = find_volume_node_by_partial_name(volume_name)
            np_array = slicer.util.arrayFromVolume(volume_node)  # (Z, Y, X)
            ijkToRasMatrix = vtk.vtkMatrix4x4()
            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"candidate_y_values: {candidate_y_values}")
                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 - CBCT_offset

                #print(f"CT offset: {CT_offset}, CBCT offset: {CBCT_offset}")
                #print(f"CT spacing: {CT_spacing} mm/voxel, CBCT spacing: {volumeNode.GetSpacing()[1]} mm/voxel")
                #print(f"Aligning CBCT with CT. Offset in mm: {alignment_offset_mm}")

                # Uporabi transformacijo
                transform = vtk.vtkTransform()
                transform.Translate(0, alignment_offset_mm, 0)

                transformNode = slicer.vtkMRMLTransformNode()
                slicer.mrmlScene.AddNode(transformNode)
                transformNode.SetAndObserveTransformToParent(transform)
                volumeNode.SetAndObserveTransformNodeID(transformNode.GetID())
                slicer.vtkSlicerTransformLogic().hardenTransform(volumeNode)
                slicer.mrmlScene.RemoveNode(transformNode)
                
                # 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 outlier: točko z največjo 3D razdaljo do robustnega centroida (median).
            Združljivo: 'points' je lahko [XYZ] ali [(XYZ, intensity)].

            :param points: list[np.array([x,y,z])] ali list[(np.array([x,y,z]), meta)]
            :return: isti list, z odstranjenim outlierjem
            """
            import numpy as np

            if not points or len(points) <= 3:
                # nič ne odstrani, potrebujemo vsaj 4 kandidate (npr. 3 markerji + 1 koža)
                return points

            # Normaliziraj vhod v Nx3 matriko koordinat + obdrži indeksno preslikavo
            coords = []
            for p in points:
                if isinstance(p, (list, tuple)) and len(p) == 2 and hasattr(p[0], "__len__"):
                    xyz = np.asarray(p[0], dtype=float)  # (centroid, intensity)
                else:
                    xyz = np.asarray(p, dtype=float)     # samo centroid
                coords.append(xyz[:3])
            A = np.vstack(coords)  # (N,3)

            # Robustni centroid
            center = np.median(A, axis=0)

            # Evklidske razdalje (kvadrati zadoščajo za max)
            d2 = np.sum((A - center) ** 2, axis=1)
            idx_out = int(np.argmax(d2))

            removed_point = A[idx_out]
            removed = points.pop(idx_out)

            print(
                f"⚠️ Odstranjen outlier (najbolj oddaljen od median-centroida): "
                f"{removed_point.round(2)} | d={float(np.sqrt(d2[idx_out])):.2f} mm | "
                f"centroid≈{center.round(2)}"
            )
            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())
                seg_node.SetAttribute("DICOM.StructureSetLabel", f"Transformiran_{cbct_date}")
                seg_node.SetAttribute("DICOM.SeriesDescription", f"Transformiran {studyName.replace('^', ' ')} ({cbct_date})")
                seg_node.SetAttribute("DICOM.ReferencedSeriesInstanceUID", ct_volume_node.GetAttribute("DICOM.SeriesInstanceUID"))
                seg_node.SetAttribute("DICOM.ReferencedStudyInstanceUID", ct_volume_node.GetAttribute("DICOM.StudyInstanceUID"))




                # 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 dicom_timestamp_from_volume(volume_node):
            """Vrne npr. '20250112_093045' iz DICOM metapodatkov ali None, če jih ni."""
            uids_attr = volume_node.GetAttribute('DICOM.instanceUIDs')
            if not uids_attr:
                return None
            uid0 = uids_attr.split()[0]
            db = slicer.dicomDatabase
            fpath = db.fileForInstance(uid0)

            def tag(tagstr):
                v = db.fileValue(fpath, tagstr)
                return v if v not in (None, "", "Unknown") else None

            # Datumi: Acquisition (0008,0022) → Series (0008,0021) → Study (0008,0020)
            date = tag("0008,0022") or tag("0008,0021") or tag("0008,0020")
            # Časi:   Acquisition (0008,0032) → Series (0008,0031) → Study (0008,0030)
            time = tag("0008,0032") or tag("0008,0031") or tag("0008,0030")

            if not date:
                return None

            # Normaliziraj: YYYYMMDD in HHMMSS
            date = date.replace("-", "").replace(".", "").strip()
            date = (date + "00000000")[:8]          # zaščita, če je prekratko
            if time:
                time = time.split(".")[0]           # odreži frakcije
                time = (time + "000000")[:6]
                return f"{date}_{time}"
            return date

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

            Če je podanih 4 točk, avtomatsko odstrani eno (ujemajočo), ki najmanj vpliva na podobnost.

            :param p1: seznam (np.array) treh ali štirih točk v 3D (npr. ct_points)
            :param p2: seznam treh ali štirih točk v 3D (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):
                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])

            # Če imamo 4 točke v vsaki množici, iščemo najboljše 3 ujemajoče
            if len(p1) == 4 and len(p2) == 4:
                best_score = np.inf
                best_idx = None
                for i in range(4):
                    tri1 = np.delete(p1, i, axis=0)
                    tri2 = np.delete(p2, i, axis=0)
                    l1 = side_lengths(tri1)
                    l2 = side_lengths(tri2)
                    a1 = angles(tri1)
                    a2 = angles(tri2)
                    score = np.mean(np.abs(l1 - l2)) + np.mean(np.abs(a1 - a2)) / 10
                    if score < best_score:
                        best_score = score
                        best_idx = i
                print(f"📐 Samodejno odstranjena točka {best_idx} za trikotniško primerjavo.")
                p1 = np.delete(p1, best_idx, axis=0)
                p2 = np.delete(p2, best_idx, axis=0)

            if len(p1) != 3 or len(p2) != 3:
                raise ValueError("Obe množici točk morata vsebovati 3 (ali 4) točke.")

            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 _pick_dose_node():
            """Izbere RTDOSE; če obstaja 'Accumulated_*', jo vzame prednostno."""
            dose_nodes = [n for n in slicer.util.getNodesByClass("vtkMRMLScalarVolumeNode")
                        if "DOSE" in (n.GetName() or "").upper()]
            if not dose_nodes:
                return None
            acc = [n for n in dose_nodes if "ACCUMULATED" in n.GetName().upper()]
            return acc[0] if acc else dose_nodes[0]

        def _pick_segmentation_node():
            """Vzemi RTSTRUCT segmentacijo (ali prvo smiselno segmentacijo)."""

            cand = slicer.util.getNodesByClass("vtkMRMLSegmentationNode")
            if not cand:
                return None
            # prednostno RTSTRUCT*
            for n in cand:
                if "RTSTRUCT" in n.GetName().upper():
                    return n
            return cand[0]

        def _segment_ids(seg_node, name_filters=None):
            """Vrne [(segmentId, name), ...] (filtrirano po imenih, če je podano)."""
            seg = seg_node.GetSegmentation()
            seg.CreateRepresentation("BinaryLabelmap")
            all_ids = vtk.vtkStringArray()
            seg.GetSegmentIDs(all_ids)
            out = []
            all_names = []
            for i in range(all_ids.GetNumberOfValues()):
                sid = all_ids.GetValue(i)
                name = seg.GetSegment(sid).GetName() or sid
                all_names.append(name)
                if (not name_filters) or any(f.lower() in name.lower() for f in name_filters):
                    out.append((sid, name))
            #print(f"[DVH/DEBUG] _segment_ids: total={all_ids.GetNumberOfValues()} names={all_names}")
            #if name_filters:
                #print(f"[DVH/DEBUG] _segment_ids: filters={name_filters} -> kept={[n for _,n in out]}")
            return out

        def compute_dvh_numpy(dose_node, seg_node, name_filters=None, out_csv_path=None):
            """
            Izračuna DVH z notraj-voxel histogramiranjem (numpy).
            Segmentacijo sproti preslika v geometrijo doze (ExportSegmentsToLabelmapNode),
            zato ni več težav z 'no overlap' ali geometrijo.
            Vrne: list [ [Name,Dmin,Dmax,Dmean,D98,D95,D50,D2], ... ]
            """
            logic = slicer.modules.segmentations.logic()
            results = []
            #print("[DVH/DEBUG] compute_dvh_numpy: dose =", dose_node.GetName(), "shape=",
              #  slicer.util.arrayFromVolume(dose_node).shape)
            #print("[DVH/DEBUG] compute_dvh_numpy: seg  =", seg_node.GetName())
            
            try:
                src_rep = seg_node.GetSegmentation().GetSourceRepresentationName()
            except Exception:
                src_rep = "?"
            #print("[DVH/DEBUG] compute_dvh_numpy: seg =", seg_node.GetName(), "sourceRep=", src_rep)



            seg_list = _segment_ids(seg_node, name_filters)
            if not seg_list:
                print("[DVH] Ni segmentov za obdelavo.")
                return results

            dose_arr = slicer.util.arrayFromVolume(dose_node).astype(np.float64)

            for sid, sname in seg_list:
                #print(f"[DVH/DEBUG] Processing segment: '{sname}' (id={sid})")

                # en segment → labelmap v geometriji doze
                lm = _labelmap_of_segment_on_dose(seg_node, dose_node, sid)
                if lm is None:
                    print(f"[DVH] Prazna maska za '{sname}'.")
                    continue
                arrLM = slicer.util.arrayFromVolume(lm)
                #print("[DVH/DEBUG]  labelmap shape:", None if arrLM is None else arrLM.shape,
                    #" nonzero:", None if arrLM is None else int((arrLM>0).sum()))
                
                seg_api = seg_node.GetSegmentation()
                #print("[DVH/DEBUG] Export using sourceRep:", seg_api.GetSourceRepresentationName(), " -> dose:", dose_node.GetName())

                mask = arrLM > 0 if arrLM is not None else None
                slicer.mrmlScene.RemoveNode(lm)
                if mask is None or not mask.any():
                    print(f"[DVH] Prazna maska za '{sname}'.")
                    continue
                vox = dose_arr[mask]
                #print(f"[DVH/DEBUG]  voxels in mask: {vox.size}, Dmean={vox.mean():.2f} Gy")


                # osnovne metrike (Gy)
                Dmin = float(np.min(vox))
                Dmax = float(np.max(vox))
                Dmean = float(np.mean(vox))
                D98 = float(np.percentile(vox, 2))
                D95 = float(np.percentile(vox, 5))
                D50 = float(np.percentile(vox, 50))
                D2  = float(np.percentile(vox, 98))


                results.append([sname, Dmin, Dmax, Dmean, D98, D95, D50, D2])
            if out_csv_path and results:
                os.makedirs(os.path.dirname(out_csv_path), exist_ok=True)
                with open(out_csv_path, "w", newline="") as f:
                    w = csv.writer(f)
                    w.writerow(["Segment","Dmin(Gy)","Dmax(Gy)","Dmean(Gy)","D98(Gy)","D95(Gy)","D50(Gy)","D2(Gy)"])
                    w.writerows(results)

            return results

        def _for_uid(volume_node):
            """Vrne DICOM FrameOfReferenceUID (0020,0052) ali None."""
            # 1) neposredno z MRML atributa (če je prisoten)
            direct = volume_node.GetAttribute('DICOM.FrameOfReferenceUID')
            if direct:
                return direct
            # 2) iz DICOM DB prek prve instance
            uids_attr = volume_node.GetAttribute('DICOM.instanceUIDs')
            if not uids_attr:
                return None
            db = slicer.dicomDatabase
            f = db.fileForInstance(uids_attr.split()[0])
            return db.fileValue(f, "0020,0052") or None

        def _pick_matching_dose_node(ct_node):
            """Vrne RTDOSE z istim FORUID kot CT (če obstaja), sicer None/first fallback."""
            ct_for = _for_uid(ct_node) if ct_node else None
            dose_nodes = [n for n in slicer.util.getNodesByClass("vtkMRMLScalarVolumeNode")
                        if "DOSE" in (n.GetName() or "").upper()]
            if not dose_nodes:
                return None
            if ct_for:
                for dn in dose_nodes:
                    if _for_uid(dn) == ct_for:
                        return dn
            # fallback: prva doza (bolje kot nič, a lahko ne prekriva)
            return dose_nodes[0]

        def _ensure_overlap_by_resample(dose_node, ct_node):
            """Če doza in CT nimata istega FORUID, dozo resampla na CT mrežo in vrne novo dozo."""
            if not dose_node or not ct_node:
                return dose_node
            if _for_uid(dose_node) == _for_uid(ct_node):
                return dose_node  # že matcha
            import slicer
            newDose = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLScalarVolumeNode",
                                                        dose_node.GetName() + "_onCT")
            params = {
                "inputVolume": dose_node.GetID(),
                "referenceVolume": ct_node.GetID(),
                "outputVolume": newDose.GetID(),
                "interpolationType": "linear",
            }
            slicer.cli.runSync(slicer.modules.resamplescalarvectordwivolume, None, params)
            print("[DVH] Fallback: resampled dose to CT grid ->", newDose.GetName())
            return newDose

        def _labelmap_of_segment_on_dose(seg_node, dose_node, segment_id):
            """Vrne TMP labelmap node za en segment v geometriji doze (ali None ob neuspehu)."""
            logic = slicer.modules.segmentations.logic()

            # 0) vedno referenciraj geometrijo na DOZO
            seg_node.SetReferenceImageGeometryParameterFromVolumeNode(dose_node)

            # 1) poskusi direktni export
            sid = vtk.vtkStringArray(); sid.InsertNextValue(segment_id)
            lm = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLLabelMapVolumeNode", f"tmpLM_{segment_id}")
            ok = logic.ExportSegmentsToLabelmapNode(seg_node, sid, lm, dose_node)
            if ok:
                arr = slicer.util.arrayFromVolume(lm)
                if arr is not None and (arr > 0).any():
                    return lm
            slicer.mrmlScene.RemoveNode(lm)

            # 2) FALLBACK: prisili pot PlanarContour -> Closed surface -> BinaryLabelmap na mreži doze
            seg = seg_node.GetSegmentation()
            seg.CreateRepresentation("Closed surface")
            seg.SetMasterRepresentationName("Closed surface")

            # kopiraj samo ta segment v novo, "čisto" segmentacijo na mreži doze
            seg_tmp = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLSegmentationNode", "Seg_onDose_TMP")
            seg_tmp.SetReferenceImageGeometryParameterFromVolumeNode(dose_node)
            seg_tmp.GetSegmentation().CopySegmentFromSegmentation(seg, segment_id, segment_id)
            seg_tmp.GetSegmentation().CreateRepresentation("BinaryLabelmap")

            lm2 = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLLabelMapVolumeNode", f"tmpLM2_{segment_id}")
            ok2 = logic.ExportSegmentsToLabelmapNode(seg_tmp, sid, lm2, dose_node)

            slicer.mrmlScene.RemoveNode(seg_tmp)
            if not ok2:
                slicer.mrmlScene.RemoveNode(lm2)
                return None

            arr2 = slicer.util.arrayFromVolume(lm2)
            if arr2 is None or not (arr2 > 0).any():
                slicer.mrmlScene.RemoveNode(lm2)
                return None
            return lm2

        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)
            
        def resample_scalar_to_reference(input_node, reference_node, interpolator="linear"):
            """Resampla input_node na geometrijo reference_node in PREPIŠE image data + IJK↔RAS + origin/spacing."""
            out = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLScalarVolumeNode", input_node.GetName() + "_tmp")
            params = {
                "inputVolume": input_node.GetID(),
                "referenceVolume": reference_node.GetID(),
                "outputVolume": out.GetID(),
                "interpolationType": "linear" if interpolator=="linear" else "nn",
            }
            slicer.cli.runSync(slicer.modules.resamplescalarvectordwivolume, None, params)

            # --- PREPIŠI ImageData ---
            input_node.SetAndObserveImageData(out.GetImageData())

            # --- PREPIŠI GEOMETRIJO (IJK↔RAS) ---
            ijkToRas = vtk.vtkMatrix4x4(); out.GetIJKToRASMatrix(ijkToRas)
            rasToIjk = vtk.vtkMatrix4x4(); out.GetRASToIJKMatrix(rasToIjk)
            input_node.SetIJKToRASMatrix(ijkToRas)
            input_node.SetRASToIJKMatrix(rasToIjk)

            # --- PREPIŠI origin/spacing (za vsak slučaj, čeprav IJK↔RAS običajno zadošča) ---
            input_node.SetOrigin(out.GetOrigin())
            input_node.SetSpacing(out.GetSpacing())

            # Po želji: prenesi še DisplayNode lastnosti ipd.
            slicer.mrmlScene.RemoveNode(out)
            
        def export_seg_or_rtstruct(ct_volume_Node, seg_node, export_dir):
            """
            Poskusi izvoziti RTSTRUCT (Planar contour). Če ne uspe, izvozi DICOM SEG.
            Zahteve:
            - seg_node je že PORAVNAN na CT (apliciran transform)
            - reference image geometry in UID-ji so nastavljeni
            """
            import slicer
            from DICOMPlugins import DicomRtImportExportPlugin
            from DICOMExportSegmentations import DICOMSegmentationExporter

            # 0) Referenčna geometrija in ujemanje v SubjectHierarchy
            seg_node.SetReferenceImageGeometryParameterFromVolumeNode(ct_volume_Node)

            sh = slicer.vtkMRMLSubjectHierarchyNode.GetSubjectHierarchyNode(slicer.mrmlScene)
            ct_item  = sh.GetItemByDataNode(ct_volume_Node)
            seg_item = sh.GetItemByDataNode(seg_node)

            # poskrbi, da sta CT in SEG pod ISTIM Study (DICOM exporter to uporablja)
            study_ct = sh.GetItemParent(ct_item)
            if study_ct and study_ct != sh.GetItemParent(seg_item):
                sh.SetItemParent(seg_item, study_ct)

            # 1) Vnesi ključne DICOM atribute iz CT na SEG (Reference Series in FORUID)
            ct_series_uid = ct_volume_Node.GetAttribute("DICOM.SeriesInstanceUID") or ""
            ct_study_uid  = ct_volume_Node.GetAttribute("DICOM.StudyInstanceUID") or ""
            ct_for_uid    = ct_volume_Node.GetAttribute("DICOM.FrameOfReferenceUID") or ""

            if ct_series_uid:
                sh.SetItemAttribute(seg_item, "DICOM.ReferencedInstanceUIDs", ct_series_uid)
            if ct_for_uid:
                seg_node.SetAttribute("DICOM.FrameOfReferenceUID", ct_for_uid)

            # 2) Poskrbi za reprezentance in nastavi "Source" (nova API namesto "Master")
            seg = seg_node.GetSegmentation()

            # Closed surface naj obstaja
            if not seg.ContainsRepresentation("Closed surface"):
                seg.CreateRepresentation("Closed surface")

            # Poskusi ustvariti Planar contour (zahteva veljavno referenčno geometrijo)
            try:
                if not seg.ContainsRepresentation("Planar contour"):
                    seg.CreateRepresentation("Planar contour")
            except Exception as e:
                print("[RTSTRUCT] CreateRepresentation('Planar contour') failed:", e)

            # Najprej probamo s Planar contour kot "Source" (nov API)
            if seg.ContainsRepresentation("Planar contour"):
                seg.SetSourceRepresentationName("Planar contour")
                source_repr = "Planar contour"
            else:
                seg.SetSourceRepresentationName("Closed surface")
                source_repr = "Closed surface"

            print(f"[RTSTRUCT] Source: {seg.GetSourceRepresentationName()}")
            print(f"[RTSTRUCT] Has Planar: {seg.ContainsRepresentation('Planar contour')}")
            print(f"[RTSTRUCT] Has Closed: {seg.ContainsRepresentation('Closed surface')}")
            print("Ref Series UID:", sh.GetItemAttribute(seg_item, "DICOM.ReferencedInstanceUIDs"))
            print("FOR UID:", seg_node.GetAttribute("DICOM.FrameOfReferenceUID"))

            # 3) Pripravi exportables in filtriraj RTSTRUCT kandidat(e)
            plugin = DicomRtImportExportPlugin.DicomRtImportExportPluginClass()

            def fill_tags(exp):
                pn = ct_volume_Node.GetAttribute("DICOM.PatientName")
                pid = ct_volume_Node.GetAttribute("DICOM.PatientID")
                study_uid = ct_volume_Node.GetAttribute("DICOM.StudyInstanceUID")
                if pn:   exp.setTag("0010,0010", pn)
                if pid:  exp.setTag("0010,0020", pid)
                if study_uid: exp.setTag("0020,000D", study_uid)
                exp.directory = export_dir

            exp_all = []
            exp_all += plugin.examineForExport(sh.GetItemByDataNode(ct_volume_Node))
            exp_all += plugin.examineForExport(sh.GetItemByDataNode(seg_node))

            def pick_rtstruct(exportables):
                cand = []
                for e in exportables:
                    name   = getattr(e, "name", "")
                    etype  = getattr(e, "exportType", "")
                    sop    = getattr(e, "SOPClassUID", "")
                    if etype == "RTSTRUCT" or "RTSTRUCT" in (name or "") or sop == "1.2.840.10008.5.1.4.1.1.481.3":
                        cand.append(e)
                return cand

            cand_rt = pick_rtstruct(exp_all)

            # Če RTSTRUCT ni našel in trenutno source ni Closed surface, poskusi še enkrat
            if not cand_rt and source_repr != "Closed surface":
                seg.SetSourceRepresentationName("Closed surface")
                print("[RTSTRUCT] Retry with source=Closed surface")
                exp_all = []
                exp_all += plugin.examineForExport(sh.GetItemByDataNode(ct_volume_Node))
                exp_all += plugin.examineForExport(sh.GetItemByDataNode(seg_node))
                cand_rt = pick_rtstruct(exp_all)

            # 4) Izvoz
            if cand_rt:
                for e in cand_rt:
                    fill_tags(e)
                ok = plugin.export(cand_rt)
                print(f"✅ RTSTRUCT export status: {ok} → {export_dir}")
                return "RTSTRUCT"
            else:
                # Fallback: DICOM SEG (OpenTPS 2.x ima odlično podporo in DVH deluje)
                options = DICOMSegmentationExporter.ExportOptions()
                options.outputDirectory = export_dir
                options.segmentationNode = seg_node
                options.referencedVolumeNode = ct_volume_Node
                DICOMSegmentationExporter.export(options)
                print(f"✅ DICOM SEG exported → {export_dir} (DVH/CCC v OpenTPS 2.x deluje tudi iz SEG)")
                return "SEG"

            
            
            
            
        # 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}")
                
                     
                
                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}")
                        tablefound = False
                        ct_table_found = False
                        cbct_table_found = False
                        CT_offset = None
                        CT_spacing = None                        
                        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)
                        table_shift_matrix = np.eye(4)
                        fixing = time.time()
                        if tablefind:
                            # --- CT table reference ---
                            yesCbct = False
                            table1_time = time.time()
                            makemarkerscheck = True
                            resCT = find_table_top_z(ct_volume_name, writefilecheck, makemarkerscheck, yesCbct)
                            if resCT is not None:
                                ct_edge_mm, _ = resCT
                                ct_table_found = True
                                # zapišemo referenco CT (ne premikamo nič)
                                CT_offset, CT_spacing = align_cbct_to_ct(ct_volume_Node, "CT", ct_edge_mm)
                            else:
                                print("⚠️ CT table top not found – skip table alignment.")
                                ct_table_found = False
                            table1end_time = time.time()

                            # --- CBCT table measure & align ---
                            yesCbct = True
                            table2_time = time.time()
                            makemarkerscheck = False
                            if ct_table_found:
                                resCBCT = find_table_top_z(cbct_volume_name, writefilecheck, makemarkerscheck, yesCbct)
                                #skupnires = resCBCT - resCT if resCBCT is not None else None
                                if resCBCT is not None:
                                    cbct_edge_mm, _ = resCBCT
                                    # Poravnava CBCT -> CT (harden transform na CBCT)
                                    align_cbct_to_ct(cbct_volume_node, "CBCT", cbct_edge_mm, CT_offset, CT_spacing)
                                    tablefound = True
                                    # ⛔️ NE dodajaj table_shift_matrix v cumulative_matrices (hardeno je že na volumnu)
                                else:
                                    print("⚠️ CBCT table top not found – skip table alignment.")
                                    tablefound = False
                            table2end_time = time.time()
                        else:
                            tablefound = False       
                        
                        resample_scalar_to_reference(cbct_volume_node, ct_volume_Node, interpolator="nn")  # Resampla CBCT na CT mrežo                    
                        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()
                        if not np.allclose(cbct_spacing, ct_spacing, atol=1e-6):
                            cbct_points = rescale_points_to_match_spacing(cbct_points, cbct_spacing, ct_spacing)
                        # sicer ne reskaliramo, ker je grid že ujemajoč

                        
                        #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 (Horn)":
                                _, 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
                            dy = mean_delta_y
                            if tablefind and tablefound:
                                dy = 0.0  # meja: ali np.clip(dy, -3.0, 3.0)
                            fine_shift = np.array([0.0, dy, 0.0])

                            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, tablefound, 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

                        

                        
                        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)

                        pre_fine_matrices = {}
                        pre_fine_matrices[(studyName, cbct_volume_name)] = cumulative_matrices[(studyName, cbct_volume_name)].copy()
                        fine_shift = np.zeros(3)
                        if(use_fine_shift):
                            # 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)
                                print(f"Fine correction shifts: ΔX={fine_shift[0]:.2f} mm, ΔY={fine_shift[1]:.2f} mm, ΔZ={fine_shift[2]:.2f} mm")

                        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, tablefound, 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] =   # 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)]
                        
                        #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}")
                        
                        M_cbct2ct_RAS = cumulative_matrices[(studyName, cbct_volume_name)].copy()

                        # RAS<->LPS pretvornik
                        R_lps = np.diag([-1.0, -1.0, 1.0, 1.0])

                        # CBCT->CT v LPS
                        M_cbct2ct_LPS = R_lps @ M_cbct2ct_RAS @ R_lps

                        # Kar želiš aplicirati na CT v Ariji/Eclipse (CT->CBCT), rigid v mm, LPS:
                        M_ct2cbct_LPS = np.linalg.inv(M_cbct2ct_LPS)
                        
                        if (resCBCT is not None) and (resCT is not None):
                            deltaY_RAS = float(resCBCT[0]) - float(resCT[0])   # samo mm komponenta (RAS)
                        else:
                            deltaY_RAS = 0.0
                        # RAS -> LPS: X in Y zamenjata predznak, Z ostane
                        deltaY_LPS = -deltaY_RAS
                        M_ct2cbct_LPS = M_ct2cbct_LPS.copy()
                        M_ct2cbct_LPS[1,3] += deltaY_LPS

                        
                        #shrani transformacijsko matriko v datoteko
                        save_transform_matrix(M_ct2cbct_LPS, studyName, cbct_volume_name)
                        
                        Tx, Ty, Tz = float(M_ct2cbct_LPS[0,3]), float(M_ct2cbct_LPS[1,3]), float(M_ct2cbct_LPS[2,3])
                        R = M_ct2cbct_LPS[:3,:3]
                        ry = np.degrees(np.arcsin(np.clip(-R[2,0], -1.0, 1.0)))
                        rx = np.degrees(np.arctan2(R[2,1], R[2,2]))
                        rz = np.degrees(np.arctan2(R[1,0], R[0,0]))
                        print(f"Eclipse 6DoF (LPS): Tx={Tx:.2f} mm, Ty={Ty:.2f} mm, Tz={Tz:.2f} mm | Rx={rx:.3f}°, Ry={ry:.3f}°, Rz={rz:.3f}°")

                        # M_ct2cbct_LPS že imaš
                        R = np.diag([-1.0, -1.0, 1.0, 1.0])
                        M_ct2cbct_RAS = R @ M_ct2cbct_LPS @ R  # pretvorba LPS->RAS

                        # 1) ustvari transform node in nastavi 4x4
                        tNode = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLLinearTransformNode", "CTmove_to_CBCT")
                        m = vtk.vtkMatrix4x4()
                        for r in range(4):
                            for c in range(4):
                                m.SetElement(r, c, float(M_ct2cbct_RAS[r, c]))
                        tNode.SetMatrixTransformToParent(m)

                        # 2) apliciraj na CT in (po želji) na strukture/segmente
                        #ctNode = find_volume_node_by_partial_name(ct_volume_name)  # ime CT-ja, ki ga že uporabljaš
                        #ctNode.SetAndObserveTransformNodeID(tNode.GetID())
                        #slicer.vtkSlicerTransformLogic().hardenTransform(ctNode)

                        # če imaš segmentation/RTSTRUCT kot segmentation node:
                        # segNode = slicer.util.getNode("Segmentation")  # prilagodi ime
                        # segNode.SetAndObserveTransformNodeID(tNode.GetID()); slicer.vtkSlicerTransformLogic().hardenTransform(segNode)

                        #print("✅ CT je bil transformiran (harden) v RAS-mm po M_ct2cbct.")
                        
                        
                        #Pridobi podatke o isocentru
                        # iso_coords = extract_isocenter_from_loaded_rtplan()
                        # if iso_coords is not None:
                        #     matrix = cumulative_matrices.get((studyName, cbct_volume_name))
                        #     if matrix is not None:
                        #         homogeneous_iso = np.append(iso_coords, 1.0)  # Doda homogeni člen
                        #         transformed_iso = np.dot(matrix, homogeneous_iso)
                        #         print(f"🧭 Transformirani izocenter: {transformed_iso[:3]} mm")
                        #         print(f"Razlika med izocentroma: ΔX={transformed_iso[0]-iso_coords[0]:.2f} mm, ΔY={transformed_iso[1]-iso_coords[1]:.2f} mm, ΔZ={transformed_iso[2]-iso_coords[2]:.2f} mm")
                        
                        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)

                        cbct_date = dicom_timestamp_from_volume(cbct_volume_node) or datetime.now().strftime("%Y%m%d_%H%M%S")
                        export_dir = os.path.join(study_dir, f"{cbct_date}_DICOM")
                        os.makedirs(export_dir, exist_ok=True)
                        
                        
                        #save_as_dicom = False 
                        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_acq = cbct_volume_node.GetAttribute("DICOM.AcquisitionDate") or "unknownDate"
                            #print(f"CBCT date (SeriesDate): {cbct_date}, AcquisitionDate: {cbct_date_acq}")   
                   
                            
                            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 in [{export_dir}]")
                            plugin.export(exportables)
                            

                            # def export_seg_as_dicom_seg(seg_node, ref_ct_node, out_dir, series_desc="AutoExport SEG"):
                            #     # 0) Referenčna geometrija in hierarhija
                            #     seg_node.SetReferenceImageGeometryParameterFromVolumeNode(ref_ct_node)
                            #     sh = slicer.vtkMRMLSubjectHierarchyNode.GetSubjectHierarchyNode(slicer.mrmlScene)
                            #     ct_item  = sh.GetItemByDataNode(ref_ct_node)
                            #     seg_item = sh.GetItemByDataNode(seg_node)
                            #     study_ct = sh.GetItemParent(ct_item)
                            #     if study_ct and study_ct != sh.GetItemParent(seg_item):
                            #         sh.SetItemParent(seg_item, study_ct)

                            #     # 1) Poskrbi za binarni labelmap kot source (DICOM SEG to pričakuje)
                            #     seg = seg_node.GetSegmentation()
                            #     if not seg.ContainsRepresentation("Binary labelmap"):
                            #         seg.CreateRepresentation("Binary labelmap")
                            #     seg.SetSourceRepresentationName("Binary labelmap")

                            #     # 2) Eksplicitno nastavi DICOM reference (pomoč, če manjka v SH atributih)
                            #     ct_series_uid = ref_ct_node.GetAttribute("DICOM.SeriesInstanceUID") or ""
                            #     ct_for_uid    = ref_ct_node.GetAttribute("DICOM.FrameOfReferenceUID") or ""
                            #     if ct_series_uid:
                            #         sh.SetItemAttribute(seg_item, "DICOM.ReferencedInstanceUIDs", ct_series_uid)
                            #     if ct_for_uid:
                            #         seg_node.SetAttribute("DICOM.FrameOfReferenceUID", ct_for_uid)

                            #     # 3) IZVOZ (brez dodatnih modulov)
                            #     logic = slicer.modules.segmentations.logic()
                            #     ok = logic.ExportAllSegmentsToDICOMSegmentation(seg_node, ref_ct_node, out_dir, series_desc)
                            #     print(f"✅ DICOM SEG export status: {ok} → {out_dir}")
                            #     return ok

                            #export_seg_as_dicom_seg(seg_node, ct_volume_Node, export_dir)




                            


                                            
                        # 📏 Izračun napake
                        errors = [np.linalg.norm(cbct - ct) for cbct, ct in zip(cbct_points, ct_points)]
                        mean_error = np.mean(errors)
                        ct_pts = np.array(ct_points)
                        cbct_pts = np.array(cbct_points)
                        ct_centroid   = ct_pts.mean(axis=0)
                        cbct_centroid = cbct_pts.mean(axis=0)
                        centroid_delta = cbct_centroid - ct_centroid
                        centroid_norm  = float(np.linalg.norm(centroid_delta))

                        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")
                        


                        if len(ct_pts) == 3 and len(cbct_pts) == 3 or len(ct_pts) == 4 and len(cbct_pts) == 4:
                            sim = triangle_similarity(ct_pts, cbct_pts)
                            print("\n Triangle similarity analysis:")
                            print(f"Centroid Δ: ΔX={centroid_delta[0]:.2f} mm, ΔY={centroid_delta[1]:.2f} mm, ΔZ={centroid_delta[2]:.2f} mm | |Δ|={centroid_norm:.2f} mm")
                            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(["Centroid (CT) X/Y/Z (mm)", f"{ct_centroid[0]:.2f}", f"{ct_centroid[1]:.2f}", f"{ct_centroid[2]:.2f}"])
                                    writer.writerow(["Centroid (CBCT) X/Y/Z (mm)", f"{cbct_centroid[0]:.2f}", f"{cbct_centroid[1]:.2f}", f"{cbct_centroid[2]:.2f}"])
                                    writer.writerow(["Centroid Δ X/Y/Z (mm) | |Δ|", f"{centroid_delta[0]:.2f}", f"{centroid_delta[1]:.2f}", f"{centroid_delta[2]:.2f}", f"{centroid_norm:.2f}"])
                                    writer.writerow([])

                                    writer.writerow([])
                                 
                                #writer.writerow(["Isocenter at: ", transformed_iso[:3]])
                            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)
                                
                           
            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 pacient: {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")