In this notebook our goal is to estimate the radius of spherical markers from an image (Cone-Beam CT volume).
We will use two approaches:
It should be noted that all of the operations, filtering and computations, are natively in 3D. This is the "magic" of ITK and SimpleITK at work.
In [ ]:
import SimpleITK as sitk
from downloaddata import fetch_data as fdata
import matplotlib.pyplot as plt
%matplotlib inline
import numpy as np
from scipy import linalg
from ipywidgets import interact, fixed
Load the volume and look at the image (visualization requires window-leveling).
In [ ]:
spherical_fiducials_image = sitk.ReadImage(fdata("spherical_fiducials.mha"))
sitk.Show(spherical_fiducials_image, "spheres")
After looking at the image you should have identified two spheres. Now select a Region Of Interest (ROI) around the sphere which you want to analyze.
In [ ]:
rois = {"ROI1":[range(280,320), range(65,90), range(8, 30)],
"ROI2":[range(200,240), range(65,100), range(15, 40)]}
mask_value = 255
def select_roi_dropdown_callback(roi_name, roi_dict):
global mask, mask_ranges
mask_ranges = roi_dict.get(roi_name)
if mask_ranges:
mask = sitk.Image(spherical_fiducials_image.GetSize(), sitk.sitkUInt8)
mask.CopyInformation(spherical_fiducials_image)
for x in mask_ranges[0]:
for y in mask_ranges[1]:
for z in mask_ranges[2]:
mask[x,y,z] = mask_value
# Use nice magic numbers for windowing the image. We need to do this as we are alpha blending the mask
# with the original image.
sitk.Show(sitk.LabelOverlay(sitk.Cast(sitk.IntensityWindowing(spherical_fiducials_image, windowMinimum=-32767,
windowMaximum=-29611), sitk.sitkUInt8),
mask, opacity=0.5))
roi_list = rois.keys()
roi_list.insert(0,'Select ROI')
interact(select_roi_dropdown_callback, roi_name=roi_list, roi_dict=fixed(rois));
In [ ]:
intensity_values = sitk.GetArrayFromImage(spherical_fiducials_image)
roi_intensity_values = intensity_values[mask_ranges[2][0]:mask_ranges[2][-1],
mask_ranges[1][0]:mask_ranges[1][-1],
mask_ranges[0][0]:mask_ranges[0][-1]].flatten()
plt.hist(roi_intensity_values, bins=100)
plt.title("Intensity Values in ROI")
plt.show()
Can you identify the region of the histogram associated with the sphere?
In our case it looks like we can automatically select a threshold separating the sphere from the background. We will use Otsu's method for threshold selection to segment the sphere and estimate its radius.
In [ ]:
# Set pixels that are in [min_intensity,otsu_threshold] to inside_value, values above otsu_threshold are
# set to outside_value. The sphere's have higher intensity values than the background, so they are outside.
inside_value = 0
outside_value = 255
number_of_histogram_bins = 100
mask_output = True
labeled_result = sitk.OtsuThreshold(spherical_fiducials_image, mask, inside_value, outside_value,
number_of_histogram_bins, mask_output, mask_value)
# Estimate the sphere radius from the segmented image using the LabelShapeStatisticsImageFilter.
label_shape_analysis = sitk.LabelShapeStatisticsImageFilter()
label_shape_analysis.SetBackgroundValue(inside_value)
label_shape_analysis.Execute(labeled_result)
print("The sphere's radius is: {0:.2f}mm".format(label_shape_analysis.GetEquivalentSphericalRadius(outside_value)))
In [ ]:
# Visually inspect the results of segmentation, just to make sure.
sitk.Show(sitk.LabelOverlay(sitk.Cast(sitk.IntensityWindowing(spherical_fiducials_image, windowMinimum=-32767, windowMaximum=-29611),
sitk.sitkUInt8), labeled_result, opacity=0.5))
Based on your visual inspection, did the automatic threshold correctly segment the sphere or did it over/under segment it?
If automatic thresholding did not provide the desired result, you can correct it by allowing the user to modify the threshold under visual inspection. Implement this approach below.
In [ ]:
# Your code here:
In this approach we will localize the sphere's edges in 3D using SimpleITK. We then compute the least squares sphere that optimally fits the 3D points using scipy/numpy. The mathematical formulation for this solution is described in this Insight Journal paper.
In [ ]:
# Create a cropped version of the original image.
sub_image = spherical_fiducials_image[mask_ranges[0][0]:mask_ranges[0][-1],
mask_ranges[1][0]:mask_ranges[1][-1],
mask_ranges[2][0]:mask_ranges[2][-1]]
# Edge detection on the sub_image with appropriate thresholds and smoothing.
edges = sitk.CannyEdgeDetection(sitk.Cast(sub_image, sitk.sitkFloat32), lowerThreshold=0.0,
upperThreshold=200.0, variance = (5.0,5.0,5.0))
Get the 3D location of the edge points and fit a sphere to them.
In [ ]:
edge_indexes = np.where(sitk.GetArrayFromImage(edges) == 1.0)
# Note the reversed order of access between SimpleITK and numpy (z,y,x)
physical_points = [edges.TransformIndexToPhysicalPoint([int(x), int(y), int(z)]) \
for z,y,x in zip(edge_indexes[0], edge_indexes[1], edge_indexes[2])]
# Setup and solve linear equation system.
A = np.ones((len(physical_points),4))
b = np.zeros(len(physical_points))
for row, point in enumerate(physical_points):
A[row,0:3] = -2*np.array(point)
b[row] = -linalg.norm(point)**2
res,_,_,_ = linalg.lstsq(A,b)
print("The sphere's radius is: {0:.2f}mm".format(np.sqrt(linalg.norm(res[0:3])**2 - res[3])))
In [ ]:
# Visually inspect the results of edge detection, just to make sure. Note that because SimpleITK is working in the
# physical world (not pixels, but mm) we can easily transfer the edges localized in the cropped image to the original.
edge_label = sitk.Image(spherical_fiducials_image.GetSize(), sitk.sitkUInt16)
edge_label.CopyInformation(spherical_fiducials_image)
e_label = 255
for point in physical_points:
edge_label[edge_label.TransformPhysicalPointToIndex(point)] = e_label
sitk.Show(sitk.LabelOverlay(sitk.Cast(sitk.IntensityWindowing(spherical_fiducials_image, windowMinimum=-32767, windowMaximum=-29611),
sitk.sitkUInt8), edge_label, opacity=0.5))
In [ ]: