Computational Geometry in Python

Posted on August 9th, 2016, by tom in Code

This post provides a list of Python functions that are very common in structural design calculations. They are not specifically optimised in any way. Therefore, you may end up using entirely different implementations in a real project. The list is a work in progress. I will try to add more examples soon…

Before continuing, note that there are several ways to do the same (or similar) things in Python. Some are faster than others, and, depending on who you ask, some are more “pythonic” than others. To illustrate this, but most certainly without making any claims about which one is “better” or faster, I include here 10 different ways to compute the length of a vector. There are probably quite a few more, but some of the variations are already a bit contrived, so I will leave it at that :).

from operator import add, pow 

l = (v[0] ** v[0] + v[1] ** v[1] + v[2] ** v[2]) **0.5
l = (v[0] **2 + v[1] **2 + v[2] **2) **0.5
l = sum([v[i]**2 for i in range(3)]) **0.5
l = sum(v[i]**2 for i in range(3)) **0.5
l = sum(x * x for x in v) **0.5
l = sum(x **2 for x in v) **0.5
l = sum(map(pow, v, [2, 2, 2])) **0.5
l = sum(map(pow, v, [2] * 3)) **0.5
l = reduce(add, map(pow, v, [2, 2, 2])) **0.5
l = reduce(lambda x, y: x + y, [x **2 for x in v]) **0.5

# v = [1., 1., 1.]
# 1.7320508075688772

Note that I have used **0.5 to take the square root. Therefore, for each variation, there is at least one other version using the math.sqrt function instead… For a speed comparison between the two, and many opinions about which one is more “pythonic”, see this StackOverflow post.

And then now, without fruther ado, the list of geometry functions…

Length

def length(v):
    return sqrt(sum(axis * axis for axis in v))

Distance

def distance(a, b):
    return sqrt(sum((a[i] - b[i]) ** 2 for i in range(3)))

Dot product

def dot(u, v):
    return sum(u[i] * v[i] for i in range(3))

Cross product

def cross(u, v):
    return [u[1] * v[2] - u[2] * v[1], 
            u[2] * v[0] - u[0] * v[2],
            u[0] * v[1] - u[1] * v[0]]

Centroid

def centroid(points):
    p = len(points)
    return [sum(axis) / p for axis in zip(*points)]

Note that in the context of vector computations, you could think of the combination of the “unpacking operator (*)” and the zip function as taking the transpose of a list of vectors.

>>> points = [[1, 2, 3], [1, 2, 3], [1, 2, 3]]
>>> zip(*points)
[[1, 1, 1], [2, 2, 2], [3, 3, 3]]

Area (triangle)

ab = [b[i] - a[i] for i in range(3)]
ac = [c[i] - a[i] for i in range(3)]

area = 0.5 * length(cross(ab, ac))

Area (polygon)

area = 0
c = centroid(polygon)

for i in range(-1, len(polygon) - 1):
    a  = polygon[i]
    b  = polygon[i + 1]
    ab = [b[j] - a[j] for j in range(3)] 
    ac = [c[j] - a[j] for j in range(3)] 
    area += 0.5 * length(cross(ab, ac))

Volume

…

Center of mass (polygon)

…

Center of mass (polyhedron)

…

BRG blog