# Generating x,y coordinates for an edge detection

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.

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.

• @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, 2014 at 12:02
• @Ö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, 2014 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, 2014 at 12: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.sstatic.net/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;


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;
Prolog -> {Red, Opacity[0.5],
constraints /. {y0_ <= y[x_] <= y1_ -> Line[{{x, y0}, {x, y1}}]}}]


Or, in 3d:

radiusInterpolation = ListInterpolation[smoothRadii];


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)];


• "smoothness = Total[Differences[vars, 2]^2];" Whoa. That is cool.
– user484
Apr 27, 2014 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 :) ). Apr 28, 2014 at 13:57
• @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. Apr 28, 2014 at 15:40

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.sstatic.net/UhLrU.png"];
img = ImageTake[img, All, {1, ImageDimensions[img][[2]]/2}];
{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.

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]


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[], [email protected], Brown,
BSplineSurface[wgpts, SplineClosed -> {False, True},
SplineDegree -> {3, 2}, SplineKnots -> {knots, circKnots},
SplineWeights -> wgwts]}, Boxed -> False]