Implementation of SAT (Separating axis theorem)

A project I was working on required the usage of the Separating Axis Theorem to detect collisions between two convex polygons in real time. So I implemented a basic class (ConvexFrame in the code) to keep track of the information about the polygon that is necessary for the SAT algorithm. I implemented it, and it works, but it is unfortunately quite slow for the project I'm working on (which is a game by the way).

For example, testing for a collision between an 8-sided shape and a 3-sided shape takes on average 0.00015 seconds of computation time, which means that in a normal game tick (1/60s = 0.016...s) I can calculate a maximum of 100 collisions between convex polygons. I need it to be faster than this, and I'm not sure how I can optimize it. Can someone help me understand where I can optimize the code?

The code is split up into 2 main files: geometry.py (imported as geo in the other file), and physical.py. geometry.py contains basic functions for vector calculations and the likes, where physical.py contains the SAT algorithm and the ConvexFrame class. I made sure that most of the functions in the geometry file were as optimized as I could get them to be, so that shouldn't be the problem, but just incase I included the average runtime of each of the functions in geo.

geometry.py:

```
import math
import maths # maths is an even simpler file containing constants and other basic functions
# there is no need to include it here.

def centroid(*points):
"""Calculate the centroid from a set of points."""
# Average time for 4 points: 1.4572602962591971e-06s
x, y = zip(*points)
_len = len(x)
return [sum(x)/_len, sum(y)/_len]

def icentroid(*points):
"""Faster than normal centroid, but returns an iterator.

Since this returns an iterator, to separate it up into an
(x, y) tuple, simply say:

>>> x, y = icentroid(*points)
"""
# Average time for 4 points: 9.622882809023352e-07s
_len = len(points)
retu

1. Vectors

Much of this code is awkward and long-winded because points and vectors are represented by plain Python tuples. A simple operation like subtracting two points requires disassembling the points into their elements, subtracting the elements, and then reassembling the result. If the code represented points using some kind of vector data structure, then a lot of it could be simplified.

For example, here are eight lines of code for computing the offset of each vertex from the origin:

self._offsets = []
orx, ory = self._origin
append_to_offsets = self._offsets.append
for vertex in coordinates:
    x, y = vertex
    offx = x - orx
    offy = y - ory
    append_to_offsets([offx, offy])

but with a class of vectors that supported subtraction, this would be one line:

self._offsets = [v - self._origin for v in coordinates]

If you have NumPy to hand, then it would make sense to use NumPy arrays as your vectors, but if not then it's easy to write such a class. For example, you might start with something simple like this:

class Vector(tuple):
    def __add__(self, other):
        return Vector(v + w for v, w in zip(self, other))

    def __sub__(self, other):
        return Vector(v - w for v, w in zip(self, other))

    def __mul__(self, s):
        return Vector(v * s for v in self)

    def __abs__(self):
        return sqrt(sum(v * v for v in self))

    def dot(self, other):
        """Return the dot product with the other vector."""
        return sum(v * w for v, w in zip(self, other))

(See my vector.py for a full-featured implementation.)

With this class, many of your geometry functions could be simplified. For example, to calculate the distance from v1 to v2 you currently have:

x1, y1 = v1
x2, y2 = v2
deltax = x2 - x1
deltay = y2 - y1
return math.sqrt(deltax * deltax + deltay * deltay)

but using the Vector class given above, this becomes so trivial that it might not be worth defining a function for it:

return abs(v1 - v2)

This approach won't speed up your code (the same operations are being carried out) but it will make it shorter, clearer, and easier to work with, and that will help you when you do come to make performance improvements.

1. Projection

The algorithm in project has a bug: there's a division by zero error if start and end have the same y-coordinate:

>>> project((1, 2), (0, 0), (2, 0)) # expecting (1, 0)
Traceback (most recent call last):
    m2 = -deltax/deltay
ZeroDivisionError: division by zero

In computer geometry you should always use vectors if possible: trying to work with the slope-and-intercept representation of lines leads into difficulty because of the exceptional cases.

The reliable way to project v onto the line from start to end is to use the dot product, like this:

w = end - start
return w * (w.dot(v) / w.dot(w))

This still has the possibility of division by zero, but only in the case where w is zero (that is, if start == end) and in that case no projection is possible.

1. Performance

-
Pure Python is always going to be a bit slow for this kind of problem, so you should look into using NumPy to speed up the calculations in collide. In this function you carry out a computation for each edge of each figure: this would be easy to vectorize so that all computations are carried out at once.

-
Because the collision test is so expensive, it's worth taking some trouble to avoid it. In particular, it would be worth storing a bounding circle with origin \$ o \$ and radius \$ r \$ for each polygon: namely, the smallest circle that contains all points in the polygon. Then, before doing the full polygon/polygon collide test, do a circle/circle test: if two polygons have bounding circles \$ (o_1, r_1) \$ and \$ (o_2, r_2) \$ then they can only collide if \$ \left| o_1 - o_ 2 \right| ≤ r_1 + r_2 \$. This should allow you to reject most collisions cheaply.

-
Consider storing your polygons in a space-partitioning data structure such as a quadtree that would let you efficiently find pairs of polygons that might collide. SciPy has an implementation in scipy.spatial.KDTree.

-
It's likely that you are going to be repeatedly testing the same set of polygons for collision (for example, in a video game you'd be doing this each frame). In that case, when you find that two polygons are separated by a particular axis, remember that axis and test it first next time. The insight is that a pair of polygons don't move very far in one time step, and so an axis that separates them at time \$ t \$ will continue to separate them at time \$ t + δ \$. This is the method of "caching witnesses" — see Rabbitz, "Fast Collision Detection of Moving Convex Polyhedra" in Graphics Gems IV.

StackExchange Code Review Q#47111, answer score: 6