How to draw Feret diameters (CaliperLength and CaliperWidth) and calculate angle between them

I would need to extract the angle between Feret diameters (and represent them) for an image such as:

Any hint?

The first step is probably always to calculate the convex hull (the caliper length/width are always on the convex hull, and it's efficient to calculate):

img = ImageTake[
ColorNegate@
Binarize@
Import["https://i.stack.imgur.com/RNUzy.jpg"], {3, -3}, {3, -3}]

boundary = PixelValuePositions[MorphologicalPerimeter[img], 1];
{hull} = MeshPrimitives[ConvexHullMesh[boundary], 2];

HighlightImage[img, hull]


If you're not too concerned about performance, you can simply calculate the distance between each pair of points on the hull, and pick the pair that maximizes the distance. That's the caliper length:

pointDist = DistanceMatrix[hull[[1]]];
len1 = Ordering[Max /@ pointDist, -1][[1]];
len2 = Ordering[pointDist[[len1]], -1][[1]];
caliperLen = hull[[1, {len1, len2}]];


For the caliper width, you need the distance between every point and every line segment between two vertices on the hull:

direction = Normalize /@ (hull[[1]] - RotateLeft[hull[[1]]]);
Abs[d.{{0, 1}, {-1, 0}}.(Transpose[hull[[1]]] - p)]
],
{hull[[1]], direction}];

width1 = Ordering[Max /@ segmentDist, 1][[1]];
width2 = Ordering[segmentDist[[width1]], -1][[1]];
caliperWidth = hull[[1, {width1, width2}]];
d = direction[[width1]];
caliperWidth[[2]] -= (caliperWidth[[2]] - caliperWidth[[1]]).d*d;

HighlightImage[img, {hull, Line[caliperLen], Line[caliperWidth]}]


• Great! although, to be honest, why do I need to SubtractFrom: caliperWidth[[2]] -= (caliperWidth[[2]] - caliperWidth[[1]]).d*d  ? Jan 4, 2018 at 14:17
• If you look at the caliper width line, it doesn't connect two vertices; Instead, it's a line normal to the convex hull border. To get this line from two vertices, I have to subtract the part that's parallel to the hull border. Does that make sense? Jan 4, 2018 at 14:41
• Absolutely yes! understood.. Thank you! Jan 4, 2018 at 14:59
im = "your image";
im = ImageTake[im, {3, -3}, {3, -3}]; (* To remove non-black pixel in the edge *)

curve = ColorReplace[ColorConvert[im, "Grayscale"], White -> Black]


The non-black pixels of the curve image indicate the boundary, but the further their pixel value is from 0.5, the less they indicate it. Hence the following weights:

which = PixelValuePositions[Binarize[curve, 1/255/2], 1];
weights = HeavisideLambda[2 Abs[1/2 - PixelValue[curve, which]]];


In order to fit a curve to the boundary, I take out a small sample of the pixels to use for an initial guess.

ini = With[{path = Module[{remain = which}, (Remove[remain]; #) &[
NestWhileList[With[{next = Nearest[remain, #, 18]}, (* 18 is case-specific tuning *)
remain = Complement[remain, next]; Last[next]] &,
First[which], Length[remain] >= 2 &]]]},
path[[2 ;; -2]]]; (* 2 and -2 are case-specific tuning *)
Show[curve, ListLinePlot[ini - 1/2, PlotStyle -> {Red}, AspectRatio -> Automatic]]


Fitting a periodic quadratic spline to the small sample:

knots = Range[-#, 1 + #, #] &[1/Length[ini]];
sz = Length[knots] - 4;
basis = Table[BSplineBasis[{2, knots}, i, Mod[t, 1]] +
BSplineBasis[{2, knots}, i, Mod[t, 1] - 1] +
BSplineBasis[{2, knots}, i, Mod[t, 1] + 1]
// PiecewiseExpand // Simplify[#, 0 <= Mod[t, 1] < 1] &, {i, 0, sz}];

mat = SparseArray[Table[basis, {t, MovingAverage[knots[[2 ;; -2]], 2]}]];
coef = LinearSolve[mat, N[ini]];

SetAttributes[spline, Listable]
spline[t_] = basis.coef;

Show[curve, ParametricPlot[spline[t] - {1/2, 1/2}, {t, 0, 1}, PlotStyle -> {Red}, AspectRatio -> Automatic]]


The spline is adjusted iteratively with weighed total least squares on the full sample.

Do[
(* Approximately project the sample points to the current spline fit *)
ts = N[Range[0, 1 - #, #] &[1/1000]];
XYs = spline[ts];
tEstimate = ts[[Nearest[XYs -> "Index", which, 1][[All, 1]]]];

(* Total least squares *)
mat = ArrayFlatten[{{#, 0}, {0, #}}] &[Sqrt[weights] SparseArray[Table[basis, {t, tEstimate}]]];
coef = Transpose[Partition[LinearSolve[Transpose[mat].mat, Transpose[mat].Flatten[Sqrt[weights] which, {2, 1}]], Length[basis]]];
spline[t_] = basis.coef;,
{6}]

(* Final fit *)
Show[curve, ParametricPlot[spline[t] - {1/2, 1/2}, {t, 0, 1}, PlotStyle -> {Red}, AspectRatio -> Automatic]]


To compute the diameter each point on the curve is projected to a line perpendicular to the curve's direction. Given some angle, the diameter is the distance between the furthest projections on the line perpendicular to that angle.

der[t_] = (spline'[t] /. {t > Floor[t] -> True}); (* Differentiating Mod causes the need for this replacement *)

proj[t_] := With[{vec = {Sin[Mod[ArcTan @@ der[t], \[Pi]]],
-Cos[Mod[ArcTan @@ der[t], \[Pi]]]}},
Norm[#] Sign[#.vec] &[Projection[spline[t], vec]]]
plot = ParametricPlot[{Cos[2 ArcTan @@ der[t]], Sin[2 ArcTan @@ der[t]]} (proj[t] + 300), {t, 0, 1}]


The plot shows the projections as the radius while the argument correspond to the angle of direction.
Below the plot data is used to create a function from the angle of direction to the diameter by linear interpolation:

pieces = Join @@ (Partition[#, 2, 1, {1, 1}] & /@ Cases[plot, Line[___], {0, ∞}][[All, 1]]);
With[{n = Map[Norm, pieces, {2}], v = Sort /@ Mod[Apply[ArcTan, pieces, {2}], 2 π, 0]},
(maxs = Pick[v, #][[All, 2]]; norms = Pick[n, #];
diffs = Pick[-Subtract @@@ v, #]; ints = Pick[Interval /@ v, #];) &[
Thread[Differences[v, {0, 1}][[All, 1]] < 1.9 π]}]]]

dia[v_?NumericQ] := With[{memQ = IntervalMemberQ[ints, 2 v]}, If[Count[memQ, True] >= 2,
{-1, 1}.MinMax[Diagonal[Pick[norms, memQ].{#, 1 - #}] &[
(Pick[maxs, memQ] - 2 v)/Pick[diffs, memQ]]], Undefined]]
Plot[dia[v], {v, 0, π}]


The angles with minimal and maximal diameter:

angles = Join[FindArgMax[dia[v], {v, 0.1, 0.2}], FindArgMin[dia[v], {v, 2, 2.2}]]
(* {0.089053929, 2.1780192} *)


Below I find those places lineT where the spline has these directions by plotting the angle of direction w.r.t. t and looking through the plot data:

lineT = Map[With[{plot = Plot[Mod[ArcTan @@ der[t] - #, π], {t, 0, 1}, Exclusions -> None]},
With[{pieces = Partition[FirstCase[plot, Line[a___] :> a, {}, {0, ∞}], 2, 1, {1, 1}]},
Mean /@ Pick[pieces[[All, All, 1]], Thread[Abs[Differences[
pieces[[All, All, 2]], {0, 1}][[All, 1]]] > 0.9 π]]]] &, angles] // Join @@ # &

Show[curve, ParametricPlot[spline[t] - {1/2, 1/2}, {t, 0, 1},
PlotStyle -> {Red}, AspectRatio -> Automatic],
Graphics[{Red, InfiniteLine[spline[#] - {1/2, 1/2}, der[#]] & /@ Delete[lineT, {5}]}]]