# Get x and z coordinate from an image and make a parametric surface

In this blog entry, this step is mentioned:

I grabbed an image of a Peugeot wine glass from the web and used the Get Coordinates function to digitize its outline, which I fed to Interpolation to get the x and z parametric functions of the curve.

While I understand that I can get the coordinates of x and y, how can I get the z coordinates, and convert that to a Parametric Curve?

-
I think you just need to replace y by z in your statement; i.e. since he will generate a surface of revolution what you call y is in fact z. – chris Oct 6 '12 at 15:12
@chris is right; think of the outline of the wine glass as being embedded in $x$-$z$ plane. – J. M. Oct 6 '12 at 15:32
I did this once for a poster, let me dig up the code – ssch Feb 9 '13 at 1:16

(I might as well.)

wineGlassPoints = {{32.75, 283.75}, {37.75, 275.75}, {43.25, 267.25}, {49.25, 256.75},
{53.75, 247.75}, {58.25, 236.25}, {61.75, 224.25}, {64.25, 211.75},
{65.25, 198.75}, {64.25, 185.75}, {61.75, 174.25}, {58.25, 165.75},
{53.25, 157.25}, {46.25, 149.25}, {38.75, 143.75}, {30.25, 139.75},
{23.25, 137.25}, {17.75, 135.25}, {13.75, 134.75}, {13.75, 128.75},
{10.25, 126.75}, {7.25, 122.25}, {5.75, 115.75}, {4.75, 109.25},
{4.25, 101.25}, {3.75, 88.25}, {3.75, 70.25}, {3.75, 53.75},
{4.75, 39.75}, {6.75, 25.75}, {8.25, 20.25}, {11.25, 16.75},
{17.25, 14.25}, {25.25, 11.25}, {33.75, 9.25}, {41.25, 7.75},
{48.25, 5.75}, {48.75, 0.25}};


With these points, we can use Eugene Lee's centripetal parametrization method, to generate parameter values corresponding to the points:

parametrizeCurve[pts_List, a : (_?NumericQ) : 1/2] :=
FoldList[Plus, 0, Normalize[(Norm /@ Differences[pts])^a, Total]] /; MatrixQ[pts, NumericQ]

tvals = parametrizeCurve[wineGlassPoints]
{0, 0.0274652, 0.0559174, 0.0870137, 0.115379, 0.146802, 0.178417, 0.210343,
0.242632, 0.27492, 0.305596, 0.332707, 0.360788, 0.389942, 0.417212, 0.44462,
0.468999, 0.490631, 0.508584, 0.530488, 0.548441, 0.569236, 0.592332, 0.615263,
0.64058, 0.672833, 0.71077, 0.747093, 0.780593, 0.814221, 0.835571, 0.85477,
0.877568, 0.903705, 0.930129, 0.954859, 0.978986, 1.}


From here, we can easily build an InterpolatingFunction[] corresponding to the wine glass's outline:

wineGlassFunction = Interpolation[Transpose[{tvals, wineGlassPoints}]];


Have a look at the outline:

ParametricPlot[wineGlassFunction[t], {t, 0, 1},
Epilog -> {AbsolutePointSize[4], Point /@ wineGlassPoints}]


It's not too hard to embed this curve in the $x$-$z$ plane; just insert a $0$ as the second ($y$) component:

ParametricPlot3D[Insert[wineGlassFunction[t], 0, 2], {t, 0, 1}]


One could certainly use RevolutionPlot3D[] to generate the corresponding surface of revolution, but I choose to use ParametricPlot3D[] and RotationTransform[] for illustrative purposes (I also take the opportunity to give the surface a little flair):

ParametricPlot3D[RotationTransform[θ, {0, 0, 1}][Insert[wineGlassFunction[t], 0, 2]],
{t, 0, 1}, {θ, -π, π}, Axes -> None, Boxed -> False,
Lighting -> "Neutral", Mesh -> False,
PlotStyle -> Opacity[1/5, ColorData["Legacy", "PowderBlue"]]]


And that's how to make a (virtual) wine glass. Filling it with (virtual) wine is left as an exercise for the interested reader.

-

A different (but related) possibility from my previous answer is to use B-splines for rendering both the outline and the corresponding surface of revolution. The strategy that follows is adapted from The NURBS Book by Piegl and Tiller.

First, we build the B-spline corresponding to the outline. Here again are the points for the wine glass contour:

wineGlassPoints = {{32.75, 283.75}, {37.75, 275.75}, {43.25, 267.25}, {49.25, 256.75},
{53.75, 247.75}, {58.25, 236.25}, {61.75, 224.25}, {64.25, 211.75},
{65.25, 198.75}, {64.25, 185.75}, {61.75, 174.25}, {58.25, 165.75},
{53.25, 157.25}, {46.25, 149.25}, {38.75, 143.75}, {30.25, 139.75},
{23.25, 137.25}, {17.75, 135.25}, {13.75, 134.75}, {13.75, 128.75},
{10.25, 126.75}, {7.25, 122.25}, {5.75, 115.75}, {4.75, 109.25},
{4.25, 101.25}, {3.75, 88.25}, {3.75, 70.25}, {3.75, 53.75},
{4.75, 39.75}, {6.75, 25.75}, {8.25, 20.25}, {11.25, 16.75},
{17.25, 14.25}, {25.25, 11.25}, {33.75, 9.25}, {41.25, 7.75},
{48.25, 5.75}, {48.75, 0.25}};


We will use, again, Eugene Lee's centripetal parametrization scheme to generate corresponding parameter values:

parametrizeCurve[pts_List, a : (_?NumericQ) : 1/2] :=
FoldList[Plus, 0, Normalize[(Norm /@ Differences[pts])^a, Total]] /; MatrixQ[pts, NumericQ]

tvals = parametrizeCurve[wineGlassPoints];


At this point, we do things a bit differently from my previous answer. One thing about B-splines is that they need the control points, as opposed to the interpolation points. We thus have to generate the control points for the (cubic) B-spline that passes through wineGlassPoints. To do that, we use a procedure suggested by Piegl and Tiller:

m = 3; (* degree of the B-spline *)
(* knots for interpolating B-spline *)
knots = Join[ConstantArray[0, m + 1], MovingAverage[ArrayPad[tvals, -1], m],
ConstantArray[1, m + 1]];
(* basis function matrix *)
bas = Table[BSplineBasis[{m, knots}, j - 1, tvals[[i]]],
{i, Length[wineGlassPoints]}, {j, Length[wineGlassPoints]}];
ctrlpts = LinearSolve[bas, wineGlassPoints];


With the control points now available, it's easy to render the B-spline corresponding to the wine glass profile. You could do something like ParametricPlot[BSplineFunction[ctrlpts, SplineDegree -> 3, SplineKnots -> knots][u] // Evaluate, {u, 0, 1}] to see the profile curve, but you should know that Mathematica has the BSplineCurve[] graphics primitive that you can use instead:

Graphics[{BSplineCurve[ctrlpts, SplineDegree -> 3, SplineKnots -> knots],
{AbsolutePointSize[4], Point[wineGlassPoints]}},
Axes -> None, Frame -> True]


For generating the corresponding surface of revolution, we could again use an appropriately constructed rotation matrix to generate the equations of the surface of revolution (just like in my previous answer). Instead, we use a NURBS representation of the circle for generating the surface of revolution. This is the approach taken in this article by Piegl and Tiller (also discussed in The NURBS Book).

First, recall the following B-spline representation of a unit circle:

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};
BSplineCurve[circPoints, SplineDegree -> 2, SplineKnots -> circKnots,
SplineWeights -> circWts] // Graphics


Here, then, is the procedure for generating the B-spline surface of revolution corresponding to the wine glass:

(* control points for surface of revolution *)
wgpts = Map[Function[pt, Append[#1 pt, #2]], circPoints] & @@@ ctrlpts;
(* weights for surface of revolution *)
wgwts = ConstantArray[circWts, Length[ctrlpts]];


We could use ParametricPlot3D[] along with BSplineFunction[] for rendering the surface of revolution, but we could instead use the BSplineSurface[] primitive. Let's generate the wine glass in the same style as my previous answer:

Graphics3D[{Directive[EdgeForm[], Opacity[1/5, ColorData["Legacy", "PowderBlue"]]],
BSplineSurface[wgpts, SplineClosed -> {False, True}, SplineDegree -> {3, 2},
SplineKnots -> {knots, circKnots}, SplineWeights -> wgwts]},
Boxed -> False, Lighting -> "Neutral"]


And we're done.

After writing up the code for this answer, some searching around shows that a similar strategy for generating surfaces of revolution has been done by Yu-Sung Chang in this Wolfram Demonstration.

-
Sadly you'll probably never read this comment but: this is wonderful! – Mr.Wizard Sep 30 '14 at 14:40
Thanks for the complement. :) – J. M. May 1 at 19:22
"Back... somewhat." Excellent! Welcome back! :-) – Mr.Wizard May 1 at 20:39

Here is some old code I used once, not highly efficient or anything, main idea is to negate the image. Go from top to bottom until a non-white pixel is found, that becomes the function value to interpolate, then I RevolutionPlot it, Image courtesy of drinkstuff.com. Works for most beer-glass images with white background you can find on google image search.

img = Import["http://www.drinkstuff.com/productimg/50540.jpg"];
img = ImageTake[img, All, {1, ImageDimensions[img][[2]]/2}]
img
{width, height} = ImageDimensions[img];
width = width - 3;
(* For given y, find first non-white pixel *)
f[y_] := Module[{i = 1},
While[ImageValue[img, {i, y}] > .99 && i < width, i = i + 1];
width - i + 1]
DiscretePlot[f[y], {y, 1, height}, AspectRatio -> Automatic]
g = Interpolation[
{#1, f[#1]} & /@ Range[1, height, 5]
, InterpolationOrder -> 7];
RevolutionPlot3D[g[y], {y, 1, height - 35}, BoxRatios -> 1,
RevolutionAxis -> {1, 0, 0}, AspectRatio -> 1,
PlotRange -> {{0, height}, {-height/2., height/2.}, {-height/2.,
height/2.}}, ColorFunction -> "GrayTones", Axes -> False,
Boxed -> False, PlotStyle -> Directive[Opacity[0.08]],
MeshStyle -> Directive[Opacity[0.07]], Background -> Transparent]


-