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