OpenCV：使用solvePnP确定单应性

Question

在过去的几周里，我尝试学习如何纠正图像，在这里大家的帮助下，我对这个过程有了更好的理解。大约一周前，我设置了一个测试示例，想要进行图像纠正（可以从上方查看图像）。这个过程运行得很好 (原始图像: http://sitedezign.net/original.jpg，纠正后的图像: http://sitedezign.net/rectified.jpg)，使用的函数是 T = cv2.getPerspectiveTransform(UV_cp, XYZ_gcp)，这里的 T 就是我们所说的单应性。

但是，当我尝试用现实世界的照片进行同样的操作时，结果却失败了，因为现实世界的坐标并不完全在一个平面上（而是大约有10个控制点，这些点在空间中有X、Y和Z坐标）。因此，我决定使用 solvePnP，希望能够创建一个可以使用的单应性。

我在测试示例上尝试了这个方法，但没有得到预期的结果：图像没有被纠正，而我用 solvePnP 计算出的单应性与使用 getPerspectiveTransform 计算的单应性不相等。

我的代码如下：

# Set UV (image) and XYZ (real life)
UV_cp = np.array([[1300.0, 2544.0], # left down
                  [1607.0, 1000.0], # left up
                  [3681.0, 2516.0], # right down
                  [3320.0, 983.0]], np.float32) # right up

# Z is on 0 plane, so Z=0.0
XYZ_gcp = np.array([[0.0, 400.0, 0.0],
                    [0.0, 0.0, 0.0],
                    [300.0, 400.0, 0.0],
                    [300.0, 0.0, 0.0]], np.float32)

rvec, tvec = cv2.solvePnP(XYZ_gcp, UV_cp, K, D)
rotM_cam = cv2.Rodrigues(rvec)[0]

# calculate camera position (= translation), in mm from 0,0,0 point
cameraPosition = -np.matrix(rotM_cam).T * np.matrix(tvec)

# 3x3 Identity matrix
I = np.identity(3)

# [I|-C]
I1_extended = np.hstack((I,-cameraPosition))

# P = K*R*I
P_cam = K.dot(rotM_cam).dot(I1_extended)

# create P2 = image from above: R = 0,0,0, translation = x, y, z = 0,0,-1000 (mm)
R_rec = matr.getR(0.0,0.0,0.0)
newZ = -1000.0
new_cameraPosition = np.array([[0.0],[0.0],[newZ]])
I2_extended = np.hstack((I,new_cameraPosition))
P_rec = K.dot(R_rec).dot(I2_extended)

# correct Homography T from getPerspectiveTransform:
T = np.array([[4.70332834e-01, 9.35182514e-02, -4.24671558e+02],
              [9.62104844e-03, 9.69462117e-01, -4.92461571e+02],
              [3.54859924e-06, 6.80081146e-04, 1.00000000e+00]])

# Homography Matrix = H = P_rect * pinv(P) => P2 * pinv(P1)
H = P_rec.dot(np.linalg.pinv(P_cam))

结果是一个变形的图像，和上面显示的纠正后的图像相差甚远。此外，应该是正确的单应性 T（来自 getPerspectiveTransform）与使用 solvePnP 结果计算的单应性 H 也相差很大。

H from solvePnP:
[[  1.01865631e+00   2.68683332e-01  -2.04519580e+03]
 [ -3.24304366e-02   6.82672680e-01  -1.15688010e+03]
 [  2.03399902e-05   1.24191993e-04  -5.41378561e-01]]

H from getPerspectiveTransform:
[[  4.70332834e-01   9.35182514e-02  -4.24671558e+02]
 [  9.62104844e-03   9.69462117e-01  -4.92461571e+02]
 [  3.54859924e-06   6.80081146e-04   1.00000000e+00]]

有没有人知道哪里出了问题？

附注：用于确定 K 矩阵和畸变系数的代码（这些值是根据我的相机 Pentax K-5 在33mm焦距下通过 Adobe Camera Raw 获取的）：

# Focal length, sensor size (mm and px)
f = 33.0 # mm
pix_width = 4928.0 # sensor size has 4928px in width
pix_height = 3624.0 # sensor size has 4928px in width
sensor_width = 23.7 # mm
sensor_height = 15.7 # mm

# set center pixel
u0 = int(pix_width / 2.0)
v0 = int(pix_height / 2.0)

# determine values of camera-matrix
mu = pix_width / sensor_width # px/mm
alpha_u = f * mu # px

mv = pix_height / sensor_height # px/mm
alpha_v = f * mv # px

# Distortion coefs 
D = np.array([[0.0, 0.0, 0.0, 0.0]])

# Camera matrix
K = np.array([[alpha_u, 0.0, u0],
              [0.0, alpha_v, v0],
              [0.0, 0.0, 1.0]])

Answer 1

你的 K 矩阵看起来是合适的，但这可能不足以在真实图像中获得良好的准确性。我认为，与其给出一些合理的值（特别是光学中心像素和镜头畸变系数），不如使用 calibrateCamera 函数来校准你的相机（文档链接, 教程）。不过，我觉得你提到的问题并不是由这个引起的。

我认为你的问题出在 P_rec 的定义上。

首先，要注意如果你使用 newZ = -1000.0，实际上你是在将相机移动1000米（而不是毫米）。

其次，你需要非常小心你考虑的3D点，以及你希望它们在图像中投影的位置：

1. 因为你在 solvePnP 函数中使用了 XYZ_gcp，这意味着你将这些坐标作为3D点使用。

2. 因为你在 getPerspectiveTransform 函数中使用了 XYZ_gcp，这意味着你也将这些作为2D坐标使用。注意，严格来说，你不能这样做，因为 getPerspectiveTransform 期望两个4x2的数组（而不是一个4x2和一个4x3），但我假设你省略了第三个坐标，它们总是为0。

因此，你的 P_rec 应该定义为 [x; y; 1] = P_rec * [x; y; 0; 1]。所以，P_rec 应该这样定义：

P_rec = [ [1 0 0 0] [0 1 0 0] [0 0 0 1] ].