Computational Geometry - Maple Programming Help

Online Help

All Products    Maple    MapleSim


Home : Support : Online Help : System : Information : Updates : Maple 2019 : updates/Maple2019/ComputationalGeometry

Computational Geometry

Maple 2019 includes a number of new commands in the ComputationalGeometry package.

 

PointOrientation

PointInCircle

PointOnSegment

SegmentsIntersect

MultiSegmentIntersect

ClosestPointPair

DelaunayTriangulation

PointOrientation

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

(1.1)

PointOrientation(a, b, d);

cw

(1.2)

PointOrientation(a, b, e);

degenerate

(1.3)

 

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 );

0.00.01.01.00.40.80.80.40.50.5

(1.4)

PointOrientation( A, [1,2,5]);

degenerate

(1.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);

ccw

(1.6)

PointOrientation(a, b, c, e);

cw

(1.7)

PointOrientation(a, b, c, f);

degenerate

(1.8)

PointInCircle

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

(2.1)

PointInCircle(a, b, c, e);

inside

(2.2)

PointInCircle(a, b, c, f);

boundary

(2.3)

PointOnSegment

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

(3.1)

PointOnSegment(d, [a, b]);

false

(3.2)

PointOnSegment(e, [a, b]);

false

(3.3)
 

SegmentsIntersect

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]]);

true

(4.1)

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) };

false

(4.2)

All of the green segments intersect L.

{ seq(SegmentsIntersect(L, K[i]), i=8..14) };

true

(4.3)

 

MultiSegmentIntersect

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" )

false

(5.1)

The fourth segment in L intersects M just once.

MultiSegmentIntersect( [L[4],  M ], mode="polygon" )

1.000000000,1.000000000

(5.2)

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

(5.3)

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.

(5.4)

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" );

segints0.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

(5.5)

 

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:

MultiSegmentIntersect:

0.398919953613456,0.801080046115181,−0.795771836740542,−0.404228163227617,0.673570378992075,−0.526429620902329,0.711020994207761,−0.530310867474975

(5.6)

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:

S1,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

MultiSegmentIntersect:

MultiSegmentIntersect:

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

(5.7)

ClosestPointPair

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

(6.1)

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:

ClosestPointPair:

ClosestPointPair:

ClosestPointPair:

ClosestPointPair:

ClosestPointPair:

ClosestPointPair:

ClosestPointPair:

ClosestPointPair:

ClosestPointPair:

0.3535533906,18,26

(6.1.1)


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

(6.2)

(d,p) := ClosestPointPair( L ); L[p[1]],L[p[2]];

d0.1732050808

p2150,4914

1.75,1.75,1.75,1.85,1.85,1.85

(6.3)

DelaunayTriangulation

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

(7.1)

t := DelaunayTriangulation(m);

t13,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

(7.2)

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);