## ## Title: intersectplot ## Created: Thu Dec 9 16:26:07 1993 ## Authors: Paul Zimmermann and Sylvain Petitjean ## ## ## Copyright (C) 1993 Paul Zimmermann and Sylvain Petitjean ## ## This program is free software; you can redistribute it and/or modify it ## under the terms of the GNU General Public License as published by the ## Free Software Foundation; either version 2 of the License, or (at your ## option) any later version. ## ## This program is distributed in the hope that it will be useful, but WITHOUT ## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or ## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for ## more details. ## ## You should have received a copy of the GNU General Public License along ## with this program; see the file COPYING. If not, write to the Free ## Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA ## 02111-1307, USA. ## ## Description: plots the intersection of two 3D surfaces given by implicit ## equation. This implementation is based on the following ## papers. ## ## [1] Tracing surface intersections, by Bajaj, Hoffmann, Lynch and Hopcroft, ## Computer Aided Geometric Design 5 (1988), pages 285-307 ## ## [2] A New Curve Tracing Algorithm and Some Applications, by David J. ## Kriegman and Jean Ponce, in Curves and Surfaces, by Laurent, ## Le Mehaute', Schumaker (eds), pages 267-270, Academic Press, 1991 alias( Solve=`intersectplot/solve`, Eps=`intersectplot/epsilon`, Maxit=`intersectplot/maxit`, Newpt=`intersectplot/newpt` ); Eps:=Float(1,-8): # precision for Newton iteration Maxit:=20: # maximum number of Newton iterations # the steps (A) ... (D) are those of the article [2] intersectplot := proc(Eq1,Eq2:algebraic, vars:[name,name,name]) local x,y,z,J,extrema,p,l,n,q,mids,invJ,dx0,pp,dp,newpp,branch,dir,s,h,detJ, eq1,eq2,Deltak,features,brpoints,feat,colorinfo; option `Copyright Paul Zimmermann and Sylvain Petitjean, 1993`; description `intersection of two surfaces`; features:=NULL; colorinfo:=NULL; brpoints:=40; # First read the different features. The color is treated # separately. The option numpoints serves as the number of # points per branch. for feat in [args[4 .. nargs]] do if op(1,feat) = 'color' or op(1,feat) = 'colour' then colorinfo := readlib('`plot/color`')(op(2,feat)); elif op(1,feat) = 'numpoints' then brpoints := op(2,feat) else features := features,feat fi od; features := readlib('`plot3d/options3d`')(features); eq1:=Eq1; if type(eq1,polynom) then eq1:=convert([seq(op(1,x),x=sqrfree(eq1)[2])],`*`) fi; eq2:=Eq2; if type(eq2,polynom) then eq2:=convert([seq(op(1,x),x=sqrfree(eq2)[2])],`*`) fi; x:=op(1,vars); y:=op(2,vars); z:=op(3,vars); # step (A) : compute all extremal points in direction x # computing the Jacobian of [eq1,eq2] with respect to [y,z] J := linalg[jacobian]([eq1,eq2],[y,z]); detJ := linalg[det](J); if detJ=0 then # try with direction y l:=x; x:=y; y:=z; z:=l; J := linalg[jacobian]([eq1,eq2],[y,z]); # now [z,x] detJ := linalg[det](J); if detJ=0 then # try with direction z l:=x; x:=y; y:=z; z:=l; J := linalg[jacobian]([eq1,eq2],[y,z]); # now [x,y] if detJ=0 then RETURN(FAIL) fi; fi; fi; extrema := Solve({eq1,eq2,detJ},{x,y,z}); if extrema={} then RETURN(PLOT3D(CURVES())) fi; # empty intersection userinfo(2,plots,`extremal points`,op(extrema)); # step (B) : compute all intersections of G with the hyperplanes # orthogonal to the x axis at the extremal points l := {}; for p in extrema do l := l union {subs(p,x)}; od; l := sort(convert(l,list)); # step (C) : intersect with the midpoints of the intervals mids := NULL; for n to nops(l)-1 do userinfo(2,plots,`considering interval`,l[n]..l[n+1]); p := (l[n]+l[n+1])/2; s := Solve(subs(x=p,{eq1,eq2}),{y,z}); mids := mids,seq([l[n]..l[n+1],q union {x=p}],q=s) od; # try towards -infinity and +infinity if nops(l)=1 then h:=1 else h:=evalf((l[nops(l)]-l[1])/(nops(l)-1)) fi; p:=l[1]-h; s := traperror(Solve(subs(x=p,{eq1,eq2}),{y,z})); if s<>lasterror then mids:=mids,seq([p-h..p+h,q union {x=p}],q=s) fi; p:=l[nops(l)]+h; s := traperror(Solve(subs(x=p,{eq1,eq2}),{y,z})); if s<>lasterror then mids:=mids,seq([p-h..p+h,q union {x=p}],q=s) fi; userinfo(2,plots,`midpoints of real branches`,mids); # step (D) : march numerically from the intersections found in step (C) # to those found in step (B) by predicting new points through # Taylor expansions and correcting through Newton iterations invJ := linalg[inverse](J); # precompute Taylor approximation dp := linalg[multiply](invJ, linalg[scalarmul]([diff(eq1,x),diff(eq2,x)],-dx0)); dp := [x+dx0,y+dp[1],z+dp[2]]; # Newton iteration Deltak := linalg[multiply](invJ,linalg[scalarmul]([eq1,eq2],-1)); Deltak := subs(x=pp[1],y=pp[2],z=pp[3],[x,y+Deltak[1],z+Deltak[2]]); s := NULL; userinfo(2,plots,`number of branchs=`,nops([mids])); for q in [mids] do # compute branch with midpoint q userinfo(3,plots,`computing branch on`,q); branch := subs(op(2,q),[x,y,z]); h := (op(2,op(1,q))-op(1,op(1,q)))/brpoints; for dir in [-1,1] do # left and right p := subs(op(2,q),[x,y,z]); while abs(op(1,p)-op((3+dir)/2,op(1,q)))>2*h/3 do # first order Taylor expansion from p with step dx0 pp := eval(subs(x=p[1],y=p[2],z=p[3],dx0=dir*h,dp)); userinfo(3,plots,`Newton iteration, starting from`,pp); for n do # Newton iteration newpp := eval(Deltak); userinfo(4,plots,`Newton iteration, new point is`,newpp); if (op(2,newpp)-op(2,pp))^2+(op(3,newpp)-op(3,pp))^2Maxit then break fi; od; if n>Maxit then break fi; if dir=-1 then branch := newpp,branch else branch:=branch,newpp fi; p := pp; od; # while ... od; # for dir ... # add branch s := s,[branch]; od; PLOT3D(CURVES(s,colorinfo),features); end: Solve := proc(sys,unks) local res,s; res := solve(sys,unks); if res=NULL or has([res],fsolve) then res := traperror(fsolve(sys,unks)); if res=lasterror then res:=NULL fi else #if has([res],RootOf) then res := seq(allvalues(s,'d'),s={res}) fi; if has([res],RootOf) then res := seq(allvalues(s,'dependent'),s={res}) fi; res := op(evalf({res})); fi; # only consider real solutions select(proc(x) not has(x,I) end,{res}); end: rimplot := proc(eq1:algebraic, vars:[name,name,name], viewpt:[constant,constant,constant]) local x,y,z,x0,y0,z0,feat,eq2,features,r,phi,theta; option `Copyright Paul Zimmermann and Sylvain Petitjean, 1993`; description `rim or contour generator of a surface, i.e. the locus of points at which the line of sight grazes the surface`; x:=op(1,vars); y:=op(2,vars); z:=op(3,vars); x0:=op(1,viewpt); y0:=op(2,viewpt); z0:=op(3,viewpt); features := NULL; for feat in [args[4 .. nargs]] do features := features,feat od; r := sqrt(x0^2+y0^2+z0^2); if r = 0 then r = 3; x0=3; y0=0; z0=0; fi; phi := arccos(z0/r); if x0 = 0 then if sign(y0)*sign(sin(phi)) > 0 then theta := 90 else theta := 270 fi else theta := 180*arctan(y0/x0)/Pi fi; phi := 180*phi/Pi; features := features,orientation=[theta,phi]; eq2:= x0 * diff(eq1,x) + y0 * diff(eq1,y) + z0 * diff(eq1,z); intersectplot(eq1,eq2,vars,features); end: outlineplot := proc(eq1:algebraic, vars:[name,name,name], viewpt:[constant,constant,constant]) local x0,y0,z0,feat,features,r,phi,theta,rim,points,points2,ptlist,colorinfo; option `Copyright Paul Zimmermann and Sylvain Petitjean, 1993`; description `outline of a surface, i.e. projection onto some image plane of the locus of points at which the line of sight grazes the surface`; x0:=op(1,viewpt); y0:=op(2,viewpt); z0:=op(3,viewpt); features := NULL; colorinfo := NULL; for feat in [args[4 .. nargs]] do if op(1,feat) = 'color' or op(1,feat) = 'colour' then colorinfo := readlib('`plot/color`')(op(2,feat)) else features := features,feat fi od; r := sqrt(x0^2+y0^2+z0^2); if r = 0 then r = 3; x0=3; y0=0; z0=0; fi; phi := arccos(z0/r); if x0 = 0 then if sign(y0)*sign(sin(phi)) > 0 then theta := 90 else theta := 270 fi else theta := 180*arctan(y0/x0)/Pi fi; phi := 180*phi/Pi; rim := rimplot(eq1,vars,viewpt,features); features := NULL; for feat in [args[4 .. nargs]] do if op(1,feat) = 'numpoints' or op(1,feat) = 'color' or op(1,feat) = 'colour' then else features := features,feat fi od; features := features,orientation=[theta,phi]; features := readlib('`plot3d/options3d`')(features); points:=op(op(rim)[1]); points2 := NULL; for ptlist in points do points2 := points2,map(proc(u,x0,y0,z0) newpt(u,x0,y0,z0) end,ptlist,x0,y0,z0); od; PLOT3D(CURVES(points2,colorinfo),features); end: newpt := proc(pt,x0,y0,z0) local l,i; l := 1-(op(1,pt)*x0+op(2,pt)*y0+op(3,pt)*z0)/(x0^2+y0^2+z0^2); [seq(op(i,pt)+l*args[i+1],i=1..3)] end: `help/text/outlineplot` := TEXT( `FUNCTION: outlineplot - outline of a surface`, ``, `CALLING SEQUENCE:`, ` outlineplot(expr,[x,y,z],pt)`, ``, `PARAMETERS:`, ` expr - an expression in x, y and x`, ` pt - a point given as [,,]`, ``, `SYNOPSIS:`, `- The function outlineplot plots the outline of a three-dimensional`, ` implicitly defined surface observed under orthogonal projection`, ` from a given point.`, ``, `- Additional arguments are passed on to plot3d, apart from the argument`, ` numpoints=m which means that the number of points in each branch`, ` is reset to m.`, ``, `- Note that the surface is assumed to be transparent.`, ``, `EXAMPLES:`, `# outline of a squash-shaped object`, `outlineplot(4*y^4+3*x*y^2-5*y^2+4*z^2+6*x^2-2*x*y+2*x+3*y-1,[x,y,z],[3,0,0]);`, `# outline of a dimpled surface`, `outlineplot((4*x^2+3*y^2)^2-4*x^2-5*y^2+4*z^2-1,[x,y,z],[0,3,0],numpoints=80);`, `# outline of a banana-shaped object`, `outlineplot(23*x^4+x^2*y^2-37*x^2*y-2*x*y^2-15*x^2-2*x*y+16*y^2+16*z^2-x+16*y,[x,y,z],[0,0,3],color=blue,numpoints=100);` ): `help/text/rimplot` := TEXT( `FUNCTION: rimplot - rim of a surface`, ``, `CALLING SEQUENCE:`, ` rimplot(expr,[x,y,z],pt)`, ``, `PARAMETERS:`, ` expr - an expression in x, y and x`, ` pt - a point given as [,,]`, ``, `SYNOPSIS:`, `- The function rimplot plots the rim or contour generator of a`, ` three-dimensional implicitly defined surface observed under`, ` orthogonal projection from a given point, i.e. the locus of points at`, ` which the line of sight grazes the surface.`, ``, `- Additional arguments are passed on to plot3d, apart from the argument`, ` numpoints=m which means that the number of points in each branch`, ` is reset to m.`, ``, `- Note that the surface is assumed to be transparent.`, ``, `EXAMPLES:`, `# rim of a squash-shaped object`, `rimplot(4*y^4+3*x*y^2-5*y^2+4*z^2+6*x^2-2*x*y+2*x+3*y-1,[x,y,z],[3,0,0]);`, `# rim of a dimpled surface`, `rimplot((4*x^2+3*y^2)^2-4*x^2-5*y^2+4*z^2-1,[x,y,z],[0,3,0],numpoints=80);`, `# rim of a banana-shaped object`, `rimplot(23*x^4+x^2*y^2-37*x^2*y-2*x*y^2-15*x^2-2*x*y+16*y^2+16*z^2-x+16*y,[x,y,z],[0,0,3],color=blue,numpoints=100);` ): `help/text/intersectplot` := TEXT( `FUNCTION: intersectplot - intersection of two surfaces`, ``, `CALLING SEQUENCE:`, ` intersectplot(expr1,expr2,[x,y,z],)`, ``, `PARAMETERS:`, ` expr1,expr2 - two expressions in x, y and z`, ``, `SYNOPSIS:`, `- The function intersectplot plots the intersection of two three-dimensional`, ` implicitly defined surfaces. It uses an algorithm that first computes the`, ` extremal or singular points of the intersection in the x-direction, then`, ` constructs all branches between extremal or singular points using Taylor`, ` expansions and Newton iteration.`, ``, `- The algorithm used is described in the following papers:`, ` [1] Tracing surface intersections, by Bajaj, Hoffmann, Lynch and Hopcroft,`, ` Computer Aided Geometric Design 5 (1988), pages 285-307.`, ` [2] A New Curve Tracing Algorithm and Some Applications,`, ` by David J. Kriegman and Jean Ponce, in Curves and Surfaces, 1991,`, ` by Laurent, Le Mehaute', Schumaker (eds), pages 267-270, Academic Press.`, ``, `- The ranges in x, y and z of the output are automatically determined from`, ` the values of the extremal or singular points.`, ``, `- Two variables are used to control the algorithm and can be reset by the`, ` user: intersectplot/epsilon is the square of the precision needed`, ` in Newton iteration (Float(1,-8) by default) and intersectplot/maxit is`, ` the maximum number of Newton iterations allowed (20 by default).`, ``, `- Remaining arguments are interpreted as options which are specified as equa-`, ` tions of the form option = value. These options are the same as those found`, ` in plot3d. See also ?plot3d[options]. In particular numpoints=m means that`, ` the number of points in each branch is m (the default value is 40). Note`, ` that to avoid possible holes, this value has to be sufficiently high.`, ``, `- Details about the steps of the algorithm as described in [2] can be`, ` obtained by setting the variable infolevel[plots] to a higher value.`, ``, `EXAMPLES:`, `# cylinder-sphere intersection`, `intersectplot(x^2+y^2-1,x^2+y^2+z^2-2,[x,y,z]);`, `# sphere-hyperboloid intersection`, `intersectplot(x^2+y^2+z^2-5,x*y*z-1,[x,y,z]);`, `# nodal singularity`, `intersectplot(z+y^2-x^3,z+x^2,[x,y,z]);`, `# tacnode singularity`, `intersectplot(z+x^4+y^4,z+y^2,[x,y,z]);`, `# tacnode and nodal singularities`, `intersectplot(z-2*x^4-y^4,z-3*x^2*y+y^2-2*y^3,[x,y,z]);`, `# cyclinder-cylindrical gaussian intersection`, `intersectplot(exp(-x^2-y^2)-z,y^2+z^2-x^2,[x,y,z],numpoints=50);`, `# cylinder-cylinder intersection with both surfaces`, `e1 := x^2+z^2+2*z:`, `e2 := y^2+z^2+4*z:`, `with(plots):`, `c1 := implicitplot3d(e1,x=-2..2,y=-3..3,z=-4..0):`, `c2 := implicitplot3d(e2,x=-2..2,y=-3..3,z=-4..0):`, `d := intersectplot(e1,e2,[x,y,z]):`, `display({c1,c2,d},style=PATCHNOGRID);` ):