我正在使用OpenCV fitLine将线拟合到3D点。计算所得拟合残差的最佳方法是什么？或者，由于除了拟合之外我还想要残差，是否有比fitLine更好的方法？
可以使用以下方法，但是必须有更好的方法。
# fit points
u, v, w, x, y, z = cv2.fitLine(points, cv2.DIST_L2, 0, 1, 0.01)
v = np.array(u[0], v[0], w[0])
p = np.array(x[0], y[0], z[0])
# rotate fit to z axis
k = np.cross(v, [0, 0, 1])
mag = np.linalg.norm(k)
R, _ = cv2.Rodrigues(k * np.arcsin(mag) / mag)
# rotate points and calculate distance to z-axis
rot_points = np.dot(R, (points-p).T)
err = rot_points[0]**2 + rot_points[1]**2
我假设
fitLine在估算行的同时计算残差err，因此不得不自己重新计算它们是一种浪费。基本上，知道我要线和残差后，是否有比fitLine更好的选择，后者仅返回线？ 最佳答案
我不知道有任何直接方法可以从cv2.fitLine本身获取残差的总和，而我将只专注于加速现有代码。现在，在以相对较高的点进行基准测试后，它表明大部分运行时都花在了最后两行，即得到rot_points和err。另外，似乎我们并没有完全使用rot_points的最后一行来计算err，因此希望我们可以通过仅切成前两行来节省一部分运行时。
让我们深入研究获得rot_points和err的有效方法。
1)rot_points = np.dot(R，(points-p).T)
此步骤涉及以points-p广播，此后通过矩阵乘法来减少。现在，broadcasting占用大量内存，如果我们分别将matrix-multiplication和R分开points的p，则可以跳过这种情况。另外，让我们引入前面讨论的 slice 的前两行。因此，我们可以像这样获得rot_points的前两行-
rot_points_row2 = np.dot(R[:2], (points.T)) - np.dot(R[:2],p[:,None])
2)err = rot_points [0] ** 2 + rot_points [1] ** 2
可以使用
np.einsum 加快第二步，以实现高效的squaring和sum-reduction，如下所示-err = np.einsum('ij,ij->j',rot_points_row2,rot_points_row2)
对于相对较少的点
2000，计算mag:mag = np.linalg.norm(k)的步骤在运行时方面也可能变得很重要。因此，为了加快速度，可以选择再次使用np.einsum，如下所示:mag = np.sqrt(np.einsum('i,i->',k,k))
运行时测试
让我们使用
2000空间中的3D点的随机数组作为输入points，并使用原始方法和最后两行建议的方法查看关联的运行时编号。In [44]: # Setup input points
...: N = 2000
...: points = np.random.rand(N,3)
...:
...: u, v, w, x, y, z = cv2.fitLine(points, cv2.DIST_L2, 0, 1, 0.01)
...: v = np.array([u[0], v[0], w[0]])
...: p = np.array([x[0], y[0], z[0]])
...:
...: # rotate fit to z axis
...: k = np.cross(v, [0, 0, 1])
...: mag = np.linalg.norm(k)
...: R, _ = cv2.Rodrigues(k * np.arcsin(mag) / mag)
...:
...: # rotate points and calculate distance to z-axis
...: rot_points = np.dot(R, (points-p).T)
...: err = rot_points[0]**2 + rot_points[1]**2
...:
让我们运行我们提出的方法，并对照原始输出验证其输出-
In [45]: rot_points_row2 = np.dot(R[:2], (points.T)) - np.dot(R[:2],p[:,None])
...: err2 = np.einsum('ij,ij->j',rot_points_row2,rot_points_row2)
...:
In [46]: np.allclose(rot_points[:2],rot_points_row2)
Out[46]: True
In [47]: np.allclose(err,err2)
Out[47]: True
最后也是最重要的一点，让我们时间介绍这些代码部分-
In [48]: %timeit np.dot(R, (points-p).T) # Original code
10000 loops, best of 3: 79.5 µs per loop
In [49]: %timeit np.dot(R[:2], (points.T)) - np.dot(R[:2],p[:,None]) # Proposed
10000 loops, best of 3: 49.7 µs per loop
In [50]: %timeit rot_points[0]**2 + rot_points[1]**2 # Original code
100000 loops, best of 3: 12.6 µs per loop
In [51]: %timeit np.einsum('ij,ij->j',rot_points_row2,rot_points_row2) # Proposed
100000 loops, best of 3: 11.7 µs per loop
当点数增加到巨大数量时，运行时看起来更有希望。有了
N = 5000000点，我们得到-In [59]: %timeit np.dot(R, (points-p).T) # Original code
1 loops, best of 3: 410 ms per loop
In [60]: %timeit np.dot(R[:2], (points.T)) - np.dot(R[:2],p[:,None]) # Proposed
1 loops, best of 3: 254 ms per loop
In [61]: %timeit rot_points[0]**2 + rot_points[1]**2 # Original code
10 loops, best of 3: 144 ms per loop
In [62]: %timeit np.einsum('ij,ij->j',rot_points_row2,rot_points_row2) # Proposed
10 loops, best of 3: 77.5 ms per loop