MathGroup Archive 2009

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: Re: Slicing a surface

  • To: mathgroup at
  • Subject: [mg97661] Re: [mg97421] Re: Slicing a surface
  • From: Daniel Lichtblau <danl at>
  • Date: Wed, 18 Mar 2009 04:53:27 -0500 (EST)
  • References: <>

SotonJames wrote:
> Thanks, but I think the real problem is that I had to solve the equations numerically, so I can't solve x,y,z in terms of s and t or anything like that, and I can't use Reduce. The three functions that make up the surface are interpolating functions and I'd just like to be able to plot the intersection of the surface and the z=constant plane and  plot it in a flattened x-y plane. It originally seemed to me that this should be quite straightforward, but it's surprising how tricky it's turned out to be.
> Thanks again for the help, no worries if you can't see how to solve it.
> James

If there is a straightforward way to do this, I don't see it. I will 
describe a couple of approaches that take a bit of work. I do not use 
interpolating functions in the example below but I believe the same 
ideas will work (both methods, though the first will require reasonable 
smoothness of said functions).

Neither method is dressed up with error detection, and the second one is 
in no way modularized. So the code is just an indication of how one 
might proceed.

Here is the example I use.

exprs = {Cos[s] + Sin[t] + s*t, s - Sin[s*t], -Cos[s]*Sin[t] + t};

ParametricPlot3D[exprs, {s, -4, 4}, {t, -4, 4}]

We'll look at the slice for z = 2.9.

Method 1: Homotopy tracking.

We first find a point on the slice, then set up and solve an ODE system 
to track the level curve. This has a drawback: if there are disconnected 
curves, you'd need to find an initial point on each component.

To find one point I use NMinimize because it is fairly robust and not in 
need of good initial guesses. I minimize the absolute value of the z 
coordinate function minus the slice value.

   z0_Real, {s_, smin_, smax_}, {t_, tmin_, tmax_}] := {s, t} /.
   Last[NMinimize[{Abs[zfunc - z0], smin <= s <= smax,
      tmin <= t <= tmax}, {s, t}]]

The level curve function will find one point using pt[] above, then 
solve for the surface parameters {s,t} in terms of a new parameter u. 
One equation is given by constancy of z on the slice curve. In order not 
to be underdetermined we need to impose a further condition on {s',t'}. 
We treat this u a bit like "arc length" parameter, in giving it an 
absolute value condition to satisy (sum of s and t components must be 
1). This can cause trouble at singularities (perhaps at crossings, say; 
I'm not actually sure). I don't use sum-of-squares because that seems to 
cause trouble, probably because it is not longer a linear system in the 

levelcurve[{xfun_, yfun_, zfun_}, zconst_,
   srange : {s_, smin_, smax_}, trange : {t_, tmin_, tmax_}] :=
   {s0, t0, pred, corr, u, endu, subst, diffeqns, paramfuncs},
   {s0, t0} = pt[zfun, zconst, srange, trange];
   subst = {s -> s[u], t -> t[u]};
   diffeqns = {(D[zfun, s] /. subst)*s'[u] + (D[zfun, t] /. subst)*
        t'[u] == 0, s'[u] + t'[u] == 1, s[0] == s0, t[0] == t0};
   paramfuncs =
    First[NDSolve[diffeqns, {s[u], t[u]}, {u, 0, Infinity}]];
   endu = paramfuncs[[1, 2, 0, 1, 1, 2]];
   ParametricPlot[({xfun, yfun} /. subst) /. paramfuncs, {u, 0, endu}]

Notice we let u run from 0 to infinity. There will be a value at which 
it ends, and we grab it from the resulting interpolating functions from 
the NDSolve solution. I suspect there is a better way to do this using 
Method->{"EvenLocator",opts} but I do not know how to go about it myself.

Anyway, this shows the level curve (slice).

levelcurve[exprs, 2.9, {s, -4, 4}, {t, -4, 4}]

Method 2: Find where the surface graphics intersects the slice plane.

First we take the polygons from the Graphics3D created by ParametricPlot.

pp = ParametricPlot3D[exprs, {s, -4, 4}, {t, -4, 4}]

This gives the raw points.

pts = First[Cases[pp, GraphicsComplex[a_, b___] :> a, Infinity]];

Now get the indices of polygons, and use them to get the vertex lists 
for the polygons.

polygons = Flatten[Cases[pp, Polygon[a_, b___] :> a, Infinity], 1];
vertexpolygons = Map[pts[[#]] &, polygons];

We now make polygons into triangles. This is the same as is done 
implicitly in plotting, except for 4-gons (we treat them the same as 
n-gons for n>4, instead of triangulating along a diagonal).

triangulate[{x_, y_, z_}] := {x, y, z}
triangulate[ll_List /; Length[ll] > 3] :=
  With[{midp = Total[ll]/Length[ll]},
   Map[Append[#, midp] &, Partition[ll, 2, 1, -1]]]

vertextriangles = Flatten[Map[triangulate, vertexpolygons], 1];

Now we need to find only the triangles that interesct our slice. The 
idea: take the ones for which there is at least one vertex z coordinate 
greater than, and at least one vertex z coordinate less than, the slice 

triangleIntersectsZ[{{_, _, z1_}, {_, _, z2_}, {_, _, z3_}}, z_] :=
  Not[z1 > z && z2 > z && z3 > z] && Not[z1 < z && z2 < z && z3 < z]
intersectingTriangles[triangles_, val_] :=
  Select[triangles, triangleIntersectsZ[#, val] &]

triangles29 = intersectingTriangles[vertextriangles, 2.9];

For each triangle there are two segments that cross the slice, unless 
one lies on it. That's a special case I really do not carefully address 
below. There is a bit of code but it is entirely untested and probably 
will require tweaking to work.

sliceIntersectionPointPairs[{p1 : {_, _, z1_}, p2 : {_, _, z2_},
    p3 : {_, _, z3_}}, z_] :=
  Select[{{p1, p2}, {p2, p3}, {p3,
     p1}}, (Last[#[[1]]] - z)*(Last[#[[2]]] - z) <= 0 &]

pointpairs29 = Map[sliceIntersectionPointPairs[#, 2.9] &, triangles29];

With the segment endpoints, we can now find where the segments hit the 
slicing plane.

sliceIntersectionSegment[{p1 : {_, _, z_}, p2 : {_, _, z_}}, z_] :=
  Map[Most, {p1, p2}]
sliceIntersectionSegment[{p1 : {_, _, z1_}, p2 : {_, _, z2_}}, z_] :=
  Map[Most, {p1 + (z1 - z)/(z1 - z2)*(p2 - p1)}]

segpts29 =
   Map[sliceIntersectionSegment[#, 2.9] &, pointpairs29, {2}];

Now make lines out of endpoints of triangles intersecting the slice plane.

segments29 =
   Apply[Line, segpts29, {1}] /. Line[{a_}, {b_}] :> Line[{a, b}];

We can now plot this. It seems to be plausible.


Further remarks:

Visually, the levelcurve code seems to do a bit betteneed work to handle 
slices with disconnected curve components, and might have trouble with 
curves that have singularities. i suspect the triangles intersection 
approach can be made fairly robust by someone more adept at programming 
with Graphics3d objects.

Daniel Lichtblau
Wolfram Research

  • Prev by Date: Re: Problem with mathlink:PLEASE HELP ME
  • Next by Date: Re: formatted table output to ascii file
  • Previous by thread: Re: Slicing a surface
  • Next by thread: Re: Slicing a surface