四顶点校正透视变换的线性方程解(转载)
设计给了两张图,一张俯视图a,一张侧视图b,要把a上的点映射到b上,b上的点映射到a上。怎么做?
只有两张图
我们可以通过工具,获取每张图4个顶点的坐标,坐标之间必然存在映射关系,pa = f(pb),用矩阵表示,pa = m*pa。我们要求m。m用3*3的矩阵。
参考:
透视变换(perspective transformation)用于解决仿射变换(affine transformation)无法改变形状内部的相对位置关系的问题。类似photoshop中的“*变换”功能,或者gimp中的“透视”功能,都可以用透视变换矩阵来实现。
现在给定2个四边形:poly1={{x1, y1}, {x2, y2}, {x3, y3}, {x4, y4}}、poly2={{u1, v1}, {u2, v2}, {u3, v3}, {u4, v4}},求做一个透视变换matrix,满足poly1的点能够变形(warp)到poly2中的点。
透视变换矩阵的形式为:
显然,我们需要求解的只是矩阵中8个未知量;已知量包括四边形的4个顶点 x 2个坐标(x, y),一共是8个方程,数量刚刚好。对于每个顶点坐标,poly1中的点总是通过下面的方程转换到poly2中:
问题的关键在于,上面的公式需要求解非线性方程。因此,在解决homography的问题上,常常用最小二乘法(least squares)求得参数估计。
其实稍加分析就不难发现,上面那个使用了分式的所谓“非线性方程”,实际上可以变形为:
它居然变成了线性方程!看来一切试图用newton迭代的算法都是没有必要的,因为我们可以用最简单的线性代数来解决:
下面是我用matlab计算的代码,mat就是我们要求的矩阵,顺便把mat的逆矩阵mat2也求出来,逆变换也有了。
没有matlab也没关系,运算不难,矩阵运算的库也不难找。
u1 = 0; v1 = 0; u2 = 1000; v2 = 0; u3 = 1000; v3 = 1000; u4 = 0; v4 = 1000; x1 = -375; y1 = 510; x2 = 254; y2 = 152; x3 = 1134; y3 = 304; x4 = 782; y4 = 828; a = [ x1, y1, 1, 0, 0, 0, -x1*u1, -y1*u1; 0, 0, 0, x1, y1, 1, -x1*v1, -y1*v1; x2, y2, 1, 0, 0, 0, -x2*u2, -y2*u2; 0, 0, 0, x2, y2, 1, -x2*v2, -y2*v2; x3, y3, 1, 0, 0, 0, -x3*u3, -y3*u3; 0, 0, 0, x3, y3, 1, -x3*v3, -y3*v3; x4, y4, 1, 0, 0, 0, -x4*u4, -y4*u4; 0, 0, 0, x4, y4, 1, -x4*v4, -y4*v4 ]; u = [u1,v1,u2,v2,u3,v3,u4,v4]'; m = a\u; mat = [m(1),m(2),m(3);m(4),m(5),m(6);m(7),m(8),1]; mat2 = inv(mat);
最后求出来结果是
mat = 0.622681800333835 -2.26554353140329 1388.93287614087 1.22602742955981 2.15410964579083 -638.835633268398 1.38327614714492e-06 0.00133144792441439 1 mat2 = 0.434597104229077 0.59516884730679 -223.411138468529 -0.177460145184904 0.0897866610409582 303.8391483172 0.000235677774164818 -0.000120369546353085 0.595763035916078
可以用mat*[x1,y1,1]来验证结果的正确性。
-------------------------------------------------------------------------------------
c#求解方法(应用math.net库)
-------------------------------------------------------------------------------------
double u1 = penvelope.minx; double v1 = penvelope.miny; double u2 = penvelope.maxx; double v2 = penvelope.miny; double u3 = penvelope.maxx; double v3 = penvelope.maxy; double u4 = penvelope.minx; double v4 = penvelope.maxy; double x1 = results[0][0]; double y1 = results[0][1]; double x2 = results[1][0]; double y2 = results[1][1]; double x3 = results[2][0]; double y3 = results[2][1]; double x4 = results[3][0]; double y4 = results[3][1]; // double[,] x = { { x1, y1, 1, 0, 0, 0, -x1*u1, -y1*u1 }, { 0, 0, 0, x1, y1, 1, -x1 * v1, -y1 * v1}, { x2, y2, 1, 0, 0, 0, -x2 * u2, -y2 * u2}, { 0, 0, 0, x2, y2, 1, -x2 * v2, -y2 * v2}, { x3, y3, 1, 0, 0, 0, -x3 * u3, -y3 * u3}, { 0, 0, 0, x3, y3, 1, -x3 * v3, -y3 * v3}, { x4, y4, 1, 0, 0, 0, -x4 * u4, -y4 * u4}, { 0, 0, 0, x4, y4, 1, -x4 * v4, -y4 * v4} }; var matrixx = densematrix.ofarray(x); var vectorb = new densevector(new[] { u1, v1, u2, v2, u3, v3, u4, v4 }); var resultx = new densevector(8); matrixx.gramschmidt().solve(vectorb, resultx); //验证 //console.writeline(x1 + ", " + y1 + ", " + x2 + ", " + y2 + ", " + x3 + ", " + y3 + ", " + x4 + ", " + y4); //console.writeline(u1 + ", " + v1 + ", " + u2 + ", " + v2 + ", " + u3 + ", " + v3 + ", " + u4 + ", " + v4); //console.writeline(resultx); //var u = (resultx[0] * x1 + resultx[1] * y1 + resultx[2]) / (resultx[6] * x1 + resultx[7] * y1 + 1); //console.writeline(u1 + " " + u); double[] result = new double[] { resultx[0], resultx[1], resultx[2], resultx[3], resultx[4], resultx[5], resultx[6], resultx[7], 1 };
推荐阅读