Instantly share code, notes, and snippets.

# VolkerH/feret_diameter.py

Created January 12, 2021 19:45
Show Gist options
• Save VolkerH/0d07d05d5cb189b56362e8ee41882abf to your computer and use it in GitHub Desktop.
Calculate minimum and maximum feret diameters for connected components in label images in python
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 # Volker Hilsenstein # BSD-3 license import numpy as np import skimage.morphology from rotating_calipers import min_max_feret def get_min_max_feret_from_labelim(label_im, labels=None): """ given a label image, calculate the oriented bounding box of each connected component with label in labels. If labels is None, all labels > 0 will be analyzed. Parameters: label_im: numpy array with labelled connected components (integer) Output: obbs: dictionary of oriented bounding boxes. The dictionary keys correspond to the respective labels """ if labels is None: labels = set(np.unique(label_im)) - {0} results = {} for label in labels: results[label] = get_min_max_feret_from_mask(label_im == label) return results def get_min_max_feret_from_mask(mask_im): """ given a binary mask, calculate the minimum and maximum feret diameter of the foreground object. This is done by calculating the outline of the object, transform the pixel coordinates of the outline into a list of points and then calling Parameters: mask_im: binary numpy array """ eroded = skimage.morphology.erosion(mask_im) outline = mask_im ^ eroded boundary_points = np.argwhere(outline > 0) # convert numpy array to a list of (x,y) tuple points boundary_point_list = list(map(list, list(boundary_points))) return min_max_feret(boundary_point_list)
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 # This file contains code taken from # http://code.activestate.com/recipes/117225-convex-hull-and-diameter-of-2d-point-sets/ # convex hull (Graham scan by x-coordinate) and diameter of a set of points # David Eppstein, UC Irvine, 7 Mar 2002 # According to that website the code is under the PSF licencse # https://en.wikipedia.org/wiki/Python_Software_Foundation_License # modifications by Volker Hilsenstein from __future__ import generators from math import sqrt def orientation(p,q,r): '''Return positive if p-q-r are clockwise, neg if ccw, zero if colinear.''' return (q[1]-p[1])*(r[0]-p[0]) - (q[0]-p[0])*(r[1]-p[1]) def hulls(Points): '''Graham scan to find upper and lower convex hulls of a set of 2d points.''' U = [] L = [] Points.sort() for p in Points: while len(U) > 1 and orientation(U[-2],U[-1],p) <= 0: U.pop() while len(L) > 1 and orientation(L[-2],L[-1],p) >= 0: L.pop() U.append(p) L.append(p) return U,L def rotatingCalipers(Points): '''Given a list of 2d points, finds all ways of sandwiching the points between two parallel lines that touch one point each, and yields the sequence of pairs of points touched by each pair of lines.''' U,L = hulls(Points) i = 0 j = len(L) - 1 while i < len(U) - 1 or j > 0: yield U[i],L[j] # if all the way through one side of hull, advance the other side if i == len(U) - 1: j -= 1 elif j == 0: i += 1 # still points left on both lists, compare slopes of next hull edges # being careful to avoid divide-by-zero in slope calculation elif (U[i+1][1]-U[i][1])*(L[j][0]-L[j-1][0]) > \ (L[j][1]-L[j-1][1])*(U[i+1][0]-U[i][0]): i += 1 else: j -= 1 def min_max_feret(Points): '''Given a list of 2d points, returns the minimum and maximum feret diameters.''' squared_distance_per_pair = [((p[0]-q[0])**2 + (p[1]-q[1])**2, (p,q)) for p,q in rotatingCalipers(Points)] min_feret_sq, min_feret_pair = min(squared_distance_per_pair) max_feret_sq, max_feret_pair = max(squared_distance_per_pair) return sqrt(min_feret_sq), sqrt(max_feret_sq) def diameter(Points): '''Given a list of 2d points, returns the pair that's farthest apart.''' diam,pair = max([((p[0]-q[0])**2 + (p[1]-q[1])**2, (p,q)) for p,q in rotatingCalipers(Points)]) return diam, pair def min_feret(Points): '''Given a list of 2d points, returns the pair that's farthest apart.''' min_feret_sq,pair = min([((p[0]-q[0])**2 + (p[1]-q[1])**2, (p,q)) for p,q in rotatingCalipers(Points)]) return min_feret_sq, pair