在MATLAB中从4个点计算2D齐次透视变换矩阵

rfbsl7qr  于 8个月前  发布在  Matlab
关注(0)|答案(3)|浏览(104)

我已经得到了在2D中形成矩形的4个点的坐标,以及应用透视变换后的坐标。

透视变换在齐次坐标中计算,并由3x3矩阵M定义。如果矩阵是未知的,我如何从给定的点计算它?
一个点的计算方法是:

| M11 M12 M13 |   | P1.x |   | w*P1'.x |
| M21 M22 M23 | * | P1.y | = | w*P1'.y |
| M31 M32 M33 |   | 1    |   | w*1     |

为了同时计算所有的点,我把它们一起写在一个矩阵A中,类似地,把变换后的点写在一个矩阵B中:

| P1.x P2.x P3.x P4.x |
A = | P1.y P2.y P3.y P4.y |
    | 1    1    1    1    |

所以方程是M*A=B,在MATLAB中可以通过M = B/AM = (A'\B')'求解M
但没那么容易。我知道变换后点的坐标,但我不知道确切的B,因为有因子w,在齐次变换后不需要1。因为在齐次坐标系中,向量的每一个倍数都是同一个点,我不知道我会得到哪个倍数。
为了考虑这些未知因素,我将方程写为M*A=B*W,其中W是一个对角矩阵,B中对角线上的每个点的因子为w1.w4。所以AB现在是完全已知的,我必须解出MW的方程。
如果我可以把方程重新排列成x*A=BA*x=B的形式,其中x将是类似于M*W的东西,我就可以解出它,知道M*W的解可能已经足够了。然而,尽管尝试了所有可能的重新安排,我还是没能做到。直到我意识到封装(M*W)是不可能的,因为一个是3x3矩阵,另一个是4x4矩阵。我被困在这里了。
M*A=B*W也没有M的单一解,因为M的每一个倍数都是相同的变换。将其写成线性方程组,可以简单地固定M中的一个条目以获得单个解。此外,可能有一些输入根本没有M的解决方案,但我们现在不要担心这个问题。
我实际上试图实现的是某种矢量图形编辑程序,用户可以拖动形状边界框的角来转换它,同时在内部计算转换矩阵。
实际上,我需要用JavaScript来实现,但如果我不能用MATLAB来解决这个问题,我就完全被卡住了。

aurhwmvo

aurhwmvo1#

OpenCV有一个很好的函数来完成这个任务,叫做getPerspectiveTransform。这个函数的源代码可以在github上找到,描述如下:

/* Calculates coefficients of perspective transformation
 * which maps (xi,yi) to (ui,vi), (i=1,2,3,4):
 *
 *      c00*xi + c01*yi + c02
 * ui = ---------------------
 *      c20*xi + c21*yi + c22
 *
 *      c10*xi + c11*yi + c12
 * vi = ---------------------
 *      c20*xi + c21*yi + c22
 *
 * Coefficients are calculated by solving linear system:
 * / x0 y0  1  0  0  0 -x0*u0 -y0*u0 \ /c00\ /u0\
 * | x1 y1  1  0  0  0 -x1*u1 -y1*u1 | |c01| |u1|
 * | x2 y2  1  0  0  0 -x2*u2 -y2*u2 | |c02| |u2|
 * | x3 y3  1  0  0  0 -x3*u3 -y3*u3 |.|c10|=|u3|,
 * |  0  0  0 x0 y0  1 -x0*v0 -y0*v0 | |c11| |v0|
 * |  0  0  0 x1 y1  1 -x1*v1 -y1*v1 | |c12| |v1|
 * |  0  0  0 x2 y2  1 -x2*v2 -y2*v2 | |c20| |v2|
 * \  0  0  0 x3 y3  1 -x3*v3 -y3*v3 / \c21/ \v3/
 *
 * where:
 *   cij - matrix coefficients, c22 = 1
 */

这个方程组更小,因为它避免了求解WM33(OpenCV称为c22)。那么它是如何工作的呢?线性系统可以通过以下步骤重新创建:
从一个点的方程开始:

| c00 c01 c02 |   | xi |   | w*ui |
| c10 c11 c12 | * | yi | = | w*vi |
| c20 c21 c22 |   | 1  |   | w*1  |

将其转换为方程组,求解uivi,并消除w。你会得到投影变换的公式:

c00*xi + c01*yi + c02
ui = ---------------------
     c20*xi + c21*yi + c22

     c10*xi + c11*yi + c12
vi = ---------------------
     c20*xi + c21*yi + c22

两边乘以分母:

(c20*xi + c21*yi + c22) * ui = c00*xi + c01*yi + c02
(c20*xi + c21*yi + c22) * vi = c10*xi + c11*yi + c12

分发uivi

c20*xi*ui + c21*yi*ui + c22*ui = c00*xi + c01*yi + c02
c20*xi*vi + c21*yi*vi + c22*vi = c10*xi + c11*yi + c12

假设c22 = 1

c20*xi*ui + c21*yi*ui + ui = c00*xi + c01*yi + c02
c20*xi*vi + c21*yi*vi + vi = c10*xi + c11*yi + c12

收集左手侧的所有cij

c00*xi + c01*yi + c02 - c20*xi*ui - c21*yi*ui = ui
c10*xi + c11*yi + c12 - c20*xi*vi - c21*yi*vi = vi

最后转换成四对点的矩阵形式:

/ x0 y0  1  0  0  0 -x0*u0 -y0*u0 \ /c00\ /u0\
| x1 y1  1  0  0  0 -x1*u1 -y1*u1 | |c01| |u1|
| x2 y2  1  0  0  0 -x2*u2 -y2*u2 | |c02| |u2|
| x3 y3  1  0  0  0 -x3*u3 -y3*u3 |.|c10|=|u3|
|  0  0  0 x0 y0  1 -x0*v0 -y0*v0 | |c11| |v0|
|  0  0  0 x1 y1  1 -x1*v1 -y1*v1 | |c12| |v1|
|  0  0  0 x2 y2  1 -x2*v2 -y2*v2 | |c20| |v2|
\  0  0  0 x3 y3  1 -x3*v3 -y3*v3 / \c21/ \v3/

这是现在的形式Ax=b和解决方案可以得到x = A\b。记住c22 = 1

ocebsuys

ocebsuys2#

这个问题应该很简单那么,如何将M*A=B*W转化为可解形式呢?这只是矩阵乘法,所以我们可以把它写成线性方程组。你知道,就像:M11*A11 + M12*A21 + M13*A31 = B11*W11 + B12*W21 + B13*W31 + B14*W41。每个线性方程组都可以写成Ax=b的形式,或者为了避免与我的问题中已经使用的变量混淆:N*x=y就这样
根据我的问题举个例子:我用已知的MW生成一些输入数据:

M = [
    1 2 3;
    4 5 6;
    7 8 1
];
A = [
    0 0 1 1;
    0 1 0 1;
    1 1 1 1
];
W = [
    4 0 0 0;
    0 3 0 0;
    0 0 2 0;
    0 0 0 1
];
B = M*A*(W^-1);

然后我忘记了MW。意味着我现在有13个变量要解。我把M*A=B*W重写成一个线性方程组,并从那里变成N*x=y的形式。在N中,每一列都有一个变量的因子:

N = [
    A(1,1) A(2,1) A(3,1)      0      0      0      0      0      0 -B(1,1)       0       0       0;
         0      0      0 A(1,1) A(2,1) A(3,1)      0      0      0 -B(2,1)       0       0       0;
         0      0      0      0      0      0 A(1,1) A(2,1) A(3,1) -B(3,1)       0       0       0;
    A(1,2) A(2,2) A(3,2)      0      0      0      0      0      0       0 -B(1,2)       0       0;
         0      0      0 A(1,2) A(2,2) A(3,2)      0      0      0       0 -B(2,2)       0       0;
         0      0      0      0      0      0 A(1,2) A(2,2) A(3,2)       0 -B(3,2)       0       0;
    A(1,3) A(2,3) A(3,3)      0      0      0      0      0      0       0       0 -B(1,3)       0;
         0      0      0 A(1,3) A(2,3) A(3,3)      0      0      0       0       0 -B(2,3)       0;
         0      0      0      0      0      0 A(1,3) A(2,3) A(3,3)       0       0 -B(3,3)       0;
    A(1,4) A(2,4) A(3,4)      0      0      0      0      0      0       0       0       0 -B(1,4);
         0      0      0 A(1,4) A(2,4) A(3,4)      0      0      0       0       0       0 -B(2,4);
         0      0      0      0      0      0 A(1,4) A(2,4) A(3,4)       0       0       0 -B(3,4);
         0      0      0      0      0      0      0      0      1       0       0       0       0
];

y是:

y = [ 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 1 ];

注意N中最后一行描述的方程,根据y,其解为1。这就是我在我的问题中提到的,你必须修复M的一个条目才能得到一个解决方案。(我们可以这样做,因为M的每一个倍数都是相同的变换。)根据这个等式,我说M33应该是1。
我们对x求解:

x = N\y

并获得:

x = [ 1.00000; 2.00000; 3.00000; 4.00000; 5.00000; 6.00000; 7.00000; 8.00000; 1.00000; 4.00000; 3.00000; 2.00000; 1.00000 ]

[ M11, M12, M13, M21, M22, M23, M31, M32, M33, w1, w2, w3, w4 ]的解决方案
M计算完成后就不需要W了。对于一般点(x, y),在求解x'y'时计算对应的w

| M11 M12 M13 |   | x |   | w * x' |
| M21 M22 M23 | * | y | = | w * y' |
| M31 M32 M33 |   | 1 |   | w * 1  |

当在JavaScript中解决这个问题时,我可以使用Numeric JavaScript库,该库具有所需的函数solve来解决Ax=b。

b1zrtrql

b1zrtrql3#

function T = calculatePerspectiveTransformMatrix(srcPts, destPts)
    % Check if the input points are valid
    if size(srcPts, 1) ~= 4 || size(destPts, 1) ~= 4
        error('Both source and destination points must have 4 rows');
    end

    % Construct the coefficient matrix A
    A = zeros(8, 8);
    b = zeros(8, 1);
    for i = 1:4
        A(i*2-1, :) = [srcPts(i, 1), srcPts(i, 2), 1, 0, 0, 0, -srcPts(i, 1)*destPts(i, 1), -srcPts(i, 2)*destPts(i, 1)];
        A(i*2, :) = [0, 0, 0, srcPts(i, 1), srcPts(i, 2), 1, -srcPts(i, 1)*destPts(i, 2), -srcPts(i, 2)*destPts(i, 2)];
        b(i*2-1:i*2) = [destPts(i, 1); destPts(i, 2)];
    end

    % Solve the linear system to obtain the perspective transform matrix
    x = A \ b;
    T = reshape([x; 1], [3, 3])';

    % Normalize the matrix
    T = T / T(3, 3);
end

要使用此函数,可以传入两组对应的点作为输入参数(srcPts和destPts)。每组点应该是一个四行两列的矩阵,其中每行代表一个点的x和y坐标。
例如,你可以像这样调用函数:

srcPts = [10, 20; 30, 40; 50, 60; 70, 80];
destPts = [100, 200; 300, 400; 500, 600; 700, 800];

T = calculatePerspectiveTransformMatrix(srcPts, destPts);

该函数将返回将源点Map到目的地点的3x3透视变换矩阵T。

相关问题