/**************************************************************************** * * MODULE: imagery library * AUTHOR(S): Original author(s) name(s) unknown - written by CERL * Written By: Brian J. Buckley * * At: The Center for Remote Sensing * Michigan State University * * PURPOSE: Image processing library * COPYRIGHT: (C) 1999, 2005 by the GRASS Development Team * * This program is free software under the GNU General Public * License (>=v2). Read the file COPYING that comes with GRASS * for details. * *****************************************************************************/ /* * Written: 12/19/91 * * Last Update: 12/26/91 Brian J. Buckley * Last Update: 1/24/92 Brian J. Buckley * Added printout of trnfile. Triggered by BDEBUG. * Last Update: 1/27/92 Brian J. Buckley * Fixed bug so that only the active control points were used. * */ #include #include #include #include #include /* STRUCTURE FOR USE INTERNALLY WITH THESE FUNCTIONS. THESE FUNCTIONS EXPECT SQUARE MATRICES SO ONLY ONE VARIABLE IS GIVEN (N) FOR THE MATRIX SIZE */ struct MATRIX { int n; /* SIZE OF THIS MATRIX (N x N) */ double *v; }; /* CALCULATE OFFSET INTO ARRAY BASED ON R/C */ #define M(row,col) m->v[(((row)-1)*(m->n))+(col)-1] #define MSUCCESS 1 /* SUCCESS */ #define MNPTERR 0 /* NOT ENOUGH POINTS */ #define MUNSOLVABLE -1 /* NOT SOLVABLE */ #define MMEMERR -2 /* NOT ENOUGH MEMORY */ #define MPARMERR -3 /* PARAMETER ERROR */ #define MINTERR -4 /* INTERNAL ERROR */ #define MAXORDER 3 /* HIGHEST SUPPORTED ORDER OF TRANSFORMATION */ /*********************************************************************** FUNCTION PROTOTYPES FOR STATIC (INTERNAL) FUNCTIONS ************************************************************************/ static int calccoef(struct Control_Points *, double *, double *, int); static int calcls(struct Control_Points *, struct MATRIX *, double *, double *, double *, double *); static int exactdet(struct Control_Points *, struct MATRIX *, double *, double *, double *, double *); static int solvemat(struct MATRIX *, double *, double *, double *, double *); static double term(int, double, double); /*********************************************************************** TRANSFORM A SINGLE COORDINATE PAIR. ************************************************************************/ int I_georef(double e1, /* EASTING TO BE TRANSFORMED */ double n1, /* NORTHING TO BE TRANSFORMED */ double *e, /* EASTING, TRANSFORMED */ double *n, /* NORTHING, TRANSFORMED */ double E[], /* EASTING COEFFICIENTS */ double N[], /* NORTHING COEFFICIENTS */ int order /* ORDER OF TRANSFORMATION TO BE PERFORMED, MUST MATCH THE ORDER USED TO CALCULATE THE COEFFICIENTS */ ) { double e3, e2n, en2, n3, e2, en, n2; switch (order) { case 1: *e = E[0] + E[1] * e1 + E[2] * n1; *n = N[0] + N[1] * e1 + N[2] * n1; break; case 2: e2 = e1 * e1; n2 = n1 * n1; en = e1 * n1; *e = E[0] + E[1] * e1 + E[2] * n1 + E[3] * e2 + E[4] * en + E[5] * n2; *n = N[0] + N[1] * e1 + N[2] * n1 + N[3] * e2 + N[4] * en + N[5] * n2; break; case 3: e2 = e1 * e1; en = e1 * n1; n2 = n1 * n1; e3 = e1 * e2; e2n = e2 * n1; en2 = e1 * n2; n3 = n1 * n2; *e = E[0] + E[1] * e1 + E[2] * n1 + E[3] * e2 + E[4] * en + E[5] * n2 + E[6] * e3 + E[7] * e2n + E[8] * en2 + E[9] * n3; *n = N[0] + N[1] * e1 + N[2] * n1 + N[3] * e2 + N[4] * en + N[5] * n2 + N[6] * e3 + N[7] * e2n + N[8] * en2 + N[9] * n3; break; default: return MPARMERR; } return MSUCCESS; } /*********************************************************************** COMPUTE THE FORWARD AND BACKWARD GEOREFFERENCING COEFFICIENTS BASED ON A SET OF CONTROL POINTS ************************************************************************/ int I_compute_georef_equations(struct Control_Points *cp, double E12[], double N12[], double E21[], double N21[], int order) { double *tempptr; int status; if (order < 1 || order > MAXORDER) return MPARMERR; /* CALCULATE THE FORWARD TRANSFORMATION COEFFICIENTS */ status = calccoef(cp, E12, N12, order); if (status != MSUCCESS) return status; /* SWITCH THE 1 AND 2 EASTING AND NORTHING ARRAYS */ tempptr = cp->e1; cp->e1 = cp->e2; cp->e2 = tempptr; tempptr = cp->n1; cp->n1 = cp->n2; cp->n2 = tempptr; /* CALCULATE THE BACKWARD TRANSFORMATION COEFFICIENTS */ status = calccoef(cp, E21, N21, order); /* SWITCH THE 1 AND 2 EASTING AND NORTHING ARRAYS BACK */ tempptr = cp->e1; cp->e1 = cp->e2; cp->e2 = tempptr; tempptr = cp->n1; cp->n1 = cp->n2; cp->n2 = tempptr; return status; } /*********************************************************************** COMPUTE THE GEOREFFERENCING COEFFICIENTS BASED ON A SET OF CONTROL POINTS ************************************************************************/ static int calccoef(struct Control_Points *cp, double E[], double N[], int order) { struct MATRIX m; double *a; double *b; int numactive; /* NUMBER OF ACTIVE CONTROL POINTS */ int status, i; /* CALCULATE THE NUMBER OF VALID CONTROL POINTS */ for (i = numactive = 0; i < cp->count; i++) { if (cp->status[i] > 0) numactive++; } /* CALCULATE THE MINIMUM NUMBER OF CONTROL POINTS NEEDED TO DETERMINE A TRANSFORMATION OF THIS ORDER */ m.n = ((order + 1) * (order + 2)) / 2; if (numactive < m.n) return MNPTERR; /* INITIALIZE MATRIX */ m.v = G_calloc(m.n * m.n, sizeof(double)); a = G_calloc(m.n, sizeof(double)); b = G_calloc(m.n, sizeof(double)); if (numactive == m.n) status = exactdet(cp, &m, a, b, E, N); else status = calcls(cp, &m, a, b, E, N); G_free(m.v); G_free(a); G_free(b); return status; } /*********************************************************************** CALCULATE THE TRANSFORMATION COEFFICIENTS WITH EXACTLY THE MINIMUM NUMBER OF CONTROL POINTS REQUIRED FOR THIS TRANSFORMATION. ************************************************************************/ static int exactdet(struct Control_Points *cp, struct MATRIX *m, double a[], double b[], double E[], /* EASTING COEFFICIENTS */ double N[] /* NORTHING COEFFICIENTS */ ) { int pntnow, currow, j; currow = 1; for (pntnow = 0; pntnow < cp->count; pntnow++) { if (cp->status[pntnow] > 0) { /* POPULATE MATRIX M */ for (j = 1; j <= m->n; j++) M(currow, j) = term(j, cp->e1[pntnow], cp->n1[pntnow]); /* POPULATE MATRIX A AND B */ a[currow - 1] = cp->e2[pntnow]; b[currow - 1] = cp->n2[pntnow]; currow++; } } if (currow - 1 != m->n) return MINTERR; return solvemat(m, a, b, E, N); } /*********************************************************************** CALCULATE THE TRANSFORMATION COEFFICIENTS WITH MORE THAN THE MINIMUM NUMBER OF CONTROL POINTS REQUIRED FOR THIS TRANSFORMATION. THIS ROUTINE USES THE LEAST SQUARES METHOD TO COMPUTE THE COEFFICIENTS. ************************************************************************/ static int calcls(struct Control_Points *cp, struct MATRIX *m, double a[], double b[], double E[], /* EASTING COEFFICIENTS */ double N[] /* NORTHING COEFFICIENTS */ ) { int i, j, n, numactive = 0; /* INITIALIZE THE UPPER HALF OF THE MATRIX AND THE TWO COLUMN VECTORS */ for (i = 1; i <= m->n; i++) { for (j = i; j <= m->n; j++) M(i, j) = 0.0; a[i - 1] = b[i - 1] = 0.0; } /* SUM THE UPPER HALF OF THE MATRIX AND THE COLUMN VECTORS ACCORDING TO THE LEAST SQUARES METHOD OF SOLVING OVER DETERMINED SYSTEMS */ for (n = 0; n < cp->count; n++) { if (cp->status[n] > 0) { numactive++; for (i = 1; i <= m->n; i++) { for (j = i; j <= m->n; j++) M(i, j) += term(i, cp->e1[n], cp->n1[n]) * term(j, cp->e1[n], cp->n1[n]); a[i - 1] += cp->e2[n] * term(i, cp->e1[n], cp->n1[n]); b[i - 1] += cp->n2[n] * term(i, cp->e1[n], cp->n1[n]); } } } if (numactive <= m->n) return MINTERR; /* TRANSPOSE VALUES IN UPPER HALF OF M TO OTHER HALF */ for (i = 2; i <= m->n; i++) for (j = 1; j < i; j++) M(i, j) = M(j, i); return solvemat(m, a, b, E, N); } /*********************************************************************** CALCULATE THE X/Y TERM BASED ON THE TERM NUMBER ORDER\TERM 1 2 3 4 5 6 7 8 9 10 1 e0n0 e1n0 e0n1 2 e0n0 e1n0 e0n1 e2n0 e1n1 e0n2 3 e0n0 e1n0 e0n1 e2n0 e1n1 e0n2 e3n0 e2n1 e1n2 e0n3 ************************************************************************/ static double term(int term, double e, double n) { switch (term) { case 1: return 1.0; case 2: return e; case 3: return n; case 4: return e * e; case 5: return e * n; case 6: return n * n; case 7: return e * e * e; case 8: return e * e * n; case 9: return e * n * n; case 10: return n * n * n; } return 0.0; } /*********************************************************************** SOLVE FOR THE 'E' AND 'N' COEFFICIENTS BY USING A SOMEWHAT MODIFIED GAUSSIAN ELIMINATION METHOD. | M11 M12 ... M1n | | E0 | | a0 | | M21 M22 ... M2n | | E1 | = | a1 | | . . . . | | . | | . | | Mn1 Mn2 ... Mnn | | En-1 | | an-1 | and | M11 M12 ... M1n | | N0 | | b0 | | M21 M22 ... M2n | | N1 | = | b1 | | . . . . | | . | | . | | Mn1 Mn2 ... Mnn | | Nn-1 | | bn-1 | ************************************************************************/ static int solvemat(struct MATRIX *m, double a[], double b[], double E[], double N[]) { int i, j, i2, j2, imark; double factor, temp; double pivot; /* ACTUAL VALUE OF THE LARGEST PIVOT CANDIDATE */ for (i = 1; i <= m->n; i++) { j = i; /* find row with largest magnitude value for pivot value */ pivot = M(i, j); imark = i; for (i2 = i + 1; i2 <= m->n; i2++) { temp = fabs(M(i2, j)); if (temp > fabs(pivot)) { pivot = M(i2, j); imark = i2; } } /* if the pivot is very small then the points are nearly co-linear */ /* co-linear points result in an undefined matrix, and nearly */ /* co-linear points results in a solution with rounding error */ if (pivot == 0.0) return MUNSOLVABLE; /* if row with highest pivot is not the current row, switch them */ if (imark != i) { for (j2 = 1; j2 <= m->n; j2++) { temp = M(imark, j2); M(imark, j2) = M(i, j2); M(i, j2) = temp; } temp = a[imark - 1]; a[imark - 1] = a[i - 1]; a[i - 1] = temp; temp = b[imark - 1]; b[imark - 1] = b[i - 1]; b[i - 1] = temp; } /* compute zeros above and below the pivot, and compute values for the rest of the row as well */ for (i2 = 1; i2 <= m->n; i2++) { if (i2 != i) { factor = M(i2, j) / pivot; for (j2 = j; j2 <= m->n; j2++) M(i2, j2) -= factor * M(i, j2); a[i2 - 1] -= factor * a[i - 1]; b[i2 - 1] -= factor * b[i - 1]; } } } /* SINCE ALL OTHER VALUES IN THE MATRIX ARE ZERO NOW, CALCULATE THE COEFFICIENTS BY DIVIDING THE COLUMN VECTORS BY THE DIAGONAL VALUES. */ for (i = 1; i <= m->n; i++) { E[i - 1] = a[i - 1] / M(i, i); N[i - 1] = b[i - 1] / M(i, i); } return MSUCCESS; }