上一个教程 ： AKAZE 和 ORB 平面跟踪
| 兼容性 | OpenCV>=3.0 |
简介
本教程将通过一些代码演示同源性的基本概念。有关理论的详细解释，请参阅计算机视觉课程或计算机视觉书籍，如《计算机视觉中的多视图几何》（Multiple View Geometry in Computer Vision）：
- 《计算机视觉中的多视图几何》，理查德-哈特利（Richard Hartley）和安德鲁-齐瑟曼（Andrew Zisserman），[107]（此处有部分章节示例，此处有 CVPR 教程）。
- 《三维视觉邀请函： 从图像到几何模型》，Yi Ma、Stefano Soatto、Jana Kosecka 和 S. Shankar Sastry，[165]（在此可获得计算机视觉书籍讲义）
- 《计算机视觉： 算法与应用》，Richard Szeliski，[244]（电子版可在此处获取）
- 深入理解基于视觉的控制的同构分解，Ezio Malis、Manuel Vargas，[168] （可在此处开放访问）
增强现实的姿势估计： 实践调查，Eric Marchand、Hideaki Uchiyama、Fabien Spindler，[170] （此处可开放访问）
教程代码可在此处找到 C++、Python 和 Java。本教程中使用的图片可在此处找到（left*.jpg）。
基本理论
什么是同构矩阵？
简而言之，平面同构关系涉及两个平面之间的变换（不超过比例因子）：
s
[
x
′
y
′
1
]
=
H
[
x
y
1
]
=
[
h
11
h
12
h
13
h
21
h
22
h
23
h
31
h
32
h
33
]
[
x
y
1
]
s \begin{bmatrix} x^{'} \\ y^{'} \\ 1 \end{bmatrix} = \mathbf{H} \begin{bmatrix} x \\ y \\ 1 \end{bmatrix} = \begin{bmatrix} h_{11} & h_{12} & h_{13} \\ h_{21} & h_{22} & h_{23} \\ h_{31} & h_{32} & h_{33} \end{bmatrix} \begin{bmatrix} x \\ y \\ 1 \end{bmatrix}
s
x′y′1
=H
xy1
=
h11h21h31h12h22h32h13h23h33
xy1
同构矩阵是一个 3x3 矩阵，但有 8 个 DoF（自由度），因为它是按比例估算的。它通常被归一化（另见 1）为
h
33
=
1
h_{33} = 1
h33=1 或
h
11
2
+
h
12
2
+
h
13
2
+
h
21
2
+
h
22
2
+
h
23
2
+
h
31
2
+
h
32
2
+
h
33
2
=
1
h_{11}^2 + h_{12}^2 + h_{13}^2 + h_{21}^2 + h_{22}^2 + h_{23}^2 + h_{31}^2 + h_{32}^2 + h_{33}^2 = 1
h112+h122+h132+h212+h222+h232+h312+h322+h332=1
下面的例子展示了不同种类的变换，但都与两个平面之间的变换有关。
-
一个平面和图像平面（图像取自 2）
-
从两个摄像机位置观察的平面（图像取自 3 和 2）
-
摄像机绕其投影轴旋转，相当于认为各点位于无限远处的平面上（图片取自 2）
同轴变换如何发挥作用？
-
例如，通过共面点估算相机姿态，用于带有标记的增强现实技术（参见前面的第一个示例）
-
透视消除/校正（见前面第二个例子）
-
全景拼接（见前面的第二和第三个示例）
演示代码
演示 1：从共面点估算姿势
注意事项
请注意，根据同轴度估计摄像机姿态的代码只是一个示例，如果您想估计平面或任意物体的摄像机姿态，请使用 cv::solvePnP。
例如，可以使用直接线性变换（DLT）算法估算同轴度（更多信息请参见第 1 节）。由于物体是平面的，因此在物体帧中表示的点与在归一化摄像机帧中表示的投影到图像平面中的点之间的变换就是同轴变换。正因为物体是平面的，所以在已知摄像机固有参数的前提下（见 2 或 4），可以从同轴变换中获取摄像机姿态。使用国际象棋棋盘对象和 findChessboardCorners() 来获取图像中的角位置，可以很容易地测试这一点。
要检测棋盘角，首先需要棋盘大小（patternSize），这里是 9x6：
vector<Point2f> corners；
bool found = findChessboardCorners(img, patternSize, corners)；
只要知道棋盘正方形的大小，就可以轻松计算出对象帧中表示的对象点：
for( int i = 0; i < boardSize.height; i++ )
for( int j = 0; j < boardSize.width; j++ )
corners.push_back(Point3f(float(j*squareSize)、
float(i*squareSize)，0))；
必须删除坐标 Z=0 以用于同构估计部分：
vector<Point3f> objectPoints；
calcChessboardCorners(patternSize, squareSize, objectPoints)；
vector<Point2f> objectPointsPlanar；
for (size_t i = 0; i < objectPoints.size(); i++)
{
objectPointsPlanar.push_back(Point2f(objectPoints[i].x, objectPoints[i].y))；
}
通过角点和使用摄像机本征和畸变系数进行反向透视变换，可以计算出归一化摄像机中表示的图像点：
FileStorage fs( samples::findFile( intrinsicsPath ), FileStorage::READ)；
Mat cameraMatrix, distCoeffs； fs["camera_matrix"] >> cameraMatrix；
fs["distortion_coefficients"] >> distCoeffs；
vector<Point2f> imagePoints；
undistortPoints(corners, imagePoints, cameraMatrix, distCoeffs)；
然后就可以用以下方法估算同构了：
Mat H = findHomography(objectPointsPlanar, imagePoints)；
cout << "H:\n" << H << endl；
从同构矩阵获取姿态的快速解决方案是（见 5）：
// 归一化以确保 ||c1|| = 1
double norm = sqrt(H.at<double>(0,0)*H.at<double>(0,0) +
H.at<double>(1,0)*H.at<double>(1,0) +
H.at<double>(2,0)*H.at<double>(2,0))；
H /= norm；
Mat c1 = H.col(0)；
Mat c2 = H.col(1)；
Mat c3 = c1.cross(c2)；
Mat tvec = H.col(2)；
Mat R(3, 3, CV_64F)；
for (int i = 0; i < 3; i++)
{
R.at<double>(i,0) = c1.at<double>(i,0)；
R.at<double>(i,1) = c2.at<double>(i,0)；
R.at<double>(i,2) = c3.at<double>(i,0)；
}
X = ( X , Y , 0 , 1 ) x = P X = K [ r 1 r 2 r 3 t ] ( X Y 0 1 ) = K [ r 1 r 2 t ] ( X Y 1 ) = H ( X Y 1 ) \begin{align*} \mathbf{X} &= \left( X, Y, 0, 1 \right ) \\ \mathbf{x} &= \mathbf{P}\mathbf{X} \\ &= \mathbf{K} \left[ \mathbf{r_1} \hspace{0.5em} \mathbf{r_2} \hspace{0.5em} \mathbf{r_3} \hspace{0.5em} \mathbf{t} \right ] \begin{pmatrix} X \\ Y \\ 0 \\ 1 \end{pmatrix} \\ &= \mathbf{K} \left[ \mathbf{r_1} \hspace{0.5em} \mathbf{r_2} \hspace{0.5em} \mathbf{t} \right ] \begin{pmatrix} X \\ Y \\ 1 \end{pmatrix} \\ &= \mathbf{H} \begin{pmatrix} X \\ Y \\ 1 \end{pmatrix} \end{align*} Xx=(X,Y,0,1)=PX=K[r1r2r3t] XY01 =K[r1r2t] XY1 =H XY1
H
=
λ
K
[
r
1
r
2
t
]
K
−
1
H
=
λ
[
r
1
r
2
t
]
P
=
K
[
r
1
r
2
(
r
1
×
r
2
)
t
]
\begin{align*} \mathbf{H} &= \lambda \mathbf{K} \left[ \mathbf{r_1} \hspace{0.5em} \mathbf{r_2} \hspace{0.5em} \mathbf{t} \right ] \\ \mathbf{K}^{-1} \mathbf{H} &= \lambda \left[ \mathbf{r_1} \hspace{0.5em} \mathbf{r_2} \hspace{0.5em} \mathbf{t} \right ] \\ \mathbf{P} &= \mathbf{K} \left[ \mathbf{r_1} \hspace{0.5em} \mathbf{r_2} \hspace{0.5em} \left( \mathbf{r_1} \times \mathbf{r_2} \right ) \hspace{0.5em} \mathbf{t} \right ] \end{align*}
HK−1HP=λK[r1r2t]=λ[r1r2t]=K[r1r2(r1×r2)t]
这是一个快速的解决方案（另见 2），因为这并不能确保所得到的旋转矩阵是正交的，而刻度大致是通过将第一列归一化为 1 来估算的。
要获得适当的旋转矩阵（具有旋转矩阵的特性），可以采用极性分解或旋转矩阵的正交化（相关信息请参见 6 或 7 或 8 或 9）：
cout << "R (before polar decomposition):\n" << R << "\ndet(R): " << determinant(R) << endl;
Mat_<double> W, U, Vt;
SVDecomp(R, W, U, Vt);
R = U*Vt;
double det = determinant(R);
if (det < 0)
{
Vt.at<double>(2,0) *= -1;
Vt.at<double>(2,1) *= -1;
Vt.at<double>(2,2) *= -1;
R = U*Vt;
}
cout << "R (after polar decomposition):\n" << R << "\ndet(R): " << determinant(R) << endl;
为了检查结果，将显示以估计的相机姿态投射到图像中的对象帧：
演示 2：透视校正
在本示例中，将通过计算将源点映射到所需点的同轴度，将源图像转换为所需的透视图。下图显示了源图像（左）和我们要转换为所需棋盘视图的棋盘视图（右）。
第一步是检测源图像和预期图像中的棋盘角：
C++
vector<Point2f> corners1, corners2;
bool found1 = findChessboardCorners(img1, patternSize, corners1);
bool found2 = findChessboardCorners(img2, patternSize, corners2);
Java
MatOfPoint2f corners1 = new MatOfPoint2f(), corners2 = new MatOfPoint2f();
boolean found1 = Calib3d.findChessboardCorners(img1, new Size(9, 6), corners1 );
boolean found2 = Calib3d.findChessboardCorners(img2, new Size(9, 6), corners2 );
Python
ret1, corners1 = cv.findChessboardCorners(img1, patternSize)
ret2, corners2 = cv.findChessboardCorners(img2, patternSize)
通过以下方法可以轻松估算出同轴度：
C++
Mat H = findHomography(corners1, corners2);
cout << "H:\n" << H << endl;
Java
Mat H = new Mat();
H = Calib3d.findHomography(corners1, corners2);
System.out.println(H.dump());
Python
H, _ = cv.findHomography(corners1, corners2)
print(H)
要将源棋盘视图翘曲成所需的棋盘视图，我们使用 cv::warpPerspective
C++
Mat img1_warp;
warpPerspective(img1, img1_warp, H, img1.size());
Java
Mat img1_warp = new Mat();
Imgproc.warpPerspective(img1, img1_warp, H, img1.size());
Python
img1_warp = cv.warpPerspective(img1, H, (img1.shape[1], img1.shape[0]))
结果图像为
计算经同色法转换的源角坐标：
C++
Mat img_draw_matches;
hconcat(img1, img2, img_draw_matches);
for (size_t i = 0; i < corners1.size(); i++)
{
Mat pt1 = (Mat_<double>(3,1) << corners1[i].x, corners1[i].y, 1);
Mat pt2 = H * pt1;
pt2 /= pt2.at<double>(2);
Point end( (int) (img1.cols + pt2.at<double>(0)), (int) pt2.at<double>(1) );
line(img_draw_matches, corners1[i], end, randomColor(rng), 2);
}
imshow("Draw matches", img_draw_matches);
waitKey();
Java
Mat img_draw_matches = new Mat();
list2.add(img1);
list2.add(img2);
Core.hconcat(list2, img_draw_matches);
Point []corners1Arr = corners1.toArray();
for (int i = 0 ; i < corners1Arr.length; i++) {
Mat pt1 = new Mat(3, 1, CvType.CV_64FC1), pt2 = new Mat();
pt1.put(0, 0, corners1Arr[i].x, corners1Arr[i].y, 1 );
Core.gemm(H, pt1, 1, new Mat(), 0, pt2);
double[] data = pt2.get(2, 0);
Core.divide(pt2, new Scalar(data[0]), pt2);
double[] data1 =pt2.get(0, 0);
double[] data2 = pt2.get(1, 0);
Point end = new Point((int)(img1.cols()+ data1[0]), (int)data2[0]);
Imgproc.line(img_draw_matches, corners1Arr[i], end, RandomColor(), 2);
}
HighGui.imshow("Draw matches", img_draw_matches);
HighGui.waitKey(0);
Python
img_draw_matches = cv.hconcat([img1, img2])
for i in range(len(corners1)):
pt1 = np.array([corners1[i][0], corners1[i][1], 1])
pt1 = pt1.reshape(3, 1)
pt2 = np.dot(H, pt1)
pt2 = pt2/pt2[2]
end = (int(img1.shape[1] + pt2[0]), int(pt2[1]))
cv.line(img_draw_matches, tuple([int(j) for j in corners1[i]]), end, randomColor(), 2)
cv.imshow("Draw matches", img_draw_matches)
cv.waitKey(0)
为了检查计算的正确性，将显示匹配的线条：
演示 3：根据摄像机位移进行同轴测量
同轴度与两个平面之间的变换有关，因此可以获取相应的摄像机位移，以便从第一个平面视图切换到第二个平面视图（更多信息请参阅 [168]）。在详细介绍如何根据摄像机位移计算同轴度之前，先回顾一下摄像机姿态和同轴变换
通过函数 cv::solvePnP，可以根据三维物体点（在物体帧中表示的点）和投影二维图像点（在图像中查看的物体点）的对应关系计算摄像机姿态。需要本征参数和畸变系数（请参阅摄像机校准过程）。
s
[
u
v
1
]
=
[
f
x
0
c
x
0
f
y
c
y
0
0
1
]
[
r
11
r
12
r
13
t
x
r
21
r
22
r
23
t
y
r
31
r
32
r
33
t
z
]
[
X
O
Y
O
Z
O
1
]
s\left[ \begin{array}{c} u\\ v\\ 1\\ \end{array} \right] =\left[ \begin{matrix} f_x& 0& c_x\\ 0& f_y& c_y\\ 0& 0& 1\\ \end{matrix} \right] \left[ \begin{matrix} r_{11}& r_{12}& \begin{matrix} r_{13}& t_x\\ \end{matrix}\\ r_{21}& r_{22}& \begin{matrix} r_{23}& t_y\\ \end{matrix}\\ r_{31}& r_{32}& \begin{matrix} r_{33}& t_z\\ \end{matrix}\\ \end{matrix} \right] \left[ \begin{array}{c} X_O\\ Y_O\\ Z_O\\ 1\\ \end{array} \right]
s
uv1
=
fx000fy0cxcy1
r11r21r31r12r22r32r13txr23tyr33tz
XOYOZO1
=
K
C
M
O
[
X
O
Y
O
Z
O
1
]
=K^CM_O\left[ \begin{array}{c} X_O\\ Y_O\\ Z_O\\ 1\\ \end{array} \right]
=KCMO
XOYOZO1
K 是本征矩阵，cMo 是摄像机姿态。cv::solvePnP 的输出正是如此：rvec 是罗德里格斯(Rodrigues)旋转矢量，tvec 是平移矢量。
cMo 可以用同质形式表示，可以将物体帧中的点转换到摄像机帧中：
[
X
C
Y
C
Z
C
1
]
=
C
M
O
[
X
O
Y
O
Z
O
1
]
\left[ \begin{array}{c} X_C\\ Y_C\\ Z_C\\ 1\\ \end{array} \right] =^CM_O\left[ \begin{array}{c} X_O\\ Y_O\\ Z_O\\ 1\\ \end{array} \right]
XCYCZC1
=CMO
XOYOZO1
=
[
C
R
O
C
t
O
0
1
×
3
1
]
[
X
O
Y
O
Z
O
1
]
=\left[ \begin{matrix} ^CR_O& ^Ct_O\\ 0_{1\times 3}& 1\\ \end{matrix} \right] \left[ \begin{array}{c} X_O\\ Y_O\\ Z_O\\ 1\\ \end{array} \right]
=[CRO01×3CtO1]
XOYOZO1
=
[
r
11
r
12
r
13
t
x
r
21
r
22
r
23
t
y
r
31
r
32
r
33
t
z
0
0
0
1
]
[
X
O
Y
O
Z
O
1
]
=\left[ \begin{matrix} r_{11}& r_{12}& r_{13}& t_x\\ r_{21}& r_{22}& r_{23}& t_y\\ r_{31}& r_{32}& r_{33}& t_z\\ 0& 0& 0& 1\\ \end{matrix} \right] \left[ \begin{array}{c} X_O\\ Y_O\\ Z_O\\ 1\\ \end{array} \right]
=
r11r21r310r12r22r320r13r23r330txtytz1
XOYOZO1
通过矩阵乘法，可以轻松地将一帧中表示的点转换到另一帧中：
- c 1 M o ^{c_1}\mathbf{M}_o c1Mo 是摄像机 1 的姿势
- c 2 M o ^{c_2}\mathbf{M}_o c2Mo 是摄像机 2 的姿态 将摄像机 1 帧中的三维点转换为摄像机 2 帧中的三维点：
C 2 M C 1 = C 2 M O ⋅ O M C 1 = C 2 M O ( C 1 M O ) − 1 = [ C 2 R O C 2 t O 0 3 × 1 1 ] ⋅ [ C 1 R O T − C 1 R O T ⋅ C 1 t O 0 1 × 3 1 ] ^{C2}M_{C1}=^{C2}M_O\cdot ^OM_{C1}=^{C2}M_O\left( ^{C1}M_O \right) ^{-1}=\left[ \begin{matrix} ^{C2}R_O& ^{C2}t_O\\ 0_{3\times 1}& 1\\ \end{matrix} \right] \cdot \left[ \begin{matrix} ^{C1}R_O^T& -^{C1}R_{O}^{T}\cdot ^{C1}t_O\\ 0_{1\times 3}& 1\\ \end{matrix} \right] C2MC1=C2MO⋅OMC1=C2MO(C1MO)−1=[C2RO03×1C2tO1]⋅[C1ROT01×3−C1ROT⋅C1tO1]
在本例中，我们将计算相对于棋盘物体的两个摄像机位置之间的摄像机位移。第一步是计算两幅图像的摄像机位置：
vector<Point2f> corners1, corners2;
bool found1 = findChessboardCorners(img1, patternSize, corners1);
bool found2 = findChessboardCorners(img2, patternSize, corners2);
if (!found1 || !found2)
{
cout << "Error, cannot find the chessboard corners in both images." << endl;
return;
}
vector<Point3f> objectPoints;
calcChessboardCorners(patternSize, squareSize, objectPoints);
FileStorage fs( samples::findFile( intrinsicsPath ), FileStorage::READ);
Mat cameraMatrix, distCoeffs;
fs["camera_matrix"] >> cameraMatrix;
fs["distortion_coefficients"] >> distCoeffs;
Mat rvec1, tvec1;
solvePnP(objectPoints, corners1, cameraMatrix, distCoeffs, rvec1, tvec1);
Mat rvec2, tvec2;
solvePnP(objectPoints, corners2, cameraMatrix, distCoeffs, rvec2, tvec2);
使用上述公式可以根据摄像机姿势计算出摄像机位移：
void computeC2MC1（const Mat &R1，const Mat &tvec1，const Mat &R2，const Mat &tvec2、
Mat &R_1to2, Mat &tvec_1to2)
{
//c2Mc1 = c2Mo * oMc1 = c2Mo * c1Mo.inv()
R_1to2 = R2 * R1.t()；
tvec_1to2 = R2 * (-R1.t()*tvec1) + tvec2；
}
根据摄像机位移计算出的与特定平面相关的同轴度为
在该图中，n 是平面的法向量，d 是摄像机画面与平面之间沿平面法线的距离。根据摄像机位移计算同轴度的方程为
2 H 1 = 2 R 1 − 2 t 1 ⋅ 1 n T 2 d ^2H_1=^2R_1-\frac{^2t_1\cdot ^1n^T}{^2d} 2H1=2R1−2d2t1⋅1nT
其中，2H1 是将第一个相机帧中的点映射到第二个相机帧中相应点的同轴矩阵，
2
R
1
=
C
2
R
O
⋅
C
1
R
O
T
^2R_1=^{C2}R_O\cdot ^{C1}R_O^T
2R1=C2RO⋅C1ROT
是表示两个相机帧之间旋转的旋转矩阵，
2
t
1
=
C
2
R
O
⋅
(
−
C
1
R
O
T
⋅
C
1
t
O
)
+
C
2
t
O
^2t_1=^{C2}R_O\cdot \left( -^{C1}R_{O}^{T}\cdot ^{C1}t_{O}^{} \right) +^{C2}t_O
2t1=C2RO⋅(−C1ROT⋅C1tO)+C2tO
是两个相机帧之间的平移向量。
这里的法线矢量 n 是在摄像机帧 1 中表示的平面法线，可以通过 2 个矢量的乘积（使用平面上的 3 个非对齐点）来计算，或者在我们的例子中直接通过以下方法计算：
Mat normal = (Mat_<double>(3,1) << 0, 0, 1)；
Mat normal1 = R1*normal；
距离 d 可以用平面法线与平面上一点的点积计算，也可以通过计算平面方程和使用 D 系数来计算：
Mat origin(3, 1, CV_64F, Scalar(0))；
Mat origin1 = R1*origin + tvec1；
double d_inv1 = 1.0 / normal1.dot(origin1)；
投影同形矩阵 G 可以通过欧几里得同形矩阵 H 和本征矩阵 K 计算得出（见 [168]），这里假设两个平面视图之间的摄像机相同：
G = γ KH K − 1 \textbf{G} = \gamma \textbf{K} \textbf{H} \textbf{K}^{-1} G=γKHK−1
Mat computeHomography(const Mat &R_1to2, const Mat &tvec_1to2, const double d_inv, const Mat &normal)
{
Mat homography = R_1to2 + d_inv * tvec_1to2*normal.t();
return homography;
}
在我们的例子中，棋盘的 Z 轴位于物体内部，而在 homography 图中，Z 轴位于物体外部。这只是符号问题：
2 H 1 = 2 R 1 + 2 t 1 ⋅ 1 n T 1 d ^2H_1=^2R_1+\frac{^2t_1\cdot ^1n^T}{^1d} 2H1=2R1+1d2t1⋅1nT
Mat homography_euclidean = computeHomography(R_1to2, t_1to2, d_inv1, normal1)；
Mat homography = cameraMatrix * homography_euclidean * cameraMatrix.inv()；
homography /= homography.at<double>(2,2)；
homography_euclidean /= homography_euclidean.at<double>(2,2)；
现在，我们将比较根据摄像机位移计算出的投影同轴度和使用 cv::findHomography 估算出的投影同轴度
findHomography H:
[0.32903393332201, -1.244138808862929, 536.4769088231476;
0.6969763913334046, -0.08935909072571542, -80.34068504082403;
0.00040511729592961, -0.001079740100565013, 0.9999999999999999]
homography from camera displacement:
[0.4160569997384721, -1.306889006892538, 553.7055461075881;
0.7917584252773352, -0.06341244158456338, -108.2770029401219;
0.0005926357240956578, -0.001020651672127799, 1]
同构矩阵相似。如果我们使用这两个同构矩阵来比较图像 1 的扭曲效果，我们可以得出以下结果
从视觉上看，根据摄像机位移计算的同轴度和使用 cv::findHomography 函数估计的同轴度所得到的图像很难区分。
练习
本演示将向您展示如何根据两个摄像机姿势计算同色变换。尝试执行相同的操作，但这次要计算 N 个同轴度。不是计算一个同轴度来直接将源图像扭曲到所需的摄像机视角，而是执行 N 个扭曲操作来查看不同变换的运行情况。
您应该会得到与下面类似的结果：
演示 4：分解同构矩阵
OpenCV 3 包含一个函数 cv::decomposeHomographyMat，可将同色矩阵分解为一组旋转、平移和平面法线。首先，我们将分解根据摄像机位移计算出的同轴矩阵：
Mat homography_euclidean = computeHomography(R_1to2, t_1to2, d_inv1, normal1)；
Mat homography = cameraMatrix * homography_euclidean * cameraMatrix.inv()；
homography /= homography.at<double>(2,2)；
homography_euclidean /= homography_euclidean.at<double>(2,2)；
cv::decomposeHomographyMat 的结果如下：
vector<Mat> Rs_decomp, ts_decomp, normals_decomp;
int solutions = decomposeHomographyMat(homography, cameraMatrix, Rs_decomp, ts_decomp, normals_decomp);
cout << "Decompose homography matrix computed from the camera displacement:" << endl << endl;
for (int i = 0; i < solutions; i++)
{
double factor_d1 = 1.0 / d_inv1;
Mat rvec_decomp;
Rodrigues(Rs_decomp[i], rvec_decomp);
cout << "Solution " << i << ":" << endl;
cout << "rvec from homography decomposition: " << rvec_decomp.t() << endl;
cout << "rvec from camera displacement: " << rvec_1to2.t() << endl;
cout << "tvec from homography decomposition: " << ts_decomp[i].t() << " and scaled by d: " << factor_d1 * ts_decomp[i].t() << endl;
cout << "tvec from camera displacement: " << t_1to2.t() << endl;
cout << "plane normal from homography decomposition: " << normals_decomp[i].t() << endl;
cout << "plane normal at camera 1 pose: " << normal1.t() << endl << endl;
}
Solution 0:
rvec from homography decomposition: [-0.0919829920641369, -0.5372581036567992, 1.310868863540717]
rvec from camera displacement: [-0.09198299206413783, -0.5372581036567995, 1.310868863540717]
tvec from homography decomposition: [-0.7747961019053186, -0.02751124463434032, -0.6791980037590677] and scaled by d: [-0.1578091561210742, -0.005603443652993778, -0.1383378976078466]
tvec from camera displacement: [0.1578091561210745, 0.005603443652993617, 0.1383378976078466]
plane normal from homography decomposition: [-0.1973513139420648, 0.6283451996579074, -0.7524857267431757]
plane normal at camera 1 pose: [0.1973513139420654, -0.6283451996579068, 0.752485726743176]
Solution 1:
rvec from homography decomposition: [-0.0919829920641369, -0.5372581036567992, 1.310868863540717]
rvec from camera displacement: [-0.09198299206413783, -0.5372581036567995, 1.310868863540717]
tvec from homography decomposition: [0.7747961019053186, 0.02751124463434032, 0.6791980037590677] and scaled by d: [0.1578091561210742, 0.005603443652993778, 0.1383378976078466]
tvec from camera displacement: [0.1578091561210745, 0.005603443652993617, 0.1383378976078466]
plane normal from homography decomposition: [0.1973513139420648, -0.6283451996579074, 0.7524857267431757]
plane normal at camera 1 pose: [0.1973513139420654, -0.6283451996579068, 0.752485726743176]
Solution 2:
rvec from homography decomposition: [0.1053487907109967, -0.1561929144786397, 1.401356552358475]
rvec from camera displacement: [-0.09198299206413783, -0.5372581036567995, 1.310868863540717]
tvec from homography decomposition: [-0.4666552552894618, 0.1050032934770042, -0.913007654671646] and scaled by d: [-0.0950475510338766, 0.02138689274867372, -0.1859598508065552]
tvec from camera displacement: [0.1578091561210745, 0.005603443652993617, 0.1383378976078466]
plane normal from homography decomposition: [-0.3131715472900788, 0.8421206145721947, -0.4390403768225507]
plane normal at camera 1 pose: [0.1973513139420654, -0.6283451996579068, 0.752485726743176]
Solution 3:
rvec from homography decomposition: [0.1053487907109967, -0.1561929144786397, 1.401356552358475]
rvec from camera displacement: [-0.09198299206413783, -0.5372581036567995, 1.310868863540717]
tvec from homography decomposition: [0.4666552552894618, -0.1050032934770042, 0.913007654671646] and scaled by d: [0.0950475510338766, -0.02138689274867372, 0.1859598508065552]
tvec from camera displacement: [0.1578091561210745, 0.005603443652993617, 0.1383378976078466]
plane normal from homography decomposition: [0.3131715472900788, -0.8421206145721947, 0.4390403768225507]
plane normal at camera 1 pose: [0.1973513139420654, -0.6283451996579068, 0.752485726743176]
同构矩阵分解的结果只能恢复到与距离 d 相对应的比例因子，因为法线是单位长度。正如您所看到的，有一种解法与计算出的摄像机位移几乎完全吻合。正如文档中所述
如果通过应用正深度约束（所有点必须位于摄像机前方）获得了点对应关系，则至少有两个解可以进一步失效。
由于分解的结果是摄像机位移，如果我们有初始摄像机姿态 c1Mo，我们可以计算当前摄像机姿态
C
2
M
O
=
C
2
M
C
1
⋅
C
1
M
O
^{C2}M_O=^{C2}M_{C1}\cdot ^{C1}M_O
C2MO=C2MC1⋅C1MO
并测试属于平面的三维物体点是否投射到摄像机前方。另一种解决方案是，如果我们知道摄像机 1 姿态下的平面法线，则保留法线最接近的解决方案。
同理，使用 cv::findHomography 估算同构矩阵。
Solution 0:
rvec from homography decomposition: [0.1552207729599141, -0.152132696119647, 1.323678695078694]
rvec from camera displacement: [-0.09198299206413783, -0.5372581036567995, 1.310868863540717]
tvec from homography decomposition: [-0.4482361704818117, 0.02485247635491922, -1.034409687207331] and scaled by d: [-0.09129598307571339, 0.005061910238634657, -0.2106868109173855]
tvec from camera displacement: [0.1578091561210745, 0.005603443652993617, 0.1383378976078466]
plane normal from homography decomposition: [-0.1384902722707529, 0.9063331452766947, -0.3992250922214516]
plane normal at camera 1 pose: [0.1973513139420654, -0.6283451996579068, 0.752485726743176]
Solution 1:
rvec from homography decomposition: [0.1552207729599141, -0.152132696119647, 1.323678695078694]
rvec from camera displacement: [-0.09198299206413783, -0.5372581036567995, 1.310868863540717]
tvec from homography decomposition: [0.4482361704818117, -0.02485247635491922, 1.034409687207331] and scaled by d: [0.09129598307571339, -0.005061910238634657, 0.2106868109173855]
tvec from camera displacement: [0.1578091561210745, 0.005603443652993617, 0.1383378976078466]
plane normal from homography decomposition: [0.1384902722707529, -0.9063331452766947, 0.3992250922214516]
plane normal at camera 1 pose: [0.1973513139420654, -0.6283451996579068, 0.752485726743176]
Solution 2:
rvec from homography decomposition: [-0.2886605671759886, -0.521049903923871, 1.381242030882511]
rvec from camera displacement: [-0.09198299206413783, -0.5372581036567995, 1.310868863540717]
tvec from homography decomposition: [-0.8705961357284295, 0.1353018038908477, -0.7037702049789747] and scaled by d: [-0.177321544550518, 0.02755804196893467, -0.1433427218822783]
tvec from camera displacement: [0.1578091561210745, 0.005603443652993617, 0.1383378976078466]
plane normal from homography decomposition: [-0.2284582117722427, 0.6009247303964522, -0.7659610393954643]
plane normal at camera 1 pose: [0.1973513139420654, -0.6283451996579068, 0.752485726743176]
Solution 3:
rvec from homography decomposition: [-0.2886605671759886, -0.521049903923871, 1.381242030882511]
rvec from camera displacement: [-0.09198299206413783, -0.5372581036567995, 1.310868863540717]
tvec from homography decomposition: [0.8705961357284295, -0.1353018038908477, 0.7037702049789747] and scaled by d: [0.177321544550518, -0.02755804196893467, 0.1433427218822783]
tvec from camera displacement: [0.1578091561210745, 0.005603443652993617, 0.1383378976078466]
plane normal from homography decomposition: [0.2284582117722427, -0.6009247303964522, 0.7659610393954643]
plane normal at camera 1 pose: [0.1973513139420654, -0.6283451996579068, 0.752485726743176]
同样，也有与计算出的摄像机位移相匹配的解决方案。
演示 5：旋转相机的基本全景拼接
注释
本示例旨在说明基于相机纯粹旋转运动的图像拼接概念，不应被用于拼接全景图像。 拼接模块 提供了完整的图像拼接流程。
同轴变换仅适用于平面结构。但在旋转摄像机的情况下（围绕摄像机投影轴的纯旋转，无平移），可以考虑任意世界（见 前文 ）。
然后就可以利用旋转变换和摄像机固有参数计算出同轴度（例如，见 10）：
s
[
x
′
y
′
1
]
=
K
R
K
−
1
[
x
y
1
]
s \begin{bmatrix} x^{'} \\ y^{'} \\ 1 \end{bmatrix} = \bf{K} \hspace{0.1em} \bf{R} \hspace{0.1em} \bf{K}^{-1} \begin{bmatrix} x \\ y \\ 1 \end{bmatrix}
s
x′y′1
=KRK−1
xy1
为了说明这一点，我们使用 Blender（一款免费开源的 3D 计算机制图软件）生成了两个摄像机视图，它们之间只进行了旋转变换。更多关于如何获取摄像机内在参数和相对于世界的 3x4 外在矩阵的信息，请参见 Blender 11（需要额外的变换来获取摄像机和物体帧之间的变换）。
下图显示了苏珊娜模型生成的两个视图，其中只有旋转变换：
利用已知的相关摄像机姿势和固有参数，可以计算出两个视图之间的相对旋转：
C++
Mat R1 = c1Mo(Range(0,3), Range(0,3));
Mat R2 = c2Mo(Range(0,3), Range(0,3));
//c1Mo * oMc2
Mat R_2to1 = R1*R2.t();
Java
Range rowRange = new Range(0,3);
Range colRange = new Range(0,3);
//c1Mo * oMc2
Mat R1 = new Mat(c1Mo, rowRange, colRange);
Mat R2 = new Mat(c2Mo, rowRange, colRange);
Mat R_2to1 = new Mat();
Core.gemm(R1, R2.t(), 1, new Mat(), 0, R_2to1 );
Python
R1 = c1Mo[0:3, 0:3]
R2 = c2Mo[0:3, 0:3]
R2 = R2.transpose()
R_2to1 = np.dot(R1,R2)
在这里，第二幅图像将相对于第一幅图像进行拼接。同轴度可以用上面的公式计算：
C++
Mat H = cameraMatrix * R_2to1 * cameraMatrix.inv();
H /= H.at<double>(2,2);
cout << "H:\n" << H << endl;
Java
Mat tmp = new Mat(), H = new Mat();
Core.gemm(cameraMatrix, R_2to1, 1, new Mat(), 0, tmp);
Core.gemm(tmp, cameraMatrix.inv(), 1, new Mat(), 0, H);
Scalar s = new Scalar(H.get(2, 2)[0]);
Core.divide(H, s, H);
System.out.println(H.dump());
Python
H = cameraMatrix.dot(R_2to1).dot(np.linalg.inv(cameraMatrix))
H = H / H[2][2]
只需使用
C++
Mat img_stitch;
warpPerspective(img2, img_stitch, H, Size(img2.cols*2, img2.rows));
Mat half = img_stitch(Rect(0, 0, img1.cols, img1.rows));
img1.copyTo(half);
Java
Mat img_stitch = new Mat();
Imgproc.warpPerspective(img2, img_stitch, H, new Size(img2.cols()*2, img2.rows()) );
Mat half = new Mat();
half = new Mat(img_stitch, new Rect(0, 0, img1.cols(), img1.rows()));
img1.copyTo(half);
Python
img_stitch = cv.warpPerspective(img2, H, (img2.shape[1]*2, img2.shape[0]))
img_stitch[0:img1.shape[0], 0:img1.shape[1]] = img1
得到的图像是
其他参考资料
-
- 第 16 讲：平面同构，罗伯特-柯林斯
-
- 二维投影变换 (同构), Christiano Gava, Gabriele Bleser
-
- 《计算机视觉： 算法与应用》，Richard Szeliski
-
- 《用于视觉跟踪和平面标记的逐步摄像机姿态估计》，Richard Szeliski
-
- 奇异值分解个人访谈，马坦-加维什（Matan Gavish
-
- 卡布希算法，最优旋转矩阵的计算
-
- 同构，格哈德-罗斯博士
扫一扫
1212
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
