In a previous article I introduced basic vector concepts. In this article, I show how apply the theory to common geometric problems.

About the Implementations

The implementations given here are written to illustrate some basic concepts. It is not production code, so it is not the best implementation in terms of maintainability, speed, or completeness. Read the tips at the end that point out some of the limitations (and possible improvements).

The code is written in Python, but should be easily readable to programmers of other languages.

Shape Representation

In most cases you will not choose your own representations of geometric figures, as these are typically determined by the library or game engine you are using. If you do roll out your own geometry primitives, however, you need to consider how you will typically get source data (that is, data formats of 3D models, hard-coded shapes, and scripted shapes), and what the typical operations on these primitives will be. In addition, you should take care that your representations work well together, otherwise your algorithms built on top will be clumsy.

The representations provided here will be used by the recipes.

Points and Vectors

As explained the previous article, there is no need to distinguish points from vectors, thus a single vector can be used to represent a point. Vectors are usually simply represented as a record, class or list of numbers – as many as the number of dimensions.

The following is a simple implementation of a 2D vector and its operations. It is pretty straightforward, so I present the entire block here without much explanation.

class Vector2D ():
  def __init__(self, x, y):
    self.x = x
    self.y = y
 
  def __str__(self):
      return '[%f, %f]' % (self.x, self.y)
 
  def neg(self):
    return Vector2D(-self.x, -self.y)
 
  def add(self, other):
    x_sum = self.x + other.x
    y_sum = self.y + other.y
 
    return Vector2D(x_sum, y_sum)
 
  def sub(self, other):
    x_dif = self.x - other.x
    y_dif = self.y - other.y
 
    return Vector2D(x_dif, y_dif)
 
  def sqr_length(self):
    return self.x * self.x + self.y * self.y
 
  def length(self):
    return sqrt(self.sqr_length())
 
  def scale(self, factor):
    return Vector2D(self.x * factor, self.y * factor)
 
  def duplicate(self):
    return Vector2D(self.x, self.y)
 
  def normalize(self):
    v = self.duplicate()
    v.normalize_inplace()
    return v
 
  def equals(self, other):
    return sqr(self.x - other.x) + sqr(self.y - other.y) < ZERO_DISTANCE_SQR
 
  def perp_op(self):
    return Vector2D(-self.y, self.x)
 
  def dot(self, other):
    return self.x * other.x + self.y * other.y
 
  def perpdot(self, other):
    return - self.y * other.x + self.x * other.y
 
  def proj_on(self, other):
    scale_factor = self.dot(other)/other.sqr_length()
    new_vector = other.scale(scale_factor)
    return new_vector
 
  def perp_on(self, other):
    scale_factor = -(self.perpdot(other)/other.sqr_length())
    new_vector = other.scale(scale_factor)
 
    return new_vector

A 3D Vector implementation is basically the same, with the obvious changes:

• An additional z component is added, and all functions take this extra component into account.
• The perp dot product is replaced with the cross product, like this:

  def cross(self, other): x = self.y*other.z – self.z*other.y y = self.z*other.x – self.x*other.z z = self.x*other.y – self.y*other.x

  def cross(self, other): x = self.y*other.z – self.z*other.y y = self.z*other.x – self.x*other.z z = self.x*other.y – self.y*other.x

• The perp_on method is replaced with this:

  def perp_op(self, other): proj_vector = self.proj_on()   return self.sub(proj_vector)

  def perp_op(self, other): proj_vector = self.proj_on() return self.sub(proj_vector)

Lines

Lines are best represented by a point on the line (that is, remember, a vector from the origin to that point), and a vector in the same direction as the line. The length of this vector does not matter, as long as it is not 0. (The direction vector is usually normalised; it is not done here so that normalisation is explicit in the algorithms that follow, which makes them slightly easier to read.)

class Line():
  def __init__(self, point, direction):
    self.point = point
    self.direction = direction

It is convenient to have constructors or static creation methods that work with other representations. For instance, we might find a method that uses two points to create a line useful:

class Line():
  ...
 
  @staticmethod
  def new_from_points(point1, point2):
    direction = point2.sub(point1)
    return Line(point1, direction)

Line segment, Triangles, and Other Polygons

Line segments, triangles and other polygons can be represented by their vertices.

class Polygon:
  def __init__(self, vertices):
    self.vertices = vertices #list of vertices as vectors.
    self.vertex_count = len(self.vertices)
 
  def get_edge_line(self, index):
    vertex1 = self.vertices[index]
    vertex2 = self.vertices[(index + 1) % self.vertex_count]
    direction = vertex2.sub(vertex1)

For example, the segment from [-1, 1] to [1, -1] will be represented like this:

segment = Polygon([Vector2D(-1, 1), Vector2D(1, -1)])

It is useful to use synonyms in languages that allow them to help our thinking. In Python, we could simply use this line after the Polygon class definition:

Segment = Polygon

This allows us to construct segments like this:

segment = Segment([Vector2D(-1, 1), Vector2D(1, -1)])

which reads better.

In C++, we could use a typedef statement, in C we could use a macro, and so on. Remember that this will not make the compiler enforce the number of vertices–the synonym is merely for making the code more readable.

This implementation of polygon has a drawback that becomes apparent when you try to see whether two polygons are the same. A line segment has two representations (one with the vertices inverted), a triangle has six representations (all the permutations of the vertices), and in general a n-agon has 2nrepresentations. So an equality check has to check all possible permutations, which is not ideal if the check is performed often. Instead, you might want to normalise the representation, so that every polygon has exactly one representation. One way of doing it is to identify a unique vertex (for instance, by using the first by lexicographical ordering – assuming there are no duplicate vertices), and then take other vertices clockwise.

This implementation also allows construction of polygons with degenerate edges (that is, where consecutive vertices are the same), polygons with redundant vertices (three vertices in a straight line), and self-intersecting polygons. Many algorithms depend on polygons not to be one of those, so you might want to check for these conditions at construction time. In addition, if you are working in 3D, you will also want to check that all the points lie on the same plane.

Circles and Spheres

Circles and spheres are usually represented by the point of the centre, and a radius.

def distance_of_point_from_line(point, line):
  line_to_point = point.sub(line.point)
  return line_to_point.perp_on(line.direction).length()

Distances

What is the distance between two points?

The distance is simply the length of the vector between the points:

def dist(v1, v2):
  return v1.sub(v2).length()

When comparing vectors, it is more efficient to calculate the squared lengths, and compare these values, which is why the Vector2D class has a sqr_length() method.

What is the distance between this point and line?

To calculate the distance from a point to a line, we simply calculate the length of the perpendicular to the line, using the line’s direction vector.

def distance_of_point_from_line(point, line):
  line_to_point = point.sub(line.point)
  return line_to_point.perp_on(line.direction).length()
 
def sqr_distance_of_point_from_line(point, line):
  line_to_point = point.sub(line.point)
  return line_to_point.perp_on(line.direction).sqr_length()

Are these two lines parallel?

Lines are parallel when their direction vectors are parallel. We can check this by comparing the normalised direction vectors. Because any line has two possible direction vectors, we need to make two comparisons, inverting the one vector in the second comparison.

def is_parallel1(line1, line2):
  v = line1.direction.normalize()
  w = line2.direction.normalize()
 
  return v.equals(w) or v.equals(w.neg())

In 2D, another way is to compute the perp dot product of the direction vectors:

def is_parallel(line1, line2):
  v = line1.direction
  w = line2.direction
 
  return abs(v.perpdot(w)) < 0.001

The last implementation avoids making a double comparison and doing expensive normalisations. It is, however, somewhat incorrect in that the result will depend on the vector lengths (even though vector lengths should not affect whether vectors are parallel or not). See the discussion under the tip Watch out for special cases below.

Are these two lines perpendicular?

Lines are perpendicular when their direction vectors are perpendicular. We can check this by computing the dot product of the direction vectors.

def is_perpendicular(line1, line2):
  v = line1.direction
  w = line2.direction
 
  return abs(v.dot(w)) < 0.001

As with checking for parallel lines, the above implementation depends on vector lengths.

Are these two points on the same side of this line?

Here we want to determine whether a point is on the same side of a line as some reference vector (which is perpendicular to the line).

We do this by computing the dot product of the reference vector, and the vector from the line point to the point under consideration. If this value is positive, the cosine of the angle between these vectors is positive, and hence the angle must lie between -180º and 180º, that is, the vectors must point to the same side of the line.

def is_point_on_same_side_of_line(line, point1, point2):
  v1 = point1.sub(line.point)
  v2 = point2.sub(line.point)
 
  if line.direction.perpdot(v1)*line.direction.perpdot(v2) > 0:
    return True

Does this line lie between these two points?

This is a variation of the previous recipe.

def is_line_between_points(line, point1, point2):
  proj1 = point1.sub(line.point).perpdot(line.direction)
  proj2 = point2.sub(line.point).perpdot(line.direction)
 
  return  proj1 * proj2 < 0

Bisecting a Line Segment

The point that lies halfway between points \(\bv v\) and \(\bv w\) is given by \((\bv v + \bv w) / 2\). With our line segment representation, this calculation becomes:

def bisect_segment(segment):
  center = segment.vertices[0].add(segment.vertices[1]).scale(0.5)

Angles

Bisecting an Angle

A vector that lies halfway between vectors \(\bv v\) and \(\bv w\) is given by \(\bv v/|\bv v| + \bv w/|\bv w|\).

def bisect_angle(v1, v2):
  return v1.normalize().add(v2.normalize())

This works as long as the angle is not 180º exactly. When the angle is 180º, the algorithm returns the zero vector. In this case, there are two possible solutions – we can select one arbitrarily, or one on the same side of an additional reference vector.

Is this object in this viewing cone?

Knowing when an object is inside some viewing cone is useful for culling. We can easily do a check by considering the angle between the vector from the cone vertex to the object, and the direction vector of the cone. Moreover, we need only compare the cosines of the angles, not the angles themselves – thus we can use the dot product for an efficient test.

This test does not take the size of objects into account, and is thus more suitable for AI culling than graphics culling. For example, large objects just outside the cone will be culled when parts of them are actually in the viewing cone.

def inside_view_cone(view_vertex_pos, view_direction, halfview_cosine, object_pos):
  object_vector = object_pos.sub(view_vertex_pos)
  object_cosine = object_vector.dot(view_direction)/(object_vector.length()*view_direction.length())
 
  return object_cosine > halfview_cosine

Polygons

The centroid of a triangle

You might recall from high school geometry that the centroid of a triangle is where the three medians intersect (a median is a point from a vertex to a point halfway between the other two vertices). The centroid has various interesting properties, (it is the centre of gravity of a triangle, for instance). The property that interests us, though, is rather mundane: it is a point that is guaranteed to lie inside the triangle that it is easy (and fast) to compute. Finding a point inside a triangle is useful, for example, in determining whether some other point is inside or outside a convex polygon – see below.

If a triangle has vertices \(\bv u\), \(\bv v\) and \(\bv w\), the centre of gravity is given by their average, \((\bv u + \bv v + \bv w)/3\).

The implementation is simply:

def triangle_centroid(triangle):
  centroid = triangle.vertices[0].add(triangle.vertices[1].add(triangle.vertices[2])).scale(1.0 / 3.0)
 
  return centroid

Is this polygon convex?

A convex polygon is one that has the following property: a line segment between any two points inside the polygon lies entirely inside the polygon. A polygon that is not convex is called concave. Algorithms designed for convex polygons are much simpler than algorithms designed to work on all polygons.

Testing whether a general polygon is convex is tricky. There are several tests that cover most cases; the one below is such a test. It does not work in cases where the polygon is self-intersecting (such as a pentagram), where three consecutive vertices lie in a straight (or very nearly straight) line, or when consecutive vertices coincide. However, if we eliminate these cases, the test is still useful. For a deeper treatment of this subject, and a robust test that works in all cases, see the book Real-Time Collision Detection.

The test below is based on the following property for convex polygons: if we walk along the edges, we will only turn left at each vertex to the next, or we will only turn right. If we do both, the polygon is not convex.

The perp dot product can be used to determine whether we are turning left or right. Thus, we traverse consecutive edge pairs, compute the perp dot product between each pair (using vectors from one vertex to the next for each edge), and make sure that the sign is the same for all of them.

def is_polygon_convex(polygon, point):
  for vertex_index in range(polygon.vertex_count):
    edge_line = polygon.get_edge_line(vertex_index)
 
    if vertex_index == polygon.vertex_count:
      next_vertex_index = 0
    else:
      next_vertex_index = vertex_index + 1
 
    next_edge_line = polygon.get_edge_line(next_vertex_index)
 
    if vertex_index == 0:
      direction = signum(edge_line.perpdot(next_edge_line))
    else:
      if signum(edge_line.perpdot(next_edge_line)) != direction_sign:
        return False
 
  return True

Is this point inside this convex polygon?

To test whether a point is inside a convex polygon, do the following: Choose any three vertices of the polygon (the first three, for convenience), and calculate the centroid. This point is guaranteed to lie inside a convex polygon.

Now for each of the edges, determine whether the given point lies on the same side of the edge as the calculated interior point. If this test passes for all edges, the point is inside.

Here is how this looks:

def point_inside_convex_polygon(polygon, point):
  interior_point = triangle_centroid(polygon)
 
  for vertex_index in range(polygon.vertex_count):
    edge_line = polygon.get_edge_line(vertex_index)
 
    if not is_point_on_same_side_of_line(edge_line, interior_point, point):
      return False
 
  return True

Intersection

To find points of intersection, we usually have to solve equations. When working with vectors, matrices provide an elegant way of doing this (but that is another article!) Fortunately, efficient tests for intersection can be performed without explicitly finding the points of intersection, as the next few examples show.

You already saw some intersection tests in the Collision Detection series. The tests below are essentially equivalent, provided here so that you can see how they look when implemented with vectors.

Do these two lines intersect?

Lines intersect when they are not parallel.

def lines_intersect(line1, line2):
   not lines_parallel(line2, line2)

Note that this is exactly the same test given in the Collision Detection series, since the direction vectors of the two lines are \([x_{A_2} – x_{A_1}, y_{A_2} – y_{A_1}]\) and \([x_{B_2} – x_{B_1}, y_{B_2} – y_{B_1}]\). Their perp dot product is the value

$$(x_{A_2} – x_{A_1})(y_{B_2} – y_{B_1}) – (y_{A_2} – y_{A_1})(x_{B_2} – x_{B_1}).$$

Which is exactly what the original function checks for.

Do these two line segments intersect?

Here we check whether the line of segment1 lies between the endpoints of segment2, and then the other way around.

def line_segments_intersects(segment1, segment2):
  if not is_line_between_points(segment1.get_edge_line(0), segment2.vertices[0], segment2.vertices[1]):
    return False  
 
  return is_line_between_points(segment2.get_edge_line(0), segment1.vertices[0], segment1.vertices[1])

Does this line and circle intersect?

This algorithm is based on the following observation: if a line intersects a circle, the length of the perpendicular from the circle centre to the line will be less than the radius.

Let w = the vector from any point on the line to the circle centre, and let v be the direction vector of the line. Then

$$|\bperp _{\bv v} \bv w| < r$$

Taking squares,

$$|\bperp _{\bv v} \bv w|^2 < r^2$$

Substituting the formula for perp:

$$\left|\frac{\bv v * \bv w}{|\bv v|^2} \bv v^\perp\right|^2 < r^2$$

Simplifying,

$$(\bv v * \bv w)^2/|\bv v|^2 < r^2$$

And finally, avoiding division, we get:

$$(\bv v * \bv w)^2 < r^2|\bv v|^2$$

This inequality is what we use for our check.

def line_intersects_circle(line, circle):
  line_to_circle = circle.center.sub(line.point)  
 
  return sqr(line.direction.perpdot(line_to_circle)) <= line.direction.sqr_length() * sqr(circle.radius)

Do these circles intersect?

Two circles through points \(\bv p\) and \(\bv q\), with radius \(r\) and \(s\) intersect if \(|\bv p – \bv q|^2 > (r + s)^2\).

def circles_intersect(circle1, circle2):
  circle_to_circle = circle2.center.sub(circle1.point)
  return circle_to_circle.square_length() <= sqr(circle1.radius + circle2.radius)

Tips for doing geometric calculations

Watch out for special cases

In all of the recipes above, I have conveniently ignored special cases, such as points “on” lines (instead of either left or right). Sometimes this is also a valid approach when programming, provided that you consider the implications carefully. I usually structure algorithms so that the special case is bundled with the “correct” action, so that I need no explicit special case testing. For example, in a piece of code where special action is taken when the camera is upside down, the special case of the camera being exactly horizontal is bundled with the case where the camera is pointing up. My test will then be

if camera_direction.z < 0:
  turn_camera_up()

Special cases might occur so infrequently that coding special actions (even when it makes sense) is overkill.

One type of special case cannot be ignored in practice, however unlikely it is to occur: special cases that lead to illegal computations (such as division by 0). Especially watch out for normalising vectors: trying to normalise a 0 vector is a nasty thing. In fact, 0 vectors themselves are nasty things. Try to get in the mindset where you question how any vector in your program might be 0, what the implications of such a zero vector are, and how it can be avoided if necessary.

The following piece of code shows a safe version of checking whether a line is between two points.

def is_point_on_same_side_of_line(point, line, reference_vector):
  if is_parallel(line.dir, reference_vector):
    throw Exception('Line and reference vector is parallel')  
 
  v = sub(point, line.point)
 
  if v.is_zero():
    throw Exception('Point is on the line')
 
  return dot(v, reference_vector) > 0

The last check can be eliminated if we change our return value to 1 if the point and reference vector lie on the same side of the line, -1 if they lie on opposite sides, and 0 if the point lies on the line.

def is_point_on_same_side_of_line(point, line, reference_vector):
  if is_parallel(line.dir, reference_vector):
    throw Exception('Line and reference vector is parallel')  
 
  v = point.sub(line.point)
 
  return signum(v.dot(reference_vector))

Object Construction

def bisect_angle(v1, v2):
  return v1.normalized_dupliacte().add(v2.normalized_dusplicate())

Construction of new objects is slow in most, if not all, languages. The piece above might be called hundreds or thousands of times – and it constructs three new vectors. Using in-place operations can reduce object construction significantly.

Here is an example of in-place addition, followed by an improved bisecting algorithm.

class Vector2D:
...
  def add_inplace(self, other):
    self.x += other.x
    self.y += other.y

def bisect_angle(v1, v2):
  length_v1 = v1.length()
  length_v2 = v2.length()
 
  new_vector = v1.duplicate()
  new_vector.add_in_place(v2)
  new_vector.scale(length_v1*length_v2)
 
  return new_vector

Robustness and Precision

Many of the computations presented here will have unexpected behaviour because of round-off errors. Implementing algorithms that are numerically stable (that is, behaves as expected despite round-off) is an advanced topic – if you are interested in some details, look at the resources on this page. It is also interesting to read the philosophy behind CGAL (a Geometry Library in C++).

In naive implementations, a lot of problems can be avoided by handling equality appropriately. Consider for example the parallel checking algorithm:

def is_parallel(line1, line2):
  v = line1.dir
  w = line2.dir
  return perpdot(v, w) == 0

The last line check performs an equality check. But there might be a slight error in either the lines’ direction vectors, or the computation of the perp dot product, which means the equality check is just about useless.

As is common with floating point numbers, we can use interval checking instead:

def is_parallel(line1, line2):
  v = line1.dir
  w = line2.dir
  return abs(perpdot(v, w)) < 0.001

This is already an improvement, but it suffers from an important defect: the longer the direction vectors are, the less likely they are deemed parallel! To see what is going on, recall that the perp dot product gives the sine of the angle between them, multiplied with the lengths of the vectors. Now if they are exactly parallel, this does not matter, because the sine between them will be exactly 0. But the lengths do affect the answer when they are only almost parallel.

To eliminate this defect, we have to normalise the vectors first:

def is_parallel(line1, line2):
  v = line1.dir.normalize_duplicate()
  w = line2.dir.normalize_duplicate()
  return abs(perpdot(v, w)) < 0.001

The above is more correct, but the expense of normalising might not be justifiable. This illustrates why storing line directions as normalised vectors is a good idea.

Another important example is to check whether two vectors are equal. Here is the unsafe check (in 2D):

def vector_equals(v1, v2)
  return v1.x == v2.x and v1.y == v2.y

We can improve this check by using intervals, like this:

def vector_equals(v1, v2)

return abs(v1.x – v2.x) < 0.001 and abs(v1.y – v2.y) < 0.001

But this still has a problem. Consider these vectors:

\(\bv u = [0, 0]\), \(\bv v = [0.0009, 0.0009]\), \(\bv w = [0.001, 0]\)

The algorithm will return true if we test \(\bv u\) and \(\bv v\), but false when we test \(\bv u\) and \(\bv w\), despite the fact that \(\bv w\)‘s endpoint is actually closer to \(\bv u\). To solve this problem, we have to check how far the endpoints of the vectors are apart. Thus, our algorithm becomes:

def vector_equals(v1, v2)
  return v1.sub(v2).sqr_length() < 0.000001

We use a much smaller threshold value here (0.001²) because we are dealing with square lengths.

Build your geometric intuition

Although algebraic fluency is necessary when working with vectors, a strong geometric intuition can help find solutions that are not algebraically clear. A good intuition also helps with identifying all the cases that an algorithm should be able to handle.

CaR is a program that is useful for exploring geometric features of a problem alla Euclid. Every construction starts with some free points. Lines are drawn through points, parallel or perpendicular to other lines. Circles and angles are similarly constructed relative to existing points and lines. When you move the free points, the whole construction moves with, keeping original relationships (intersection, parallelism, etc.) intact.

You will be amazed that after solving problems with CaR for a few days you will be able to do geometry in your head. (If you are looking for good problems, check out the resources on this page. The Notes on Euclidean Geometry (pdf) is a gem.)