# Replace matrix @ vector list comprehensions with something more efficient

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.

• Even though your problem takes place in Blender, it's a generic math problem and there is a StackExchange for math. I imagine that you will have an easier time finding linear algebra guru's over yonder ;) There will be some here too of course, no offense intended :) math.stackexchange.com May 1 '19 at 10:07
• Maybe yes, but I don't speak their language :) and felt this should be a common enough (or at east useful) issue in Blender circles. May 1 '19 at 11:44
• I'm voting to close this question as off-topic because it is a general Python problem which isn't specific to Blender. May 1 '19 at 12:03
• I think the question is fine. I would suggest editing title, such that it is more akin to the question "replace matrix @ vector list comprehensions with something more efficient" May 1 '19 at 12:23

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)


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.

• Thank ambi, I appreciate you chiming in! I'm also using empty and reshape. I've also just noticed that float64 seems to be faster than float32. I still need float32(need coords for gpu drawing) at the end, so I just change it at the end, and it's still faster. As for, np.dot() I'm actually getting a different result from it than via the einsum() method. I need to study this some more. Thanks again, much appreciated! May 2 '19 at 9:13

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):

'''


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.

• Thanks, that's really interesting to use mesh.transform(). It is indeed incredibly fast and I've started to use it in some other areas now. May 2 '19 at 8:25