# How to make a smooth revolution surface plot

I'm trying to plot a 3d revolution plot from a set of 2d points. These data points form a 2d curve, then we rotate that curve around y axis and get a 3d surface. @J. M. has a well explained and very helpful post at here which deals exactly the problem I have. However, I tried to use the method, and get a 3d surface that is very rough and not smooth.

Here is the 2d data points:

and this is how it looks like

Graphics[Line[points], Frame -> True]

and then I use J. M.'s method(code copied and modified from here) (it will take about 10 seconds to run)

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

tvals = parametrizeCurve[points];

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[points]}, {j, Length[points]}];
ctrlpts = LinearSolve[bas, points];

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

and then generate the 3d surface

Graphics3D[{Directive[EdgeForm[]],
BSplineSurface[wgpts, SplineClosed -> {False, True},
SplineDegree -> {3, 2}, SplineKnots -> {knots, circKnots},
SplineWeights -> wgwts]}, Boxed -> False]

We can see that the surface is not smooth. It looks like the surface is composed by flat rings. So how can we make the surface smooth?

Edit:

I think this unsmoothness may come from my data rather than the rotation process, so I tried to smooth my data using something like

pointsSmooth=ExponentialMovingAverage[points, 1/20];

then I get the a smoother surface, but ExponentialMovingAverage seems to have removed the end points and there is a hole on the surface, which I don't want.

Also smoothing using a smooth constant like 1/20 largely modified the original data:

Graphics[{Red, Line[ExponentialMovingAverage[points, 1/20]], Blue,
Line[points]}, Frame -> True, AspectRatio -> 1]

So is it possible to smooth the data while keep the general shape so that it will give a better smooth surface? Or there are other ways to contract a smooth surface from the data?

• Now that you found the source of the lack of smoothness lies in your data, it is not clear what your question is now. – m_goldberg May 29 '14 at 4:21
• @m_goldberg Still, my goal is get a smooth surface :) – xslittlegrass May 29 '14 at 4:37
• Why is it so important to avoid the small deviation between the actual data and approximate surface? – Jens May 29 '14 at 16:57
• @Jens What did you refer to deviation, the roughness (discontinuity in the derivative) at the surface or the change of the curve after ExponentialMovingAverage? 1) If it's the later. I'm trying to make a movie of the wave function evolution by plotting the isosurface, when an atom is experiencing a laser field. since it's a physical simulation, I would like to make it looks good but also keep the data as unmodified as possible. 2) About the roughness, I generate the data and then rendered in povray. This roughness is more visible in the rendered image and is very annoying. – xslittlegrass May 29 '14 at 17:17

Sounds like what you actually need after your edit is a way to smooth a list of data while keeping the endpoints fixed. Here's a dumb approach that will work with any "symmetrical" smoothing filter, including GaussianFilter, MeanFilter, even MedianFilter. It won't work with ExponentialMovingAverage, though, because that's not symmetrical, although it should if you average the results from ExponentialMovingAverage and Reverse@ExponentialMovingAverage@Reverse.

smooth[list_, filter_] :=
Take[filter[Join[
(2 First@list - #) & /@ Reverse@Rest@list,
list,
(2 Last@list - #) & /@ Reverse@Most@list]],
{Length@list, 2 Length@list - 1}]

All it does is it extends the data in a "flipped" form about both endpoints -- for example, $[1,2,10]$ will become $[\color{grey}{-8,0},1,2,10,\color{grey}{18,19}]$ -- then smooths that, and drops the extra entries.

For example:

smooth[{a, b, c, d, e}, GaussianFilter[#, 1] &] // Simplify
{a,
(b BesselI[0, 1/4] + (a + c) BesselI[1, 1/4])/(BesselI[0, 1/4] + 2 BesselI[1, 1/4]),
(c BesselI[0, 1/4] + (b + d) BesselI[1, 1/4])/(BesselI[0, 1/4] + 2 BesselI[1, 1/4]),
(d BesselI[0, 1/4] + (c + e) BesselI[1, 1/4])/(BesselI[0, 1/4] + 2 BesselI[1, 1/4]),
e}

:)

Turns out GaussianSmooth smooths across all dimensions of the array by default, so the $x$-coordinate gets averaged with the $y$-coordinate and vice versa. Oops.

{xs, ys} = smooth[#, GaussianFilter[#, 5] &] & /@ Transpose[points];
{fx, fy} = Interpolation /@ {xs, ys};
RevolutionPlot3D[{fx[i], fy[i]}, {i, 1, Length[points]}]

You can set smoothedPoints = Transpose@{xs, ys} if you want to use the whole parametrizeCurve stuff instead.

• I don't know where I messed up but I can't make it work using your method. The smoothed curve seems dramatically changes. Could you give more details and show a plot? – xslittlegrass May 29 '14 at 19:02
• @xslittlegrass: Fixed. – Rahul May 30 '14 at 3:11
• Still, downsampling the data through a spline fit, as in your answer, is probably a more sensible approach since you don't really need to retain all those hundreds of degrees of freedom. – Rahul May 30 '14 at 4:13
• Thanks a lot for your answer! Your method is simpler and works very well. – xslittlegrass May 30 '14 at 4:39
• It looks more expedient to use ArrayPad[]: smooth[list_, filter_] := With[{n = Length[list]}, Take[filter[ArrayPad[list, n - 1, "Extrapolated", InterpolationOrder -> 1]], {n, 2 n - 1}]] – J. M. is away Mar 12 '18 at 12:57

The the roughness on the surface is due to the noise in the original data, so that the first derivative is not smooth.

This shows the first derivative of the x and y components using original data, we can see it's very noisy

tvals = parametrizeCurve[points];
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[points]}, {j, Length[points]}];
ctrlpts = LinearSolve[bas, points];

df1[u_] =
Evaluate@D[
BSplineFunction[ctrlpts, SplineDegree -> 3, SplineKnots -> knots][
u], u];
Plot[{Evaluate[df1[u]], Evaluate[df1[u]][[2]]}, {u, 0, 1}]

We need to remove that noise in order to get a smooth surface. This can be achieved using linear fit. Andy Ross presented very effective a spline model fit at here. I modified his code and create this function to apply to my 2d data.

dataSmoothSplineMethod[ls_?(ArrayQ[#, 1] &), knotNum_Integer] :=
Module[{SplineModel, lth = Length[ls], knots, data},
data = Transpose[{Range[1, lth], ls}];
SplineModel[data_, deg_, knots_] :=
Block[{basis, allKnots, n, kmin, kmax}, n = Length[knots] + 2;
kmin = 0;
kmax = Ceiling[Max[data[[All, 1]]]] + 1;
basis =
Array[\[FormalX]^# &, deg + 1, 0]~Join~
Table[BSplineBasis[{deg, knots}, i, \[FormalX]], {i, 0,
Length[knots] - deg - 2}];
LinearModelFit[data, basis, \[FormalX]]];
knots = {1}~Join~Range[2, lth - 1, Round[lth/(knotNum - 2)]]~
Join~{lth};
SplineModel[data, 3, knots]
]

dataSmoothSplineMethod[ls_?(ArrayQ[#, 2] &), knotNum_Integer] :=
Module[{SplineModel, lth = Length[ls], lsx, lsy},
lsx = ls[[All, 1]];
lsy = ls[[All, 2]];
{dataSmoothSplineMethod[lsx, knotNum],
dataSmoothSplineMethod[lsy, knotNum]}
]

f2 = dataSmoothSplineMethod[points, 20];

and here shows the comparison of the smoothed data and original data, we can the changes are almost not visible

Legended[Show[
ParametricPlot[Through@f2[x], {x, 1, Length[points]},
PlotPoints -> 100, PlotStyle -> Red],
Graphics[{Blue, Line[points]}]],
Placed[LineLegend[{Red, Blue}, {"spline fit", "original data"}],
ImageScaled@{0.5, 0.5}]]

but the first derivative is much more smooth

df2[u_] = {D[f2[[1]][u], u], D[f2[[2]][u], u]};

Plot[Evaluate[df2[u]], {u, 1, Length[points]}]

Using this smoothed data, we can create a much more smooth 3d surface

pointsSmoothed = Table[Through@f2[x], {x, 1, Length[points]}];
parametrizeCurve[pts_List, a : (_?NumericQ) : 1/2] :=
FoldList[Plus, 0, Normalize[(Norm /@ Differences[pts])^a, Total]] /;
MatrixQ[pts, NumericQ]
tvals = parametrizeCurve[pointsSmoothed];
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[pointsSmoothed]}, {j, Length[pointsSmoothed]}];
ctrlpts = LinearSolve[bas, pointsSmoothed];
circpointsSmoothed = {{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]], circpointsSmoothed] & @@@
ctrlpts;
wgwts = ConstantArray[circWts, Length[ctrlpts]];

Graphics3D[{Directive[EdgeForm[]],
BSplineSurface[wgpts, SplineClosed -> {False, True},
SplineDegree -> {3, 2}, SplineKnots -> {knots, circKnots},
SplineWeights -> wgwts]}, Boxed -> False]

You have a lot of points, so generating two interpolation functions from them, one for the bottom of the surface and the other for the top, should give a smooth surface of revolution. If the default plot isn't smooth enough, you can always increase PlotPoints.

max = Max[First /@ points];
ii = Position[points, {max, _}][[1, 1]];
btm = Interpolation[points[[;; ii]]];
top = Interpolation[points[[ii ;;]]];
RevolutionPlot3D[{{t, btm[t]}, {t, top[t]}}, {t, 0., max}]

• Thanks, but it doesn't seems to work well even I increased the PlotPoints to 100. Here is the result I get imgur.com/t61hEW8 – xslittlegrass May 29 '14 at 15:24