相机纯横摇（Pan）旋转时两图像间单应性矩阵H的求解疑问

最终问题 #

我最近在研究相机纯横摇旋转时，两张图像之间的单应性矩阵H的计算方法，遇到了一个困惑：H是3×3矩阵，其中H[0,0]、H[0,1]、H[1,0]、H[1,1]看起来可以通过旋转矩阵推导得到，但H[0,2]和H[1,2]这两个元素该怎么计算呢？

理论背景 #

先给大家梳理下我目前用到的理论基础：
给定相机内参矩阵（f为焦距）：

K = np.diag([f, f, 1])

相机中心在世界坐标系中的位置：

Tw = np.array([tx, ty, tz])

构建矩阵ITw：

ITw = np.column_stack((np.eye(3), Tw))

可以计算投影矩阵P：

P = K @ Rot @ ITw

其中Rot是3×3的旋转矩阵。对于世界坐标系中Z=0平面上的点，对应的单应性矩阵可以通过下式计算：

Hcalc = P[:,[0,1,3]]

旋转矩阵Rot的定义方式如下（用scipy基于欧拉角计算）：

from scipy.spatial.transform import Rotation as R
pan = 10
tilt = -30
roll = 4

pan *= np.pi / 180
tilt *= np.pi / 180
roll *= np.pi / 180

r_pan = R.from_euler('Z', pan).as_matrix()
r_tilt = R.from_euler('X', tilt).as_matrix()
r_roll = R.from_euler('Z', roll).as_matrix()

Rot = r_roll @ r_tilt @ r_pan

示例验证 #

我用一个具体例子验证了前面的理论：
取焦距f=400，相机中心位置Tw=[10, -5, 6]，地面上（Z=0）的一组源点：

src_pts = np.array([ 
       [ 0.  , -7.62],
       [ 0.  , -6.35],
       [ 0.  ,  0.  ],
       [ 0.  ,  1.27],
       [ 0.  ,  2.54],
       [ 0.  ,  3.81],
       [20.  , -7.62],
       [20.  , -6.35],
       [20.  ,  0.  ],
       [20.  ,  1.27],
       [20.  ,  2.54],
       [20.  ,  3.81]])

我通过Hcalc计算出这些点在图像上的对应点pts10，并用OpenCV的findHomography验证：

h10, status = cv2.findHomography(src_pts.reshape(-1,1,2), pts10.reshape(-1,1,2))

得到的结果和Hcalc完全一致。之后我把pan设为15°，重复同样步骤得到h15和pts15。

我的困惑点 #

现在我需要计算从第一张图像（pan=10°）到第二张图像（pan=15°）的单应性矩阵，用OpenCV得到的结果是：

H
Out[187]: 
array([[ 9.97147963e-01, -7.56172873e-02,  1.73592176e+01],
       [ 7.54846993e-02,  9.98090883e-01,  1.87521001e+00],
       [-1.09070429e-04, -3.49361263e-06,  1.00000000e+00]])

我分别计算了pan=10°和pan=15°对应的旋转矩阵rot10和rot15，然后计算：

rot15 @ rot10.T / (rot15 @ rot10.T)[2,2]
Out[191]: 
array([[ 0.99714794, -0.07561722,  0.04339806],
       [ 0.0754847 ,  0.9980909 ,  0.00468803],
       [-0.04362816, -0.00139744,  1.        ]])

从这里我找到了H部分元素的计算规律：
令rotNorm = rot15 @ rot10.T / (rot15 @ rot10.T)[2,2]，那么对于i,j < 2，有H[i,j] = RotNorm[i,j]；
另外H的第三行前两列满足：H[2,j] = RotNorm[2,j] / f（j < 2）。

但是H[0,2]和H[1,2]这两个元素该怎么推导计算呢？

非常感谢大家能看到这里，希望有人能帮我解答这个疑问！

备注：内容来源于stack exchange，提问作者fenaux