(* Convex Hull in 3D *) (* Wouter Meeussen, 25/04/1998, 09/05/2002. *) (* About *) (* Very public domain, no rights to be reserved whatsoever. Feed back, corrections, enhancements, speed-up and the like are appreciated at : Wouter Meeussen, Testing, suggestions and encouragements by Russell Towle . Package design along Roman E. Maeder's "Template.nb". *) (* Limitations *) (* Points should be unique. Duplicate points are removed. HullVolume, HullArea and HullShow were adapted to operate on the original non-unique point set. Line added 09/05/2002 to handle a coplanar point set as a flattened 3D hull. Input points can be either Reals or Integers, or expressions that evaluate to numbers. The calculation converts its input to Reals, but gives output as indices pointing to the input-point list. For exact input, the HullVolume and HullArea need be FullSimplified to produce a standard mathematical form. *) (* Package *) (* This part declares the publicly visible functions, options, and values. *) (* Set up the package context, including public imports *) BeginPackage["DiscreteMath`ConvexHull3D`"] (* Usage messages for the exported functions and the context itself *) ConvexHull3D::usage = "ConvexHull3D.m is a package that calculates the 3D convex hull of a set of points in 3D. It returns the indices of the points ( their position in the input list) that form the 3D polygons that together form the polytope spanned by these points. Non spanning points ( along the edges or in the faces) are optional." ConvexHull3D::usage = "ConvexHull3D[pointlist:{{x,y,z}..}] calculates the 3D convex hull of a set of points in 3D. It returns the indices of the points ( their position in the input list) that form plane 3D polygons that together make up the polytope spanned by these points." HullVolume::usage = "HullVolume[ pointlist:{{x,y,z}..}, hull:{{_Integer..}..} ] calculates the volume of the convex hull around pointlist as defined by hull." HullArea::usage = "HullArea[ pointlist:{{x,y,z}..}, hull:{{_Integer..}..} ] calculates the area of the convex hull around pointlist as defined by hull." HullShow::usage = "HullShow[ pointlist:{{x,y,z}..}, hull:{{_Integer..}..},options_:{} ] displays the convex hull." (* Error messages for the exported objects *) ConvexHull3D::badarg = "function `1` was called with argument `2`, a list of length-3 lists of numerical data was expected." HullVolume::badarg = "You twit, you called `1` with argument `2`!" HullArea::badarg = "You twit, you called `1` with argument `2`!" (* Begin the private context (implementation part) *) Begin["`Private`"] (* Definition of auxiliary functions *) myOrdering[z_List]:=(Sort[Transpose[{z,Range[Length[z]]}]]//Transpose) [[2]] norm[vec_List]:=Sqrt[vec.vec] bary[pts_]:=1/Length[pts] Plus@@pts uniq[li_List]:=Block[{f},f[g_]:=(f[g]=Sequence[];g);f/@li] (* For points "vi" in 3D space, given as vi : {Real_,Real_,Real_}, the following function defines the distance from v to the plane v1_v2_v3: *) (* plane[{v1_,v2_,v3_},v_]:=Det[{(v-v1),(v1-v2), (v2-v3)}] is faster with compile *) plane=Compile[{{u,_Real,2},{v,_Real,1}},Det[{(v-u[[1]]),(u[[1]]-u[[2]]), (u[[2]]-u[[3]])}]] (* this provides the angle between (v2-v1) and (v3-v2), usefull only in case of 4 or more coplanar points: *) ang2[{v1_List,v2_List},v2_]:=1. ;ang2[{v1_List,v2_List},v1_]:=-1. ang2[{v1_List,v2_List},v3_List]:=(v2-v1).(v3-v2)/norm[v2-v1]/norm[v3-v2] ang[{v1_List,v2_List},v2_|v1_]:=1. ang[{v1_List,v2_List},v3_List]:=(v2-v1).(v3-v2)/norm[v2-v1]/norm[v3-v2] (* the function ch2D handles coplanarity of 4 or more points: *) ch2D[indices_List]:=Module[{pts,cen,dis,idx,xfar,far,xsec,hull}, pts=points[[indices]];cen=bary[pts];dis=norm[#-cen]& /@ pts;idx=myOrdering[ dis];xfar=Last[idx];far=pts[[xfar]]; xsec=Sort[{100+ ang2[{cen,far},pts[[#]] ]-100. , norm[(far-pts[[#]])] ,# }&/@idx ][[-2,-1]]; hull={xfar,xsec}; FixedPoint[Last[AppendTo[hull, Sort[ {100+ang[ pts[[Take[hull,-2] ]] , pts[[#]] ]-100., norm[(pts[[ Last[ hull] ]] -pts[[#]])] ,# } &/@idx] [[-3,-1]] ] ]&,0,SameTest->(First[hull]==#2 &)]; indices[[Drop[hull,-1]]]] (* The plane has a sign, which we will need. All distances from points v on the same side of the plane {v1,v2,v3} have the same sign. It can be "+" or "-", dependent on the right or left circulation of the set {v1,v2,v3} as seen from the origin. *) test[tria:{_Integer,_Integer,_Integer}]:= Block[{},res={}; If[Unequal@@tria && Chop[norm[ Cross[points[[tria[[1]] ]]-points[[tria[[2]] ]], points[[tria[[1]] ]] -points[[tria[[3]] ]]]]]=!=0 , Block[{i=1,oneside=True}, While[(oneside= Or[FreeQ[res,-1],FreeQ[res,1]])&&i<=n , res={res, If[FreeQ[tria,i], Sign@Chop[ plane[points[[tria]],points[[i]] ] ],0] };i++]; oneside],False]] (* Definition of the exported functions *) (* ConvexHull3D[] *) ConvexHull3D[inputpoints:{{_?NumericQ,_?NumericQ,_?NumericQ}..}]:= Module[{cen,cenpoints,distx,far,crosx,nix,it,faces,oldlegs,legs,disordered, testing,testable,circulation}, time=TimeUsed[]; points=N[uniq@inputpoints]; n=Length[points]; cen=Plus@@points /n; cenpoints=(#-cen)&/@points; distx=Reverse[myOrdering[(#.#)&/@cenpoints]]; far=distx[[1]]; (* crosx=Reverse[myOrdering[(Abs[Cross[points[[far]]-cen,#]]&/@cenpoints)* ( (#.#)&/@cenpoints ) ]]; *) crosx= Reverse[myOrdering[1/(1+2 norm/@((points[[far]]-#)& /@ points)) * (norm[#-cen]& /@ points + norm/@( Cross[(points[[far]]-cen) ,(#-cen) ]& /@ points )) ]]; nix=True;i=1;badones=disordered={}; While[nix&&i<=n-1,j=1; While[ j<=n-1 && (nix=Not[test[it={far,crosx[[i]],crosx[[j]]}]]), j++]; i++]; If[Count[res,0,-1]>3, it=ch2D[disordered=Flatten[Position[Flatten[res],0]] ] ; badones=Complement[disordered,it]; ]; faces={it}; (* Print[TimeUsed[]-time,"ruff=",disordered,"faces=",faces,"bad=",badones]; *) If[n == Length[it] + Length[badones], Return[{it,FixedPoint[RotateLeft,#,SameTest->(First[#2]==Min[#2]&)]&@(Reverse@it)}]];(* added 09/05/2002 to handle a coplanar point set *) oldlegs=legs={}; faces=Union@Flatten[{faces,{it}},1]; legs=Split[Sort[Sort/@Flatten[Partition[#~Append~First[#],2,1]&/@faces,1]]]; legs=Flatten[DeleteCases[legs,{_List,_List}],1]; While[legs=!={}&&legs=!=oldlegs, oldlegs=legs; legs=Split[Sort[Sort/@Flatten[Partition[#~Append~First[#],2,1]&/@faces,1]]]; legs=Flatten[DeleteCases[legs,{q_List,q_List..}],1]; testing=First@ Select[faces, (!FreeQ[#, legs[[1,1]] ]&& !FreeQ[#,legs[[1,2]] ])&]; nix=True; j=1; testable= DeleteCases[distx,Alternatives@@Join[testing,badones]]; While[ j<=Length[testable] && (nix=Not[test[it=Flatten[{legs[[1]],testable[[j]]}] ] ]), j++];disordered={}; If[Count[res,0,-1]>3, it=ch2D[disordered=Flatten[Position[Flatten[res],0]] ] ; badones=Union@Flatten[{badones,Complement[disordered,it]}] ; ]; (* Print["now",legs[[1]],testing,"ruff",disordered, "it",it,"bad= ",badones] ; *) faces= Flatten[{faces,{it}},1]; legs=Split[Sort[Sort/@Flatten[Partition[#~Append~First[#],2,1]&/@faces,1]]]; legs=Flatten[DeleteCases[legs,{q_List,q_List..}],1]; (* Print["2:legs", legs,"\nfaces=",faces ] ; *) ];(*endwhile*) circulation=Sign[Chop[plane[Take[#,3]/.k_Integer:>points[[k]],cen]]]&/@faces; faces=MapThread[List,{faces,circulation}]/.{poly_List,1}:>Reverse[poly]/.{ poly_List,-1}:>poly; faces=Sort[FixedPoint[RotateLeft,#,SameTest->(First[#2]==Min[#2]&)]&/@faces]; If[legs==={},faces,{{}}] ] (* HullVolume[] *) HullVolume[ pointlist:{{_?NumericQ,_?NumericQ,_?NumericQ}..}, hull:{{_Integer..}..} ]:=Block[{cen,upointlist,n},upointlist=uniq[pointlist];n=Length[upointlist];cen=Plus@@upointlist/n; 1/6* Plus@@ Apply[(upointlist[[#1]]-cen).Cross[upointlist[[#2]]-cen,upointlist[[#3]]-cen]&, Flatten[ Thread[ z[First[#],Partition[Rest[#],2,1] ]] &/@ hull ,1]/.z[q__]:> Flatten[{q}] ,1] ] (* /;MatchQ[pointlist,{{_?NumericQ,_?NumericQ,_?NumericQ}..}]&&MatchQ[ hull,{{_Integer..}..}] || Message[ConvexHull3D::badarg, HullVolume, pointlist,hull] *) (* HullArea[] *) HullArea[ pointlist:{{_?NumericQ,_?NumericQ,_?NumericQ}..}, hull:{{_Integer..}..} ]:= Block[{upointlist},upointlist=uniq[pointlist];1/2*Plus@@Apply[norm[Cross[(upointlist[[#2]]-upointlist[[#1]]),(upointlist[[#3]]- upointlist[[#2]])]]&, Flatten[ Thread[ z[First[#],Partition[Rest[#],2,1] ]] &/@ hull ,1]/.z[q__]:> Flatten[{q}] ,1]] (*/;MatchQ[pointlist,{{_?NumericQ,_?NumericQ,_?NumericQ}..}]&&MatchQ[ hull,{{_Integer..}..}]|| Message[ConvexHull3D::badarg, HullArea, pointlist,hull] *) (* HullShow[] *) HullShow[ pointlist:{{_?NumericQ,_?NumericQ,_?NumericQ}..}, hull:{{_Integer..}..},options_:{} ]:= Block[{upointlist},upointlist=uniq[pointlist];Show[Graphics3D[{Map[Point, pointlist,{1}]}], Graphics3D[Polygon/@(upointlist[[#]]& /@hull)],Join[options,{Axes->True, AspectRatio -> Automatic ,BoxRatios->{1,1,1}}]] ] (* /;MatchQ[pointlist,{{_?NumericQ,_?NumericQ,_?NumericQ}..}]&&MatchQ[ hull,{{_Integer..}..}] || Message[ConvexHull3D::badarg, HullArea, pointlist,hull] *) (* End the private context *) End[ ] EndPackage[ ]