5

I'd like to be able to process a number of graphics that have edge contours and output a list of $x$, $y$ coordinates to be used for generating a spline in another program, or in Mathematica, to make a SOR, prism, etc.

Essentially, given a photo of a turned spindle, output the radius as a function of length.

enter image description here

I was working with code from Get x and z coordinate from an image and make a parametric surface. I've got to the point where I can output a graph of the "height" vs length, but I need to get a list of coordinates.

I suppose being able to draw a line to measure from in the case of tricky graphics or slanted [eventual] centers of rotation would be a helpful feature.

Bald Eagle
  • 51
  • 1
  • 3
  • @kuba, are you planning on answering this question? I believe that it's a duplicate of the Q&A linked in the question. The OP is "only" asking for the coordinates, which can be roughly given by ssch' answer. – Öskå Apr 26 '14 at 12:02
  • 1
    @Öskå No, I'm not. I just though it may be of interest of future visitors and someone who like image processing will answer so I've edited information provided by OP. – Kuba Apr 26 '14 at 12:04
  • @Kuba Well, IMO it's a duplicate. Combining f[y] from ssch's answer and using J.M methods give a decent 3D plot. – Öskå Apr 26 '14 at 12:15

2 Answers2

15

For a simple image like that (high-contrast object in front of a monochrome background), calling EdgeDetect is actually enough to find the edges:

img = Import["https://i.stack.imgur.com/UhLrU.png"];    
edges = EdgeDetect[img];

Then we can find the leftmost and rightmost points on every line to get an array of radii:

radii = Map[Max[#[[All, 1]]] - Min[#[[All, 1]]] &, 
    GatherBy[PixelValuePositions[edges, 1], Last]]/2;
ListLinePlot[radii]

enter image description here

This line is "jagged", because every edge position was essentially "rounded" to the nearest integer. What we really want is a line that is

  • as smooth as possible
  • but nowhere more than 0.5px off these radii.

That's an optimization problem:

n = Length[radii];
vars = Array[y, n];
constraints = Array[radii[[#]] - 0.5 <= y[#] <= radii[[#]] + 0.5 &, n];
smoothness = Total[Differences[vars, 2]^2];

{fit, sol} = FindMinimum[{smoothness, constraints}, vars];

We can plot these smoothed values against the original "jagged" values:

smoothRadii = vars /. sol;    
ListLinePlot[smoothRadii, 
 Prolog -> {Red, Opacity[0.5], 
   constraints /. {y0_ <= y[x_] <= y1_ -> Line[{{x, y0}, {x, y1}}]}}]

enter image description here

Or, in 3d:

radiusInterpolation = ListInterpolation[smoothRadii];
RevolutionPlot3D[{radiusInterpolation[u], u}, {u, 1, 
  Length[smoothRadii]}, PlotPoints -> 50]

enter image description here

We can even do (a bit) better using a smoothness measure that doesn't penalize "kinks" as hard:

smoothness = Total[Clip[#, {0, .5}] & /@ (Differences[vars, 2]^2)];

enter image description here

Niki Estner
  • 36,101
  • 3
  • 92
  • 152
  • 1
    "smoothness = Total[Differences[vars, 2]^2];" Whoa. That is cool. –  Apr 27 '14 at 18:33
  • @nikie Overall, nice approach, but I think you need to clean shadows on the left hand side of the image to get it right size (min max includes the shadow :) ). – s.s.o Apr 28 '14 at 13:57
  • 2
    @s.s.o: Oops, you're right; I didn't realize there was a shadow. On the other hand, I just answered this question because I didn't notice it was months old (and probably abandoned by the OP). Adding "image processing cleanup code" that will only work for this single image anyway seems kind of pointless now. – Niki Estner Apr 28 '14 at 15:40
7

Just to show that this question is a duplicate of the link given in its body:

By using ssch's answer one can get the coordinates of the shape:

img = Import["https://i.stack.imgur.com/UhLrU.png"];
img = ImageTake[img, All, {1, ImageDimensions[img][[2]]/2}];
img = ColorNegate[GradientFilter[img, 4]];
{width, height} = ImageDimensions[img];
width = width - 3;
f[y_] := Module[{i = 1}, 
  While[ImageValue[img, {i, y}] > .99 && i < width, i = i + 1]; width - i + 1]

Where f[y] gives the radius in terms of the height. Some quick operations can be done in order to obtain a good looking shape:

coor[n_] := Select[Reverse /@ Sort@({#1, f[#1]} & /@ Range[1,height,n]), 5 < First@# < 50 &]
shapePos = coor@10

where n would represent the sharpness/precision of the shape.

Then by using J.M.'s answer:

parametrizeCurve[pts_List, a : (_?NumericQ) : 1/2] := 
  FoldList[Plus, 0, Normalize[(Norm /@ Differences[pts])^a, Total]] /;MatrixQ[pts, NumericQ]
tvals = N@parametrizeCurve[shapePos];
m = 3;
knots = Join[ConstantArray[0, m + 1], MovingAverage[ArrayPad[tvals, -1], m], ConstantArray[1, m + 1]];
bas = Table[BSplineBasis[{m, knots}, j - 1, tvals[[i]]],
  {i,Length[shapePos]}, {j, Length[shapePos]}];
ctrlpts = LinearSolve[bas, shapePos];
Graphics[{BSplineCurve[ctrlpts, SplineDegree -> 3, SplineKnots -> knots],
  {AbsolutePointSize[4], Point[shapePos]}}, Axes -> None, Frame -> True]

enter image description here

circPoints = {{1, 0}, {1, 1}, {-1, 1}, {-1, 0}, {-1, -1}, {1, -1}, {1, 0}};
circKnots = {0, 0, 0, 1/4, 1/2, 1/2, 3/4, 1, 1, 1};
circWts = {1, 1/2, 1/2, 1, 1/2, 1/2, 1};
wgpts = Map[Function[pt, Append[#1 pt, #2]], circPoints] & @@@ ctrlpts;

wgwts = ConstantArray[circWts, Length[ctrlpts]];
Graphics3D[{EdgeForm[], Opacity@.75, Brown, 
  BSplineSurface[wgpts, SplineClosed -> {False, True}, 
  SplineDegree -> {3, 2}, SplineKnots -> {knots, circKnots}, 
  SplineWeights -> wgwts]}, Boxed -> False]

enter image description here

Öskå
  • 8,587
  • 4
  • 30
  • 49