In [32]:
from ICGChapter00 import *
In [33]:
# Uncomment the line below for better graphics if you have a retina display
%config InlineBackend.figure_format = 'retina'

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. For example, some algorithms only work if no two points have the same x or y. We will indicate any such additional constraints as we encounter them.



We will represent a point as a 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 that makes it easy to explore the characteristics of the computational geometry algorithms including 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 [35]:
# 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

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

    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)

    def __add__(self, p):
        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:

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:
    • Addition: __add__
    • 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 [36]:
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))
    x, y = Point(5,7)
    print("Destructuring: the expression x, y = Point(5,7) gives x = {}, y = {}".format(x,y))
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


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 [37]:
def distance(p1, p2):
    """Calculate the distance between two Point objects."""
    return math.hypot(p2.x - p1.x, p2.y - p1.y)

Here is a simple demo program. For ease of identification, and to separate them from the computational geometry algorithms, all demo programs and illustrations are surrounded by comments thus: #Demo ... #Demo. Their cell also has the special DEMO tag.

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

We will generally adopt functional programming style, and implement computational geometry algorithms as functions that operate on simple entity classes such as Point. An entity class is simply a class that has very little behaviour - it is really just a data structure.

A more object oriented approach would mean having large classes with computational geometry algorithms as methods. This has some advantages, but we find that full object orientation doesn't work so well in this case for a couple of reasons:

  1. The functional programming style lends itself very well to mathematical programming. After all, most mathematics is about functions in one form or another, and aggregating these functions into classes really adds very little value.

  2. We need to present the material in this book in a logical, sequential fashion. This means we can't create an all-singing-all-dancing Point class with a full set of computational geometry algorithms, because we haven't introduced those algorithms yet.

We are actually expert object-oriented programmers (see our books at but in this case, the functional style is a better fit. One of the joys of Python is that it is a multi-paradigm language, so we can use whatever programming paradigm is best for a particular problem.

Lists of Points


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

In [39]:
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])
            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."""
    def appendAll(self, plist):
        """Modify the PointList object by appending a list of Point objects."""
        for p in plist:

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

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 [40]:
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))
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 PointList containing the appropriate elements. For example:

In [41]:
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

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 [42]:
def _demoAppendAndRemove():
    pl = PointList((4,4), (5,2), (6,3), (7,4), (8,5))
    pl.appendAll((Point(1,1), Point(2,2)))

    pl.removeAll((Point(1,1), Point(2,2)))

Ideally, we would have liked our PointList class to be immutable. 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.

Be aware that when you append or remove an item from the PointList it changes the PointList! If you don't want that, then use the copy(...) method first to create a deep copy of the PointList. 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


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 [43]:
def _demoIntegerPointList():
    pl = PointList((2.0, 3.141), (5.23, 4.29), (6.7, 8.9))
    pli = pl.roundPoints()

Ordered Point sets


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 Python3 is to use the facilities of the OrderedDict class. We create an OrderedDict using the Points as keys. This automatically removes duplicates. We then convert this back to a list.

In [44]:
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 [45]:
print(OrderedPointSet((1,2), (3,4), (5,6), (1,2)))

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


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 [46]:

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 [47]:
extremePoints([Point(-1,-5), Point(2,-2), Point(-4,-4), Point(2,1), Point(2,2)])

Random sets of Points not in general position


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 [48]:
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 [49]:
[Point(3,0), Point(3,-1), Point(-1,-4), Point(-3,1), Point(0,3)]

Displaying Points and other geometrical objects


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 on which to draw our computational geometry objects. 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.

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:

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

In [50]:
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.fig.add_subplot(111, aspect='equal')
        plt.xticks(range(p1.x, p2.x))
        plt.yticks(range(p1.y, p2.y))[p1.x, p2.x])[p1.y, p2.y])

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

This is what it looks like:

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

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 uses to mark an interactive (movable) point.
  • text(...) - text.
  • annotation(...) - an arrow.
  • lineFromTo(...) - a line between two points.
  • path(...) - a closed polygonal shape that is usually filled.

Here are the implementations:

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

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.""", pt.y), radius, fill=fill, **kwargs))
    if dot: circle(cgc, pt, radius = 0.15, fill=True, **kwargs) 
def text(cgc, pt, text, **kwargs):
    """Draw text.""", pt.y, text, **kwargs)

def annotation(cgc, pt1, pt2, arrowstyle="->", text="", **kwargs):
    """Draw an annotation comprising a line with an arrow head.""", 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.""", 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 = [Path.MOVETO] + pathCodes + [Path.CLOSEPOLY], pathCodes), **kwargs))

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 [53]:
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 [54]:
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 [55]:
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:

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

In [56]:
_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)
    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)

Lines, rays, and line segments


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 treat all of these as being directed. The direction is from the first point to the last point.

Representing line segments


We will represent a line segment as an LineSegment class that is a subclass of PointList that holds exactly two Points.

In [57]:
class LineSegment(PointList):
    def __init__(self, p1, p2):
        """A line segment between two Point objects, p1 and p2."""
        PointList.__init__(self, p1, p2)

    def p1(self):
        """Return Point p1."""
        return self.points[0]

    def p2(self):
        """Return Point p2."""
        return self.points[1]

The geometry of line segments


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 [58]:
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$')

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, sin)$

  • 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 [59]:
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


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

  • Lines pass through $p_1$ and $p_2$ and extend to infinity in both directions.

  • Rays begin at $p_1$ and pass through $p_2$ to infinity.

  • Line segment begin at $p_1$ and end 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 [60]:
def lineSegment(cgc, p1, p2, **kwargs):
    """Draw a line segment from Point p1 to Point p2.""", 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 [61]:
@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))



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