from ICGChapter02 import *
# Comment out the line below if the graphics do not display
%config InlineBackend.figure_format = 'retina'
%%javascript
IPython.OutputArea.prototype._should_scroll = function(lines) {return false;}
In this chapter, we will start to work with triangles, and see how they are an essential tool in computational geometry.
displayTOC("ICGChapter03.ipynb")
Because triangles are three sided polygons, we need to consider a few basic principles of polygons first. We will consider polygons in much more detail later. Here are some definitions:
Considering the definition of a polygon given above:
If non-consecutive edges do not intersect, then it is a simple closed curve and therefore a simple polygon: $$e_i \cap e_j = 0 \: for \: j \ne i+1$$
A convex polygon has an interior that is a convex set. In other words for a convex polygon, every edge between two vertexes must be inside the polygon. For this to be the case, every internal angle must be $\le 180^{\circ}$. We won't provide a proof of this, but it means that convex polygons are always simple. If a polygon is not convex, then it is concave. Concave polygons are not necessarily simple.
As we said above, a triangle is a polygon with three sides. It is simple and convex, because intersections and concavities are impossible when there are only three sides.
CONVEX_HEXAGON_PLIST = PointList((4, 0), (2, 4), (-2, 4), (-4, 0), (-2, -4), (2, -4))
CONCAVE_HEXAGON_PLIST = PointList((4, 0), (2, 4), (-2, 4), (0, 0), (-2, -4), (2, -4))
SQUARE_PLIST = PointList((2, -2), (2, 2), (-2, 2), (-2, -2))
This is what they look like:
#Demo
@interact(poly=["Convex hexagon", "Concave hexagon", "Square"])
def _demoPolygons(poly="Convex hexagon"):
# Geometry
if poly == "Convex hexagon":
pol = CONVEX_HEXAGON_PLIST
if poly == "Concave hexagon":
pol = CONCAVE_HEXAGON_PLIST
if poly == "Square":
pol = SQUARE_PLIST
# Graphics
cgc = CGCanvas("Example closed polygonal chains")
polygon(cgc, pol, linewidth=2)
labelPointList(cgc, pol, template="$v_{}$")
cgc.show()
#Demo
We will adopt the convention that polygon vertexes are always listed in an anti-clockwise order. In other words, if you step from one vertex to the next, the interior of the polygon is always on your left hand side. This anti-clockwise convention is very important, as we will see shortly when we come to calculate polygon areas.
Now that we know how to represent polygons we can begin to investigate their properties. We will start with the simplest polygon - the triangle.
def triangle(cgc, plist, labels=("a", "b", "c"), dot=False, colors=("red", "green", "blue"), **kwargs):
"""Draw a triangle and label its vertexes."""
if dot: circles(cgc, plist)
a,b,c = plist
lineSegment(cgc, a, b, color=colors[0], **kwargs)
lineSegment(cgc, b, c, color=colors[1], **kwargs)
lineSegment(cgc, c, a, color=colors[2], **kwargs)
text(cgc, offsetPoint(a), text=labels[0], color=colors[0])
text(cgc, offsetPoint(b), text=labels[1], color=colors[1])
text(cgc, offsetPoint(c), text=labels[2], color=colors[2])
For example:
#Demo
def _demoTriangle():
# Geometry
t = PointList((0,3), (0,0), (4,0))
# Graphics
cgc = CGCanvas("Triangle")
triangle(cgc, t, dot=False)
cgc.show()
_demoTriangle()
#Demo
A key function in computational geometry is the function that gives the signed area of a 2D triangle that is positive if the vertexes are given in an anti-clockwise order, and negative if they are given in a clockwise order.
Because this function is so important, we will spend some time developing it. We will start in 3D, because that is the basis of the 2D function, and then define a flexible and efficient 2D function.
Given a triangle, $T$, in 3D space with vertexes $a$, $b$ and $c$, we can convert the edges $(a,b)$ and $(a,c)$ to edge vectors simply by subtracting $a$. This translates the triangle to put $a$ at the origin. So we have:
$\overrightarrow{v} = b - a$
$\overrightarrow{w} = c - a$
The area of the triangle is then given by:
$$A(T_{3D}) = \frac{|\overrightarrow{v} \times \overrightarrow{w}|}{2}$$
A proof of this is beyond this book, but any text on coordinate geometry should explain it.
We need a Point3D class to help us with this. This class is the pretty much the same code as our standard Point class but with an extra $z$ coordinate. We don't need all of this code for our purposes here, but we have included it for completeness.
class Point3D:
"""A class to represent a point in 3D space with x, y and z coordinates. Point3D objects are immutable."""
def __init__(self, x=0, y=0, z=0):
self.__x = x
self.__y = y
self.__z = z
@property
def x(self):
"""Return the x coordinate."""
return self.__x
@property
def y(self):
"""Return the y coordinate."""
return self.__y
@property
def z(self):
"""Return the z coordinate."""
return self.__z
def __eq__(self, p):
return (self.x == p.x) and (self.y == p.y) and (self.z == p.z)
def __ne__(self, p):
return (self.x != p.x) or (self.y != p.y) or (self.z != p.z)
def __add__(self, p):
return Point(self.x + p.x, self.y + p.y, self.z + p.z)
def __sub__(self, p):
return Point3D(self.x - p.x, self.y - p.y, self.z - p.z)
def __mul__(self, p):
return Point(self.x * p.x, self.y * p.y, self.z * p.z)
def __truediv__(self, p):
return Point(self.x/p.x, self.y/p.y, self.z/p.z)
def __hash__(self):
return hash((self.x, self.y, self.z))
def __iter__(self):
return iter((self.x, self.y, self.z))
def __str__(self):
return "({},{},{})".format(self.x, self.y, self.z)
def __repr__(self):
return "Point({!r},{!r},{}!r)".format(self.x, self.y, self.z)
We can now define a function for the area of a triangle in 3D as follows:
def triangleArea3D(a, b, c):
"""Calculate the signed area of the triangle defined by the Point3D objects a, b, and c. Positive if points ant-clockwise."""
v = b - a
w = c - a
cross = np.cross((v.x, v.y, v.z), (w.x, w.y, w.z))
return np.linalg.norm(cross)*0.5
Of course, we don't have to use built-in functions, we can just multiply everything out as follows:
$$\overrightarrow{v} \times \overrightarrow{w} = -a_z b_y + a_y b_z + a_z c_y-b_z c_y - a_y c_z + b_y c_z, a_z b_x - a_x b_z - a_z c_x + b_z c_x + a_x c_z - b_x c_z, -a_y b_x + a_x b_y + a_y c_x - b_y c_x - a_x c_y + b_x c_y$$
No - I didn't do that by hand, I used a symbolic algebra system - Mathematica to be precise.
This gives:
$$A(T_{3D}) = \frac{1}{2}\sqrt{(-a_y b_x + a_x b_y + a_y c_x - b_y c_x - a_x c_y + b_x c_y)^2 + (a_z b_x - a_x b_z - a_z c_x + b_z c_x + a_x c_z - b_x c_z)^2 + (-a_z b_y + a_y b_z + a_z c_y - b_z c_y - a_y c_z + b_y c_z)^2}$$
For 2D triangles, T, where the z coordinate is 0, this simplifies quite a bit because we lose all of the terms involving z:
$$A(T) = \frac{1}{2}\sqrt{(-a_y b_x + a_x b_y + a_y c_x - b_y c_x - a_x c_y + b_x c_y)^2} = \frac{1}{2}(-a_y b_x + a_x b_y + a_y c_x - b_y c_x - a_x c_y + b_x c_y)$$
We can implement this function as follows:
def triangleArea(a,b,c):
"""Calculate the signed area of the triangle defined by the Points a, b, and c. Positive if points ant-clockwise."""
return 0.5* (-a.y * b.x + a.x * b.y + a.y * c.x - b.y * c.x - a.x * c.y + b.x * c.y)
Another way to express the $A(T_{3D})$ is as half the determinant of a matrix comprising the $x$, $y$ and $z$ coordinates of the points:
$$A(T_{3D}) = \frac{1}{2} \left | \begin{bmatrix} a _x & a_y & a_z\\ b_x & b_y & b_z \\ c_x & c_y & c_z \end{bmatrix} \right |$$
For 2D points, the $z$ coordinates are set to 1:
$$A(T) = \frac{1}{2} \left | \begin{bmatrix} a _x & a_y & 1\\ b_x & b_y & 1 \\ c_x & c_y & 1 \end{bmatrix} \right |$$
Using the NumPy library this is trivial to implement:
#Demo
def triangleAreaDet3D(a, b, c):
return np.linalg.det([[a.x, a.y, a.z], [b.x, b.y, b.z], [c.x, c.y, c.z]])
def triangleAreaDet(a, b, c):
return np.linalg.det([[a.x, a.y, 1], [b.x, b.y, 1], [c.x, c.y, 1]])
#Demo
But which version is faster? We will look at the 2D case only, because that is our main interest here, but the 3D case will be similar. We can ask Python to profile the functions for us as follows:
# %prun for i in range(10000) : triangleArea(Point(0, 0), Point(1, 0), Point(1, 1)) # time = 0.120
# %prun for i in range(10000) : triangleAreaDet(Point(0, 0), Point(1, 0), Point(1, 1)) # time = 0.514
If you want to see the results on your own system, you can uncomment the two cells above and run them.
The result is that triangleArea(...) is just over four times faster than triangleAreaDet(...). This isn't surprising, because even though triangleAreaDet(...) uses the hugely efficient NumPy mathematics library, it still makes a lot of expensive function calls, whereas triangleArea(...) only uses basic arithmetic.
Below is an interactive demo program for triangleArea(...). Notice how the triangle area is positive when the points are ordered anti-clockwise and negative when they are ordered clockwise.
#Demo
@interact(ax=(-5,5), ay=(-5,5), bx=(-5,5), by=(-5,5), cx=(-5,5), cy=(-5,5))
def _demoTriangleArea(ax=0, ay=3, bx=0, by=0, cx=4, cy=0):
# Geometry
a = Point(ax,ay)
b = Point(bx, by)
c = Point(cx, cy)
area = triangleArea(a, b, c)
# Graphics
cgc = CGCanvas("Triangle area, a to b to c = {}".format(area))
triangle(cgc, PointList(a,b,c))
cgc.show()
#Demo
Now that we have a function to determine the signed area of a triangle, we can use this function to develop some fundamental algorithms of computational geometry - those that determine the relationships of points to lines.
If we treat a triangle as a line segment $(a, b)$ and a point $c$, then the signed triangle area tells us about the relationship between point $c$ and the line $L$ determined by the line segment:
If triangle $(a, b, c)$ area is $\gt 0$, then point $c$ is to the left of line $L$ as you travel from $a$ to $b$. Otherwise it is to the right.
If triangle $(a, b, c)$ area is $\ge 0$ then point $c$ is to the left of the line $L$ or collinear with it.
If triangle $(a, b, c)$ area is $= 0$ then point $c$ is collinear with the line $L$.
You can easily verify this for yourself by playing about with the, "Triangle area", demonstration in the previous section.
def isCollinear(a, b, c):
"""Return True if Points a, b and c are collinear."""
return triangleArea(a, b, c) == 0
def isLeft(a, b, c):
"""Return True if Point c is to the left of line (a, b) but not on it."""
return triangleArea(a, b, c) > 0
def isLeftOn(a, b, c):
"""Return True if Point c is to the left of or on the line (a, b)."""
return triangleArea(a, b, c) >= 0
def isRight(a, b, c):
"""Return True if Point c is to the right of line (a, b) but not on it."""
return not isLeftOn(a, b, c)
def isRightOn(a, b, c):
"""Return True if Point c is to the right of or on the line (a, b)."""
return isLeft(a, b, c)
There are a couple of things to note:
We clearly distinguish between the line $L$ defined by the line line segment $(a, b)$ and the line segment itself. Lines have infinite length, so point $c$ does not need to be on the line segment between $a$ and $b$ for these functions to return True - they work over the infinite length of the line.
It is crucial to understand what "left" actually means in this context: if you are traveling along the line from $a$ to $b$, then "left" means on your left hand side.
Another important relationship we need to test for is "between" - is point $c$ is between $a$ and $b$? By "between", we mean that $c$ is somewhere on the line segment $(a, b)$. This breaks down into two considerations:
For $c$ to be between $a$ and $b$ it is necessary (but not sufficient) that it is collinear with $a$ and $b$.
The point $c$ must not be off either end of the line segment $(a, b)$. In other words, $c_x$ must be between $a_x$ and $b_x$ AND $c_y$ must be between $a_y$ and $b_y$.
We can express this second constraint as follows:
$((a_x \le c_x \le b_x) \lor (b_x \le c_x \le a_x)) \land ((a_y \le c_y \le b_y) \lor (b_y \le c_y \le a_y))$
Where $\land$ means logical and, and $\lor$ means logical or.
This translates directly into Python code as follows:
def isBetween(a, b, c):
"""Return True if Point c is somewhere on the line segment (a, b)."""
if not isCollinear(a, b, c):
return False
return ((a.x <= c.x <= b.x) or (b.x <= c.x <= a.x)) and ((a.y <= c.y <= b.y) or (b.y <= c.y <= a.y))
#Demo
@interact(cx=(-5,5), cy=(-5,5))
def _demoPointCBetweenAAndB(cx=0, cy=3):
# Geometry
c = Point(cx, cy)
a = Point(-4,0)
b = Point(4, 0)
between = isBetween(a, b, c)
collinear = isCollinear(a, b, c)
leftOn = isLeftOn(a, b, c)
left = isLeft(a, b, c)
# Graphics
cgc = CGCanvas("isBetween = {}, isCollinear = {}, isLeftOn = {}, isLeft ={}".format(between, collinear, leftOn, left))
if between:
col = cols(1)
elif left:
col = cols(2)
elif leftOn:
col = cols(3)
elif collinear:
col = cols(4)
else:
col = cols(5)
circles(cgc, [a, b])
text(cgc, pt=a + Point(0,-0.7), text="a")
text(cgc, pt=b + Point(0,-0.7), text="b")
interactor(cgc, c, color=col)
cgc.show()
#Demo
Now that we have a test for collinearity we can write a test to see if a set of points is in general position. The algorithm is simple: take all subsets of the point set that have cardinality 3, and if none of these is collinear, then the point set is in general position.
In Python, we can generate all subsets of a set of $n$ points of size 3 as follows:
#Demo
def subsets3(n):
subsets = []
for i in range(n):
for j in range(i+1, n):
for k in range(j+1, n):
subsets.append([i, j, k])
return subsets
subsets3(4)
#Demo
However, we won't use this approach, because Python has the itertools library. This contains itertools.combinations(...), which does the same thing, but faster:
#Demo
list(itertools.combinations(range(4), 3))
#Demo
We can use this to write a function to check a point set for general position:
def isGeneralPosition(plist):
"""Return True if a list of Points is in general position."""
perms = itertools.combinations(plist, 3)
for p in perms:
if isCollinear(p[0], p[1], p[2]):
return False
return True
Here is a demo program. Move the points around and see under what conditions they go into and out of general position:
#Demo
@interact(p0x=(-5,5), p0y=(-5,5), p1x=(-5,5), p1y=(-5,5), p2x=(-5,5), p2y=(-5,5), p3x=(-5,5), p3y=(-5,5))
def _demoPointsInGeneralPosition(p0x=0, p0y=3, p1x=0, p1y=0, p2x=4, p2y=0, p3x=3, p3y=0):
# Geometry
plist = PointList((p0x, p0y), (p1x, p1y), (p2x, p2y), (p3x, p3y))
gp = isGeneralPosition(plist)
# Graphics
cgc = CGCanvas("Points in general position is {}".format(gp))
col = "green" if gp else "red"
circles(cgc, plist, color=col)
labelPointList(cgc, plist)
cgc.show()
#Demo
At this point it is interesting to revisit our randomPointSetNonGP(...) function with the question, what is the probability of generating a random point set in general position? If we were using floating point numbers, then the answer would be almost certainty, because there are so many possibilities to choose from. If we are using integers on a grid, then the answer depends on the size of the grid, and how many points we ask for.
We can perform an experiment: If the size of the grid is $n$ x $n$ and we ask for a set of $m$ random points, we can calculate the probability they are in general position by running randomPointSetNonGP(...) many times and recording what proportion of the runs give a result in general position. Here is a function to do this:
def GPProbability(n, m, runs=500):
"""Return the probability of a Point set generated with randomPointSetNonGP(...) is in general position."""
gp = 0
for r in range(runs):
if isGeneralPosition(randomPointSetNonGP(m, start=1, end=n)):
gp += 1
return gp/runs
For example, if we have $n$ = 11 (our standard grid) and ask for 10 points, the probability that these 10 points will be in general position is:
#Demo
GPProbability(11, 10)
#Demo
This is very low! Let's plot probability of general position vs. $\frac{m}{n}$ for $n = 11$, to see what is going on. Here is a function to generate a probability curve for $m = 0$ to $m = mMax$:
def GPProbabilityCurve(n, mMax, runs=500):
"""Probability curve for a set of m points chosen from a square grid of n points being in general position."""
x = []
y = []
for m in range(mMax):
x.append(m/n)
y.append(GPProbability(n, m, runs=500))
return [x, y, n]
Here is the curve for $n = 11$ and $m=0$ to $m=15$:
#Demo
x, y, n = GPProbabilityCurve(11, 15)
plt.plot(x, y)
plt.ylabel("Probability of GP for n = {}".format(n))
plt.xlabel(r"Ratio $\frac{m}{n}$")
plt.show()
#Demo
You can see from the above graph for $n = 11$, that as $m$ approaches $n$, and $\frac{m}{n}$ approaches 1, the probability of the point set being in general position gets very low. Rather than take values off the above graph, we can calculate $m$ for a given probability $P(n)$:
def mForPn(P, n):
"""Calculate the size of a Point set drawn from a square grid of n Points that has a probability P of being in general position."""
m = 0
while GPProbability(n, m) > P:
m += 1
return m
Use the demo program below to find $m$ for a given $P(n)$. You can see that for $n = 11$ and a $P(n)$ of 0.5 (a 50:50 chance) you can ask for no more than 7 points:
#Demo
@interact(P=[0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9], n=range(3, 15))
def _demoProbs(P=0.5, n=11):
print(r"For n = {}, the number of points, m, for probability of {} = {}".format(n, P, mForPn(P, n)))
#Demo
We can also explore the effect of increasing $n$:
#Demo
curves = [GPProbabilityCurve(n, 2*n, runs=500) for n in range(1, 14, 2)]
for c in curves:
x, y, n = c
plt.plot(x, y, label="n={}".format(n))
plt.ylabel("Probability of GP")
plt.xlabel(r"Ratio $\frac{m}{n}$")
plt.legend()
plt.show()
#Demo
The above graph shows that as $n$ increases, the ratio $\frac{m}{n}$ for which you can get any given probability of a general position set decreases.
You can see that this simple experiment has allowed us to get a much clearer idea of how likely the function RandomPointSetNonGP(...) will give a point set in general position. It is often necessary to perform experiments such as this in order to understand the behavior of certain functions.
In this section we will use the function isGeneralPosition(...), to create a simple constructive algorithm to generate random sets of points in general position. As usual, we will generate points on an $n$ x $n$ grid, so the maximum number of points we can have is obviously $n^2$, in which case there is a point on each intersection of the grid.
Firstly, we need to work out the maximum number of points it is possible to have in general position for a given grid.
The no-three-in-line problem asks for the maximum number of points that can be placed on an $n$ x $n$ grid so that no three points are collinear. This number is at most $2n$, since if $2n + 1$ points are placed in the grid, then by the pigeonhole principle some row and some column must contain three points. Given this result, on our usual $(-5, -5)$ to $(5, 5)$ grid, which is 11 x 11, the maximum number of general position points is 11 x 2 = 22.
Here is a simple constructive algorithm to generate a point set $P$ of $m$ random points in general position where $m \lt 2n$ and $n$ is the side of the grid:
Let $P = \emptyset$ (the empty set)
While $|P| \lt m$
2.1 Generate a random point, $p$.
2.2 If $P \cup \{p\}$ is in general position then $P = P \cup \{p\}$
The algorithm generates random points and adds them to the set if they are in general position until we have the number we want. To implement this, we first of all need a function to tell us if a point is in general position with respect to a point set:
def isPointInGP(pts, p):
"""Return True if p is in general position with respect to Point set pts."""
ptsGP = copy.deepcopy(pts)
ptsGP.append(p)
return (p not in pts) and isGeneralPosition(ptsGP)
Next, we simply add points in general position to get our point set:
def randomPointSetGP(m, start=-5, end=5):
"""Return a random set of m points in general position."""
pts = []
m = min(m, int((end - start + 1) + (end - start + 1)/2))
while len(pts) < m:
p = Point(random.randint(start, end), random.randint(start, end))
if isPointInGP(pts, p):
pts.append(p)
return pts
Notice that for a given $n$ x $n$ grid we have set the maximum number of points in general position to be the integer value of $n + \frac{n}{2}$. This is a purely pragmatic value based on our experience with the function. The problem is that if we allowed the maximum possible number of points ($2n$) we might never be able to get there because of how the early random points are positioned, and this would lead to an infinite loop. The value of $n + \frac{n}{2}$ is nice compromise that is low enough that the algorithm always terminates (note: we do not have a proof for this) and we get a good random scattering of points that doesn't look too constrained.
Here is a demo program. Notice that the function limits the number of points as we just described:
#Demo
@interact(m=(0,20))
def _demoRandomPointSet(m):
# Geometry
pts = randomPointSetGP(m)
# Graphics
cgc = CGCanvas("Random point set in general position for m = {}".format(len(pts)))
circles(cgc, pts)
cgc.show()
#Demo
Sometimes all we need to know is what side of a line a set of points is on. Another way of saying this is whether a set of points all lie in the same open half plane (points on the line not included) or closed half plane (points on the line included). These functions are trivial given the functions we developed above:
def isSameSide(a, b, plist):
"""Return True if all Points in plist are on the same side (but not on) line (a,b)."""
left = right = True
for p in plist:
left = left and isLeft(a, b, p)
right = right and isRight(a, b, p)
return left or right
def isSameSideOn(a, b, plist):
"""Return True if all Points in plist are on the same side or on line (a,b)."""
leftOn = rightOn = True
for p in plist:
leftOn = leftOn and isLeftOn(a, b, p)
rightOn = rightOn and isRightOn(a, b, p)
return leftOn or rightOn
Here is a demo program:
#Demo
@interact(p0x=(-6, 6), p0y=(-6,6), p1x=(-6, 6), p1y=(-6,6), p2x=(-6, 6), p2y=(-6,6))
def _demoIsSameSide(p0x=5, p0y=-3, p1x=4, p1y=3, p2x=-3, p2y=5):
# Geometry
plist = PointList((p0x, p0y), (p1x, p1y), (p2x, p2y))
a = Point(-5, 5)
b = Point(5, -5)
sameSide = isSameSide(a, b, plist)
sameSideOn = isSameSideOn(a, b, plist)
# Graphics
cgc = CGCanvas("Same side = {}, same side on = {}".format(sameSide, sameSideOn))
line(cgc, a, b)
if sameSide:
col = "green"
elif sameSideOn:
col = "blue"
else:
col = "red"
for p in plist:
interactor(cgc, p, color=col)
cgc.show()
#Demo
Triangles are the key to computational geometry. They allow us to discover relationships between points and lines by forming them into triangles and considering the signed triangle areas. These relationships are the building blocks for the more complex algorithms we will introduce in subsequent chapters. Always remember that for our algorithms to work, we need to adhere to our convention that polygon vertexes always go in an anticlockwise order.