This code localizes fiducials given corresponding fiducial detection results.

localize_fiducials.py

import nibabel as nib
import pcl
import numpy as np
from pcl import IterativeClosestPoint
from pcl import pcl_visualization


# used for icp rigid alignemnt of fiducial candidate with synthetic fiducial for lits localization
def generate_synthetic_fiducial():
   
    # init array 
    synthetic_fiducial = [] 
    global synthetic_fiducial_pcl_cloud 
    synthetic_fiducial_pcl_cloud = pcl.PointCloud() 

    # config of fiducial
    # reference: http://izimed.com/all-products/multi-modality-fiducial-markers-mm30xx.html
    global height_fid_mm 
    height_fid_mm = 3.5 
    inner_rad_fid_mm = 2.0 # not specified
    outer_rad_fid_mm = 7.5 
    # resolution at which to generate representation
    global header 
    header = segmentation_result.header
    x_resolution = header.get_zooms()[1]
    y_resolution = header.get_zooms()[2]
    z_resolution = header.get_zooms()[0]
    point_id = 0 
    # loop to generate synthetic fiducial model 
    for x in np.arange(-1.0 * outer_rad_fid_mm, outer_rad_fid_mm, x_resolution /2 ):
        for z in np.arange(-1.0 * outer_rad_fid_mm, outer_rad_fid_mm , z_resolution/ 2):
            for y in np.arange(-1.0 * height_fid_mm  / 2.0, height_fid_mm / 2.0 , y_resolution/2 ):
                if (( inner_rad_fid_mm * inner_rad_fid_mm <= x * x + z * z <= outer_rad_fid_mm * outer_rad_fid_mm) and (y == height_fid_mm / 2.0 or y == -1.0 * height_fid_mm / 2.0)): # on the top/botton surface
                   synthetic_fiducial.append([x, y,  z]) #why linear translation required, TODO

    # write the point cloud back
    synthetic_fiducial_pcl_cloud.from_array(np.asarray(synthetic_fiducial, dtype=np.float32))
    pcl.save(synthetic_fiducial_pcl_cloud, 'synthetic_fiducial.pcd')


# localizes fiducials segmented by NN solution
def localize_fiducials():

    # get segmentation data as numpy array
    segmentation_result_np_array = segmentation_result.get_fdata()
    print ('shape of segmentation result array: ', segmentation_result_np_array.shape) # debug
    x_dim = segmentation_result_np_array.shape[0]
    y_dim = segmentation_result_np_array.shape[1]
    z_dim = segmentation_result_np_array.shape[2]

    segmentation_points_of_interest =  [(x, y, z) for x in range(x_dim) for y in range(y_dim) for z in range(z_dim)  if segmentation_result_np_array[x][y][z] == 1]
    
    # compute connected components:
    segmentation_result_point_cloud = pcl.PointCloud(segmentation_points_of_interest) 


    print ('shape of array of point of interest: ', np.asarray(segmentation_points_of_interest).shape)
     
    pcl.save(segmentation_result_point_cloud, 'only_fiducial_segmentations.pcd') 
    
    # noise filtering and outlier removal TODO

    # euclidean clustering 
    vg = segmentation_result_point_cloud.make_voxel_grid_filter()
    vg.set_leaf_size(1, 1, 1)
    cloud_filtered = vg.filter()
    pcl.save(cloud_filtered, 'after_voxel_grid_filter.pcd')
    tree = cloud_filtered.make_kdtree()

    connected_components = cloud_filtered.make_EuclideanClusterExtraction()
    connected_components.set_ClusterTolerance(0.02)
    connected_components.set_MinClusterSize(5)
    connected_components.set_MaxClusterSize(15000)
    connected_components.set_SearchMethod(tree)
    cluster_indices = connected_components.Extract()
    print ('number of clusters', len(cluster_indices))
    

    for j, indices in enumerate(cluster_indices):
        cloud_cluster = pcl.PointCloud()
        print ('indices = ', str(len(indices)))
        points = np.zeros((len(indices), 3), dtype = np.float32)
 
        for i, indice in enumerate(indices):
            points[i][0] = cloud_filtered[indice][0]
            points[i][1] = cloud_filtered[indice][1]
            points[i][2] = cloud_filtered[indice][2]
        cloud_cluster.from_array(points) # collects point-clouds corresponding to each cluster
        # add original segmentation
        pcl.save(cloud_cluster, str(j) + '_cluster.pcd')
        icp = cloud_cluster.make_IterativeClosestPoint()
        # registration process
        converged, transf, estimate, fitness = icp.icp(cloud_cluster, synthetic_fiducial_pcl_cloud, max_iter=1000)
        # localization
        print('fitness of localization: ', fitness)
        print('icp converged?: ', converged)

        point_of_interest_on_synthetic_fiducial = [0, height_fid_mm / 2.0, 0, 1]
        coordinates_of_point_of_interest_fiducial_image_coordinates = np.dot (transf, point_of_interest_on_synthetic_fiducial) # in image coordinate system

        for i in range(2):
            coordinates_of_point_of_interest_fiducial_image_coordinates[i] /= coordinates_of_point_of_interest_fiducial_image_coordinates[2]

        coordinates_of_point_of_interest_fiducial_patient_coordinates = np.dot(segmentation_result.affine, coordinates_of_point_of_interest_fiducial_image_coordinates) # in patient coordinate system
        for i in range(2):
            coordinates_of_point_of_interest_fiducial_patient_coordinates[i] /= coordinates_of_point_of_interest_fiducial_patient_coordinates[2]

        # convert from RAS+ to LPS+
        for i in range(2):
            coordinates_of_point_of_interest_fiducial_patient_coordinates[i] *= -1.0 * coordinates_of_point_of_interest_fiducial_patient_coordinates[i]
        # z-axis stays the same.    

       
        # save per scan/segmentation   
        np.savetxt('point_of_interest_coordinates_' + str(j) + '.txt', coordinates_of_point_of_interest_fiducial_patient_coordinates)   
        
'''        viewer = pcl.pcl_visualization.PCLVisualizering('ground_cluster' + str(i))
        viewer.AddPointCloud(segmentation_result_point_cloud)
        viewer.AddPointCloud(cloud_cluster)
'''

if __name__ == '__main__':
    
    # load segmentation results
    global segmentation_result 
    segmentation_result = nib.load('test.nii.gz')

    # DEBUG
    print('image header: ', segmentation_result.header)
    print('check version: ', segmentation_result) # to test whether nii1 or nii2
    print('affine matrix: ', segmentation_result.affine)

    generate_synthetic_fiducial()
    localize_fiducials()