(* :Title: ImplicitPlot3D *)

(* :Author: Steven Wilkinson, Northern Kentucky University *)

(* :Summary:  Plot3D[ ] plots a function of two variables.
Many simple surfaces (e.g., ellipsoids, hyperboloids etc.)
are not the graphs of functions of two variables.
ImplicitPlot3D[ ] allows the user to plot surfaces implicitly
defined by an equation of three variables. *)

(* :Context: Graphics`ImplicitPlot3D` *)

(* :Package Version: 2.0 *)

(* :History:
        V2.0 by Steven Wilkinson, July 1992
*)

(* :Keywords: solution set, 3D graphics,
        implicitly defined surfaces *)

(* :Warning: The algorithm used by ImplicitPlot3D is
essentially cubic with respect to the desired resolution.
Therefore, it can be time and memory intensive. *)

(* :Mathematica Version: 2.0 *)

(* :Limitations:  The package works by taking the box given by
the range of the three variables and dividing it into a number
of subboxes as specified by the PlotPoints option. Using the
intermediate value theorem, the package checks each of the
subboxes to see if the surface passes through any of its edges.
                The algorithm is actually a little faster than cubic
because of the following modification. The user can have the
package make a number of prepasses. The first pass divides the
original box coarsely. The package retains only those subboxes
in which it can detect the surface. These subboxes are then
repeatedly subdivided in the same manner retaining only those
new subboxes containing the surface.
                The limitations of this package (which are
especially emphasized by the modification) is that if the
surface intersects an edge of a subbox an even number of times,
the package may not be able to discern that the surface has
entered that subbox. The coarser the subdivisions, the more
likely this is to happen.*)

BeginPackage[
        "Graphics`ImplicitPlot3D`","Utilities`FilterOptions`"
]

ImplicitPlot3D::usage =
        "ImplicitPlot3D[eqn,{x,a,b},{y,c,d},{z,e,f}] draws a graph
        of the set of points with coordinates (x,y,z) in the box
        a2x2b, c2y2d, e2z2f that satisfy the equation eqn.
        ImplicitPlot3D[func,{x,a,b},{y,c,d},{z,e,f}] draws a graph
        of the set of points with coordinates (x,y,z) in the box
        a2x2b, c2y2d, e2z2f that satisfy the equation func = 0.";

Passes::usage =
        "Passes is an option for ImplicitPlot3D. It is used to
        specify how many times to subdivide the plot range in
        order to reach the stated value of PlotPoints. This option
        is ignored if the PlotPoints option as a list of lists."

Options[ImplicitPlot3D] =
        {Compiled -> True,
        PlotPoints -> {{5,5,5},{3,3,3}},
        Passes -> Automatic};

Begin["`Private`"]

        vector = {{1,0,0},{0,1,0},{0,0,1},{0,0,0}};

        vertex =
                Flatten[Table[{i,j,k},{i,0,1},{j,0,1},{k,0,1}],2];

        edges =
                Select[
                        Flatten[
                                Table[{vertex[[i]],vector[[j]]},{i,8},{j,4}]
                        ,1],
                        (Max[#[[1]]+#[[2]]]<2)&
                ];

        prepassEdges =
                        Select[
                                Flatten[
                                        Table[{vertex[[i]],vector[[j]]},{i,7},{j,3}]
                                        ,1],
                                (Max[#[[1]]+#[[2]]]<2)&
                        ];


ImplicitPlot3D[
        lhs_ == rhs_, xr:{_,_,_},yr:{_,_,_},zr:{_,_,_},opts___
] :=
    ImplicitPlot3D[lhs - rhs,xr,yr,zr,opts]

ImplicitPlot3D[
        func_,{x_,xMin_,xMax_},{y_,yMin_,yMax_},{z_,zMin_,zMax_},
        opts___
] :=
        Module[{cubes,cubePointer,f,boxes},

                {cubes,cubePointer,divisions} =
                        ReleaseHold[
                                Hold[
                                        makeCubeList[PlotPoints,Passes]
                                ]/. {opts} /. Options[ImplicitPlot3D]
                        ];

                f =
                        If[Compiled /. {opts} /. Options[ImplicitPlot3D],
                                Compile[{x,y,z},func],
                                Function[{x,y,z},func]
                        ];

                boxes = {{{xMin,xMax},{yMin,yMax},{zMin,zMax}}};

                If[Length[cubePointer] > 1,
                        boxes =
                                Fold[
                                        prepassBoxes[
                                                f,#1,#2,cubes,divisions,cubePointer
                                        ]&,
                                        boxes,
                                        Range[Length[divisions] - 1]
                                ]
                ];

                boxes =
                                Flatten[
                                        finalBoxes[
                                                f,#,divisions[[-1]],
                                                cubes[[ cubePointer[[-1]] ]]
                                        ]& /@ boxes
                                ,1];

                Show[
                        Graphics3D[boxes],FilterOptions[Graphics3D,opts]
                ]
        ]



(* Functions makeCubeList, standardizePlotPoints, checkErr,
factorInteger2, adjustLength, and reduceList are used to take
the options PlotPoints and Passes and put them in a standard
form. *)

makeCubeList[plotPoints_,passes_]:=
        Module[{plotPts,boxPtr,plotPtsReduced},
                plotPts=standardizePlotPoints[plotPoints,passes];
                If[plotPts===$Failed,Return[$Failed] ];
                {plotPtsReduced,boxPtr} = reduceList[plotPts];
                {       Apply[
                                        Flatten[Array[{#1-1,#2-1,#3-1}&,{#1,#2,#3}],2]&,
                                        plotPtsReduced,
                                        1
                        ],
                        boxPtr,
                        plotPts
                }
        ]

standardizePlotPoints[plotPts_,pss_] :=
        checkErr /@ plotPts /; Depth[plotPts] == 3

standardizePlotPoints[plotPts_,pss_] :=
        standardizePlotPoints[{plotPts},pss] /; Depth[plotPts] == 1

standardizePlotPoints[plotPts_,pss_] :=
        Module[{},
                Message[ImplicitPlot3D::ppta,plotPoints];
                $Failed
        ]  /; Depth[plotPts] != 2

ImplicitPlot3D::ppta = "`1` is not a list of numbers or a list of lists of numbers.";

standardizePlotPoints[plotPts_,passes_] :=
        Module[{plotPtsFactored,pss},
                plotPtsFactored = factorInteger2 /@ checkErr[plotPts];
                pss =
                        Which[
                                passes===Automatic,Max[Length /@ plotPtsFactored
   ],
                                passes > 0, Ceiling[passes],
                                True,1
                        ];
                Reverse[
                        Select[
                                Transpose[
                                        adjustLength[#,pss]& /@ plotPtsFactored
                                ],
                                # != {1,1,1}&
                        ]
                ]
        ]

checkErr[{p_Integer}] := {p,p,p} /; p > 0

checkErr[{p_Integer,q_Integer,r_Integer}] :=
        {p,q,r} /; p > 0 && q > 0 && r > 0

checkErr[p___] :=
        Message[ImplicitPlot3D::pptb,p]

ImplicitPlot3D::pptb = "`1` is not a list of either 1 or 3 positive integers.";


factorInteger2[n_] :=
        Flatten[
                Table[#[[1]],{#[[2]]}]& /@ FactorInteger[n]
        ]


adjustLength[factoredList_,n_] :=
        factoredList /; Length[factoredList] == n

adjustLength[factoredList_,n_] :=
        Join[   Table[1,{n-Length[factoredList]} ],
                        factoredList ] /;       Length[factoredList] < n

adjustLength[factoredList_,n_] :=
        Nest[
                Sort[
                        Rest[
                                ReplacePart[#,#[[1]]*#[[2]],2]
                        ]
                ]&,
                factoredList,
                Length[factoredList] - n
        ]

reduceList[x_List] :=
        Module[{reducedx},
                reducedx = Union[x];
                {reducedx, Flatten[Position[reducedx,#]& /@ x]}
        ]



(* Functions subtractAltered, sameFaceQ, orderVertices,
and reorderedVertices are used to find the planar convex hull
of the given points. *)

subtractAltered[x_,y_] := x-y /; IntegerQ[x]&&IntegerQ[y]

subtractAltered[x_,y_] := z

SetAttributes[subtractAltered,Listable];

sameFaceQ[{v1_,e1_},{v2_,e2_}]:=
        Module[{z},
                Apply[Times,subtractAltered[v1+z*e1,v2+z*e2] ] == 0
        ]

orderVertices[{ordered_,unordered_}] :=
        {
                Join[ordered,unordered],{}
        } /; Length[unordered] < 2

orderVertices[{ordered_,unordered_}] :=
        {
                Join[ordered,{unordered[[1]]} ], Rest[unordered]
        } /; sameFaceQ[ordered[[-1,2]],unordered[[1,2]] ]

orderVertices[{ordered_,unordered_}] :=
        {ordered,RotateLeft[unordered]}

reorderedVertices[vertices_] :=
        {} /; Length[vertices] < 2

reorderedVertices[vertices_] :=
        Transpose[vertices][[1]] /; Length[vertices] < 4

reorderedVertices[vertices_] :=
        Transpose[
                First[
                        FixedPoint[
                                orderVertices,
                                {{vertices[[1]]},Rest[vertices]},
                                Length[vertices](Length[vertices]-1)/2
                        ]
                ]
        ][[1]]



(* Functions locusQ, prepassCubeQ, prepassLocusQ, prepassBoxes,
ppBox, finalBoxes are used to determine if and where the
surface intersects a given box. *)


locusQ[h_,{i1_,j1_,k1_},{{i2_,j2_,k2_},{0,0,0}}] :=
        Module[
                {value =
                        Abs[h[i1+i2,j1+j2,k1+k2]]},
                        NumberQ[value] && value == 0.
                ]

locusQ[h_,{i1_,j1_,k1_},{{i2_,j2_,k2_},{i3_,j3_,k3_}}] :=
        Module[
                {value =
                        h[i1+i2,j1+j2,k1+k2] h[i1+i2+i3,j1+j2+j3,k1+k2+k3]
                },
                NumberQ[value] && value < 0.
        ]


prepassCubeQ[g__] :=
        Apply[Or,Map[prepassLocusQ[g,#]&,prepassEdges] ]

prepassLocusQ[h_,{i_,j_,k_},{{i1_,j1_,k1_},{i2_,j2_,k2_}}] :=
        Module[
                {value =
                        h[i+i1,j+j1,k+k1]*h[i+i1+i2,j+j1+j2,k+k1+k2]
                },
                NumberQ[value] && value <= 0.
        ]

prepassBoxes[func_,boxes_,n_,cb_,divis_,ptr_] :=
        Flatten[
                Map[
                        ppBox[func,#,n,cb,divis,ptr]&,
                        boxes
                ],
                1
        ]

ppBox[
        func_,{{u1_,u2_},{v1_,v2_},{w1_,w2_}},n_,cb_,divis_,ptr_
] :=
        Module[{f,du,dv,dw,nu,nv,nw,cubes},
                {nu,nv,nw} = divis[[n]];
                cubes = cb[[ptr[[n]] ]];
                du=(u2-u1)/nu;
                dv=(v2-v1)/nv;
                dw=(w2-w1)/nw;
                f[i_,j_,k_] := f[i,j,k] =
                        func[i*du+u1,j*dv+v1,k*dw+w1];
                Return[
                                Apply[{{#1*du+u1,(1+#1)*du+u1},
                                                {#2*dv+v1,(1+#2)*dv+v1},
                                                {#3*dw+w1,(1+#3)*dw+w1}}&,
                                                Select[cubes,prepassCubeQ[f,#]&],{1}
                                ]
                ]
        ]

finalBoxes[
        func_,{{u1_,u2_},{v1_,v2_},{w1_,w2_}},{nu_,nv_,nw_},cb_
] :=
        Module[{f,du,dv,dw,coord,plate},
                du=(u2-u1)/nu;
                dv=(v2-v1)/nv;
                dw=(w2-w1)/nw;
                f[i_,j_,k_] := f[i,j,k] =
                        func[i*du+u1,j*dv+v1,k*dw+w1];
                coord[{i1_,j1_,k1_},{{i2_,j2_,k2_},{0,0,0}}] :=
                        {u1+(i1+i2)*du,v1+(j1+j2)*dv,w1+(k1+k2)*dw};
                coord[{i1_,j1_,k1_},{{i2_,j2_,k2_},{i3_,j3_,k3_}}] :=
                        Module[{f1=f[i1+i2,j1+j2,k1+k2],f2},
                                f2=f1/(f[i1+i2+i3,j1+j2+j3,k1+k2+k3]-f1);
                                {u1 + du(i1+i2-i3*f2),
                                        v1 + dv(j1+j2-j3*f2),
                                        w1 + dw(k1+k2-k3*f2)}
                        ];

                plate[i_,j_,k_] :=
                        Map[{coord[{i,j,k},#],#}&,
                                Select[edges,locusQ[f,{i,j,k},#]&]
                        ];

                Return[
                        Polygon /@
                                        reorderedVertices /@
                                                Select[Apply[plate,cb,{1}], (# != {})&]
                        ]
        ]

 End[] (* "`Private`" *)

EndPackage[]  (* "Graphics`ImplicitPlot3D`" *)
