# Plotting implicitly-defined space curves

It is known that space curves can either be defined parametrically,

\begin{align*}x&=f(t)\\y&=g(t)\\z&=h(t)\end{align*}

or as the intersection of two surfaces,

\begin{align*}F(x,y,z)&=0\\G(x,y,z)&=0\end{align*}

Curves represented parametrically can of course be plotted in Mathematica using ParametricPlot3D[]. Though implicitly-defined plane curves can be plotted with ContourPlot[], and implicitly-defined surfaces can be plotted with ContourPlot3D[], no facilities exist for plotting space curves like the intersection of the torus $(x^2+y^2+z^2+8)^2=36(x^2+y^2)$ and the cylinder $y^2+(z-2)^2=4$:

Sometimes, one might be lucky and manage to find a parametrization for the intersection of two algebraic surfaces, but these situations are few and far between, especially if the two surfaces are of sufficiently high degree. The situation is worse if at least one of the surfaces is transcendental.

How might one write a routine that plots space curves defined as the intersection of two implicitly-defined surfaces?

It would be preferable if the routine returns only Line[] objects representing the space curve. A routine that handles only algebraic surfaces would be an acceptable answer, but it would be nice if your routine can handle transcendental surfaces as well. A bonus feature for the routine might be the ability to determine if the two surfaces given do not have a space curve intersection, or intersect only at a single point, or other such degeneracies.

-
Finally, I'll be able to properly intersect the columns with the roof –  rm -rf May 23 '12 at 20:18
@R.M, Funny, I thought getting the intersection of two parametrically-defined surfaces was an even tougher problem... so I decided not to ask for it in the question. :) That would be nice, though. –  Ｊ. Ｍ. May 23 '12 at 20:25
@J. M.♦ Can we extend the question a bit by considering intersection of two BSpline surfaces that has no explicit analytical form? This will be very effective for 3D modeling. –  PlatoManiac May 23 '12 at 21:28
@Plato: that might be best done as a separate question... (that would be covered by the "two parametrically-defined surfaces" case, as BSplineSurface[] objects can be transformed into BSplineFunction[]s). –  Ｊ. Ｍ. May 23 '12 at 21:36
@rcollyer Regarding the equation: if you accept a differential equation then it's easy: take the cross product of the gradients of the implicit functions, and perhaps normalize it (call it t[r], where r is a 3D point). Then the curve is described by D[r[u],u] == t[r[u]] for curve parameter u. This is because t must be tangent to both surfaces for all points r on the intersection curve. You "just" have to integrate that equation with a starting point r0 on the intersection. –  Jens May 24 '12 at 3:26

I take zero credit for this. It is a method I learned from Maxim Rytin.

ContourPlot3D[{(x^2 + y^2 + z^2 + 8)^2 - 36 (x^2 + y^2),
y^2 + (z - 2)^2 - 4}, {x, -4, 4}, {y, -4, 4}, {z, -2, 2},
Contours -> {0}, ContourStyle -> Opacity[0], Mesh -> None,
BoundaryStyle -> {1 -> None, 2 -> None, {1, 2} -> {{Green, Tube[.03]}}}, Boxed -> False]


-
+1 Magical ;-) {1,2}->Green means intersection (boundary) between surface 1 and 2 will be green. Here is a minimal set of options to make it work: ContourPlot3D[{(x^2 + y^2 + z^2 + 8)^2 - 36 (x^2 + y^2), y^2 + (z - 2)^2 - 4}, {x, -4, 4}, {y, -4, 4}, {z, -2, 2}, ContourStyle -> None, Mesh -> None, BoundaryStyle -> {1 -> None, 2 -> None, {1, 2} -> Green}] –  Vitaliy Kaurov May 23 '12 at 22:28
Nice! I find it interesting, though, that the docs do not say that the Filling format specification can be applied to BoundaryStyle. Very useful. –  rcollyer May 24 '12 at 3:01
I'll just note that if you take a look at the InputForm[] of the output of Daniel's version, the Polygon[] objects representing the isosurfaces are still there, just transparent. One could of course do something like ContourPlot3D[{(x^2 + y^2 + z^2 + 8)^2 - 36 (x^2 + y^2), y^2 + (z - 2)^2 - 4}, {x, -4, 4}, {y, -4, 4}, {z, -2, 2}, BoundaryStyle -> {1 -> None, 2 -> None, {1, 2} -> {}}, ContourStyle -> None, Mesh -> None] (as Vitaliy comments) if one just wants the curves themselves. So many unused points, though... –  Ｊ. Ｍ. May 24 '12 at 5:49
@R.M Thanks for the updates. I seem to be getting an awful lot of upvotes for code that was never of my own devising, and not notably better (that I can tell) than that in the other response. Embarrassed am I. –  Daniel Lichtblau May 24 '12 at 14:22
As rightly you should be keeping all the geometric goodies to yourself. There should be a Spelunking tag... –  Yves Klett Sep 17 '12 at 8:33

This is admittedly a hack, and does not work as well with PlotPoints / MaxRecursion as a dedicated function would (you will notice that I needed to increase PlotPoints for a good result), but it's a good way to make a quick and dirty plot:

ContourPlot3D[(x^2 + y^2 + z^2 + 8)^2 == 36 (x^2 + y^2),
{x, -4, 4}, {y, -4, 4}, {z, -2, 2},
MeshFunctions -> {Function[{x, y, z}, y^2 + (z - 2)^2 - 4]},
Mesh -> {{0}},
ContourStyle -> None,
PlotPoints -> 30, BoxRatios -> Automatic]


Per @whuber's comment below, the same thing can be achieved using RegionFunction:

ContourPlot3D[(x^2 + y^2 + z^2 + 8)^2 == 36 (x^2 + y^2),
{x, -4, 4}, {y, -4, 4}, {z, -2, 2},
RegionFunction -> Function[{x, y, z}, y^2 + (z - 2)^2 < 4],
ContourStyle -> None, Mesh -> None, BoxRatios -> Automatic]


The style can be adjusted using BoundaryStyle.

I believe that this will give good quality results even without increasing PlotPoints, but I'll need to do some more testing before I can be certain about that...

-
What made this problem hard back in the days of old Mathematica was that the version of ContourPlot3D[] available in the standard packages didn't have the MeshFunctions option, and thus it took a bit of labor to look at which polygons were intersecting. This is nice, even if it is quick-and-dirty. :) –  Ｊ. Ｍ. May 23 '12 at 20:29
+1 With suitable graphics settings, the RegionFunction option will accomplish something similar. –  whuber May 23 '12 at 20:32
@J.M. I think Mathematica's plotting framework is really powerful, and it uses adaptive sampling at every stage. I am not sure if mesh lines are computed with a higher resolution than the surface or not, but the system definitely has the capability to do this for the intersection of two surfaces: it does it when computing exclusions. Try e.g. Plot3D[Boole[x^2 + y^2 < 1], {x, -1.1, 1.1}, {y, -1.1, 1.1}] and not how smooth and sharp the circle boundary is. If you turn off exclusions (Exclusions -> None), it'll be a lot more jagged. –  Szabolcs May 23 '12 at 20:33
@whuber And you just answered the question that didn't fit in the comment I wrote above! RegionFunction will sample the intersection curve very smoothly, similarly to how Exclusions are computed! (The Plot example above.) This should give a better result than MeshFunctions, I think, but I'd need to test it. –  Szabolcs May 23 '12 at 20:34