Homography matrix(单应性矩阵)在广告投放中的实践

全世界只有不到2%的人关注了视觉IMAX​

你真是个特别的人

前言

由于近期在研究相机与投影仪的标定程序时，需要将结构光图片与灰点相机拍摄得到的图片中，找出角点之间的对应性，使用了如下一条代码：

Mat HomoMatrix = findHomography( imgPoints , projPoints,  CV_RANSAC);

如果看过【视觉IMAX】中之前发表过的两篇文章：1）图像处理的仿射变换与透视变换  2）透视变换进阶，应该对findHomography()函数不会感到陌生。

但是前面的两篇文章中，对于单应性矩阵并未做太多讲解。恰巧，今天【视觉IMAX】知识星球中的一个小伙伴也对单应性矩阵进行了发问。

本文我将尽我所能，和大家聊一聊单应性矩阵的基本概念、作用及其的一个简单应用案例，如有不到之处，还望批评指正。

一 单应性矩阵概念

 对于单应性矩阵的概念，此处结合着《Learning OpenCV》，对其进行简单介绍。

在计算机视觉中，平面的单应性被定义为从一个平面到另一个平面的投影映射（小注：术语「单应性」在不同学科上有各种不同的含义。例如，在数学上，它有更通用的意思。在计算机视觉中，对单应性最感兴趣的部分只是其他意义中的一个子集）。

因此，一个二维平面上的点映射到摄像机成像仪上的映射就是平面单应性的例子。如果对点Q到成像仪上的点q的映射使用齐次坐标，这种映射可以用矩阵相乘的方式表示。若有以下定义：

则可以将单应性简单表示为：

这里引入参数s,它是一个任意尺度比例（目的是使得单应性被定义到该尺度比例）。通常习惯放在H的外面。

稍微利用一点几何和矩阵代数的知识，便可以求解这个变换矩阵。最重要的是，H有两部分：用于定位观察的物体平面的物理变换和使用摄像机内参数矩阵的投影。

图1.1：用单应性来描述平面物体的观测：从物体平面到图像平面的映射，同时表征了这两个平面的相对位置和摄像机投影矩阵

物理变换部分是与观测到的图像平面相关的部分旋转R和部分平移t的影响之和。因为使用齐次坐标，我们可以把它们组合到一个单一矩阵中，如下所示：（注：这里W=[R t]是一个3x4矩阵，前三列包含R的9个元素，最后一列由拥有三个元素的向量t组成。）

然后，通过乘以，得到摄像机矩阵M（用来表示影射坐标），即：

这看起来已经完成了。但实际上，我们的关注点不是表征所有空间的坐标，而只是定义我们所寻找的平面的坐标。这需要简化。

考虑到一般性，我们可以选择定义这个物体平面，使得Z=0。这样做的原因是如果把旋转矩阵也分解为3个3x1向量（即R=[r1 r2 r3]），那么其中的一个列向量就不需要了。具体如下：

映射目标点到成像仪的单应性矩阵H可以完全用H=sM[r1 r2 t]表述，其中：

注意，H现在是3x3矩阵。

OpenCV使用上述公式来计算单应性矩阵。它使用同一物体的多个图像来计算每个视场的旋转和平移，同时也计算摄像机的内参数（对所有视场不变）。如我们所讨论的，旋转和平移分别用三个角度和三个偏移量定义，因此对每个视场，有6个未知量。这没有问题，因为一个已知平面物体（如棋盘）能够提供8个方程，即映射一个正方形到四边形可以用4个(x,y)点来描述。每个新的图像帧为计算新的6个未知量提供8个方程，因此若给定足够图像，我们能够计算出未知内参数的任何值（或者更多）。

通过下面的简单方程，单应性矩阵H把源图像平面上的点集位置与目标图像平面（通常为成像仪平面）上的点集位置联系起来：

注意到，我们可以在不知道摄像机内参数的情况下计算H。事实上，OpenCV正是利用从多个视场计算多个单应性矩阵的方法来求解摄像机内参数，如下文所示。

OpenCV提供了一个方便的函数findHomography()，以对应点序列作为输入，返回最佳描述这些对应点的单应性矩阵。我们至少需要4个点来求解H，但是如果有，通常提供更多的点（自然，要保证精确解我们只需要4个对应点。如果提供更多，我们的计算结果便是最小二乘误差意义上的最优解）。使用更多解的好处是尽可能减少噪声和其他不确定因素所带来的干扰。

输入数组srcPoints和dstPoints，即可以是Nx2，也可以是Nx3矩阵。对于前者，点是一像素坐标表示，对于后者，希望是其次坐标。返回的变量是Mat矩阵，其值为3x3矩阵，以保证反向投影误差最小。至于拟合方法，可以选择八点法、ransac、lmeds等，对于噪音比较多的情况，ransac和lmeds会取得比较好的效果。文章开头部分，也是选用的ransac方法。

在单应性矩阵中只有8个独立参数，我们选择归一化，使得=1。但通常的方法是对整个单应性矩阵乘以一个尺度比例。

二 单应性矩阵的作用

由上面的分析可知，单应性矩阵主要用来解决两个问题：

1） 表述真实世界中一个平面与对应它图像的透视变换

2） 通过透视变换实现图像从一种视图变换到另外一种视图

除了概念的理解之外，最重要的，我们更关心它的应用场景，下面就从应用层面简单说说单应性矩阵的应用。

a)用来实现图像拼接时的对齐问题

b)可以用于计算机图形学中的纹理渲染与计算平面阴影

c) 解决拍照时候图像扭曲问题。这可以见文章开篇介绍的两篇文章。

三 一个简单的应用案例

如果我们现在希望在路边的广告牌中，将广告牌中的内容替换为我们自己的宣传内容（当然是虚拟的）。街拍图如下所示：

接下来，我想将我的公众号宣传图投放到红框中，该如何操作呢？

首先，准备好投放的广告素材(本公众号的二维码)，如下图所示。

接下来是最重要的步骤：在原图中获取到对应的广告牌位置，得到广告牌的四个角坐标，通过计算单应性矩阵实现内容替换，最终生成的效果图如下：

附上实现的代码，如下：

#include <opencv.hpp>
#include <iostream>
#include <math.h>
using namespace cv;
using namespace std;
int main(int argc, char** argv) {
    // load images

    Mat src = imread("E:/8-LabPic/imgTest/srcImg.jpg");
    if (!src.data) {
        printf("could not load image...\n");
        return -1;
    }
    // show images
    namedWindow("input image", CV_WINDOW_AUTOSIZE);
    imshow("input image", src);
    Mat replaceImg =imread("E:/
     8-LabPic/imgTest/replaceImg.jpg");
    imshow("adv content", replaceImg);
    // 定义两个平面上四个角坐标
    vector<Point> src_corners(4);
    vector<Point> dst_corners(4);
    // 原图像平面四点坐标
    src_corners[0] = Point(0, 0);
    src_corners[1] = Point(replaceImg.cols, 0);
    src_corners[2] = Point(0, replaceImg.rows);
    src_corners[3] = Point(replaceImg.cols, 
     replaceImg.rows);
    // 目标平面四个角坐标
    dst_corners[0] = Point(218,50);
    dst_corners[1] = Point(557, 30);
    dst_corners[2] = Point(273, 236);
    dst_corners[3] = Point(559, 215);
    // 计算单应性矩阵与透视变换
    Mat h = findHomography(src_corners, dst_corners);
    Mat output_img;
    warpPerspective(replaceImg, output_img, h, 
     src.size());
    // create mask
    Mat m1 = Mat::zeros(replaceImg.size(), CV_8UC1);
    m1 = Scalar(255);
    Mat mask_output;
    warpPerspective(m1, mask_output, h, src.size());
    imshow("Mask Result", mask_output);
    // use mask
    Mat result1;
    add(output_img, output_img, result1, mask_output);
    Mat result2;
    bitwise_not(mask_output, mask_output);
    add(src, result1, result2, mask_output);
    // put them together
    Mat result;
    add(result1, result2, result);
    imshow("Final Result", result);
    imwrite("E:/8-LabPic/imgTest/result.jpg", result);
    waitKey(0);
    return 0;
}

四 跋

顺便交代一下：如果已知图像平面的角点坐标，同时也求出了图像平面与投影平面之间的单应性矩阵，那么，对应点在投影仪中的坐标该如何求呢？

Mat HomoMatrix = findHomography( imgPoints, projPoints,  CV_RANSAC);

Point3d transformedCorner = Point3d( Mat(HomoMatrix * Mat( Point3d(cornerInImage.x, cornerInImage.y, 1.0))));

cornerInProjector=Point2f(transformedCorner.x/ transformedCorner.z,transformedCorner.y/ transformedCorner.z);

文章中对于单应性的理解，讲解不到位的地方，还望多多交流，批评指正。

​