Computational Geometry
Maple 2019 includes a number of new commands in the ComputationalGeometry package.
PointOrientation
PointInCircle
PointOnSegment
SegmentsIntersect
MultiSegmentIntersect
ClosestPointPair
DelaunayTriangulation
This command gives a very efficient test of the orientation of three points in two dimensions or four points in three dimensions. Either "cw" or "ccw" for points in clockwise or counterclockwise order or "degenerate" if they are colinear or coplanar.
with(ComputationalGeometry):
a := [0,0]: b := [1,1]: c := [0.4,0.8]: d := [0.8,0.4]: e:=[0.5,0.5]:
PointOrientation(a, b, c);
ccw
PointOrientation(a, b, d);
cw
PointOrientation(a, b, e);
degenerate
The most efficient way to use this command is on a large array of point coordinates. PointOrientation works on these in place without having to copy the points into a new memory location making it efficient when doing many of these calculations at a time.
A := Array( [a,b,c,d,e], datatype=float[8], order=C_order );
A≔0.0.1.1.0.4000000000000000.8000000000000000.8000000000000000.4000000000000000.5000000000000000.500000000000000
PointOrientation( A, [1,2,5]);
a := [0,0,0]: b := [1,1,0]: c := [0.4,0.8,0]: d := [0.8,0.4,-1]: e := [0.5,0.5,1]: f := [0.5,0.5,0]:
PointOrientation(a, b, c, d);
PointOrientation(a, b, c, e);
PointOrientation(a, b, c, f);
This command gives a very efficient test to determine if a fourth point is inside a circle defined by three points. Like PointOrientation, the PointInCircle command also has an efficient mode for arrays of points.
a := [0,0]: b := [1,1]: c := [0,1]: d := [0.5,1.25]: e:=[0.25,0.75]: f:=[1,0]:
plots:-display(plottools:-point([a,b,c,d,e,f],symbolsize=20), plots:-textplot([d[], "d"], align=["above", "left"]), plots:-textplot([e[], "e"], align=["above", "left"]), plots:-textplot([f[], "f"], align=["above", "left"]), plottools:-circle([1/2,1/2], sqrt(2)/2, style=line), axes=box );
PointInCircle(a, b, c, d);
outside
PointInCircle(a, b, c, e);
inside
PointInCircle(a, b, c, f);
boundary
This command gives an accurate test to determine if a third point is on a line segment defined by two points. Like PointOrientation, the PointOnSegment command also has an efficient mode for arrays of points.
a := [0,0]: b := [1,1]: c := [0.4,0.4]: d := [1.1,1.1]: e:=[0.5,0.55]:
plots:-display(plottools:-point([a,b,c,d,e],symbolsize=20), plots:-textplot([c[], "c"], align=["above", "left"]), plots:-textplot([d[], "d"], align=["below", "left"]), plots:-textplot([e[], "e"], align=["above", "left"]), plottools:-line(a, b), axes=box );
PointOnSegment(c, [a, b]);
true
PointOnSegment(d, [a, b]);
false
PointOnSegment(e, [a, b]);
This command gives an efficient test to determine if two line segments intersect. Like PointOrientation, the SegmentsIntersect command also has an efficient mode for arrays of points.
SegmentsIntersect([[0,0], [1,1]], [[0,1], [1,0]]);
L := [[1,0], [1,5]]:
K := [ [[1,-2], [1,-1]], [[0,1], [0,3]], [[2,1], [2,3]], [[1,6.5], [1,8]], [[1,6], [1.5,8]], [[0.5,4], [1.5,7]], [[0.5, 0.5], [1, -0.75]], [[1,5], [1.5,6]], [[1,4.5], [1.5,5]], [[.5,3.75], [1.5,3.5]], [[1,2.5], [1,3.5]], [[0.5,1.5], [1,2]], [[1,0], [1.5,-2]], [[1, 0.5], [1, -0.5]] ]:
plots:-display( plottools:-line(L[], legend="L", color="Dalton 2", thickness=3), seq(plottools:-line(K[i][], legend="K " || i, color="Blue"), i=1..7), seq(plottools:-line(K[i][], legend="K " || i, color="Green"), i=8..14), axes=boxed, style=pointline);
None of the blue segments intersect L (orange).
{ seq(SegmentsIntersect(L, K[i]), i=1..7) };
All of the green segments intersect L.
{ seq(SegmentsIntersect(L, K[i]), i=8..14) };
The MultiSegmentIntersect command uses a sweep line algorithm to determine if the collections of line segments have intersections. It can be used, to name a few cases, to test if a line segment intersects a polygon, test if a collection of segments has any intersections, or find all the intersections in a collection of polygons.
Here L is a collection of line segments, and M is a list of the vertices of a polygon.
L := [ [[-.5, -1], [-.5,2]], [[.5,-1], [.5,2]], [[1.5,-1], [1.5,2]], [[0.75,1.5], [1.25,0.5]], [[-.25,0], [1.25, 0]], [[.2,1], [.4,1]], [[-.25,.5], [.25,.25]], [[0, 0.5], [0, 1.5]] ]: M := [[0,0], [0,1], [1,1], [1,0]]:
plots:-display( seq(plottools:-line(L[i][], legend="L "||i, color="Red",style=pointline),i=1..nops(L)), plottools:-polygon(M,style=pointline, legend="M"), axes=none);
MultiSegmentIntersect( [L[1], M ], mode="polygon", output="truefalse" )
The fourth segment in L intersects M just once.
MultiSegmentIntersect( [L[4], M ], mode="polygon" )
1.000000000,1.000000000
The second segment in L intersects M twice, but we can use the 'lazy' option to find just one of them.
MultiSegmentIntersect( [L[2], M ], mode="polygon", lazy )
0.5000000000,1.000000000
We can also find all the intersections of L and M at once by flattening L into a list of points to be treated as point pairs.
MultiSegmentIntersect( [map(op,L), M ], mode=["pairs","polygon"] )
0.5000000000,0.,0.5000000000,1.000000000,1.000000000,1.000000000,0,0,1,0,0.2,1,0.4,1,0.,1.000000000,0.,0.,0.,0.3750000000,1.000000000,1.000000000,1.000000000,−0.
P, Q, and R are the vertices of three irregular polygons.
P := op(1,plottools:-polygonbyname("8gon", inbox=[1.2,1.2])): Q := op(1,plottools:-polygonbyname("3gon", irregular, star, radius=1/2)): R := op(1,plottools:-polygonbyname("4gon", irregular, star=3, radius=3/4)):
theplot := plots:-display( plottools:-polygon(P, style=line, color="Niagara 6"), plottools:-polygon(Q, style=line, color="Niagara 5"), plottools:-polygon(R, style=line, color="Niagara 9" ), axes=none );
segints := MultiSegmentIntersect( [ P, Q, R], mode="polygon" );
segints≔0.748398218595573,0.451601781343994,0.398919953613456,0.801080046115181,0.0318673784065071,0.848528137100000,−0.136783930483585,0.848528137100000,−0.351314818812679,0.848528137100000,−0.848528137400000,0.318862801207573,−0.783470937602132,−0.416529062358602,−0.795771836740542,−0.404228163227617,−0.805408563800242,−0.394591436173733,−0.310621779662489,−0.848528137100000,0.673570378992075,−0.526429620902329,0.754216232835860,−0.445783767107218,0.711020994207761,−0.530310867474975,0.777977217493130,−0.498309822607008,−0.760087396244286,−0.417676105328972,−0.793448721745016,−0.400683269457881
plots:-display( theplot, plottools:-point( segints, symbol="circle", 'symbolsize'=15, 'color'="Black"), axes=none );
As with several other new ComputationalGeometry commands setting the infolevel of CompGeomPlot to 1 or greater will give some graphical representation of the steps of the algorithm. In this example, we use the default mode="pairs" mode which treats the collections of points as unconnected line segments. An infolevel value of 1 plots the input and, at the end, the intersections found. Values of 2 or 4 will show the steps of the algorithm as it runs.
infolevel[CompGeomPlot] := 1: MultiSegmentIntersect( [ P, Q, R] );
MultiSegmentIntersect:
0.398919953613456,0.801080046115181,−0.795771836740542,−0.404228163227617,0.673570378992075,−0.526429620902329,0.711020994207761,−0.530310867474975
The mode option of MultiSegmentIntersect is meant to allow for convenience. The command can also take a single rtable of point coordinates together with a list of lists of index pairs specifying line segment endpoints in the point rtable.
M := Matrix([[P],[Q],[R]]): S := [ [seq([i,i+1], i=1..8-1), [8,1]], [seq([8+i,8+i+1], i=1..6-1), [8+6,8+1]], [seq([14+i,14+i+1], i=1..8-1), [14+8,14+1]] ]; infolevel[CompGeomPlot] := 1: MultiSegmentIntersect(M, S); infolevel[CompGeomPlot] := 0:
S≔1,2,2,3,3,4,4,5,5,6,6,7,7,8,8,1,9,10,10,11,11,12,12,13,13,14,14,9,15,16,16,17,17,18,18,19,19,20,20,21,21,22,22,15
0.748398218595573,0.451601781343994,0.398919953613456,0.801080046115181,0.0318673784065071,0.848528137100000,−0.136783930483585,0.848528137100000,−0.351314818812679,0.848528137100000,−0.848528137400000,0.318862801207573,−0.783470937602132,−0.416529062358602,−0.795771836740542,−0.404228163227617,−0.805408563800242,−0.394591436173733,−0.310621779662489,−0.848528137100000,0.673570378992075,−0.526429620902329,0.754216232835860,−0.445783767107218,0.711020994207761,−0.530310867474975,0.777977217493130,−0.498309822607008,−0.760087396244286,−0.417676105328972,−0.793448721745016,−0.400683269457881
Use the ClosestPointPair command to find a closest pair of points in a collection of points in any dimension. It uses a classic divide and conquer approach.
ClosestPointPair( [[0,1], [2,2], [1,1], [4,4]] );
1,1,3
Verbose Demonstration
This command can be used to demonstrate the algorithm visually by setting the infolevel of CompGeomPlot to 4.
infolevel[CompGeomPlot]:=4: ClosestPointPair( [seq(seq(seq([i, j], i = 0 .. 2, .5), j = 0 .. 2, .5)), [1.25, 1.25]] ); infolevel[CompGeomPlot]:=1:
ClosestPointPair:
0.3535533906,18,26
Here is a slightly larger example in three dimensions.
L := [seq(seq(seq([i, j, k], i = 0 .. 4, .25), j = 0 .. 4, .25), k = 0 .. 4, .25), [1.85, 1.85, 1.85]]: numelems(L);
4914
(d,p) := ClosestPointPair( L ); L[p[1]],L[p[2]];
d≔0.1732050808
p≔2150,4914
1.75,1.75,1.75,1.85,1.85,1.85
The DelaunayTriangulation command introduced in Maple 2018 now supports computing triangularizations in dimensions up to 10.
m := RandomTools:-Generate(list(list(float(range=-95..95),3),15));
m≔−22.70990285,−84.26614411,−3.86641598,−4.49508400,57.52946750,−83.20308836,9.35541031,16.74985062,46.67376825,−66.65183214,30.97968824,−42.87742988,6.90092535,77.44513031,−53.51950332,−11.84027187,60.31680656,21.73755375,8.97960333,61.76051242,27.40775176,60.28912701,8.33727108,−74.22489097,27.09842774,−83.58721340,−0.78447362,−5.17789092,−34.43051994,13.07930916,15.60891286,−31.91703837,−65.32248447,22.19030206,22.85228936,78.14996416,−42.78547388,−56.49945055,79.84896994,79.69509190,−40.48209840,−93.75987033,−7.53357838,20.82785949,74.72960972
t := DelaunayTriangulation(m);
t≔13,10,1,4,10,11,1,4,10,9,13,1,11,9,1,14,11,9,10,1,6,11,10,4,6,11,7,10,6,5,11,4,5,6,11,7,11,8,7,10,8,9,11,14,9,8,11,10,8,5,11,7,9,12,10,13,8,12,9,14,5,2,11,4,2,8,5,11,2,8,11,14,15,10,13,4,15,6,10,4,12,3,10,13,3,15,10,13,15,3,12,13,3,15,12,7,8,3,7,10,3,12,8,7,3,6,7,10,3,15,6,10,15,3,6,7,3,8,9,10,12,3,9,10,12,3,8,9
We can collect the outputs to form the polygon faces of the tetrahedrons in t and plot them.
tetrahedron := A->plottools:-polygon([[A[1], A[2], A[3]], [A[1], A[2], A[4]], [A[1], A[3], A[4]], [A[2], A[3], A[4]]], style=line, color="Black"): plots:-display( plottools:-point(m, symbolsize=20, color="Red"), map(x->tetrahedron(m[x]),t),axes=none);
Download Help Document
