No Subject
- To: mathgroup@smc.vnet.net
- From: Wouter Meeussen <w.meeussen.vdmcc@vandemoortele.be>
- Date: Tue, 31 Mar 1998 02:28:44 -0500
hi,
in the add-on DiscreteMath`ComputationalGeometry` lives the program
ConvexHull2D. It finds the points on the perimeter (fence posts) of a
given set of points.
For those who would like to have ConvexHull3D, finding the bounding
polyhedron in 3D, here it comes. It gives surface and volume too (as a
side-kick).
(* By the way, the Area[Dodecahedron] in Geometry`Polytopes` is wrong:
it is given as a single pentagon, in stead of twelve of'm. *)
Hoping not to infuriating anyone by sending 17 Kb as attachment,
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[ 12981, 401]*)
(*NotebookOutlinePosition[ 13685, 426]*) (*
CellTagsIndexPosition[ 13641, 422]*) (*WindowFrame->Normal*)
Notebook[{
Cell[CellGroupData[{
Cell["Convex Hull in 3D", "Subtitle"],
Cell["\<\
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\ \>", "Text"],
Cell[CellGroupData[{
Cell["Initialisations", "Subsection"],
Cell[BoxData[
\(myOrdering[z_List] :=
\((Sort[Transpose[{z, Range[Length[z]]}]] // Transpose)\)\
[\([2]\)]\)],
"Input",
InitializationCell->True],
Cell["norm[vec_List]:=Sqrt[vec.vec]", "Input",
InitializationCell->True],
Cell["bary[pts_]:=1/Length[pts] Plus@@pts", "Input",
InitializationCell->True],
Cell["\<\
For points \"vi\" in 3D space, given as vi : {Real_,Real_,Real_}, the \
following function defines a plane :\ \>", "Text"],
Cell[BoxData[
\(plane[{v1_, v2_, v3_}, v_] :=
Det[{\((v - v1)\), \((v1 - v2)\), \ \((v2 - v3)\)}]\)], "Input",
InitializationCell->True],
Cell["\<\
this provides the angle between (v2-v1) and (v3-v2), usefull only in
case of \ 4 or more coplanar points:\
\>", "Text"],
Cell["\<\
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]\
\>", "Input",
InitializationCell->True],
Cell["\<\
ch2D[indices_List]:=Module[{idx,cen,pts,xfar,far,hull},
pts=points[[indices]];cen=bary[pts];dis=norm[#-cen]& /@ \
pts;idx=myOrdering[dis];xfar=Last[idx];far=pts[[xfar]];xsec=idx[[myOrdering[\
ang[{cen,far},pts[[#]] ]&/@idx][[-2]] ]]; hull={xfar,xsec};
FixedPoint[Last[AppendTo[hull,idx[[myOrdering[ang[pts[[Take[hull,-2] \
]],pts[[#]] ]&/@idx]
[[-3]] ]] ] ]&,0,SameTest->(First[hull]==#2
&)];indices[[Drop[hull,-1]]]]\ \>", "Input",
InitializationCell->True],
Cell["\<\
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.\
\>", "Text"],
Cell[BoxData[
\(test[tria : {_Integer, _Integer, _Integer}] := \n
If[Unequal@@tria, \t\n
Block[{i = 1, oneside = True}, res = {}; \n\t
While[\((oneside = \ Or[FreeQ[res, \(-1\)], FreeQ[res, 1]])\)
&&
i <= n\ , \n\t\t\t\t
res = {res, \n\t\t\t\t\t
If[FreeQ[tria, i],
Sign@\(Chop@plane[points[\([tria]\)], points[\([i]\)]\
]\),
0]\n\t\t\t\t\t\t}; \(i++\)]; \noneside], False]\)],
"Input",\
InitializationCell->True],
Cell[BoxData[
\(Remove[ConvexHull3D]\)], "Input"],
Cell[BoxData[
RowBox[{
RowBox[{
\(ConvexHull3D[points : {{_?NumericQ, _?NumericQ,
_?NumericQ}..}]\),
":=", "\n", "\t",
RowBox[{"Module", "[",
RowBox[{
\({distx, far, crosx, nix, it, faces, legs, cen}\), ",", "\n",
RowBox[{
\(n = Length[points]\), ";", "\n", \(cen = Plus@@points\
/n\),
";", "\n", \(cenpoints = \(\((# - cen)\)&\)/@points\),
";",
"\n", \(distx =
Reverse[myOrdering[\(\((#.#)\)&\)/@cenpoints]]\), ";",
"\n",
\(far = distx[\([1]\)]\), ";", "\n",
\(crosx =
Reverse[
myOrdering[
\((\(Abs[Cross[points[\([far]\)] - cen,
#]]&\)/@cenpoints)
\)*\n\t\t\((\ \(\((#.#)\)&\)/@cenpoints\
)\)\t]]\),
";", "\n", " ", \(nix = True\), ";", \(i = 1\), ";", "\n",
StyleBox[
\(While[nix && i <= n - 1, j = 1; \n\t\t\t\t
While[\ j <= n - 1\ && \n\t\t\t
\((nix =
Not[test[
it = {far, crosx[\([i]\)],
crosx[\([j]\)]}]])
\), \n\t\t\(j++\)]; \n\t\(i++\)]\),
CellFrame->True,
FontColor->RGBColor[1, 0, 0],
Background->None],
StyleBox[";",
CellFrame->True,
FontColor->RGBColor[1, 0, 0],
Background->None], "\n",
\(If[Count[res, 0, \(-1\)] > 3,
it = ch2D[Flatten[Position[Flatten[res], 0]]\ \ ]\ \
]\), ";",
"\n", \(faces = {it}\), ";", "\n", \(legs = {}\), ";",
"\n",
\(faces = Union@Flatten[{faces, {it}}, 1]\), ";", "\n",
\(legs =
Split[Sort[
Sort/@Flatten[
\(Partition[#~Append~First[#], 2, 1]&\)/@faces,
1]]]
\), ";", "\n",
\(legs = Flatten[DeleteCases[legs, {_List, _List}], 1]\),
";",
"\n", \(While[legs =!= {}, \n
faces = Union@Flatten[{faces, {it}}, 1]; \n
legs = Split[
Sort[Sort/@
Flatten[
\(Partition[#~Append~First[#], 2,
1]&\)/@faces,
1]]]; \n
legs = Flatten[DeleteCases[legs, {q_List, q_List..}],
1]; \n
testing =
First@Select[faces, \
\((\(! FreeQ[#, \ legs[\([1, 1]\)]\ ]\) && \
\(! FreeQ[#, legs[\([1, 2]\)]\ ]\))\)&]; \n
nix = True; \nj = 1; \
testable = \tDeleteCases[crosx, Alternatives@@testing];
\n
\t\t\t\tWhile[\
j <= Length[testable]\ && \n\t\t\t
\((nix =
Not[test[
it = Flatten[{legs[\([1]\)],
\n\t\t\t\t\t\t\t
testable[\([j]\)]}]\ ]\ ])\), \n\t\t
\(j++\)]; \n
If[Count[res, 0, \(-1\)] > 3,
it = ch2D[Flatten[Position[Flatten[res], 0]]\ \ ]\ \
]; \n
\t\t (*\t\ Print[{testing, legs\ , it}]\ ; \ *) \n\t
faces = Union@Flatten[{faces, {it}}, 1]; \n
legs = Split[
Sort[Sort/@
Flatten[
\(Partition[#~Append~First[#], 2,
1]&\)/@faces,
1]]]; \n
legs = Flatten[DeleteCases[legs, {q_List, q_List..}],
1]; \ ]
\), ";", "\n",
\(circulation =
\(Sign[Chop[
plane[Take[#, 3] /. k_Integer :>
points[\([k]\)],
cen]]]&\)/@faces\), ";", "\t\t", "\n",
\(faces =
\(MapThread[List, {faces, circulation}] /. {poly_List,
1} :>
Reverse[poly]\) /. {poly_List, \(-1\)} :> poly\),
";",
"\n", \(Sort[
\(FixedPoint[RotateLeft, #,
SameTest -> \((First[#2] ==
Min[#2]&)\)]&\)/@faces]
\)}]}], "\n", " ", "]"}]}], "\n", "\t\t"}]],
"Input",
InitializationCell->True]
}, Closed]],
Cell[CellGroupData[{
Cell["Example", "Subsection"],
Cell["\<\
Create a random set of points :
(* about 100 random points would take approx 60 sec on a pentium 90 MHz,
\ 32Mbyte Ram, winNT3.51 *)\
\>", "Text"],
Cell[BoxData[
\(n = 100; \n
points = \ \
Table[Plus@@Table[Random[Real, {\(-1\), 1}], {8}], {n}, {3}]; \)],
"Input"],
Cell["or get them from a regular solid :", "Text"],
Cell[BoxData[
\(<< Geometry`Polytopes`\)], "Input"],
Cell[BoxData[
\(points = Vertices[Dodecahedron]; \ n = Length[points]\)],
"Input"],
Cell["\<\
or even from a truncated, stellated, or other form from the \
Graphics`Polyhedra` add-on,\
\>", "Text"],
Cell[BoxData[
\(<< Graphics`Polyhedra`\)], "Input"],
Cell["\<\
points=Union[Chop[(First[Truncate[ Polyhedron[Dodecahedron], .4]]
/.Polygon->Sequence)~Flatten~1 ],SameTest->(Max@@(#1-#2)< 10^-7&)];\
\>", "Input"],
Cell[BoxData[
\(n = Length[points]\)], "Input"],
Cell["the number of planes (in absence of coplanarity) is:", "Text"],
Cell[BoxData[
\(Binomial[n, 3]\)], "Input"],
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["this are the indices of the points on the convex hull:", "Text"],
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"]
}, Closed]],
Cell[CellGroupData[{
Cell["Testing", "Subsection"],
Cell["\<\
here we check the circulations around each polygon : should be outward \
pointing :\
\>", "Text"],
Cell["\<\
cen=bary[points]
circulation=Sign[Chop[plane[Take[#,3]/.k_Integer:>points[[k]],cen]]]&/@ch3D\
\>", "Input"],
Cell["\<\
If we want volume and surface, then polygons larger than triangles
should be \ devided into coplanar triangles :\
\>", "Text"],
Cell["\<\
pieces=Flatten[ Thread[ z[First[#],Partition[Rest[#],2,1] ]] &/@ ch3D \
,1]/.z[q__]:>Flatten[{q}]\
\>", "Input"],
Cell["\<\
now we find the volume, either measured from the barycenter (all \
contributions positive)\
\>", "Text"],
Cell["\<\
volume=1/6* Plus@@ \
Apply[(points[[#1]]-cen).Cross[points[[#2]]-cen,points[[#3]]-cen]&,pieces,1]
\
\>", "Input"],
Cell["\<\
or from the coordinate zero, with some positive and some negative terms,
but it comes to the same result :\
\>", "Text"],
Cell["\<\
volume= 1/6* Plus@@ Apply[(points[[#1]] ).Cross[points[[#2]]
,points[[#3]] \ ]&,pieces,1]\
\>", "Input"],
Cell["for the surface (or area), we do not need the barycenter :",
"Text"],
Cell["\<\
surface=1/2*Plus@@Apply[norm[Cross[(points[[#2]]-points[[#1]]),(points[[#3]]-\
points[[#2]])]]&,pieces,1]\
\>", "Input"],
Cell[TextData[StyleBox[
"If we choose the dodecahedron for testing volume & area, we should take
the \ length of our sides into account :", "Text"]], "Text"],
Cell["\<\
?Volume
Volume[Dodecahedron]
volume==(%//N)*norm[points[[1]]-points[[2]]]^3\ \>", "Input"],
Cell["Area[Dodecahedron]/Area[Pentagon]//N", "Input"],
Cell["For the area, we find : ", "Text"],
Cell["\<\
?Area
Area[Dodecahedron]
surface== 12 Area[Dodecahedron]*norm[points[[1]]-points[[2]]]^2\
\>", "Input"],
Cell["\<\
A spot of trouble here : the Area[Dodecahedron] is in error since it is
the \ same as the area of a regular pentagon, in stead of 12 times that
area. (Oops).\ \>", "Text"]
}, Closed]],
Cell[CellGroupData[{
Cell["Dispaying", "Subsection"],
Cell[BoxData[
\(\(polis = Polygon/@\((\(points[\([#]\)]&\)\ /@\ ch3D)\); \)\)],
"Input"],
Cell[BoxData[
RowBox[{"still", "=",
RowBox[{"Show", "[",
RowBox[{
\(Graphics3D[{PointSize[0.05], Map[Point, \ points, {1}]}]\),
",",
\(Graphics3D[polis]\), ",", \(Axes -> True\), ",",
RowBox[{
StyleBox["AspectRatio",
"MR"], " ",
StyleBox["->",
"MR"], " ",
StyleBox["Automatic",
"MR"]}], " ", ",", \(BoxRatios -> {1, 1, 1}\)}], "]"}]}]],
"Input"],
Cell["this shows the polygon in rotation :", "Text"],
Cell[BoxData[
\(<< Graphics`Animation`\)], "Input"],
Cell[BoxData[
\(SpinShow[still]\)], "Input"]
}, Closed]]
}, Open ]]
},
FrontEndVersion->"Microsoft Windows 3.0", ScreenRectangle->{{0, 1024},
{0, 740}}, AutoGeneratedPackage->None,
WindowToolbars->"EditBar",
WindowSize->{726, 559},
WindowMargins->{{0, Automatic}, {Automatic, 6}} ]
(***********************************************************************
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[CellGroupData[{
Cell[1731, 51, 37, 0, 54, "Subtitle"], Cell[1771, 53, 206, 5, 71,
"Text"],
Cell[CellGroupData[{
Cell[2002, 62, 37, 0, 47, "Subsection"], Cell[2042, 64, 161, 4, 30,
"Input",
InitializationCell->True],
Cell[2206, 70, 74, 1, 30, "Input",
InitializationCell->True],
Cell[2283, 73, 80, 1, 30, "Input",
InitializationCell->True],
Cell[2366, 76, 132, 3, 33, "Text"],
Cell[2501, 81, 149, 3, 30, "Input",
InitializationCell->True],
Cell[2653, 86, 129, 3, 33, "Text"],
Cell[2785, 91, 159, 4, 48, "Input",
InitializationCell->True],
Cell[2947, 97, 470, 10, 120, "Input",
InitializationCell->True],
Cell[3420, 109, 276, 5, 52, "Text"], Cell[3699, 116, 524, 11, 170,
"Input",
InitializationCell->True],
Cell[4226, 129, 53, 1, 30, "Input"], Cell[4282, 132, 4318, 88, 870,
"Input",
InitializationCell->True]
}, Closed]],
Cell[CellGroupData[{
Cell[8637, 225, 29, 0, 31, "Subsection"], Cell[8669, 227, 156, 4, 52,
"Text"], Cell[8828, 233, 134, 4, 50, "Input"], Cell[8965, 239, 50, 0,
33, "Text"],
Cell[9018, 241, 55, 1, 30, "Input"], Cell[9076, 244, 86, 1, 30,
"Input"], Cell[9165, 247, 113, 3, 33, "Text"], Cell[9281, 252, 55, 1,
30, "Input"], Cell[9339, 255, 159, 3, 48, "Input"], Cell[9501, 260, 51,
1, 30, "Input"], Cell[9555, 263, 68, 0, 33, "Text"],
Cell[9626, 265, 47, 1, 30, "Input"], Cell[9676, 268, 140, 3, 33,
"Text"], Cell[9819, 273, 68, 1, 30, "Input"], Cell[9890, 276, 70, 0,
33, "Text"],
Cell[9963, 278, 53, 1, 30, "Input"], Cell[10019, 281, 72, 0, 33,
"Text"], Cell[10094, 283, 69, 1, 30, "Input"], Cell[10166, 286, 48, 1,
30, "Input"] }, Closed]],
Cell[CellGroupData[{
Cell[10251, 292, 29, 0, 31, "Subsection"], Cell[10283, 294, 107, 3, 33,
"Text"], Cell[10393, 299, 117, 3, 48, "Input"], Cell[10513, 304, 136,
3, 33, "Text"], Cell[10652, 309, 122, 3, 48, "Input"], Cell[10777, 314,
114, 3, 33, "Text"], Cell[10894, 319, 124, 4, 66, "Input"], Cell[11021,
325, 130, 3, 52, "Text"], Cell[11154, 330, 115, 3, 30, "Input"],
Cell[11272, 335, 74, 0, 33, "Text"], Cell[11349, 337, 130, 3, 48,
"Input"], Cell[11482, 342, 157, 2, 33, "Text"], Cell[11642, 346, 100,
4, 66, "Input"], Cell[11745, 352, 53, 0, 30, "Input"], Cell[11801, 354,
40, 0, 33, "Text"], Cell[11844, 356, 117, 4, 66, "Input"], Cell[11964,
362, 181, 4, 52, "Text"] }, Closed]],
Cell[CellGroupData[{
Cell[12182, 371, 31, 0, 31, "Subsection"], Cell[12216, 373, 92, 1, 30,
"Input"], Cell[12311, 376, 478, 13, 50, "Input"], Cell[12792, 391, 52,
0, 33, "Text"], Cell[12847, 393, 55, 1, 30, "Input"], Cell[12905, 396,
48, 1, 30, "Input"] }, Closed]]
}, Open ]]
}
]
*)
(***********************************************************************
End of Mathematica Notebook file.
***********************************************************************)
NV Vandemoortele Coordination Center Oils & Fats Applied Research
Prins Albertlaan 79
Postbus 40
B-8870 Izegem (Belgium)
Tel: +/32/51/33 21 11
Fax: +/32/51/33 21 75
vdmcc@vandemoortele.be