创建转换矩阵的numpy数组的正确方法是什么答案

【问题标题】：What is the proper way to create a numpy array of transformation matrices创建转换矩阵的numpy数组的正确方法是什么
【发布时间】：2015-12-10 23:45:35
【问题描述】：

给定一个旋转角度列表（让我们说一下 X 轴）：

import numpy as np
x_axis_rotations = np.radians([0,10,32,44,165])

我可以通过这样做创建一个匹配这些角度的矩阵数组：

matrices = []
for angle in x_axis_rotations:
    matrices.append(np.asarray([[1 , 0 , 0],[0, np.cos(angle), -np.sin(angle)], [0, np.sin(angle), np.cos(angle)]]))
matrices = np.array(matrices)

这可行，但它没有利用 numpy 处理大型数组的优势......所以如果我的角度数组有数百万，这样做不会很快。

有没有更好（更快）的方法从输入数组创建变换矩阵数组？

【问题讨论】：

如果我理解正确的话，你想要一个 3D 数组吗？输入来自哪里？文件或您在运行时生成的东西？
您能解释一下为什么要存储矩阵而不是在使用时计算它们的原因吗？
@Simpom 我想实现this paper 中概述的点到三角形技术。这个想法是，通过为网格的每个三角形预先计算一个变换矩阵，可以将点到三角形的 3D 距离问题减少为更简单的 2D 问题，这可以比通过重心计算更快地计算。如果我想要的只是单点查询，那么是的，不需要计算矩阵。但是，如果我有一百万个点要查询，那么这应该很有用。
@ventsyv 我想创建一个 4x4 矩阵的数组。阅读我刚刚写给 Simpom 的评论，了解我想要的特定用例。输入来自外部程序，形式为描述三角形网格的顶点位置数组，以及用于查询网格上最近点的 3d 点数组。

标签： python arrays performance numpy matrix

【解决方案1】：

这是一个直接而简单的方法：

c = np.cos(x_axis_rotations)
s = np.sin(x_axis_rotations)
matrices = np.zeros((len(x_axis_rotations), 3, 3))
matrices[:, 0, 0] =  1
matrices[:, 1, 1] =  c
matrices[:, 1, 2] = -s
matrices[:, 2, 1] =  s
matrices[:, 2, 2] =  c

时间，对于好奇的人：

In [30]: angles = 2 * np.pi * np.random.rand(1000)

In [31]: timeit OP(angles)
100 loops, best of 3: 5.46 ms per loop

In [32]: timeit askewchan(angles)
10000 loops, best of 3: 39.6 µs per loop

In [33]: timeit divakar(angles)
10000 loops, best of 3: 93.8 µs per loop

In [34]: timeit divakar_oneline(angles)
10000 loops, best of 3: 56.1 µs per loop

In [35]: timeit divakar_combine(angles)
10000 loops, best of 3: 43.9 µs per loop

所有都比你的循环快得多，所以使用你最喜欢的一个:)

【讨论】：

简洁明了！谢谢！

【解决方案2】：

你可以使用linear indexing来帮忙，像这样-

# Get cosine and sine values in one-go
cosv = np.cos(x_axis_rotations)
sinv = np.sin(x_axis_rotations)

# Get size parameter
N = x_axis_rotations.size

# Initialize output array
out = np.zeros((N,3,3))

# Set the first element in each 3D slice as 1
out[:,0,0] = 1

# Calculate the first of positions where cosine valued elements are to be put
idx1 = 4 + 9*np.arange(N)[:,None]

# One by one put those 4 values in 2x2 blocks across all 3D slices
out.ravel()[idx1] = cosv
out.ravel()[idx1+1] = -sinv

out.ravel()[idx1+3] = sinv
out.ravel()[idx1+4] = cosv

或者，您可以在使用zeros 初始化输出数组并将每个切片中的第一个元素设置为1 后一次性设置元素，就像这样 -

out.reshape(N,-1)[:,[4,5,7,8]] = np.column_stack((cosv,-sinv,sinv,cosv))

在上面提到的两种方法之间，可能会演变出另外两种middleground方法，在用零初始化并将每个 3D 切片中的第一个元素设置为1 之后再次放置，就像这样 -

out.reshape(N,-1)[:,[4,8]] = cosv[:,None]
out.reshape(N,-1)[:,[5,7]] = np.column_stack((-sinv[:,None],sinv[:,None]))

最后一个是 -

out.reshape(N,-1)[:,[4,8]] = cosv[:,None]
out.reshape(N,-1)[:,5] = -sinv
out.reshape(N,-1)[:,7] = sinv

【讨论】：

我喜欢这些优雅的答案，但我认为索引最终只会减慢速度。
@askewchan 好吧，我很害怕 column_stacks，因为众所周知它们很慢。我想最好把所有这些都计时。
你最后的解决方案不应该有最后两行的新轴，因为 LHS 上的切片只有一个维度。
@askewchan 谢谢，确实是这样，已编辑。我得到了和你一样的运行时结果。
@Divakar 谢谢！我会纯粹因为它是多么简洁而接受 askewchan 的回答，但你的回答紧随其后！