(*---------------------------start of LiveGraphics3D.m-------------------*)


(*:Mathematica Version: 4.0 *)

(*:Package Version: 1.15 *)

(*:Title: LiveGraphics3D Toolbox *)

(* :Context: LiveGraphics3D` *)

(* :Author: Martin Kraus (mk) (Martin.Kraus@informatik.uni-stuttgart.de)
   some routines modified by Eric W. Weisstein (eww) (eww@wolfram.com)
*)

(* :Summary:  
   This package implements several functions for three-dimensional
   graphics, which are especially useful in preparing graphics
   for the LiveGraphics3D applet.

Master version available from:
  http://mathworld.wolfram.com/notebooks/LiveGraphics3D.m

*)


(*:History:
  v1.00 (1998-11-08): Written by Martin Kraus (mk)
  v1.01 (1998-11-11): Dir added (eww)
  v1.02 (1998-11-12): Minor mods to WriteLiveAnimationForm and reduced
                      number of sides for circle and disk (eww)
  v1.03 (1998-11-15): improved NonIntersectingPolygons,
                      reprogrammed DropInnerPolygons, 
                      removed old version of DropInnerPolygons and
                      PlotCentersOfMasses (mk)
  v1.04 (1998-11-16): corrected ClipAtPlotRange (did not work in v1.03) (mk)
  v1.05 (1998-11-16): documentation updated slightly,
                      tweaked circle, disk (eww)
  v1.06 (1998-11-25): new function UnifyPolygons and helpers 
                      (MergePolygons, PolygonInSamePlane) (mk)
  v1.07 (1998-11-28): NonIntersectingLines (mk)
  v1.08 (1998-12-05): documentation updated slightly (eww)
  v1.09 (1998-12-11): PolyhedronPolygon (eww+mk)
  v1.10 (1999-01-08): LinelessFaces (eww+mk)
  v1.11 (1999-01-17): OverlappingEdgesQ now more tolerant (mk)
  v1.12 (1999-03-30): Improved DropInnerPolygons to make the front face
                      of the starting polygon facing away from the origin. (mk)
  v1.13 (1999-07-17): Added BoundingBox[{}]={} to stop occasional errors in
                      working with duals of uniform polyhedra.  At some point,
                      I should figure out why this is happening and
                      actually fix it. (eww)
  v1.14 (1999-08-18): Added workaround for writing only 5 significant
                      numbers to files since NumberForm[InputForm[N[g]],5]
                      does not work properly under Mathematica 4.0.
  v1.15 (2000-01-20): Updated emails and URL (eww)
  
  (c) 1998-2000 Martin Kraus
*)

(*:Keywords:
  	Graphics3D, LiveGraphics3D
*)

(*:Sources:
   None. 
*)

(* :Limitations: 
	Several functions work only with flat, convex polygons.
	    However, PolygonTriangulation routines can be
	    used to break up concave polygons.
	Several functions deal only with polygons while they
		should also deal with lines.
	Many functions have a certain tolerance for numerical
		deviations; however, this tolerance assumes a certain
		scale (coordinates not much greater than 10 and not
		much smaller than 1/10).
	The output file consists of a single long line with no
		newlines, and therefore is difficult to view in
		Mathematica and some text editors.
	Nothing has been seriously tested.
*)

(* :Discussion: 
	This package is still under construction.
	
	Martin Kraus's polygon triangulation routines available from
	
	  http://library.wolfram.com/packages/polygontriangulation/
	  
	may be used in addition to or instead of some of the
	routines included here.
	
	Please send bug reports, comments and suggestions to 
	Martin.Kraus@informatik.uni-stuttgart.de
*)

(*:References: *)


BeginPackage["LiveGraphics3D`",
	{"Graphics`Animation`","Graphics`Graphics3D`","Graphics`Shapes`"}]


CenterAtCenterOfMass::usage =
"CenterAtCenterOfMass[graphics3d] resets the plot range such
that the center of mass is in the center of the bounding
box. (All polygons are assumed to have equal weights.)"

ClipAtPlotRange::usage =
"ClipAtPlotRange[graphics3d] clips all polygons at the bounding box.";

Dir::usage =
"Dir->directory is an option to WriteLiveForm and WriteLiveAnimationForm
which specified in which directory the output file should be placed."

DropInnerPolygons::usage =
"DropInnerPolygons[graphics3d] tries to drop all polygons which are
not part of the surface of the one solid in graphics3d.  It
also reverses the order of points in some polygons such that the
same face (BACK or FRONT) will be on the surface of the solid.
Therefore, DropInnerPolygons is also useful for solids without any
inner polygons (like eww's obelisk.m example).  The dropping of inner 
polygons should reduce the size of certain large graphics dramatically.
DropInnerPolygons must often be called using the syntax
DropInnerPolygons[NonIntersectingPolygons[graphics3D]], and often
best results are obtained with applet parameters VISIBLE_FACES=BACK
or VISIBLE_FACES=FRONT.  DropInnerPolygons returns some progress information
as it is running."

Graphics2DTo3D::usage =
"Graphics2DTo3D[graphics2D] gives a flat 3D object from a 2D object.
Now that Mathematica 4.0 supports animiated GIF exporting, this function
isn't needed so much."

LinelessFaces::usage =
"LinelessFaces[graphics3D] returns the solid with only exterior edges
included."

LiveAnimationForm::usage =
"LiveAnimationForm[list[graphics3d]] converts the list of
graphics3d objects to a form appropriate for an animation
to be displayed by the LiveGraphics3D applet."

LiveForm::usage =
"LiveForm[graphics3d] converts g to a form appropriate for a rotatable 
3-D solid to be displayed by the LiveGraphics3D applet."

NonIntersectingPolygons::usage =
"NonIntersectingPolygons[graphics3D] breaks up all intersecting polygons 
in such a manner that the resulting polygons do not intersect each other.
NonIntersectingPolygons can take a long time to run (i.e., hours) for 
complicated solids such as complex uniform polyhedra."

NonIntersectingLines::usage =
"NonIntersectingLines[graphics3D] breaks up all lines which intersect
polygons."

PlotCentersOfMasses::usage =
"PlotCentersOfMasses[graphics3d, p] shows a plot of points 
which shows the distance of the center of masses of all polygons
to p. (The y coordinate of the points is random.)"

PolyhedronPolygon::usage =
"PolyhedronPolygon[polyhedron,polygon list] returns a Graphics3D
object representing polyhedron with a line on the surface.
The polygon list is a pure list of vertices of the intersecting plane,
which is assumed to contain and be symmetric with respect to the origin."

SplitLines::usage =
"SplitLines[graphics3d,n] splits all lines into smaller lines.
Each original straight line is broken into n pieces."

Triangulate::usage =
"Triangulate[g,n] takes a Graphics3D polyhedron and triangulates all 
faces with n sides."

UnifyPolygons::usage = 
"UnifyPolygons[graphics3d] unifies neighbouring polygons which are
in the same plane.  Example: 
DropInnerPolygons[NonIntersectingPolygons[UnifyPolygons[Graphics3D[NumberedPolyhedron[21]]]]]."

WriteLiveAnimationForm::usage =
"WriteLiveAnimationForm[file,g,<opts>] takes a list of Graphics3D
objects and writes an animation to file which can be displayed by
the LiveGraphics3D applet.  Any 2D graphics must first be converted
using Graphics2DTo3D."

WriteLiveForm::usage =
"WriteLiveForm[file,g,<opts>] takes a Graphics3D object or list of such 
objects and writes an animation to file which can be displayed by
the LiveGraphics3D applet.  Any 2D graphics must first be converted
using Graphics2DTo3D."


(************* set your default directory for output files here *************)

Options[LiveForm]={
  Dir->""
};


Begin["`Private`"]

addn[n_,k_,m_]:=Mod[n+k-1,m]+1;


(* Under certain circumstances, e.g.,
 NonIntersectingPolygons[Graphics3D[Dual[NumberedPolyhedron[18]]]]][[1]]],
 a couple BoundingBox[{}]'s seem to get thrown in
 *)

BoundingBox[{}]:={}

BoundingBox[Polygon[l_List]]:=
  Module[{xs,ys,zs},{xs,ys,zs}=
    Transpose[l];{{Min[xs],Max[xs]},{Min[ys],Max[ys]},{Min[zs],Max[zs]}}];

BoundingBox[Line[l_List]]:=
  Module[{xs,ys,zs},{xs,ys,zs}=
    Transpose[l];{{Min[xs],Max[xs]},{Min[ys],Max[ys]},{Min[zs],Max[zs]}}];

BreakLineWithPolygon[{Line[l1_List],bb1_List},{p_Polygon,bb_List}]:=
  Module[{i,breakpoint,broken={},part,rest},
    i=1;
    rest=l1;
    While[i<Length[rest],
      breakpoint=PointOfIntersection[rest[[i]],rest[[i+1]],{p,bb}];
      If[{}!=breakpoint,
        part=Append[Take[rest,{1,i}],breakpoint];
        broken=Append[broken,{Line[part],BoundingBox[Line[part]]}];
        rest=Prepend[Take[rest,{i+1,Length[rest]}],breakpoint];
        i=0;
      ];
      i=i+1;
    ];
    broken=Append[broken,{Line[rest],BoundingBox[Line[rest]]}];
  	 Return[broken]];

BreakLineWithPolygons[{l1_Line,bb1_List},pl_List]:=
  Module[{broken={{l1,bb1}},i=1,l},
    While[i<=Length[pl],
      broken=broken/.{{Line[l_],box_List}:>
        BreakLineWithPolygon[{Line[l],box},pl[[i]]]};
      i=i+1];Flatten[broken/.{{Line[l_],_List}:>Line[l]}]];

BreakPolygonWithPolygon[{Polygon[l1_List],bb1_List},{Polygon[l2_List],
  bb2_List}]:=Module[{line,s1,s2,p1,p2,i1,i2,h,fulll1,result},
  	If[!OverlappingBoxesQ[bb1,bb2],Return[{{Polygon[l1],bb1}}]];
    fulll1=Flatten[{l1,l1},1];
    line=LineOfIntersection[Polygon[l1],Polygon[l2]];
    If[line=={},Return[{{Polygon[l1],bb1}}]];
    s1=ParametersOfClosestPoints[line,Polygon[l1]];
    s2=ParametersOfClosestPoints[line,Polygon[l2]];
    If[s1=={}||s2=={},Return[{{Polygon[l1],bb1}}]];
    If[s1[[2,1]]<s2[[1,1]]||s2[[2,1]]<s1[[1,1]],Return[{{Polygon[l1],bb1}}]];
    p1=line[[1]]+s1[[1,1]]*line[[2]];p2=line[[1]]+s1[[2,1]]*line[[2]];
    i1=Mod[s1[[1,2]]-1,Length[l1]]+1;i2=Mod[s1[[2,2]]-1,Length[l1]]+1;
    If[i1>i2,h=i1;i1=i2;i2=h;h=p1;p1=p2;p2=h];
	result={DropDoubledPoints[Polygon[Flatten[
		{{p1},Take[fulll1,{i1+1,i2}],{p2}}
	  ,1]]],DropDoubledPoints[Polygon[Flatten[
      	{{p2},Take[fulll1,{i2+1,i1+Length[l1]}],{p1}}
      ,1]]]};
    {{result[[1]],BoundingBox[result[[1]]]},
    {result[[2]],BoundingBox[result[[2]]]}}];

BreakPolygonWithPolygons[{p_Polygon,bb_List},pl_List]:=
  Module[{broken={{p,bb}},i=1,l},
    While[i<=Length[pl],
      broken=broken/.{{Polygon[l_],box_List}:>
        BreakPolygonWithPolygon[{Polygon[l],box},pl[[i]]]};
      i=i+1];Flatten[broken/.{{Polygon[l_],_List}:>Polygon[l]}]];

CenterOfMass[Polygon[l_]]:=Plus@@N[l]/Length[l];

CenterOfMass[g_Graphics3D]:=
  Module[{com,p},
    com=(Cases[g[[1]],_Polygon,Infinity]/.p_Polygon:>CenterOfMass[p]);
    Plus@@com/Length[com]];
    
CenterAtCenterOfMass[g_Graphics3D,opts___]:=
  Module[{com,pr,oldcenter,p,i},
    com=(Cases[g[[1]],_Polygon,Infinity]/.p_Polygon:>CenterOfMass[p]);
    com=Plus@@com/Length[com];pr=FullOptions[g,PlotRange];
    oldcenter=Apply[Plus,pr,2]/2.;i=1;
    While[i<=3,
      If[oldcenter[[i]]<com[[i]],
        pr[[i,2]]=pr[[i,2]]+2(com[[i]]-oldcenter[[i]]),
        pr[[i,1]]=pr[[i,1]]-2(oldcenter[[i]]-com[[i]])];i=i+1];
    Show[g,opts,PlotRange->pr]
]

centroid[l_List]:=(Plus@@l)/Length[l];

circle[{x_,y_},{a_,b_}]:=Module[{i,n=25},
	Line[Table[{x+a Cos[2.Pi i/n],y+b Sin[2.Pi i/n]},{i,0,n}]]
]

circle[{x_,y_},r_]:=circle[{x,y},{r,r}]

ClipAtPlotRange[SurfaceGraphics[gl_,opts___]]:=
	ClipAtPlotRange[Graphics3D[SurfaceGraphics[gl,opts]]]

ClipAtPlotRange[Graphics3D[gl_,opts___]]:=Module[
	{g,l,xmin,xmax,ymin,ymax,zmin,zmax,polys,pr,polysbbs},
	Off[Solve::svars];
    Off[RowReduce::luc];
    pr={{xmin,xmax},{ymin,ymax},{zmin,zmax}}=
        FullOptions[Graphics3D[gl,opts],PlotRange];
    polys={Polygon[{{xmin,ymin,zmin},{xmin,ymin,zmax},{xmin,ymax,zmax},{xmin,
              ymax,zmin}}],
        Polygon[{{xmin,ymin,zmin},{xmin,ymin,zmax},{xmax,ymin,zmax},{xmax,
              ymin,zmin}}],
        Polygon[{{xmin,ymin,zmin},{xmin,ymax,zmin},{xmax,ymax,zmin},{xmax,
              ymin,zmin}}],
        Polygon[{{xmax,ymin,zmin},{xmax,ymin,zmax},{xmax,ymax,zmax},{xmax,
              ymax,zmin}}],
        Polygon[{{xmin,ymax,zmin},{xmin,ymax,zmax},{xmax,ymax,zmax},{xmax,
              ymax,zmin}}],
        Polygon[{{xmin,ymin,zmax},{xmin,ymax,zmax},{xmax,ymax,zmax},{xmax,
              ymin,zmax}}]};
    polysbbs=polys/.{p_Polygon:>{p,BoundingBox[p]}};
    g=gl/.{p_Polygon:>BreakPolygonWithPolygons[{p,BoundingBox[p]},
          polysbbs]};
    grownpr=Transpose[{
          Transpose[pr][[1]]-0.001*(Transpose[pr][[2]]-Transpose[pr][[1]]),
          Transpose[pr][[2]]+0.001*(Transpose[pr][[2]]-Transpose[pr][[1]])}];
    g=g/.{p_Polygon:>If[IsPointOfPolygonOutside[p,grownpr],{},p]};
    On[Solve::svars];
    On[RowReduce::luc];
    Graphics3D[g,opts]];

CosAngleBetweenPolygons[normal1_,sign1_,dir1_,normal2_,dir2_]:=
	Module[{normal2p,into1,cosine},
		If[dir1.dir2<0.,
			normal2p=normal2,
			normal2p=-normal2];
		into1=Cross[normal1,dir1];
		cosine=N[Re[ArcCos[-normal1.normal2p]]];
		If[normal2p.into1<0.,cosine=2.*Pi-cosine];
		If[sign1<0,cosine=2.*Pi-cosine];
		cosine];

disk[{x_,y_},{a_,b_}]:=Module[{i,n=25},
	{EdgeForm[],Polygon[Table[{x+a Cos[2.Pi i/n],y+b Sin[2.Pi i/n]},{i,n}]]}
]

disk[{x_,y_},r_]:=disk[{x,y},{r,r}]

Distance[p1_List,p2_List]:=N[Sqrt[(p1-p2).(p1-p2)]];

DropDoubledPoints[Polygon[l_]]:=
  Module[{maxdist,newl={Last[l]},i=1},
    maxdist=10.^-8*
        Max[Map[((#[[1]]-#[[2]]).(#[[1]]-#[[2]]))&,
            Transpose[{l,RotateRight[l]}]]];
    While[i<=Length[l],
      If[(l[[i]]-Last[newl]).(l[[i]]-Last[newl])>maxdist,
        newl=Append[newl,l[[i]]]];i=i+1];
    newl=Drop[newl,1]; If[Length[newl]>2,Polygon[newl],{}]];

DropInnerPolygons[g_Graphics3D]:=
  Module[{startpolys,xcoors,minx,current,polys,normals,edges,signs,mainpolygon,
      mainsign,surfaces,checkedsurfaces,i,gp},polys=Cases[g,_Polygon,Infinity];
    normals=Map[NormalVector,polys];edges=Map[Edges,polys];
    boxes=Map[BoundingBox,polys];
		signs=Table[1,{Length[polys]}];
		checkedsurfaces={};
		Print["There are a total of ",Length[polys]," polygons."];
		minx=Min[Transpose[Transpose[boxes][[1]]][[1]]]+10.^-4;
		i=1;startpolys={};
    While[i<Length[polys],
      If[boxes[[i,1,1]]<=minx,xcoors=Transpose[polys[[i,1]]][[1]];
        If[Length[xcoors]>2,
          startpolys=
            Append[startpolys,{Sort[xcoors][[2]],-Abs[normals[[i,1]]],i}]]];
      i=i+1];
		surfaces={Sort[startpolys][[1,3]]};
		i=0;
		While[Length[surfaces]>0,
			current=First[surfaces];
			i=i+1;
			If[i<4||Mod[i,10]==0,
			Print["processing polygon no.",current, " (",Length[checkedsurfaces], 
          " processed, ",Length[surfaces]," in buffer)"],
				If[i==4,Print["Reporting only every 10th."]]];
			{newsurfs,signs}=
        NeighbouringSurfaces[current,polys,normals,boxes,edges,signs,surfaces,
          checkedsurfaces];
			checkedsurfaces=Append[checkedsurfaces,First[surfaces]];
			surfaces=Join[Drop[surfaces,1],newsurfs];
				(*Show[Graphics3D[{Extract[polys,Partition[checkedsurfaces,1]]}],
          Axes->True];*)
			];
		Print["Returning ",Length[checkedsurfaces]," of ",Length[polys],
      " polygons."];
      mainpolygon=polys[[First[checkedsurfaces]]];
      mainsign=Sign[CenterOfMass[mainpolygon].NormalVector[mainpolygon]];
		i=1;gp=g;
		While[i<=Length[polys],
			pos=Position[gp,polys[[i]]];
			If[!MemberQ[checkedsurfaces,i],
				gp=ReplacePart[gp,{},pos],
				If[signs[[i]]*mainsign<0,
					gp=ReplacePart[gp,Polygon[Reverse[polys[[i,1]]]],pos]]];
			i=i+1];
		gp];
		
Edges[p_Polygon]:=
  Transpose[{Transpose[{p[[1]],RotateLeft[p[[1]]]}],
      Map[Normalize,RotateLeft[p[[1]]]-p[[1]]]}];

Graphics2DTo3D[g_Graphics,r_:1,opts___]:=Module[{},
	Off[StackGraphics::bdpts];
	Show[StackGraphics[{g} /. {Circle->circle,Disk->disk}],
		opts,ViewPoint->{0,-3,0},Boxed->False,Axes->False,Ticks->None,
		BoxRatios->{1,1,r},PlotRange->All]
]

Graphics2DTo3D[g_List,r_:1,opts___]:=Module[{i},
	Table[Graphics2DTo3D[g[[i]],r,opts],{i,Length[g]}]
]

Graphics2DTo3D[gl_GraphicsArray,r_:1,opts___]:=Module[{i,g=Flatten[gl[[1]]]},
	Table[Graphics2DTo3D[g[[i]],r,opts],{i,Length[g]}]
]

IsPointOutside[{x_,y_,z_},{{xmin_,xmax_},{ymin_,ymax_},{zmin_,zmax_}}]:=
  If[x<xmin||x>xmax||y<ymin||y>ymax||z<zmin||z>zmax,True,False];

IsPointOfPolygonOutside[Polygon[l_List],plotrange_List]:=
  Or@@Map[IsPointOutside[#,plotrange]&,l];
  
line[l_List]:=Line[Join[l,{l[[1]]}]]

LinelessFaces[g_]:=Module[
	{g1=Show[g,DisplayFunction->Identity],
	g2,g3,g4
	},
	Print["Step 1 of 4. Splitting intersecting and dropping inner polygons..."];
	g2=DropInnerPolygons[NonIntersectingPolygons[g]];
	Print["Step 2 of 4. Re-unifying polygons..."];
	g3=UnifyPolygons[g2];
	Print["Step 3 of 4. Determining Edges..."];
	g4=Show[SplitLines[g3/.Polygon[l_]:>Line[Append[l,First[l]]],1],
		DisplayFunction->Identity];
	Print["Step 4 of 4. Combining faces and edges..."];
	Graphics3D[{EdgeForm[],g2[[1]],g4[[1]]}]
]

LineOfIntersection[Polygon[l1_List],Polygon[l2_List]]:=
  Module[{n1=NormalVector[Polygon[l1]],n2=NormalVector[Polygon[l2]],n,p1,p2,
      p3},n=Cross[n1,n2];If[n.n<10.^-5,Return[{}]];n=n/Sqrt[n.n];
    sol=Solve[{n1.({p1,p2,p3}-l1[[1]])==0,n2.({p1,p2,p3}-l2[[1]])==0,
          n.({p1,p2,p3}-l1[[1]])==0},{p1,p2,p3}];
    If[Length[sol]=!=1||Length[sol[[1]]]=!=3,
      Return[{}]];{{p1,p2,p3}/.sol[[1]],n}];

LiveAnimationForm[gl_]:=
	HoldForm[ShowAnimation][LiveForm[gl]]

LiveAnimationForm[gl_,opts___]:=
	HoldForm[ShowAnimation][LiveForm[gl],LiveForm[opts]]

LiveForm[gr_]:=Module[
	{g=Switch[gr,
		_List,gr,
		_Graphics3D,gr,
		_ContourGraphics,Graphics3D[SurfaceGraphics[gr]],
		_DensityGraphics,Graphics3D[SurfaceGraphics[gr]],
		_SurfaceGraphics,Graphics3D[gr],
		_,gr
	]},
	NumberForm[InputForm[N[g]],5]
]

MergePolygons[p1_Polygon,e1_Integer,p2_Polygon,e2_Integer]:=
  Module[{newpol,dif,i,d,maxdist=10.^-8},
    newpol = DropDoubledPoints[
    Polygon[Flatten[{Take[p1[[1]],{1,e1}],
          Take[p2[[1]],{e2+1,Length[p2[[1]]]}],Take[p2[[1]],{1,e2}],
          Take[p1[[1]],{e1+1,Length[p1[[1]]]}]},1]]];
  i=1;
  While[i <= Length[newpol[[1]]],
    dif=newpol[[1,Mod[i-1-1,Length[newpol[[1]]]]+1]]-
    newpol[[1,Mod[i+1-1,Length[newpol[[1]]]]+1]];
    If[dif.dif<maxdist,
      newpol=Polygon[Delete[newpol[[1]],{{i},{Mod[i+1-1,Length[newpol[[1]]]]+1}}]];
      i=Mod[i-1-1,Length[newpol[[1]]]]+1,
      i=i+1]
  ];
  newpol
]    

NeighbouringSurfaces[pindex_,polys_,normals_,boxes_,edges_,signs_,surfaces_,
    checkedsurfaces_]:=
	Module[{newsigns=signs,result={},opolys,oedges,i},
		opolys=Position[boxes,box_/;OverlappingBoxesQ[boxes[[pindex]],box,10.^-4],
        1];
		i=1;
		While[i<=Length[edges[[pindex]]],
			currentedge=edges[[pindex,i]];
			oedges=
        Map[({#[[1]],
                OverlappingEdgesPositions[edges[[#[[1]]]],currentedge]})&,
          opolys];
			(* oedges has now elements of the form 
					{index of polygon,{positions of overlapping edges}}*)
		 oedges=
        Map[(If[#[[2]]=={}||#[[1]]==pindex,{},{{
                    CosAngleBetweenPolygons[normals[[pindex]],signs[[pindex]],
                      currentedge[[2]],normals[[#[[1]]]],
                      edges[[#[[1]],#[[2,1,1]],2]]],#[[1]],
                    currentedge[[2]].edges[[#[[1]],#[[2,1,1]],2]]}}])&,
          oedges];
			oedges=Sort[Flatten[oedges,1]];
			If[Length[oedges]>0&&!MemberQ[surfaces,oedges[[1,2]]]&&!
            MemberQ[checkedsurfaces,oedges[[1,2]]],
        result=Append[result,oedges[[1,2]]];
        newsigns[[oedges[[1,2]]]]=-signs[[pindex]]*Sign[oedges[[1,3]]]];
			i=i+1];
		{result,newsigns}];

NonIntersectingPolygons[Graphics3D[gl_]]:=
	NonIntersectingPolygons[Graphics3D[gl,{}]]

NonIntersectingPolygons[Graphics3D[gl_,opts___]]:=Module[{g,l,polys},
	Off[Solve::svars];
	Off[RowReduce::luc];
    polys=Cases[gl,_Polygon,Infinity];
    polys=polys/.Polygon[l_]:>{Polygon[l],BoundingBox[Polygon[l]]};
    g=Graphics3D[gl/.{Polygon[l_]:>BreakPolygonWithPolygons[
    	{Polygon[l],BoundingBox[Polygon[l]]},polys]},
      opts];
	On[Solve::svars];
	On[RowReduce::luc];
	g
];

NonIntersectingLines[Graphics3D[gl_,opts___]]:=Module[{g,l,polys},
	Off[Solve::svars];
	Off[RowReduce::luc];
    polys=Cases[gl,_Polygon,Infinity];
    polys=polys/.Polygon[l_]:>{Polygon[l],BoundingBox[Polygon[l]]};
    g=Graphics3D[gl/.{Line[l_]:>
      BreakLineWithPolygons[{Line[l],BoundingBox[Line[l]]},polys]},
      opts];
	On[Solve::svars];
	On[RowReduce::luc];
	g
];

Normalize[v_List]:=
  Module[{length=N[Sqrt[v.v]]},If[length<10.^-20,v,v/length]];
  
NormalVector[Polygon[l_List]]:=
  Module[{i=1,k=1,v},
    While[i<Length[l]&&k<3,
      If[Distance[l[[i]],l[[i+1]]]>10.^-5,v[k]=l[[i+1]]-l[[i]];k=k+1];i=i+1];
    If[k!=3,{1,1,1},Normalize[Cross[v[1],v[2]]]]];

OverlappingBoxesQ[{{xmin1_,xmax1_},{ymin1_,ymax1_},{zmin1_,zmax1_}},{{xmin2_,
        xmax2_},{ymin2_,ymax2_},{zmin2_,zmax2_}},tol_:0]:=
  If[(xmin1>xmax2+tol||xmax1+tol < xmin2)||(
        ymin1>ymax2+tol||ymax1+tol < ymin2)||(
        zmin1>zmax2+tol||zmax1+tol < zmin2),False,True];

OverlappingEdgesPositions[edges_,{{p1a_,p1b_},dir1_}]:=
  Position[edges,{{p2a_List,p2b_List},dir2_List}/;
      OverlappingEdgesQ[{{p1a,p1b},dir1},{{p2a,p2b},dir2}]];
      
OverlappingEdgesQ[{{p1a_,p1b_},dir1_},{{p2a_,p2b_},dir2_}]:=
	Module[{dot},
		If[Plus@@Abs[dir1-dir2]<10.^-4&&
        Plus@@Abs[(p2a-p1a)-((p2a-p1a).dir1)*dir1]<10.^-4,
			dot=(p1b-p1a).(p2a-p1a);
			If[dot>0,
				If[dot >(p1b-p1a).(p1b-p1a)-10.^-4,
					False,
					True],
				dot=(p1b-p1a).(p2b-p1a);
				If[dot>10.^-4,
					True,
					False]],
			If[Plus@@Abs[dir1+dir2]<10^-4&&
          	Plus@@Abs[(p2a-p1a)-((p2a-p1a).dir1)*dir1]<10.^-4,
				dot=(p1b-p1a).(p2b-p1a);
				If[dot>0,
					If[dot>(p1b-p1a).(p1b-p1a)-10.^-4,
						False,
						True],
					dot=(p1b-p1a).(p2a-p1a);
					If[dot >10.^-4,
						True,
						False]],
				False]]];      

ParameterOfClosestPoint[{p1_List,v1_List},{p2_List,v2_List}]:=
  Module[{},-(Cross[v1,v2].Cross[p1-p2,v2])/(Cross[v1,v2].Cross[v1,v2])];

ParameterOfClosestPointToPoint[{p1_List,v1_List},p_List]:=(p-p1).v1/Sqrt[v1.v1];

ParametersOfClosestPoints[{p_List,v_List},Polygon[l_List]]:=
  Module[{n=NormalVector[Polygon[l]],i,dir,fulll=Flatten[{l,l},1],oldsign,
      newsign,index,parameter},dir=Cross[v,n];i=2;k=1;
    oldsign=TolerantSign[(fulll[[i]]-p).dir];
    While[i<=Length[fulll]&&k<3,newsign=TolerantSign[(fulll[[i]]-p).dir];
      If[newsign!=0&&oldsign!=0&&newsign!=oldsign,index[k]=i-1;
        parameter[k]=
          ParameterOfClosestPoint[{p,v},{fulll[[i-1]],
              fulll[[i]]-fulll[[i-1]]}];k=k+1];If[newsign!=0,oldsign=newsign];
      i=i+1];If[k!=3,{},
      Sort[{{parameter[1],index[1]},{parameter[2],index[2]}}]]];

PlotCentersOfMasses[g_Graphics3D,origin_]:=
  Module[{pl,p},
    pl=Cases[g[[1]],_Polygon,Infinity]/.p_Polygon:>CenterOfMass[p];
    pl=Map[(Distance[origin,#])&,pl];
		ListPlot[Transpose[{pl,Table[Random[],{Length[pl]}]}],Axes->{True,False}]];

PointInBoundingBoxQ[{x_,y_,z_}, {{xmin_,xmax_},{ymin_,ymax_},{zmin_,zmax_}},tol_:0]:=
  If[(xmin-tol<=x&&x<=xmax+tol)&&(ymin-tol<=y&&y<=ymax+tol)&&(zmin-tol<=z&&z<=zmax+tol),True,False];

PointOfIntersection[l1_,l2_,{poly_Polygon,bbox_}]:=
  Module[{l,p=poly[[1,1]],normal=NormalVector[poly],denom},
    denom=((l2-l1).normal);
    If[Abs[denom]<0.000001,Return[{}]];
    l=((p-l1).normal)/denom;
    point=l1+l*(l2-l1);
    If[l<=0.000001||l>=0.999999||!PointInBoundingBoxQ[point,bbox],{},point]
  ];

PolygonInSamePlane[pindex_,polys_,normals_,boxes_,edges_,checkedsurfaces_]:=
  Module[{result={},opolys,oedges,i},
		opolys=Position[boxes,box_/;OverlappingBoxesQ[boxes[[pindex]],box,10.^-4],
        1];
		opolys=Intersection[opolys,
        Position[normals, 
          normal_/;Abs[Abs[normal.normals[[pindex]]]-1]<10.^-4,1]];
		If[Length[opolys]==0,Return[{}]];
		i=1;
		oedges={};
		While[i<=Length[edges[[pindex]]]&&Length[oedges]==0,
			currentedge=edges[[pindex,i]];
			oedges=
        Map[({#[[1]],
                OverlappingEdgesPositions[edges[[#[[1]]]],currentedge]})&,
          opolys];
			(* oedges has now elements of the form 
					{index of polygon,{positions of overlapping edges}}*)
		oedges=
        Map[(
              If[#[[2]]=={}||#[[1]]==pindex||
                  MemberQ[checkedsurfaces,#[[1]]],{},{{
                    i,#[[1]],#[[2,1,1]]}}])&, oedges];
		oedges=Flatten[oedges,1];
			i=i+1];
		If[Length[oedges]>0,oedges[[1]],{}]];

PolyhedronPolygon[solid_,polygon_List]:=Module[
	{g=Show[Graphics3D[{solid,Polygon[polygon]}],DisplayFunction->$DisplayFunction]},
	Graphics3D[
		NonIntersectingPolygons[Show[{DropInnerPolygons[Graphics3D[g[[1,1]]]],
		AffineShape[Graphics3D[g[[1,2]]],{1.1,1.1,1.1}]}]][[1,1]]
	]
]

SplitPoly[p_]:=Module[{i,c=centroid[p],n=Length[p],f=1},
	{Table[Polygon[{c+f(p[[i]]-c),
	c+f(p[[addn[i,1,n]]]-c),c}],{i,n}] (*,line[p] *)}
]
	
SplitLines[g_,n_Integer]:=
  Module[{gl},
    gl=g/.Line[l_List]:>Map[Line[#]&,Transpose[{Drop[l,-1],Drop[l,1]}]];
    If[n>1,gl=gl/.Line[{a_,b_}]:>Line[Table[N[a+i*(b-a)],{i,0,1,1/n}]];
      gl=gl/.Line[l_List]:>Map[Line[#]&,Transpose[{Drop[l,-1],Drop[l,1]}]];];
    gl];

TolerantSign[x_]:=If[Abs[x]<10.^-4,0,Sign[x]];

Triangulate[g_,s_]:=Module[{poly=g[[1]] /.Polygon->Identity,base,sides},
	{base=Select[poly,Length[#]==s&],sides=Select[poly,Length[#]!=s&]};
	Graphics3D[{{EdgeForm[],SplitPoly/@base},Polygon/@sides}]
]

UnifyPolygons[g_Graphics3D]:=
  Module[{current,polys,oldpolys,normals,edges,signs,surfaces,checkedsurfaces,
      mergedsurfaces,i,gp,mergethis},polys=Cases[g,_Polygon,Infinity];
    oldpolys=polys;normals=Map[NormalVector,polys];edges=Map[Edges,polys];
    boxes=Map[BoundingBox,polys];
		checkedsurfaces={};
		mergedsurfaces={};
		Print["There are a total of ",Length[polys]," polygons."];
		surfaces=Table[i,{i,Length[polys]}];
		i=1;
		While[Length[surfaces]>0,
			current=First[surfaces];
			(*If[i<4||Mod[i,10]==0,
			Print["processing polygon no.",current, " (",Length[checkedsurfaces], 
          " processed, ",Length[surfaces]," in buffer)"],
				If[i==4,Print["Reporting only every 10th."]]];*)
			If[MemberQ[checkedsurfaces,current],mergethis={},mergethis=
       PolygonInSamePlane[current,polys,normals,boxes,edges, 
            checkedsurfaces]];
			If[Length[mergethis]==0,
				checkedsurfaces=Append[checkedsurfaces,current];
				surfaces=Drop[surfaces,1];
				i=i+1,
				{edgeindex,pindex,pedgeindex}=mergethis;
				Print["Merging polygon no.",current," (edge ",edgeindex,") with no.",
          pindex, " (edge ",pedgeindex,")"];
				polys[[current]]=
          MergePolygons[polys[[current]],edgeindex,polys[[pindex]],
            pedgeindex];
				edges[[current]]=Edges[polys[[current]]];
				boxes[[current]]=BoundingBox[polys[[current]]];
				checkedsurfaces=Append[checkedsurfaces,pindex];
				mergedsurfaces=Append[mergedsurfaces,pindex]]
				(*Show[Graphics3D[{Extract[polys,Partition[checkedsurfaces,1]]}],
          Axes->True];*)
			];
		i=1;gp=g;
		While[i<=Length[polys],
			pos=Position[gp,oldpolys[[i]]];
			If[MemberQ[mergedsurfaces,i],
				gp=Delete[gp,pos],
				gp=ReplacePart[gp,polys[[i]],pos]];
			i=i+1];
		gp];

WriteLiveAnimationForm[filename_,gl_,opts___]:=Module[
         {outfile,dir=Dir/.{opts}/.Options[LiveForm]},
         Unprotect[Real];
         Format[x_Real, InputForm] := OutputForm[NumberForm[x, 5,
            NumberFormat -> (If[#3 == "", #1, SequenceForm[#1, "*^", #3]] &)]];
         Protect[Real];
         outfile=dir<>filename;
         WriteString[outfile,
         ToString[HoldForm[ShowAnimation][LiveForm[gl],LiveForm[opts]]]];
         Close[outfile];
         Unprotect[Real]; Format[x_Real, InputForm] =.; Protect[Real];
]

WriteLiveAnimationForm[filename_,gl_]:=Module[
         {outfile,dir=Dir /. Options[LiveForm]},
         Unprotect[Real];
         Format[x_Real, InputForm] := OutputForm[NumberForm[x, 5,
            NumberFormat -> (If[#3 == "", #1, SequenceForm[#1, "*^", #3]] &)]];
         Protect[Real];
         outfile=dir<>filename;
         WriteString[outfile,
         ToString[HoldForm[ShowAnimation][LiveForm[gl]]]];
         Close[outfile];
         Unprotect[Real]; Format[x_Real, InputForm] =.; Protect[Real];
]

WriteLiveForm[filename_,g_,opts___]:=Module[
         {outfile,dir=Dir/.{opts}/.Options[LiveForm]},
         Unprotect[Real];
         Format[x_Real, InputForm] := OutputForm[NumberForm[x, 5,
            NumberFormat -> (If[#3 == "", #1, SequenceForm[#1, "*^", #3]] &)]];
         Protect[Real];
         outfile=dir<>filename;
         WriteString[outfile,ToString[LiveForm[g]]];
         Close[outfile];
         Unprotect[Real]; Format[x_Real, InputForm] =.; Protect[Real];
]


End[]

(* Protect[  ] *)

EndPackage[]

(*---------------------end of LiveGraphics3D.m-----------------------*)
