beta testing ConvexHull3D
- To: mathgroup@smc.vnet.net
- Subject: [mg12165] beta testing ConvexHull3D
- From: Wouter Meeussen <eu000949@pophost.eunet.be>
- Date: Fri, 1 May 1998 03:08:18 -0400
hi all,
Here is the package version of ConvexHull3D. I would appreciate your
beta-testing input before I commit it to mathsource.
Attached are a short notebook hulltester.nb and ConvexHull3D.M
enjoy,
wouter.
(***********************************************************************
Mathematica-Compatible Notebook
This notebook can be used on any computer system with Mathematica 3.0,
MathReader 3.0, or any compatible application. The data for the
notebook starts with the line of stars above.
To get the notebook into a Mathematica-compatible application, do one of
the following:
* Save the data starting with the line of stars above into a file
with a name ending in .nb, then open the file inside the application;
* Copy the data starting with the line of stars above to the
clipboard, then use the Paste menu command inside the application.
Data for notebooks contains only printable 7-bit ASCII and can be sent
directly in email or through ftp in text mode. Newlines can be CR, LF
or CRLF (Unix, Macintosh or MS-DOS style).
NOTE: If you modify the data for this notebook not in a Mathematica-
compatible application, you must delete the line below containing the
word CacheID, otherwise Mathematica-compatible applications may try to
use invalid cache data.
For more information on notebooks and Mathematica-compatible
applications, contact Wolfram Research:
web: http://www.wolfram.com
email: info@wolfram.com
phone: +1-217-398-0700 (U.S.)
Notebook reader applications are available free of charge from Wolfram
Research.
***********************************************************************)
(*CacheID: 232*)
(*NotebookFileLineBreakTest
NotebookFileLineBreakTest*)
(*NotebookOptionsPosition[ 3483, 124]*)
(*NotebookOutlinePosition[ 4132, 147]*) (*
CellTagsIndexPosition[ 4088, 143]*) (*WindowFrame->Normal*)
Notebook[{
Cell[BoxData[
\(\( (*\ SeedRandom@12\ *) \n
points = \
Union[\ Table[
Plus@@Table[\ Random[Integer, {\(-1\), 1}], {2}], {24},
{3}]]\)\)],
"Input"],
Cell["\<\
the cell below is set to non-evaluatable : you should first adapt the
correspondence between it and your storage \ directory.
You might run into \"path-problems\" here. I did!. So I copied the *.m
file \ into the DiscreteMath directory.\
\>", "Text"],
Cell[BoxData[
\(Needs["\<DiscreteMath`ConvexHull3D`\>"]\)], "Input",
Evaluatable->False],
Cell[BoxData[
\(\(?ConvexHull3D\)\)], "Input"],
Cell[CellGroupData[{
Cell["the ConvexHull3D:", "Subsubsection"],
Cell["\<\
the convex hull is the following set of polygons, the points given as \
indices to their positions in \"points\" :\ \>", "Text"],
Cell[BoxData[
\(Timing[ch3D = ConvexHull3D[points]]\)], "Input"],
Cell["\<\
these points are on the hull, but not in vertex positions (mid-edge or
plane \ center for instance) :
here you could use Contexts[ ] first to find the full path name.\ \>",
"Text"],
Cell["DiscreteMath`ConvexHull3D`Private`badones", "Input"],
Cell[BoxData[
\(Union[Flatten[ch3D]]\)], "Input"],
Cell["and these are the indices of the points inside the hull:",
"Text"],
Cell[BoxData[
\(Complement[Range[Length[points]], %]\)], "Input"],
Cell[BoxData[
\(Length/@{%%, %}\)], "Input"],
Cell[BoxData[
\(HullVolume[points, ch3D]\)], "Input"],
Cell[BoxData[
\(HullArea[points, ch3D]\)], "Input"],
Cell[BoxData[
\(% // N\)], "Input"],
Cell[BoxData[
\(still = HullShow[points, ch3D]\)], "Input"],
Cell[BoxData[
\(<< Graphics`Animation`\)], "Input"],
Cell[BoxData[
\(SpinShow[still]\)], "Input"],
Cell[BoxData[
\(\(?Hull*\)\)], "Input"]
}, Open ]]
},
FrontEndVersion->"Microsoft Windows 3.0", ScreenRectangle->{{0, 1024},
{0, 712}}, WindowSize->{527, 614},
WindowMargins->{{0, Automatic}, {Automatic, 5}} ]
(***********************************************************************
Cached data follows. If you edit this Notebook file directly, not
using Mathematica, you must remove the line containing CacheID at the
top of the file. The cache data will then be recreated when you save
this file from within Mathematica.
***********************************************************************)
(*CellTagsOutline
CellTagsIndex->{}
*)
(*CellTagsIndex
CellTagsIndex->{}
*)
(*NotebookFileOutline
Notebook[{
Cell[1709, 49, 174, 5, 70, "Input"], Cell[1886, 56, 260, 6, 90, "Text"],
Cell[2149, 64, 94, 2, 30, "Input"],
Cell[2246, 68, 50, 1, 30, "Input"],
Cell[CellGroupData[{
Cell[2321, 73, 42, 0, 43, "Subsubsection"], Cell[2366, 75, 140, 3, 52,
"Text"],
Cell[2509, 80, 68, 1, 30, "Input"],
Cell[2580, 83, 191, 4, 52, "Text"],
Cell[2774, 89, 58, 0, 30, "Input"],
Cell[2835, 91, 53, 1, 30, "Input"],
Cell[2891, 94, 72, 0, 33, "Text"],
Cell[2966, 96, 69, 1, 30, "Input"],
Cell[3038, 99, 48, 1, 30, "Input"],
Cell[3089, 102, 57, 1, 30, "Input"], Cell[3149, 105, 55, 1, 30,
"Input"], Cell[3207, 108, 39, 1, 30, "Input"], Cell[3249, 111, 63, 1,
30, "Input"], Cell[3315, 114, 55, 1, 30, "Input"], Cell[3373, 117, 48,
1, 30, "Input"], Cell[3424, 120, 43, 1, 30, "Input"]
}, Open ]]
}
]
*)
(***********************************************************************
End of Mathematica Notebook file.
***********************************************************************)
--=====================_893619951==_ Content-Type: text/plain;
charset="us-ascii"
(*
Convex Hull in 3D
*)
(*
Wouter Meeussen, 25/04/1998.
*)
(*
About
*)
(*
Very public domain, no rights to be reserved whatsoever. Feed back,
corrections, enhancements, speed-up and the like are appreciated at
:
Wouter Meeussen, <eu000949@pophost.eunet.be> Testing, suggestions and
encouragements by Russell Towle <rustybel@ foothill.net>.
Package design along Roman E. Maeder's "Template.nb". *)
(*
Limitations
*)
(*
Points should be unique. Identical points lead into trouble. If need be,
a routine should be written to eliminate them. 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..}..} ] 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
(*
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)}] (*
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[{10+ ang2[{cen,far},pts[[#]] ]-10. , norm[(far-pts[[#]])]
,#
}&/@idx ][[-2,-1]];
hull={xfar,xsec};
FixedPoint[Last[AppendTo[hull,
Sort[ {10+ang[ pts[[Take[hull,-2] ]] , pts[[#]] ]-10., 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[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];
*)
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..}..} ]:= 1/6* Plus@@ Apply[(pointlist[[#1]]
).Cross[pointlist[[#2]] ,pointlist[[#3]]
]&,
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..}..} ]:=
1/2*Plus@@Apply[norm[Cross[(pointlist[[#2]]-pointlist[[#1]]),(pointlist[[#3]]-
pointlist[[#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..}..} ]:= Show[Graphics3D[{PointSize[0.04],Map[Point,
pointlist,{1}]}],
Graphics3D[Polygon/@(pointlist[[#]]& /@hull)],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[ ]