**Computational Geometry** is a field of mathematics that seeks the development of efficient algorithms to solve problems described in terms of basic geometrical objects. We differentiate between **Combinatorial Computational Geometry** and **Numerical Computational Geometry**.

Combinatorial Computational Geometry deals with interaction of basic geometrical objects: points, segments, lines, polygons, and polyhedra. In this setting we have three categories of problems:

**Static Problems**: The construction of a known target object is required from a set of input geometric objects.**Geometric Query Problems**: Given a set of known objects (the**search space**) and a sought property (the**query**), these problems deal with the search of objects that satisfy the query.**Dynamic Problems**: Similar to problems from the previous two categories, with the added challenge that the input is not known in advance, and objects are inserted or deleted between queries/constructions.

Numerical Computational Geometry deals mostly with representation of objects in space described by means of curves, surfaces, and regions in space bounded by those.

Before we proceed to the development and analysis of the different algorithms in those two settings, it pays off to explore the basic background: Plane Geometry.

The basic Geometry capabilities are usually treated through the Geometry module of the `sympy`

libraries. Rather than giving an academical description of all objects and properties in that module, we discover the most useful ones through a series of small self-explanatory `python`

sessions.

We start with the concepts of **point** and **segment**. The aim is to illustrate how easily we can check for collinearity, compute lengths, midpoints, or slopes of segments, for example. We also show how to quickly compute the angle between two segments, as well as deciding whether a given point belongs to a segment or not. The next diagram illustrates an example, which we follow up with code.

In [1]:

```
from sympy.geometry import *
P1 = Point(0, 0)
P2 = Point(3, 4)
P3 = Point(2, -1)
P4 = Point(-1, 5)
S1 = Segment(P1, P2)
S2 = Segment(P3, P4)
```

In [2]:

```
Point.is_collinear(P1, P2, P3)
```

Out[2]:

In [3]:

```
S1.length
```

Out[3]:

In [4]:

```
S2.midpoint
```

Out[4]:

In [5]:

```
S1.slope
```

Out[5]:

In [6]:

```
S1.intersection(S2)
```

Out[6]:

In [7]:

```
Segment.angle_between(S1, S2)
```

Out[7]:

In [8]:

```
S1.contains(P3)
```

Out[8]:

The next logical geometrical concept is the **line**. We can perform more interesting operations with lines, and to that effect we have a few more constructors. We can find their equations, compute the distance between a point and a line, and many other operations.

In [9]:

```
L1 = Line(P1, P2)
L2 = L1.perpendicular_line(P3) # perpendicular line to L1
```

In [10]:

```
L2.arbitrary_point() # parametric equation of L2
```

Out[10]:

In [11]:

```
L2.equation() # algebraic equation of L2
```

Out[11]:

In [12]:

```
L2.contains(P4) # is P4 in L2?
```

Out[12]:

In [13]:

```
L2.distance(P4) # distance from P4 to L2
```

Out[13]:

In [14]:

```
L1.is_parallel(S2) # is S2 parallel to L1?
```

Out[14]:

The next geometrical concept we are to explore is the **circle**. We may define a circle by its center and radius, or by three points on it. We can easily compute all of its properties.

In [15]:

```
C1 = Circle(P1, 3)
C2 = Circle(P2, P3, P4)
```

In [16]:

```
C2.area
```

Out[16]:

In [17]:

```
C2.radius
```

Out[17]:

In [18]:

```
C2.equation()
```

Out[18]:

In [19]:

```
C2.center
```

Out[19]:

In [20]:

```
C2.circumference
```

Out[20]:

Computing intersections with other objects, checking whether a line is tangent to a circle, or finding the tangent lines through an non-interior point, are simple tasks too:

In [21]:

```
C2.intersection(C1)
```

Out[21]:

In [22]:

```
C2.intersection(S1)
```

Out[22]:

In [23]:

```
C2.is_tangent(L2)
```

Out[23]:

In [24]:

```
C1.tangent_lines(P4)
```

Out[24]:

One of the most useful basic geometric concept is the **triangle**. We need robust and fast algorithms to manipulate and extract information from them. Let us show first the definition of one, together with a series of queries to describe its properties:

In [25]:

```
T = Triangle(P1, P2, P3)
```

In [26]:

```
T.area # Note it gives a signed area
```

Out[26]:

In [27]:

```
T.angles
```

Out[27]:

In [28]:

```
T.sides
```

Out[28]:

In [29]:

```
T.perimeter
```

Out[29]:

In [30]:

```
T.is_right() # Is T a right triangle?
```

Out[30]:

In [31]:

```
T.is_equilateral() # Is T equilateral?
```

Out[31]:

In [32]:

```
T.is_scalene() # Is T scalene?
```

Out[32]:

In [33]:

```
T.is_isosceles() # Is T isosceles?
```

Out[33]:

Next, note how easily we can obtain representation of the different segments, centers, and circles associated with triangles, as well as the *medial triangle* (the triangle with vertices at the midpoints of the segments).

In [34]:

```
T.altitudes
```

Out[34]:

In [35]:

```
T.orthocenter # Intersection of the altitudes
```

Out[35]:

In [36]:

```
T.bisectors() # Angle bisectors
```

Out[36]:

In [37]:

```
T.incenter # Intersection of angle bisectors
```

Out[37]:

In [38]:

```
T.incircle
```

Out[38]:

In [39]:

```
T.inradius
```

Out[39]:

In [40]:

```
T.medians
```

Out[40]:

In [41]:

```
T.centroid # Intersection of the medians
```

Out[41]:

In [42]:

```
T.circumcenter # Intersection of perpendicular bisectors
```

Out[42]:

In [43]:

```
T.circumcircle
```

Out[43]:

In [44]:

```
T.circumradius
```

Out[44]:

In [45]:

```
T.medial
```

Out[45]:

Some other interesting operations with triangles:

- Intersection with other objects
- Computation of the minimum distance from a point to each of the segments.
- Checking whether two triangles are similar.

In [46]:

```
T.intersection(C1)
```

Out[46]:

In [47]:

```
T.distance(T.circumcenter)
```

Out[47]:

In [48]:

```
T.is_similar(Triangle(P1, P2, P4))
```

Out[48]:

The other basic geometrical objects currently coded in the Geometry module are

`LinearEntity`

. This is a superclass (which we never use directly), with three subclasses`Segment`

,`Line`

and`Ray`

.`LinearEntity`

enjoys the following basic methods:`are_concurrent(o1, o2, ..., on)`

`are_parallel(o1, o2)`

`are_perpendicular(o1, o2)`

`parallel_line(self, Point)`

`perpendicular_line(self, Point)`

`perpendicular_segment(self, Point)`

`Ellipse`

. An object with a center, together with horizontal and vertical radii.`Circle`

is, as a matter of fact, a subclass of`Ellipse`

with both radii equal.`Polygon`

. A superclass we can instantiate by listing a set of vertices. Triangles are a subclass of`Polygon`

, for example. The basic methods of`Polygon`

are`area`

`perimeter`

`centroid`

`sides`

`vertices`

`RegularPolygon`

. This is a subclass of`Polygon`

, with extra attributes:`apothem`

`center`

`circumcircle`

`exterior_angle`

`incircle`

`interior_angle`

`radius`

For more information about this module, refer to the official `sympy`

documentation at docs.sympy.org/dev/modules/geometry.html

There is also a *non-basic* geometric object: A `Curve`

, which we define by providing parametric equations, together with the interval of definition of the parameter. It currently does not have many useful methods, other than those describing its constructors. Let us illustrate how to deal with these objects. For example, a three-quarters arc of an ellipse could be coded as follows:

In [49]:

```
from sympy import var, pi, sin, cos
var('t', real=True)
Arc = Curve((3*cos(t), 4*sin(t)), (t, 0, 3*pi/4))
```

To end the exposition on basic objects from the geometry module in the `sympy`

library, we must mention that we may apply any basic affine transformations to any of the previous objects. This is done by combination of the methods `reflect`

, `rotate`

, `translate`

and `scale`

.

In [50]:

```
T.reflect(L1)
```

Out[50]:

In [51]:

```
T.rotate(pi/2, P2)
```

Out[51]:

In [52]:

```
T.translate(5,4)
```

Out[52]:

In [53]:

```
T.scale(9)
```

Out[53]:

In [54]:

```
Arc.rotate(pi/2, P3).translate(pi,pi).scale(0.5)
```

Out[54]:

With these basic definitions and operations, we are ready to address more complex situations. Let us explore these new challenges next.

Also called **algorithmic geometry**, the applications of this field are plenty: in robotics, these are used to solve visibility problems, and motion planning, for instance. Similar applications are employed to design route planning or search algorithms in **geographic information systems** (**GIS**).

Let us describe the different categories of problems, making emphasis on the tools for solving them that are available in the `scipy`

stack.

The fundamental problems in this category are the following:

**Convex hulls**: Given a set of points in space, find the smallest convex polyhedron containing them.**Voronoi diagrams**: Given a set of points in space (the*seeds*), compute a partition in regions consisting of all points closer to each seed.**Triangulations**: Partition the plane with triangles, in a way that two triangles are either disjoint, or otherwise they share an edge or a vertex. There are different triangulations depending on the input objects, or constraints on the properties of the triangles.**Shortest paths**: Given a set of obstacles in a space, and two points, find the shortest path between the points that does not intersect any of the obstacles.

The problems of computation of convex hulls, basic triangulations, and Voronoi diagrams are intimately linked. The theory that explains this beautiful topic is explained in detail in a monograph in Computer Science titled

Computational Geometry, written by Franco Preparata and Michael Shamos. It was published by Springer-Verlag in 1985.

While it is possible to compute the convex hull of a reasonably large set of points in the plane through the geometry module of the library `sympy`

, this is not recommendable. A much faster and reliable code is available in the module `scipy.spatial`

through the routine `ConvexHull`

, which is in turn a wrapper to `qconvex`

, from the `Qhull`

libraries www.qhull.org. This routine also allows the computation of convex hulls in higher dimensions. Let us compare both methods, with the famous *Lake Superior* polygon, `superior.poly`

.

`poly`

files representplanar straight line graphs---a simple list of vertices and edges, together with information about holes and concavities, in some cases. The running example can be downloaded from www.math.sc.edu/~blanco/superior.poly. It contains a polygonal description of the coastline of Lake Superior, with 7 holes (for the islands), 518 vertices, and 518 edges.For a complete description of the

`poly`

format, refer to www.cs.cmu.edu/~quake/triangle.poly.html. With that information, we can write a simple reader without much effort. This is an example:

In [104]:

```
from numpy import array
def read_poly(file_name):
"""
Simple poly-file reader, that creates a python dictionary
with information about vertices, edges and holes.
It assumes that vertices have no attributes or boundary markers.
It assumes that edges have no boundary markers.
No regional attributes or area constraints are parsed.
"""
output = {'vertices': None, 'holes': None, 'segments': None}
# open file and store lines in a list
file = open(file_name, 'r')
lines = file.readlines()
file.close()
lines = [x.strip('\n').split() for x in lines]
# Store vertices
vertices= []
N_vertices, dimension, attr, bdry_markers = [int(x) for x in lines[0]]
# We assume attr = bdrt_markers = 0
for k in range(N_vertices):
label, x, y = [items for items in lines[k+1]]
vertices.append([float(x), float(y)])
output['vertices']=array(vertices)
# Store segments
segments = []
N_segments, bdry_markers = [int(x) for x in lines[N_vertices+1]]
for k in range(N_segments):
label, pointer_1, pointer_2 = [items for items in lines[N_vertices+k+2]]
segments.append([int(pointer_1)-1, int(pointer_2)-1])
output['segments'] = array(segments)
# Store holes
N_holes = int(lines[N_segments+N_vertices+2][0])
holes = []
for k in range(N_holes):
label, x, y = [items for items in lines[N_segments + N_vertices + 3 + k]]
holes.append([float(x), float(y)])
output['holes'] = array(holes)
return output
```

In [56]:

```
import numpy as np
from scipy.spatial import ConvexHull
import matplotlib.pyplot as plt
%matplotlib inline
```

In [57]:

```
lake_superior = read_poly("/Users/francisco.blanco.silva/Dropbox/Documents/Books/Mastering/chapter6/superior.poly")
vertices_ls = lake_superior['vertices']
```

In [58]:

```
%time hull = ConvexHull(vertices_ls)
```

In [80]:

```
plt.figure(figsize=(18, 18))
plt.xlim(vertices_ls[:,0].min()-0.01, vertices_ls[:,0].max()+0.01)
plt.ylim(vertices_ls[:,1].min()-0.01, vertices_ls[:,1].max()+0.01)
plt.axis('off')
plt.axes().set_aspect('equal')
plt.plot(vertices_ls[:,0], vertices_ls[:,1], 'b.')
for simplex in hull.simplices:
plt.plot(vertices_ls[simplex, 0], vertices_ls[simplex, 1], 'r-')
plt.show()
```

Let us now illustrate a few advanced uses of `ConvexHull`

. First, the computation of the convex hull of a random set of points in the 3D space. For visualization, we will use the `mayavi`

libraries.

In [ ]:

```
from mayavi import mlab
```

In [ ]:

```
points = np.random.rand(320, 3)
hull = ConvexHull(points)
X = hull.points[:, 0]
Y = hull.points[:, 1]
Z = hull.points[:, 2]
mlab.triangular_mesh(X, Y, X, hull.simplices, colormap='gray',
opacity=0.5, representation='wireframe')
mlab.show()
```

Computing the Voronoi diagram of a set of vertices (our *seeds*) can be done with the routine `Voronoi`

(and its companion `voronoi_plot_2d`

for visualization) from the module `scipy.spatial`

. The routine `Voronoi`

is in turn a wrapper to the function `qvoronoi`

from the `Qhull`

libraries, with the following default `qvoronoi`

controls: `qhull_option='Qbb Qc Qz Qx'`

if the dimension of the points is greater than 4, and `qhull_options='Qbb Qc Qz'`

otherwise. For the computation of the **furthest-site Voronoi diagram**, instead of the nearest-site, we would use the extra control `Qu`

.

In [60]:

```
from scipy.spatial import Voronoi, voronoi_plot_2d
```

In [78]:

```
vor = Voronoi(vertices_ls)
plt.figure(figsize=(10, 10))
ax = plt.subplot(111, aspect='equal')
voronoi_plot_2d(vor, ax=ax)
plt.xlim( 0.45, 0.50)
plt.ylim(-0.40, -0.35)
plt.show()
```

- The small dots are the original seeds with
`x`

-coordinates between`0.45`

and`0.50`

, and`y`

-coordinates between`-0.40`

and`-0.35`

. We access those values either from the original list`vertices_ls`

, or from`vor.points`

. - The plane gets partitioned into different regions (the Voronoi cells), one for each seed. These regions contain all points in the plane which are closest to its seed. Each region receives an index, which is not necessarily the same index as the index of its seed in the
`vor.points`

list. To access the corresponding region to a given seed, we use`vor.point_region`

.

In [62]:

```
vor.point_region
```

Out[62]:

- Each Voronoi cell is defined by its delimiting vertices and edges (also known as
*ridges*in Voronoi jargon). The list with the coordinates of the computed vertices of the Voronoi diagram can be obtained with`vor.vertices`

. These vertices were represented as bigger dots in the previous image, and are easily identifiable because they are always at the intersection of at least two edges --- while the seeds have no incoming edges.

In [63]:

```
vor.vertices
```

Out[63]:

- For each of the regions, we can access the set of delimiting vertices with
`vor.regions`

. For instance, to obtain the coordinates of the vertices that delimit the region around the 4th seed, we could issue

In [64]:

```
[vor.vertices[x] for x in vor.regions[vor.point_region[4]]]
```

Out[64]:

Care must be taken with the previous step: Some of the vertices of the Voronoi cells are not actual vertices, but lie *at infinity*. When this is the case, they are identified with the index `-1`

. In this situation, to provide an accurate representation of a ridge of these characteristics we must use the knowledge of the two seeds whose contiguous Voronoi cells intersect on said ridge---since the ridge is perpendicular to the segment defined by those two seeds. We obtain the information about those seeds with `vor.ridge_points`

In [65]:

```
vor.ridge_points
```

Out[65]:

The first entry of `vor.ridge_points`

can be read as follows: There is a ridge perpendicular to both the first and second seeds.

There are other attributes of the object `vor`

that we may use to inquire properties of the Voronoi diagram, but the ones we have described should be enough to replicate the previous image. We leave this as a nice exercise:

- Gather the indices of the seeds from
`vor.points`

that have their`x`

- and`y`

-coordinates in the required window. Plot them. - For each of those seeds, gather information about vertices of their corresponding Voronoi cells. Plot those vertices not at the infity with a different style as the seeds.
- Gather information about the ridges of each relevant region, and plot them as simple thin segments. Some of the ridges cannot be represented by their two vertices. In that case, we use the information about the seeds that determine them.

A **triangulation** of a set of vertices in the plane is a division of the convex hull of the vertices into triangles, satisfying one important condition. Any two given triangles:

- must be disjoint, or
- must intersect only at one common vertex, or
- must share one common edge.

These plain triangulations have not much computational value, since some of its triangles might be too skinny --- this leads to uncomfortable rounding errors, or computation or erroneous areas, centers, etc. Among all possible triangulations, we always seek one where the properties of the triangles are somehow balanced.

With this purpose in mind, we have the **Delaunay triangulation** of a set of vertices. This triangulation satisfies an extra condition: None of the vertices lies in the interior of the circumcircle of any triangle. We refer to triangles with this property as **Delaunay triangles**.

For this simpler setting, in the module `scipy.spatial`

, we have the routine `Delaunay`

, which is in turn a wrapper to the function `qdelaunay`

from the `Qhull`

libraries, with the `qdelaunay`

controls set exactly as for the Voronoi diagram computations.

In [135]:

```
from scipy.spatial import Delaunay, delaunay_plot_2d
```

In [138]:

```
tri = Delaunay(vertices_ls)
plt.figure(figsize=(18, 18))
ax = plt.subplot(111, aspect='equal')
delaunay_plot_2d(tri, ax=ax)
plt.show()
```

It is possible to generate triangulations with imposed edges too. Given a collection of vertices and edges, a **constrained Delaunay triangulation** is a division of the space into triangles with those prescribed features. The triangles in this triangulation are not necessarily Delaunay.

We can accomplish this extra condition sometimes by subdivision of each of the imposed edges. We call this triangulation **conforming Delaunay**, and the new (artificial) vertices needed to subdivide the edges are called **Steiner points**.

A **constrained conforming Delaunay triangulation** of an imposed set of vertices and edges satisfies a few more conditions, usually setting thresholds on the values of angles or areas of the triangles. This is achieved by introducing a new set of **Steiner points**, which are allowed anywhere, not only on edges.

To achieve these high-level triangulations, we need to step outside of the `scipy`

stack. We have a `python`

wrapper to the amazing implementation of mesh generators www.cs.cmu.edu/~quake/triangle.html by Richard Shewchuck . This wrapper, together with examples and other related functions, can be installed by issuing

In [ ]:

```
!pip install triangle
```

For more information on this module, refer to the documentation online from its author, Dzhelil Rufat, at dzhelil.info/triangle/index.html

Let us compute those different triangulations for our running example. We use once again the `poly`

file with the features of Lake Superior, which we read into a dictionary with all the information about vertices, segments and holes. The first example is that of the constrained Delaunay triangulation (`cndt`

). We accomplish this task with the flag `p`

(indicating that the source is a *planar straight line graph*, rather than a set of vertices).

In [69]:

```
from triangle import triangulate, plot as tplot
```

In [76]:

```
cndt = triangulate(lake_superior, 'p')
plt.figure(figsize=(18, 18))
ax = plt.subplot(111, aspect='equal')
tplot.plot(ax, **cndt)
plt.show()
```

The next step is the computation of a conforming Delaunay triangulation (`cfdt`

). We enforce Steiner points on some segments to ensure as many Delaunay triangles as possible. We achieve this with extra flag `D`

.

In [71]:

```
cfdt = triangulate(lake_superior, 'pD')
```

But slight or no improvements with respect to the previous diagram can be observed in this case. The real improvement arises when we further impose constraints in the values of minimum angles on triangles (with the flag `q`

), or in the maximum values of the areas of triangles (with the flag 'a'). For instance, if we require a constrained conforming Delaunay triangulation (`cncfdt`

) in which all triangles have a minimum angle of at least 20 degrees, we issue the following command

In [75]:

```
cncfq20dt = triangulate(lake_superior, 'pq20D')
plt.figure(figsize=(18,18))
ax = plt.subplot(111, aspect='equal')
tplot.plot(ax, **cncfq20dt)
plt.show()
```