SimpleITK_OrientedBBoxExample.py

import numpy as np
import matplotlib.pyplot as plt
import SimpleITK as sitk
from skimage.draw import ellipsoid
from math import ceil

# create our label shape
ell = ellipsoid(6, 10, 16).astype(np.int32)
ell = np.pad(ell, (16,0),'constant')
labels = sitk.GetImageFromArray(ell)

# get label shape stats
shape_stats = sitk.LabelShapeStatisticsImageFilter()
shape_stats.SetComputeOrientedBoundingBox(True)
shape_stats.Execute(labels)

# show that we can use the oriented bounding box output correctly
resampler = sitk.ResampleImageFilter()
aligned_image_spacing = [1,1,1]
aligned_image_size = [int(ceil(shape_stats.GetOrientedBoundingBoxSize(1)[i] / aligned_image_spacing[i])) for i in
                      range(3)]
direction_mat = shape_stats.GetOrientedBoundingBoxDirection(1)
aligned_image_direction = [direction_mat[0], direction_mat[3], direction_mat[6],
                           direction_mat[1], direction_mat[4], direction_mat[7],
                           direction_mat[2], direction_mat[5], direction_mat[8]]
resampler.SetOutputDirection(aligned_image_direction)
resampler.SetOutputOrigin(shape_stats.GetOrientedBoundingBoxOrigin(1))
resampler.SetOutputSpacing(aligned_image_spacing)
resampler.SetSize(aligned_image_size)

bounded_img = resampler.Execute(labels)
bounded_array = sitk.GetArrayFromImage(bounded_img)
fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(20, 10))
center = np.array(bounded_array.shape) // 2
ax1.imshow(bounded_array[center[0], :, :])
ax2.imshow(bounded_array[:, center[1], :])
ax3.imshow(bounded_array[:, :, center[2]])


# Show that we can correctly plot centroid in real world coordinates
c = shape_stats.GetCentroid(1)
c = np.flip(np.array(c)).astype(np.int32)
fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(20, 10))
center = c
ax1.imshow(ell[center[0], :, :])
ax2.imshow(ell[:, center[1], :])
ax3.imshow(ell[:, :, center[2]])

ax1.scatter(c[2], c[1], s=100, marker='+', c='green')  # c=annot_c, cmap=cmap)
ax2.scatter(c[2], c[0], s=100, marker='+', c='green')
ax3.scatter(c[1], c[0], s=100, marker='+', c='green')

# Show that the oriented bounding box vertices are incorrect as in real world coordinates
oriented_bbox_vertices = shape_stats.GetOrientedBoundingBoxVertices(1)
oriented_bbox_vertices = np.array(oriented_bbox_vertices).reshape(8,3)
oriented_bbox_vertices_in_voxels = np.zeros_like(oriented_bbox_vertices)
for i,o in enumerate(oriented_bbox_vertices):
    oriented_bbox_vertices_in_voxels[i] = labels.TransformPhysicalPointToContinuousIndex(o)

for o in oriented_bbox_vertices_in_voxels:
    o = np.flip(o)
    ax1.scatter(o[2], o[1], s=200, marker='o', c='red')  # c=annot_c, cmap=cmap)
    ax2.scatter(o[2], o[0], s=200, marker='o', c='red')
    ax3.scatter(o[1], o[0], s=200, marker='o', c='red')

# now transform using the direction matrix, see that the coordinates look great
T = sitk.Similarity3DTransform()
T.SetMatrix(aligned_image_direction)
T.SetCenter(shape_stats.GetOrientedBoundingBoxOrigin(1))

oriented_bbox_vertices_rotated = np.zeros_like(oriented_bbox_vertices)
for i,o in enumerate(oriented_bbox_vertices):
    oriented_bbox_vertices_rotated[i] = T.TransformPoint(o)

oriented_bbox_vertices_rotated_in_voxels = np.zeros_like(oriented_bbox_vertices_rotated)
for i, o in enumerate(oriented_bbox_vertices_rotated):
    oriented_bbox_vertices_rotated_in_voxels[i] = labels.TransformPhysicalPointToContinuousIndex(o)

for o in oriented_bbox_vertices_rotated_in_voxels:
    o = np.flip(o)
    ax1.scatter(o[2], o[1], s=200, marker='o', c='green')  # c=annot_c, cmap=cmap)
    ax2.scatter(o[2], o[0], s=200, marker='o', c='green')
    ax3.scatter(o[1], o[0], s=200, marker='o', c='green')