On Fri, Nov 21, 2008 at 2:57 AM, Sébastien Barthélemy barthelemy@crans.orgwrote:

Hello,

I would like to port a matlab library which provides functions for rigid body mechanics such as operations on homogeneous matrices (in SE(3)), twists (in se(3)) and so.

In matlab, the library worked on 3d matrices: n homoneous matrices were stacked along the 3d dimension. This speeded up computations on multiple matrices, at the price that native matrix opertators such H1

• H2 did not work (it is not defined for 3d matrices).

In this spirit, in numpy a set of rotation matrices could be built in the following way:

def rotx(theta): """ SE(3) matrices corresponding to a rotation around x-axis. Theta is a 1-d array """ costheta = np.cos(theta) sintheta = np.sin(theta) H = np.zeros((theta.size,4,4)) H[:,0,0] = 1 H[:,3,3] = 1 H[:,1,1] = costheta H[:,2,2] = costheta H[:,2,1] = sintheta H[:,1,2] = sintheta return H

I'm now seeking advices regarding an implementation with numpy (it's my first time with numpy).

• Is there any difference between a 3d array and a 1-d array of 2-d

arrays (seems not)

• would you use a 3d-array or a list of 2d arrays or anything else to

store these matrices ?

• Is there a way to work easily on multiple matrices of the same shape

? (in single-instruction/multiple-data spirit) for instance if A.shape==(3,2,4) and B.shape=(3,4,1), what is the best way to compute C for i = (0,1,2) C[i,:,:] = np.dot(A[i,:,:],B[i,:,:])

There is work going on along these lines called generalized ufuncs, but it isn't there yet. The easiest current workaround is to use lists of matrices or arrays. Or you can get clever and do the matrix multiplications yourself using the current numpy operators, i.e.

C = (A[:,:,:,newaxis]*B[:,newaxis,:,:]).sum(axis=-2)

But this uses more memory than the list of matrices approach.

Chuck