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 itertools
from DICOMLib import DICOMUtils
from collections import deque, Counter
import vtk
from slicer.ScriptedLoadableModule import *
import qt
import matplotlib.pyplot as plt
import csv
import time

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

cumulative_matrices = {}

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

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

        # Dropdown menu za izbiro metode
        self.rotationMethodComboBox = qt.QComboBox()
        self.rotationMethodComboBox.addItems(["Kabsch", "Horn", "Iterative Closest Point (Horn)"])
        self.layout.addWidget(self.rotationMethodComboBox)

        # Checkboxi za transformacije
        self.scalingCheckBox = qt.QCheckBox("Scaling")
        self.scalingCheckBox.setChecked(True)
        self.layout.addWidget(self.scalingCheckBox)
        
        self.rotationCheckBox = qt.QCheckBox("Rotation")
        self.rotationCheckBox.setChecked(True)
        self.layout.addWidget(self.rotationCheckBox)

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

        self.markersCheckBox = qt.QCheckBox("Place control points for detected markers")
        self.markersCheckBox.setChecked(True)
        self.layout.addWidget(self.markersCheckBox)
        
        self.writefileCheckBox = qt.QCheckBox("Write distances to csv file")
        self.writefileCheckBox.setChecked(True)
        self.layout.addWidget(self.writefileCheckBox)

        # 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):
        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()

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


class MyTransformModuleLogic(ScriptedLoadableModuleLogic):
    """
    Core logic of the module.
    """
    
    
    def run(self, selectedMethod, applyRotation, applyTranslation, applyScaling, applymarkers, writefilecheck):
        start_time = time.time()
        print("Calculating...")
        
        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=True, fallback_if_worse=True, normalize_lengths=True, normalize_angles=False, min_distance=5):

            def side_lengths(points):
                lengths = [
                    np.linalg.norm(points[0] - points[1]),
                    np.linalg.norm(points[1] - points[2]),
                    np.linalg.norm(points[2] - points[0])
                ]
                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):
                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))

                return w_len * score_len + w_ang * score_ang

            cbct_points = list(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.5
                weight_angle = 0.5

            best_perm = None
            best_score = float('inf')
            #print(f"CT points: {ct_points}")
            for perm in itertools.permutations(cbct_points):
                perm = np.array(perm)
                score = permutation_score(perm, ct_lengths, ct_angles, weight_length, weight_angle)

                #print(f"Perm: {perm.tolist()} -> Score: {score:.4f}")

                if score < best_score:
                    best_score = score
                    best_perm = perm

            #print("CT centroid:", np.mean(ct_points, axis=0))
            #print("CBCT centroid (best perm):", np.mean(best_perm, axis=0))
            
            if fallback_if_worse:
                original_score = permutation_score(np.array(cbct_points), ct_lengths, ct_angles, weight_length, weight_angle)
                #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 matrix_to_array(vtk_transform):
            """
            Converts a vtkTransform to a 4x4 numpy array.
            """
            vtk_matrix = vtk.vtkMatrix4x4()
            vtk_transform.GetMatrix(vtk_matrix)
            np_matrix = np.zeros((4, 4))
            for i in range(4):
                for j in range(4):
                    np_matrix[i, j] = vtk_matrix.GetElement(i, j)
            return np_matrix
        
        def save_transform_matrix(matrix, study_name, cbct_volume_name, ct_table_found):
            """
            Appends the given 4x4 matrix to a text file under the given study folder.
            """
            base_folder = 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
            if ct_table_found:
                safe_cbct_name += "_MizaPoravnava"
            else:
                safe_cbct_name += "_NIMizaPoravnana"
            filename = os.path.join(study_folder, f"{safe_cbct_name}.txt")
            with open(filename, "a") as f:
                #f.write("Transformacija:\n")
                for row in matrix:
                    f.write(" ".join(f"{elem:.6f}" for elem in row) + "\n")
                f.write("\n")  # Dodaj prazen vrstico med transformacijami
            print(f"Transform matrix saved to {filename}")
                
        def detect_points_region_growing(volume_name, yesCbct, create_marker, intensity_threshold=3000, x_min=90, x_max=380, y_min=190, y_max=380, z_min=50, z_max=140, max_distance=9, centroid_merge_threshold=5):
            volume_node = find_volume_node_by_partial_name(volume_name)
            if not volume_node:
                raise RuntimeError(f"Volume {volume_name} not found.")
            
            image_data = volume_node.GetImageData()
            matrix = vtk.vtkMatrix4x4()
            volume_node.GetIJKToRASMatrix(matrix)

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

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

            visited = set()

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

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

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

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

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

            unique_centroids = []
            for centroid, intensity in centroids:
                if not any(np.linalg.norm(centroid - existing_centroid) < centroid_merge_threshold for existing_centroid, _ in unique_centroids):
                    unique_centroids.append((centroid, intensity))
            
            if create_marker:        
                markups_node = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLMarkupsFiducialNode", f"Markers_{volume_name}")
                for centroid, intensity in unique_centroids:
                    markups_node.AddControlPoint(*centroid)
                    markups_node.SetDisplayVisibility(False)
                    #print(f"Detected Centroid (RAS): {centroid}, Max Intensity: {intensity}")

            return unique_centroids

        def find_table_top_z(ct_volume_name, writefilecheck, makemarkerscheck, yesCbct):
            """
            Najde višino zgornjega roba mize v CT/CBCT volumnu in po želji doda marker v sceno.

            Args:
                ct_volume_name (str): Ime volumna.
                writefilecheck (bool): Ali naj se rezultat shrani v .csv.
                makemarkerscheck (bool): Ali naj se doda marker v 3D Slicer.
                yesCbct (bool): True, če je CBCT; False, če je CT.

            Returns:
                (float, int): Z komponenta v RAS prostoru, in Y indeks v slicerjevem volumnu.
            """
            # --- Pridobi volume node ---
            ct_volume_node = find_volume_node_by_partial_name(ct_volume_name)
            np_array = slicer.util.arrayFromVolume(ct_volume_node)  # (Z, Y, X)
            ijkToRasMatrix = vtk.vtkMatrix4x4()
            ct_volume_node.GetIJKToRASMatrix(ijkToRasMatrix)

            # --- Določimo lokacijo stolpca ---
            z_index = np_array.shape[0] // 2  # srednji slice
            y_size = np_array.shape[1]
            # x_index = int(np_array.shape[2] * 0.15)
            # x_index = max(0, min(x_index, np_array.shape[2] - 1))
            

            # --- Izračun spodnje tretjine (spodnji del slike) ---
            y_start = int(y_size * 2 / 3)
            slice_data = np_array[z_index, :, :]  # (Y, X)
            y_end = y_size  # Do dna slike
            #column_values = slice_data[y_start:y_end, x_index]  # (Y)

            # --- Parametri za rob ---
            threshold_high = -300 if yesCbct else -100
            threshold_low = -700 if yesCbct else -350
            min_jump = 100 if yesCbct else 100
            window_size = 4  # število voxelov nad/pod
            #previous_value = column_values[-1]
            table_top_y = None

            # --- Več stolpcev okoli x_index ---
            x_center = np_array.shape[2] // 2
            x_offset = 30  # 30 levo od sredine
            x_index_base = max(0, x_center - x_offset)
            candidate_y_values = []

            search_range = range(-5, 6)  # od -5 do +5 stolpcev
            
            for dx in search_range:
                x_index = x_index_base + dx
                if x_index < 0 or x_index >= np_array.shape[2]:
                    continue

                column_values = slice_data[y_start:y_end, x_index]

                for i in range(window_size, len(column_values) - window_size):
                    curr = column_values[i]
                    above_avg = np.mean(column_values[i - window_size:i])
                    below_avg = np.mean(column_values[i + 1:i + 1 + window_size])

                    if (threshold_low < curr < threshold_high
                        and (above_avg - below_avg) > min_jump
                        and below_avg < -400
                        and above_avg > -300):
                        y_found = y_start + i
                        candidate_y_values.append(y_found)
                        break  # samo prvi zadetek v stolpcu
            if candidate_y_values:
                most_common_y, _ = Counter(candidate_y_values).most_common(1)[0]
                table_top_y = most_common_y
                print(f"✅ Rob mize (najpogostejši Y): {table_top_y}, pojavitev: {candidate_y_values.count(table_top_y)}×")

            
            """ # --- 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 točko, ki je najnižja po dani osi (default Y-os v RAS prostoru).
            """
            if len(points) <= 3:
                return points  # Ni potrebe po brisanju

            arr = np.array(points)
            lowest_index = np.argmin(arr[:, axis])
            removed = points.pop(lowest_index)
            print(f"⚠️ Odstranjena najnižja točka (os {axis}): {removed}")
            return points
       
        def update_timing_csv(timing_data, study_name):
            file_path = os.path.join(os.path.dirname(__file__), "timing_summary.csv")
            file_exists = os.path.isfile(file_path)

            with open(file_path, mode='a', newline='') as csvfile:
                fieldnames = ["Study", "IO", "Fixing", "Scaling", "CentroidAlign", "Rotation", "Translation", "Transform", "FileSave", "Total"]
                writer = csv.DictWriter(csvfile, fieldnames=fieldnames)

                if not file_exists:
                    writer.writeheader()

                row = {"Study": study_name}
                row.update(timing_data)
                writer.writerow(row)
                
        def find_volume_node_by_partial_name(partial_name):
            for node in slicer.util.getNodesByClass("vtkMRMLScalarVolumeNode"):
                if partial_name in node.GetName():
                    return node
            raise RuntimeError(f"❌ Volume with name containing '{partial_name}' not found.")
                
        # Globalni seznami za končno statistiko
        prostate_size_est = []
        ctcbct_distance = []
        table_z_values = {}

        # Pridobimo SubjectHierarchyNode
        shNode = slicer.vtkMRMLSubjectHierarchyNode.GetSubjectHierarchyNode(slicer.mrmlScene)
        
        studyItems = vtk.vtkIdList()
        shNode.GetItemChildren(shNode.GetSceneItemID(), studyItems)
        
        for i in range(studyItems.GetNumberOfIds()):
            study_start_time = time.time()
            start_io = time.time()
            studyItem = studyItems.GetId(i)
            studyName = shNode.GetItemName(studyItem)
            print(f"\nProcessing study: {studyName}")
            # **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)
                    #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:
                ct_volume_name = ct_list[0]  # Uporabi prvi CT kot referenco
                ct_volume_Node = find_volume_node_by_partial_name(ct_volume_name)
                print(f"\nProcessing CT: {ct_volume_name}")
                yesCbct = False
                makemarkerscheck = True                               
                result = find_table_top_z(ct_volume_name, writefilecheck, makemarkerscheck, yesCbct)
                if result is not None:
                    mm_offset, pixel_offset = result
                    ct_table_found = True
                    #print(f"✔️ Poravnava z višino mize: ΔY = {mm_offset:.2f} mm")
                    # Dodaj ΔY k translaciji ali transformaciji po potrebi
                else:
                    print("⚠️ Table top not found – continue without correction on Y axis.")
                    mm_offset = 0.0  # ali None, če želiš eksplicitno ignorirati
                    ct_table_found = False
                CT_offset, CT_spacing = align_cbct_to_ct(ct_volume_Node, "CT", mm_offset)                                    
                     
                
                ct_points = [centroid for centroid, _ in volume_points_dict[("CT", ct_volume_name)]]
                ct_points = remove_lowest_marker(ct_points) #odstrani marker v riti, če obstaja
                print(f"CT points: {ct_points}")                
                if len(ct_points) < 3:
                    print(f"CT volume {ct_volume_name} doesn't have enough points for registration. Points: {len(ct_points)}")
                    continue
                else:
                    for cbct_volume_name in cbct_list:
                        print(f"\nProcessing CBCT Volume: {cbct_volume_name}")                        
                        yesCbct = True
                        scan_type = "CBCT" #redundant but here we are
                        cbct_volume_node = find_volume_node_by_partial_name(cbct_volume_name)
                        
                        key = (studyName, cbct_volume_name)
                        if key not in cumulative_matrices:
                            cumulative_matrices[key] = np.eye(4)
                        
                        fixing = time.time()
                        makemarkerscheck = False  # Ustvari CBCT miza markerje pred poravnavo
                        if(ct_table_found):
                            result = find_table_top_z(cbct_volume_name, writefilecheck, makemarkerscheck, yesCbct) #!!!!!!!!!!!!!???????????? ct_volume_name
                            if result is not None:
                                mm_offset, pixel_offset = result
                                #print(f"✔️ Poravnava z višino mize: ΔY = {mm_offset:.2f} mm")
                                skupni_offset = align_cbct_to_ct(cbct_volume_node, scan_type, mm_offset, CT_offset, CT_spacing) #poravna CBCT in sporoči skupni offset                                
                                
                                table_shift_matrix = np.eye(4)
                                table_shift_matrix[1, 3] = skupni_offset  # Premik po Y

                                # Pomnoži obstoječo kumulativno matriko
                                key = (studyName, cbct_volume_name)
                                if key not in cumulative_matrices:
                                    cumulative_matrices[key] = np.eye(4)

                                cumulative_matrices[key] = np.dot(cumulative_matrices[key], table_shift_matrix)
                                
                            else:
                                print("⚠️ Table top not found – continue without correction on Y axis.")
                                mm_offset = 0.0             
                                            
                        cbct_points = [centroid for centroid, _ in volume_points_dict[("CBCT", cbct_volume_name)]] #zastareli podatki
                        cbct_points_array = np.array(cbct_points)  # Pretvorba v numpy array
                        
                        # print_orientation(ct_volume_name)
                        # print_orientation(cbct_volume_name)
                        
                        #for i, (cb, ct) in enumerate(zip(cbct_points, ct_points)):
                        #    print(f"Pair {i}: CBCT {cb}, CT {ct}, diff: {np.linalg.norm(cb - ct):.2f}")
                        
                        
                        
                        ustvari_marker = False  # Ustvari markerje
                        cbct_points = [centroid for centroid, _ in detect_points_region_growing(cbct_volume_name, yesCbct, ustvari_marker, intensity_threshold=3000)]                     
                        #cbct_points = detect_points_region_growing(cbct_volume_name, yesCbct, intensity_threshold=3000)   
                        #cbct_points = [centroid for centroid, _ in volume_points_dict[("CBCT", cbct_volume_name)]] #zastareli podatki
                                                
                                                
                        if len(cbct_points) < 3:
                            print(f"CBCT Volume '{cbct_volume_name}' doesn't have enough points for registration. Points: {len(cbct_points)}")
                            continue
                        
                        cbct_spacing = cbct_volume_node.GetSpacing()
                        ct_spacing = ct_volume_Node.GetSpacing()
                        cbct_points = rescale_points_to_match_spacing(cbct_points, cbct_spacing, ct_spacing)
                        
                        #Sortiramo točke po X/Y/Z da se izognemo težavam pri poravnavi
                        cbct_points = match_points(cbct_points, ct_points)
                        fixing_end = time.time()
                        #visualize_point_matches_in_slicer(cbct_points, ct_points, studyName) #poveže pare markerjev
        
                        if writefilecheck:
                            # Shranjevanje razdalj
                            distances_ct_cbct = []
                            distances_internal = {"A-B": [], "B-C": [], "C-A": []}

                            cbct_volume_node = find_volume_node_by_partial_name(cbct_volume_name)
                            
                            # Sortiramo točke po Z-koordinati (ali X/Y, če raje uporabljaš drugo os)
                            cbct_points_sorted = cbct_points_array[np.argsort(cbct_points_array[:, 2])]

                            # Razdalje med CT in CBCT (SORTIRANE točke!)
                            d_ct_cbct = np.linalg.norm(cbct_points_sorted - ct_points, axis=1)
                            distances_ct_cbct.append(d_ct_cbct)

                            # Razdalje med točkami znotraj SORTIRANIH cbct_points
                            d_ab = np.linalg.norm(cbct_points_sorted[0] - cbct_points_sorted[1])
                            d_bc = np.linalg.norm(cbct_points_sorted[1] - cbct_points_sorted[2])
                            d_ca = np.linalg.norm(cbct_points_sorted[2] - cbct_points_sorted[0])

                            # Sortiramo razdalje po velikosti, da so vedno v enakem vrstnem redu
                            sorted_distances = sorted([d_ab, d_bc, d_ca])

                            distances_internal["A-B"].append(sorted_distances[0])
                            distances_internal["B-C"].append(sorted_distances[1])
                            distances_internal["C-A"].append(sorted_distances[2])
                            
                            # Dodamo ime študije za v statistiko
                            studyName = shNode.GetItemName(studyItem)
                            
                            # **Shrani razdalje v globalne sezname**
                            prostate_size_est.append({"Study": studyName, "Distances": sorted_distances})
                            ctcbct_distance.append({"Study": studyName, "Distances": list(distances_ct_cbct[-1])})  # Pretvorimo v seznam

                        # Izberi metodo glede na uporabnikov izbor
                        chosen_rotation_matrix = np.eye(3)
                        chosen_translation_vector = np.zeros(3)
                        #print("Markerji pred transformacijo:", cbct_points, ct_points)
                        start_scaling = time.time()
                        scaling_factors = None
                        if applyScaling:
                            scaling_factors = compute_scaling_stddev(cbct_points, ct_points)
                            #print("Scaling factors: ", scaling_factors)
                            cbct_points = compute_scaling(cbct_points, scaling_factors)
                            
                        
                        
                            
                            
                        end_scaling = time.time()
                        start_align = time.time()
                        initial_error = np.mean(np.linalg.norm(np.array(cbct_points) - np.array(ct_points), axis=1))
                        if initial_error > 30:
                            #print("⚠️ Initial distance too large, applying centroid prealignment.")
                            cbct_points, transvector = prealign_by_centroid(cbct_points, ct_points)
                        else:  
                            transvector = np.zeros(3)   
                        end_align = time.time()
                        
                        start_rotation = time.time()
                        if applyRotation:
                            if selectedMethod == "Kabsch":
                                chosen_rotation_matrix = compute_Kabsch_rotation(cbct_points, ct_points)
                            elif selectedMethod == "Horn":
                                chosen_rotation_matrix = compute_Horn_rotation(cbct_points, ct_points)
                            elif selectedMethod == "Iterative Closest Point (Kabsch)":
                                _, chosen_rotation_matrix, _ = icp_algorithm(cbct_points, ct_points)
                            #print("Rotation Matrix:\n", chosen_rotation_matrix)
                        end_rotation = time.time()

                        
                        start_translation = time.time()
                        fine_shift = np.zeros(3)  # Inicializiraj fine premike
                        if applyTranslation:
                            chosen_translation_vector, cbct_points_transformed, method_used = choose_best_translation(
                                cbct_points, ct_points, chosen_rotation_matrix)
                            
                            # Sistematična razlika (signed shift)
                            rotated_cbct = np.dot(cbct_points, chosen_rotation_matrix.T)
                            translated_cbct = rotated_cbct + chosen_translation_vector
                            delta_y_list = [ct[1] - cbct[1] for ct, cbct in zip(ct_points, translated_cbct)]
                            mean_delta_y = np.mean(delta_y_list)


                            # Uporabi sistematični shift za dodatno poravnavo v y-osi
                            fine_shift = np.array([0.0, mean_delta_y, 0.0])  # samo Y-os
                            cbct_points_transformed += fine_shift

                        end_translation = time.time()

                        start_transform = time.time()
                        # ✅ Kombinirana transformacija
                        total_translation = chosen_translation_vector + fine_shift
                        chosen_translation_vector = total_translation
                        vtk_transform = create_vtk_transform(chosen_rotation_matrix, chosen_translation_vector, studyName, cbct_volume_name, scaling_factors)
                        
                        
                        
                        
                        
                        combined_matrix = np.eye(4)
                        
                        if scaling_factors is not None:
                            scaling_matrix = np.diag([scaling_factors[0], scaling_factors[1], scaling_factors[2], 1.0])
                            combined_matrix = np.dot(scaling_matrix, combined_matrix)
                        
                        if chosen_translation_vector is not None:
                            combined_matrix[:3, 3] = chosen_translation_vector + transvector # glavna poravnava, sistematični y shift in groba poravnava
                        if chosen_rotation_matrix is not None:
                            combined_matrix[:3, :3] = chosen_rotation_matrix
                        
                        
                        cumulative_matrices[(studyName, cbct_volume_name)] = np.dot(cumulative_matrices[(studyName, cbct_volume_name)], combined_matrix)



                        # 🔄 Pripni transformacijo
                        imeTransformNoda = cbct_volume_name + " Transform"
                        transform_node = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLTransformNode", imeTransformNoda)
                        transform_node.SetAndObserveTransformToParent(vtk_transform)

                        cbct_volume_node = 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)
                        
                        #popravek v x osi
                        delta_x_list = [ct[0] - cbct[0] for ct, cbct in zip(ct_points, cbct_points)]
                        mean_delta_x = np.mean(delta_x_list)
                        #popravek v y osi
                        delta_y_list = [ct[1] - cbct[1] for ct, cbct in zip(ct_points, cbct_points)]
                        mean_delta_y = np.mean(delta_y_list)
                        #popravek v z osi
                        delta_z_list = [ct[2] - cbct[2] for ct, cbct in zip(ct_points, cbct_points)]
                        mean_delta_z = np.mean(delta_z_list)

                        # Uporabi sistematični shift za dodatno poravnavo
                        fine_shift = np.array([mean_delta_x, mean_delta_y, mean_delta_z])
                        #cbct_points_transformed += fine_shift
                        
                        if fine_shift is not None:
                            shift_matrix = np.eye(4)
                            shift_matrix[:3, 3] = fine_shift
                            cumulative_matrices[(studyName, cbct_volume_name)] = np.dot(cumulative_matrices[(studyName, cbct_volume_name)], shift_matrix)
                 
                    
                        chosen_rotation_matrix = np.eye(3) #tokrat brez rotacije
                        vtk_transform = create_vtk_transform(chosen_rotation_matrix, fine_shift, studyName, cbct_volume_name) #Tukaj se tudi izpiše transformacijska matrika

                        # 🔄 Pripni transformacijo
                        imeTransformNoda = cbct_volume_name + " Transform"
                        transform_node = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLTransformNode", imeTransformNoda)
                        transform_node.SetAndObserveTransformToParent(vtk_transform)

                        cbct_volume_node = find_volume_node_by_partial_name(cbct_volume_name)
                        cbct_volume_node.SetAndObserveTransformNodeID(transform_node.GetID())

                        # 🔨 Uporabi (ali shrani transformacijo kasneje)
                        slicer.vtkSlicerTransformLogic().hardenTransform(cbct_volume_node)
                        slicer.mrmlScene.RemoveNode(transform_node)

                        
                        ustvari_marker = True
                        cbct_points = [centroid for centroid, _ in detect_points_region_growing(cbct_volume_name, yesCbct, ustvari_marker, intensity_threshold=3000)]
                        print(f"Fine correction shifts: ΔX={fine_shift[0]:.2f} mm, ΔY={fine_shift[1]:.2f} mm, ΔZ={fine_shift[2]:.2f} mm")

                        #shrani transformacijsko matriko v datoteko
                        save_transform_matrix(cumulative_matrices[(studyName, cbct_volume_name)], studyName, cbct_volume_name, ct_table_found)
                        
                        
                        # 📏 Izračun napake
                        errors = [np.linalg.norm(cbct - ct) for cbct, ct in zip(cbct_points, ct_points)]
                        mean_error = np.mean(errors)

                        print("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")
                        
                        
                
            else:
                print(f"Study {studyItem} doesn't have any appropriate CT or CBCT volumes.")
                continue
            
            study_end_time = time.time()    
            timing_data = {
                "IO": end_io - start_io,
                "Fixing": fixing_end - fixing,
                "Scaling": end_scaling - start_scaling,
                "CentroidAlign": end_align - start_align,
                "Rotation": end_rotation - start_rotation,
                "Translation": end_translation - start_translation,
                "Transform": end_transform - start_transform,
                "Total": study_end_time - study_start_time}
            update_timing_csv(timing_data, studyName)
            print(f"Timing data for {studyName}: {timing_data}")

        # Izpis globalne statistike
        
        start_save = time.time()
        if writefilecheck:
            #print("Distances between CT & CBCT markers: ", ctcbct_distance)
            #print("Distances between pairs of markers for each volume: ", prostate_size_est)
            
            # 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")