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: Problem with mathlink:PLEASE HELP ME**

**Re: formatted table output to ascii file**

**Re: Slicing a surface**

**Re: Slicing a surface**