如何将图片中的一个任意四边形区域的图像转化为矩形【附源码】

前段时间导师给了一个任务，任务中包含一个功能，将一幅图像中的任意四边形区域映射为矩形区域。比如，我们有时拍照片，因为角度的问题，本来是矩形的广告牌，在照片中可能表现为一个任意的四边形。事实上，在现实中本来是矩形的物体，在照片中往往呈现为任意四边形。将任意四边形的图片转化为矩形的应用现在在我们身边也时常见到：安卓端的adobe reader最近新推出的扫描功能，因为拍照片时不可能完全正对我们的文档，就包含了该功能；同期，WPS office 也推出了扫描功能；除此之外，当我们扫描二维码时，也可以利用该变换，来实现不正对图片进行扫码。该映射学术名称为中心投影变换，Perspective Mapping。

先来看一下一个示例：（示例只是为了说明情况，也没有找到很好的源图，只是为了给大家一个直观的理解，所以并没有实现精确的变换，显示效果不是很好，大家见谅）

-----------------------------------------------------分割线-----------------------------------------------------------

这是源图，我们的目的是把里边的显示屏显示的部分提取出来，并且映射为矩形

图1. 源图

这是处理后的效果图

图2. 效果图

看了这两张图，大家应该知道算法的目的了。

-----------------------------------------------------分割线-----------------------------------------------------------

原文地址：https://www.geometrictools.com/Documentation/PerspectiveMappings.pdf

大家可以看看原文，自己按照里边的内容在草稿纸上推，很快就能推出来。但是里边是全中文的，且有些术语可能有些难理解，所以我在这里对该文章中用的中心投影变换简单解释一下。

先直接上原理图，大家才好有一个直观的认识

图3. 中心投影变换原理图  (a) 中心投影变换图示  (b)任意四边形  (c)矩形

假设E点是发光源，任意四边形 q00 q01 q11 q10 可以在一个平面上投影为 r00 r01 r11 r10 。如图3(a)所示。q00投影为r00，q01投影为r01，q11投影为r11，q10投影为r10。利用线性代数的知识，只要知道了图3(a)中的各个顶点的坐标，那么可以很简单地将任意四边形中的q的坐标映射到矩形中的映射点r的坐标，我们假设这些点的坐标都已知，接下俩我们需要找到q到r或者r到q的映射关系。

注：

1） q是Quadrilaterals（四边形）的第一个字母，r是rectangle（矩形）的第一个字母；

2） 为了简化问题，我们将q00和r00重合为一点，设为o点；

3） 其实我们可以证明，一个任意四边形在图3(a)的投影方式下一定可以投影为矩形，后面的算法分析过程可以证明这一点，所以这里不讨论这个问题。事实上还可以证明，最终的变换与E点的位置无关。

问题转化为将图3(b)中的任意四边形变换为图3(c)中的矩形。(c)中的矩形是很容易“放”到(a)中去的，至少这个过程能够比较直观地想象出来。但是将(b)中的任意四边形“放"到(a)中，这个过程可能不容易想象。那么我们用线性代数或者高中学的空间几何知识来模型化这个过程。

事实上，我们只需要(b)中的任意四边形Q00 Q01 Q11 Q10映射到(a)中的任意四边形q00 q01q11 q10；将(c)中的矩形R00 R01 R11 R10映射到(a)中的矩形r00 r01 r11 r10，然后根据r到q的映射关系，就能解决该问题了。

进下来进入公式推导环节。为了和原文一致，建立以下坐标系

注：将矩形的四个顶点定义为o, (1 0 0), (0 1 0), (1 1 0)不影响结果，虽然它已经变成正方形了。原因是，我们关心的任意的点r的坐标是通过向量or10和向量or01的线性组合的系数，即

r = x0r10 + x1r01 (1)

所以无论坐标系怎么设置，x0和x1的值并不会变，总是在（0~1  0~1  0）区间内，所以先不要纠结r10和r01的坐标问题，直接往下继续公式推导。

注：r表示or向量，r10表示or10向量，下文也是同样的标注方式，以加粗的方式表示向量。

为了让q00(o点） q01 q11 q10四点共面，最好的方法就是将向量oq11 表示为向量 oq01 和向量 oq10的线性组合，记为

q11 = a0q10 + a1q01(2)

只要任意四边形给定，a0和a1两个值可以直接求出，因此可以看作是常量。

同样，任意四边形中的任意一点q可以表示为

q = y0q10 + y1q01(3)

因为点E，q，r三点共线，所以r点的坐标可以表示为

r = E + t(q-E)(4)

其中t为乘法系数

因此可以得到一下四个方程

r00 = (0,0,0) = E + t00(q00 -E)(5)

r10 = (1,0,0) = E + t10(q10 - E)(6)

r01 = (0,1,0) = E + t01(q01 - E)(7)

r11 = (1,1,0) = E + t11(q11 - E)(8)

易知，t00 = 1

假设平面oq01q11q10存在法向量N，则上面四个方程左右边同时点乘N = (n0,n1,n2)，可得

(1-t10)N·E = n0, (1-t01)N·E = n1, (1-t11)N·E = n0+n1;

将前两个等式等号左边部分相加，得到第三个公式左边部分，即得到

t11 = t10 + t01 -1 (9)

结合公式(9)与(2),(5),(6),(7),(8)可以解出t00,t10,t01,t11

t00 = 1, t10 = a0 / (a0+a1-1), t01 = a1 / (a0+a1-1), t11 = 1 / (a0+a1-1)(10)

注：a0和a1理解为常数，因为一旦任意四边形确定，a0和a1可以马上算出来，体现在matlab中就是一个矩阵左除运算。

根据公式

(1) r = x0r10 + x1r01

(3) q = y0q10 + y1q01

(4) r = E + t(q-E)

将(1)(3)带入(4)，并将r表示为(x0,x1,0)，可得

(x0, x1, 0) = E + t(q - E) =E+ t(y0q10 + y1q01 - E)(11)

根据公式(6),(7)以及已经计算出来的r00,t10,t01,t11值，可以将q10和q01用E和(1,0,0)和(0,1,0)表示出来，然后带入公式(11)，可得

考虑到E和(1,0,0)和(0,1,0)三个向量线性无关，而他们的线性组合等于0向量，所以三个系数都为0

于是解得

该公式的物理意义我解释一下：

1）我们知道了矩形区域内的一个像素点的坐标，假设为(height,width) = (100,100)，同时假设我们想要输出的矩形图像的大小为200*300(height*width)，那么x0为1/2，x1为1/3；也就是说x0和x1为将矩形视为单位正方形的情况下，点的坐标，这也是上面可以用单位长度来表征矩形的原因。

2）对于任意四边形中的一点，y0和y1的值也可以相似地理解，只不过这时候的基向量不再相互垂直，可以在matlab中用左除矩阵的运算方式很快地求出这两个值。

3）该公式可以完成这样一种功能，假如加入我们要求矩形中的一点的像素值，我们可以求出该点对应的x0,x1的值，然后通过该公式找到y0和y1的值，再利用y0和y1和公式(3)，找到任意四边形中一个具体的坐标

注：在任意四边形中找到的具体坐标通常不是整数，所以需要采用一定的插值策略，我的代码中用的最近邻域插值，采用其他高级插值方式本文不讨论。

同时，我们也可以找到反变换，物理意义同上。

以下是matlab源码

%函数功能：中心投影变换。输入源图，源图中的任意四边形的4个点的坐标（左上，右上，左下，右下），以及输出图像的大小（高，宽）
function Imgback = m_PerspectiveTransformation(imgIn,pointLT,pointRT,pointLB,pointRB,outHeitht,outWidth)
    [imgInHeight,imgInWidth,imgInDimension] = size(imgIn);
    %为了中心投影变换，需要将4个点转化为三个向量，具体看参考文献
    vector10 = pointLB - pointLT;
    vector01 = pointRT - pointLT;
    vector11 = pointRB - pointLT;
    %把vector11表示成vector10和vector01的线性组合，以使三个向量共面
    A = [vector10' , vector01'];
    B = vector11' ;
    S = A\B;
    a0 = S(1);
    a1 = S(2);
    
    
    %输出矩形
    Imgback = uint8(zeros(outHeitht,outWidth,imgInDimension));
    
    
    %利用循环操作来对每个像素点赋值
    for heightLoop = 1:outHeitht
        for widthLoop = 1:outWidth
            %以下算法为参考文献中的公式表示
            x0 = heightLoop/outHeitht;
            x1 = widthLoop/outWidth;
            FenMu = a0+a1-1+(1-a1)*x0+(1-a0)*x1;            %分母
            y0 = a0*x0/FenMu;
            y1 = a1*x1/FenMu;
            
            %根据得到的参数找到对应的源图像中的坐标位置，并赋值
            coordInOri = y0*vector10 + y1*vector01 + pointLT;
            heightC = round(coordInOri(1));
            widthC = round(coordInOri(2));
                if (heightC > imgInHeight || heightC <= 0 || widthC >imgInWidth || widthC <=0 )
                    disp(['m_PerspectiveTransformation超出范围' num2str(heightC) num2str(widthC)]);
                    pause();
                    return;
                end
            for dimentionLoop = 1:imgInDimension
                %使用最近邻域插值，使用高级插值方法效果会更好
                Imgback(heightLoop,widthLoop,dimentionLoop) = imgIn(heightC,widthC,dimentionLoop);
            end
        end
    end
    
    
%     figure; imshow(Imgback); title('投影变换的结果');