射影几何(Projective Geometry)

1. 点

点\(P(x,y)\)二维坐标表示为 \[ P(x,y)\sim\begin{pmatrix} x \\ y \\ 1 \end{pmatrix} = \lambda \begin{pmatrix} x \\ y \\ 1 \end{pmatrix} = \begin{pmatrix} x' \\ y' \\ \lambda \end{pmatrix}\sim \begin{pmatrix} x \\ y \\ z \end{pmatrix} \] 若\(z = 0\) 那么，它表示在\(O(0,0)到P(x,y)\)的直线上的最远端的无穷点

2. 直线

一个三位向量可以表示一条直线，我们通常使用法向量。其中，\(A,B,C \neq0\) \[ Ax+By+C=0 \sim \begin{pmatrix} A \\ B \\ C \end{pmatrix} = \lambda\begin{pmatrix} A \\ B \\ C \end{pmatrix} \] \(\begin{pmatrix} 0 \\ 0 \\ 1 \end{pmatrix}\)代表无穷远处的直线。

2.1 点在直线上

\[ 直线：l=\begin{pmatrix} a \\ b \\ c \end{pmatrix},点：P=\begin{pmatrix} x \\ y \\ z \end{pmatrix},P在直线上,l \cdot P = 0 \]

2.2 过两点的直线

\[ P_{1}=\begin{pmatrix} x_{1} \\ y_{1} \\ z_{1} \end{pmatrix},P_{2}=\begin{pmatrix} x_{2} \\ y_{2} \\ z_{2} \end{pmatrix} \Rightarrow l = P_{1} \times P_{2} \]

2.3 两条直线的交点

\[ l_{1}=\begin{pmatrix} x_{1} \\ y_{1} \\ z_{1} \end{pmatrix},l_{2}=\begin{pmatrix} x_{2} \\ y_{2} \\ z_{2} \end{pmatrix} \Rightarrow P = l_{1} \times l_{2} \]

如果两条平行线相交，那么结果是一个无穷远的点。

3. 射影变换

通过一个矩阵可以实现任意射影变换, \(P'\)是点\(P\)变换后的点，那么有以下关系。 \[ P' = M \times P \]

image-20200116150606601

3.1 求变换矩阵M

我们需要4个点以及射影后的坐标，\(a,b,c,d,a',b',c',d'\)。其中,a,b,c,d坐标自己定，一般可以定成下面的 \[ a =\begin{pmatrix} 0 \\ 0 \\ 1 \end{pmatrix},b =\begin{pmatrix} 1 \\ 0 \\ 1 \end{pmatrix},c =\begin{pmatrix} 1 \\ 1 \\ 1 \end{pmatrix},d =\begin{pmatrix} 0 \\ 1 \\ 1 \end{pmatrix} \\ a' =\begin{pmatrix} a_{x} \\ a_{y} \\ 1 \end{pmatrix},b' =\begin{pmatrix} b_{x} \\ b_{y} \\ 1 \end{pmatrix},c' =\begin{pmatrix} c_{x} \\ c_{y} \\ 1 \end{pmatrix},d' =\begin{pmatrix} d_{x} \\ d_{y} \\ 1 \end{pmatrix} \] 以上都是已知数，可以列出下面4个关系式，因为\(z\)是1，所以要加\(\lambda\)这个参数 \[ M\times a=\lambda _{a}a' \\ M\times b=\lambda _{b}b' \\ M\times c=\lambda _{c}c' \\ M\times d=\lambda _{d}d' , \\ \lambda_{d}=1 \] 以a点举例，展开这个式子 \[ M=\begin{pmatrix} m_{1}& m_{2} &m_{3} \\ m_{4}& m_{5} &m_{6} \\ m_{7}& m_{8} & m_{9} \end{pmatrix} \times \begin{pmatrix} a_{x} \\ a_{y}\\ 1 \end{pmatrix}= \begin{pmatrix} m_{1}a_{x}+m_{2}a_{y}+m_{3} \\ m_{4}a_{x}+m_{5}a_{y}+m_{6}\\ m_{7}a_{x}+m_{8}a_{y}+m_{9} \end{pmatrix} =\begin{pmatrix} \lambda_{a}a_{x}' \\ \lambda_{a}a_{y}'\\ \lambda_{a} \end{pmatrix} \] 把以上的式子全部展开，可以得到一个12个未知数的线性方程组（矩阵9个+\(\lambda\)变量3个）。使用高斯消元法即可解出M。

3.2 还原射影图像的坐标

我们可以通过变换矩阵的逆矩阵，来反向变换，实现点的还原。 \[ M \times M^{-1}=I \\ a' \times M^{-1}=a \] 所以先求M，再求M的逆矩阵，然后相乘得到原来的点。

4. 代码实现（C++）

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
const int maxn = 1010;
const double ZERO = 1e-10;

namespace PG { // projective geometry

	struct Vector {
		int n;
		double data[maxn];
		
		Vector(int x=13) {
			n = x;
			for (int i = 1;i <= n; ++i) {
				data[i] = 0;
			}
		}
		
	};

	struct Matrix {
		int n, m;
		double data[maxn][maxn];
		Matrix(int x, int y) {
			n = x; m = y;
		}

		void multline(int line, double value) {
			for (int i = 1;i <= m; ++i) data[line][i] *= value;
		}

		void addwith(int line, int other) {
			for (int i = 1;i <= m; ++i) data[line][i] += data[other][i];
		}

		void appendline(Vector v) {
			for (int i = 1;i <= v.n; ++i) {
				data[n+1][i] = v.data[i];
			}
			n++;
		}

		void mul(double val) {
			for (int i = 1;i <= n; ++i) {
				for (int j = 1;j <= m; ++j) {
					data[i][j] *= val;
				}
			}
		}

		void add(Vector v) {
			int cnt = 1;
			for (int i = 1;i <= n; ++i) {
				for (int j = 1;j <= m; ++j) {
					data[i][j] = v.data[cnt++];
				}
			}
		}

		Matrix& operator = (const Matrix& other) {
			n = other.n; m = other.m;
			for (int i = 1;i <= n; ++i) {
				for (int j = 1;j <= m; ++j) {
					data[i][j] = other.data[i][j];
				}
			}
			return *this;
		}

		Vector gauss() {
			int vis[n+1], used[n+1];
			memset(vis,0,sizeof(vis)); memset(used,0,sizeof(used));

			for (int i = 1;i <= n; ++i) { // 选择主元
				int flag = 0;
				for (int j = 1;j <= n; ++j) { // 选择方程
					if (data[j][i] != 0 && !used[j]) {
						flag = 1;
						vis[i] = j;
						used[j] = 1;
						multline(j, 1.0/data[j][i]);
						for (int k = 1;k <= n; ++k) {
							if (used[k] || abs(data[k][i]) < ZERO) continue; // 不选择j
							multline(k, -1.0/data[k][i]);
							addwith(k, j); //消元i
						}
						break;
					}
				}
				if (!flag) return Vector(-1);
			}

			Vector ret(n);

			for (int i = n; i >= 1; --i) {
				int now = vis[i];
				double v = data[now][m];
				for (int j = 1;j <= n; ++j) {
					if (j == now) continue;
					data[j][m] -= data[j][i] * v;
					data[j][i] = 0;
				}
				ret.data[i] = v;
			}

			return ret;
		}

		void gauss_jordan() {
			for(int i = 1,r;i <= n; ++i){	
				r = i;
				for(int j=i+1;j<=n;++j) 
					if(abs(data[j][i]) > abs(data[r][i])) r = j;
				if(abs(data[r][i]) < ZERO) {
					n = -1; return; // no solution
				}
				if(i != r) swap(data[i],data[r]);
				
				for(int k = 1;k <= n; ++k){
					if(k == i || abs(data[k][i]) < ZERO) continue;
				
					multline(k, -data[i][i]/data[k][i]);
					addwith(k, i);
				
				} 
			}	
			for(int i = 1;i <= n; ++i) multline(i, 1.0/data[i][i]);
		}

		Matrix getI(int size) {
			Matrix ret(size, size);
			memset(ret.data,0,sizeof(ret.data));
			for (int i = 1;i <= size; ++i) ret.data[i][i] = 1;
			return ret;
		}

		void concat(Matrix other) {
			for (int i = 1;i <= other.n; ++i) {
				for (int j = 1;j <= other.m; ++j) {
					data[i][j+m] = other.data[i][j];
				}
			}
			m = m + other.m;
		}

		void spilt(int a, int b) { //分割位a-b列
			for (int i = 1;i <= n; ++i) {
				for (int j = 1;j <= b-a+1; ++j) {
					data[i][j] = data[i][a+j-1];
				}
			}
			m = b-a+1;
		}

		void print() {
			for (int i = 1;i <= n; ++i) {
				for (int j = 1;j <= m; ++j) {
					cout << data[i][j] << " ";
				}
				cout << endl;
			}
			cout << endl;
		}

		Matrix inverse() { // 求逆矩阵，高斯乔丹方法
			Matrix now(n,m); now = *this;
			now.concat(getI(n));
			now.gauss_jordan();
			now.spilt(now.n+1,now.m);
			return now;
		}

		Matrix operator *(Matrix other) {
			Matrix ret(n, other.m);
			for (int i = 1;i <= n; ++i) {
				for (int j = 1;j <= other.m; ++j) {
					ret.data[i][j] = 0;
					for (int k = 1;k <= m; ++k) ret.data[i][j] += data[i][k]*other.data[k][j];
				}
			}
			return ret;
		}
	};
}