How to get the long & short axis OR get length of mask at a point orthogonal to the feret diameter?

I'm trying to measure a polygon mask's longest diameter (feret diameter) while also getting the length of the orthogonal line from the center of the feret diameter. Diagram of what I'm trying to do can be found here: https://radiologyassistant.nl/assets/fleischner-2017-guideline-for-pulmonary-nodules/a59385c3993fd3_TAB-measurement.png

I cannot use the major_axis and minor_axis properties of skimage.measure.regionprops as that fits an ellipse over my mask and is found to be wildly unreliable for the mask shapes I have.

Are there any functions or other ways I can do this? I've only found skimage's methods to be most applicable for what I am trying to do. I don't know if there are other terms for what I'm trying to do. Any insight is helpful.

This is the basis of what I'm doing:

mask_ps = np.transpose(skimage.draw.polygon2mask(im2d_ps.shape, this_ps_seg_points)) #turn into a polygon mask

mask_ps_label = skimage.measure.label(mask_ps)
props_ps = skimage.measure.regionprops(mask_ps_label)
props_ps.sort(key=lambda x: x.area)
print("props_ps.area", props_ps[-1].area) #max area
this_feret = props_ps[-1].feret_diameter_max

The mask that I'm using (mask_ps) is quite large so I can't insert it here by copy/paste.

I really like the idea that Christoph suggested in the comments.

Here is some code implementing that in Python, using SimpleITK.

import SimpleITK as sitk
import numpy as np

import matplotlib.pyplot as plt

# create an empty array
arr = np.zeros((512,512), dtype=bool)
xx,yy = np.indices(arr.shape)


# add a bunch of overlapping circles
cond = (xx-256)**2 + (yy-256)**2 < 50**2
arr[cond] = 1

cond = (xx-300)**2 + (yy-300)**2 < 50**2
arr[cond] = 1

cond = (xx-210)**2 + (yy-200)**2 < 80**2
arr[cond] = 1

cond = (xx-300)**2 + (yy-190)**2 < 90**2
arr[cond] = 1

# display
fig, ax = plt.subplots(1,1,figsize=(5,5))
ax.imshow(arr, interpolation=None, cmap=plt.cm.Greys_r)
fig.show()

# create a SimpleITK image
img = sitk.GetImageFromArray(arr.astype(int))

# generate label 
filter_label = sitk.LabelShapeStatisticsImageFilter()
filter_label.SetComputeFeretDiameter(True)
filter_label.Execute(img)

# compute the Feret diameter
# the 1 means we are computing for the label with value 1
filter_label.GetFeretDiameter(1)

Returns:

264.25177388240934

# we have to get a bit smarter for the principal moments
pc1_x, pc1_y, pc2_x, pc2_y = filter_label.GetPrincipalAxes(1)

# get the center of mass
com_y, com_x = filter_label.GetCentroid(1)

# now trace the distance from the centroid to the edge along the principal axes
# we use some linear algebra

# get the position of each point in the image
v_x, v_y = np.where(arr)

# convert these positions to a vector from the centroid
v_pts = np.array((v_x - com_x, v_y - com_y)).T

# project along the first principal component
distances_pc1 = np.dot(v_pts, np.array((pc1_x, pc1_y)))

# get the extent
dmax_1 = distances_pc1.max()
dmin_1 = distances_pc1.min()

# project along the second principal component
distances_pc2 = np.dot(v_pts, np.array((pc2_x, pc2_y)))

# get the extent
dmax_2 = distances_pc2.max()
dmin_2 = distances_pc2.min()

# the total diameter is the difference in these distances
print("Distance along major axis:", dmax_1 - dmin_1)
print("Distance along minor axis:", dmax_2 - dmin_2)

Distance along major axis: 258.5497770889507

Distance along minor axis: 251.97123773093475

# display
fig, ax = plt.subplots(1,1,figsize=(5,5))
ax.imshow(arr, interpolation=None, cmap=plt.cm.Greys_r)

ax.scatter(com_y, com_x, c="g", marker="o", s=50, zorder=99, label="COM")

ax.plot((com_y, com_y+dmax_1*pc1_y), (com_x, com_x+dmax_1*pc1_x), lw=2, c='b')
ax.plot((com_y, com_y+dmin_1*pc1_y), (com_x, com_x+dmin_1*pc1_x), lw=2, c='b', label="Major axis")

ax.plot((com_y, com_y+dmax_2*pc2_y), (com_x, com_x+dmax_2*pc2_x), lw=2, c='r')
ax.plot((com_y, com_y+dmin_2*pc2_y), (com_x, com_x+dmin_2*pc2_x), lw=2, c='r', label="Minor axis")

ax.legend()

fig.show()

You'll see that the projected distances along the major/minor axes are not necessarily equal to the maximum (Feret) diameter. If you are truly only interested in these, then I'd recommend checking out the feret Python package. In any case, I hope this code gives you something helpful! Also, SimpleITK is actually really great for lots of image processing tasks, and is a nice supplement to SciKit-Image.