Replace matrix @ vector list comprehensions with something more efficient

Question

I'm looking to replace list comprehensions like the following with something more efficient and I was wondering if numpy can be used.

[obj.matrix_world @ v.co for v in obj.data.vertices]

I can do this, which is significantly faster, but still lacks the matrix multiplication, that needs to happen for each coordinate.

coords = np.empty((len(obj.data.vertices), 3), 'f')
obj.data.vertices.foreach_get('co', np.reshape(coords, len(obj.data.vertices) * 3))

I was hoping it would just be

np.matrix(obj.matrix_world) @ coords

But that fails with (FWIW, It's just a cube, hence (8, 3))

ValueError: shapes (4,4) and (8,3) not aligned: 4 (dim 1) != 8 (dim 0)

Similarly, doing

for co in coords:
    np.matrix(obj.matrix_world) @ co

fails with

ValueError: shapes (4,4) and (3,) not aligned: 4 (dim 1) != 3 (dim 0)

So I've hit a wall now. Numpy can't multiply a 4x4 matrix with a 3d vector, but Blender can. How do I solve this in a numpy way?

---

edit:

Based on @JoseConseco's excellent reply below, this is what I'm doing now

mx = active.matrix_world
verts = active.data.vertices
vert_count = len(verts)

coords = np.empty((vert_count, 3), 'f')
verts.foreach_get('co', np.reshape(coords, vert_count * 3))

coords_4d = np.ones((vert_count, 4), 'f')
coords_4d[:, :-1] = coords

coords = np.einsum('ij,aj->ai', mx, coords_4d)[:, :-1]

Note, that there's no need to convert the (mathutils) Matrix into an np array. It works fine just like that.

Based on my tests, doing this via numpy is ~3 times faster than the list comprehension.

Answer 1

Quickest way AFAIK:

mat = np.array(self.curveObj.matrix_world)

verts_co = np.zeros((vertCount*3), dtype=np.float)
mesh.vertices.foreach_get("co", verts_co)
verts_co.shape = (vertCount, 3)
verts_co_4d = np.ones(shape=(vertCount, 4), dtype=np.float)
verts_co_4d[:, :-1] = verts_co  # cos v (x,y,z,1) - point,   v(x,y,z,0)- vector
splinePoints_world = np.einsum('ij,aj->ai', mat, splinePoints)

Answer 2

In addition to JoseConseco answer, changing

np.einsum('ij,aj->ai', mx, coords_4d)[:, :-1]

into

np.dot(coords_4d, mx.inverted())[:, :-1]

Should improve performance.

The entire script:

vlen = len(obj.data.vertices)

vco = np.empty(vlen * 3, 'f')
obj.data.vertices.foreach_get('co', vco)

coords = np.empty((vlen, 4), 'f')
coords[::4] = 1.0
coords[:,:-1] = vco.reshape((vlen, 3))

coords = np.dot(coords, obj.matrix_world)
res = coords[:,:-1].reshape((vlen*3))

obj.data.vertices.foreach_set('co', res)

reshape is basically free operation in Numpy, so should be used liberally. np.empty also just allocates memory, doesn't write it so maybe saves some time.

Answer 3

Mesh.transform(Matrix)

As an alternative to using numpy

If the desired result is to transform the vertex coordinates by a matrix then the mesh transform method does exactly this, internally with the passed matrix.

ob = context.object
mw = ob.matrix_world
me = ob.data
me.transform(mw) # transforms all verts by matrix

Timing it.

Timing the methods outlined. Have included mesh transform, a method using mesh transform but returning a list of transformed coordinates, and the two numpy methods outlined above.

Note it is my opinion (test code to put empties at result commented out) that using dot product requires the matrix to be transposed, rather than inverted.

import bpy
import numpy as np
from mathutils import Matrix

import time

    
REPS = 1000

def time_it(func):
    def wrapper(*arg, **kw):
        t1 = time.time()
        for i in range(REPS):
            func(*arg, **kw)
        t2 = time.time()
        print(func.__name__, (t2 - t1))
    return wrapper


def vert_coords(me):
    n = len(me.vertices)
    coords = np.empty((n * 3), dtype=float)
    me.vertices.foreach_get("co", coords)
    return np.reshape(coords, (n, 3))


@time_it
def transform00(me, matrix=Matrix()):
    me.transform(matrix)

@time_it
def transform0(me, matrix=Matrix()):
    copy = me.copy()
    copy.transform(matrix)
    coords = vert_coords(copy)
    bpy.data.meshes.remove(copy)
    return coords

@time_it        
def transform1(coords, matrix=Matrix()):
    n = len(coords)
    coords4d = np.empty(shape=(n, 4), dtype=float)
    coords4d[::-1] = 1
    coords4d[:, :-1] = coords
    return np.dot(coords4d, matrix.transposed())[:,:-1]

@time_it    
def transform2(coords, matrix=Matrix()):
    n = len(coords)
    coords4d = np.empty(shape=(n, 4), dtype=float)
    coords4d[::-1] = 1
    coords4d[:, :-1] = coords
    return np.einsum('ij,aj->ai', matrix,  coords4d)[:,:-1]

   
ob = bpy.context.object
mw = ob.matrix_world
mesh = ob.data
copy = ob.data.copy()
coords = vert_coords(mesh)
transform00(copy, mw)
bpy.data.meshes.remove(copy) # clean up
transform0(mesh, mw)
transform1(coords, mw)
transform2(coords, mw)

'''
# add empty at coordinates to test
transform = transform1 # choose a method
for p in transform(vert_coords(mesh), ob.matrix_world):
    bpy.ops.object.empty_add(location=p)
    
'''

Some test runs on default primitives

Cube verts: 8
transform00 0.004887104034423828
transform0 0.18675923347473145
transform1 0.08976554870605469
transform2 0.0721430778503418
Sphere verts: 482
transform00 0.014344215393066406
transform0 0.40707826614379883
transform1 0.12292027473449707
transform2 0.12022829055786133
Suzanne verts: 507
transform00 0.012979984283447266
transform0 0.398883581161499
transform1 0.1199808120727539
transform2 0.11667346954345703

Findings.

If the outcome is transforming the vertices then Mesh.transform is way quicker.

The einsum method comes out a squeak quicker than dot product

Trying to leverage off transform, to get an array of transformed coords by copying a mesh, transforming and foreach getting the result, has to much overhead and is slowest.

Iterating over each vertex via a list comprehension (results not shown) is slower again.

Answer 4

This is not an answer, but a consecutive question. I hope that is ok.

I'm working on a modal operator and want to quickly get vertex coordinates to pass them to a shader.

The code works fine when working with a list comprehension like this:

self.vert_data = np.array([mat @ v.co for v in verts], np.float32)

Here is how it looks like in the viewport, when working as intended:

For some reason, I'm not able to get implement any of the above methods. Whichever I use, the result turns out like this:

The vertices blow up and cover the entire scene (or so it seems). I'm pretty sure that I'm missing something very obvious here, but I can't find. Here is the code that I wrote following the answers above:

coords = np.zeros(vert_count * 3, dtype=np.float32)
verts.foreach_get('co', coords.ravel())
coords = coords.reshape(vert_count, 3)
coords_4d = np.ones((vert_count, 4), 'f')
coords_4d[:, :-1] = coords
# self.vert_data = np.dot(coords_4d, mat.inverted())[:, :-1]
self.vert_data = np.einsum('ij,aj->ai', mat, coords_4d)[:, :-1]

Many thanks!