In [15]:
from ICGChapter00 import *

In [16]:
# Comment out the line below if the graphics do not display
%config InlineBackend.figure_format = 'retina'

In [17]:
%%javascript
IPython.OutputArea.prototype._should_scroll = function(lines) {return false;}


# Chapter 1: Points and lines¶

In this chapter we are going to consider points, and the different types of line. Let's start with some definitions:

Euclidean plane: a flat two dimensional surface.

Point: a location.

General position for points: a set of points is said to be in general position if no three points are collinear.

The Euclidean plane is our work surface on which we will construct our geometry. This is 2 dimensional, so henceforth, when we refer to "point" we will mean "2D point". We will be mostly working in two dimensions.

General position for points is quite an important concept, because, as we will see shortly, computational geometry algorithms that involve sorting points often fail if the points are not in general position. We will indicate these as we encounter them. It is worth noting that there are other, more constrained, definitions for general position than the one given above. For example, some algorithms only work if no two points have the same x or y. We will indicate any such additional constraints as they arise.

In [18]:
displayTOC("ICGChapter01.ipynb")


# Points¶

Contents

We will represent a point as a Python class with two attributes, x and y, that will usually be integers. Why integers? Not for speed reasons (we are only dealing with very small examples in this text) but because integers allow exact placement of points in our interactive demonstrations. This makes it easy to explore the characteristics of the computational geometry algorithms including any special cases. In fact the difference in speed between processing integers and floating point numbers is very small for modern computers, and in some cases, floating point operations may actually be faster.

Here is our Point class:

In [28]:
# An immutable point
class Point:
"""A simple Point class. It comprises two numbers, x and y."""
def __init__(self, x=0, y=0):
self.__x = x
self.__y = y

@property
def x(self):
"""Return the x value of the point."""
return self.__x

@property
def y(self):
"""Return the y value of the point."""
return self.__y

def __eq__(self, p):
return (self.x == p.x) and (self.y == p.y)

def __ne__(self, p):
return (self.x != p.x) or (self.y != p.y)

return Point(self.x + p.x, self.y + p.y)

def __sub__(self, p):
return Point(self.x - p.x, self.y - p.y)

def __mul__(self, p):
return Point(self.x * p.x, self.y * p.y)

def __truediv__(self, p):
return Point(self.x/p.x, self.y/p.y)

def __hash__(self):
return hash((self.x, self.y))

def __iter__(self):
return iter((self.x, self.y))

def __str__(self):
return "({},{})".format(self.x, self.y)

def __repr__(self):
return "Point({!r},{!r})".format(self.x, self.y)

def roundXY(self):
"""Round the x an y values to make them integers."""
return Point(round(self.x), round(self.y))


A key aim of our Point class is to make Points as "Pythonic" as possible. What does this mean? It means to ensure that Point objects integrate as seamlessly as possible into the rest of the Python language. Python has a lot of built-in facilities defined in its data model (the Python term for its object model), and we want the Point class to implement as many of these as are appropriate. The methods whose names are surrounded by double underscores are special methods that allow Point objects to integrate with Python. They have nothing to do with computational geometry per se, but, as we will see, they make the implementation of computational geometry algorithms very simple indeed.

The best possible reference for Python special functions is, "Fluent Python", by Luciano Ramalho:

http://shop.oreilly.com/product/0636920032519.do

Every serious Python programmer needs to read this book!

The facilities the Point class offers are:

• Holding the x and y values of the Point: variables x and y.
• Test for equality between two Points: (__eq__ and __ne__)
• Basic arithmetic between two Points:
• Subtraction: __sub__
• Multiplication: __mul__
• Division: __truediv__
• Identity based on the x and y attribute pair: __hash__ returns a unique identifier for a Point object that is based on the values of its x and y attributes. This method allows us to store Points in Python dictionaries and sets. We assume that two Points with the same x and y are the same point.
• Destructuring: __iter__ allows us to access the x and y attributes of a Point object e.g. a, b = Point(x,y).
• Printable representation: __str__
• Code representation: __repr__
• Conversion to integer: roundXY(...) returns a new Point object with integer x and y values by rounding the current Point object x and y values.

Points are immutable - once the $x$ and $y$ values of the Point have been set when the Point is created, they can't be changed. In that respect they behave much like the integers. For example, the value of the integer 5 is always 5, similarly the value of the Point(1,1) is always (1,1).

Here are some simple examples of what you can do with Points:

In [29]:
#Demo
def _demoPoints():
p1 = Point(1,1)
p2 = Point(2,2)
p3 = Point(3,4)
print("p3 == p1 is {}".format(p3 == p1))
print("p3 != p1 is {}".format(p3 != p1))
print("p1 + p2 = {}".format(p1 + p2))
print("p1 - p2 = {}".format(p1 - p2))
print("p2 * p3 = {}".format(p2 * p3))
print("p2 / p3 = {}".format(p2 / p3))
print("(p2 / p3).roundXY() = {}".format((p2 / p3).roundXY()))
print("p1 as string = {!s}".format(p1))
print("p3 as code = {!r}".format(p3))
print()
x, y = Point(5,7)
print("Destructuring: the expression x, y = Point(5,7) gives x = {}, y = {}".format(x,y))

_demoPoints()
#Demo

p3 == p1 is False
p3 != p1 is True
p1 + p2 = (3,3)
p1 - p2 = (-1,-1)
p2 * p3 = (6,8)
p2 / p3 = (0.6666666666666666,0.5)
(p2 / p3).roundXY() = (1,0)
p1 as string = (1,1)
p3 as code = Point(3,4)

Destructuring: the expression x, y = Point(5,7) gives x = 5, y = 7


As you can see, Point behaves very much like a built-in Python class.

## The distance between two points¶

Contents

The distance between two Points $(x_1,y_1)$ and $(x_2, y_2)$ is given from Pythagoras's Theorem by:

$$d = \sqrt{(x_2 - x_1)^2 + (y_2 - y_1)^2}$$

We will create a function to calculate this distance as follows:

In [30]:
def distance(p1, p2):
"""Calculate the distance between two Point objects."""
return math.hypot(p2.x - p1.x, p2.y - p1.y)


The function math.hypot(...) takes two arguments, x and y, and returns their Euclidean norm, sqrt(xx + yy). This is the length of the vector from the origin to point (x, y).

Here is a simple demo program.

In [31]:
#Demo
def _demoDistance():
p1 = Point(1,1)
p2 = Point(2,2)
print("Distance between p1 and p2 = {}".format(distance(p1,p2)))

_demoDistance()
#Demo

Distance between p1 and p2 = 1.4142135623730951


# Lists of Points - the PointList class¶

Contents

We will define a PointList class that holds a list of Points as follows:

In [32]:
class PointList:
def __init__(self, *args):
"""A simple class to hold a list of Point objects. Can accept points in several formats."""
# args is a list of Points
if type(args[0]) is Point:
self.points = list(args)
# args is a list of tuples, (0,1) etc.
elif type(args[0]) is tuple:
self.points = [Point(t[0], t[1]) for t in args]
# args is a list containing a list of Points [Point(...), Point(...)...]
elif type(args[0]) is list:
self.points = list(args[0])
else:
raise ValueError("PointList constructor - args contains wrong type.")

def copy(self):
"""Return a deep copy of the PointList object."""
return copy.deepcopy(self)

def __len__(self):
return len(self.points)

def __eq__(self, pl):
if len(self.points) != len(pl.points):
return False
for i,p in enumerate(self.points):
if p != pl.points[i]:
return False
return True

def __iter__(self):
return iter(self.points)

def __getitem__(self, i):
if isinstance(i, slice):
return type(self)(self.points[i])
elif isinstance(i, int):
return self.points[i]

def __contains__(self, p):
for pt in self.points:
if pt == p: return True
return False

def __str__(self):
return "(" + ",".join(["{!s}".format(p) for p in self.points]) + ")"

def __repr__(self):
return type(self).__name__ + '(' + ",".join(["{!r}".format(p) for p in self.points]) + ")"

# Returns a new PointList object
def roundPoints(self):
"""Return a copy of the PointList in which all Point objects have x and y values rounded to integers."""
return PointList([p.roundXY() for p in self.points])

# These methods modify the PointList object
def append(self, p):
"""Modify the PointList object by appending a Point object."""
self.points.append(p)

def appendAll(self, plist):
"""Modify the PointList object by appending a list of Point objects."""
for p in plist:
self.append(p)

def remove(self, p):
"""Modify the PointList object by removing a Point object."""
self.points.remove(p)

def removeAll(self, plist):
"""Modify the PointList object by removing a list of Point objects."""
for p in plist:
self.remove(p)


This is a very Pythonic class, and you can generally use it interchangeably with list. In fact, we usually prefer native lists, and we only use PointList when we need to take advantage of its special features such as its flexible __init__ method.

The __init__ method is the constructor for the class, and this can take several different types of arguments as follows:

In [33]:
#Demo
def _demoPointListInit():
pl1 = PointList((1,1), (2,3), (4,6), (7,9))            # From tuples
pl2 = PointList(Point(4,5), Point(10,8), Point(1,7))   # From Points
pl3 = PointList([Point(4,5), Point(10,8), Point(1,7)]) # From a list of Points
print("pl1 = {}".format(pl1))
print("pl2 = {}".format(pl2))
print("pl3 = {}".format(pl3))

_demoPointListInit()
#Demo

pl1 = ((1,1),(2,3),(4,6),(7,9))
pl2 = ((4,5),(10,8),(1,7))
pl3 = ((4,5),(10,8),(1,7))


Typically, you will hardly ever use the "from list of Points" variant, but it is useful for implementation reasons.

The __getitem__ method needs some consideration. This method allows you to use Python indexing and slicing on objects of the PointList class. If the method is passed an integer, then it returns the Point at that index. If, however, it is passed a slice, then it returns a new sub list of the PointList containing the appropriate elements. For example:

In [34]:
#Demo
def _demoPointListGetitem():
pl = PointList((4,4), (5,2), (6,3), (7,4), (8,5))
print(pl[1:3])   # A simple index
print(pl[0:6:2]) # A slice from Point 0 to 6 step 2

_demoPointListGetitem()
#Demo

((5,2),(6,3))
((4,4),(6,3),(8,5))


In some of the more advanced computational geometry algorithms, it is useful to be able to add and remove Points from the PointList as follows:

In [35]:
#Demo
def _demoAppendAndRemove():
pl = PointList((4,4), (5,2), (6,3), (7,4), (8,5))
print(pl)

pl.append(Point(0,0))
print(pl)
pl.appendAll((Point(1,1), Point(2,2)))
print(pl)

pl.remove(Point(0,0))
print(pl)
pl.removeAll((Point(1,1), Point(2,2)))
print(pl)

_demoAppendAndRemove()
#Demo

((4,4),(5,2),(6,3),(7,4),(8,5))
((4,4),(5,2),(6,3),(7,4),(8,5),(0,0))
((4,4),(5,2),(6,3),(7,4),(8,5),(0,0),(1,1),(2,2))
((4,4),(5,2),(6,3),(7,4),(8,5),(1,1),(2,2))
((4,4),(5,2),(6,3),(7,4),(8,5))


Ideally, we would have liked our PointList class to be immutable like our Point class. That would mean that whenever you added or removed a Point, you got back a new PointList rather than just updating the existing PointList in place. However, this proved to be more trouble than it was worth for the following reasons:

1. The Python list is mutable, so in order for our PointList to be as compatible with the Python list as possible it also has to be mutable. Otherwise append(...) and remove(...) don't work as expected.

2. Python doesn't have immutability baked-in like languages such as Clojure, so there can be significant performance issues resulting from peforming deep copies of the PointList on each append or remove.

PointList objects are mutable, so when you append or remove an item from a the object, it object changes! If you don't want to change the original, then use the copy(...) method to create a deep copy of it. You will see that several algorithms we present later are destructive, so they always make a copy of what is passed in to them first.

## Integer point lists¶

Contents

Sometimes we will get list of points that have floating point values (for example, every regular polygon apart from the square must have one or more non-integer vertexes), but in this book we generally want to work on an integer grid. The PointList __int__(...) method allows us to round the values to make them integers. Of course, this changes the positions some or all of the Points in the PointList by $\pm 0.5$.

In [36]:
#Demo
def _demoRoundPoints():
pl = PointList((2.0, 3.141), (5.23, 4.29), (6.7, 8.9))
pli = pl.roundPoints()
print(pl)
print(pli)

_demoRoundPoints()
#Demo

((2.0,3.141),(5.23,4.29),(6.7,8.9))
((2,3),(5,4),(7,9))


## Ordered Point sets¶

Contents

We sometimes need an ordered set of Points. This is a list of Points in which there are no duplicates. The recommended way to do this in Python 3 is to use the facilities of the OrderedDict class. We create an OrderedDict using the Points as keys, and this automatically removes duplicates. We then convert this back to a list. Here is a simple OrderedPointSet class that uses this technique:

In [37]:
class OrderedPointSet(PointList):
def __init__(self, *args):
"""A mathematical set of Point objects with no duplicates but with order preserved."""
PointList.__init__(self, *args)
# Turn the points attribute into a set
self.points = list(OrderedDict.fromkeys(self.points))


The OrderedPointSet removes duplicates automatically:

In [38]:
#Demo
print(OrderedPointSet((1,2), (3,4), (5,6), (1,2)))
#Demo

((1,2),(3,4),(5,6))


Python has a built-in set class that is unordered, but which supports the full range of set operations: union, intersection etc. If we are not concerned with preserving order, then will will use a native set of Points rather than OrderedPointSet.

# Extreme Points of a PointList¶

Contents

The first true computational geometry algorithm we will look at is finding the extreme Points of a PointList. Here is a definition of extreme point:

Extreme point: an extreme point of a point list is a point that is the furthest along a given direction.

A very useful function for working with point lists is one that finds the extreme points in the X and Y directions. On the X axis, these are the rightmost and leftmost points, and on the Y axis, the topmost and bottommost points. The extreme points are often unique, but not always. Any one of them can overlap with another one e.g. a point which is rightmost AND bottommost. There must be at least two unique extreme points (in which case the other points are positioned on a line between them) and, obviously, there can be no more than four. If there are $n_e$ unique extreme points then:

$$2 <= n_e <= 4$$

When there are two (or more) candidate points for a given extreme (for example two or more points with the same extreme Y coordinate), the algorithm we use here gives you the last one it finds:

In [39]:
RIGHTMOST = 0
TOPMOST = 1
LEFTMOST = 2
BOTTOMMOST = 3

def extremePoints(plist):
"""Return a list of the rightmost, topmost, leftmost and bottommost Points of a Point list."""
rightmst = leftmst = topmst = bottommst = plist[0];
for p in plist:
if (p.x > rightmst.x): rightmst = p
if (p.y > topmst.y): topmst = p
if (p.x < leftmst.x): leftmst = p
if (p.y < bottommst.y): bottommst = p
return PointList([rightmst, topmst, leftmst, bottommst])


For example:

In [40]:
#Demo
extremePoints([Point(-1,-5), Point(2,-2), Point(-4,-4), Point(2,1), Point(2,2)])
#Demo

Out[40]:
PointList(Point(2,-2),Point(2,2),Point(-4,-4),Point(-1,-5))

# Random sets of Points not in general position¶

Contents

We often need a set of random points as input to computational geometry algorithms. For convenience, we will generate points on an $n$ x $n$ grid, so the maximum number of points we can have is obviously $n^2$, when there is a point on each intersection of the grid.

We will generally be working in a 11 x 11 grid that goes from -5 to +5 in both X and Y, so in this case, the maximum number of points is 121.

Here is a very simple constructive algorithm for generating a set, $P$, of $m$ random points on an $n$ x $n$ grid that ranges from $min$ to $max$ in both $X$ and $Y$. Essentially, we just keep adding random points to a set until we get the number we want:

1. Let $P = \emptyset$ (the empty set)

2. Repeat while $|P| < m$:

2.1 Generate a random point, $p = (x, y)$

2.2 $P = P \cup \{p\}$

Where $x \in \mathbb{Z}$, $y \in \mathbb{Z}$, $min \le x \le max$, $min \le y \le max$ and $m \le n^2$.

This algorithm has the huge advantage that it is very, very fast so we can easily generate sets of thousands of points. However, the points are not guaranteed to be in general position. We will find out how to generate random point sets in general position later.

In [41]:
def randomPointSetNonGP(numberOfPoints, start=-5, end=5):
"""Return a random set of Points with no duplicates, but not guaranteed to be in general position."""
maxNumberOfPoints = (end - start + 1)**2
numberOfPoints = min(maxNumberOfPoints, numberOfPoints)
randomPts = set()
while len(randomPts) < numberOfPoints:
p = Point(random.randint(start, end), random.randint(start, end))
return list(randomPts)


Here is a demo program:

In [42]:
#Demo
randomPointSetNonGP(5)
#Demo

Out[42]:
[Point(1,2), Point(0,-5), Point(3,3), Point(-3,-3), Point(-1,-2)]

# Displaying Points and other geometrical objects¶

Contents

PyPlot is the standard way to display plots and other graphics in Jupyter notebooks (such as this one), but it is complex. Our goal in this section is to create a canvas on which we can graphically demonstrate computational geometry using very short and simple programs. To do this, we will create a class called CGCanvas that provides a blank canvas on which to draw. We will then define a set of drawing functions on an as-needed basis to draw our computational geometry objects on the canvas. These functions will provide an easy to use interface to the relevant parts of PyPlot. The combination of CGCanvas and associated drawing functions provides a mini graphics language tailored to our needs. This is known as a domain specific language or DSL.

PyPlot is very flexible and each drawing function has many options expressed as Python keyword arguments. Our graphics language will provide lightweight wrapper that passes through keyword arguments (**kwargs) to the underlying PyPlot functions. You can find out what options are available here:

https://matplotlib.org/users/beginner.html

Our wrapper always provides you with sensible defaults, so that you typically don't have to worry about the various options if you don't want to. However, you can still access many of them using the keyword argument pass-though mechanism.

Here is the code for the blank canvas on which we will draw:

In [43]:
class CGCanvas:
"""A canvas class on which we can draaw graphics."""
def __init__(self, title, xlabel='X', ylabel='Y', p1=Point(-7, -7), p2=Point(7, 7)):
self.fig = plt.figure()
self.fig.set_size_inches(CANVAS_WIDTH, CANVAS_HEIGHT)
self.ax = self.fig.add_subplot(111, aspect='equal')
plt.title(title)
plt.xlabel(xlabel)
plt.ylabel(ylabel)
plt.xticks(range(p1.x, p2.x))
plt.yticks(range(p1.y, p2.y))
self.ax.grid(True)
self.ax.set_xlim([p1.x, p2.x])
self.ax.set_ylim([p1.y, p2.y])

def show(self):
"""Show the canvas, displaying any graphics drawn on it."""
plt.show()
plt.clf()
plt.cla()
plt.close()


This is what it looks like:

In [44]:
#Demo
def _demoCGCanvas():
CGCanvas(title="This is our standard blank canvas").show()

_demoCGCanvas()
#Demo


We can now define some functions to draw on this:

• circle(...) - draws a filled circle.
• circles(...) - draws a sequence of filled circles, one for each point in the point list parameter.
• interactor(...) - draws an empty circle used to mark an interactive (movable) point.
• text(...) - text, including mathematical notation in LaTeX format.
• annotation(...) - an arrow.
• lineFromTo(...) - a line between two points.
• path(...) - a closed polygonal shape that is usually filled.

Here are the implementations:

In [45]:
def circle(cgc, pt, radius=0.15, **kwargs):
"""Draw a circle."""

def circles(cgc, plist, radius=0.15, **kwargs):
"""Draw circles at each point in plist."""
for p in plist:
circle(cgc, p, radius, **kwargs)

def interactor(cgc, pt, radius=0.5, dot=True, fill=None, **kwargs):
"""Draw an interactor, which is a holow circle around an optional disk."""
if dot: circle(cgc, pt, radius = 0.15, fill=True, **kwargs)

def text(cgc, pt, text, **kwargs):
"""Draw text."""
cgc.ax.text(pt.x, pt.y, text, **kwargs)

def annotation(cgc, pt1, pt2, arrowstyle="->", text="", **kwargs):
"""Draw an annotation comprising a line with an arrow head."""
cgc.ax.annotate(text, xy=(pt1.x, pt1.y), xycoords='data',
xytext=(pt2.x, pt2.y), textcoords='data',
arrowprops=dict(arrowstyle=arrowstyle, connectionstyle="arc3"), **kwargs)

def lineFromTo(cgc, p1, p2, **kwargs):
"""Draw a line from Point object p1 to Point object p2."""
cgc.ax.add_patch(patches.FancyArrow(p1.x, p1.y, p2.x - p1.x, p2.y - p1.y, **kwargs))

def path(cgc, poly, **kwargs):
"""Draw a polyline between the Point objects in poly."""
pathVertexes = [(p.x, p.y) for p in poly] + [(poly[0].x, poly[0].y)]
pathCodes = []
for i in range(2, len(pathVertexes)):
pathCodes.append(Path.LINETO)
pathCodes = [Path.MOVETO] + pathCodes + [Path.CLOSEPOLY]


We often need to annotate a point with text, and for clarity, the text needs to be offset from the point by some amount. Typically, the offsets for the x and y coordinates of a point should have the same sign as x and y: so positive x will move more to the right, negative x more to the left, positive y more above, and negative y more below the original point. This will usually put the offset point in a reasonable position relative to the original point.

In [46]:
def offsetPoint(pt, amount=Point(0.5, 0.5)):
"""Offset a point left, right, up or down, according to the quadrant it is in."""
return pt + amount*Point(np.sign(pt.x),np.sign(pt.y))


Lists of points are usually labeled by their index in the list as $(0, 1, 2, ... n)$. The following function has a default template that labels points 0, 1, 2, 3, ... etc. If you give it a template such as "p{}", this will label the points p1, p2, p3 ... and so on. You can also pass in a special offset function for the labels. The default is offsetPoint(...) described above.

In [47]:
def labelPointList(cgc, plist, template = "{}", offset=offsetPoint, **kwargs):
"""Label a list of Points. The template parameter can be used to create complex labels."""
for i, p in enumerate(plist):
text(cgc, offset(p), text=template.format(i), **kwargs)


We often need to use an indexed set of colors. The function col(...) gives access to the matplotlib, "cyclic", colors, C1 to C9 via an integer index that wraps around in the range 0..9:

In [48]:
def cols(i):
"""Return one of the ten Matplotlib cyclic colors."""
return "C" + str(i%10)


Color handling in matplotlib is actually quite complex, and you can find out more about it here: https://matplotlib.org/users/colors.html

We can now use our graphics functions to demonstrate our randomPointSetNonGP(...) and extremePoints(...) functions. Just move the interactor with the sliders:

In [49]:
#Demo
_rps1 = randomPointSetNonGP(numberOfPoints=10, start=-4, end=4)

@interact(x=(-5,5), y=(-5,5))
def _demoRandomPointSetAndExtremePoints(x=0, y=0):
# Geometry
rps = _rps1.copy()
p = Point(x,y)
rps.append(p)
exts = extremePoints(rps)

# Graphics
cgc = CGCanvas("randomPointSetNonGP(...) and extremePoints(...)")
# All points
circles(cgc, rps)
# Extreme points
circles(cgc, exts, facecolor="red")
# Label the extreme points
text(cgc, exts[RIGHTMOST]+Point(0,0), text='Rightmost    ', horizontalalignment='right')
text(cgc, exts[TOPMOST]+Point(0,0.5), text='Topmost')
text(cgc, exts[LEFTMOST]+Point(0,0), text='    Leftmost', horizontalalignment='left')
text(cgc, exts[BOTTOMMOST]+Point(0,-0.5), text='Bottommost')
# Interactive point
interactor(cgc, p, dot=False, fill=None)

cgc.show()
#Demo


# Lines, rays, and line segments¶

Contents

There are quite a few definitions related to lines which we have indented to indicate taxonomic relationships:

• Line: a straight one-dimensional figure having no thickness and extending infinitely in both directions.
• Ray: a straight one-dimensional figure having no thickness, beginning at a point and extending infinitely in one direction.
• Line segment: a straight one-dimensional figure having no thickness, beginning at a point and ending at a point.
• Euclidean Plane: A flat, two dimensional surface.
• Half plane: a planar region consisting of all points on one side of an infinite straight line, and no points on the other side.
• Closed half plane: a half plane that includes the points on the line.
• Open half plane: a half plane that does not include the points on the line.
• General position for lines: a set of lines is said to be in general position if no three lines are concurrent (intersect at the same point).
• Vertex: The common endpoint of two or more rays or line segments.

For most computational geometry algorithms, the direction of a line, ray or line segment is important, so we will generally treat all of these as being directed where the direction is from the first point to the last point.

## The geometry of line segments¶

Contents

Once we have 2 points, $p_1$ and $p_2$, forming a line segment, $(p_1, p_2)$, interesting geometry begins to emerge. Consider the diagram below:

In [50]:
#Demo
def _demoGeometryOfLineSegments():
# Geometry
p1 = Point(1,1)
p2 = Point(5,4)
offset = Point(0.4, 0.4)

# Graphics
cgc = CGCanvas("Geometry of line segments", p1=Point(0,0), p2=Point(7,6))

circle(cgc, p1)
circle(cgc, p2)

annotation(cgc,p1,p2, arrowstyle="<-")

lineFromTo(cgc,p1,Point(p2.x, p1.y), color="red")
lineFromTo(cgc,p2,Point(p2.x, p1.y), color="red")

text(cgc, p1 - offset, text="p1")
text(cgc, p2 + offset, text="p2")
text(cgc, Point(5.5,2.5), text="o")
text(cgc, Point(3,0.5), text="a")
text(cgc, Point(3,3), text="h")
text(cgc, Point(2,1.2), text='$\\theta$')

cgc.show()

_demoGeometryOfLineSegments()
#Demo


From the right angled triangle formed by $p_1$ and $p_2$ we have the adjacent, $a$, the opposite, $o$, and the hypotenuse, $h$.

• The angle $\theta$, that the line segment makes with the X axis can be calculated, but generally speaking, the sin and cos of $\theta$ are more useful:

• $sin \theta = \frac{o}{h}$

• $cos \theta = \frac{a}{h}$

• We can calculate the length of $a$ and $o$ by subtraction and/or by projecting $h$ onto the X and Y axies:

• $a$ = p2.x - p1.x = $h \: cos \: \theta$

• $o$ = p2.y - p1.y = $h \: sin \: \theta$

• Pythagoras gives us the length of the hypotenuse: $h = \sqrt{a^2+o^2}$ = length of the line segment

• The midpoint of line segment (p1, p2) is just the average of the two points.

• Subtracting p1 from p2 makes p2 equivalent to a vector from the origin.

• The line segment as a unit vector is just $(cos \: \theta, sin \: \theta)$

• If we extend the length, l of the line segment to L, the new end point, P is given by: P = ( L $cos \: \theta$ + p1.x, L $sin \: \theta$ + p1.y )

Here are some functions to implement this geometry:

In [51]:
def length(p):
"""The length of Point p taken as a vector i.e. the distance from the origin to p."""
return distance(Point(0,0), p)

def vector(p1, p2):
"""The vector from Point p1 to Point p2."""
return p2 - p1

def pSin(p1, p2):
"""The sine of the angle between Point p1 and Point p2."""
return (p2.y - p1.y) / distance(p1, p2)

def pCos(p1, p2):
"""The cosine of the angle between Point p1 and Point p2."""
return (p2.x - p1.x) / distance(p1, p2)

def unitVector(p1, p2):
"""The unit vector from Point p1 to Point p2."""
return Point(pCos(p1, p2), pSin(p1, p2)) if p1 != p2 else Point(0,0)

def midpoint(p1, p2):
"""The Point that is the midpoint between Point p1 and Point p2."""
return p1 + (p2 - p1)/Point(2,2)

def extendLength(p1, p2, l):
"""Extend the line from Point p1 to Point p2 by length l."""
x = l * pCos(p1, p2) + p1.x
y = l * pSin(p1, p2) + p1.y
return [p1, Point(x,y)]


# Displaying lines, rays and line segments¶

Contents

Lines, rays and line segments are defined by exactly two points, $p_1$ and $p_2$:

• A line passes through $p_1$ and $p_2$ and extends to infinity in both directions.

• A rays begins at $p_1$ and passes through $p_2$ to infinity.

• A line segment begins at $p_1$ and ends at $p_2$.

Provided we take "infinity" to mean "outside our plotting area", we can create rays and lines just by extending the length of a line segment in one or both directions outside of our plot. This is very easy to implement:

In [52]:
def lineSegment(cgc, p1, p2, **kwargs):
"""Draw a line segment from Point p1 to Point p2."""
cgc.ax.add_patch(patches.FancyArrow(p1.x, p1.y, p2.x - p1.x, p2.y - p1.y, **kwargs))

def ray(cgc, p1, p2, **kwargs):
"""Draw a ray from Point p1 through Point p2 to infinity (and beyond!)."""
if p1 != p2:
posInf1, posInf2 = extendLength(p1, p2, 1000)
lineSegment(cgc, posInf1, posInf2, **kwargs)

def line(cgc, p1, p2, **kwargs):
"""Draw a line from infinity though Point p1 and Point p2 to infinity."""
if p1 != p2:
posInf1, posInf2 = extendLength(p1, p2, 1000)
negInf1, negInf2 = extendLength(p1, p2, -1000)
lineSegment(cgc, posInf1, posInf2, **kwargs)
lineSegment(cgc, negInf1, negInf2, **kwargs)


Here is a simple demo program:

In [57]:
#Demo
@interact(option = ["Line segment", "Ray", "Line"],x1=(-5,5), y1=(-5,5),x2=(-5,5), y2=(-5,5))
def _demoLines( option = [], x1=1, y1=1, x2=3, y2=4):
# Geometry
p1 = Point(x1,y1)
p2 = Point(x2,y2)
d = distance(p1, p2)
l1 = length(p1)
l2 = length(p2)
vp1p2 = vector(p1, p2)
uvp1p2 = unitVector(p1, p2)

# Graphics
cgc = CGCanvas("Line segments, rays and lines")

text(cgc, offsetPoint(p1), text="p1")
interactor(cgc, p1)

text(cgc, offsetPoint(p2), text="p2")
interactor(cgc, p2)

circle(cgc, midpoint(p1, p2), color="red")

if option == "Line segment":
lineSegment(cgc, p1, p2)
if option == "Ray":
ray(cgc, p1, p2)
if option == "Line":
line(cgc, p1, p2)

text(cgc, Point(-6,-3), text="distance = {}".format(d))
text(cgc, Point(-6,-4), text="length p1 = {}, length p2 ={}".format(l1, l2))
text(cgc, Point(-6,-5), text="vector = {}".format(vp1p2))
text(cgc, Point(-6,-6), text="unit vector = {}".format(uvp1p2))

cgc.show()
#Demo


# Summary¶

Contents

In this chapter we have put in place the basic concepts, objects and tools to allow us to explore computational geometry through simple interactive programs. We will build on this base in subsequent chapters, constructing algorithms, more complex objects and interactive demonstrations.

It is interesting how much work we have put into constructing Point and PointList classes. Most of this has been to make the classes work as seamlessly with Python as possible. In fact, we will see that it often doesn't matter whether we use a PointList object, a Python list or some other type of collection. The code just works. As long as the object in question implements the right interface, it doesn't really matter what it is. This is often known as, "duck typing", based on the phrase, "if it looks like a duck, swims like a duck, and quacks like a duck, then it probably is a duck".