In [1]:

```
from ICGChapter02 import *
```

In [2]:

```
# Uncomment the line below for better graphics if you have a retina display
%config InlineBackend.figure_format = 'retina'
```

In this chapter, we will start to work with triangles, and see how they are an essential tool in computational geometry.

In [3]:

```
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:

**Polygon**: the region of a plane bounded by a closed curve. A polygon is a closed planar polygonal chain.**Simple polygon**: the region of a plane bounded by a simple (non self-intersecting closed curve).**Convex polygon**: a simple polygon whose interior is a convex set. A convex polygon is monotone with respect to any straight line.**Triangle**: a three-sided simple, convex polygon.

**Monotone polygon**: a polygon is monotone with respect to a straight line L if every line perpendicular to L intersects the polygon at most twice.

**Convex set**: a set of points, $p_1, p_2, ...p_n$ such that all edges, $(p_i, p_j)$ where $i \le n$, $j \le n$ and $i \ne j$ lie within the set.

Considering the definition of a polygon given above:

- A polygon has $n$ consecutive vertexes $v_1, v_2, ... v_n$. The edges of the polygon $e_1, e_2, ... e_n$ are the line segments $(v_1, v_2), (v_2, v_3), ... (v_n, v_1)$
- The edges are connected end to end to form a curve. In other words, the intersection of $e_i$ and $e_{i+1}$ is the vertex $v_i$: $$e_i \cap e_{i+1} = v_{i+1} \: for \: i=1..n$$
- The last vertex of the last edge is connected back to the first vertex of the first edge (a closed curve). An alternate definition (that we don't use in this text) is that the first and last vertexes are identical. $$e_n \cap e_1 = v_1$$

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.

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.

In [4]:

```
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:

In [5]:

```
#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.

In [6]:

```
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:

In [7]:

```
#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 moves 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 (REF).

We need a Point3D class to help us with this. This class is 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.

In [8]:

```
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:

In [9]:

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

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:

In [10]:

```
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:

In [11]:

```
#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:

In [12]:

```
# %prun for i in range(10000) : triangleArea(Point(0, 0), Point(1, 0), Point(1, 1)) # time = 0.120
```

In [13]:

```
# %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.

Here is an interactive demo program for triangleArea:

In [14]:

```
#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 very fundamental algorithms of computational geometry - those that determine the relationships of points to lines.

If we treat a triangle as an edge $(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 edge:

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.

In [15]:

```
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 edge $(a, b)$ and the edge itself. Because lines have infinite length, point $c$ does not need to be 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 edge $(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)$ OR $(b_x \le c_x \le a_x))$ AND $((a_y \le c_y \le b_y)$ OR $(b_y \le c_y \le a_y))$

This translates directly into Python code as follows:

In [16]:

```
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))
```

In [17]:

```
#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 so we can write a test to see if a set of points are 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:

In [18]:

```
#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
```

Out[18]:

However, we won't use this approach, becaue Python has the itertools library. This contains itertools.combinations(...), which does the same thing, but faster:

In [19]:

```
#Demo
list(itertools.combinations(range(4), 3))
#Demo
```

Out[19]:

We can use this to write a function to check a point set for general position:

In [20]:

```
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 to test if a set of points is in general position:

In [21]:

```
#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.

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 general position result. Here is a function to do this:

In [22]:

```
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:

In [23]:

```
#Demo
GPProbability(11, 10)
#Demo
```

Out[23]:

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$:

In [24]:

```
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..15$:

In [25]:

```
#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)$:

In [26]:

```
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 50:50 you can ask for no more than 7 points:

In [27]:

```
#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$:

In [28]:

```
#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.

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$, where there is a point on each intersection of the grid. We need to consider the maximum number of points it is possible to have in general position.

The no-three-in-line problem asks for the maximum number of points that can be placed in the $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:

In [29]:

```
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:

In [30]:

```
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:

In [31]:

```
#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:

In [32]:

```
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:

In [33]:

```
#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
```

There is no doubt that triangles are key to computational geometry. They allow us to discover relationships between points by forming them into triangles and considering the signed triangle areas. These relationships are the building blocks for more the complex algorithms we will introduce later. We should also point out that in order for all this to work, we need to adhere to our convention that polygon vertexes always go in an anticlockwise order.