[GRASS-SVN] r49846 - grass/trunk/misc/m.transform
svn_grass at osgeo.org
svn_grass at osgeo.org
Tue Dec 20 10:10:56 EST 2011
Author: mmetz
Date: 2011-12-20 07:10:56 -0800 (Tue, 20 Dec 2011)
New Revision: 49846
Removed:
grass/trunk/misc/m.transform/crs.c
grass/trunk/misc/m.transform/crs.h
Modified:
grass/trunk/misc/m.transform/main.c
Log:
use imagery lib for transformation
Deleted: grass/trunk/misc/m.transform/crs.c
===================================================================
--- grass/trunk/misc/m.transform/crs.c 2011-12-20 15:04:03 UTC (rev 49845)
+++ grass/trunk/misc/m.transform/crs.c 2011-12-20 15:10:56 UTC (rev 49846)
@@ -1,462 +0,0 @@
-
-/***************************************************************************/
-
-/***************************************************************************/
-/*
- CRS.C - Center for Remote Sensing rectification routines
-
- Written By: Brian J. Buckley
-
- At: The Center for Remote Sensing
- Michigan State University
- 302 Berkey Hall
- East Lansing, MI 48824
- (517)353-7195
-
- 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 <stdio.h>
-#include <string.h>
-#include <stdlib.h>
-#include <math.h>
-#include <limits.h>
-
-#include <grass/gis.h>
-#include <grass/imagery.h>
-
-#include "crs.h"
-
-/* 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 */
-
-/***************************************************************************/
-/*
- 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 CRS_georef(double e1, /* EASTINGS TO BE TRANSFORMED */
- double n1, /* NORTHINGS TO BE TRANSFORMED */
- double *e, /* EASTINGS TO BE TRANSFORMED */
- double *n, /* NORTHINGS TO BE 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 GEOREFFERENCING COEFFICIENTS BASED ON A SET OF CONTROL POINTS
- */
-
-/***************************************************************************/
-
-int CRS_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;
-}
Deleted: grass/trunk/misc/m.transform/crs.h
===================================================================
--- grass/trunk/misc/m.transform/crs.h 2011-12-20 15:04:03 UTC (rev 49845)
+++ grass/trunk/misc/m.transform/crs.h 2011-12-20 15:10:56 UTC (rev 49846)
@@ -1,30 +0,0 @@
-
-/***************************************************************************/
-
-/***************************************************************************/
-/*
- CRS.H - Center for Remote Sensing rectification routines
-
- Written By: Brian J. Buckley
-
- At: The Center for Remote Sensing
- Michigan State University
- 302 Berkey Hall
- East Lansing, MI 48824
- (517)353-7195
-
- Written: 12/19/91
-
- Last Update: 12/26/91 Brian J. Buckley
- */
-
-/***************************************************************************/
-
-/***************************************************************************/
-
-#define MAXORDER 3
-
-extern int CRS_compute_georef_equations(struct Control_Points *, double *,
- double *, double *, double *, int);
-extern int CRS_georef(double, double, double *, double *, double *, double *,
- int);
Modified: grass/trunk/misc/m.transform/main.c
===================================================================
--- grass/trunk/misc/m.transform/main.c 2011-12-20 15:04:03 UTC (rev 49845)
+++ grass/trunk/misc/m.transform/main.c 2011-12-20 15:10:56 UTC (rev 49846)
@@ -25,11 +25,6 @@
#include <grass/imagery.h>
#include <grass/glocale.h>
-extern int CRS_compute_georef_equations(struct Control_Points *, double *,
- double *, double *, double *, int);
-extern int CRS_georef(double, double, double *, double *, double *, double *,
- int);
-
struct Max
{
int idx;
@@ -91,7 +86,7 @@
int n, i;
equation_stat =
- CRS_compute_georef_equations(&points, E12, N12, E21, N21, order);
+ I_compute_georef_equations(&points, E12, N12, E21, N21, order);
if (equation_stat == 0)
G_fatal_error(_("Not enough points, %d are required"),
@@ -113,7 +108,7 @@
count++;
if (need_fwd) {
- CRS_georef(points.e1[n], points.n1[n], &e2, &n2, E12, N12, order);
+ I_georef(points.e1[n], points.n1[n], &e2, &n2, E12, N12, order);
fx = fabs(e2 - points.e2[n]);
fy = fabs(n2 - points.n2[n]);
@@ -126,7 +121,7 @@
}
if (need_rev) {
- CRS_georef(points.e2[n], points.n2[n], &e1, &n1, E21, N21, order);
+ I_georef(points.e2[n], points.n2[n], &e1, &n1, E21, N21, order);
rx = fabs(e1 - points.e1[n]);
ry = fabs(n1 - points.n1[n]);
@@ -260,9 +255,9 @@
double xe, xn;
if(forward)
- CRS_georef(east, north, &xe, &xn, E12, N12, order);
+ I_georef(east, north, &xe, &xn, E12, N12, order);
else
- CRS_georef(east, north, &xe, &xn, E21, N21, order);
+ I_georef(east, north, &xe, &xn, E21, N21, order);
fprintf(stdout, "%.15g %.15g\n", xe, xn);
}
More information about the grass-commit
mailing list