New Features in Maple 2018 - Computational Geometry - Maplesoft

What's New in Maple 2018

Computational Geometry

Maple 2018 contains a new Computational Geometry package, based on Qhull (http://www.qhull.org). It contains functionality for programmers who intend to apply computational methods to polygons and clouds of points.  Computational geometry problems occur in many applications involving points in two- or higher-dimensional spaces, such as feature recognition, predicting vapor-liquid phase diagrams, delineating closely related regions for scattered data, and more..



with(ComputationalGeometry) 

[ConvexHull, DelaunayTriangulation, PolygonTriangulation, VoronoiDiagram]
 

The VoronoiDiagram command computes and plots the Voronoi diagram of a set of input points. 

m := LinearAlgebra:-RandomMatrix(40, 2); -1 

VoronoiDiagram(m, showpoints, symbol = solidcircle, symbolsize = 7, colorregions = ColorTools:-GetPalette( 

Plot_2d
 

 

The ConvexHull command can also compute with sets of points in higher dimensions. Here we pick some random points in 3-D, where the x-coordinates, y-coordinates, and z-coordinates are chosen according to different probability distributions: the x-coordinates are the absolute values of samples from Student's t-distribution; the y-coordinates are sampled from a standard normal distribution; and the z-coordinates come from a geometric distribution. 

n := 500; -1 

s := Statistics:-Sample; -1 

points := `<|>`(`~`[abs](`^`(s(StudentT(1), n), %T)), `^`(s(Normal(0, 1), n), %T), `^`(s(Geometric(`/`(1, 10)), n), %T)) 

Typesetting:-mprintslash([points := RTABLE(18446744074101631086, float[8], Matrix, rectangular, Fortran_order, [], 2, 1 .. 500, 1 .. 3)], [Matrix(%id = 18446744074101631086)])
 

hull := ConvexHull(points); -1; hull := map(proc (face) options operator, arrow; points[face] end proc, hull); -1 

plots:-display(plots:-pointplot3d(points, color = black), map(plottools:-polygon, hull, style = polygon), axes = none, transparency = .2) 

Plot_2d
 

The following example uses polygon data from a section of the built-in world map. This finds all land areas that intersect with the given viewport. 

m := DataSets:-Builtin:-WorldMap(); -1 

SetView(m, -150, 19, -50, 110); -1 

d := Display(m, view = [DEFAULT, DEFAULT]); -1; d; 1 

Plot_2d
 

data := plottools:-getdata(d); -1 

polygons := map(proc (x) options operator, arrow; x[3] end proc, [data]); -1 

numelems(polygons); 1 

41
 

There are 41 polygons representing various elements of this map. Such polygons will be displayed throughout this example, using the procedure in the Code Edit Region below, along with additional code to help with colors. 

 

Embedded component

 

displayPolygons(polygons); 1 

Plot_2d
 

Remove the polygon representing the background (and the oceans). 

polygons := select(proc (p) options operator, arrow; `<`(4, upperbound(p, 1)) end proc, polygons); -1 

numelems(polygons); 1 

40
 

In this example, every other polygon has more than four vertices, so that is an easy criterion. 

colors := randomColors(numelems(polygons)); -1 

pd := displayPolygons(polygons, colors); -1; pd; 1 

Plot_2d
 

This first demonstration uses the ConvexHull command. 

ConvexHull(polygons[1]); 1 

[192, 147, 146, 115, 112, 111, 86, 85, 84, 82, 74, 12, 11, 268, 267, 228, 217, 216, 215, 195]
 

It returns the indices of the points occurring in a convex hull. You can easily show these convex hulls as follows: 

convexHulls := map(proc (p) options operator, arrow; p[ConvexHull(p)] end proc, polygons); -1 

hd := displayPolygons(convexHulls, colors, transparency = .3); -1; hd; 1 

Plot_2d
 

plots:-display(pd, hd); 1 

Plot_2d
 

For the rest of this document, you will work with a single polygon. Use the one for which the area of the convex hull is largest: 

mainlandCanada := sort(polygons, key = (proc (p) options operator, arrow; ConvexHull(p, volume) end proc))[-1] 

Typesetting:-mprintslash([mainlandCanada := RTABLE(18446744074141058582, float[8], Matrix, rectangular, Fortran_order, [], 2, 1 .. 272, 1 .. 2)], [Matrix(%id = 18446744074141058582)])
 

The DelaunayTriangulation of a set of points subdivides the convex hull of these points into triangles in a way that avoids long and thin triangles as much as possible. 

tris := DelaunayTriangulation(mainlandCanada); -1 

dt := displayPolygons(map(proc (t) options operator, arrow; mainlandCanada[t] end proc, tris), style = polygon); -1; dt; 1 

Plot_2d
 

This overlay shows that some of the triangles fall inside the given polygon, and some fall outside. This is because the input is interpreted as a point cloud, not an ordered polygon. 

plots:-display(plottools:-polygon(mainlandCanada, transparency = .6, color = white), dt); 1 

Plot_2d
 

A VoronoiDiagram is a diagram that shows the dual of a Delaunay triangulation. 

VoronoiDiagram(mainlandCanada); 1 

Plot_2d
 

plots:-display(VoronoiDiagram(mainlandCanada, colorregions = false, transparency = .75), plots:-pointplot(mainlandCanada), dt); 1 

Plot_2d
 

The PolygonTriangulation command interprets its argument as a polygon, that is, an ordered sequence of points, in contrast with the other commands mentioned so far, which each interprets its argument as an unordered set of points. It subdivides this polygon into triangles, trying to avoid triangles that are long and narrow. 

triangulation := PolygonTriangulation(mainlandCanada); -1 

displayPolygons(map(proc (t) options operator, arrow; mainlandCanada[t] end proc, triangulation), style = polygon); 1 

Plot_2d