AvtomatskaPoravnavaMarkerjev/FirstTryLinear.txt

import os
import numpy as np
import scipy
from scipy.spatial.distance import cdist
from scipy.spatial.transform import Rotation as R
import slicer
from DICOMLib import DICOMUtils
from collections import deque

# Define a threshold for grouping nearby points (in voxel space)
#distance_threshold = 4  # This can be adjusted based on your dataset

# Function to group points that are close to each other
def group_points(points, threshold):
    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 detect_points_region_growing(volume_name, intensity_threshold=3000, x_min=90, x_max=380, y_min=190, y_max=380, z_min=80, z_max=120, 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 = []

    # Check if it's CT or CBCT
    is_cbct = "cbct" in volume_name.lower()

    if is_cbct:
        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))
            
    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}")

    return unique_centroids

def compute_rigid_transform(moving_points, fixed_points):
    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)
    R_optimal = np.dot(Vt.T, U.T)

    # Correct improper rotation (reflection)
    if np.linalg.det(R_optimal) < 0:
        Vt[-1, :] *= -1
        R_optimal = np.dot(Vt.T, U.T)

    # Compute translation
    translation = centroid_fixed - np.dot(centroid_moving, R_optimal)

    return R_optimal, translation

def apply_transform(points, rotation_matrix, translation_vector):
    points = np.array(points)
    transformed_points = np.dot(points, rotation_matrix.T) + translation_vector
    return transformed_points




# Initialize lists and dictionary
cbct_list = []
ct_list = []
volume_points_dict = {}

# Process loaded volumes
for volumeNode in slicer.util.getNodesByClass("vtkMRMLScalarVolumeNode"):
    volumeName = volumeNode.GetName()
    shNode = slicer.vtkMRMLSubjectHierarchyNode.GetSubjectHierarchyNode(slicer.mrmlScene)
    imageItem = shNode.GetItemByDataNode(volumeNode)
    
    modality = shNode.GetItemAttribute(imageItem, 'DICOM.Modality')
    #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
    
    # Determine scan type
    if "cbct" in volumeName.lower():
        cbct_list.append(volumeName)
        scan_type = "CBCT"
    else:
        ct_list.append(volumeName)
        scan_type = "CT"
    
    # Detect points using region growing
    grouped_points = detect_points_region_growing(volumeName, intensity_threshold=3000)
    volume_points_dict[(scan_type, volumeName)] = grouped_points

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


cbct_points = [centroid for centroid, _ in volume_points_dict[("CBCT", "CBCT1")]]  # Extract only centroids (coordinates)
ct_points = [centroid for centroid, _ in volume_points_dict[("CT", "CT")]]  # Extract only centroids (coordinates)
print("CBCT points: ", np.array(cbct_points))
# Ensure we have enough points for registration
if len(cbct_points) >= 3 and len(ct_points) >= 3:
    rotation_matrix, translation_vector = compute_rigid_transform(cbct_points, ct_points)
    transformed_cbct_points = apply_transform(cbct_points, rotation_matrix, translation_vector)

    print("Optimal Rotation Matrix:\n", rotation_matrix)
    print("Translation Vector:\n", translation_vector)
    print("Transformed CBCT Points:\n", transformed_cbct_points)
else:
    print("Not enough points for registration.")