我正在嘗試加速一些依賴于每個相機對之間基本矩陣計算的多相機系統。
請注意以下是偽代碼。@表示矩陣乘法，|表示串聯。
我有代碼來計算F每一對calculate_f(camera_matrix1_3x4, camera_matrix1_3x4)，而天真的解決方案是
for c1 in cameras:
for c2 in cameras:
if c1 != c2:
f = calculate_f(c1.proj_matrix, c2.proj_matrix)
這很慢，我想加快速度。我有大約 5000 個攝像頭。
我已經預先計算了每對相機之間的所有旋轉和平移（在世界坐標中）以及內部引數k，因此對于每個相機c，它都認為c.matrix = c.k @ (c.rot | c.t)
我可以使用這些引數r, t來幫助加速以下計算F嗎？
在數學形式，對3個不同的攝像頭c1，c2，c3我有
f12=(c1.proj_matrix, c2.proj_matrix)，我想f23=(c2.proj_matrix, c3.proj_matrix)，f13=(c1.proj_matrix, c3.proj_matrix)有一些功能f23, f13 = fast_f(f12, c1.r, c1.t, c2.r, c2.t, c3.r, c3.t)？
用于計算 numpy 中基本矩陣的作業函式：
def fundamental_3x3_from_projections(p_left_3x4: np.array, p_right__3x4: np.array) -> np.array:
# The following is based on OpenCv-contrib's c implementation.
# see https://github.com/opencv/opencv_contrib/blob/master/modules/sfm/src/fundamental.cpp#L109
# see https://sourishghosh.com/2016/fundamental-matrix-from-camera-matrices/
# see https://answers.opencv.org/question/131017/how-do-i-compute-the-fundamental-matrix-from-2-projection-matrices/
f_3x3 = np.zeros((3, 3))
p1, p2 = p_left_3x4, p_right__3x4
x = np.empty((3, 2, 4), dtype=np.float)
x[0, :, :] = np.vstack([p1[1, :], p1[2, :]])
x[1, :, :] = np.vstack([p1[2, :], p1[0, :]])
x[2, :, :] = np.vstack([p1[0, :], p1[1, :]])
y = np.empty((3, 2, 4), dtype=np.float)
y[0, :, :] = np.vstack([p2[1, :], p2[2, :]])
y[1, :, :] = np.vstack([p2[2, :], p2[0, :]])
y[2, :, :] = np.vstack([p2[0, :], p2[1, :]])
for i in range(3):
for j in range(3):
xy = np.vstack([x[j, :], y[i, :]])
f_3x3[i, j] = np.linalg.det(xy)
return f_3x3
uj5u.com熱心網友回復：
Numpy 顯然沒有針對小矩陣進行優化。CPython 輸入物件的決議、內部檢查和函式呼叫引入了顯著的開銷，遠大于執行實際計算所需的執行時間。更不用說創建許多臨時陣列也很昂貴。解決此問題的一種方法是使用 Numba 或 Cython。
此外，行列式的計算可以優化很多，因為您知道矩陣的確切大小并且矩陣的一部分并不總是改變。事實上，使用基本代數運算式4×4行列式幫助編譯器優化了很多整體計算得益于公共子運算式洗掉（不CPython的解釋執行）和去除復雜的回圈/條件句np.linalg.det。
這是結果代碼：
import numba as nb
@nb.njit('float64(float64[:,::1])')
def det_4x4(mat):
a, b, c, d = mat[0,0], mat[0,1], mat[0,2], mat[0,3]
e, f, g, h = mat[1,0], mat[1,1], mat[1,2], mat[1,3]
i, j, k, l = mat[2,0], mat[2,1], mat[2,2], mat[2,3]
m, n, o, p = mat[3,0], mat[3,1], mat[3,2], mat[3,3]
return a * (f * (k*p - l*o) g * (l*n - j*p) h * (j*o - k*n)) \
b * (e * (l*o - k*p) g * (i*p - l*m) h * (k*m - i*o)) \
c * (e * (j*p - l*n) f * (l*m - i*p) h * (i*n - j*m)) \
d * (e * (k*n - j*o) f * (i*o - k*m) g * (j*m - i*n))
@nb.njit('float64[:,::1](float64[:,::1], float64[:,::1])')
def fundamental_3x3_from_projections(p_left_3x4, p_right_3x4):
f_3x3 = np.empty((3, 3))
p1, p2 = p_left_3x4, p_right_3x4
x = np.empty((3, 2, 4), dtype=np.float64)
x[0, 0, :] = p1[1, :]
x[0, 1, :] = p1[2, :]
x[1, 0, :] = p1[2, :]
x[1, 1, :] = p1[0, :]
x[2, 0, :] = p1[0, :]
x[2, 1, :] = p1[1, :]
y = np.empty((3, 2, 4), dtype=np.float64)
y[0, 0, :] = p2[1, :]
y[0, 1, :] = p2[2, :]
y[1, 0, :] = p2[2, :]
y[1, 1, :] = p2[0, :]
y[2, 0, :] = p2[0, :]
y[2, 1, :] = p2[1, :]
xy = np.empty((4, 4), dtype=np.float64)
for i in range(3):
xy[2:4, :] = y[i, :, :]
for j in range(3):
xy[0:2, :] = x[j, :, :]
f_3x3[i, j] = det_4x4(xy)
return f_3x3
這在我的機器上快了130 倍（85.6 us VS 0.66 us）。
如果應用的函式是可交換的（即f(c1, c2) == f(c2, c1)），您可以將程序加快兩倍。如果是這樣，您只能計算上半部分。事實證明，您的函式具有一些有趣的屬性，因為它f(c1, c2) == f(c2, c1).T似乎總是正確的。另一種可能的優化是并行運行回圈。
通過所有這些優化，生成的程式應該快大約3 個數量級。
分析方法的準確性
提供的精度似乎與原始精度相似。關于輸入矩陣，結果有時比 Numpy 方法更準確，有時更不準確。這特別是由于行列式的計算。使用 24 位小數，與Wolphram Alpha的可靠結果相比，沒有明顯的錯誤。這表明該方法是正確的，由于數值穩定性細節，結果不一樣。以下是用于測驗方法準確性的代碼：
# Imports
from decimal import Decimal
import numba as nb
# Definitions
def det_4x4(mat):
a, b, c, d = mat[0,0], mat[0,1], mat[0,2], mat[0,3]
e, f, g, h = mat[1,0], mat[1,1], mat[1,2], mat[1,3]
i, j, k, l = mat[2,0], mat[2,1], mat[2,2], mat[2,3]
m, n, o, p = mat[3,0], mat[3,1], mat[3,2], mat[3,3]
return a * (f * (k*p - l*o) g * (l*n - j*p) h * (j*o - k*n)) \
b * (e * (l*o - k*p) g * (i*p - l*m) h * (k*m - i*o)) \
c * (e * (j*p - l*n) f * (l*m - i*p) h * (i*n - j*m)) \
d * (e * (k*n - j*o) f * (i*o - k*m) g * (j*m - i*n))
@nb.njit('float64(float64[:,::1])')
def det_4x4_numba(mat):
a, b, c, d = mat[0,0], mat[0,1], mat[0,2], mat[0,3]
e, f, g, h = mat[1,0], mat[1,1], mat[1,2], mat[1,3]
i, j, k, l = mat[2,0], mat[2,1], mat[2,2], mat[2,3]
m, n, o, p = mat[3,0], mat[3,1], mat[3,2], mat[3,3]
return a * (f * (k*p - l*o) g * (l*n - j*p) h * (j*o - k*n)) \
b * (e * (l*o - k*p) g * (i*p - l*m) h * (k*m - i*o)) \
c * (e * (j*p - l*n) f * (l*m - i*p) h * (i*n - j*m)) \
d * (e * (k*n - j*o) f * (i*o - k*m) g * (j*m - i*n))
# Example matrix
precise_xy = np.array(
[[Decimal('42'),Decimal('-6248'),Decimal('4060'),Decimal('845')],
[Decimal('-0.00992'),Decimal('-0.704'),Decimal('-0.71173298417'),Decimal('300.532')],
[Decimal('-8.94274'),Decimal('-7554.39'),Decimal('604.57'),Decimal('706282')],
[Decimal('-0.0132'),Decimal('-0.2757'),Decimal('-0.961'),Decimal('247.65')]]
)
xy = precise_xy.astype(np.float64)
res_numpy = Decimal(np.linalg.det(xy))
res_numba = Decimal(det_4x4_numba(xy))
res_precise = det_4x4(precise_xy)
# The Wolphram Alpha expression used is:
# det({{42,-6248,4060,845},
# {-0.00992,-0.704,-0.71173298417,300.532},
# {-8.94274,-7554.39,604.57,706282},
# {-0.0132,-0.2757,-0.961,247.65}})
res_wolframalpha = Decimal('-323312.2164828991329828243')
# The result got from Wolfram-Alpha have a 25-digit precision
# and is exactly the same than the one of det_4x4 using 24-digit decimals.
assert res_precise == res_wolframalpha
print(abs((res_numpy-res_precise)/res_precise)) # 1.7E-14
print(abs((res_numba-res_precise)/res_precise)) # 3.1E-14
# => Similar relative error (Numba slightly less accurate
# but both are not close to the 1e-16 relative epsilon)
轉載請註明出處，本文鏈接：https://www.uj5u.com/gongcheng/402193.html
上一篇：理解中的兩個命令