import os import numpy as np import scipy import re from scipy.spatial.distance import cdist from scipy.spatial.transform import Rotation as R import slicer import slicer.util import itertools import DICOMLib from DICOMLib import DICOMUtils import DicomRtImportExportPlugin from collections import deque, Counter import vtk from slicer.ScriptedLoadableModule import * import qt from datetime import datetime import csv import time import logging import matplotlib matplotlib.use('Agg') # << to dodaš ZGORAJ, da omogoči PNG zapis brez GUI import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D #exec(open("C:/Users/lkomar/Documents/Prostata/FirstTryRegister.py").read()) cumulative_matrices = {} class SeekTransformModule(ScriptedLoadableModule): """ Module description shown in the module panel. """ def __init__(self, parent): ScriptedLoadableModule.__init__(self, parent) self.parent.title = "Seek Transform module" self.parent.categories = ["Image Processing"] self.parent.contributors = ["Luka Komar (Onkološki Inštitut Ljubljana, Fakulteta za Matematiko in Fiziko Ljubljana)"] self.parent.helpText = "This module applies rigid transformations to CBCT volumes based on reference CT volumes." self.parent.acknowledgementText = "Supported by doc. Primož Peterlin & prof. Andrej Studen" class SeekTransformModuleWidget(ScriptedLoadableModuleWidget): """ GUI of the module. """ def setup(self): ScriptedLoadableModuleWidget.setup(self) # Dropdown menu za izbiro metode self.rotationMethodComboBox = qt.QComboBox() self.rotationMethodComboBox.addItems(["Kabsch", "Horn", "Iterative Closest Point (Horn)"]) self.layout.addWidget(self.rotationMethodComboBox) # Checkboxi za transformacije self.scalingCheckBox = qt.QCheckBox("Scaling") self.scalingCheckBox.setChecked(False) self.layout.addWidget(self.scalingCheckBox) self.rotationCheckBox = qt.QCheckBox("Rotation") self.rotationCheckBox.setChecked(True) self.layout.addWidget(self.rotationCheckBox) self.translationCheckBox = qt.QCheckBox("Translation") self.translationCheckBox.setChecked(True) self.layout.addWidget(self.translationCheckBox) self.markersCheckBox = qt.QCheckBox("Place control points for detected markers") self.markersCheckBox.setChecked(False) self.layout.addWidget(self.markersCheckBox) self.writefileCheckBox = qt.QCheckBox("Write data to csv file") self.writefileCheckBox.setChecked(True) self.layout.addWidget(self.writefileCheckBox) self.tableCheckBox = qt.QCheckBox("Find top of the table and match height") self.tableCheckBox.setChecked(True) self.layout.addWidget(self.tableCheckBox) self.FineShiftCheckBox = qt.QCheckBox("Use extra fine shift correction") self.FineShiftCheckBox.setChecked(True) self.layout.addWidget(self.FineShiftCheckBox) self.save_as_dicomCheckBox = qt.QCheckBox("Save transformed CT and segmentations as DICOM files") self.save_as_dicomCheckBox.setChecked(False) self.layout.addWidget(self.save_as_dicomCheckBox) # Load button self.applyButton = qt.QPushButton("Find markers and transform") self.applyButton.toolTip = "Finds markers, computes optimal rigid transform and applies it to CBCT volumes." self.applyButton.enabled = True self.layout.addWidget(self.applyButton) # Connect button to logic self.applyButton.connect('clicked(bool)', self.onApplyButton) self.layout.addStretch(1) def onApplyButton(self): # Nastavi globalni logger log_file_path = os.path.join("C:/Users/lkomar/Documents/Prostata", "seektransform_log.txt") logging.basicConfig(filename=log_file_path, level=logging.INFO, format='%(asctime)s - %(message)s') try: logging.info("▶️ onApplyButton pressed.") except Exception as e: print("❌ Logging setup failed:", e) logic = MyTransformModuleLogic() selectedMethod = self.rotationMethodComboBox.currentText # izberi metodo izračuna rotacije # Preberi stanje checkboxov applyRotation = self.rotationCheckBox.isChecked() applyTranslation = self.translationCheckBox.isChecked() applyScaling = self.scalingCheckBox.isChecked() applyMarkers = self.markersCheckBox.isChecked() writefilecheck = self.writefileCheckBox.isChecked() tablefind = self.tableCheckBox.isChecked() use_fine_shift = self.FineShiftCheckBox.isChecked() save_as_dicom = self.save_as_dicomCheckBox.isChecked() # Pokliči logiko z izbranimi nastavitvami logic.run(selectedMethod, applyRotation, applyTranslation, applyScaling, applyMarkers, writefilecheck, tablefind, use_fine_shift, save_as_dicom) class MyTransformModuleLogic(ScriptedLoadableModuleLogic): """ Core logic of the module. """ def run(self, selectedMethod, applyRotation, applyTranslation, applyScaling, applymarkers, writefilecheck, tablefind, use_fine_shift, save_as_dicom): start_time = time.time() print("Calculating...") #slicer.util.delayDisplay(f"Starting", 1000) def group_points(points, threshold): # Function to group points that are close to each other grouped_points = [] while points: point = points.pop() # Take one point from the list group = [point] # Start a new group # Find all points close to this one distances = cdist([point], points) # Calculate distances from this point to others close_points = [i for i, dist in enumerate(distances[0]) if dist < threshold] # Add the close points to the group group.extend([points[i] for i in close_points]) # Remove the grouped points from the list points = [point for i, point in enumerate(points) if i not in close_points] # Add the group to the result grouped_points.append(group) return grouped_points def region_growing(image_data, seed, intensity_threshold, max_distance): dimensions = image_data.GetDimensions() visited = set() region = [] queue = deque([seed]) while queue: x, y, z = queue.popleft() if (x, y, z) in visited: continue visited.add((x, y, z)) voxel_value = image_data.GetScalarComponentAsDouble(x, y, z, 0) if voxel_value >= intensity_threshold: region.append((x, y, z)) # Add neighbors within bounds for dx, dy, dz in [(1, 0, 0), (-1, 0, 0), (0, 1, 0), (0, -1, 0), (0, 0, 1), (0, 0, -1)]: nx, ny, nz = x + dx, y + dy, z + dz if 0 <= nx < dimensions[0] and 0 <= ny < dimensions[1] and 0 <= nz < dimensions[2]: if (nx, ny, nz) not in visited: queue.append((nx, ny, nz)) return region def compute_scaling_stddev(moving_points, fixed_points): moving = np.array(moving_points) fixed = np.array(fixed_points) # Standard deviation around centroid, po osi scaling_factors = np.std(fixed, axis=0) / np.std(moving, axis=0) return tuple(scaling_factors) def compute_scaling(cbct_points, scaling_factors): """Applies non-uniform scaling to CBCT points. Args: cbct_points (list of lists): List of (x, y, z) points. scaling_factors (tuple): Scaling factors (sx, sy, sz) for each axis. Returns: np.ndarray: Scaled CBCT points. """ sx, sy, sz = scaling_factors # Extract scaling factors scaling_matrix = np.diag([sx, sy, sz]) # Create diagonal scaling matrix cbct_points_np = np.array(cbct_points) # Convert to numpy array scaled_points = cbct_points_np @ scaling_matrix.T # Apply scaling scaling_4x4 = np.eye(4) scaling_4x4[0, 0] = sx scaling_4x4[1, 1] = sy scaling_4x4[2, 2] = sz return scaled_points.tolist() # Convert back to list def compute_Kabsch_rotation(moving_points, fixed_points): """ Computes the optimal rotation matrix to align moving_points to fixed_points. Parameters: moving_points (list or ndarray): List of points to be rotated CBCT fixed_points (list or ndarray): List of reference points CT Returns: ndarray: Optimal rotation matrix. """ assert len(moving_points) == len(fixed_points), "Point lists must be the same length." # Convert to numpy arrays moving = np.array(moving_points) fixed = np.array(fixed_points) # Compute centroids centroid_moving = np.mean(moving, axis=0) centroid_fixed = np.mean(fixed, axis=0) # Center the points moving_centered = moving - centroid_moving fixed_centered = fixed - centroid_fixed # Compute covariance matrix H = np.dot(moving_centered.T, fixed_centered) # SVD decomposition U, _, Vt = np.linalg.svd(H) Rotate_optimal = np.dot(Vt.T, U.T) # Correct improper rotation (reflection) if np.linalg.det(Rotate_optimal) < 0: Vt[-1, :] *= -1 Rotate_optimal = np.dot(Vt.T, U.T) return Rotate_optimal def compute_Horn_rotation(moving_points, fixed_points): """ Computes the optimal rotation matrix using quaternions. Parameters: moving_points (list or ndarray): List of points to be rotated. fixed_points (list or ndarray): List of reference points. Returns: ndarray: Optimal rotation matrix. """ assert len(moving_points) == len(fixed_points), "Point lists must be the same length." moving = np.array(moving_points) fixed = np.array(fixed_points) # Compute centroids centroid_moving = np.mean(moving, axis=0) centroid_fixed = np.mean(fixed, axis=0) # Center the points moving_centered = moving - centroid_moving fixed_centered = fixed - centroid_fixed # Construct the cross-dispersion matrix M = np.dot(moving_centered.T, fixed_centered) # Construct the N matrix for quaternion solution A = M - M.T delta = np.array([A[1, 2], A[2, 0], A[0, 1]]) trace = np.trace(M) N = np.zeros((4, 4)) N[0, 0] = trace N[1:, 0] = delta N[0, 1:] = delta N[1:, 1:] = M + M.T - np.eye(3) * trace # Compute the eigenvector corresponding to the maximum eigenvalue eigvals, eigvecs = np.linalg.eigh(N) q_optimal = eigvecs[:, np.argmax(eigvals)] # Optimal quaternion # Convert quaternion to rotation matrix w, x, y, z = q_optimal R = np.array([ [1 - 2*(y**2 + z**2), 2*(x*y - z*w), 2*(x*z + y*w)], [2*(x*y + z*w), 1 - 2*(x**2 + z**2), 2*(y*z - x*w)], [2*(x*z - y*w), 2*(y*z + x*w), 1 - 2*(x**2 + y**2)] ]) return R def icp_algorithm(moving_points, fixed_points, max_iterations=100, tolerance=1e-5): """ Iterative Closest Point (ICP) algorithm to align moving_points to fixed_points. Parameters: moving_points (list or ndarray): List of points to be aligned. fixed_points (list or ndarray): List of reference points. max_iterations (int): Maximum number of iterations. tolerance (float): Convergence tolerance. Returns: ndarray: Transformed moving points. ndarray: Optimal rotation matrix. ndarray: Optimal translation vector. """ # Convert to numpy arrays moving = np.array(moving_points) fixed = np.array(fixed_points) # Initialize transformation R = np.eye(3) # Identity matrix for rotation t = np.zeros(3) # Zero vector for translation prev_error = np.inf # Initialize previous error to a large value for iteration in range(max_iterations): # Step 1: Find the nearest neighbors (correspondences) distances = np.linalg.norm(moving[:, np.newaxis] - fixed, axis=2) nearest_indices = np.argmin(distances, axis=1) nearest_points = fixed[nearest_indices] # Step 2: Compute the optimal rotation and translation R_new = compute_Horn_rotation(moving, nearest_points) centroid_moving = np.mean(moving, axis=0) centroid_fixed = np.mean(nearest_points, axis=0) t_new = centroid_fixed - np.dot(R_new, centroid_moving) # Step 3: Apply the transformation moving = np.dot(moving, R_new.T) + t_new # Update the cumulative transformation R = np.dot(R_new, R) t = np.dot(R_new, t) + t_new # Step 4: Check for convergence mean_error = np.mean(np.linalg.norm(moving - nearest_points, axis=1)) if np.abs(prev_error - mean_error) < tolerance: print(f"ICP converged after {iteration + 1} iterations.") break prev_error = mean_error else: print(f"ICP reached maximum iterations ({max_iterations}).") return moving, R, t def match_points(cbct_points, ct_points, auto_weights=False, fallback_if_worse=False, normalize_lengths=True, normalize_angles=False, min_distance=5, w_order=1.0): def side_lengths(points): lengths = [ np.linalg.norm(points[0] - points[1]), np.linalg.norm(points[1] - points[2]), np.linalg.norm(points[2] - points[0]) ] return lengths def triangle_angles(points): a = np.linalg.norm(points[1] - points[2]) b = np.linalg.norm(points[0] - points[2]) c = np.linalg.norm(points[0] - points[1]) angle_A = np.arccos(np.clip((b**2 + c**2 - a**2) / (2 * b * c), -1.0, 1.0)) angle_B = np.arccos(np.clip((a**2 + c**2 - b**2) / (2 * a * c), -1.0, 1.0)) angle_C = np.pi - angle_A - angle_B return [angle_A, angle_B, angle_C] def normalize(vec): norm = np.linalg.norm(vec) return [v / norm for v in vec] if norm > 0 else vec def permutation_score(perm, ct_lengths, ct_angles, w_len, w_ang, penalty_angle_thresh=np.deg2rad(10)): perm_lengths = side_lengths(perm) perm_angles = triangle_angles(perm) # Filter za minimum razdalje if min(perm_lengths) < min_distance: return float('inf') lengths_1 = normalize(perm_lengths) if normalize_lengths else perm_lengths lengths_2 = normalize(ct_lengths) if normalize_lengths else ct_lengths angles_1 = normalize(perm_angles) if normalize_angles else perm_angles angles_2 = normalize(ct_angles) if normalize_angles else ct_angles score_len = sum(abs(a - b) for a, b in zip(lengths_1, lengths_2)) score_ang = sum(abs(a - b) for a, b in zip(angles_1, angles_2)) order_penalty = order_mismatch_penalty(ct_points, perm, axis='z') return w_len * score_len + w_ang * score_ang + order_penalty + w_order * order_penalty def smart_sort_cbct_points(cbct_points, z_threshold=5.0): """ Sortira točke tako, da poskusi najprej po Z. Če so razlike po Z manjše od praga, sortira po (Y, X), sicer sortira po Z. """ z_values = [pt[2] for pt in cbct_points] z_range = max(z_values) - min(z_values) if z_range < z_threshold: # Sortiraj po Y, nato X (če so točke v isti ravnini po Z) return sorted(cbct_points, key=lambda pt: (pt[1], pt[0])) else: # Sortiraj po Z, nato Y, nato X return sorted(cbct_points, key=lambda pt: (pt[2], pt[1], pt[0])) def order_mismatch_penalty(ct_points, perm, axis='z'): axis_idx = {'x': 0, 'y': 1, 'z': 2}[axis] ct_sorted = np.argsort([pt[axis_idx] for pt in ct_points]) perm_sorted = np.argsort([pt[axis_idx] for pt in perm]) return sum(1 for a, b in zip(ct_sorted, perm_sorted) if a != b) cbct_points = list(cbct_points) print("CBCT points:", cbct_points) ct_lengths = side_lengths(np.array(ct_points)) ct_angles = triangle_angles(np.array(ct_points)) if auto_weights: var_len = np.var(ct_lengths) var_ang = np.var(ct_angles) total_var = var_len + var_ang + 1e-6 weight_length = (1 - var_len / total_var) weight_angle = (1 - var_ang / total_var) else: weight_length = 0.8 weight_angle = 0.2 cbct_sorted = smart_sort_cbct_points(cbct_points) original_score = permutation_score(np.array(cbct_sorted), ct_lengths, ct_angles, weight_length, weight_angle) # Če je ta rezultat dovolj dober, uporabi best_score = float('inf') best_perm = None if original_score < float('inf'): # lahko dodaš prag če želiš best_score = original_score best_perm = np.array(cbct_sorted) # Nato preveri vse permutacije (vključno s prvotnim vrstnim redom, če fallback_if_worse=True) for perm in itertools.permutations(cbct_points): perm = np.array(perm) score = permutation_score(perm, ct_lengths, ct_angles, weight_length, weight_angle) if score < best_score: best_score = score best_perm = perm print(f"New best permutation found with perm: {perm}") #print("CT centroid:", np.mean(ct_points, axis=0)) #print("CBCT centroid (best perm):", np.mean(best_perm, axis=0)) if fallback_if_worse: #original_score = permutation_score(np.array(cbct_points), ct_lengths, ct_angles, weight_length, weight_angle) print("Original score: ", original_score) if original_score <= best_score: print("Fallback to original points due to worse score of the permutation.") return list(cbct_points) return list(best_perm) def compute_translation(moving_points, fixed_points, rotation_matrix): """ Computes the translation vector to align moving_points to fixed_points given a rotation matrix. Parameters: moving_points (list or ndarray): List of points to be translated. fixed_points (list or ndarray): List of reference points. rotation_matrix (ndarray): Rotation matrix. Returns: ndarray: Translation vector. """ # Convert to numpy arrays moving = np.array(moving_points) fixed = np.array(fixed_points) # Compute centroids centroid_moving = np.mean(moving, axis=0) centroid_fixed = np.mean(fixed, axis=0) # Compute translation translation = centroid_fixed - np.dot(centroid_moving, rotation_matrix) return translation def create_vtk_transform(rotation_matrix, translation_vector, tablefound, study_name=None, cbct_volume_name=None, scaling_factors=None): """ Creates a vtkTransform from scaling, rotation, and translation. Shrani tudi kumulativno matriko v globalni slovar cumulative_matrices. """ # ----- Inicializacija ----- global cumulative_matrices transform = vtk.vtkTransform() # ----- 1. Skaliranje ----- if scaling_factors is not None: sx, sy, sz = scaling_factors transform.Scale(sx, sy, sz) # ----- 2. Rotacija ----- # Rotacijsko matriko in translacijo pretvori v homogeno matriko affine_matrix = np.eye(4) affine_matrix[:3, :3] = rotation_matrix affine_matrix[:3, 3] = translation_vector # Vstavi v vtkMatrix4x4 vtk_matrix = vtk.vtkMatrix4x4() for i in range(4): for j in range(4): vtk_matrix.SetElement(i, j, affine_matrix[i, j]) transform.Concatenate(vtk_matrix) # # ----- 3. Debug izpis ----- # print("Transform matrix:") # for i in range(4): # print(" ".join(f"{vtk_matrix.GetElement(i, j):.6f}" for j in range(4))) # ----- 4. Shrani v kumulativni matriki ----- if study_name and cbct_volume_name: key = (study_name, cbct_volume_name) if key not in cumulative_matrices: cumulative_matrices[key] = np.eye(4) cumulative_matrices[key] = np.dot(cumulative_matrices[key], affine_matrix) return transform def save_transform_matrix(matrix, study_name, cbct_volume_name): """ Appends the given 4x4 matrix to a text file under the given study folder. """ base_folder = os.path.join(os.path.dirname(__file__), "Transformacijske matrike") study_folder = os.path.join(base_folder, study_name) os.makedirs(study_folder, exist_ok=True) # Create folders if they don't exist safe_cbct_name = re.sub(r'[<>:"/\\|?*]', '_', cbct_volume_name) # Preveri ali je CT miza najdena filename = os.path.join(study_folder, f"{safe_cbct_name}.txt") with open(filename, "w") as f: #f.write("Transformacija:\n") for row in matrix: f.write(" ".join(f"{elem:.6f}" for elem in row) + "\n") print(" ".join(f"{elem:.6f}" for elem in row) + "\n") f.write("\n") # Dodaj prazen vrstico med transformacijami #print(f"Transform matrix saved to {filename}") def detect_points_region_growing(volume_name, yesCbct, create_marker, intensity_threshold=3000, x_min=90, x_max=380, y_min=200, y_max=380, z_min=25, z_max=140, max_distance=9, centroid_merge_threshold=5): volume_node = find_volume_node_by_partial_name(volume_name) if not volume_node: raise RuntimeError(f"Volume {volume_name} not found.") image_data = volume_node.GetImageData() matrix = vtk.vtkMatrix4x4() volume_node.GetIJKToRASMatrix(matrix) dimensions = image_data.GetDimensions() #detected_regions = [] if yesCbct: #je cbct ali ct? valid_x_min, valid_x_max = 0, dimensions[0] - 1 valid_y_min, valid_y_max = 0, dimensions[1] - 1 valid_z_min, valid_z_max = 0, dimensions[2] - 1 else: valid_x_min, valid_x_max = max(x_min, 0), min(x_max, dimensions[0] - 1) valid_y_min, valid_y_max = max(y_min, 0), min(y_max, dimensions[1] - 1) valid_z_min, valid_z_max = max(z_min, 0), min(z_max, dimensions[2] - 1) visited = set() def grow_region(x, y, z): if (x, y, z) in visited: return None voxel_value = image_data.GetScalarComponentAsDouble(x, y, z, 0) if voxel_value < intensity_threshold: return None region = region_growing(image_data, (x, y, z), intensity_threshold, max_distance=max_distance) if region: for point in region: visited.add(tuple(point)) return region return None regions = [] for z in range(valid_z_min, valid_z_max + 1): for y in range(valid_y_min, valid_y_max + 1): for x in range(valid_x_min, valid_x_max + 1): region = grow_region(x, y, z) if region: regions.append(region) # Collect centroids using intensity-weighted average centroids = [] for region in regions: points = np.array([matrix.MultiplyPoint([*point, 1])[:3] for point in region]) intensities = np.array([image_data.GetScalarComponentAsDouble(*point, 0) for point in region]) if intensities.sum() > 0: weighted_centroid = np.average(points, axis=0, weights=intensities) max_intensity = intensities.max() centroids.append((np.round(weighted_centroid, 2), max_intensity)) unique_centroids = [] for centroid, intensity in centroids: if not any(np.linalg.norm(centroid - existing_centroid) < centroid_merge_threshold for existing_centroid, _ in unique_centroids): unique_centroids.append((centroid, intensity)) if create_marker: markups_node = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLMarkupsFiducialNode", f"Markers_{volume_name}") for centroid, intensity in unique_centroids: markups_node.AddControlPoint(*centroid) markups_node.SetDisplayVisibility(False) #print(f"Detected Centroid (RAS): {centroid}, Max Intensity: {intensity}") return unique_centroids def find_table_top_z(volume_name, writefilecheck, makemarkerscheck, yesCbct): """ Najde višino zgornjega roba mize v CT/CBCT volumnu in po želji doda marker v sceno. Args: ct_volume_name (str): Ime volumna. writefilecheck (bool): Ali naj se rezultat shrani v .csv. makemarkerscheck (bool): Ali naj se doda marker v 3D Slicer. yesCbct (bool): True, če je CBCT; False, če je CT. Returns: (float, int): Z komponenta v RAS prostoru, in Y indeks v slicerjevem volumnu. """ # --- Pridobi volume node --- volume_node = find_volume_node_by_partial_name(volume_name) np_array = slicer.util.arrayFromVolume(volume_node) # (Z, Y, X) ijkToRasMatrix = vtk.vtkMatrix4x4() volume_node.GetIJKToRASMatrix(ijkToRasMatrix) # --- Določimo lokacijo stolpca --- z_index = np_array.shape[0] // 2 # srednji slice y_size = np_array.shape[1] # x_index = int(np_array.shape[2] * 0.15) # x_index = max(0, min(x_index, np_array.shape[2] - 1)) # --- Izračun spodnje tretjine (spodnji del slike) --- y_start = int(y_size * 2 / 3) slice_data = np_array[z_index, :, :] # (Y, X) y_end = y_size # Do dna slike #column_values = slice_data[y_start:y_end, x_index] # (Y) # --- Parametri za rob --- threshold_high = -300 if yesCbct else -100 threshold_low = -700 if yesCbct else -350 min_jump = 100 if yesCbct else 100 window_size = 4 # število voxelov nad/pod #previous_value = column_values[-1] table_top_y = None # --- Več stolpcev okoli x_index --- x_center = np_array.shape[2] // 2 x_offset = 30 # 30 levo od sredine x_index_base = max(0, x_center - x_offset) candidate_y_values = [] search_range = range(-5, 6) # od -5 do +5 stolpcev for dx in search_range: x_index = x_index_base + dx if x_index < 0 or x_index >= np_array.shape[2]: continue column_values = slice_data[y_start:y_end, x_index] for i in range(window_size, len(column_values) - window_size): curr = column_values[i] above_avg = np.mean(column_values[i - window_size:i]) below_avg = np.mean(column_values[i + 1:i + 1 + window_size]) if (threshold_low < curr < threshold_high and (above_avg - below_avg) > min_jump and below_avg < -400 and above_avg > -300): y_found = y_start + i candidate_y_values.append(y_found) break # samo prvi zadetek v stolpcu if candidate_y_values: most_common_y, _ = Counter(candidate_y_values).most_common(1)[0] table_top_y = most_common_y print(f"candidate_y_values: {candidate_y_values}") print(f"✅ Rob mize (najpogostejši Y): {table_top_y}, pojavitev: {candidate_y_values.count(table_top_y)}/11") """ # --- Poišči skok navzdol pod prag (od spodaj navzgor) --- for i in range(len(column_values) - 2, -1, -1): # od spodaj proti vrhu intensity = column_values[i] if (intensity - previous_value) > min_jump and intensity < thresholdhigh and intensity > thresholdlow: table_top_y = y_start + i - 1 print(f"✅ Rob mize najden pri Y = {table_top_y}, intenziteta = {intensity}") print("Column values (partial):", column_values.tolist()) break previous_value = intensity """ if table_top_y is None: print(f"⚠️ Rob mize ni bil najden (X = {x_index})") print("Column values (partial):", column_values.tolist()) return None # --- Pretvorba v RAS koordinato --- table_ijk = [x_index, table_top_y, z_index] table_ras = np.array(ijkToRasMatrix.MultiplyPoint([*table_ijk, 1]))[:3] # --- Marker --- if makemarkerscheck: table_node = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLMarkupsFiducialNode", f"VišinaMize_{ct_volume_name}") table_node.AddControlPoint(table_ras) table_node.SetDisplayVisibility(False) # --- Shrani v CSV --- if writefilecheck: height_file = os.path.join(os.path.dirname(__file__), "heightdata.csv") with open(height_file, mode='a', newline='') as file: writer = csv.writer(file) modality = "CBCT" if yesCbct else "CT" writer.writerow([modality, ct_volume_name, f"Upper table edge at Z = {table_ras[1]:.2f} mm"]) return table_ras[1], table_top_y def align_cbct_to_ct(volumeNode, scan_type, offset, CT_offset=None, CT_spacing=None): """ Aligns CBCT volume to CT volume based on height offset. Args: volumeNode (vtkMRMLScalarVolumeNode): The volume node to be aligned. scan_type (str): The type of scan ("CT" or "CBCT"). offset (float): The height offset of the current volume from the center in mm. CT_offset (float, optional): The height offset of the CT volume from the center. Required for CBCT alignment. CT_spacing (float, optional): The voxel spacing of the CT volume in mm (for scaling the offset). Returns: float: The alignment offset applied to the CBCT volume (if applicable). """ if scan_type == "CT": CT_offset = offset CT_spacing = volumeNode.GetSpacing()[1] #print(f"CT offset set to: {CT_offset}, CT spacing: {CT_spacing} mm/voxel") return CT_offset, CT_spacing else: if CT_offset is None or CT_spacing is None: raise ValueError("CT_offset and CT_spacing must be provided to align CBCT to CT.") CBCT_offset = offset # Razlika v mm brez skaliranja na CBCT_spacing alignment_offset_mm = CT_offset - CBCT_offset #print(f"CT offset: {CT_offset}, CBCT offset: {CBCT_offset}") #print(f"CT spacing: {CT_spacing} mm/voxel, CBCT spacing: {volumeNode.GetSpacing()[1]} mm/voxel") #print(f"Aligning CBCT with CT. Offset in mm: {alignment_offset_mm}") # Uporabi transformacijo transform = vtk.vtkTransform() transform.Translate(0, alignment_offset_mm, 0) transformNode = slicer.vtkMRMLTransformNode() slicer.mrmlScene.AddNode(transformNode) transformNode.SetAndObserveTransformToParent(transform) volumeNode.SetAndObserveTransformNodeID(transformNode.GetID()) slicer.vtkSlicerTransformLogic().hardenTransform(volumeNode) slicer.mrmlScene.RemoveNode(transformNode) # Poskusi najti ustrezen marker in ga premakniti marker_name = f"VišinaMize_{volumeNode.GetName()}" # Robustno iskanje markerja po imenu table_node = None for node in slicer.util.getNodesByClass("vtkMRMLMarkupsFiducialNode"): if node.GetName() == marker_name: table_node = node break if table_node is not None: current_point = [0, 0, 0] table_node.GetNthControlPointPosition(0, current_point) moved_point = [ current_point[0], current_point[1] + alignment_offset_mm, current_point[2] ] table_node.SetNthControlPointPosition(0, *moved_point) return alignment_offset_mm def print_orientation(volume_name): node = find_volume_node_by_partial_name(volume_name) matrix = vtk.vtkMatrix4x4() node.GetIJKToRASMatrix(matrix) print(f"{volume_name} IJK→RAS:") for i in range(3): print([matrix.GetElement(i, j) for j in range(3)]) def prealign_by_centroid(cbct_points, ct_points): """ Predporavna CBCT markerje na CT markerje glede na centrične točke. Args: cbct_points: List ali ndarray točk iz CBCT. ct_points: List ali ndarray točk iz CT. Returns: List: CBCT točke premaknjene tako, da so centrične točke usklajene. """ cbct_points = np.array(cbct_points) ct_points = np.array(ct_points) cbct_centroid = np.mean(cbct_points, axis=0) ct_centroid = np.mean(ct_points, axis=0) translation_vector = ct_centroid - cbct_centroid aligned_cbct = cbct_points + translation_vector return aligned_cbct, translation_vector def choose_best_translation(cbct_points, ct_points, rotation_matrix): """ Izbere boljšo translacijo: centroidno ali povprečno po rotaciji (retranslation). Args: cbct_points (array-like): Točke iz CBCT (še ne rotirane). ct_points (array-like): Ciljne CT točke. rotation_matrix (ndarray): Rotacijska matrika. Returns: tuple: (best_translation_vector, transformed_cbct_points, used_method) """ cbct_points = np.array(cbct_points) ct_points = np.array(ct_points) # 1. Rotiraj CBCT točke rotated_cbct = np.dot(cbct_points, rotation_matrix.T) # 2. Centroid translacija centroid_moving = np.mean(cbct_points, axis=0) centroid_fixed = np.mean(ct_points, axis=0) translation_centroid = centroid_fixed - np.dot(centroid_moving, rotation_matrix) transformed_centroid = rotated_cbct + translation_centroid error_centroid = np.mean(np.linalg.norm(transformed_centroid - ct_points, axis=1)) # 3. Retranslacija (srednja razlika) translation_recomputed = np.mean(ct_points - rotated_cbct, axis=0) transformed_recomputed = rotated_cbct + translation_recomputed error_recomputed = np.mean(np.linalg.norm(transformed_recomputed - ct_points, axis=1)) # 4. Izberi boljšo if error_recomputed < error_centroid: #print(f"✅ Using retranslation (error: {error_recomputed:.2f} mm)") return translation_recomputed, transformed_recomputed, "retranslation" else: #print(f"✅ Using centroid-based translation (error: {error_centroid:.2f} mm)") return translation_centroid, transformed_centroid, "centroid" def rescale_points_to_match_spacing(points, source_spacing, target_spacing): scale_factors = np.array(target_spacing) / np.array(source_spacing) return np.array(points) * scale_factors def visualize_point_matches_in_slicer(cbct_points, ct_points, study_name="MatchVisualization"): assert len(cbct_points) == len(ct_points), "Mora biti enako število točk!" # Ustvari markups za CBCT cbct_node = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLMarkupsFiducialNode", f"{study_name}_CBCT") cbct_node.GetDisplayNode().SetSelectedColor(0, 0, 1) # modra # Ustvari markups za CT ct_node = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLMarkupsFiducialNode", f"{study_name}_CT") ct_node.GetDisplayNode().SetSelectedColor(1, 0, 0) # rdeča # Dodaj točke for i, (cbct, ct) in enumerate(zip(cbct_points, ct_points)): cbct_node.AddControlPoint(*cbct, f"CBCT_{i}") ct_node.AddControlPoint(*ct, f"CT_{i}") # Ustvari model z linijami med pari points = vtk.vtkPoints() lines = vtk.vtkCellArray() for i, (p1, p2) in enumerate(zip(cbct_points, ct_points)): id1 = points.InsertNextPoint(p1) id2 = points.InsertNextPoint(p2) line = vtk.vtkLine() line.GetPointIds().SetId(0, id1) line.GetPointIds().SetId(1, id2) lines.InsertNextCell(line) polyData = vtk.vtkPolyData() polyData.SetPoints(points) polyData.SetLines(lines) # Model node modelNode = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLModelNode", f"{study_name}_Connections") modelNode.SetAndObservePolyData(polyData) modelDisplay = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLModelDisplayNode") modelDisplay.SetColor(0, 0, 0) # črna modelDisplay.SetLineWidth(2) modelDisplay.SetVisibility(True) modelNode.SetAndObserveDisplayNodeID(modelDisplay.GetID()) modelNode.SetAndObservePolyData(polyData) print(f"✅ Vizualizacija dodana za {study_name} (točke + povezave)") def remove_lowest_marker(points, axis=1): """ Odstrani outlier: točko z največjo 3D razdaljo do robustnega centroida (median). Združljivo: 'points' je lahko [XYZ] ali [(XYZ, intensity)]. :param points: list[np.array([x,y,z])] ali list[(np.array([x,y,z]), meta)] :return: isti list, z odstranjenim outlierjem """ import numpy as np if not points or len(points) <= 3: # nič ne odstrani, potrebujemo vsaj 4 kandidate (npr. 3 markerji + 1 koža) return points # Normaliziraj vhod v Nx3 matriko koordinat + obdrži indeksno preslikavo coords = [] for p in points: if isinstance(p, (list, tuple)) and len(p) == 2 and hasattr(p[0], "__len__"): xyz = np.asarray(p[0], dtype=float) # (centroid, intensity) else: xyz = np.asarray(p, dtype=float) # samo centroid coords.append(xyz[:3]) A = np.vstack(coords) # (N,3) # Robustni centroid center = np.median(A, axis=0) # Evklidske razdalje (kvadrati zadoščajo za max) d2 = np.sum((A - center) ** 2, axis=1) idx_out = int(np.argmax(d2)) removed_point = A[idx_out] removed = points.pop(idx_out) print( f"⚠️ Odstranjen outlier (najbolj oddaljen od median-centroida): " f"{removed_point.round(2)} | d={float(np.sqrt(d2[idx_out])):.2f} mm | " f"centroid≈{center.round(2)}" ) return points def update_timing_csv(timing_data, study_name): file_path = os.path.join(os.path.dirname(__file__), "timing_summary.csv") file_exists = os.path.isfile(file_path) with open(file_path, mode='a', newline='') as csvfile: fieldnames = ["Study", "IO", "Fixing", "Table", "Scaling", "CentroidAlign", "Rotation", "Translation", "Transform", "FileSave", "Total"] writer = csv.DictWriter(csvfile, fieldnames=fieldnames) if not file_exists: writer.writeheader() row = {"Study": study_name} row.update(timing_data) writer.writerow(row) def find_volume_node_by_partial_name(partial_name): for node in slicer.util.getNodesByClass("vtkMRMLScalarVolumeNode"): if partial_name in node.GetName(): return node raise RuntimeError(f"❌ Volume with name containing '{partial_name}' not found.") def convert_rtstruct_to_segmentation_nodes(ct_volume_node): shNode = slicer.vtkMRMLSubjectHierarchyNode.GetSubjectHierarchyNode(slicer.mrmlScene) segmentation_nodes = [] # ✅ Najdi vse SegmentationNode-e, ki imajo segmente for seg_node in slicer.util.getNodesByClass("vtkMRMLSegmentationNode"): num_segments = seg_node.GetSegmentation().GetNumberOfSegments() if num_segments == 0: continue # Če še nima referenceVolume, jo nastavimo if not seg_node.GetNodeReferenceID("referenceVolume"): seg_node.SetReferenceImageGeometryParameterFromVolumeNode(ct_volume_node) seg_node.SetNodeReferenceID("referenceVolume", ct_volume_node.GetID()) print(f"🔗 Nastavljena referenca na CT za: {seg_node.GetName()}") segmentation_nodes.append(seg_node) print(f"📎 Najdena segmentacija z vsebino: {seg_node.GetName()}, segmentov: {num_segments}") if segmentation_nodes: return segmentation_nodes for item_id in range(shNode.GetNumberOfItems()): modality = shNode.GetItemAttribute(item_id, "DICOM.Modality") if modality != "RTSTRUCT": continue rtstruct_node = shNode.GetItemDataNode(item_id) if not rtstruct_node: continue print(f"📎 Najden RTSTRUCT: {rtstruct_node.GetName()}") # Ustvari nov SegmentationNode seg_node = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLSegmentationNode", f"Seg_{rtstruct_node.GetName()}") seg_node.SetReferenceImageGeometryParameterFromVolumeNode(ct_volume_node) seg_node.SetNodeReferenceID("referenceVolume", ct_volume_node.GetID()) seg_node.SetAttribute("DICOM.StructureSetLabel", f"Transformiran_{cbct_date}") seg_node.SetAttribute("DICOM.SeriesDescription", f"Transformiran {studyName.replace('^', ' ')} ({cbct_date})") seg_node.SetAttribute("DICOM.ReferencedSeriesInstanceUID", ct_volume_node.GetAttribute("DICOM.SeriesInstanceUID")) seg_node.SetAttribute("DICOM.ReferencedStudyInstanceUID", ct_volume_node.GetAttribute("DICOM.StudyInstanceUID")) # Najdi child elemente v SubjectHierarchy (strukture) rtstruct_item_id = shNode.GetItemByDataNode(rtstruct_node) structure_ids = vtk.vtkIdList() shNode.GetItemChildren(rtstruct_item_id, structure_ids) segment_count = 0 for i in range(structure_ids.GetNumberOfIds()): structure_item_id = structure_ids.GetId(i) name = shNode.GetItemName(structure_item_id) associated_node = shNode.GetItemDataNode(structure_item_id) if associated_node and associated_node.IsA("vtkMRMLModelNode"): print(f" ➕ Dodajam strukturo: {name}") slicer.modules.segmentations.logic().ImportModelToSegmentationNode(associated_node, seg_node) seg_node.GetSegmentation().GetSegment(seg_node.GetSegmentation().GetNumberOfSegments() - 1).SetName(name) segment_count += 1 # Če ni bilo nič dodano iz SH: poskusi uvoziti modele iz scene if segment_count == 0: print("⚠️ RTSTRUCT nima struktur v SubjectHierarchy – poskus uvoza vseh modelov iz scene.") for model_node in slicer.util.getNodesByClass("vtkMRMLModelNode"): if "RTSTRUCT" in model_node.GetName().upper() or model_node.GetName().startswith("Model"): print(f" ➕ [fallback] Uvoz modela: {model_node.GetName()}") slicer.modules.segmentations.logic().ImportModelToSegmentationNode(model_node, seg_node) seg_node.GetSegmentation().GetSegment(seg_node.GetSegmentation().GetNumberOfSegments() - 1).SetName(model_node.GetName()) print(f"📊 Segmentov v {seg_node.GetName()}: {seg_node.GetSegmentation().GetNumberOfSegments()}") segmentation_nodes.append(seg_node) return segmentation_nodes def apply_cumulative_transform_to_segmentation(segmentation_node, matrix): transform = vtk.vtkTransform() vtk_matrix = vtk.vtkMatrix4x4() for i in range(4): for j in range(4): vtk_matrix.SetElement(i, j, matrix[i, j]) transform.SetMatrix(vtk_matrix) transform_node = slicer.vtkMRMLTransformNode() slicer.mrmlScene.AddNode(transform_node) transform_node.SetAndObserveTransformToParent(transform) segmentation_node.SetAndObserveTransformNodeID(transform_node.GetID()) slicer.vtkSlicerTransformLogic().hardenTransform(segmentation_node) slicer.mrmlScene.RemoveNode(transform_node) def convert_all_models_to_segmentation(reference_volume_name: str, prefix: str = "Imported_"): """ Pretvori vse modele (ModelNode) v sceni v enoten vtkMRMLSegmentationNode. 📥 VHODI: ---------- reference_volume_name : str Ime obstoječega CT volumna (npr. "CT_1"), ki določa geometrijo za segmentacijo. Ta volumen mora biti že naložen v sceni. prefix : str Predpona za ime novega segmentation noda (npr. "Imported_"). Ime novega noda bo nekaj kot: "Imported_Segmentation". 📤 IZHOD: ---------- segmentation_node : vtkMRMLSegmentationNode Nov nod, ki vsebuje en segment za vsak najden model v sceni. Ta segmentacijski nod je pripravljen za transformacijo in DICOM export (vsebuje BinaryLabelmap). """ import slicer # Pridobi referenčni volumen (CT) reference_volume = slicer.util.getNode(reference_volume_name) # Ustvari nov segmentacijski nod segmentation_node = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLSegmentationNode", prefix + "Segmentation") segmentation_node.SetReferenceImageGeometryParameterFromVolumeNode(reference_volume) segmentation_node.SetNodeReferenceID("referenceVolume", reference_volume.GetID()) # Najdi vse modele model_nodes = slicer.util.getNodesByClass("vtkMRMLModelNode") #print(f"📦 Najdenih modelov: {len(model_nodes)}") for model_node in model_nodes: name = model_node.GetName() #print(f"🔍 Model: {name}") skip = ( name.lower().startswith("segmentation") or name.lower().startswith("surface") or name.lower() in ["red volume slice", "green volume slice", "yellow volume slice"] or "rtstruct" in name.lower() ) if skip: continue # Uvozi model kot segment success = slicer.modules.segmentations.logic().ImportModelToSegmentationNode(model_node, segmentation_node) if not success: #print(f"✅ Model '{name}' uvožen kot segment.") print(f"❌ Napaka pri uvozu modela: {name}") # Ustvari BinaryLabelmap reprezentacijo (nujno za DICOM export) created = segmentation_node.GetSegmentation().CreateRepresentation("BinaryLabelmap") if created: print("✅ BinaryLabelmap reprezentacija uspešno ustvarjena.") segmentation_node.GetSegmentation().SetMasterRepresentationName("BinaryLabelmap") #else: #print("❌ Pretvorba v BinaryLabelmap ni uspela.") return segmentation_node def apply_cumulative_transform_to_volume(volume_node, matrix): transform = vtk.vtkTransform() vtk_matrix = vtk.vtkMatrix4x4() for i in range(4): for j in range(4): vtk_matrix.SetElement(i, j, matrix[i, j]) transform.SetMatrix(vtk_matrix) transform_node = slicer.vtkMRMLTransformNode() slicer.mrmlScene.AddNode(transform_node) transform_node.SetAndObserveTransformToParent(transform) volume_node.SetAndObserveTransformNodeID(transform_node.GetID()) slicer.vtkSlicerTransformLogic().hardenTransform(volume_node) slicer.mrmlScene.RemoveNode(transform_node) def dicom_timestamp_from_volume(volume_node): """Vrne npr. '20250112_093045' iz DICOM metapodatkov ali None, če jih ni.""" uids_attr = volume_node.GetAttribute('DICOM.instanceUIDs') if not uids_attr: return None uid0 = uids_attr.split()[0] db = slicer.dicomDatabase fpath = db.fileForInstance(uid0) def tag(tagstr): v = db.fileValue(fpath, tagstr) return v if v not in (None, "", "Unknown") else None # Datumi: Acquisition (0008,0022) → Series (0008,0021) → Study (0008,0020) date = tag("0008,0022") or tag("0008,0021") or tag("0008,0020") # Časi: Acquisition (0008,0032) → Series (0008,0031) → Study (0008,0030) time = tag("0008,0032") or tag("0008,0031") or tag("0008,0030") if not date: return None # Normaliziraj: YYYYMMDD in HHMMSS date = date.replace("-", "").replace(".", "").strip() date = (date + "00000000")[:8] # zaščita, če je prekratko if time: time = time.split(".")[0] # odreži frakcije time = (time + "000000")[:6] return f"{date}_{time}" return date def triangle_similarity(p1, p2): """ Primerja dva trikotnika (vsak definiran s 3 točkami v 3D) na podlagi: - razmerij dolžin - razlik med koti Če je podanih 4 točk, avtomatsko odstrani eno (ujemajočo), ki najmanj vpliva na podobnost. :param p1: seznam (np.array) treh ali štirih točk v 3D (npr. ct_points) :param p2: seznam treh ali štirih točk v 3D (npr. cbct_points) :return: dict z razliko dolžin, razliko kotov, povprečno napako """ def side_lengths(pts): a = np.linalg.norm(pts[1] - pts[0]) b = np.linalg.norm(pts[2] - pts[1]) c = np.linalg.norm(pts[0] - pts[2]) return np.array([a, b, c]) def angles(pts): a, b, c = side_lengths(pts) angle_A = np.arccos((b**2 + c**2 - a**2) / (2 * b * c)) angle_B = np.arccos((a**2 + c**2 - b**2) / (2 * a * c)) angle_C = np.pi - angle_A - angle_B return np.degrees([angle_A, angle_B, angle_C]) # Če imamo 4 točke v vsaki množici, iščemo najboljše 3 ujemajoče if len(p1) == 4 and len(p2) == 4: best_score = np.inf best_idx = None for i in range(4): tri1 = np.delete(p1, i, axis=0) tri2 = np.delete(p2, i, axis=0) l1 = side_lengths(tri1) l2 = side_lengths(tri2) a1 = angles(tri1) a2 = angles(tri2) score = np.mean(np.abs(l1 - l2)) + np.mean(np.abs(a1 - a2)) / 10 if score < best_score: best_score = score best_idx = i print(f"📐 Samodejno odstranjena točka {best_idx} za trikotniško primerjavo.") p1 = np.delete(p1, best_idx, axis=0) p2 = np.delete(p2, best_idx, axis=0) if len(p1) != 3 or len(p2) != 3: raise ValueError("Obe množici točk morata vsebovati 3 (ali 4) točke.") l1 = side_lengths(p1) l2 = side_lengths(p2) a1 = angles(p1) a2 = angles(p2) length_diff = np.abs(l1 - l2) angle_diff = np.abs(a1 - a2) return { "side_length_diff_mm": length_diff, "angle_diff_deg": angle_diff, "mean_length_error_mm": np.mean(length_diff), "mean_angle_error_deg": np.mean(angle_diff), "TSR": 1 / (1 + np.mean(length_diff) + np.mean(angle_diff) / 10) } def _pick_dose_node(): """Izbere RTDOSE; če obstaja 'Accumulated_*', jo vzame prednostno.""" dose_nodes = [n for n in slicer.util.getNodesByClass("vtkMRMLScalarVolumeNode") if "DOSE" in (n.GetName() or "").upper()] if not dose_nodes: return None acc = [n for n in dose_nodes if "ACCUMULATED" in n.GetName().upper()] return acc[0] if acc else dose_nodes[0] def _pick_segmentation_node(): """Vzemi RTSTRUCT segmentacijo (ali prvo smiselno segmentacijo).""" cand = slicer.util.getNodesByClass("vtkMRMLSegmentationNode") if not cand: return None # prednostno RTSTRUCT* for n in cand: if "RTSTRUCT" in n.GetName().upper(): return n return cand[0] def _segment_ids(seg_node, name_filters=None): """Vrne [(segmentId, name), ...] (filtrirano po imenih, če je podano).""" seg = seg_node.GetSegmentation() seg.CreateRepresentation("BinaryLabelmap") all_ids = vtk.vtkStringArray() seg.GetSegmentIDs(all_ids) out = [] all_names = [] for i in range(all_ids.GetNumberOfValues()): sid = all_ids.GetValue(i) name = seg.GetSegment(sid).GetName() or sid all_names.append(name) if (not name_filters) or any(f.lower() in name.lower() for f in name_filters): out.append((sid, name)) #print(f"[DVH/DEBUG] _segment_ids: total={all_ids.GetNumberOfValues()} names={all_names}") #if name_filters: #print(f"[DVH/DEBUG] _segment_ids: filters={name_filters} -> kept={[n for _,n in out]}") return out def compute_dvh_numpy(dose_node, seg_node, name_filters=None, out_csv_path=None): """ Izračuna DVH z notraj-voxel histogramiranjem (numpy). Segmentacijo sproti preslika v geometrijo doze (ExportSegmentsToLabelmapNode), zato ni več težav z 'no overlap' ali geometrijo. Vrne: list [ [Name,Dmin,Dmax,Dmean,D98,D95,D50,D2], ... ] """ logic = slicer.modules.segmentations.logic() results = [] #print("[DVH/DEBUG] compute_dvh_numpy: dose =", dose_node.GetName(), "shape=", # slicer.util.arrayFromVolume(dose_node).shape) #print("[DVH/DEBUG] compute_dvh_numpy: seg =", seg_node.GetName()) try: src_rep = seg_node.GetSegmentation().GetSourceRepresentationName() except Exception: src_rep = "?" #print("[DVH/DEBUG] compute_dvh_numpy: seg =", seg_node.GetName(), "sourceRep=", src_rep) seg_list = _segment_ids(seg_node, name_filters) if not seg_list: print("[DVH] Ni segmentov za obdelavo.") return results dose_arr = slicer.util.arrayFromVolume(dose_node).astype(np.float64) for sid, sname in seg_list: #print(f"[DVH/DEBUG] Processing segment: '{sname}' (id={sid})") # en segment → labelmap v geometriji doze lm = _labelmap_of_segment_on_dose(seg_node, dose_node, sid) if lm is None: print(f"[DVH] Prazna maska za '{sname}'.") continue arrLM = slicer.util.arrayFromVolume(lm) #print("[DVH/DEBUG] labelmap shape:", None if arrLM is None else arrLM.shape, #" nonzero:", None if arrLM is None else int((arrLM>0).sum())) seg_api = seg_node.GetSegmentation() #print("[DVH/DEBUG] Export using sourceRep:", seg_api.GetSourceRepresentationName(), " -> dose:", dose_node.GetName()) mask = arrLM > 0 if arrLM is not None else None slicer.mrmlScene.RemoveNode(lm) if mask is None or not mask.any(): print(f"[DVH] Prazna maska za '{sname}'.") continue vox = dose_arr[mask] #print(f"[DVH/DEBUG] voxels in mask: {vox.size}, Dmean={vox.mean():.2f} Gy") # osnovne metrike (Gy) Dmin = float(np.min(vox)) Dmax = float(np.max(vox)) Dmean = float(np.mean(vox)) D98 = float(np.percentile(vox, 2)) D95 = float(np.percentile(vox, 5)) D50 = float(np.percentile(vox, 50)) D2 = float(np.percentile(vox, 98)) results.append([sname, Dmin, Dmax, Dmean, D98, D95, D50, D2]) if out_csv_path and results: os.makedirs(os.path.dirname(out_csv_path), exist_ok=True) with open(out_csv_path, "w", newline="") as f: w = csv.writer(f) w.writerow(["Segment","Dmin(Gy)","Dmax(Gy)","Dmean(Gy)","D98(Gy)","D95(Gy)","D50(Gy)","D2(Gy)"]) w.writerows(results) return results def _for_uid(volume_node): """Vrne DICOM FrameOfReferenceUID (0020,0052) ali None.""" # 1) neposredno z MRML atributa (če je prisoten) direct = volume_node.GetAttribute('DICOM.FrameOfReferenceUID') if direct: return direct # 2) iz DICOM DB prek prve instance uids_attr = volume_node.GetAttribute('DICOM.instanceUIDs') if not uids_attr: return None db = slicer.dicomDatabase f = db.fileForInstance(uids_attr.split()[0]) return db.fileValue(f, "0020,0052") or None def _pick_matching_dose_node(ct_node): """Vrne RTDOSE z istim FORUID kot CT (če obstaja), sicer None/first fallback.""" ct_for = _for_uid(ct_node) if ct_node else None dose_nodes = [n for n in slicer.util.getNodesByClass("vtkMRMLScalarVolumeNode") if "DOSE" in (n.GetName() or "").upper()] if not dose_nodes: return None if ct_for: for dn in dose_nodes: if _for_uid(dn) == ct_for: return dn # fallback: prva doza (bolje kot nič, a lahko ne prekriva) return dose_nodes[0] def _ensure_overlap_by_resample(dose_node, ct_node): """Če doza in CT nimata istega FORUID, dozo resampla na CT mrežo in vrne novo dozo.""" if not dose_node or not ct_node: return dose_node if _for_uid(dose_node) == _for_uid(ct_node): return dose_node # že matcha import slicer newDose = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLScalarVolumeNode", dose_node.GetName() + "_onCT") params = { "inputVolume": dose_node.GetID(), "referenceVolume": ct_node.GetID(), "outputVolume": newDose.GetID(), "interpolationType": "linear", } slicer.cli.runSync(slicer.modules.resamplescalarvectordwivolume, None, params) print("[DVH] Fallback: resampled dose to CT grid ->", newDose.GetName()) return newDose def _labelmap_of_segment_on_dose(seg_node, dose_node, segment_id): """Vrne TMP labelmap node za en segment v geometriji doze (ali None ob neuspehu).""" logic = slicer.modules.segmentations.logic() # 0) vedno referenciraj geometrijo na DOZO seg_node.SetReferenceImageGeometryParameterFromVolumeNode(dose_node) # 1) poskusi direktni export sid = vtk.vtkStringArray(); sid.InsertNextValue(segment_id) lm = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLLabelMapVolumeNode", f"tmpLM_{segment_id}") ok = logic.ExportSegmentsToLabelmapNode(seg_node, sid, lm, dose_node) if ok: arr = slicer.util.arrayFromVolume(lm) if arr is not None and (arr > 0).any(): return lm slicer.mrmlScene.RemoveNode(lm) # 2) FALLBACK: prisili pot PlanarContour -> Closed surface -> BinaryLabelmap na mreži doze seg = seg_node.GetSegmentation() seg.CreateRepresentation("Closed surface") seg.SetMasterRepresentationName("Closed surface") # kopiraj samo ta segment v novo, "čisto" segmentacijo na mreži doze seg_tmp = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLSegmentationNode", "Seg_onDose_TMP") seg_tmp.SetReferenceImageGeometryParameterFromVolumeNode(dose_node) seg_tmp.GetSegmentation().CopySegmentFromSegmentation(seg, segment_id, segment_id) seg_tmp.GetSegmentation().CreateRepresentation("BinaryLabelmap") lm2 = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLLabelMapVolumeNode", f"tmpLM2_{segment_id}") ok2 = logic.ExportSegmentsToLabelmapNode(seg_tmp, sid, lm2, dose_node) slicer.mrmlScene.RemoveNode(seg_tmp) if not ok2: slicer.mrmlScene.RemoveNode(lm2) return None arr2 = slicer.util.arrayFromVolume(lm2) if arr2 is None or not (arr2 > 0).any(): slicer.mrmlScene.RemoveNode(lm2) return None return lm2 def save_triangle_visualization(ct_pts, cbct_pts, outpath): fig = plt.figure(figsize=(10, 8)) ax = fig.add_subplot(111, projection='3d') ct_closed = np.vstack([ct_pts, ct_pts[0]]) cbct_closed = np.vstack([cbct_pts, cbct_pts[0]]) ax.plot(ct_closed[:, 0], ct_closed[:, 1], ct_closed[:, 2], 'b-', label='CT trikotnik') ax.scatter(ct_pts[:, 0], ct_pts[:, 1], ct_pts[:, 2], c='blue') ax.plot(cbct_closed[:, 0], cbct_closed[:, 1], cbct_closed[:, 2], 'r--', label='CBCT trikotnik') ax.scatter(cbct_pts[:, 0], cbct_pts[:, 1], cbct_pts[:, 2], c='red') ct_centroid = np.mean(ct_pts, axis=0) cbct_centroid = np.mean(cbct_pts, axis=0) ax.scatter(*ct_centroid, c='blue', marker='x', s=60) ax.scatter(*cbct_centroid, c='red', marker='x', s=60) ax.set_title("Trikotnik markerjev: CT (modro) vs CBCT (rdeče)") ax.set_xlabel('X (mm)') ax.set_ylabel('Y (mm)') ax.set_zlabel('Z (mm)') ax.legend() ax.view_init(elev=20, azim=30) plt.tight_layout() try: fig.savefig(outpath) #print(f"🖼 Triangle visualization saved to: {outpath}") except Exception as e: print(f"❌ Failed to save triangle visualization: {e}") finally: plt.close(fig) def resample_scalar_to_reference(input_node, reference_node, interpolator="linear"): """Resampla input_node na geometrijo reference_node in PREPIŠE image data + IJK↔RAS + origin/spacing.""" out = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLScalarVolumeNode", input_node.GetName() + "_tmp") params = { "inputVolume": input_node.GetID(), "referenceVolume": reference_node.GetID(), "outputVolume": out.GetID(), "interpolationType": "linear" if interpolator=="linear" else "nn", } slicer.cli.runSync(slicer.modules.resamplescalarvectordwivolume, None, params) # --- PREPIŠI ImageData --- input_node.SetAndObserveImageData(out.GetImageData()) # --- PREPIŠI GEOMETRIJO (IJK↔RAS) --- ijkToRas = vtk.vtkMatrix4x4(); out.GetIJKToRASMatrix(ijkToRas) rasToIjk = vtk.vtkMatrix4x4(); out.GetRASToIJKMatrix(rasToIjk) input_node.SetIJKToRASMatrix(ijkToRas) input_node.SetRASToIJKMatrix(rasToIjk) # --- PREPIŠI origin/spacing (za vsak slučaj, čeprav IJK↔RAS običajno zadošča) --- input_node.SetOrigin(out.GetOrigin()) input_node.SetSpacing(out.GetSpacing()) # Po želji: prenesi še DisplayNode lastnosti ipd. slicer.mrmlScene.RemoveNode(out) def export_seg_or_rtstruct(ct_volume_Node, seg_node, export_dir): """ Poskusi izvoziti RTSTRUCT (Planar contour). Če ne uspe, izvozi DICOM SEG. Zahteve: - seg_node je že PORAVNAN na CT (apliciran transform) - reference image geometry in UID-ji so nastavljeni """ import slicer from DICOMPlugins import DicomRtImportExportPlugin from DICOMExportSegmentations import DICOMSegmentationExporter # 0) Referenčna geometrija in ujemanje v SubjectHierarchy seg_node.SetReferenceImageGeometryParameterFromVolumeNode(ct_volume_Node) sh = slicer.vtkMRMLSubjectHierarchyNode.GetSubjectHierarchyNode(slicer.mrmlScene) ct_item = sh.GetItemByDataNode(ct_volume_Node) seg_item = sh.GetItemByDataNode(seg_node) # poskrbi, da sta CT in SEG pod ISTIM Study (DICOM exporter to uporablja) study_ct = sh.GetItemParent(ct_item) if study_ct and study_ct != sh.GetItemParent(seg_item): sh.SetItemParent(seg_item, study_ct) # 1) Vnesi ključne DICOM atribute iz CT na SEG (Reference Series in FORUID) ct_series_uid = ct_volume_Node.GetAttribute("DICOM.SeriesInstanceUID") or "" ct_study_uid = ct_volume_Node.GetAttribute("DICOM.StudyInstanceUID") or "" ct_for_uid = ct_volume_Node.GetAttribute("DICOM.FrameOfReferenceUID") or "" if ct_series_uid: sh.SetItemAttribute(seg_item, "DICOM.ReferencedInstanceUIDs", ct_series_uid) if ct_for_uid: seg_node.SetAttribute("DICOM.FrameOfReferenceUID", ct_for_uid) # 2) Poskrbi za reprezentance in nastavi "Source" (nova API namesto "Master") seg = seg_node.GetSegmentation() # Closed surface naj obstaja if not seg.ContainsRepresentation("Closed surface"): seg.CreateRepresentation("Closed surface") # Poskusi ustvariti Planar contour (zahteva veljavno referenčno geometrijo) try: if not seg.ContainsRepresentation("Planar contour"): seg.CreateRepresentation("Planar contour") except Exception as e: print("[RTSTRUCT] CreateRepresentation('Planar contour') failed:", e) # Najprej probamo s Planar contour kot "Source" (nov API) if seg.ContainsRepresentation("Planar contour"): seg.SetSourceRepresentationName("Planar contour") source_repr = "Planar contour" else: seg.SetSourceRepresentationName("Closed surface") source_repr = "Closed surface" print(f"[RTSTRUCT] Source: {seg.GetSourceRepresentationName()}") print(f"[RTSTRUCT] Has Planar: {seg.ContainsRepresentation('Planar contour')}") print(f"[RTSTRUCT] Has Closed: {seg.ContainsRepresentation('Closed surface')}") print("Ref Series UID:", sh.GetItemAttribute(seg_item, "DICOM.ReferencedInstanceUIDs")) print("FOR UID:", seg_node.GetAttribute("DICOM.FrameOfReferenceUID")) # 3) Pripravi exportables in filtriraj RTSTRUCT kandidat(e) plugin = DicomRtImportExportPlugin.DicomRtImportExportPluginClass() def fill_tags(exp): pn = ct_volume_Node.GetAttribute("DICOM.PatientName") pid = ct_volume_Node.GetAttribute("DICOM.PatientID") study_uid = ct_volume_Node.GetAttribute("DICOM.StudyInstanceUID") if pn: exp.setTag("0010,0010", pn) if pid: exp.setTag("0010,0020", pid) if study_uid: exp.setTag("0020,000D", study_uid) exp.directory = export_dir exp_all = [] exp_all += plugin.examineForExport(sh.GetItemByDataNode(ct_volume_Node)) exp_all += plugin.examineForExport(sh.GetItemByDataNode(seg_node)) def pick_rtstruct(exportables): cand = [] for e in exportables: name = getattr(e, "name", "") etype = getattr(e, "exportType", "") sop = getattr(e, "SOPClassUID", "") if etype == "RTSTRUCT" or "RTSTRUCT" in (name or "") or sop == "1.2.840.10008.5.1.4.1.1.481.3": cand.append(e) return cand cand_rt = pick_rtstruct(exp_all) # Če RTSTRUCT ni našel in trenutno source ni Closed surface, poskusi še enkrat if not cand_rt and source_repr != "Closed surface": seg.SetSourceRepresentationName("Closed surface") print("[RTSTRUCT] Retry with source=Closed surface") exp_all = [] exp_all += plugin.examineForExport(sh.GetItemByDataNode(ct_volume_Node)) exp_all += plugin.examineForExport(sh.GetItemByDataNode(seg_node)) cand_rt = pick_rtstruct(exp_all) # 4) Izvoz if cand_rt: for e in cand_rt: fill_tags(e) ok = plugin.export(cand_rt) print(f"✅ RTSTRUCT export status: {ok} → {export_dir}") return "RTSTRUCT" else: # Fallback: DICOM SEG (OpenTPS 2.x ima odlično podporo in DVH deluje) options = DICOMSegmentationExporter.ExportOptions() options.outputDirectory = export_dir options.segmentationNode = seg_node options.referencedVolumeNode = ct_volume_Node DICOMSegmentationExporter.export(options) print(f"✅ DICOM SEG exported → {export_dir} (DVH/CCC v OpenTPS 2.x deluje tudi iz SEG)") return "SEG" # Globalni seznami za končno statistiko prostate_size_est = [] ctcbct_distance = [] table_z_values = {} # Pridobimo SubjectHierarchyNode shNode = slicer.vtkMRMLSubjectHierarchyNode.GetSubjectHierarchyNode(slicer.mrmlScene) studyItems = vtk.vtkIdList() shNode.GetItemChildren(shNode.GetSceneItemID(), studyItems) #slicer.util.delayDisplay(f"[DEBUG] Starting the loops", 1000) for i in range(studyItems.GetNumberOfIds()): study_start_time = time.time() start_io = time.time() studyItem = studyItems.GetId(i) studyName = shNode.GetItemName(studyItem) print(f"\nProcessing study: {studyName}") # **LOKALNI** seznami, resetirajo se pri vsakem study-ju cbct_list = [] ct_list = [] volume_points_dict = {} CT_offset = 0 # Get child items of the study item volumeItems = vtk.vtkIdList() shNode.GetItemChildren(studyItem, volumeItems) # Iteracija čez vse volumne v posameznem studyju for j in range(volumeItems.GetNumberOfIds()): intermediateItem = volumeItems.GetId(j) finalVolumeItems = vtk.vtkIdList() shNode.GetItemChildren(intermediateItem, finalVolumeItems) # Išči globlje! for k in range(finalVolumeItems.GetNumberOfIds()): volumeItem = finalVolumeItems.GetId(k) volumeNode = shNode.GetItemDataNode(volumeItem) try: dicomUIDs = volumeNode.GetAttribute("DICOM.instanceUIDs") except AttributeError: print(f"⚠️ Volume node '{volumeNode}' not found or no attribute 'DICOM.instanceUIDs'. Skip.") dicomUIDs = None continue # Preskoči, če ni veljaven volume if not dicomUIDs: print("❌ This is an NRRD volume!") continue # Preskoči, če ni DICOM volume volumeName = volumeNode.GetName() imageItem = shNode.GetItemByDataNode(volumeNode) modality = shNode.GetItemAttribute(imageItem, "DICOM.Modality") #deluje! #dimensions = volumeNode.GetImageData().GetDimensions() #spacing = volumeNode.GetSpacing() #print(f"Volume {volumeNode.GetName()} - Dimenzije: {dimensions}, Spacing: {spacing}") if modality != "CT": print("Not a CT") continue # Preskoči, če ni CT # Preveri, če volume obstaja v sceni if not slicer.mrmlScene.IsNodePresent(volumeNode): print(f"Volume {volumeName} not present in the scene.") continue # Preverimo proizvajalca (DICOM metapodatki) manufacturer = shNode.GetItemAttribute(imageItem, 'DICOM.Manufacturer') #manufacturer = volumeNode.GetAttribute("DICOM.Manufacturer") #manufacturer = slicer.dicomDatabase.fileValue(uid, "0008,0070") #print(manufacturer) # Določimo, ali gre za CBCT ali CT if "varian" in manufacturer.lower() or "elekta" in manufacturer.lower(): cbct_list.append(volumeName) scan_type = "CBCT" yesCbct = True else: # Siemens ali Philips ct_list.append(volumeName) scan_type = "CT" yesCbct = False if volumeNode and volumeNode.IsA("vtkMRMLScalarVolumeNode"): print(f"✔️ {scan_type} {volumeNode.GetName()} (ID: {volumeItem})") if not volumeNode or not volumeNode.IsA("vtkMRMLScalarVolumeNode"): print("Can't find volumeNode") #continue # Preskoči, če ni veljaven volume # Detekcija točk v volumnu ustvari_marker = not yesCbct # Ustvari markerje pred poravnavo na mizo grouped_points = detect_points_region_growing(volumeName, yesCbct, ustvari_marker, intensity_threshold=3000) if grouped_points is None or len(grouped_points) < 3: print(f"⚠️ Volume {volumeName} doesn't have enough points for registration. Points: {len(grouped_points)}") continue if not yesCbct: # loči koordinate in intenzitete coords_only = [pt for pt, _ in grouped_points] intensities = [intensity for _, intensity in grouped_points] # permutiraj koordinate (npr. zaradi boljšega ujemanja) coords_sorted = match_points(coords_only, coords_only) # ponovno sestavi pare (točka, intenziteta) grouped_points = list(zip(coords_sorted, intensities)) #print(f"Populating volume_points_dict with key ('{scan_type}', '{volumeName}')") volume_points_dict[(scan_type, volumeName)] = grouped_points #print(volume_points_dict) # Check if the key is correctly added # Če imamo oba tipa volumna (CBCT in CT) **znotraj istega studyja** end_io = time.time() if cbct_list and ct_list: fixing = fixing_end = 0 table1_time = table1end_time = 0 table2_time = table2end_time = 0 start_scaling = end_scaling = 0 start_align = end_align = 0 start_rotation = end_rotation = 0 start_translation = end_translation = 0 start_transform = end_transform = 0 study_start_time = study_end_time = 0 ct_volume_name = ct_list[0] # Uporabi prvi CT kot referenco ct_volume_Node = find_volume_node_by_partial_name(ct_volume_name) print(f"\nProcessing CT: {ct_volume_name}") ct_points = [centroid for centroid, _ in volume_points_dict[("CT", ct_volume_name)]] print(f"CT points: {ct_points}") if len(ct_points) < 3: print(f"CT volume {ct_volume_name} doesn't have enough points for registration. Points: {len(ct_points)}") continue else: # if len(ct_points) == 4: # ct_points = remove_lowest_marker(ct_points) #odstrani marker v riti, če obstaja for cbct_volume_name in cbct_list: print(f"\nProcessing CBCT Volume: {cbct_volume_name}") tablefound = False ct_table_found = False cbct_table_found = False CT_offset = None CT_spacing = None yesCbct = True scan_type = "CBCT" #redundant but here we are cbct_volume_node = find_volume_node_by_partial_name(cbct_volume_name) key = (studyName, cbct_volume_name) if key not in cumulative_matrices: cumulative_matrices[key] = np.eye(4) table_shift_matrix = np.eye(4) fixing = time.time() if tablefind: # --- CT table reference --- yesCbct = False table1_time = time.time() makemarkerscheck = True resCT = find_table_top_z(ct_volume_name, writefilecheck, makemarkerscheck, yesCbct) if resCT is not None: ct_edge_mm, _ = resCT ct_table_found = True # zapišemo referenco CT (ne premikamo nič) CT_offset, CT_spacing = align_cbct_to_ct(ct_volume_Node, "CT", ct_edge_mm) else: print("⚠️ CT table top not found – skip table alignment.") ct_table_found = False table1end_time = time.time() # --- CBCT table measure & align --- yesCbct = True table2_time = time.time() makemarkerscheck = False if ct_table_found: resCBCT = find_table_top_z(cbct_volume_name, writefilecheck, makemarkerscheck, yesCbct) #skupnires = resCBCT - resCT if resCBCT is not None else None if resCBCT is not None: cbct_edge_mm, _ = resCBCT # Poravnava CBCT -> CT (harden transform na CBCT) align_cbct_to_ct(cbct_volume_node, "CBCT", cbct_edge_mm, CT_offset, CT_spacing) tablefound = True # ⛔️ NE dodajaj table_shift_matrix v cumulative_matrices (hardeno je že na volumnu) else: print("⚠️ CBCT table top not found – skip table alignment.") tablefound = False table2end_time = time.time() else: tablefound = False resample_scalar_to_reference(cbct_volume_node, ct_volume_Node, interpolator="nn") # Resampla CBCT na CT mrežo cbct_points = [centroid for centroid, _ in volume_points_dict[("CBCT", cbct_volume_name)]] #zastareli podatki cbct_points_array = np.array(cbct_points) # Pretvorba v numpy array if len(cbct_points) != len(ct_points): if(len(cbct_points) + 1 == len(ct_points)): ct_points = remove_lowest_marker(ct_points) #odstrani marker v riti, če obstaja else: print(f"Neujemajoče število točk! CBCT: {len(cbct_points)}, CT: {len(ct_points)}") continue # print_orientation(ct_volume_name) # print_orientation(cbct_volume_name) #for i, (cb, ct) in enumerate(zip(cbct_points, ct_points)): # print(f"Pair {i}: CBCT {cb}, CT {ct}, diff: {np.linalg.norm(cb - ct):.2f}") ustvari_marker = False # Ustvari markerje cbct_points = [centroid for centroid, _ in detect_points_region_growing(cbct_volume_name, yesCbct, ustvari_marker, intensity_threshold=3000)] #cbct_points = detect_points_region_growing(cbct_volume_name, yesCbct, intensity_threshold=3000) #cbct_points = [centroid for centroid, _ in volume_points_dict[("CBCT", cbct_volume_name)]] #zastareli podatki if len(cbct_points) < 3: print(f"CBCT Volume '{cbct_volume_name}' doesn't have enough points for registration. Points: {len(cbct_points)}") continue cbct_spacing = cbct_volume_node.GetSpacing() ct_spacing = ct_volume_Node.GetSpacing() if not np.allclose(cbct_spacing, ct_spacing, atol=1e-6): cbct_points = rescale_points_to_match_spacing(cbct_points, cbct_spacing, ct_spacing) # sicer ne reskaliramo, ker je grid že ujemajoč #Sortiramo točke po X/Y/Z da se izognemo težavam pri poravnavi cbct_points = match_points(cbct_points, ct_points) fixing_end = time.time() #visualize_point_matches_in_slicer(cbct_points, ct_points, studyName) #poveže pare markerjev if writefilecheck: # Shranjevanje razdalj distances_ct_cbct = [] distances_internal = {"A-B": [], "B-C": [], "C-A": []} cbct_volume_node = find_volume_node_by_partial_name(cbct_volume_name) # Sortiramo točke po Z-koordinati (ali X/Y, če raje uporabljaš drugo os) cbct_points_sorted = cbct_points_array[np.argsort(cbct_points_array[:, 2])] # Razdalje med CT in CBCT (SORTIRANE točke!) if cbct_points_sorted.shape[0] != len(ct_points): print(f"⚠️ Število točk CBCT ({cbct_points_sorted.shape[0]}) != CT ({len(ct_points)}), preskakujem izračun.") else: d_ct_cbct = np.linalg.norm(cbct_points_sorted - ct_points, axis=1) distances_ct_cbct.append(d_ct_cbct) # Razdalje med točkami znotraj SORTIRANIH cbct_points d_ab = np.linalg.norm(cbct_points_sorted[0] - cbct_points_sorted[1]) d_bc = np.linalg.norm(cbct_points_sorted[1] - cbct_points_sorted[2]) d_ca = np.linalg.norm(cbct_points_sorted[2] - cbct_points_sorted[0]) # Sortiramo razdalje po velikosti, da so vedno v enakem vrstnem redu sorted_distances = sorted([d_ab, d_bc, d_ca]) distances_internal["A-B"].append(sorted_distances[0]) distances_internal["B-C"].append(sorted_distances[1]) distances_internal["C-A"].append(sorted_distances[2]) # Dodamo ime študije za v statistiko studyName = shNode.GetItemName(studyItem) # **Shrani razdalje v globalne sezname** prostate_size_est.append({"Study": studyName, "Distances": sorted_distances}) ctcbct_distance.append({"Study": studyName, "Distances": list(distances_ct_cbct[-1])}) # Pretvorimo v seznam # Izberi metodo glede na uporabnikov izbor chosen_rotation_matrix = np.eye(3) chosen_translation_vector = np.zeros(3) #print("Markerji pred transformacijo:", cbct_points, ct_points) start_scaling = time.time() scaling_factors = None if applyScaling: scaling_factors = compute_scaling_stddev(cbct_points, ct_points) #print("Scaling factors: ", scaling_factors) cbct_points = compute_scaling(cbct_points, scaling_factors) end_scaling = time.time() start_align = time.time() initial_error = np.mean(np.linalg.norm(np.array(cbct_points) - np.array(ct_points), axis=1)) if initial_error > 30: #print("⚠️ Initial distance too large, applying centroid prealignment.") cbct_points, transvector = prealign_by_centroid(cbct_points, ct_points) else: transvector = np.zeros(3) end_align = time.time() start_rotation = time.time() if applyRotation: if selectedMethod == "Kabsch": chosen_rotation_matrix = compute_Kabsch_rotation(cbct_points, ct_points) elif selectedMethod == "Horn": chosen_rotation_matrix = compute_Horn_rotation(cbct_points, ct_points) elif selectedMethod == "Iterative Closest Point (Horn)": _, chosen_rotation_matrix, _ = icp_algorithm(cbct_points, ct_points) #print("Rotation Matrix:\n", chosen_rotation_matrix) end_rotation = time.time() start_translation = time.time() fine_shift = np.zeros(3) # Inicializiraj fine premike if applyTranslation: chosen_translation_vector, cbct_points_transformed, method_used = choose_best_translation( cbct_points, ct_points, chosen_rotation_matrix) # Sistematična razlika (signed shift) rotated_cbct = np.dot(cbct_points, chosen_rotation_matrix.T) translated_cbct = rotated_cbct + chosen_translation_vector delta_y_list = [ct[1] - cbct[1] for ct, cbct in zip(ct_points, translated_cbct)] mean_delta_y = np.mean(delta_y_list) # Uporabi sistematični shift za dodatno poravnavo v y-osi dy = mean_delta_y if tablefind and tablefound: dy = 0.0 # meja: ali np.clip(dy, -3.0, 3.0) fine_shift = np.array([0.0, dy, 0.0]) cbct_points_transformed += fine_shift end_translation = time.time() start_transform = time.time() # ✅ Kombinirana transformacija total_translation = chosen_translation_vector + fine_shift chosen_translation_vector = total_translation vtk_transform = create_vtk_transform(chosen_rotation_matrix, chosen_translation_vector, tablefound, studyName, cbct_volume_name, scaling_factors) combined_matrix = np.eye(4) # 1. Rotacija if chosen_rotation_matrix is not None: combined_matrix[:3, :3] = chosen_rotation_matrix # 2. Skaliranje if scaling_factors is not None: # Pomnoži rotacijo s skalirnim faktorjem po vsaki osi scaled_rotation = combined_matrix[:3, :3] * scaling_factors # broadcasting vsake vrstice combined_matrix[:3, :3] = scaled_rotation # 3. Translacija if chosen_translation_vector is not None: combined_matrix[:3, 3] = chosen_translation_vector + transvector # združena translacija cumulative_matrices[(studyName, cbct_volume_name)] = np.dot(cumulative_matrices[(studyName, cbct_volume_name)], combined_matrix) # 🔄 Pripni transformacijo imeTransformNoda = cbct_volume_name + " Transform" transform_node = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLTransformNode", imeTransformNoda) transform_node.SetAndObserveTransformToParent(vtk_transform) cbct_volume_node = find_volume_node_by_partial_name(cbct_volume_name) cbct_volume_node.SetAndObserveTransformNodeID(transform_node.GetID()) # 🔨 Uporabi (ali shrani transformacijo kasneje) slicer.vtkSlicerTransformLogic().hardenTransform(cbct_volume_node) slicer.mrmlScene.RemoveNode(transform_node) end_transform = time.time() # 📍 Detekcija markerjev po transformaciji ustvari_marker = False cbct_points = [centroid for centroid, _ in detect_points_region_growing(cbct_volume_name, yesCbct, ustvari_marker, intensity_threshold=3000)] #print("Markerji po transformaciji:\n", cbct_points, ct_points) cbct_points = match_points(cbct_points, ct_points) #popravek v x osi delta_x_list = [ct[0] - cbct[0] for ct, cbct in zip(ct_points, cbct_points)] mean_delta_x = np.mean(delta_x_list) #popravek v y osi delta_y_list = [ct[1] - cbct[1] for ct, cbct in zip(ct_points, cbct_points)] mean_delta_y = np.mean(delta_y_list) #popravek v z osi delta_z_list = [ct[2] - cbct[2] for ct, cbct in zip(ct_points, cbct_points)] mean_delta_z = np.mean(delta_z_list) pre_fine_matrices = {} pre_fine_matrices[(studyName, cbct_volume_name)] = cumulative_matrices[(studyName, cbct_volume_name)].copy() fine_shift = np.zeros(3) if(use_fine_shift): # Uporabi sistematični shift za dodatno poravnavo fine_shift = np.array([mean_delta_x, mean_delta_y, mean_delta_z]) #cbct_points_transformed += fine_shift if fine_shift is not None: shift_matrix = np.eye(4) shift_matrix[:3, 3] = fine_shift cumulative_matrices[(studyName, cbct_volume_name)] = np.dot(cumulative_matrices[(studyName, cbct_volume_name)], shift_matrix) print(f"Fine correction shifts: ΔX={fine_shift[0]:.2f} mm, ΔY={fine_shift[1]:.2f} mm, ΔZ={fine_shift[2]:.2f} mm") chosen_rotation_matrix = np.eye(3) #tokrat brez rotacije #### TEST ROTACIJA ######## # angle_deg = 0 # angle_rad = np.deg2rad(angle_deg) # chosen_rotation_matrix = np.array([ # [np.cos(angle_rad), -np.sin(angle_rad), 0], # [np.sin(angle_rad), np.cos(angle_rad), 0], # [0, 0, 1] # ]) ###KONEC TESTA### vtk_transform = create_vtk_transform(chosen_rotation_matrix, fine_shift, tablefound, studyName, cbct_volume_name) #Tukaj se tudi izpiše transformacijska matrika # 🔄 Pripni transformacijo imeTransformNoda = cbct_volume_name + " Transform" transform_node = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLTransformNode", imeTransformNoda) transform_node.SetAndObserveTransformToParent(vtk_transform) cbct_volume_node = find_volume_node_by_partial_name(cbct_volume_name) cbct_volume_node.SetAndObserveTransformNodeID(transform_node.GetID()) # 🔨 Uporabi (ali shrani transformacijo kasneje) slicer.vtkSlicerTransformLogic().hardenTransform(cbct_volume_node) slicer.mrmlScene.RemoveNode(transform_node) # table_shift_matrix_ct = np.eye(4) # table_shift_matrix_ct[1, 3] = # ali skupni_offset, če je potreben simetričen premik # cumulative_matrices[(studyName, cbct_volume_name)] = np.dot( # table_shift_matrix, # cumulative_matrices[(studyName, cbct_volume_name)] # ) #ustvari_marker = True cbct_points = [centroid for centroid, _ in detect_points_region_growing(cbct_volume_name, yesCbct, applymarkers, intensity_threshold=3000)] #preverjanje determinante rot_part = combined_matrix[:3, :3] det = np.linalg.det(rot_part) if not np.isclose(det, 1.0, atol=0.01): print(f"⚠️ Neortogonalna rotacija! Determinanta: {det}") M_cbct2ct_RAS = cumulative_matrices[(studyName, cbct_volume_name)].copy() # RAS<->LPS pretvornik R_lps = np.diag([-1.0, -1.0, 1.0, 1.0]) # CBCT->CT v LPS M_cbct2ct_LPS = R_lps @ M_cbct2ct_RAS @ R_lps # Kar želiš aplicirati na CT v Ariji/Eclipse (CT->CBCT), rigid v mm, LPS: M_ct2cbct_LPS = np.linalg.inv(M_cbct2ct_LPS) if (resCBCT is not None) and (resCT is not None): deltaY_RAS = float(resCBCT[0]) - float(resCT[0]) # samo mm komponenta (RAS) else: deltaY_RAS = 0.0 # RAS -> LPS: X in Y zamenjata predznak, Z ostane deltaY_LPS = -deltaY_RAS M_ct2cbct_LPS = M_ct2cbct_LPS.copy() M_ct2cbct_LPS[1,3] += deltaY_LPS #shrani transformacijsko matriko v datoteko save_transform_matrix(M_ct2cbct_LPS, studyName, cbct_volume_name) Tx, Ty, Tz = float(M_ct2cbct_LPS[0,3]), float(M_ct2cbct_LPS[1,3]), float(M_ct2cbct_LPS[2,3]) R = M_ct2cbct_LPS[:3,:3] ry = np.degrees(np.arcsin(np.clip(-R[2,0], -1.0, 1.0))) rx = np.degrees(np.arctan2(R[2,1], R[2,2])) rz = np.degrees(np.arctan2(R[1,0], R[0,0])) print(f"Eclipse 6DoF (LPS): Tx={Tx:.2f} mm, Ty={Ty:.2f} mm, Tz={Tz:.2f} mm | Rx={rx:.3f}°, Ry={ry:.3f}°, Rz={rz:.3f}°") # M_ct2cbct_LPS že imaš R = np.diag([-1.0, -1.0, 1.0, 1.0]) M_ct2cbct_RAS = R @ M_ct2cbct_LPS @ R # pretvorba LPS->RAS # 1) ustvari transform node in nastavi 4x4 tNode = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLLinearTransformNode", "CTmove_to_CBCT") m = vtk.vtkMatrix4x4() for r in range(4): for c in range(4): m.SetElement(r, c, float(M_ct2cbct_RAS[r, c])) tNode.SetMatrixTransformToParent(m) # 2) apliciraj na CT in (po želji) na strukture/segmente #ctNode = find_volume_node_by_partial_name(ct_volume_name) # ime CT-ja, ki ga že uporabljaš #ctNode.SetAndObserveTransformNodeID(tNode.GetID()) #slicer.vtkSlicerTransformLogic().hardenTransform(ctNode) # če imaš segmentation/RTSTRUCT kot segmentation node: # segNode = slicer.util.getNode("Segmentation") # prilagodi ime # segNode.SetAndObserveTransformNodeID(tNode.GetID()); slicer.vtkSlicerTransformLogic().hardenTransform(segNode) #print("✅ CT je bil transformiran (harden) v RAS-mm po M_ct2cbct.") #Pridobi podatke o isocentru # iso_coords = extract_isocenter_from_loaded_rtplan() # if iso_coords is not None: # matrix = cumulative_matrices.get((studyName, cbct_volume_name)) # if matrix is not None: # homogeneous_iso = np.append(iso_coords, 1.0) # Doda homogeni člen # transformed_iso = np.dot(matrix, homogeneous_iso) # print(f"🧭 Transformirani izocenter: {transformed_iso[:3]} mm") # print(f"Razlika med izocentroma: ΔX={transformed_iso[0]-iso_coords[0]:.2f} mm, ΔY={transformed_iso[1]-iso_coords[1]:.2f} mm, ΔZ={transformed_iso[2]-iso_coords[2]:.2f} mm") results_base = os.path.join(os.path.dirname(__file__), "Rezultati") study_folder = studyName.replace('^', '_') study_dir = os.path.join(results_base, study_folder) os.makedirs(study_dir, exist_ok=True) cbct_date = dicom_timestamp_from_volume(cbct_volume_node) or datetime.now().strftime("%Y%m%d_%H%M%S") export_dir = os.path.join(study_dir, f"{cbct_date}_DICOM") os.makedirs(export_dir, exist_ok=True) #save_as_dicom = False if save_as_dicom: # Apply transform to CT logging.info(f"[{studyName}] Start applying transform to CT volume.") slicer.util.delayDisplay(f"[DEBUG] Start transform on CT", 1000) apply_cumulative_transform_to_volume(ct_volume_Node, cumulative_matrices[(studyName, cbct_volume_name)]) # Convert RTSTRUCT to SegmentationNode if needed logging.info(f"[{studyName}] Getting segmentation nodes from RTSTRUCT.") #slicer.util.delayDisplay(f"[DEBUG] Start segmentation conversion", 1000) #seg_nodes = convert_rtstruct_to_segmentation_nodes(ct_volume_Node) seg_node = convert_all_models_to_segmentation(ct_volume_name, prefix="Imported_") #new_seg_node = slicer.util.getNode(new_seg_name) apply_cumulative_transform_to_segmentation(seg_node, cumulative_matrices[(studyName, cbct_volume_name)]) plugin = DicomRtImportExportPlugin.DicomRtImportExportPluginClass() shNode = slicer.vtkMRMLSubjectHierarchyNode.GetSubjectHierarchyNode(slicer.mrmlScene) ct_itemID = shNode.GetItemByDataNode(ct_volume_Node) seg_itemID = shNode.GetItemByDataNode(seg_node) cbct_item = shNode.GetItemByDataNode(cbct_volume_node) cbct_date = shNode.GetItemAttribute(cbct_item, "DICOM.SeriesDate") or "unknownDate" #cbct_date_acq = cbct_volume_node.GetAttribute("DICOM.AcquisitionDate") or "unknownDate" #print(f"CBCT date (SeriesDate): {cbct_date}, AcquisitionDate: {cbct_date_acq}") exportables = [] exportables += plugin.examineForExport(ct_itemID) exportables += plugin.examineForExport(seg_itemID) for exp in exportables: # Kopiraj podatke iz volumetričnega noda v exportable if ct_volume_Node.GetAttribute("DICOM.PatientName"): exp.setTag("0010,0010", ct_volume_Node.GetAttribute("DICOM.PatientName")) # PatientName if ct_volume_Node.GetAttribute("DICOM.PatientID"): exp.setTag("0010,0020", ct_volume_Node.GetAttribute("DICOM.PatientID")) # PatientID if ct_volume_Node.GetAttribute("DICOM.StudyInstanceUID"): exp.setTag("0020,000D", ct_volume_Node.GetAttribute("DICOM.StudyInstanceUID")) # StudyInstanceUID exp.directory = export_dir print(f"[{studyName}] ✅ Export successful in [{export_dir}]") plugin.export(exportables) # def export_seg_as_dicom_seg(seg_node, ref_ct_node, out_dir, series_desc="AutoExport SEG"): # # 0) Referenčna geometrija in hierarhija # seg_node.SetReferenceImageGeometryParameterFromVolumeNode(ref_ct_node) # sh = slicer.vtkMRMLSubjectHierarchyNode.GetSubjectHierarchyNode(slicer.mrmlScene) # ct_item = sh.GetItemByDataNode(ref_ct_node) # seg_item = sh.GetItemByDataNode(seg_node) # study_ct = sh.GetItemParent(ct_item) # if study_ct and study_ct != sh.GetItemParent(seg_item): # sh.SetItemParent(seg_item, study_ct) # # 1) Poskrbi za binarni labelmap kot source (DICOM SEG to pričakuje) # seg = seg_node.GetSegmentation() # if not seg.ContainsRepresentation("Binary labelmap"): # seg.CreateRepresentation("Binary labelmap") # seg.SetSourceRepresentationName("Binary labelmap") # # 2) Eksplicitno nastavi DICOM reference (pomoč, če manjka v SH atributih) # ct_series_uid = ref_ct_node.GetAttribute("DICOM.SeriesInstanceUID") or "" # ct_for_uid = ref_ct_node.GetAttribute("DICOM.FrameOfReferenceUID") or "" # if ct_series_uid: # sh.SetItemAttribute(seg_item, "DICOM.ReferencedInstanceUIDs", ct_series_uid) # if ct_for_uid: # seg_node.SetAttribute("DICOM.FrameOfReferenceUID", ct_for_uid) # # 3) IZVOZ (brez dodatnih modulov) # logic = slicer.modules.segmentations.logic() # ok = logic.ExportAllSegmentsToDICOMSegmentation(seg_node, ref_ct_node, out_dir, series_desc) # print(f"✅ DICOM SEG export status: {ok} → {out_dir}") # return ok #export_seg_as_dicom_seg(seg_node, ct_volume_Node, export_dir) # 📏 Izračun napake errors = [np.linalg.norm(cbct - ct) for cbct, ct in zip(cbct_points, ct_points)] mean_error = np.mean(errors) ct_pts = np.array(ct_points) cbct_pts = np.array(cbct_points) ct_centroid = ct_pts.mean(axis=0) cbct_centroid = cbct_pts.mean(axis=0) centroid_delta = cbct_centroid - ct_centroid centroid_norm = float(np.linalg.norm(centroid_delta)) print("Total Individual errors:", errors) print("Average error: {:.2f} mm".format(mean_error)) for i, (cbct, ct) in enumerate(zip(cbct_points, ct_points)): diff = np.array(cbct) - np.array(ct) print(f"Specific marker errors {i+1}: ΔX={diff[0]:.2f} mm, ΔY={diff[1]:.2f} mm, ΔZ={diff[2]:.2f} mm") if len(ct_pts) == 3 and len(cbct_pts) == 3 or len(ct_pts) == 4 and len(cbct_pts) == 4: sim = triangle_similarity(ct_pts, cbct_pts) print("\n Triangle similarity analysis:") print(f"Centroid Δ: ΔX={centroid_delta[0]:.2f} mm, ΔY={centroid_delta[1]:.2f} mm, ΔZ={centroid_delta[2]:.2f} mm | |Δ|={centroid_norm:.2f} mm") print(f"Side length differences (mm): {sim['side_length_diff_mm']}") print(f"Angle differences (deg): {sim['angle_diff_deg']}") print(f"Mean length error: {sim['mean_length_error_mm']:.2f} mm") print(f"Mean angle error: {sim['mean_angle_error_deg']:.2f}°") print(f"Triangle Similarity Ratio (TSR): {sim['TSR']:.3f}") else: print("⚠️ Za primerjavo trikotnikov potrebne točno 3 točke.") if writefilecheck: errorsfile = os.path.join(study_dir, f"{cbct_date}_errors.csv") with open(errorsfile, mode='a', newline='') as file: writer = csv.writer(file) # Glava za študijo writer.writerow(["Study", studyName]) writer.writerow(["Total Individual Errors (mm)"]) for error in errors: writer.writerow(["", f"{error:.2f}"]) writer.writerow(["Average Error (mm)", f"{mean_error:.2f}"]) # Specifične napake markerjev writer.writerow(["Specific Marker Errors"]) writer.writerow(["Marker", "ΔX (mm)", "ΔY (mm)", "ΔZ (mm)"]) for i, (cbct, ct) in enumerate(zip(cbct_points, ct_points)): diff = np.array(cbct) - np.array(ct) writer.writerow([f"Marker {i+1}", f"{diff[0]:.2f}", f"{diff[1]:.2f}", f"{diff[2]:.2f}"]) if len(cbct_points) == 3 and len(ct_points) == 3: writer.writerow([]) writer.writerow(["Triangle Similarity"]) writer.writerow(["Side Length Diff (mm)", *[f"{v:.2f}" for v in sim["side_length_diff_mm"]]]) writer.writerow(["Angle Diff (deg)", *[f"{v:.2f}" for v in sim["angle_diff_deg"]]]) writer.writerow(["Mean Length Error (mm)", f"{sim['mean_length_error_mm']:.2f}"]) writer.writerow(["Mean Angle Error (deg)", f"{sim['mean_angle_error_deg']:.2f}"]) writer.writerow(["Triangle Similarity Ratio (TSR)", f"{sim['TSR']:.3f}"]) writer.writerow(["Centroid (CT) X/Y/Z (mm)", f"{ct_centroid[0]:.2f}", f"{ct_centroid[1]:.2f}", f"{ct_centroid[2]:.2f}"]) writer.writerow(["Centroid (CBCT) X/Y/Z (mm)", f"{cbct_centroid[0]:.2f}", f"{cbct_centroid[1]:.2f}", f"{cbct_centroid[2]:.2f}"]) writer.writerow(["Centroid Δ X/Y/Z (mm) | |Δ|", f"{centroid_delta[0]:.2f}", f"{centroid_delta[1]:.2f}", f"{centroid_delta[2]:.2f}", f"{centroid_norm:.2f}"]) writer.writerow([]) writer.writerow([]) #writer.writerow(["Isocenter at: ", transformed_iso[:3]]) graphs_dir = os.path.join(results_base, study_folder) os.makedirs(graphs_dir, exist_ok=True) pngfile = os.path.join(graphs_dir, f"{cbct_date}_trikotnika.png") save_triangle_visualization(np.array(ct_points), np.array(cbct_points), pngfile) else: print(f"Study {studyItem} doesn't have any appropriate CT or CBCT volumes.") continue study_end_time = time.time() timing_data = { "IO": end_io - start_io, "Fixing": fixing_end - fixing, "Table": ((table1end_time - table1_time)+(table2end_time - table2_time)) if tablefind else 0, "Scaling": end_scaling - start_scaling, "CentroidAlign": end_align - start_align, "Rotation": end_rotation - start_rotation, "Translation": end_translation - start_translation, "Transform": end_transform - start_transform, "Total": study_end_time - study_start_time} update_timing_csv(timing_data, studyName) print(f"Timing data for pacient: {timing_data}") # Izpis globalne statistike start_save = time.time() if writefilecheck: #print("Distances between CT & CBCT markers: ", ctcbct_distance) #print("Distances between pairs of markers for each volume: ", prostate_size_est) # Define file paths prostate_size_file = os.path.join(os.path.dirname(__file__), "prostate_size.csv") ctcbct_distance_file = os.path.join(os.path.dirname(__file__), "ct_cbct_distance.csv") # Write prostate size data with open(prostate_size_file, mode='w', newline='') as file: writer = csv.writer(file) writer.writerow(["Prostate Size"]) for size in prostate_size_est: writer.writerow([size]) #print("Prostate size file written at ", prostate_size_file) # Write CT-CBCT distance data with open(ctcbct_distance_file, mode='w', newline='') as file: writer = csv.writer(file) writer.writerow(["CT-CBCT Distance"]) for distance in ctcbct_distance: writer.writerow([distance]) #print("CT-CBCT distance file written at ", ctcbct_distance_file) end_save = time.time() print(f"Saving time: {end_save - start_save:.2f} seconds") end_time = time.time() # Calculate and print elapsed time elapsed_time = end_time - start_time # Convert to minutes and seconds minutes = int(elapsed_time // 60) seconds = elapsed_time % 60 print(f"Execution time: {minutes} minutes and {seconds:.6f} seconds")