43
$\begingroup$

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:

tube

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?

$\endgroup$
  • 4
    $\begingroup$ Your picture is 3D - so do you mean a set of {x,y,z} coordinates? $\endgroup$ – Vitaliy Kaurov Mar 15 '12 at 20:37
  • $\begingroup$ My data is 2D, I just did {x, y, 0} so that it would be in the x-y plane. $\endgroup$ – s0rce Mar 16 '12 at 18:22
45
$\begingroup$

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]

Mathematica graphics

$\endgroup$
  • $\begingroup$ Thanks very much, I eventually found the Tube command but I had no idea you could replace the line in the parametric plot with a tube so easily! You get an upvote, you would also get one from my girlfriend who's work I was doing but she doesn't have an account. $\endgroup$ – s0rce Mar 16 '12 at 18:22
31
$\begingroup$

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]

B-spline tube with B-spline cross section

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

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

B-spline tube with circular cross section

$\endgroup$
25
$\begingroup$

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

tube through random points

$\endgroup$
  • $\begingroup$ @BrettChampion Thanks for the code edit ;-) $\endgroup$ – Vitaliy Kaurov Mar 16 '12 at 0:16
20
$\begingroup$

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

Mathematica graphics

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

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

Mathematica graphics

With BSplineCurve you would get

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

Mathematica graphics

$\endgroup$
19
$\begingroup$

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
]

Trefoil knot extruded with a deltoid

$\endgroup$
17
$\begingroup$

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]

tube

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]
$\endgroup$
4
$\begingroup$

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

perturbed clothoid

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

tube with gaps

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

Frenet-Serret tube

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

tube, Bishop-style

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]

square knot with threefold-symmetric cross section

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.