Drawing curves on surfaces

By Gianluca Gorni, University of Udine, Italy. Last revised: November 22, 2001.

Description: drawing a parametric curve on a parametric surface, so that the curve does not "submarines" under the surface here and there. This is achieved by adapting the curve to the triangulation of the surface.

Instructions: load the package by clicking on the button below, then evaluate the examples.

Button[
      StyleForm["Click here to load package!", 
        FontColor -> RGBColor[1, 0, 0], FontSize -> 14], 
      ButtonFunction :> 
        FrontEndExecute[FrontEndToken["EvaluateInitialization"]], 
      ButtonEvaluator -> None, Active -> True]

Mathematica package (to load into the kernel)

preamble

\(\(BeginPackage["CurveOnSurface`", "Utilities`FilterOptions`"];\)\)

Unprotect["`*"];
ClearAll["`*"];

\(\(CurveOnSurface::usage = 
"CurveOnSurface[surface,curve1,curve2,...,options] draws the curves on the surface adding plot points to the curve as to avoid submarineing of the curves under the surface. The surface can be argument for either Plot3D or ParametricPlot3D. The curves are expressed in the parameter space of the surface. Example: CurveOnSurface[{x^2+y^2,{x,-1,1},{y,-1,1}},{{Cos[t],Sin[t]},{t,0,2Pi},PlotPoints->41}].";\)\)

SetAttributes[CurveOnSurface, HoldAll];

Begin["`Private`"];

insertion of points along a line

The function toSplitQ[1,pt1,pt2] decides whether the line between the two points pt1 and pt2 crosses over a vertical line with integer abscissa, and gives a list of the integer abscissas that are crossed. SplitQ[2,pt1,pt2] does the same with crossing with horizontal lines.

toSplitQ[n_Integer /; n < 3, pt0_List /; NumericQ[pt0[[1]]], 
            pt1_List] /; 
          MatrixQ[{pt0, 
              pt1}] := 
Module[{coord = Transpose[{pt0, 
                    pt1}][[n]], med1, 
            med2}, 
{med1, 
              med2} = {Floor[Min[#]] + 1, Ceiling[Max[#]] - 1} &[
              coord]; 
{med1 ≤ 
              med2, If[OrderedQ[coord], Identity, Reverse][
              Range[med1, med2]]}];

The function makeSplit[1][{pt1,pt2}] inserts new points along the line segment between pt1 and pt2 at all crossings with vertical lines of integer abscissa. So does makeSplit[2][{pt1,pt2}] with horizontal lines.

makeSplit[
            n_Integer /; n < 3][{pt0_List /; NumericQ[pt0[[1]]], 
            pt1_List}] := 
Module[{}, 
{answer, divisions} = 
            toSplitQ[n, pt0, pt1]; 
If[
            answer, 
Join[{pt0}, 
              t*pt1 + (1 - t)*pt0 /. 
                Map[Solve[(t*pt1 + (1 - t)*
                                  pt0)[[n]] == #, 
                        t][[1]] &, divisions], {pt1}], {pt0, 
              pt1}]];

The function toSplitQ[3,pt0,pt1] decides whether the line segment between the two points pt0 and pt1 crosses the diagonal of the integer cell the two points are in, and gives the coordinates of the crossing point.

toSplitQ[3, pt0_List, pt1_List] := 
      Module[{v, c}, 
v = 
          Floor[Map[Min, 
              Transpose[{pt0, pt1}]]]; 
c = {pt0 - v, 
            pt1 - v}; 
{Subtract @@ c[[1]]*
              Subtract @@ c[[2]] < 0, c}]

The function makeSplit[3][{pt0,pt1}] adds to the line between the points pt0 and pt1 the crossing (if any) between the line and the diagonal of the integer cell the two points are in.

makeSplit[3][{pt0_List /; NumericQ[pt0[[1]]], pt1_List}] := 
        Module[{answer, c0, c1}, 
{answer, {c0, c1}} = 
            toSplitQ[3, pt0, pt1]; 
If[answer, 
            Join[{pt0}, 
              t*pt1 + (1 - t) pt0 /. 
                Solve[Equal @@ (t*c1 + (1 - t) c0), t], {pt1}], {pt0, 
              pt1}]];

The function divideLine[ptsList] takes the list of points ptsList and inserts points at all crossings with vertical or horizontal or diagonal lines.

divideLine[ptsList_List] := 
        Fold[Flatten[Map[makeSplit[#2], Partition[#1, 2, 1]], 1] &, 
            ptsList, {1, 2, 3}] //. {a___, b_List, b_, 
              c___} :> {a, b, c};

making the 3D curves

The function weight[pt] expresses the point pt as a convex combination of the four vertices of the integer cell the point is in, and gives the coefficients.

weight[pt_List /; NumericQ[pt[[1]]]] := 
        Module[{vert, loc, belowDiagQ}, 
vert = 
            Floor[pt]; 
loc = 
            pt - vert; 
belowDiagQ = (Subtract @@ 
                  loc ≥ 0); 
Map[Flatten, 
            Cases[{{vert, w00}, {vert + {1, 0}, w10}, {vert + {1, 1}, 
                      w11}, {vert + {0, 1}, w01}} /. Solve[
                      Thread[loc == 
                            w00 {0, 0} + 
                              If[belowDiagQ, w10 {1, 0}, w01 {0, 1}] + 
                              w11 {1, 1}]~
                        Join~{w00 + If[belowDiagQ, w10, w01] + w11 == 
                            1}, {w00, If[belowDiagQ, w10, w01], 
                        w11}][[1]] /. {w10 -> 0, 
                  w01 -> 0}, {_, w_ /; w =!= 0}]]];

The function polygonVerticesCoordinates makes a matrix of the coordinates of the vertices of the surface.

SetAttributes[polygonVerticesCoordinates, HoldAll];

polygonVerticesCoordinates[
          f_ /; Not[Head[f] === List], {u_, uMin_, uMax_}, {v_, vMin_, 
            vMax_}, optsS___?OptionQ] := 
        Transpose[{Plot3D[u, {u, uMin, uMax}, {v, vMin, vMax}, 
                DisplayFunction -> Identity, optsS][[1]], Plot3D[
                v, {u, uMin, uMax}, {v, vMin, vMax}, 
                DisplayFunction -> Identity, optsS][[1]], Plot3D[
                f, {u, uMin, uMax}, {v, vMin, vMax}, 
                DisplayFunction -> Identity, optsS][[1]]}, {3, 1, 
            2}];

polygonVerticesCoordinates[{fx_, fy_, fz_}, {u_, uMin_, uMax_}, {v_, 
            vMin_, vMax_}, optsS___?OptionQ] := 
        Transpose[{Plot3D[fx, {u, uMin, uMax}, {v, vMin, vMax}, 
                DisplayFunction -> Identity, optsS][[1]], Plot3D[
                fy, {u, uMin, uMax}, {v, vMin, vMax}, 
                DisplayFunction -> Identity, optsS][[1]], Plot3D[
                fz, {u, uMin, uMax}, {v, vMin, vMax}, 
                DisplayFunction -> Identity, optsS][[1]]}, {3, 1, 
            2}];

The function normalizeCoordinates[pt,xMin,passoX,yMin,passoY] takes the grid with vertex at {xMin,yMin} and takes it to the integer grid starting at the origin.

normalizeCoordinates[pt_List /; NumericQ[pt[[1]]], xMin_, 
          passoX_, yMin_, passoY_] := 
        Rationalize[(pt - {xMin, yMin})/{passoX, 
              passoY},  .00001];

The function originalCoordinates[nrmPt,xMin,passoX,yMin,passoY] takes back a point in the integer grid to the original grid.

originalCoordinates[nrmPt_List /; NumericQ[nrmPt[[1]]], xMin_, 
          passoX_, yMin_, passoY_] := 
        nrmPt*{passoX, passoY} + {xMin, yMin};

The function basicCurvePoints[{gu,gv},{t,tMin,tMax}] effectively calculates the values of {gu,gv} for t from tMin to tMax at uniform increments of t. It uses ParametricPlot3D, and not Table or such, because we want the evaluation of our function to behave like the built-in Plot3D and ParametricPlot3D.

SetAttributes[basicCurvePoints, HoldAll];

basicCurvePoints[{gu_, gv_}, {t_, tMin_, tMax_}, optsC___?OptionQ] := 
        Transpose[
            Drop[Transpose[ParametricPlot3D[{gu, gv, 0}, {t, tMin, tMax}, 
                      optsC, PlotPoints -> 15, 
                      DisplayFunction -> Identity][[1]] /. 
                  Line[{a__}] :> 
                    Sequence @@ {a}], -1]] //. {x___, y_List, y_, 
              z___} :> {x, y, z};

The function internalCurveOnSurface coordinates a parametric curve with the grid of a surface. There is a version for parametric surfaces and one for cartesian surfaces.

SetAttributes[internalCurveOnSurface, HoldAll];

internalCurveOnSurface[{f_, {u_, uMin_, uMax_}, {v_, vMin_, vMax_}, 
            optsS___?OptionQ}, {{gu_, gv_, sty2___}, {t_, tMin_, tMax_}, 
            optsC___?OptionQ}] := 
        Module[{xMin, xMax, yMin, yMax, ptNo, n1, n2, passoX, passoY, 
            gridInfo, cur0, normalizedCur0, curveXMin, curveXMax, curveYMin, 
            curveYMax, subdividedCur, weightsList, verticesData, 
            weightedCoord, 
            cur3D}, 
{xMin, xMax, yMin, yMax} = 
            N[Evaluate[{uMin, uMax, vMin, 
                  vMax}]]; 
ptNo = PlotPoints /. 
                Evaluate[{optsS}] /. 
              PlotPoints -> {15, 15}; 
{n1, n2} = 
            If[NumericQ[ptNo], {ptNo, ptNo}, 
              ptNo]; 
{passoX, 
              passoY} = ({xMax, yMax} - {xMin, yMin})/({n1, n2} - 
                  1); 
gridInfo = {xMin, passoX, yMin, 
              passoY}; 
cur0 = 
            basicCurvePoints[{gu, gv}, {t, tMin, tMax}, 
              optsC]; 
normalizedCur0 = 
            Map[normalizeCoordinates[#, Sequence @@ gridInfo] &, 
              cur0]; 
{{curveXMin, curveXMax}, {curveYMin, 
                curveYMax}} = 
            Map[{Min[Floor[#]], Max[Ceiling[#]]} &, 
              Transpose[
                normalizedCur0~
                  Join~{{xMin, yMin}, {xMax, 
                      yMax}}]]; 
subdividedCur = 
            divideLine[normalizedCur0]; 
weightsList = 
            Map[weight, subdividedCur]; 
verticesData = 
            polygonVerticesCoordinates[
              f, {u, uMin + curveXMin*passoX, uMin + curveXMax*passoX}, {v, 
                vMin + curveYMin*passoY, vMin + curveYMax*passoY}, 
              PlotPoints -> {curveXMax - curveXMin + 1, 
                  curveYMax - curveYMin + 
                    1}]; 
weightedCoord = 
            Apply[#3*
                  verticesData[[#2 - curveYMin + 1, #1 - curveXMin + 
                        1]] &, 
              weightsList, {2}]; 
cur3D = 
            Line[Apply[Plus, weightedCoord, 
                1]]; 
Graphics3D[{sty2, 
              cur3D}]];

Plotting the surface separately:

SetAttributes[internalSurface, HoldAll];

internalSurface[{f_ /; Not[Head[f] === List], {u_, uMin_, 
              uMax_}, {v_, vMin_, vMax_}, optsS___?OptionQ}] := 
        Graphics3D[
            Plot3D[f, {u, uMin, uMax}, {v, vMin, vMax}, optsS, 
              DisplayFunction -> Identity, PlotPoints -> {15, 15}, 
              Mesh -> False]] /. 
          Polygon[pts_List] :> 
            Polygon[RotateLeft[pts]];

internalSurface[{{fx_, fy_, fz_, sty1_:  EdgeForm[]}, {u_, uMin_, 
              uMax_}, {v_, vMin_, vMax_}, optsS___?OptionQ}] := 
        ParametricPlot3D[{fx, fy, fz, sty1}, {u, uMin, uMax}, {v, vMin, 
            vMax}, optsS, DisplayFunction -> Identity, 
          PlotPoints -> {15, 15}];

Combining the surface and the curve(s):

CurveOnSurface[
        surf : {f_, {u_, uMin_, uMax_}, {v_, vMin_, vMax_}, 
            optsS___?OptionQ}, 
        curves : ({{_, _, ___}, {_Symbol, _, _}, ___?OptionQ} .. ), 
        optsGen___?OptionQ] := 
      Show[{internalSurface[surf], 
          Sequence @@ Map[internalCurveOnSurface[surf, #] &, {curves}]}, 
        optsGen, DisplayFunction -> $DisplayFunction]

End[];

Protect["`*"];

EndPackage[];

Examples

The naïve way, with the curve that submarines behind the surface:

Show[{ParametricPlot3D[{4  Cos[u] Cos[v], 3  Sin[u] Cos[v], 
              2  Sin[v], EdgeForm[]}, {u, 0, 2  Pi}, {v, -Pi/2, Pi/2}, 
            PlotPoints -> {21, 11}, DisplayFunction -> Identity], 
          ParametricPlot3D[
            Evaluate[{4  Cos[u] Cos[v], 3  Sin[u] Cos[v], 
                  2  Sin[v]} /. {u -> t*Cos[30  t] - Pi/2, 
                  v -> t*Sin[30  t] + Pi/6}], {t, 0, 1}, 
            PlotPoints -> 50, DisplayFunction -> Identity]}, 
        DisplayFunction -> $DisplayFunction, 
        Ticks -> None];

With CurveOnSurface the curve shows all the way (it takes a while to compute):

CurveOnSurface[{{4  Cos[u] Cos[v], 3  Sin[u] Cos[v], 2  Sin[v]}, {u, 0, 
                2  Pi}, {v, -Pi/2, Pi/2}, PlotPoints -> {21, 11}},{{t*Cos[30  t] - Pi/2, t*Sin[30  t] + Pi/6}, {t, 0, 1}, 
              PlotPoints -> 50},Ticks -> None];

The syntax is:

CurveOnSurface[
{content of what you would write into Plot3D or ParametricPlot3D to get the surface only},{content of what you would write into ParametricPlot to describe the curve in the 2 dimensional parameter space of the surface},
optional Graphics3D options]

You can do it the Plot3D style:

f = 
        1/(1 + (x + 1)^2 + y^2) - 
          1/(1 + (x - 1)^2 + y^2);

CurveOnSurface[{f, {x, -2, 2}, {y, -2, 2}},{{t*Cos[10  t], t*Sin[10  t]}, {t, 0, 2}, 
              PlotPoints -> 50}];

You can also draw more than one curve on the same surface:

CurveOnSurface[{{(3 + Cos[u]) Cos[v], (3 + Cos[u]) Sin[v], 
                Sin[u]}, {u, 0, 2  Pi}, {v, 0, 2  Pi}, 
              PlotPoints -> {15, 41}},Evaluate[Sequence@@Table[{{Cos[θ] t + Pi/4, 
                        Sin[θ] t - Pi/4}, {t, -4, 4}, 
                      PlotPoints -> 5},{θ, 0, Pi, Pi/6}]],Boxed -> False,Axes -> False];