知源图四点坐标和目标图四点坐标，求透视变换矩阵

　　最近在搞图像处理，碰到了透视变换的问题。

　　同事给我一些代码，里边有误差，挺严重，让我帮他想想哪里出错了。捣鼓了很久，我猜测肯定是透视变换矩阵求错了，然后我的透视变换之旅就开始了。

　　后来问题解决了，发现他的矩阵和我求得矩阵一摸一样，他的代码并没有错误，是他采用的图片在做广角变换的时候有误差，导致程序结果误差。

　　首先感谢xiaowei_cqu，是她的一篇博客（http://blog.csdn.net/xiaowei_cqu/article/details/26471527）教会了我变换原理。

　　后来去用opencv来验证，发现函数调来调去太麻烦了，代码量也不小。

　　可是我只想要变换矩阵，自己写个不就好了！然后根据opencv源码里的原理，配上高斯消元自己写了个qwq。

　　talk is cheap,show me the code!!!

//Calculates coefficients of perspective transformation
/* 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
*/

void Gauss(double A[][9], int equ, int var, double* ans) { //epu:A‘s row  var:A‘s col-1
	int row, col;
	for (row = 0, col = 0; col<var&&row<equ; col++, row++) {
		int max_r = row;
		for (int i = row + 1; i<equ; i++) {
			if ((1e-12)<fabs(A[i][col]) - fabs(A[max_r][col])) {
				max_r = i;
			}
		}
		if (max_r != row)
			for (int j = 0; j<var + 1; j++)
				swap(A[row][j], A[max_r][j]);
		for (int i = row + 1; i<equ; i++) {
			if (fabs(A[i][col])<(1e-12))
				continue;
			double tmp = -A[i][col] / A[row][col];
			for (int j = col; j<var + 1; j++) {
				A[i][j] += tmp*A[row][j];
			}
		}

	}
	for (int i = var - 1; i >= 0; i--) { //计算唯一解。
		double tmp = 0;
		for (int j = i + 1; j<var; j++) {
			tmp += A[i][j] * (*(ans + j));
		}
		ans[i] = (A[i][var] - tmp) / A[i][i];
	}
}

void byx_getPerspectiveTransform(PointD * src, PointD * dst, double* ret){
	double x0 = src[0].x, x1 = src[1].x, x2 = src[3].x, x3 = src[2].x;
	double y0 = src[0].y, y1 = src[1].y, y2 = src[3].y, y3 = src[2].y;
	double u0 = dst[0].x, u1 = dst[1].x, u2 = dst[3].x, u3 = dst[2].x;
	double v0 = dst[0].y, v1 = dst[1].y, v2 = dst[3].y, v3 = dst[2].y;
	double A[8][9] = {
			{ x0, y0, 1, 0, 0, 0, -x0*u0, -y0*u0, u0 },
			{ x1, y1, 1, 0, 0, 0, -x1*u1, -y1*u1, u1 },
			{ x2, y2, 1, 0, 0, 0, -x2*u2, -y2*u2, u2 },
			{ x3, y3, 1, 0, 0, 0, -x3*u3, -y3*u3, u3 },
			{ 0, 0, 0, x0, y0, 1, -x0*v0, -y0*v0, v0 },
			{ 0, 0, 0, x1, y1, 1, -x1*v1, -y1*v1, v1 },
			{ 0, 0, 0, x2, y2, 1, -x2*v2, -y2*v2, v2 },
			{ 0, 0, 0, x3, y3, 1, -x3*v3, -y3*v3, v3 },
	};
	double ans[8] = { 0 };
	Gauss(A, 8, 8, ans);
	*(ret + 0) = ans[0];    *(ret + 1) = ans[1];    *(ret + 2) = 0;    *(ret + 3) = ans[2];
	*(ret + 4) = ans[3];    *(ret + 5) = ans[4];    *(ret + 6) = 0;    *(ret + 7) = ans[5];
	*(ret + 8) = 0;         *(ret + 9) = 0;         *(ret + 10) = 1;   *(ret + 11) = 0;
	*(ret + 12) = ans[6];   *(ret + 13) = ans[7];   *(ret + 14) = 0;   *(ret + 15) = 1;
}

PointD byx_Transform(PointD p,double * mat){
    PointD ret;
    double D = p.x*mat[12]+p.y*mat[13]+mat[15];
    ret.x=(p.x*mat[0]+p.y*mat[1]+mat[3])/D;
    ret.y=(p.x*mat[4]+p.y*mat[5]+mat[7])/D;
    return ret;
}

　　ps:变换矩阵是3*3的，我同事需要4*4的（第三行和第三列自行忽略）