最近在搞图像处理,碰到了透视变换的问题。
我只想要变换矩阵,根据opencv源码里的原理,配上高斯消元自己写了个...qwq。
talk is cheap,show me the code!!!
1 /* Calculates coefficients of perspective transformation 2 * which maps (xi,yi) to (ui,vi), (i=1,2,3,4): 3 * 4 * c00*xi + c01*yi + c02 5 * ui = --------------------- 6 * c20*xi + c21*yi + c22 7 * 8 * c10*xi + c11*yi + c12 9 * vi = --------------------- 10 * c20*xi + c21*yi + c22 11 * 12 * Coefficients are calculated by solving linear system: 13 * / x0 y0 1 0 0 0 -x0*u0 -y0*u0 \ /c00\ /u0\ 14 * | x1 y1 1 0 0 0 -x1*u1 -y1*u1 | |c01| |u1| 15 * | x2 y2 1 0 0 0 -x2*u2 -y2*u2 | |c02| |u2| 16 * | x3 y3 1 0 0 0 -x3*u3 -y3*u3 |.|c10|=|u3|, 17 * | 0 0 0 x0 y0 1 -x0*v0 -y0*v0 | |c11| |v0| 18 * | 0 0 0 x1 y1 1 -x1*v1 -y1*v1 | |c12| |v1| 19 * | 0 0 0 x2 y2 1 -x2*v2 -y2*v2 | |c20| |v2| 20 * \ 0 0 0 x3 y3 1 -x3*v3 -y3*v3 / \c21/ \v3/ 21 * 22 * where: 23 * cij - matrix coefficients, c22 = 1 24 */ 25 26 void Gauss(double A[][9], int equ, int var, double* ans) { //epu:A's row var:A's col-1 27 int row, col; 28 for (row = 0, col = 0; col<var&&row<equ; col++, row++) { 29 int max_r = row; 30 for (int i = row + 1; i<equ; i++) { 31 if ((1e-12)<fabs(A[i][col]) - fabs(A[max_r][col])) { 32 max_r = i; 33 } 34 } 35 if (max_r != row) 36 for (int j = 0; j<var + 1; j++) 37 swap(A[row][j], A[max_r][j]); 38 for (int i = row + 1; i<equ; i++) { 39 if (fabs(A[i][col])<(1e-12)) 40 continue; 41 double tmp = -A[i][col] / A[row][col]; 42 for (int j = col; j<var + 1; j++) { 43 A[i][j] += tmp*A[row][j]; 44 } 45 } 46 47 } 48 for (int i = var - 1; i >= 0; i--) { //计算唯一解。 49 double tmp = 0; 50 for (int j = i + 1; j<var; j++) { 51 tmp += A[i][j] * (*(ans + j)); 52 } 53 ans[i] = (A[i][var] - tmp) / A[i][i]; 54 } 55 } 56 57 void byx_getPerspectiveTransform(Point2f * src, Point2f * dst, double* ret){ 58 double x0 = src[0].x, x1 = src[1].x, x2 = src[3].x, x3 = src[2].x; 59 double y0 = src[0].y, y1 = src[1].y, y2 = src[3].y, y3 = src[2].y; 60 double u0 = dst[0].x, u1 = dst[1].x, u2 = dst[3].x, u3 = dst[2].x; 61 double v0 = dst[0].y, v1 = dst[1].y, v2 = dst[3].y, v3 = dst[2].y; 62 double A[8][9] = { 63 { x0, y0, 1, 0, 0, 0, -x0*u0, -y0*u0, u0 }, 64 { x1, y1, 1, 0, 0, 0, -x1*u1, -y1*u1, u1 }, 65 { x2, y2, 1, 0, 0, 0, -x2*u2, -y2*u2, u2 }, 66 { x3, y3, 1, 0, 0, 0, -x3*u3, -y3*u3, u3 }, 67 { 0, 0, 0, x0, y0, 1, -x0*v0, -y0*v0, v0 }, 68 { 0, 0, 0, x1, y1, 1, -x1*v1, -y1*v1, v1 }, 69 { 0, 0, 0, x2, y2, 1, -x2*v2, -y2*v2, v2 }, 70 { 0, 0, 0, x3, y3, 1, -x3*v3, -y3*v3, v3 }, 71 }; 72 Gauss(A, 8, 8, ret); 73 *(ret + 8) = 1; 74 } 75 76 Point2f byx_Transform(Point2f p, double* mat){ 77 Point2f ret; 78 double D = p.x*mat[6] + p.y*mat[7] + mat[8]; 79 ret.x = (p.x*mat[0] + p.y*mat[1] + mat[2]) / D; 80 ret.y = (p.x*mat[3] + p.y*mat[4] + mat[5]) / D; 81 return ret; 82 }