4
$\begingroup$

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

sample image

Any hint?

$\endgroup$

2 Answers 2

6
$\begingroup$

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]

enter image description here

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]]]);
segmentDist = MapThread[Function[{p, d},
    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]}]

enter image description here

$\endgroup$
3
  • $\begingroup$ Great! although, to be honest, why do I need to SubtractFrom: caliperWidth[[2]] -= (caliperWidth[[2]] - caliperWidth[[1]]).d*d ? $\endgroup$
    – Stefano
    Jan 4, 2018 at 14:17
  • 1
    $\begingroup$ 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? $\endgroup$ Jan 4, 2018 at 14:41
  • $\begingroup$ Absolutely yes! understood.. Thank you! $\endgroup$
    – Stefano
    Jan 4, 2018 at 14:59
3
$\begingroup$
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, #];) &[
 MapThread[And, {Thread[-Subtract @@@ v != 0],
                 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}]}]]

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.