Colinearity tests in Blender meshes using Numpy

I re-implemented the algorithm used in the select_colinear_edges add-on to select all edges that are co-linear with already selected edges, and I thought a little write-up with some details could be useful for some people.

Warning! Long read! (≈ 20 min)

The challenge

If we want to select a path of co-linear edges all we have to do is start from any already selected edge, check if its neighbor is co-linear and if it is, select it and proceed from there. If we are careful not to examine any edges more than once, this algorithm will be quite fast and the time will depend on the number of directly connected edges that prove to be co-linear. And even in a large mesh this is likely to be a small number.

But what if we do not require those connected edges to form an unbroken path?

Then for all initially selected edges we would have to test all other edges in the mesh for co-linearity, something that can take a very long time if the mesh contains millions of vertices and everything is implemented in Python using just the mathutils module.

The algorithm

How do you determine if two edges are co-linear?

The first step is to see if they are parallel. This is done by calculating the dot product of the two normalized direction functions. If this product is very close to 1 or -1 we consider them parallel.

(The dot product of two normalized vectors is the cosine of the angle between them)

Being parallel is a necessary condition but not a sufficient one to determine if two edges are co-linear. We also need to check if they are on the same line. This is done by first calculating the vector from anyone of the two vertices in one edge to any one of the vertices in the other edge.

E3 is parallel to E1 but the light blue between vector is not parallel to E1

If the length of this vector is zero, then the chosen vertices coincide, and the edges are co-linear. If not we check the angle between this vector and the direction vector of one of the edges, and if this is is very close to 1 or -1, the edges are co-linear.

This means that for all edges we need to calculate the normalized direction vector and for all initially selected edges we need to calculate this between vector for all other edges.

A numpy based solution

Numpy can work efficiently on vast arrays of numbers and is bundled with Blender. By using Numpy we can avoid two notoriously slow things: Python loops and calling functions.

Our function looks like this (See the function colinear_edges() in this file):

def colinear_edges(selected: np.ndarray, indices, coords, threshold):
    colinear = np.zeros_like(selected)

    # calculate direction vectors for each edge
    edge_dirs = coords[indices[:, 1]] - coords[indices[:, 0]]
    edge_dirs = edge_dirs / np.linalg.norm(edge_dirs, axis=1)[:, np.newaxis]

    for e in selected.nonzero()[0]:
        # get the direction vector of the selected edge
        dir1 = edge_dirs[e]
        # check all other edges for colinearity
        angles = np.arccos(np.clip(np.dot(dir1, edge_dirs.T), -1.0, 1.0))
        parallel = (angles < threshold) | (np.abs(angles - np.pi) < threshold)
        v1 = coords[indices[e, 0]]
        w1 = coords[indices[:, 0]]
        # vector between start points
        between = w1 - v1
        # if the vector between start points is zero, they share a vertex, so colinear
        between_length = np.linalg.norm(between, axis=1)
        connected = between_length < 1e-6
        angles_between = np.abs(
            np.arccos(
                np.clip(
                    np.dot(dir1, (between / between_length[:, np.newaxis]).T), -1.0, 1.0
                )
            )
        )
        bparallel = (angles_between < threshold) | (
            np.abs(angles_between - np.pi) < threshold
        )
        # colinear if they are parallel and either share a vertex or the angle between the direction vector and the vector between start points is less than the threshold
        colinear |= (connected | bparallel) & parallel

    return colinear

Lets explain a few important steps.

The function is called with 4 arguments, a boolean array that indicates which edges are currently selected, an array with indices (2 for each edge) that indexes the third argument an array with vertex coordinates, and a threshold value we'll discuss later. All those arrays come from a Blender Mesh object and we will see how later in this article.

Line 5+6: Here we calculate all direction vectors between the edge indices in one go, and then normalize them in a single statement by dividing each vector by its norm (i.e. length).

Line 8-10: We loop over each selected edge and get its direction vector.

Line 12: Then we calculate the angles with all other vectors. This is done by calculating the dot product between the direction vector and all other direction vectors in on go (note that we need to transpose the array of vectors for this to work). We clip the dot products between -1 and 1 to guard against any floating point inaccuracies and then use the arccos() function to calculate the angle (Remember that the dot product represents the cosine of the angle between two vectors)

Line 13: then the angle is checked against the threshold and if smaller (or very close to π, because we don´t care in which direction the vectors are aligned) we deem them parallel.

Line 14-17: then we take a vertex v1 from the first edge and all vertices w1 for each other edge, and calculate the between vector.

Line 19+20: we calculate the length all those between vectors, and for all of them determine if this length is so small we consider the vertices coincident.

Line 21-27: then we calculate all angle between the direction vector and the between vectors in the same way we did before.

Line 28-30: we then determine if the between vectors are parallel with the direction vector (or anti-parallel, because we don´t care about that)

Line 32: Finally we combine the logic and say two edges are co-linear if the are parallel AND their chosen vertices are coincident OR the angle is near zero. The result is OR-ed into the colinear array because we do this for each edge that was initially selected and want to return the combined set.

Calling colinear_edges()

If we have a Blender Mesh object we can access obj.data.edges and obj.data.vertices. The select-colinear() function takes reference to those properties and uses then to efficiently retrieving all the indices, selected status and vertex coordinates with the foreach_get() method (Line 7-9). It stores them in arrays we have created first (Line 4-6).

foreach_get() expects flat arrays, so we reshape them into their expected shape where needed (Line 10+11), before we call the colinear_edges() function discussed earlier (Line 13).

The result is a flat array with the new selected status of each edge which we store in the select attribute of the mesh edges with the foreach_set() method (Line 14).

And finally we return the number of selected edges by counting all non-zero values (True is considered non-zero too, so this works fine for arrays of booleans).

def select_colinear(edges, vertices, threshold):
    n_edges = len(edges)
    n_vertices = len(vertices)
    indices = np.empty(2 * n_edges, dtype=int)
    coords = np.empty(3 * n_vertices, dtype=float)
    selected = np.zeros(n_edges, dtype=bool)
    edges.foreach_get("vertices", indices)
    edges.foreach_get("select", selected)
    vertices.foreach_get("co", coords)
    coords = coords.reshape((n_vertices, 3))
    indices = indices.reshape((n_edges, 2))

    colinear = colinear_edges(selected, indices, coords, threshold)
    edges.foreach_set("select", colinear)
    return np.count_nonzero(colinear)

Summary

Using the foreach_get() / foreach_set() lets us easily get access to mesh properties in bulk, which allows us to use Numpy to implement an algorithm that calculates co-linearity without Python loops (except for the loop of all initially selected edges).

In exchange for a modest increase in complexity we gain a lot of performance: Although your mileage may vary of course, I could easily (in < 0.1 second) select all co-linear edges when picking one edge in a default cube that was subdivided 500 times (= around 3.5 million edges). Fast enough for me 😀