2

I have a numpy array of an image:

enter image description here

The red is the binary mask over a specific part and my goal is to calculate the biggest diameter length in this mask. I have tried using a bounding box and measuring the diameter of the box (from one end to the other) but it is not as accurate as I want it to be, especially if the shape of the mask is close to a circle.

My next idea is to use a bounding circle instead of a box, and measure the radius which would be very accurate. The problem is that I don't know of any packages or functions to give me a bounding circle. Even if there is such a function, I don't know how to measure the radius of the given circle as there is no sharp points in a circle to go end to end.

Appreciate any ideas.

Raein Hashemi
  • 3,346
  • 4
  • 22
  • 33
  • You're looking to compute the maximum Feret diameter (maximum projection length) of the shape. I explain an algorithm here: https://www.crisluengo.net/archives/408 – Cris Luengo Sep 18 '19 at 20:04
  • probably duplicated with this question, https://stackoverflow.com/questions/13240615/width-of-an-arbitrary-polygon – zihaozhihao Sep 18 '19 at 20:50
  • I think given the relatively small size of the perimeters (in number of pixels) the most efficient solution might end up being a brute force one, i.e. compute some ordering on the perimeter pixels p_1, ..., p_n and compare all pairs. Worst case time O(n^2) where n is the total size of the region (long skinny region = all pixels are on the perimeter) – Sean Sep 20 '19 at 16:49

2 Answers2

2

enter image description here

Some notes on my solution:

(0) It accepts a binary image as input.

(1) It finds it for all regions in the image.

(2) It finds the greatest diameter of the convex hull of the region. I think this is a reasonable thing to do, but you can adjust the implementation if you like.

(3) I use the skimage.data "coins" image so you can reproduce yourself.

import sys

# To find the diameters
from skimage.measure import regionprops, label
from sklearn.metrics import pairwise_distances
from scipy import ndimage as ndi
import numpy as np

# To generate test data
from skimage import data
from skimage.filters import sobel
from skimage.morphology import watershed

# For visualization
import matplotlib.pyplot as plt

STREL_4 = np.array([[0, 1, 0],
                    [1, 1, 1],
                    [0, 1, 0]], dtype=np.bool)


def get_border_image(region):
    convex_hull_mask = region.convex_image
    eroded_image = ndi.binary_erosion(convex_hull_mask, STREL_4, border_value=0)
    border_image = np.logical_xor(convex_hull_mask, eroded_image)
    return border_image

def get_region_diameters(img):

    assert img.dtype == np.bool and len(img.shape) == 2

    label_img = label(img, connectivity=img.ndim)

    for region in regionprops(label_img):
        border_image = get_border_image(region)
        perimeter_coordinates = np.transpose(np.nonzero(border_image))
        pairwise_distances_matrix = pairwise_distances(perimeter_coordinates)
        i, j = np.unravel_index(np.argmax(pairwise_distances_matrix), pairwise_distances_matrix.shape)
        ptA, ptB = perimeter_coordinates[i], perimeter_coordinates[j]
        region_offset = np.asarray([region.bbox[0], region.bbox[1]])
        ptA += region_offset
        ptB += region_offset
        yield pairwise_distances_matrix[i, j], ptA, ptB

if __name__ == "__main__":

    # Create a segmentation of the coins image, for testing purposes.
    # You should create a binary image
    coins = data.coins()
    elevation_map = sobel(coins)
    markers = np.zeros_like(coins)
    markers[coins < 30] = 1
    markers[coins > 150] = 2
    segmentation = (watershed(elevation_map, markers) > 1)

    for distance, ptA, ptB in get_region_diameters(segmentation):
        plt.imshow(segmentation)
        x1, x2, y1, y2 = ptA[1], ptB[1], ptA[0], ptB[0]
        plt.plot([x1, x2], [y1, y2], color='k', linestyle='-', linewidth=2)
        print(distance, ptA, ptB)

    plt.show()
Kyle Pena
  • 517
  • 2
  • 8
  • Some notes on my solution: (1) It finds it for all regions in the image. (2) It finds the greatest diameter of the convex hull of the region. I think this is a reasonable thing to do, but you can adjust the implementation if you like. – Kyle Pena Sep 20 '19 at 16:42
  • Edited as @CrisLuengo suggests – Kyle Pena Sep 20 '19 at 16:45
  • Thank you for making the effort, appreciate it. Your logic seems accurate but there is too much code for a simple task like this. Especially when I saw you use binary erosion and XOR just to get the border which seems unnecessary. I finally ended up using openCV (cv2) as @Dr.H.Lecter suggested. Getting the contour with "contours, hierarchy = cv2.findContours(image, 1, 2)" and getting the radius of the bounding circle using "(x,y),radius = cv2.minEnclosingCircle(contours[0])" – Raein Hashemi Sep 20 '19 at 17:21
0

If you are willing to accept an approximation use this method and use the property "major_axis_length": https://scikit-image.org/docs/dev/api/skimage.measure.html#skimage.measure.regionprops

Kyle Pena
  • 517
  • 2
  • 8