Re: Re: Slicing a surface
- To: mathgroup at smc.vnet.net
- Subject: [mg97661] Re: [mg97421] Re: Slicing a surface
- From: Daniel Lichtblau <danl at wolfram.com>
- Date: Wed, 18 Mar 2009 04:53:27 -0500 (EST)
- References: <200903130947.EAA20062@smc.vnet.net>
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.
pt[zfunc_,
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
derivatives.
levelcurve[{xfun_, yfun_, zfun_}, zconst_,
srange : {s_, smin_, smax_}, trange : {t_, tmin_, tmax_}] :=
Module[
{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
value.
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.
Graphics[segments29]
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
- References:
- Re: Slicing a surface
- From: SotonJames <james.french@soton.ac.uk>
- Re: Slicing a surface