geometry.py

def degree_rep(angle):
    from math import degrees
    return '∠ %.2f°' % degrees(angle)

def position_vector(vertex):
    assert isinstance(vertex, Vertex)
    return Vector(*vertex.coordinates)

def angle_magnitude_key(angle):
    '''
        Returns how great an angle is.
    '''
    from math import pi as π
    normalized_angle = abs(angle) % (2 * π)
    if normalized_angle > π:
        normalized_angle = (2 * π) - normalized_angle
    return normalized_angle

class Vertex:
    def __init__(self, x, y, z):
        if not all(is_real(coordinate) for coordinate in (x, y, z)):
            raise ValueError(str.format('Bad vertex coordinates: {!r}',
                (x, y, z),
            ))
        self.x, self.y, self.z = x, y, z
    def __repr__(self):
        return str.format('Vertex({!r}, {!r}, {!r})', self.x, self.y, self.z)
    def distance_to(self, other):
        # can't use hypot because of 3 arguments
        from math import sqrt
        return sqrt(
            (self.x - other.x) ** 2 +
            (self.y - other.y) ** 2 +
            (self.z - other.z) ** 2
        )
    @property
    def coordinates(self):
        return self.x, self.y, self.z
    def __add__(self, other):
        assert not (type(self) == type(other) == Vertex)
        if type(self) == Vector and type(other) == Vector:
            result_type = Vector
        else:
            result_type = Vertex
        return result_type(self.x + other.x, self.y + other.y, self.z + other.z)
    def __radd__(self, other):
        return self + other
    def __neg__(self):
        return type(self)(-self.x, -self.y, -self.z)
    def __sub__(self, other):
        result = self + (-position_vector(other))
        if isinstance(other, Vertex):
            return Vector(*result.coordinates)
        else:
            return result
    def __mul__(self, scalar):
        assert is_real(scalar)
        return type(self)(self.x * scalar, self.y * scalar, self.z * scalar)
    def __rmul__(self, other):
        return self * other
    def __truediv__(self, scalar):
        assert is_real(scalar)
        return type(self)(self.x / scalar, self.y / scalar, self.z / scalar)
    def __floordiv__(self, scalar):
        assert is_real(scalar)
        return type(self)(self.x // scalar, self.y // scalar, self.z // scalar)
    def __matmul__(self, transformation_matrix):
        return transform(self, transformation_matrix)
    def __eq__(self, other):
        return self.coordinates == other.coordinates
    def __lt__(self, other):
        return self.coordinates < other.coordinates
    def __hash__(self):
        return hash(self.coordinates)

class VertexOp:
    def __init__(self, callable):
        self.callable = callable
    def __rmul__(self, coordinate):
        return self.callable(coordinate)

class LineSegment:
    def __init__(self, v1, v2):
        self.v1, self.v2 = v1, v2
    def __repr__(self):
        return str.format('LineSegment({!r}, {!r})', self.v1, self.v2)
    @property
    def length(self):
        return self.v1.distance_to(self.v2)
    @property
    def vertices(self):
        return self.v1, self.v2
    def __len__(self):
        return 2
    def __getitem__(self, index):
        return self.vertices[index]

class IndexRing:
    def __init__(self, container):
        self.container = container
    def __getitem__(self, index):
        return self.container[index % len(self.container)]

class Vector(Vertex):
    def __init__(self, x, y, z):
        super().__init__(x, y, z)
    def __repr__(self):
        return str.format('Vector({!r}, {!r}, {!r})', self.x, self.y, self.z)
    def length(self):
        from math import sqrt
        return sqrt(self.x ** 2 + self.y ** 2 + self.z ** 2)
    def normalized(self):
        length = self.length()
        if length > 0:
            return Vector(self.x / length, self.y / length, self.z / length)
        else:
            return Vector(0, 0, 0)

def dot_product(v1, v2):
    ''' Returns the dot product v1 · v2. '''
    return v1.x * v2.x + v1.y * v2.y + v1.z * v2.z

def cross_product(v1, v2):
    ''' Returns the cross product v1 × v2. '''
    return Vector(
        x = v1.y * v2.z - v1.z * v2.y,
        y = v1.z * v2.x - v1.x * v2.z,
        z = v1.x * v2.y - v1.y * v2.x,
    )

def unit_normal(a, b, c):
    x = Matrix3x3([
        1, a.y, a.z,
        1, b.y, b.z,
        1, c.y, c.z,
    ]).determinant()
    y = Matrix3x3([
        a.x, 1, a.z,
        b.x, 1, b.z,
        c.x, 1, c.z,
    ]).determinant()
    z = Matrix3x3([
        a.x, a.y, 1,
        b.x, b.y, 1,
        c.x, c.y, 1,
    ]).determinant()
    return Vector(x, y, z).normalized()

def pairs(iterable, *, count = 2):
    '''
        Iterates the given iterable and returns tuples containing `count`
        sequential values in the iterable.

        Formally, for a vector A with indices 0…n and pair count k, this
        function yields:
            (A₀, A₁, …, Aₖ),
            (A₁, A₂, …, Aₖ₊₁),
            …,
            (Aₙ₋₂, Aₙ₋₁, …, Aₖ₋₂),
            (Aₙ₋₁, A₀, …, Aₖ₋₁).
    '''
    iterable_count = len(iterable)
    for index in range(iterable_count):
        yield tuple(
            iterable[(index + offset) % iterable_count]
            for offset in range(count)
        )

class Polygon:
    def __init__(self, vertices):
        self.vertices = vertices
    def __repr__(self):
        return str.format('Polygon({!r})', self.vertices)
    def area(self):
        total = Vector(0, 0, 0)
        for v1, v2 in pairs(self.vertices):
            total += cross_product(v1, v2)
        return dot_product(total, unit_normal(self.vertices[0], self.vertices[1], self.vertices[2]))
    def centroid(self):
        ...
        raise NotImplementedError
    def __getitem__(self, index):
        return self.vertices[index % len(self.vertices)]
    def __len__(self):
        return len(self.vertices)
    @property
    def perimeter_lines(self):
        for v1, v2 in pairs(self.vertices):
            yield LineSegment(v1, v2)
    @property
    def smallest_angle(self):
        from math import acos, inf
        min_angle = inf
        for v1, v2, v3 in pairs(self.vertices, count = 3):
            vec1 = (position_vector(v3) - position_vector(v2)).normalized()
            vec2 = (position_vector(v1) - position_vector(v2)).normalized()
            cosine = dot_product(vec1, vec2) / vec1.length() / vec2.length()
            min_angle = min(min_angle, angle_magnitude_key(acos(cosine)))
        return min_angle
    @property
    def hairline_ratio(self):
        lengths = [line.length for line in self.perimeter_lines]
        return max(lengths) / min(lengths)

def is_real(number):
    return isinstance(number, int) or isinstance(number, float)

class Matrix3x3:
    '''
        A 3×3 matrix forming the top-left portion of a full 4×4 transformation
        matrix.
    '''
    def __init__(self, values):
        assert(all(is_real(x) for x in values))
        assert len(values) == 9
        self.values = values
    def __repr__(self):
        return str.format('Matrix3x3({!r})', self.values)
    def __getitem__(self, index):
        return self.values[index]
    def determinant(self):
        v = self.values
        return sum([
            +(v[0] * v[4] * v[8]),
            +(v[1] * v[5] * v[6]),
            +(v[2] * v[3] * v[7]),
            -(v[2] * v[4] * v[6]),
            -(v[1] * v[3] * v[8]),
            -(v[0] * v[5] * v[7]),
        ])

def complete_matrix(matrix, anchor):
    '''
        Combines a 3×3 matrix and an anchor vertex into a full 4×4
        transformation matrix.
    '''
    return [
        matrix[0], matrix[1], matrix[2], anchor.x,
        matrix[3], matrix[4], matrix[5], anchor.y,
        matrix[6], matrix[7], matrix[8], anchor.z,
        0, 0, 0, 1,
    ]

def transform(vertex, transformation_matrix):
    '''
        Transforms a vertex by a 4×4 transformation matrix.
    '''
    u = transformation_matrix[0] * vertex.x \
        + transformation_matrix[1] * vertex.y \
        + transformation_matrix[2] * vertex.z \
        + transformation_matrix[3]
    v = transformation_matrix[4] * vertex.x \
        + transformation_matrix[5] * vertex.y \
        + transformation_matrix[6] * vertex.z \
        + transformation_matrix[7]
    w = transformation_matrix[8] * vertex.x \
        + transformation_matrix[9] * vertex.y \
        + transformation_matrix[10] * vertex.z \
        + transformation_matrix[11]
    return Vertex(u, v, w)