feq[x_,y_,eps_:10^-10]:= Plus @@ Abs[x-y]0, Return[{}] ]; If[sign1==sign2==0, Return[Indeterminate] ]; {A,B,C}=Coefficient[plane[x,y,z],{x,y,z}]; D=plane[x,y,z]-(A x+B y+C z); pt=LinearSolve[ {{A,B,C,0},{1,0,0,x1-x2},{0,1,0,y1-y2},{0,0,1,z1-z2}}, {-D,x1,y1,z1}]; pt=Take[pt,3]; If[feq[pt,pt1], Return[pt1] ]; If[feq[pt,pt2], Return[pt2] ]; pt] cut[polys_,plane_]:=Module[ {cutpoly,newfaces, newaps,newpolys,chpolys}, cutpoly[poly_]:=Module[ {sides,inters,newpoly}, sides=Transpose[{poly,RotateLeft[poly]}]; inters=inter[#,plane]& /@ sides; newaps=Join[newaps,inters]; Do[ If[inters[[i]]==Indeterminate, newaps=Join[newaps,sides[[i]]] ], {i,Length[poly]} ]; newpoly=Flatten[ Transpose[{poly,inters}],1]; newpoly=DeleteCases[newpoly, {}|Indeterminate| (p_/;fsign[plane @@ p]>0)]; If[newpoly=={}, Return[{}] ]; newpoly=#[[1]]& /@ Split[newpoly]; If[newpoly[[1]]==newpoly[[-1]], newpoly=Drop[newpoly,-1] ]; If[!MatchQ[inters,{{}..}], AppendTo[chpolys,newpoly] ]; newpoly]; newfaces[newaps_,chpolys_]:=Module[ {nextap, k,n,ap,res,face}, nextap[ap1_:{Infinity,0,0},ap2_]:=Module[ {adjfaces,apexes}, adjfaces=Select[chpolys, !MemberQ[#,ap1]&&MemberQ[#,ap2]&]; apexes=Flatten[adjfaces,1]; apexes=Select[apexes, #!=ap2&&MemberQ[newaps,#]&,1]; If[apexes=={}, {}, apexes[[1]] ] ]; res={}; n=Length[newaps]; face=Take[newaps,1]; AppendTo[face, nextap[face[[-1]] ] ]; k=2; While[True, ap=nextap[face[[-2]],face[[-1]] ]; If[ap=={}, Return[res] ]; AppendTo[face,ap]; k++; If[face[[1]]!=face[[-1]], Continue[] ]; AppendTo[res, Drop[face,-1] ]; If[--k==n, Return[res] ]; face=Select[newaps, !MemberQ[res,#,{2}]&,1]; AppendTo[face, nextap[face[[-1]] ] ]; k+=2] ]; newaps=chpolys={}; newpolys=cutpoly /@ polys; newaps=DeleteCases[newaps, {}|Indeterminate]; newaps=Union[newaps]; newpolys=DeleteCases[newpolys,{}]; Join[newpolys, newfaces[newaps,chpolys] ] ] {feq,fsign,inter,cut} (* Copyright (C) 1998, Maxim Rytin *)