OpenCV-C++实现单应性矩阵的求解

        单应性矩阵(H)(Homography)约束了同一3D空间点在两个像素平面的2D齐次坐标。

[left[ begin{matrix}   {{u}_{a}}  \   {{v}_{a}}  \   1  \end{matrix} right]=left[ begin{matrix}
   {{H}_{1}} & {{H}_{2}} & {{H}_{3}}  \   {{H}_{4}} & {{H}_{5}} & {{H}_{6}}  \   {{H}_{7}} & {{H}_{8}} & {{H}_{9}}  \end{matrix} right]left[ begin{matrix}   {{u}_{b}}  \   {{v}_{b}}  \   1  \
end{matrix} right]]

        单应性矩阵(H)具有8个自由度，已知A和B两张图像上的四对点，即可列出八个方程来求解出单应性矩阵(H)。

        单应性矩阵(H)(Homography)可以理解为描述物体在世界坐标系和像素坐标系之间的位置映射关系。

        以棋盘格相对位姿估计为例，世界坐标系定在棋盘格上，即z=0。因此世界坐标系和像素坐标系之间的关系如下：

[left[ begin{matrix}   u  \   v  \   1  \end{matrix} right]=sKleft[ begin{matrix}   R & T  \
end{matrix} right]left[ begin{matrix}   x  \   y  \   z  \   1  \end{matrix} right]=sKleft[ begin{matrix}   {{r}_{1}} & {{r}_{2}} & {{r}_{3}} & T  \end{matrix} right]left[ begin{matrix}   x  \   y  \   0  \   1  \end{matrix} right]=sKleft[ begin{matrix}   {{r}_{1}} & {{r}_{2}} & T  \end{matrix} right]left[ begin{matrix}   x  \   y  \   1  \end{matrix} right]]

        因此，(H=sKleft[ begin{matrix}   {{r}_{1}} & {{r}_{2}} & T  \end{matrix} right])即为成像平面与标定板平面之间的单应性矩阵。其中，(K)为相机内参矩阵，(s)是任意尺度的比例（目的是使得单应性定义到该尺度比例）。

        首先抛出结论：单应性矩阵具有8个自由度，需要4对点求解。

        首先，我们假设两张图像中的对应点对齐次坐标为(x',y',1)和(x,y,1)，则有：

[left[ begin{matrix}   {{x}'}  \   {{y}'}  \   1  \end{matrix} right]=left[ begin{matrix}   {{h}_{11}} & {{h}_{12}} & {{h}_{13}}  \   {{h}_{21}} & {{h}_{22}} & {{h}_{23}}  \   {{h}_{31}} & {{h}_{32}} & {{h}_{33}}  \end{matrix} right]left[ begin{matrix}   x  \   y  \   1  \end{matrix} right]]

        矩阵展开后有3个等式，将第3个等式代入前两个等式中可得：

[begin{align}  & {x}'=frac{{{h}_{11}}x+{{h}_{12}}y+{{h}_{13}}}{{{h}_{31}}x+{{h}_{32}}y+{{h}_{33}}} \  & {y}'=frac{{{h}_{21}}x+{{h}_{22}}y+{{h}_{23}}}{{{h}_{31}}x+{{h}_{32}}y+{{h}_{33}}} \ end{align}]

        也就是说，一个点对对应两个等式。

        由于这里使用的是齐次坐标系，也就是说可以进行任意尺度的缩放。比如我们把({{h}_{ij}})乘以任意一个非零常数k并不改变等式结果，所以实际上单应矩阵(H)只有8个自由度。8自由度下(H)计算过程有两种方法。

第一种方法：直接设置 ({{h}_{33}}=1)

第二种方法：将(H)添加约束条件，将(H)矩阵模变为1

        接下来由四个对应点的像素列出8个方程，根据最优化的方法（此处略）可以求解出结果。

        首先随意拍摄一张图片（3072×4096），命名为src.jpg：

        目的是要将书本校正，因此先使用图像处理软件（如全能扫描王，华为自带相册等）获得校正后的图像（2070×2500），命名为dst.jpg：

        使用Photoshop已经预先获取了src.jpg中书本四个角对应的像素点：

左上：(889,1032)

右上：(2471,1079)

左下：(276,2868)

右下：(2829,2979)

注意：x轴方向是水平向右，y轴方向是水平向下，坐标原点是图像左上角第一个像素点。

        另外，也可使用图像分隔等方法检测目标物体的四个特征点。由于不是本篇文章的重点，这里不作详细说明。

        dst.jpg中四个角对应的像素点即图像四个角的像素坐标：

左上：(0,0)

右上：(2069,0)

左下：(0,2499)

右下：(2069,2499)

        主要函数：findHomography与warpPerspective

        C++实现代码如下：

#include 
#include 
#include 

using namespace cv;
using namespace std;

Mat myResize(Mat &srcImag, double width_ratio, double height_ratio);

int main()
{
	Mat im_src,im_dst,im_out,H;
	im_src=imread("../src.jpg");
	im_dst=imread("../dst.jpg");
	Size size=im_dst.size();
	//根据源图片和目标图片相对于的四对点（左上，右上，左下，右下）来计算单应性矩阵
	vector pts_src{Point2f(889,1032),Point2f(2471,1079),Point2f(276,2868),Point2f(2829,2979)};
	vector pts_dst{Point2f(0,0),Point2f(2069,0),Point2f(0,2499),Point2f(2069,2499)};
	//计算单应性矩阵H
	H = findHomography(pts_src,pts_dst);
	cout<< H << endl;
	//根据计算出来的单应性矩阵H，对图像进行透视变换
	warpPerspective(im_src, im_out, H, size);
	//将输出图片按0.3的比例缩小
	im_out=myResize(im_out, 0.3, 0.3);
	imshow("result", im_out);
	
	waitKey();
	destroyAllWindows();
	return 0;
}

//调整图片大小
Mat myResize(Mat &srcImag, double width_ratio, double height_ratio)
{
	Mat src=srcImag;
	float scaleW = width_ratio;
    float scaleH = height_ratio;
    int width = int(src.cols * scaleW);
    int height = int(src.rows * scaleH);
	resize(src, src, Size(width, height));
	return src;
}

        计算出的单应性矩阵如下：

        对src.jpg进行单应性变换的结果如下：

        结果与dst.jpg基本一致，由此证明过程无误。