一,理论
略
二,计算
已知从原图位置
(
u
,
v
,
w
)
(u,v,w)
(u,v,w)到目标图像位置
(
x
′
,
y
′
,
z
′
)
(x',y',z')
(x′,y′,z′)的变换矩阵为
A
A
A
A
=
[
a
11
a
12
a
13
a
21
a
22
a
23
a
31
a
32
a
33
]
A=\left[ \begin{array}{ccc} a_{11} & a_{12} & a_{13}\\\\ a_{21} & a_{22} & a_{23}\\\\ a_{31} & a_{32} & a_{33}\\\\ \end{array} \right]
A=⎣⎢⎢⎢⎢⎢⎢⎡a11a21a31a12a22a32a13a23a33⎦⎥⎥⎥⎥⎥⎥⎤
即
(
x
′
,
y
′
,
z
′
)
=
(
u
,
v
,
w
)
∗
A
=
(
a
11
u
+
a
21
v
+
a
31
w
,
a
12
u
+
a
22
v
+
a
32
w
,
a
13
u
+
a
23
v
+
a
33
w
)
(x',y',z')=(u,v,w)*A=(a_{11}u+a_{21}v+a_{31}w,a_{12}u+a_{22}v+a_{32}w,a_{13}u+a_{23}v+a_{33}w)
(x′,y′,z′)=(u,v,w)∗A=(a11u+a21v+a31w,a12u+a22v+a32w,a13u+a23v+a33w)
a
33
=
1
,
w
=
1
a_{33}=1,w=1
a33=1,w=1,从而有
x = x ′ z ′ = a 11 u + a 21 v + a 31 a 13 u + a 23 v + 1 x=\frac{x'}{z'}=\frac{a_{11}u+a_{21}v+a_{31}}{a_{13}u+a_{23}v+1} x=z′x′=a13u+a23v+1a11u+a21v+a31,
y = y ′ z ′ = a 12 u + a 22 v + a 32 a 13 u + a 23 v + 1 y=\frac{y'}{z'}=\frac{a_{12}u+a_{22}v+a_{32}}{a_{13}u+a_{23}v+1} y=z′y′=a13u+a23v+1a12u+a22v+a32.---------(1)
一个点有两个方程,解
A
A
A的8个参数至少需要4个点,设原图的四个点为
(
u
1
,
v
1
)
,
(
u
2
,
v
2
)
,
(
u
3
,
v
3
)
,
(
u
4
,
v
4
)
(u_1,v_1),(u_2,v_2),(u_3,v_3),(u_4,v_4)
(u1,v1),(u2,v2),(u3,v3),(u4,v4),目标图像四个点为
(
x
1
,
y
1
)
,
(
x
2
,
y
2
)
,
(
x
3
,
y
3
)
,
(
x
4
,
y
4
)
(x_1,y_1),(x_2,y_2),(x_3,y_3),(x_4,y_4)
(x1,y1),(x2,y2),(x3,y3),(x4,y4),则有
x
i
=
a
11
u
i
+
a
21
v
i
+
a
31
a
13
u
i
+
a
23
v
i
+
1
x_i=\frac{a_{11}u_i+a_{21}v_i+a_{31}}{a_{13}u_i+a_{23}v_i+1}
xi=a13ui+a23vi+1a11ui+a21vi+a31
y i = a 12 u i + a 22 v i + a 32 a 13 u i + a 23 v i + 1 , i = 1 , 2 , 3 , 4 y_i=\frac{a_{12}u_i+a_{22}v_i+a_{32}}{a_{13}u_i+a_{23}v_i+1},i=1,2,3,4 yi=a13ui+a23vi+1a12ui+a22vi+a32,i=1,2,3,4。展开得到
a
11
u
i
+
a
12
0
−
a
13
u
i
x
i
+
a
21
v
i
+
a
22
0
−
a
23
v
i
x
i
+
a
31
+
a
32
0
=
x
i
a_{11}u_i+a_{12}0-a_{13}u_ix_i+a_{21}v_i+a_{22}0-a_{23}v_ix_i+a_{31}+a_{32}0=x_i
a11ui+a120−a13uixi+a21vi+a220−a23vixi+a31+a320=xi
a
11
0
+
a
12
u
i
−
a
13
u
i
y
i
+
a
21
0
+
a
22
v
i
−
a
23
v
i
y
i
+
a
31
0
+
a
32
=
y
i
,
i
=
1
,
2
,
3
,
4
a_{11}0+a_{12}u_i-a_{13}u_iy_i+a_{21}0+a_{22}v_i-a_{23}v_iy_i+a_{31}0+a_{32}=y_i,i=1,2,3,4
a110+a12ui−a13uiyi+a210+a22vi−a23viyi+a310+a32=yi,i=1,2,3,4。即
[
u
1
0
−
u
1
x
1
v
1
0
−
v
1
x
1
1
0
0
u
1
−
u
1
y
1
0
v
1
−
v
1
y
1
0
1
u
2
0
−
u
2
x
2
v
2
0
−
v
2
x
2
1
0
0
u
2
−
u
2
y
2
0
v
2
−
v
2
y
2
0
1
u
3
0
−
u
3
x
3
v
3
0
−
v
3
x
3
1
0
0
u
3
−
u
3
y
3
0
v
3
−
v
3
y
3
0
1
u
4
0
−
u
4
x
4
v
4
0
−
v
4
x
4
1
0
0
u
4
−
u
4
y
4
0
v
4
−
v
4
y
4
0
1
]
.
[
a
11
a
12
a
13
a
21
a
22
a
23
a
31
a
32
]
=
[
x
1
y
1
x
2
y
2
x
1
y
1
x
2
y
2
]
\left[ \begin{array}{ccc} u_1 &0&-u_1x_1 &v_1&0&-v_1x_1&1&0 \\ 0&u_1&-u_1y_1&0&v_1&-v_1y_1&0&1\\ u_2 &0&-u_2x_2 &v_2&0&-v_2x_2&1&0 \\ 0&u_2&-u_2y_2&0&v_2&-v_2y_2&0&1\\ u_3 &0&-u_3x_3 &v_3&0&-v_3x_3&1&0 \\ 0&u_3&-u_3y_3&0&v_3&-v_3y_3&0&1\\ u_4 &0&-u_4x_4 &v_4&0&-v_4x_4&1&0 \\ 0&u_4&-u_4y_4&0&v_4&-v_4y_4&0&1\\ \end{array} \right].\left[ \begin{array}{ccc} a_{11}\\ a_{12}\\ a_{13}\\ a_{21}\\ a_{22}\\ a_{23}\\ a_{31}\\ a_{32}\\ \end{array} \right]=\left[ \begin{array}{ccc} x_1\\ y_1\\ x_2\\ y_2\\ x_1\\ y_1\\ x_2\\ y_2\\ \end{array} \right]
⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡u10u20u30u400u10u20u30u4−u1x1−u1y1−u2x2−u2y2−u3x3−u3y3−u4x4−u4y4v10v20v30v400v10v20v30v4−v1x1−v1y1−v2x2−v2y2−v3x3−v3y3−v4x4−v4y41010101001010101⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤.⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡a11a12a13a21a22a23a31a32⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤=⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡x1y1x2y2x1y1x2y2⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤
从而有,
[
a
11
a
12
a
13
a
21
a
22
a
23
a
31
a
32
]
=
[
u
1
0
−
u
1
x
1
v
1
0
−
v
1
x
1
1
0
0
u
1
−
u
1
y
1
0
v
1
−
v
1
y
1
0
1
u
2
0
−
u
2
x
2
v
2
0
−
v
2
x
2
1
0
0
u
2
−
u
2
y
2
0
v
2
−
v
2
y
2
0
1
u
3
0
−
u
3
x
3
v
3
0
−
v
3
x
3
1
0
0
u
3
−
u
3
y
3
0
v
3
−
v
3
y
3
0
1
u
4
0
−
u
4
x
4
v
4
0
−
v
4
x
4
1
0
0
u
4
−
u
4
y
4
0
v
4
−
v
4
y
4
0
1
]
−
1
.
[
x
1
y
1
x
2
y
2
x
1
y
1
x
2
y
2
]
\left[ \begin{array}{ccc} a_{11}\\ a_{12}\\ a_{13}\\ a_{21}\\ a_{22}\\ a_{23}\\ a_{31}\\ a_{32}\\ \end{array} \right]=\left[ \begin{array}{ccc} u_1 &0&-u_1x_1 &v_1&0&-v_1x_1&1&0 \\ 0&u_1&-u_1y_1&0&v_1&-v_1y_1&0&1\\ u_2 &0&-u_2x_2 &v_2&0&-v_2x_2&1&0 \\ 0&u_2&-u_2y_2&0&v_2&-v_2y_2&0&1\\ u_3 &0&-u_3x_3 &v_3&0&-v_3x_3&1&0 \\ 0&u_3&-u_3y_3&0&v_3&-v_3y_3&0&1\\ u_4 &0&-u_4x_4 &v_4&0&-v_4x_4&1&0 \\ 0&u_4&-u_4y_4&0&v_4&-v_4y_4&0&1\\ \end{array} \right]^{-1}.\left[ \begin{array}{ccc} x_1\\ y_1\\ x_2\\ y_2\\ x_1\\ y_1\\ x_2\\ y_2\\ \end{array} \right]
⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡a11a12a13a21a22a23a31a32⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤=⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡u10u20u30u400u10u20u30u4−u1x1−u1y1−u2x2−u2y2−u3x3−u3y3−u4x4−u4y4v10v20v30v400v10v20v30v4−v1x1−v1y1−v2x2−v2y2−v3x3−v3y3−v4x4−v4y41010101001010101⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤−1.⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡x1y1x2y2x1y1x2y2⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤
最终得到变换矩阵,不过这个八阶矩阵求逆直接劝退我自己写代码的想法了,后来看到这个思路不错,其思路是两个四边形变换之间用一个小的正方形过度一下,这个正方形的坐标不是0就是1大大简化了计算,后面的代码也是基于这个思路做的。
如上图所示,从甲到到乙直接地变换是
A
A
A,甲到丙的变换是
A
1
A_1
A1,丙到乙的变换是
A
2
A_2
A2。先考虑从正方形丙到四边形乙的变换。其中,正方形四个点的坐标为
1
:
(
0
,
1
)
,
2
:
(
1
,
1
)
,
3
:
(
1
,
0
)
,
4
:
(
0
,
0
)
1:(0,1),2:(1,1),3:(1,0),4:(0,0)
1:(0,1),2:(1,1),3:(1,0),4:(0,0),乙的四个点的坐标为
1
:
(
x
1
,
y
1
)
,
2
:
(
x
2
,
y
2
)
,
3
:
(
x
3
,
y
3
)
,
4
:
(
x
4
,
y
4
)
1:(x_1,y_1),2:(x_2,y_2),3:(x_3,y_3),4:(x_4,y_4)
1:(x1,y1),2:(x2,y2),3:(x3,y3),4:(x4,y4),两边的四个点分别对应,由公式(1),有
x
1
=
a
21
+
a
31
a
23
+
1
x_1=\frac{a_{21} + a_{31}}{a_{23}+1}
x1=a23+1a21+a31
y
1
=
a
22
+
a
32
a
23
+
1
y_1=\frac{a_{22}+a_{32}}{a_{23}+1}
y1=a23+1a22+a32
x
2
=
a
11
+
a
21
+
a
31
a
13
+
a
23
+
1
x_2=\frac{a_{11}+a_{21}+a_{31}}{a_{13}+a_{23}+1}
x2=a13+a23+1a11+a21+a31
y
2
=
a
12
+
a
22
+
a
32
a
13
+
a
23
+
1
y_2=\frac{a_{12}+a_{22}+a_{32}}{a_{13}+a_{23}+1}
y2=a13+a23+1a12+a22+a32
x
3
=
a
11
+
a
31
a
13
+
1
x_3=\frac{a_{11}+a_{31}}{a_{13}+1}
x3=a13+1a11+a31
y
3
=
a
12
+
a
32
a
13
+
1
y_3=\frac{a_{12} + a_{32}}{a_{13}+1}
y3=a13+1a12+a32
x
4
=
a
31
x_4=a_{31}
x4=a31
y
4
=
a
32
y_4=a_{32}
y4=a32
将
x
4
,
y
4
x_4,y_4
x4,y4带入并化简得
a
21
=
a
23
x
1
+
x
1
−
x
4
a_{21}=a_{23}x_1+x_1-x_4
a21=a23x1+x1−x4 ----(2)
a
22
=
a
23
y
1
+
y
1
−
y
4
a_{22}=a_{23}y_1+y_1-y_4
a22=a23y1+y1−y4 ----(3)
a
11
=
a
13
x
3
+
x
3
−
x
4
a_{11}=a_{13}x_3+x_3-x_4
a11=a13x3+x3−x4 ----(4)
a
12
=
a
13
y
3
+
y
3
−
y
4
a_{12}=a_{13}y_3+y_3-y_4
a12=a13y3+y3−y4 ----(5)
a
11
+
a
21
−
a
13
x
2
−
a
23
x
2
=
x
2
−
x
4
a_{11}+a_{21}-a_{13}x_2-a_{23}x_2=x_2-x_4
a11+a21−a13x2−a23x2=x2−x4 ----(6)
a
12
+
a
22
−
a
13
y
2
−
a
23
y
2
=
y
2
−
y
4
a_{12}+a_{22}-a_{13}y_2-a_{23}y_2=y_2-y_4
a12+a22−a13y2−a23y2=y2−y4 ----(7)
将上面(2、3、4、5)带入到(6、7)进行消元得到关于
a
13
,
a
23
a_{13},a_{23}
a13,a23的等式
a
13
(
x
3
−
x
2
)
+
a
23
(
x
1
−
x
2
)
=
x
2
−
x
1
+
x
4
−
x
3
a
13
(
y
3
−
y
2
)
+
a
23
(
y
1
−
y
2
)
=
y
2
−
y
1
+
y
4
−
y
3
a_{13}(x_3-x_2)+a_{23}(x_1-x_2)=x_2-x_1+x_4-x_3\\ a_{13}(y_3-y_2)+a_{23}(y_1-y_2)=y_2-y_1+y_4-y_3
a13(x3−x2)+a23(x1−x2)=x2−x1+x4−x3a13(y3−y2)+a23(y1−y2)=y2−y1+y4−y3
即
[
x
3
−
x
2
x
1
−
x
2
y
3
−
y
2
y
1
−
y
2
]
.
[
a
13
a
23
]
=
[
x
2
−
x
1
−
x
3
+
x
4
y
2
−
y
1
+
y
4
−
y
3
]
\left[ \begin{array}{ccc} x_3-x_2 &x_1-x_2 \\ y_3-y_2&y_1-y_2\\ \end{array} \right].\left[ \begin{array}{ccc} a_{13}\\ a_{23}\\ \end{array} \right]=\left[ \begin{array}{ccc} x_2-x_1-x_3+x_4\\ y_2-y_1+y_4-y_3\\ \end{array} \right]
[x3−x2y3−y2x1−x2y1−y2].[a13a23]=[x2−x1−x3+x4y2−y1+y4−y3]
当最左边矩阵可逆时
[
a
13
a
23
]
=
[
x
3
−
x
2
x
1
−
x
2
y
3
−
y
2
y
1
−
y
2
]
−
1
.
[
x
2
−
x
1
−
x
3
+
x
4
y
2
−
y
1
+
y
4
−
y
3
]
\left[ \begin{array}{ccc} a_{13}\\ a_{23}\\ \end{array} \right]=\left[ \begin{array}{ccc} x_3-x_2 &x_1-x_2 \\ y_3-y_2&y_1-y_2\\ \end{array} \right]^{-1}.\left[ \begin{array}{ccc} x_2-x_1-x_3+x_4\\ y_2-y_1+y_4-y_3\\ \end{array} \right]
[a13a23]=[x3−x2y3−y2x1−x2y1−y2]−1.[x2−x1−x3+x4y2−y1+y4−y3]
左边矩阵不可逆时容易计算得到1,2,3三点共线,将透视矩阵代入可得四点共线,也就是该映射退化成将平面映射成一条线了,自己用的时候基本上用不到这种情况,这里不考虑该情况了,得到
a
13
,
a
23
a_{13},a_{23}
a13,a23后带入(2),(3),(4),(5),可将所有参数计算得出,求出
A
2
A_2
A2后,再求丙到甲的变换矩阵
A
1
−
1
A_1^{-1}
A1−1,最终得到
A
=
A
1
∗
A
2
A=A_1*A_2
A=A1∗A2
三,代码
#include <iostream>
#include <opencv2/opencv.hpp>
typedef struct _M_POINT_ {
int x;
int y;
}M_POINT;
typedef struct _POINTS_ {
M_POINT p1;
M_POINT p2;
M_POINT p3;
M_POINT p4;
}POINTS;
class Persperctive_transform {
public:
Persperctive_transform() {};
~Persperctive_transform() {};
void persperctive_transformation(POINTS p_ori, POINTS p_goal, float T[9])
{
float A1[9];
float A2[9];
float A1_1[9];
transformT(p_ori, A1);
transformT(p_goal, A2);
mat_inverse_3d(A1, A1_1);
Mat_multy(A1_1, A2, T);
}
private:
void transformT(POINTS p, float A[9])
{
//float A[9];
float det = (p.p3.x - p.p2.x) * (p.p1.y - p.p2.y) - (p.p1.x - p.p2.x) * (p.p3.y - p.p2.y);
A[2] = (p.p1.y - p.p2.y) * (p.p2.x - p.p1.x + p.p4.x - p.p3.x) + (p.p2.x - p.p1.x) * (p.p2.y - p.p1.y + p.p4.y - p.p3.y);
A[2] /= det;
A[0] = A[2] * p.p3.x + p.p3.x - p.p4.x;
A[1] = A[2] * p.p3.y + p.p3.y - p.p4.y;
A[5] = (p.p2.y - p.p3.y) * (p.p2.x - p.p1.x + p.p4.x - p.p3.x) + (p.p3.x - p.p2.x) * (p.p2.y - p.p1.y + p.p4.y - p.p3.y);
A[5] /= det;
A[3] = A[5] * p.p1.x + p.p1.x - p.p4.x;
A[4] = A[5] * p.p1.y + p.p1.y - p.p4.y;
A[6] = p.p4.x;
A[7] = p.p4.y;
A[8] = 1;
}int pow(int x, int n)
{
int t = 1;
for (int i = 0; i < n; i++)
{
t *= x;
}
return t;
}
void mat_inverse_3d(float ori[9], float res[9])
{
float det = 0;
float A11, A12, A13, A21, A22, A23, A31, A32, A33;
A11 = pow(-1, 2) * (ori[4] * ori[8] - ori[5] * ori[7]);
A12 = pow(-1, 3) * (ori[3] * ori[8] - ori[5] * ori[6]);
A13 = pow(-1, 4) * (ori[3] * ori[7] - ori[4] * ori[6]);
A21 = pow(-1, 3) * (ori[1] * ori[8] - ori[2] * ori[7]);
A22 = pow(-1, 4) * (ori[0] * ori[8] - ori[2] * ori[6]);
A23 = pow(-1, 5) * (ori[0] * ori[7] - ori[1] * ori[6]);
A31 = pow(-1, 4) * (ori[1] * ori[5] - ori[2] * ori[4]);
A32 = pow(-1, 5) * (ori[0] * ori[5] - ori[2] * ori[3]);
A33 = pow(-1, 6) * (ori[0] * ori[4] - ori[1] * ori[3]);
det += pow(-1, 2) * ori[0] * (ori[4] * ori[8] - ori[7] * ori[5]);
det += pow(-1, 3) * ori[3] * (ori[1] * ori[8] - ori[7] * ori[2]);
det += pow(-1, 4) * ori[6] * (ori[1] * ori[5] - ori[2] * ori[4]);
if (abs(det) < 0.0001)
{
std::cout << "no inverse\n";
return;
}
res[0] = A11 / det;
res[1] = A21 / det;
res[2] = A31 / det;
res[3] = A12 / det;
res[4] = A22 / det;
res[5] = A32 / det;
res[6] = A13 / det;
res[7] = A23 / det;
res[8] = A33 / det;
}
void Mat_multy(float left[9], float right[9], float res[9])
{
for (int i = 0; i < 9; i++)
{
int row = i / 3;
int col = i % 3;
float value = 0;
for (int i = 0; i < 3; i++)
{
value += left[row * 3 + i] * right[i * 3 + col];
}
if (abs(value) < 0.0001) value = 0;
res[i] = value;
}
}
};
int main()
{
Persperctive_transform PT;
POINTS p_ori;
POINTS p_goal;
cv::Mat img=cv::imread("C:\\Users\\user\\Desktop\\tes\\lena.jpg");
cv::Mat res = cv::Mat::zeros(img.rows, img.cols, CV_8UC3);
p_ori.p1.x = 0; p_ori.p1.y = 0;
p_ori.p2.x = 0; p_ori.p2.y = img.cols;
p_ori.p3.x = img.rows; p_ori.p3.y = img.cols;
p_ori.p4.x = img.rows; p_ori.p4.y = 0;
p_goal.p1.x = 150; p_goal.p1.y = 50;
p_goal.p2.x = 0; p_goal.p2.y = img.cols;
p_goal.p3.x = img.rows-100; p_goal.p3.y = img.cols-100;
p_goal.p4.x = img.rows; p_goal.p4.y = 0+0;
float T[9];
PT.persperctive_transformation(p_ori, p_goal, T);
for (int i = 0; i < img.rows; i++)
{
for (int j = 0; j < img.cols; j++)
{
float w = T[2] * i + T[5] * j + T[8];
float row = T[0] * i + T[3] * j + T[6];
float col = T[1] * i + T[4] * j + T[7];
int R = row / w;
int C = col / w;
if (R < 0 || C < 0 || R>img.rows || C>img.cols) continue;
res.at<cv::Vec3b>(R, C) = img.at<cv::Vec3b>(i, j);
}
}
cv::Mat res2 = cv::Mat::zeros(img.rows, img.cols, CV_8UC3);
p_ori.p1.x = 150; p_ori.p1.y = 50;
p_ori.p2.x = 0; p_ori.p2.y = img.cols;
p_ori.p3.x = img.rows-100; p_ori.p3.y = img.cols-100;
p_ori.p4.x = img.rows; p_ori.p4.y = 0;
p_goal.p1.x = 0; p_goal.p1.y = 0;
p_goal.p2.x = 0; p_goal.p2.y = img.cols;
p_goal.p3.x = img.rows - 0; p_goal.p3.y = img.cols - 0;
p_goal.p4.x = img.rows; p_goal.p4.y = 0 + 0;
PT.persperctive_transformation(p_ori, p_goal, T);
for (int i = 0; i < img.rows; i++)
{
for (int j = 0; j < img.cols; j++)
{
float w = T[2] * i + T[5] * j + T[8];
float row = T[0] * i + T[3] * j + T[6];
float col = T[1] * i + T[4] * j + T[7];
int R = row / w;
int C = col / w;
if (R < 0 || C < 0 || R>img.rows || C>img.cols) continue;
res2.at<cv::Vec3b>(R, C) = res.at<cv::Vec3b>(i, j);
}
}
return 0;
}