Extruding along a path

Question

I'm trying to render a 3D image of a path by extruding a circular cross-section along the path, to create a "snake-like" path.

Here is an image I found to illustrate:

I can't seem to figure out if there is a way to do this?

I just found the Tube command; now I need to find a way to turn a set of points into a curve. Basically, I need to find the necessary control points for a BezierCurve or a BSplineCurve to fit a set of {x, y} coordinates.

Can the Interpolation with the Spline method be used to generate a spline from a set of points?

Answer 1

In Mathematica 7 or 8, you can just use Tube. Please see the docs for many, many examples.

Example:

Show[ParametricPlot3D[{Cos[x], Sin[x], x/5}, {x, 0, 15}] /. 
  Line -> (Tube[#, 0.2] &), PlotRange -> All]

Answer 2

This question has been nicely addressed by the previous answers, so I'll just write about a method you can use if you want your tube to have custom cross-sections; or, like in this question, you need to have the tube as a bunch of Polygon[]s.

(* Pixar method; http://jcgt.org/published/0006/01/01/ *)
orthogonalDirections[{p1_?VectorQ, p2_?VectorQ}] := Module[{s, w, w1, xx, yy, zz},
     {xx, yy, zz} = Normalize[p2 - p1];
     s = 2 UnitStep[zz] - 1; w = -1/(s + zz); w1 = xx yy w;
     {{1 + s w xx^2, s w1, -s xx}, {w1, s + w yy^2, -yy}}]

orthogonalDirections[{p1_?VectorQ, p2_?VectorQ, p3_?VectorQ}] := Module[{d, u, v},
     {u, v} = Normalize /@ {p3 - p2, p1 - p2};
     If[Chop[Norm[u - v] Norm[u + v]] != 0,
        d = (u + v)/2; Normalize /@ {d, Cross[u, d]},
        orthogonalDirections[{p1, p2}]]]

extend[cs_, q_, d_, nrms_] :=
cs + Outer[Times, First[LinearSolve[Transpose[Prepend[-nrms, d]],
                                    q - Transpose[cs]]], d]

(* for custom cross-sections *)
crossSection[pointList_?MatrixQ, r_, csList_?MatrixQ] := Module[{p1, p2},
     {p1, p2} = Take[pointList, 2];
     (p1 + #) & /@ (r csList.orthogonalDirections[{p1, p2}])] /; 
Last[Dimensions[pointList]] == 3 && Last[Dimensions[csList]] == 2

(* for circular cross-sections *)
crossSection[pointList_?MatrixQ, r_, n_Integer] := 
 crossSection[pointList, r, Composition[Through, {Cos, Sin}] /@ Range[0, 2 Pi, 2 Pi/n]]

(* approximate vertex normals, for a smooth appearance *)
vertNormals[vl_ /; ArrayQ[vl, 3, NumericQ]] := Block[{mdu, mdv, msh},
    msh = ArrayPad[#, {{1, 1}, {1, 1}}, "Extrapolated", InterpolationOrder -> 2] & /@
          Transpose[vl, {2, 3, 1}];
    mdu = ListCorrelate[{{1, 0, -1}}/2, #, {{-2, 1}, {2, -1}}, 0] & /@ msh;
    mdv = ListCorrelate[{{-1}, {0}, {1}}/2, #, {{1, -2}, {-1, 2}}, 0] & /@ msh;
    MapThread[Composition[Normalize, Cross], Transpose[{mdu, mdv}, {1, 4, 2, 3}], 2]]

MakePolygons[vl_ /; ArrayQ[vl, 3, NumericQ], OptionsPattern[{"Normals" -> True}]] :=
Module[{dims = Most[Dimensions[vl]]}, 
       GraphicsComplex[Apply[Join, vl], 
                       Polygon[Flatten[Apply[Join[Reverse[#1], #2] &, 
                       Partition[Partition[Range[Times @@ dims], Last[dims]],
                                 {2, 2}, {1, 1}], {2}], 1]], 
                       If[TrueQ[OptionValue["Normals"] /. Automatic -> True], 
                          VertexNormals -> Apply[Join, vertNormals[vl]], 
                          Unevaluated[]]]]

Options[TubePolygons] = {"Normals" -> True, "Scale" -> 1.};
TubePolygons[path_?MatrixQ, cs : (_Integer | _?MatrixQ), OptionsPattern[]] :=
With[{p3 = PadRight[path, {Automatic, 3}]},
     MakePolygons[FoldList[Function[{p, t}, 
                                    extend[p, t[[2]], t[[2]] - t[[1]],
                                           orthogonalDirections[t]]], 
                           crossSection[p3, OptionValue["Scale"], cs], 
                           Partition[p3, 3, 1, {1, 2}, {}]], 
                  "Normals" -> OptionValue["Normals"]]]

Try it out:

path = First @ Cases[ParametricPlot3D[BSplineFunction[
       {{0, 0, 0}, {1, 1, 1}, {2, -1, -1}, {3, 0, 1}, {4, 1, -1}}][u] // Evaluate,
     {u, 0, 1}, MaxRecursion -> 1], Line[l_] :> l, Infinity];

cs = First @ Cases[ParametricPlot[
     BSplineFunction[{{0., 0.}, {0.25, 0.}, {0.5, 0.125}, {0.25, 0.25}, {0., 0.25}},
                     SplineClosed -> True][u] // Evaluate,
                 {u, 0, 1}, MaxRecursion -> 1], Line[l_] :> l, Infinity];

Graphics3D[{EdgeForm[], TubePolygons[path, cs]}, Boxed -> False]

Of course, you can elect to have a circular cross section, as is usual:

Graphics3D[{EdgeForm[], TubePolygons[path, 20, "Scale" -> .2]}, Boxed -> False]

Answer 3

I see you mentioned splines in your question. Your picture is 3D, but you used 2D {x,y} coordinates in the question. This little example uses a random set of control points and emphasizes the 3D nature of the Tube and {x,y,z} coordinates:

points = RandomReal[1, {20, 3}]; 
Export["tube.gif",
   Table[
      Graphics3D[
         {Orange, Specularity[White, 100], Tube[BSplineCurve[points], .03]},
         Boxed -> False, SphericalRegion -> True,
         ViewAngle -> .25, 
         ViewPoint -> RotationTransform[a, {0, 0, 1}][{3, 0, 3}]], 
      {a, 0, 2 Pi, .1}
   ]]

Answer 4

The disadvantage of using BSplineCurve is that generally speaking, the curve doesn't go through the control points. If you want the curve to be both smooth and go through the test points, you could Interpolation with Method -> "Spline" instead. For example

pts = RandomReal[1, {30, 3}];
interp = Interpolation[MapIndexed[{#2[[1]], #1} &, pts], Method -> "Spline"]

You could then use Szabolcs's method to plot interp and replace the line with a tube:

pl = ParametricPlot3D[interp[x], {x, 1, Length[pts]}, 
  PlotPoints -> 2 Length[pts],
  PlotStyle -> Green, PlotRange -> All] /. {Line[a_] :> Tube[a, .05]}

To show that the control points actually lie on the curve:

Show[pl, Graphics3D[{White, Sphere[#, .08] & /@ pts}]]

With BSplineCurve you would get

Graphics3D[{{Darker[Blue], Tube[BSplineCurve[pts], .05]}, {White, 
   Sphere[#, .08] & /@ pts}}]

Answer 5

Not directly related to question, but rather in addition to J.M.'s answer, here's a simple method to extrude any cross section (not just a circle) when you have a differentiable path:

(* First get the curve for the path, in this case the trefoil knot. *)
path = KnotData["Trefoil", "SpaceCurve"]

(* Now define the shape of the cross section, in this case a deltoid. *)
crosssection = {2 Cos[#] + Cos[2 #], 2 Sin[#] - Sin[2 #], 0} &

(* Plot the extruded path. First rotate the cross section to be 
   perpendicular to the tangent of the path, and then translate. *)
ParametricPlot3D[
  path[u] + 1/8 RotationMatrix[{{0, 0, 1}, path'[u]}].crosssection[v],
  {u, 0, 2 Pi}, {v, 0, 2 Pi},
  Boxed -> False, Axes -> False, Mesh -> False, MaxRecursion -> 3
]

Answer 6

Tube is the right thing to use here. Increase n to smoothly interpolate your curve, and r to control the thickness of your circle.

n = 280;
r = .15;
data = Table[{Cos[u] Sin[u^.3], Sin[u] Cos[u^.2], u/10}, {u, 0, 45,45/(n - 1)}];
Graphics3D[{CapForm[None], 
Tube[BSplineCurve[data], Table[r, {i, Length@data}]]},
Boxed -> False, Axes -> False]

If you want a function:

Extrude[{x_, y_, z_}, {t_, Start_, End_}, Discretization_, Radius_] := 
Module[{data3d},
       data3d =Table[{x, y, z}, {t, Start,End, (End - Start)/(Discretization - 1)}];
       Graphics3D[{CapForm[None],Tube[BSplineCurve[data3d],
       Table[Radius, {i, Length@data3d}]]}, Boxed -> False, Axes -> False]
      ];

The following will generate the same extrusion as the graphics:

Extrude[{Cos[u] Sin[u^.3], Sin[u] Cos[u^.2], u/10}, {u, 0, 45}, 280, .1]

Answer 7

I am writing this answer as a separate supplement to my previous one, and also as an expansion of Teake's answer.

Teake's answer features one possible solution to the problem of constructing a tube from a given path and cross section. It works for a lot of space curves, but can fail in some cases.

For instance, consider this perturbed clothoid:

path[t_] := With[{b = 1/2},
                 {FresnelS[t/Sqrt[1 + b^2]], FresnelC[t/Sqrt[1 + b^2]], b/Sqrt[1 + b^2] t}]

and say we want to make a tube with a circular cross section out of it:

circ[u_] := {Cos[u], Sin[u]}/15;

Let's try Teake's solution (where I had rearranged it into an equivalent form):

ParametricPlot3D[path[u] + circ[v].Most[RotationMatrix[{path'[u], {0, 0, 1}}]] // Evaluate,
                 {u, -5, 5}, {v, 0, 2 π}, Mesh -> False, PlotPoints -> {75, 45}]

Oh no! All those breaks!

One could try adding Exclusions -> None, but the tube still looks somewhat pinched at those places.

One could also try the more conventional method of using the Frenet-Serret frame to make the tube:

nbv[t_] = Rest[Last[FrenetSerretSystem[path[t], t]]];

ParametricPlot3D[path[t] + circ[u].nbv[t] // Evaluate, {t, -5, 5}, {u, 0, 2 π}, 
                 Mesh -> False, PlotPoints -> {75, 45}]

but that also produces a tube with a kink.

What to do?

The solution is to use the Bishop frame to set up the required local coordinate system around the space curve. Unlike the other approaches, the vectors produced by the Bishop frame are guaranteed not to collapse, which would lead to the pinches and gaps in the previously produced tubes.

Adapting Henrik's code from this answer to this situation, we have the following routine:

plotTube[f_, cs_, t0_: 0, {t_, tra__}, {u_, ura__},
         opts : OptionsPattern[{NDSolve, ParametricPlot3D}]] := 
    Module[{bf, bi, b0, nf, no, n0, s, ta, ts, w, w1, xx, yy, zz, T, γ, κ},
           γ = Function[t, f];
           (* unit tangent vector *)
           T = Function[t, Evaluate[Normalize[γ'[t]]]];
           (* curvature vector *)
           κ = Function[t, Evaluate[(γ''[t] - (γ''[t].T[t]) T[t])/Norm[γ'[t]]^2]];
           (* Pixar method; http://jcgt.org/published/0006/01/01/ *)
           ts = T[t0]; {xx, yy, zz} = ts;
           s = 2 UnitStep[zz] - 1; w = -1/(s + zz); w1 = xx yy w;
           n0 = {1 + s w xx^2, s w1, -s xx}; b0 = {w1, s + w yy^2, -yy};
           (* Bishop frame *)
           {nf, bf} = NDSolveValue[{ta'[t] == Norm[γ'[t]]
                                    {no[t].κ[t], bi[t].κ[t]}.{no[t], bi[t]},
                                    no'[t] == -Norm[γ'[t]] ta[t] no[t].κ[t],
                                    bi'[t] == -Norm[γ'[t]] ta[t] bi[t].κ[t],
                                    ta[t0] == ts, no[t0] == n0, bi[t0] == b0},
                                   {no, bi}, {t, tra},
                                   Method -> {"OrthogonalProjection",
                                              Dimensions -> {3, 3}},
                                   FilterRules[{opts}, Options[NDSolveValue]]];
           ParametricPlot3D[γ[t] + cs.{nf[t], bf[t]}, {t, tra}, {u, ura},
                            Method -> Automatic, 
                            Evaluate[FilterRules[{opts}, Options[ParametricPlot3D]]], 
                            Mesh -> False]]

where I use orthogonal projection to enforce the orthogonality of the Bishop frame.

This gives a pretty nice tube:

plotTube[path[t], circ[u], {t, -5, 5}, {u, 0, 2 π}, PlotPoints -> {75, 45}]

A more elaborate example:

plotTube[{3 Cos[t]/4 - Cos[3 t]/2, Sin[t]/4 + Sin[3 t]/2, Sin[5 t]/4},
         (2 + Cos[3 u])/3 {Cos[u], Sin[u]}/30, {t, 0, 2 π}, {u, 0, 2 π},
         PlotPoints -> 55]