Computational Geometry - New Features in Maple 2019 - Maplesoft

What's New in Maple 2019

Computational Geometry

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



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
 

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

rtable(1 .. 5, 1 .. 2, [[0., 0.], [1.0, 1.0], [.4, .8], [.8, .4], [.5, .5]], subtype = Matrix)
 

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

degenerate
 

> 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
 

> PointOrientation(a, b, c, e);
 

cw
 

> PointOrientation(a, b, c, f);
 

degenerate
 

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

Plot_2d
 

> PointInCircle(a, b, c, d);
 

outside
 

> PointInCircle(a, b, c, e);
 

inside
 

> PointInCircle(a, b, c, f);
 

boundary
 

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

Plot_2d
 

> PointOnSegment(c, [a, b]);
 

true
 

> PointOnSegment(d, [a, b]);
 

false
 

> PointOnSegment(e, [a, b]);
 

false
 

 

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
 

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

Plot_2d
 

None of the blue segments intersect L (orange). 

> { seq(SegmentsIntersect(L, K[i]), i=1..7) };
 

{false}
 

All of the green segments intersect L. 

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

{true}
 

 

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

Plot_2d
 

> MultiSegmentIntersect( [L[1],  M ], mode="polygon", output="truefalse" )
 

false
 

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 )
 

[[.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"] )
 

[[.5000000000, 0.], [.5000000000, 1.000000000], [1.000000000, 1.000000000], [0, 0], [1, 0], [.2, 1], [.4, 1], [0., 1.000000000], [0., 0.], [0., .3750000000], [1.000000000, 1.000000000], [1.000000000, ...
[[.5000000000, 0.], [.5000000000, 1.000000000], [1.000000000, 1.000000000], [0, 0], [1, 0], [.2, 1], [.4, 1], [0., 1.000000000], [0., 0.], [0., .3750000000], [1.000000000, 1.000000000], [1.000000000, ...
 

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

Plot_2d
 

 

> segints := MultiSegmentIntersect( [ P, Q, R], mode="polygon" );
 

Typesetting:-mprintslash([segints := [[.748398218595573, .451601781343994], [.398919953613456, .801080046115181], [0.318673784065071e-1, .848528137100000], [-.136783930483585, .848528137100000], [-.35...
Typesetting:-mprintslash([segints := [[.748398218595573, .451601781343994], [.398919953613456, .801080046115181], [0.318673784065071e-1, .848528137100000], [-.136783930483585, .848528137100000], [-.35...
Typesetting:-mprintslash([segints := [[.748398218595573, .451601781343994], [.398919953613456, .801080046115181], [0.318673784065071e-1, .848528137100000], [-.136783930483585, .848528137100000], [-.35...
Typesetting:-mprintslash([segints := [[.748398218595573, .451601781343994], [.398919953613456, .801080046115181], [0.318673784065071e-1, .848528137100000], [-.136783930483585, .848528137100000], [-.35...
Typesetting:-mprintslash([segints := [[.748398218595573, .451601781343994], [.398919953613456, .801080046115181], [0.318673784065071e-1, .848528137100000], [-.136783930483585, .848528137100000], [-.35...
Typesetting:-mprintslash([segints := [[.748398218595573, .451601781343994], [.398919953613456, .801080046115181], [0.318673784065071e-1, .848528137100000], [-.136783930483585, .848528137100000], [-.35...
 

 

> plots:-display( theplot,
plottools:-point( segints,
              symbol="circle", 'symbolsize'=15, 'color'="Black"), axes=none );
 

Plot_2d
 

 

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:
Plot_2d
MultiSegmentIntersect:
Plot_2d
Typesetting:-mprintslash([[[.398919953613456, .801080046115181], [-.795771836740542, -.404228163227617], [.673570378992075, -.526429620902329], [.711020994207761, -.530310867474975]]], [[[HFloat(0.398...
Typesetting:-mprintslash([[[.398919953613456, .801080046115181], [-.795771836740542, -.404228163227617], [.673570378992075, -.526429620902329], [.711020994207761, -.530310867474975]]], [[[HFloat(0.398...
 

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:
 

 

 

 

 

 

Typesetting:-mprintslash([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], [...
Typesetting:-mprintslash([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], [...
MultiSegmentIntersect:
Plot_2d
MultiSegmentIntersect:
Plot_2d
Typesetting:-mprintslash([[[.748398218595573, .451601781343994], [.398919953613456, .801080046115181], [0.318673784065071e-1, .848528137100000], [-.136783930483585, .848528137100000], [-.3513148188126...
Typesetting:-mprintslash([[[.748398218595573, .451601781343994], [.398919953613456, .801080046115181], [0.318673784065071e-1, .848528137100000], [-.136783930483585, .848528137100000], [-.3513148188126...
Typesetting:-mprintslash([[[.748398218595573, .451601781343994], [.398919953613456, .801080046115181], [0.318673784065071e-1, .848528137100000], [-.136783930483585, .848528137100000], [-.3513148188126...
Typesetting:-mprintslash([[[.748398218595573, .451601781343994], [.398919953613456, .801080046115181], [0.318673784065071e-1, .848528137100000], [-.136783930483585, .848528137100000], [-.3513148188126...
Typesetting:-mprintslash([[[.748398218595573, .451601781343994], [.398919953613456, .801080046115181], [0.318673784065071e-1, .848528137100000], [-.136783930483585, .848528137100000], [-.3513148188126...
Typesetting:-mprintslash([[[.748398218595573, .451601781343994], [.398919953613456, .801080046115181], [0.318673784065071e-1, .848528137100000], [-.136783930483585, .848528137100000], [-.3513148188126...
 

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]
 

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:
Plot_2d

 

ClosestPointPair:

Plot_2d

 

ClosestPointPair:

Plot_2d

 

ClosestPointPair:

Plot_2d

 

ClosestPointPair:

Plot_2d

 

ClosestPointPair:

Plot_2d

 

ClosestPointPair:

Plot_2d

 

ClosestPointPair:

Plot_2d

 

ClosestPointPair:

Plot_2d

 

ClosestPointPair:

Plot_2d
.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]];
 

 

 

Typesetting:-mprintslash([d := .1732050808], [.1732050808])
Typesetting:-mprintslash([p := [2150, 4914]], [[2150, 4914]])
[1.75, 1.75, 1.75], [1.85, 1.85, 1.85]
 

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

Typesetting:-mprintslash([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], ...
Typesetting:-mprintslash([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], ...
Typesetting:-mprintslash([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], ...
Typesetting:-mprintslash([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], ...
Typesetting:-mprintslash([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], ...
 

> t := DelaunayTriangulation(m);
 

Typesetting:-mprintslash([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,...
Typesetting:-mprintslash([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,...
Typesetting:-mprintslash([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,...
Typesetting:-mprintslash([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,...
 

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

Plot_2d
 



Bereit für den nächsten Schritt?

Sprechen Sie mit unseren Produktspezialisten über eine kostenlose Demoversion von Maple

*Die Maple-Evaluation ist für Schüler und Studenten bzw. die private Nutzung zurzeit nicht verfügbar.