|
|
@@ -9,6 +9,8 @@ from collections import deque
|
|
|
import vtk
|
|
|
from slicer.ScriptedLoadableModule import *
|
|
|
import qt
|
|
|
+import matplotlib.pyplot as plt
|
|
|
+import csv
|
|
|
|
|
|
#exec(open("C:/Users/lkomar/Documents/Prostata/FirstTryRegister.py").read())
|
|
|
|
|
|
@@ -35,7 +37,20 @@ class SeekTransformModuleWidget(ScriptedLoadableModuleWidget):
|
|
|
self.rotationMethodComboBox = qt.QComboBox()
|
|
|
self.rotationMethodComboBox.addItems(["Kabsch", "Horn", "Iterative Closest Point (Kabsch)"])
|
|
|
self.layout.addWidget(self.rotationMethodComboBox)
|
|
|
-
|
|
|
+
|
|
|
+ # Checkboxi za transformacije
|
|
|
+ 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.scalingCheckBox = qt.QCheckBox("Scaling")
|
|
|
+ self.scalingCheckBox.setChecked(True)
|
|
|
+ self.layout.addWidget(self.scalingCheckBox)
|
|
|
+
|
|
|
# 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."
|
|
|
@@ -49,15 +64,26 @@ class SeekTransformModuleWidget(ScriptedLoadableModuleWidget):
|
|
|
|
|
|
def onApplyButton(self):
|
|
|
logic = MyTransformModuleLogic()
|
|
|
- selectedMethod = self.rotationMethodComboBox.currentText #izberi metodo izračuna rotacije
|
|
|
- logic.run(selectedMethod)
|
|
|
+ selectedMethod = self.rotationMethodComboBox.currentText # izberi metodo izračuna rotacije
|
|
|
+
|
|
|
+ # Preberi stanje checkboxov
|
|
|
+ applyRotation = self.rotationCheckBox.isChecked()
|
|
|
+ applyTranslation = self.translationCheckBox.isChecked()
|
|
|
+ applyScaling = self.scalingCheckBox.isChecked()
|
|
|
+
|
|
|
+ # Pokliči logiko z izbranimi nastavitvami
|
|
|
+ logic.run(selectedMethod, applyRotation, applyTranslation, applyScaling)
|
|
|
+
|
|
|
|
|
|
class MyTransformModuleLogic(ScriptedLoadableModuleLogic):
|
|
|
"""
|
|
|
Core logic of the module.
|
|
|
"""
|
|
|
- def run(self, selectedMethod):
|
|
|
+
|
|
|
+
|
|
|
+ def run(self, selectedMethod, applyRotation, applyTranslation, applyScaling):
|
|
|
print("Calculating...")
|
|
|
+
|
|
|
def group_points(points, threshold):
|
|
|
# Function to group points that are close to each other
|
|
|
grouped_points = []
|
|
|
@@ -105,74 +131,50 @@ class MyTransformModuleLogic(ScriptedLoadableModuleLogic):
|
|
|
|
|
|
return region
|
|
|
|
|
|
- def detect_points_region_growing(volume_name, yesCbct, intensity_threshold=3000, x_min=90, x_max=380, y_min=190, y_max=380, z_min=80, z_max=140, max_distance=9, centroid_merge_threshold=5):
|
|
|
- volume_node = slicer.util.getNode(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)
|
|
|
+
|
|
|
+ def compute_optimal_scaling_per_axis(moving_points, fixed_points):
|
|
|
+ """Computes optimal scaling factors for each axis (X, Y, Z) to align moving points (CBCT) to fixed points (CT).
|
|
|
|
|
|
- dimensions = image_data.GetDimensions()
|
|
|
- #detected_regions = []
|
|
|
+ Args:
|
|
|
+ moving_points (list of lists): List of (x, y, z) moving points (CBCT).
|
|
|
+ fixed_points (list of lists): List of (x, y, z) fixed points (CT).
|
|
|
|
|
|
- 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)
|
|
|
+ Returns:
|
|
|
+ tuple: Scaling factors (sx, sy, sz).
|
|
|
+ """
|
|
|
+ moving_points_np = np.array(moving_points)
|
|
|
+ fixed_points_np = np.array(fixed_points)
|
|
|
|
|
|
- visited = set()
|
|
|
+ # Compute centroids
|
|
|
+ centroid_moving = np.mean(moving_points_np, axis=0)
|
|
|
+ centroid_fixed = np.mean(fixed_points_np, axis=0)
|
|
|
|
|
|
- def grow_region(x, y, z):
|
|
|
- if (x, y, z) in visited:
|
|
|
- return None
|
|
|
+ # Compute absolute distances of each point from its centroid along each axis
|
|
|
+ distances_moving = np.abs(moving_points_np - centroid_moving)
|
|
|
+ distances_fixed = np.abs(fixed_points_np - centroid_fixed)
|
|
|
|
|
|
- voxel_value = image_data.GetScalarComponentAsDouble(x, y, z, 0)
|
|
|
- if voxel_value < intensity_threshold:
|
|
|
- return None
|
|
|
+ # Compute scaling factors as the ratio of mean absolute distances per axis
|
|
|
+ scale_factors = np.mean(distances_fixed, axis=0) / np.mean(distances_moving, axis=0)
|
|
|
|
|
|
- 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
|
|
|
+ return tuple(scale_factors)
|
|
|
|
|
|
- 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)
|
|
|
+ def compute_scaling(cbct_points, scaling_factors):
|
|
|
+ """Applies non-uniform scaling to CBCT points.
|
|
|
|
|
|
- # 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))
|
|
|
+ Args:
|
|
|
+ cbct_points (list of lists): List of (x, y, z) points.
|
|
|
+ scaling_factors (tuple): Scaling factors (sx, sy, sz) for each axis.
|
|
|
|
|
|
- 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))
|
|
|
-
|
|
|
- markups_node = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLMarkupsFiducialNode", f"Markers_{volume_name}")
|
|
|
- for centroid, intensity in unique_centroids:
|
|
|
- markups_node.AddControlPoint(*centroid)
|
|
|
- #print(f"Detected Centroid (RAS): {centroid}, Max Intensity: {intensity}")
|
|
|
+ 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
|
|
|
|
|
|
- return unique_centroids
|
|
|
+ cbct_points_np = np.array(cbct_points) # Convert to numpy array
|
|
|
+ scaled_points = cbct_points_np @ scaling_matrix.T # Apply scaling
|
|
|
+
|
|
|
+ return scaled_points.tolist() # Convert back to list
|
|
|
|
|
|
def compute_Kabsch_rotation(moving_points, fixed_points):
|
|
|
"""
|
|
|
@@ -348,7 +350,7 @@ class MyTransformModuleLogic(ScriptedLoadableModuleLogic):
|
|
|
|
|
|
return translation
|
|
|
|
|
|
- def create_vtk_transform(rotation_matrix, translation_vector):
|
|
|
+ def create_vtk_transform(rotation_matrix, translation_vector):
|
|
|
"""
|
|
|
Creates a vtkTransform from a rotation matrix and a translation vector.
|
|
|
"""
|
|
|
@@ -369,130 +371,299 @@ class MyTransformModuleLogic(ScriptedLoadableModuleLogic):
|
|
|
transform = vtk.vtkTransform()
|
|
|
transform.SetMatrix(vtk_matrix)
|
|
|
return transform
|
|
|
+
|
|
|
+
|
|
|
+ def detect_points_region_growing(volume_name, yesCbct, intensity_threshold=3000, x_min=90, x_max=380, y_min=190, y_max=380, z_min=80, z_max=140, max_distance=9, centroid_merge_threshold=5):
|
|
|
+ volume_node = slicer.util.getNode(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 = []
|
|
|
|
|
|
- # Initialize lists and dictionary
|
|
|
- cbct_list = []
|
|
|
- ct_list = []
|
|
|
- volume_points_dict = {}
|
|
|
-
|
|
|
- # Process loaded volumes
|
|
|
- for volumeNode in slicer.util.getNodesByClass("vtkMRMLScalarVolumeNode"):
|
|
|
- volumeName = volumeNode.GetName()
|
|
|
- #print(volumeName)
|
|
|
- shNode = slicer.vtkMRMLSubjectHierarchyNode.GetSubjectHierarchyNode(slicer.mrmlScene)
|
|
|
- imageItem = shNode.GetItemByDataNode(volumeNode)
|
|
|
-
|
|
|
- # Pridobi vse atribute za ta element
|
|
|
- #attributeNames = shNode.GetItemAttributeNames(imageItem)
|
|
|
-
|
|
|
- #modality = shNode.GetItemAttribute(imageItem, 'DICOM.Modality')
|
|
|
- #manufacturer = shNode.GetItemAttribute(imageItem, 'DICOM.Manufacturer')
|
|
|
- #print(modality)
|
|
|
-
|
|
|
- # Check if the volume is loaded into the scene
|
|
|
- if not slicer.mrmlScene.IsNodePresent(volumeNode):
|
|
|
- print(f"Volume {volumeName} not present in the scene.")
|
|
|
- continue
|
|
|
-
|
|
|
- manufacturer = shNode.GetItemAttribute(imageItem, 'DICOM.Manufacturer')
|
|
|
- #print(manufacturer.lower())
|
|
|
- # Determine scan type
|
|
|
- if "varian" in manufacturer.lower() or "elekta" in manufacturer.lower():
|
|
|
- cbct_list.append(volumeName)
|
|
|
- scan_type = "CBCT"
|
|
|
- yesCbct = True
|
|
|
- else: #Siemens or phillips imamo CTje
|
|
|
- ct_list.append(volumeName)
|
|
|
- scan_type = "CT"
|
|
|
- yesCbct = False
|
|
|
-
|
|
|
- # Detect points using region growing
|
|
|
- grouped_points = detect_points_region_growing(volumeName, yesCbct, intensity_threshold=3000)
|
|
|
- volume_points_dict[(scan_type, volumeName)] = grouped_points
|
|
|
+ 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)
|
|
|
|
|
|
- # Print the results
|
|
|
- # print(f"\nCBCT Volumes: {cbct_list}")
|
|
|
- # print(f"CT Volumes: {ct_list}")
|
|
|
- # print("\nDetected Points by Volume:")
|
|
|
- # for (scan_type, vol_name), points in volume_points_dict.items():
|
|
|
- # print(f"{scan_type} Volume '{vol_name}': {len(points)} points detected.")
|
|
|
+ visited = set()
|
|
|
|
|
|
+ def grow_region(x, y, z):
|
|
|
+ if (x, y, z) in visited:
|
|
|
+ return None
|
|
|
|
|
|
- if cbct_list and ct_list:
|
|
|
- # Izberi prvi CT volumen kot referenco
|
|
|
- ct_volume_name = ct_list[0]
|
|
|
- ct_points = [centroid for centroid, _ in volume_points_dict[("CT", ct_volume_name)]]
|
|
|
+ voxel_value = image_data.GetScalarComponentAsDouble(x, y, z, 0)
|
|
|
+ if voxel_value < intensity_threshold:
|
|
|
+ return None
|
|
|
|
|
|
- if len(ct_points) < 3:
|
|
|
- print("CT volumen nima dovolj točk za registracijo.")
|
|
|
- else:
|
|
|
- #print("CT points: ", np.array(ct_points))
|
|
|
-
|
|
|
- for cbct_volume_name in cbct_list:
|
|
|
- # Izvleci točke za trenutni CBCT volumen
|
|
|
- cbct_points = [centroid for centroid, _ in volume_points_dict[("CBCT", cbct_volume_name)]]
|
|
|
+ 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
|
|
|
|
|
|
- print(f"\nProcessing CBCT Volume: {cbct_volume_name}")
|
|
|
- if len(cbct_points) < 3:
|
|
|
- print(f"CBCT Volume '{cbct_volume_name}' nima dovolj točk za registracijo.")
|
|
|
- continue
|
|
|
+ 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)
|
|
|
|
|
|
- #print("CBCT points: ", np.array(cbct_points))
|
|
|
+ # 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))
|
|
|
|
|
|
- # Display the results for the current CBCT volume
|
|
|
- # print("\nSVD Method:")
|
|
|
- # print("Rotation Matrix:\n", svd_rotation_matrix)
|
|
|
- # print("Translation Vector:\n", svd_translation_vector)
|
|
|
+ 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))
|
|
|
+
|
|
|
+ markups_node = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLMarkupsFiducialNode", f"Markers_{volume_name}")
|
|
|
+ for centroid, intensity in unique_centroids:
|
|
|
+ markups_node.AddControlPoint(*centroid)
|
|
|
+ #print(f"Detected Centroid (RAS): {centroid}, Max Intensity: {intensity}")
|
|
|
|
|
|
- # print("\nHorn Method:")
|
|
|
- # print("Rotation Matrix:\n", horn_rotation_matrix)
|
|
|
- # print("Translation Vector:\n", horn_translation_vector)
|
|
|
+ return unique_centroids
|
|
|
+
|
|
|
+
|
|
|
+
|
|
|
+ # Globalni seznami za končno statistiko
|
|
|
+ prostate_size_est = []
|
|
|
+ ctcbct_distance = []
|
|
|
|
|
|
- # print("\nQuaternion Method:")
|
|
|
- # print("Rotation Matrix:\n", quaternion_rotation_matrix)
|
|
|
- # print("Translation Vector:\n", quaternion_translation_vector)
|
|
|
+ # Pridobimo SubjectHierarchyNode
|
|
|
+ shNode = slicer.vtkMRMLSubjectHierarchyNode.GetSubjectHierarchyNode(slicer.mrmlScene)
|
|
|
+
|
|
|
+ studyItems = vtk.vtkIdList()
|
|
|
+ shNode.GetItemChildren(shNode.GetSceneItemID(), studyItems)
|
|
|
|
|
|
+ for i in range(studyItems.GetNumberOfIds()):
|
|
|
+ studyItem = studyItems.GetId(i)
|
|
|
+
|
|
|
+ # Dodamo ime študije v statistiko
|
|
|
+ studyName = shNode.GetItemName(studyItem)
|
|
|
+ prostate_size_est.append({"Study": studyName})
|
|
|
+ ctcbct_distance.append({"Study": studyName})
|
|
|
+
|
|
|
+ # **LOKALNI** seznami, resetirajo se pri vsakem study-ju
|
|
|
+ cbct_list = []
|
|
|
+ ct_list = []
|
|
|
+ volume_points_dict = {}
|
|
|
+
|
|
|
+ # 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)
|
|
|
+
|
|
|
+ # Preveri, ali je to dejanska skupina volumnov (npr. "No study description")
|
|
|
+ intermediateName = shNode.GetItemName(intermediateItem)
|
|
|
+ #print(f"Checking intermediate item: {intermediateName}")
|
|
|
+
|
|
|
+ finalVolumeItems = vtk.vtkIdList()
|
|
|
+ shNode.GetItemChildren(intermediateItem, finalVolumeItems) # Išči globlje!
|
|
|
+
|
|
|
+ for k in range(finalVolumeItems.GetNumberOfIds()):
|
|
|
+ volumeItem = finalVolumeItems.GetId(k)
|
|
|
+ volumeNode = shNode.GetItemDataNode(volumeItem)
|
|
|
|
|
|
+ dicomUIDs = volumeNode.GetAttribute("DICOM.instanceUIDs")
|
|
|
+ if not dicomUIDs:
|
|
|
+ print("❌ To je NRRD volume!")
|
|
|
+ continue # Preskoči, če ni DICOM volume
|
|
|
+
|
|
|
+ if volumeNode and volumeNode.IsA("vtkMRMLScalarVolumeNode"):
|
|
|
+ print(f"✔️ Najden volume: {volumeNode.GetName()} (ID: {volumeItem})")
|
|
|
+
|
|
|
+ if not volumeNode or not volumeNode.IsA("vtkMRMLScalarVolumeNode"):
|
|
|
+ print("ne najdem volumeNode")
|
|
|
+ continue # Preskoči, če ni veljaven volume
|
|
|
+
|
|
|
+ # Preveri, če volume ima StorageNode (drugače `.GetFileName()` vrže napako)
|
|
|
+ storageNode = volumeNode.GetStorageNode()
|
|
|
|
|
|
- # Izberi metodo glede na uporabnikov izbor
|
|
|
- if selectedMethod == "Kabsch":
|
|
|
- chosen_rotation_matrix = compute_Kabsch_rotation(cbct_points, ct_points)
|
|
|
- chosen_translation_vector = compute_translation(cbct_points, ct_points, chosen_rotation_matrix)
|
|
|
- print("\nKabsch Method:")
|
|
|
- print("Rotation Matrix:\n", chosen_rotation_matrix)
|
|
|
- print("Translation Vector:\n", chosen_translation_vector)
|
|
|
- elif selectedMethod == "Horn":
|
|
|
- chosen_rotation_matrix = compute_Horn_rotation(cbct_points, ct_points)
|
|
|
- chosen_translation_vector = compute_translation(cbct_points, ct_points, chosen_rotation_matrix)
|
|
|
- print("\nHorn Method:")
|
|
|
- print("Rotation Matrix:\n", chosen_rotation_matrix)
|
|
|
- print("Translation Vector:\n", chosen_translation_vector)
|
|
|
- elif selectedMethod == "Iterative Closest Point (Kabsch)":
|
|
|
- new_points, chosen_rotation_matrix, chosen_translation_vector = icp_algorithm(cbct_points, ct_points)
|
|
|
- #chosen_translation_vector = compute_translation(cbct_points, ct_points, chosen_rotation_matrix)
|
|
|
- print("\Iterative Closest Point Method:")
|
|
|
- print("Rotation Matrix:\n", chosen_rotation_matrix)
|
|
|
- print("Translation Vector:\n", chosen_translation_vector)
|
|
|
-
|
|
|
- imeTransformNoda = cbct_volume_name + " Transform"
|
|
|
- # Ustvari vtkTransformNode in ga poveži z CBCT volumenom
|
|
|
- transform_node = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLTransformNode", imeTransformNoda)
|
|
|
- # Kliči funkcijo, ki uporabi matriki
|
|
|
- vtk_transform = create_vtk_transform(chosen_rotation_matrix, chosen_translation_vector)
|
|
|
- # Dodaj transform v node
|
|
|
- transform_node.SetAndObserveTransformToParent(vtk_transform)
|
|
|
-
|
|
|
- # Pridobi CBCT volumen in aplikacijo transformacije
|
|
|
- cbct_volume_node = slicer.util.getNode(cbct_volume_name)
|
|
|
- cbct_volume_node.SetAndObserveTransformNodeID(transform_node.GetID()) # Pripni transform node na volumen
|
|
|
-
|
|
|
- # Uporabi transformacijo na volumnu (fizična aplikacija)
|
|
|
- slicer.vtkSlicerTransformLogic().hardenTransform(cbct_volume_node) # Uporabi transform na volumen
|
|
|
- print("Transform uspešen na", cbct_volume_name)
|
|
|
|
|
|
+ if not storageNode:
|
|
|
+ print("ne najdem storageNode")
|
|
|
+ continue # Preskoči, če volume nima shranjenih DICOM podatkov
|
|
|
+ volumeName = volumeNode.GetName()
|
|
|
+ #print(volumeName)
|
|
|
+ imageItem = shNode.GetItemByDataNode(volumeNode)
|
|
|
+ #print(imageItem)
|
|
|
+ #dicomUIDsList = volumeNode.GetAttribute("DICOM.instanceUIDs").split()
|
|
|
+ # Preverimo modaliteto volumna (DICOM metapodatki)
|
|
|
+ #modality = slicer.dicomDatabase.fileValue(storageNode.GetFileName(), "0008,0060") #prazen
|
|
|
+ #modality = volumeNode.GetAttribute("DICOM.Modality") #None
|
|
|
+ #modality = slicer.dicomDatabase.fileValue(uid, "0008,0060") # Modality #prazen
|
|
|
+ #modality = slicer.dicomDatabase.fileValue(dicomUIDsList[0], "0008,0060") #prazen
|
|
|
+ modality = shNode.GetItemAttribute(imageItem, "DICOM.Modality") #deluje!
|
|
|
+ #print(modality)
|
|
|
+
|
|
|
+ dimensions = volumeNode.GetImageData().GetDimensions()
|
|
|
+ spacing = volumeNode.GetSpacing()
|
|
|
+ #print(f"Volume {volumeNode.GetName()} - Dimenzije: {dimensions}, Spacing: {spacing}")
|
|
|
+
|
|
|
+ if modality != "CT":
|
|
|
+ print("Ni CT slika")
|
|
|
+ continue # Preskoči, če ni CT
|
|
|
+
|
|
|
|
|
|
- #transformed_cbct_image = create_vtk_transform(cbct_image_sitk, chosen_rotation_matrix, chosen_translation_vector)
|
|
|
+ # 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
|
|
|
+ print("CBCT")
|
|
|
+ else: # Siemens ali Philips
|
|
|
+ ct_list.append(volumeName)
|
|
|
+ scan_type = "CT"
|
|
|
+ yesCbct = False
|
|
|
+ print("CT")
|
|
|
+
|
|
|
+ # Detekcija točk v volumnu
|
|
|
+ grouped_points = detect_points_region_growing(volumeName, yesCbct, 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**
|
|
|
+ if cbct_list and ct_list:
|
|
|
+ ct_volume_name = ct_list[0] # Uporabi prvi CT kot referenco
|
|
|
+ ct_points = [centroid for centroid, _ in volume_points_dict[("CT", ct_volume_name)]]
|
|
|
+
|
|
|
+ if len(ct_points) < 3:
|
|
|
+ print(f"CT volumen {ct_volume_name} nima dovolj točk za registracijo.")
|
|
|
+ else:
|
|
|
+ for cbct_volume_name in cbct_list:
|
|
|
+ cbct_points = [centroid for centroid, _ in volume_points_dict[("CBCT", cbct_volume_name)]]
|
|
|
+
|
|
|
+ print(f"\nProcessing CBCT Volume: {cbct_volume_name}")
|
|
|
+ if len(cbct_points) < 3:
|
|
|
+ print(f"CBCT Volume '{cbct_volume_name}' nima dovolj točk za registracijo.")
|
|
|
+ continue
|
|
|
+
|
|
|
+ # Shranjevanje razdalj
|
|
|
+ distances_ct_cbct = []
|
|
|
+ distances_internal = {"A-B": [], "B-C": [], "C-A": []}
|
|
|
+
|
|
|
+ cbct_points_array = np.array(cbct_points) # Pretvorba v numpy array
|
|
|
+
|
|
|
+ ct_volume_node = slicer.util.getNode(ct_volume_name)
|
|
|
+ cbct_volume_node = slicer.util.getNode(cbct_volume_name)
|
|
|
+ ct_spacing = ct_volume_node.GetSpacing() # (x_spacing, y_spacing, z_spacing)
|
|
|
+ cbct_spacing = cbct_volume_node.GetSpacing() # (x_spacing, y_spacing, z_spacing)
|
|
|
+
|
|
|
+ ct_scale_factor = np.array(ct_spacing) # Spacing za CT (x, y, z)
|
|
|
+ cbct_scale_factor = np.array(cbct_spacing) # Spacing za CBCT (x, y, z)
|
|
|
+ print(ct_scale_factor, cbct_scale_factor)
|
|
|
+
|
|
|
+
|
|
|
+ # 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])
|
|
|
+
|
|
|
+ # **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)
|
|
|
+
|
|
|
+ if applyScaling:
|
|
|
+ scaling_factors = compute_optimal_scaling_per_axis(cbct_points, ct_points)
|
|
|
+ print("Skalirni faktorji: ", scaling_factors)
|
|
|
+ cbct_points = compute_scaling(cbct_points, scaling_factors)
|
|
|
+
|
|
|
+ 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)
|
|
|
+
|
|
|
+ if applyTranslation:
|
|
|
+ chosen_translation_vector = compute_translation(cbct_points, ct_points, chosen_rotation_matrix)
|
|
|
+ print("Translation Vector:\n", chosen_translation_vector)
|
|
|
+
|
|
|
+ # Ustvari vtkTransformNode in ga poveži z CBCT volumenom
|
|
|
+ imeTransformNoda = cbct_volume_name + " Transform"
|
|
|
+ transform_node = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLTransformNode", imeTransformNoda)
|
|
|
+
|
|
|
+ # Kreiraj transformacijo in jo uporabi
|
|
|
+ vtk_transform = create_vtk_transform(chosen_rotation_matrix, chosen_translation_vector)
|
|
|
+ transform_node.SetAndObserveTransformToParent(vtk_transform)
|
|
|
+
|
|
|
+ # Pridobi CBCT volumen in aplikacijo transformacije
|
|
|
+ cbct_volume_node = slicer.util.getNode(cbct_volume_name)
|
|
|
+ cbct_volume_node.SetAndObserveTransformNodeID(transform_node.GetID())
|
|
|
+
|
|
|
+ # Uporabi transformacijo na volumnu (fizična aplikacija)
|
|
|
+ slicer.vtkSlicerTransformLogic().hardenTransform(cbct_volume_node)
|
|
|
+ print("Transform uspešen na", cbct_volume_name)
|
|
|
+
|
|
|
+ else:
|
|
|
+ print(f"Study {studyItem} nima ustreznih CBCT in CT volumnov.")
|
|
|
+
|
|
|
+ # Izpis globalne statistike
|
|
|
+ print("Razdalje med CT in CBCTji: ", ctcbct_distance)
|
|
|
+ print("Razdalje med markerji: ", prostate_size_est)
|
|
|
+
|
|
|
+ # Define the file path for the CSV file
|
|
|
+ file_path = os.path.join(os.path.dirname(__file__), "study_data.csv")
|
|
|
+
|
|
|
+ # Write lists to the CSV file
|
|
|
+ with open(file_path, mode='w', newline='') as file: #w za write, a za append
|
|
|
+ writer = csv.writer(file)
|
|
|
+ # Write headers
|
|
|
+ writer.writerow(["Prostate Size", "CT-CBCT Distance"])
|
|
|
+ # Write data rows
|
|
|
+ for i in range(len(prostate_size_est)):
|
|
|
+ writer.writerow([prostate_size_est[i], ctcbct_distance[i]])
|
|
|
|
|
|
- else:
|
|
|
- print("CBCT ali CT volumen ni bil najden.")
|
Luka