前段时间导师给了一个任务，任务中包含一个功能，将一幅图像中的任意四边形区域映射为矩形区域。比如，我们有时拍照片，因为角度的问题，本来是矩形的广告牌，在照片中可能表现为一个任意的四边形。事实上，在现实中本来是矩形的物体，在照片中往往呈现为任意四边形，因此很多情况下有必要将这些任意四边形变换回矩形。将任意四边形的图片转化为矩形的应用现在在我们身边也时常见到：安卓端的adobe reader最近（2017年）新推出的扫描功能，因为拍照片时不可能完全正对我们的文档，就包含了该功能；同期，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)的投影方式下一定可以投影为矩形，后面的算法分析过程可以证明这一点。事实上还可以证明，r点和q点的坐标变换公式与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的值并不会变，即长方形和正方形原理是一样的。
注：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;
将前两个等式等号左边部分相加，与第三个等式比较，发现等号右边都为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）对于任意四边形中的一个点，y1的值类似地理解，只不过这时候的基向量不再相互垂直，可以在matlab中用左除矩阵的方式很快地求出这两个值。
3）该公式可以完成这样一种功能：假如需要找到矩形中一个点对应的任意四边形中的点，我们可以求出该点对应的x0,x1的值，然后通过该公式找到y0和y1的值，再利用y0和y1和公式(3)，找到任意四边形中一个具体的点的坐标。
注：在任意四边形中找到的具体坐标通常不是整数，所以需要采用一定的插值策略，下面的代码中用的最近邻域插值。采用其他高级插值方式可能效果更好，不过这不是本文的中心，本文不讨论。
同时，我们也可以找到反变换，物理意义同上。
以下是matlab源码
注：这只是核心函数，外围测试函数请参考我下面给出的链接中的资源。
%函数功能：中心投影变换。输入源图，源图中的任意四边形的4个点的坐标（左上，右上，左下，右下），对应于上面示例图中显示器四个角在源图中的四个坐标点，以及输出图像的大小（高，宽）。注意四个点的顺序不是左上，右上，右下，左下。如果有必要，请自行调整顺序。
%例如 imgBack = m_PerspectiveTransformation(imgSrc,[831 1542],[275 3096],[1727 1394],[1601 3151],900,1600);
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('投影变换的结果');
示例程序可以在这里下载（matlab代码）
596
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
