OpenCV:从decomposeHomographyMat函数得到的旋转和平移矩阵奇怪

5

我正在使用OpenCV库中的findHomography函数尝试解决平面到平面的投影问题。作为一个玩具示例,我有一组R2中的点P和另一组R2中的点Q,其中Qx = Px+50,Qy = Py。这意味着我将x坐标偏移了50个单位。现在运行以下代码:

Mat projectionMatrix = findHomography(Q, P);
vector<Point2f> projectedPoints(objectCoordinates.size());
perspectiveTransform(Q, projectedPoints, projectionMatrix);

这个程序给了我P,很不错。但是,我现在想要旋转和平移矩阵R&T,这就让我感到困惑了。OpenCV 3有一个叫做decomposeHomographyMat的函数,可以返回R和T的最多4个解(也返回法向量,但我不存储它们)。我按以下方式运行它:

vector<Mat> Rs;
vector<Mat> Ts;
decomposeHomographyMat(projectionMatrix, cameraMatrix, Rs, Ts, noArray());

我使用的 cameraMatrix 是从先前的实验中经过试验和测试的。好的,我得到了四个结果。看着它们,我注意到对于所有R,我得到了恒等矩阵,这很好。但是,所有的转换向量都是[0,0,0]T,而我希望至少有一个是[-50,0,0]T。我需要对decomposeHomographyMat的结果进行某些操作才能获得期望的行为吗?
谢谢

你是否验证了计算得出的单应性矩阵不为零? - SpamBot
@SpamBot 是的,计算出的单应性矩阵是(进行了一些非常小的四舍五入处理):[1 0 -50; 0 1 0; 0 0 1]。 - IlliteratePhD
这只是一个猜测,也许decomposeHomographyMat无法处理仿射单应性的特殊情况? - SpamBot
2个回答

15

我发现我的某些观点是错误的,因此我决定重写这个答案。

简而言之 - 您会得到奇怪的结果,是因为内在参数矩阵不正确。

使用论文“Malis、E和Vargas、M,”基于视觉控制的单应性分解的深入理解“(OpenCV中的单应性分解基于该论文的术语),透视变换用H表示,并称为欧几里得单应性矩阵,其归一化结果G = K^-1 * H * K(其中K是相机的校准矩阵)称为单应性矩阵

cv :: findHomography()cv :: decomposeHomographyMat()都使用欧几里得单应性矩阵。 但是,为了将其分解为平移和旋转,cv :: decomposeHomographyMat()欧几里得单应性矩阵归一化以获得单应性矩阵。 它依靠用户提供的K来执行这种归一化。

关于估计K的问题,我认为这超出了这个问题的范围。这个问题被称为相机自动标定,以下是来自维基百科文章的相关引用:
因此,至少需要三个视角进行完整校准,视角之间的内部参数固定。优质的现代成像传感器和光学器件还可以提供关于校准的更多先验约束条件,如零偏移(正交像素网格)和单位纵横比(正方形像素)。集成这些先验条件将减少所需图像的最小数量至两张。
看起来,您可以在零偏移和正方形像素的假设下,从同一台相机的两帧图像对中提取K。但是,我对这个主题不熟悉,无法给出更多建议。
因此,为了检查我的理解是否正确,我制作了一个小例子,在其中将一些点在3D平面上投影到2个虚拟相机上,找到单应性矩阵并对其进行分解,然后将此分解与真实的旋转和平移向量进行比较。这比使用真实输入更好,因为这样我们可以精确地知道K,并且可以将其估计误差与Rt的误差分离开来。对于我检查过的输入,它能够正确地估计旋转和平移向量,尽管由于某种原因平移向量始终比真实值小10倍。也许这个分解只定义了一个比例(现在我不确定),但有趣的是,它与真实值之间的关系是由一个固定系数确定的。
以下是源代码:
#include <opencv2/opencv.hpp>
#include <iostream>
#include <vector>


int main() {
  // set up a virtual camera
  float f = 100, w = 640, h = 480;

  cv::Mat1f K = (cv::Mat1f(3, 3) <<
      f, 0, w/2,
      0, f, h/2,
      0, 0,   1);

  // set transformation from 1st to 2nd camera (assume K is unchanged)
  cv::Mat1f rvecDeg = (cv::Mat1f(3, 1) << 45, 12, 66);
  cv::Mat1f t = (cv::Mat1f(3, 1) << 100, 200, 300);

  std::cout << "-------------------------------------------\n";
  std::cout << "Ground truth:\n";

  std::cout << "K = \n" << K << std::endl << std::endl;
  std::cout << "rvec = \n" << rvecDeg << std::endl << std::endl;
  std::cout << "t = \n" << t << std::endl << std::endl;

  // set up points on a plane
  std::vector<cv::Point3f> p3d{{0, 0, 10},
                               {100, 0, 10},
                               {0, 100, 10},
                               {100, 100, 10}};

  // project on both cameras
  std::vector<cv::Point2f> Q, P, S;

  cv::projectPoints(p3d,
                    cv::Mat1d::zeros(3, 1),
                    cv::Mat1d::zeros(3, 1),
                    K,
                    cv::Mat(),
                    Q);

  cv::projectPoints(p3d,
                    rvecDeg*CV_PI/180,
                    t,
                    K,
                    cv::Mat(),
                    P);

  // find homography
  cv::Mat H = cv::findHomography(Q, P);

  std::cout << "-------------------------------------------\n";
  std::cout << "Estimated H = \n" << H << std::endl << std::endl;


  // check by reprojection
  std::vector<cv::Point2f> P_(P.size());
  cv::perspectiveTransform(Q, P_, H);

  float sumError = 0;

  for (size_t i = 0; i < P.size(); i++) {
    sumError += cv::norm(P[i] - P_[i]);
  }

  std::cout << "-------------------------------------------\n";
  std::cout << "Average reprojection error = "
      << sumError/P.size() << std::endl << std::endl;


  // decompose using identity as internal parameters matrix
  std::vector<cv::Mat> Rs, Ts;
  cv::decomposeHomographyMat(H,
                             K,
                             Rs, Ts,
                             cv::noArray());

  std::cout << "-------------------------------------------\n";
  std::cout << "Estimated decomposition:\n\n";
  std::cout << "rvec = " << std::endl;
  for (auto R_ : Rs) {
    cv::Mat1d rvec;
    cv::Rodrigues(R_, rvec);
    std::cout << rvec*180/CV_PI << std::endl << std::endl;
  }

  std::cout << std::endl;

  std::cout << "t = " << std::endl;
  for (auto t_ : Ts) {
    std::cout << t_ << std::endl << std::endl;
  }

  return 0;
}

这是在我的机器上的输出:

-------------------------------------------
Ground truth:
K =
[100, 0, 320;
0, 100, 240;
0, 0, 1]

rvec =
[45;
12;
66]

t =
[100;
200;
300]

-------------------------------------------
Estimated H =
[0.04136041220427821, 0.04748763375951008, 358.5557917287962;
0.05074854454707714, 0.06137211243830468, 297.4585754092336;
8.294458769850147e-05, 0.0002294875562580223, 1]

-------------------------------------------
Average reprojection error = 0

-------------------------------------------
Estimated decomposition:

rvec =
[-73.21470385654712;
56.64668212487194;
82.09114210289061]

[-73.21470385654712;
56.64668212487194;
82.09114210289061]

[45.00005330430893;
12.00000697952995;
65.99998380038915]

[45.00005330430893;
12.00000697952995;
65.99998380038915]


t =
[10.76993852870029;
18.60689642878277;
30.62344129378435]

[-10.76993852870029;
-18.60689642878277;
-30.62344129378435]

[10.00001378255982;
20.00002581449634;
30.0000336510648]

[-10.00001378255982;
-20.00002581449634;
-30.0000336510648]

正如您所见,假设中存在旋转向量的正确估计,并且存在比例正确的平移估计。

1
谢谢您提供详尽和有用的答案!我已经阅读了您提供的一些参考资料,我认为我理解了一般的概念。但是,我仍然不确定K的作用,因为它似乎被完全忽略了。然而,当我将上述工作流应用于相机旋转30度和标记点的两张图片时,我的结果很糟糕。同样,使用估计的H进行透视变换可以得到完美的结果。但是分解却产生了糟糕的结果。我猜测在应用R和t时我可能错过了某种归一化点的方法,但无法找出具体原因。 - IlliteratePhD
我稍微解释一下。我已经在两个场景的重叠部分手动放置了12个标记,摄像机围绕y轴旋转了30度。由于摄像机只围绕一个轴旋转,我原本以为会得到一个展示旋转的R和t = 0。然而,我看到的是相当大的平移和毫无意义的旋转。 - IlliteratePhD
现在,我不太确定当您使用cv :: findHomography()来估计转换时,K矩阵是否应该是单位矩阵。坐标确实在幕后进行了归一化,但归一化以矩阵形式存储,并且在返回结果之前将其逆应用于估计的变换。因此,cv :: findHomography()似乎返回欧几里得单应性矩阵。但不幸的是,我无法弄清楚如何在这种情况下找到相应的K矩阵。 - alexisrozhkov
我已经编辑了答案,虽然它并没有完全解决你的问题,但我希望它仍然可能会有所帮助。 - alexisrozhkov
非常感谢您分享这个信息丰富的答案。“有一个按比例正确估计翻译”的说法,我想问一下是否可能获得真实的转换比例? - Humam Helfawi
显示剩余5条评论

2
您应该在每个平面上使用solvePnP,然后从两个相机矩阵中获取相对平移和旋转(假设您至少有4个点)。decomposeHomographyMat的问题是你可能会得到多达4种解决方案...您可以删除超出图像FOV的两个解决方案,但这仍然是一个问题。
  • @alexisrozhkov- 我已经找到了您代码中的错误——文章和函数都假定Z平面在1处(而不是像您一样的z=10)。以下是您上面代码的正确实现(修正并用Python编写):
import cv2
import numpy as np


# set up a virtual camera
f = 100
w = 640
h = 480

K = np.array([[f, 0, w/2],
              [0, f, h/2],
              [0, 0,   1]])
dist_coffs = np.zeros(5)
# set transformation from 1st to 2nd camera (assume K is unchanged)
rvecDeg = np.array([[45, 12, 66]]).T
t = np.array([[100.0, 200, 300]]).T

print("-------------------------------------------\n")
print("Ground truth:\n")

print("K = \n" + str(K))
print("rvec = \n" + str(rvecDeg))
print("t = \n" + str(t))

# set up points on a plane
p3d = np.array([[0, 0, 1],
                [100, 0, 1],
                [0, 100, 1],
                [100, 100, 1]], dtype=np.float64)

# project on both cameras

Q, _ = cv2.projectPoints(p3d,
                         np.zeros((3, 1)),
                         np.zeros((3, 1)),
                         K,
                         dist_coffs)

P, _ = cv2.projectPoints(p3d,
                         rvecDeg*np.pi/180,
                         t,
                         K,
                         dist_coffs)

# find homography
H, _ = cv2.findHomography(Q, P)

print("-------------------------------------------\n")
print("Estimated H = \n" + str(H))


# check by reprojection
P_ = cv2.perspectiveTransform(Q, H)

sumError = 0

for i in range(P.shape[0]):
    sumError += np.linalg.norm(P[i] - P_[i])


print("-------------------------------------------\n")
print("Average reprojection error = "+str(sumError/P.shape[0]))


# decompose using identity as internal parameters matrix
num_res, Rs, ts, n = cv2.decomposeHomographyMat(H, K)

print("-------------------------------------------\n")
print("Estimated decomposition:\n\n")
for i, Rt in enumerate(zip(Rs, ts)):
    R, t = Rt
    print("option " + str(i+1))
    print("rvec = ")
    rvec, _ = cv2.Rodrigues(R)
    print(rvec*180/np.pi)
    print("t = ")
    print(t)

这是结果(选项3被确定为正确答案):
-------------------------------------------

Ground truth:

K = 
[[100.   0. 320.]
 [  0. 100. 240.]
 [  0.   0.   1.]]
rvec = 
[[45]
 [12]
 [66]]
t = 
[[100.]
 [200.]
 [300.]]
-------------------------------------------

Estimated H = 
[[3.93678513e-03 4.51998690e-03 3.53830416e+02]
 [4.83037187e-03 5.84154353e-03 3.05790229e+02]
 [7.89487379e-06 2.18431675e-05 1.00000000e+00]]
-------------------------------------------

Average reprojection error = 1.1324020050061362e-05
-------------------------------------------

Estimated decomposition:


option 1
rvec = 
[[-78.28555877]
 [ 58.01301837]
 [ 81.41634175]]
t = 
[[100.79816988]
 [198.59277542]
 [300.6675498 ]]
option 2
rvec = 
[[-78.28555877]
 [ 58.01301837]
 [ 81.41634175]]
t = 
[[-100.79816988]
 [-198.59277542]
 [-300.6675498 ]]
option 3
rvec = 
[[45.0000205 ]
 [12.00005488]
 [65.99999505]]
t = 
[[100.00010826]
 [200.00026791]
 [300.00034698]]
option 4
rvec = 
[[45.0000205 ]
 [12.00005488]
 [65.99999505]]
t = 
[[-100.00010826]
 [-200.00026791]
 [-300.00034698]]

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接