add new module v.rectify

git-svn-id: https://svn.osgeo.org/grass/grass/trunk@49825 15284696-431f-4ddb-bdfa-cd5b030d7da7

vector/Makefile

@@ -70,6 +70,7 @@ SUBDIRS = \
 	v.qcount \
 	v.random \
 	v.reclass \
+	v.rectify \
 	v.sample \
 	v.segment \
 	v.select \

vector/v.rectify/Makefile

@@ -0,0 +1,10 @@
+MODULE_TOPDIR = ../..
+
+PGM = v.rectify
+
+LIBES = $(IMAGERYLIB) $(VECTORLIB) $(GISLIB)
+DEPENDENCIES = $(IMAGERYDEP) $(VECTORDEP) $(GISDEP)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+default: cmd

vector/v.rectify/cp.c

@@ -0,0 +1,361 @@
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <unistd.h>
+#include <fcntl.h>
+#include <grass/imagery.h>
+#include <grass/glocale.h>
+#include "global.h"
+#include "crs.h"
+
+struct Stats
+{
+    double x, y, z, g;
+    double sum2, rms;
+};
+
+static void update_stats(struct Stats *st, int n,
+                         double dx, double dy, double *dz,
+			 double dg, double d2)
+{
+    st->x += dx;
+    st->y += dy;
+    if (dz)
+	st->z += *dz;
+    st->g += dg;
+    st->sum2 += d2;
+}
+
+static void diagonal(double *dg, double *d2, double dx, double dy, double *dz)
+{
+    *d2 = dx * dx + dy * dy;
+    if (dz)
+	*d2 += *dz * *dz;
+    *dg = sqrt(*d2);
+}
+
+
+static void compute_rms(struct Control_Points *cp, struct Control_Points_3D *cp3,
+                        int order, int use3d, char *sep)
+{
+    int n;
+    int count, npoints;
+    struct Stats fwd, rev;
+
+    fwd.sum2 = fwd.rms = rev.sum2 = rev.rms = 0.0;
+    
+    fwd.x = fwd.y = fwd.z = fwd.g = 0.0;
+    rev.x = rev.y = rev.z = rev.g = 0.0;
+
+    count = 0;
+    
+    /* print index, forward difference, backward difference, 
+     * forward rms, backward rms */
+    if (use3d)
+	printf("index%sfwd_dx%sfwd_dy%sfwd_dz%sback_dx%sback_dy%sback_dz%sfwd_RMS%sback_RMS",
+               sep, sep, sep, sep, sep, sep, sep, sep);
+    else
+	printf("index%sfwd_dx%sfwd_dy%sback_dx%sback_dy%sfwd_RMS%sback_RMS",
+	       sep, sep, sep, sep, sep, sep);
+
+    printf("\n");
+    
+    if (use3d)
+	npoints = cp3->count;
+    else
+	npoints = cp->count;
+
+    for (n = 0; n < npoints; n++) {
+	double e1, n1, z1, e2, n2, z2;
+	double fx, fy, fz, fd, fd2;
+	double rx, ry, rz, rd, rd2;
+
+	if (use3d) {
+	    if (cp3->status[n] <= 0)
+		continue;
+	}
+	else {
+	    if (cp->status[n] <= 0)
+		continue;
+	}
+
+	count++;
+
+	/* forward: source -> target */
+	if (use3d) {
+	    CRS_georef_3d(cp3->e1[n], cp3->n1[n], cp3->z1[n],
+			  &e2, &n2, &z2, 
+			  E12, N12, Z12, 
+			  order);
+
+	    fx = fabs(e2 - cp3->e2[n]);
+	    fy = fabs(n2 - cp3->n2[n]);
+	    fz = fabs(z2 - cp3->z2[n]);
+
+	    diagonal(&fd, &fd2, fx, fy, &fz);
+
+	    update_stats(&fwd, n, fx, fy, &fz, fd, fd2);
+	}
+	else {
+	    CRS_georef(cp->e1[n], cp->n1[n], &e2, &n2, E12, N12, order);
+
+	    fx = fabs(e2 - cp->e2[n]);
+	    fy = fabs(n2 - cp->n2[n]);
+
+	    diagonal(&fd, &fd2, fx, fy, NULL);
+
+	    update_stats(&fwd, n, fx, fy, NULL, fd, fd2);
+	}
+
+	/* backward: target -> source */
+	if (use3d) {
+	    CRS_georef_3d(cp3->e2[n], cp3->n2[n], cp3->z2[n],
+			  &e1, &n1, &z1,
+			  E21, N21, Z21,
+			  order);
+
+	    rx = fabs(e1 - cp3->e1[n]);
+	    ry = fabs(n1 - cp3->n1[n]);
+	    rz = fabs(z1 - cp3->z1[n]);
+
+	    diagonal(&rd, &rd2, rx, ry, &rz);
+
+	    update_stats(&rev, n, rx, ry, &rz, rd, rd2);
+	}
+	else {
+	    CRS_georef(cp->e2[n], cp->n2[n], &e1, &n1, E21, N21, order);
+
+	    rx = fabs(e1 - cp->e1[n]);
+	    ry = fabs(n1 - cp->n1[n]);
+
+	    diagonal(&rd, &rd2, rx, ry, NULL);
+
+	    update_stats(&rev, n, rx, ry, NULL, rd, rd2);
+	}
+
+	/* print index, forward difference, backward difference, 
+	 * forward rms, backward rms */
+	printf("%d", n + 1);
+	printf("%s%f%s%f", sep, fx, sep, fy);
+	if (use3d)
+	    printf("%s%f", sep, fz);
+	printf("%s%f%s%f", sep, rx, sep, ry);
+	if (use3d)
+	    printf("%s%f", sep, rz);
+	printf("%s%.4f", sep, fd);
+	printf("%s%.4f", sep, rd);
+
+	printf("\n");
+    }
+
+    if (count > 0) {
+	fwd.x /= count;
+	fwd.y /= count;
+	fwd.g /= count;
+	rev.x /= count;
+	rev.y /= count;
+	rev.g /= count;
+	if (use3d) {
+	    fwd.z /= count;
+	    rev.z /= count;
+	}
+	fwd.rms = sqrt(fwd.sum2 / count);
+	rev.rms = sqrt(rev.sum2 / count);
+    }
+    printf("%d", count);
+    printf("%s%f%s%f", sep, fwd.x, sep, fwd.y);
+    if (use3d)
+	printf("%s%f", sep, fwd.z);
+    printf("%s%f%s%f", sep, rev.x, sep, rev.y);
+    if (use3d)
+	printf("%s%f", sep, rev.z);
+    printf("%s%.4f", sep, fwd.rms);
+    printf("%s%.4f", sep, rev.rms);
+
+    printf("\n");
+}
+
+
+int new_control_point_3d(struct Control_Points_3D *cp,
+			double e1, double n1, double z1,
+			double e2, double n2, double z2,
+			int status)
+{
+    int i;
+    unsigned int size;
+
+    if (status < 0)
+	return 1;
+    i = (cp->count)++;
+    size = cp->count * sizeof(double);
+    cp->e1 = (double *)G_realloc(cp->e1, size);
+    cp->e2 = (double *)G_realloc(cp->e2, size);
+    cp->n1 = (double *)G_realloc(cp->n1, size);
+    cp->n2 = (double *)G_realloc(cp->n2, size);
+    cp->z1 = (double *)G_realloc(cp->z1, size);
+    cp->z2 = (double *)G_realloc(cp->z2, size);
+    size = cp->count * sizeof(int);
+    cp->status = (int *)G_realloc(cp->status, size);
+
+    cp->e1[i] = e1;
+    cp->e2[i] = e2;
+    cp->n1[i] = n1;
+    cp->n2[i] = n2;
+    cp->z1[i] = z1;
+    cp->z2[i] = z2;
+    cp->status[i] = status;
+
+    return 0;
+}
+
+static int read_control_points(FILE * fd, struct Control_Points *cp)
+{
+    char buf[1000];
+    double e1, e2, n1, n2;
+    int status;
+
+    cp->count = 0;
+
+    /* read the control point lines. format is:
+       image_east image_north  target_east target_north  status
+     */
+    cp->e1 = NULL;
+    cp->e2 = NULL;
+    cp->n1 = NULL;
+    cp->n2 = NULL;
+    cp->status = NULL;
+
+    while (G_getl2(buf, sizeof buf, fd)) {
+	G_strip(buf);
+	if (*buf == '#' || *buf == 0)
+	    continue;
+	if (sscanf(buf, "%lf%lf%lf%lf%d", &e1, &n1, &e2, &n2, &status) == 5)
+	    I_new_control_point(cp, e1, n1, e2, n2, status);
+	else
+	    return -4;
+    }
+
+    return 1;
+}
+
+static int read_control_points_3d(FILE * fd, struct Control_Points_3D *cp)
+{
+    char buf[1000];
+    double e1, e2, n1, n2, z1, z2;
+    int status;
+
+    cp->count = 0;
+
+    /* read the control point lines. format is:
+       source_east source_north source_height  target_east target_north target_height  status
+     */
+    cp->e1 = NULL;
+    cp->e2 = NULL;
+    cp->n1 = NULL;
+    cp->n2 = NULL;
+    cp->z1 = NULL;
+    cp->z2 = NULL;
+    cp->status = NULL;
+
+    while (G_getl2(buf, sizeof buf, fd)) {
+	G_strip(buf);
+	if (*buf == '#' || *buf == 0)
+	    continue;
+	if (sscanf(buf, "%lf%lf%lf%lf%lf%lf%d", &e1, &n1, &z1, &e2, &n2, &z2, &status) == 7)
+	    new_control_point_3d(cp, e1, n1, z1, e2, n2, z2, status);
+	else
+	    return -4;
+    }
+
+    return 1;
+}
+
+
+int get_control_points(char *group, char *pfile, int order, int use3d, int rms, char *sep)
+{
+    char msg[200];
+    struct Control_Points cp;
+    struct Control_Points_3D cp3;
+    int ret = 0;
+    int order_pnts[2][3] = {{ 3, 6, 10 }, { 4, 10, 20 }};
+
+    if (use3d) {
+	/* read 3D GCPs from points file */
+	FILE *fp;
+	int fd, stat;
+	
+	if ((fd = open(pfile, 0)) < 0)
+	    G_fatal_error(_("Can not open file <%s>"), pfile);
+	    
+	fp = fdopen(fd, "r");
+
+	stat = read_control_points_3d(fp, &cp3);
+	fclose(fp);
+	if (stat < 0) {
+	    G_fatal_error(_("Bad format in control point file <%s>"),
+		      pfile);
+	    return 0;
+	}
+
+	ret = CRS_compute_georef_equations_3d(&cp3, E12, N12, Z12, E21, N21, Z21, order);
+    }
+    else if (pfile) {
+	/* read 2D GCPs from points file */
+	FILE *fp;
+	int fd, stat;
+	
+	if ((fd = open(pfile, 0)) < 0)
+	    G_fatal_error(_("Can not open file <%s>"), pfile);
+	    
+	fp = fdopen(fd, "r");
+
+	stat = read_control_points(fp, &cp);
+	fclose(fp);
+	if (stat < 0) {
+	    G_fatal_error(_("Bad format in control point file <%s>"),
+		      pfile);
+	    return 0;
+	}
+
+	ret = CRS_compute_georef_equations(&cp, E12, N12, E21, N21, order);
+    }
+    else {
+	/* read group control points */
+	if (!I_get_control_points(group, &cp))
+	    exit(0);
+
+	sprintf(msg, _("Control Point file for group <%s@%s> - "),
+		group, G_mapset());
+		
+	ret = CRS_compute_georef_equations(&cp, E12, N12, E21, N21, order);
+    }
+
+    switch (ret) {
+    case 0:
+	sprintf(&msg[strlen(msg)],
+		_("Not enough active control points for current order, %d are required."),
+		order_pnts[use3d != 0][order]);
+	break;
+    case -1:
+	strcat(msg, _("Poorly placed control points."));
+	strcat(msg, _(" Can not generate the transformation equation."));
+	break;
+    case -2:
+	strcat(msg, _("Not enough memory to solve for transformation equation"));
+	break;
+    case -3:
+	strcat(msg, _("Invalid order"));
+	break;
+    default:
+	break;
+    }
+    if (ret != 1)
+	G_fatal_error(msg);
+	
+    if (rms) {
+	compute_rms(&cp, &cp3, order, use3d, sep);
+    }
+	
+
+    return 1;
+}

vector/v.rectify/crs.c

@@ -0,0 +1,446 @@
+
+/***********************************************************************
+
+   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,	/* 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 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;
+}

vector/v.rectify/crs.h

@@ -0,0 +1,54 @@
+
+/***************************************************************************/
+
+/***************************************************************************/
+/*
+   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
+
+struct Control_Points_3D
+{
+    int count;
+    double *e1;
+    double *n1;
+    double *z1;
+    double *e2;
+    double *n2;
+    double *z2;
+    int *status;
+};
+
+
+/* crs.c */
+int CRS_compute_georef_equations(struct Control_Points *, double *,
+					double *, double *, double *, int);
+int CRS_georef(double, double, double *, double *, double *, double *,
+		      int);
+
+/* crs3d.c */
+int CRS_compute_georef_equations_3d(struct Control_Points_3D *,
+                                    double *, double *, double *,
+				    double *, double *, double *,
+				    int);
+int CRS_georef_3d(double, double, double,
+                  double *, double *, double *,
+		  double *, double *, double *,
+		  int);

vector/v.rectify/crs3d.c

@@ -0,0 +1,530 @@
+
+/***********************************************************************
+
+   crs3d.c
+
+   written by: Markus Metz
+
+   based on crs.c - Center for Remote Sensing rectification routines
+
+************************************************************************/
+
+#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_3D *, double *, double *, double *, int);
+static int calcls(struct Control_Points_3D *, struct MATRIX *,
+                  double *, double *, double *,
+		  double *, double *, double *);
+static int exactdet(struct Control_Points_3D *, struct MATRIX *,
+                    double *, double *, double *,
+		    double *, double *, double *);
+static int solvemat(struct MATRIX *, double *, double *, double *,
+                    double *, double *, double *);
+static double term(int, double, double, double);
+
+/***********************************************************************
+
+  TRANSFORM A SINGLE COORDINATE PAIR.
+
+************************************************************************/
+
+int CRS_georef_3d(double e1,	/* EASTING TO BE TRANSFORMED */
+	          double n1,	/* NORTHING TO BE TRANSFORMED */
+	          double z1,	/* HEIGHT TO BE TRANSFORMED */
+	          double *e,	/* EASTING, TRANSFORMED */
+	          double *n,	/* NORTHING, TRANSFORMED */
+	          double *z,	/* HEIGHT, TRANSFORMED */
+	          double E[],	/* EASTING COEFFICIENTS */
+	          double N[],	/* NORTHING COEFFICIENTS */
+	          double Z[],	/* HEIGHT COEFFICIENTS */
+		  int order     /* ORDER OF TRANSFORMATION TO BE PERFORMED, MUST MATCH THE
+				   ORDER USED TO CALCULATE THE COEFFICIENTS */
+    )
+{
+    double e2, n2, z2, en, ez, nz,
+           e3, n3, z3, e2n, e2z, en2, ez2, n2z, nz2, enz;
+
+    switch (order) {
+    case 1:
+	*e = E[0] + E[1] * e1 + E[2] * n1 + E[3] * z1;
+	*n = N[0] + N[1] * e1 + N[2] * n1 + N[3] * z1;
+	*z = Z[0] + Z[1] * e1 + Z[2] * n1 + Z[3] * z1;
+	break;
+
+    case 2:
+	e2 = e1 * e1;
+	en = e1 * n1;
+	ez = e1 * z1;
+	n2 = n1 * n1;
+	nz = n1 * z1;
+	z2 = z1 * z1;
+
+	*e = E[0] + E[1] * e1 + E[2] * n1 + E[3] * z1 +
+	     E[4] * e2 + E[5] * en + E[6] * ez + E[7] * n2 + E[8] * nz + E[9] * z2;
+	*n = N[0] + N[1] * e1 + N[2] * n1 + N[3] * z1 +
+	     N[4] * e2 + N[5] * en + N[6] * ez + N[7] * n2 + N[8] * nz + N[9] * z2;
+	*z = Z[0] + Z[1] * e1 + Z[2] * n1 + Z[3] * z1 +
+	     Z[4] * e2 + Z[5] * en + Z[6] * ez + Z[7] * n2 + Z[8] * nz + Z[9] * z2;
+	break;
+
+    case 3:
+	e2 = e1 * e1;
+	en = e1 * n1;
+	ez = e1 * z1;
+	n2 = n1 * n1;
+	nz = n1 * z1;
+	z2 = z1 * z1;
+
+	e3 = e1 * e2;
+	e2n = e2 * n1;
+	e2z = e2 * z1;
+	en2 = e1 * n2;
+	enz = e1 * n1 * z1;
+	ez2 = e1 * z2;
+	n3 = n1 * n2;
+	n2z = n2 * z1;
+	nz2 = n1 * z2;
+	z3 = z1 * z2;
+
+	*e = E[0] + E[1] * e1 + E[2] * n1 + E[3] * z1 +
+	     E[4] * e2 + E[5] * en + E[6] * ez + E[7] * n2 + E[8] * nz + E[9] * z2 +
+	     E[10] * e3 + E[11] * e2n + E[12] * e2z + E[13] * en2 + E[14] * enz +
+	     E[15] * ez2 + E[16] * n3 + E[17] * n2z + E[18] * nz2 + E[19] * z3;
+	*n = N[0] + N[1] * e1 + N[2] * n1 + N[3] * z1 +
+	     N[4] * e2 + N[5] * en + N[6] * ez + N[7] * n2 + N[8] * nz + N[9] * z2 +
+	     N[10] * e3 + N[11] * e2n + N[12] * e2z + N[13] * en2 + N[14] * enz +
+	     N[15] * ez2 + N[16] * n3 + N[17] * n2z + N[18] * nz2 + N[19] * z3;
+	*z = Z[0] + Z[1] * e1 + Z[2] * n1 + Z[3] * z1 +
+	     Z[4] * e2 + Z[5] * en + Z[6] * ez + Z[7] * n2 + Z[8] * nz + Z[9] * z2 +
+	     Z[10] * e3 + Z[11] * e2n + Z[12] * e2z + Z[13] * en2 + Z[14] * enz +
+	     Z[15] * ez2 + Z[16] * n3 + Z[17] * n2z + Z[18] * nz2 + Z[19] * z3;
+	break;
+
+    default:
+	return MPARMERR;
+    }
+
+    return MSUCCESS;
+}
+
+/***********************************************************************
+
+  COMPUTE THE FORWARD AND BACKWARD GEOREFFERENCING COEFFICIENTS
+  BASED ON A SET OF CONTROL POINTS
+
+************************************************************************/
+
+int CRS_compute_georef_equations_3d(struct Control_Points_3D *cp,
+                                    double E12[], double N12[], double Z12[],
+				    double E21[], double N21[], double Z21[],
+				    int order)
+{
+    double *tempptr;
+    int status;
+
+    /*
+    if (order < 1 || order > MAXORDER)
+	return MPARMERR;
+    */
+
+    /* CALCULATE THE FORWARD TRANSFORMATION COEFFICIENTS */
+
+    status = calccoef(cp, E12, N12, Z12, order);
+
+    if (status != MSUCCESS)
+	return status;
+
+    /* SWITCH THE 1 AND 2 EASTING, NORTHING, AND HEIGHT ARRAYS */
+
+    tempptr = cp->e1;
+    cp->e1 = cp->e2;
+    cp->e2 = tempptr;
+    tempptr = cp->n1;
+    cp->n1 = cp->n2;
+    cp->n2 = tempptr;
+    tempptr = cp->z1;
+    cp->z1 = cp->z2;
+    cp->z2 = tempptr;
+
+    /* CALCULATE THE BACKWARD TRANSFORMATION COEFFICIENTS */
+
+    status = calccoef(cp, E21, N21, Z21, order);
+
+    /* SWITCH THE 1 AND 2 EASTING, NORTHING, AND HEIGHT ARRAYS BACK */
+
+    tempptr = cp->e1;
+    cp->e1 = cp->e2;
+    cp->e2 = tempptr;
+    tempptr = cp->n1;
+    cp->n1 = cp->n2;
+    cp->n2 = tempptr;
+    tempptr = cp->z1;
+    cp->z1 = cp->z2;
+    cp->z2 = tempptr;
+
+    return status;
+}
+
+/***********************************************************************
+
+  COMPUTE THE GEOREFFERENCING COEFFICIENTS
+  BASED ON A SET OF CONTROL POINTS
+
+************************************************************************/
+
+static int calccoef(struct Control_Points_3D *cp,
+                    double E[], double N[], double Z[],
+		    int order)
+{
+    struct MATRIX m;
+    double *a;
+    double *b;
+    double *c;
+    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 */
+
+    /*
+		    2D    3D
+       1st order:    3     4
+       2nd order:    6    10
+       3rd order:   10	  20
+    */
+
+    if (order == 1)
+	m.n = 4;
+    else if (order == 2)
+	m.n = 10;
+    else if (order == 3)
+	m.n = 20;
+    
+    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));
+    c = G_calloc(m.n, sizeof(double));
+
+    if (numactive == m.n)
+	status = exactdet(cp, &m, a, b, c, E, N, Z);
+    else
+	status = calcls(cp, &m, a, b, c, E, N, Z);
+
+    G_free(m.v);
+    G_free(a);
+    G_free(b);
+    G_free(c);
+
+    return status;
+}
+
+/***********************************************************************
+
+  CALCULATE THE TRANSFORMATION COEFFICIENTS WITH EXACTLY THE MINIMUM
+  NUMBER OF CONTROL POINTS REQUIRED FOR THIS TRANSFORMATION.
+
+************************************************************************/
+
+static int exactdet(struct Control_Points_3D *cp, struct MATRIX *m,
+                    double a[], double b[], double c[],
+		    double E[],	/* EASTING COEFFICIENTS */
+		    double N[],	/* NORTHING COEFFICIENTS */
+		    double Z[]	/* HEIGHT 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], cp->z1[pntnow]);
+
+	    /* POPULATE MATRIX A AND B */
+
+	    a[currow - 1] = cp->e2[pntnow];
+	    b[currow - 1] = cp->n2[pntnow];
+	    c[currow - 1] = cp->z2[pntnow];
+
+	    currow++;
+	}
+    }
+
+    if (currow - 1 != m->n)
+	return MINTERR;
+
+    return solvemat(m, a, b, c, E, N, Z);
+}
+
+/***********************************************************************
+
+  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_3D *cp, struct MATRIX *m,
+                  double a[], double b[], double c[],
+		  double E[],	/* EASTING COEFFICIENTS */
+		  double N[],	/* NORTHING COEFFICIENTS */
+		  double Z[]	/* HEIGHT 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] = c[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], cp->z1[n]) *
+		               term(j, cp->e1[n], cp->n1[n], cp->z1[n]);
+
+		a[i - 1] += cp->e2[n] * term(i, cp->e1[n], cp->n1[n], cp->z1[n]);
+		b[i - 1] += cp->n2[n] * term(i, cp->e1[n], cp->n1[n], cp->z1[n]);
+		c[i - 1] += cp->z2[n] * term(i, cp->e1[n], cp->n1[n], cp->z1[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, c, E, N, Z);
+}
+
+/***********************************************************************
+
+  CALCULATE THE X/Y TERM BASED ON THE TERM NUMBER
+
+  ORDER\TERM   1      2      3      4      5      6      7      8      9      10
+  1            e0n0z0 e1n0z0 e0n1z0 e0n0z1
+  2            e0n0z0 e1n0z0 e0n1z0 e0n0z1 e2n0z0 e1n1z0 e1n0z1 e0n2z0 e0n1z1 e0n0z2
+  3            e0n0z0 e1n0z0 e0n1z0 e0n0z1 e2n0z0 e1n1z0 e1n0z1 e0n2z0 e0n1z1 e0n0z2
+
+  ORDER\TERM   11      12     13     14     15     16     17     18     19     20
+  3            e3n0z0  e2n1z0 e2n0z1 e1n2z0 e1n1z1 e1n0z2 e0n3z0 e0n2z1 e0n1z2 e0n0z3
+
+************************************************************************/
+
+static double term(int term, double e, double n, double z)
+{
+    switch (term) {
+    /* 1st order */
+    case 1:
+	return 1.0;
+    case 2:
+	return e;
+    case 3:
+	return n;
+    case 4:
+	return z;
+    /* 2nd order */
+    case 5:
+	return e * e;
+    case 6:
+	return e * n;
+    case 7:
+	return e * z;
+    case 8:
+	return n * n;
+    case 9:
+	return n * z;
+    case 10:
+	return z * z;
+    /* 3rd order */
+    case 11:
+	return e * e * e;
+    case 12:
+	return e * e * n;
+    case 13:
+	return e * e * z;
+    case 14:
+	return e * n * n;
+    case 15:
+	return e * n * z;
+    case 16:
+	return e * z * z;
+    case 17:
+	return n * n * n;
+    case 18:
+	return n * n * z;
+    case 19:
+	return n * z * z;
+    case 20:
+	return z * z * z;
+    }
+
+    return 0.0;
+}
+
+/***********************************************************************
+
+  SOLVE FOR THE 'E', 'N' AND 'Z' 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 |
+
+  ,
+
+  | M11 M12 ... M1n | | N0   |   | b0   |
+  | M21 M22 ... M2n | | N1   | = | b1   |
+  |  .   .   .   .  | | .    |   | .    |
+  | Mn1 Mn2 ... Mnn | | Nn-1 |   | bn-1 |
+
+  and
+
+  | M11 M12 ... M1n | | Z0   |   | c0   |
+  | M21 M22 ... M2n | | Z1   | = | c1   |
+  |  .   .   .   .  | | .    |   | .    |
+  | Mn1 Mn2 ... Mnn | | Zn-1 |   | cn-1 |
+
+************************************************************************/
+
+static int solvemat(struct MATRIX *m, double a[], double b[], double c[],
+                    double E[], double N[], double Z[])
+{
+    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;
+
+	    temp = c[imark - 1];
+	    c[imark - 1] = c[i - 1];
+	    c[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];
+		c[i2 - 1] -= factor * c[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);
+	Z[i - 1] = c[i - 1] / M(i, i);
+    }
+
+    return MSUCCESS;
+}

vector/v.rectify/env.c

@@ -0,0 +1,42 @@
+#include <unistd.h>
+#include "global.h"
+static int which_env = -1;	/* 0 = cur, 1 = target */
+
+int select_current_env(void)
+{
+    if (which_env < 0) {
+	G__create_alt_env();
+	which_env = 0;
+    }
+    if (which_env != 0) {
+	G__switch_env();
+	which_env = 0;
+    }
+
+    return 0;
+}
+
+int select_target_env(void)
+{
+    if (which_env < 0) {
+	G__create_alt_env();
+	which_env = 1;
+    }
+    if (which_env != 1) {
+	G__switch_env();
+	which_env = 1;
+    }
+
+    return 0;
+}
+
+int show_env(void)
+{
+    fprintf(stderr, "env(%d) switch to LOCATION %s, MAPSET %s\n", which_env,
+	    G__getenv("LOCATION_NAME") ==
+	    NULL ? "?" : G__getenv("LOCATION_NAME"),
+	    G__getenv("MAPSET") == NULL ? "?" : G__getenv("MAPSET"));
+    G_sleep(2);
+
+    return 0;
+}

vector/v.rectify/global.h

@@ -0,0 +1,19 @@
+#include <grass/gis.h>
+#include <grass/imagery.h>
+
+
+/* georef coefficients */
+
+extern double E12[20], N12[20], Z12[20];
+extern double E21[20], N21[20], Z21[20];
+
+/* cp.c */
+int get_control_points(char *, char *, int, int, int, char *);
+
+/* env.c */
+int select_current_env(void);
+int select_target_env(void);
+int show_env(void);
+
+/* target.c */
+int get_target(char *);

vector/v.rectify/main.c

@@ -0,0 +1,255 @@
+
+/****************************************************************************
+ *
+ * MODULE:       v.rectify
+ * AUTHOR(S):    Markus Metz
+ *               based on i.rectify
+ * PURPOSE:      calculate a transformation matrix and then convert x,y(,z) 
+ *               coordinates to standard map coordinates for all objects in 
+ *               the vector
+ *               control points can come from i.points or i.vpoints or 
+ *               a user-given text file
+ * COPYRIGHT:    (C) 2002-2011 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.
+ *
+ *****************************************************************************/
+
+#include <stdlib.h>
+#include <string.h>
+
+#include <grass/vector.h>
+#include <grass/glocale.h>
+
+#include "global.h"
+#include "crs.h"
+
+
+/* georef coefficients */
+
+double E12[20], N12[20], Z12[20];
+double E21[20], N21[20], Z21[20];
+
+
+int main(int argc, char *argv[])
+{
+    char group[INAME_LEN];
+    int order;			/* ADDED WITH CRS MODIFICATIONS */
+    int n, i, nlines, type;
+    int target_overwrite = 0;
+    char *points_file, *overstr, *rms_sep;
+    struct Map_info In, Out;
+    struct line_pnts *Points, *OPoints;
+    struct line_cats *Cats;
+    double x, y, z;
+    int use3d;
+
+    struct Option *grp,         /* imagery group */
+     *val,                      /* transformation order */
+     *in_opt,			/* input vector name */
+     *out_opt,			/* output vector name */
+     *pfile,			/* text file with GCPs */
+     *sep;			/* field separator for RMS report */
+
+    struct Flag *flag_use3d, *no_topo, *print_rms;
+
+    struct GModule *module;
+
+    G_gisinit(argv[0]);
+
+    module = G_define_module();
+    G_add_keyword(_("vector"));
+    G_add_keyword(_("rectify"));
+    module->description =
+	_("Rectifies a vector by computing a coordinate "
+	  "transformation for each object in the vector based on the "
+	  "control points.");
+
+    in_opt = G_define_standard_option(G_OPT_V_INPUT);
+    in_opt->required = YES;
+
+    out_opt = G_define_standard_option(G_OPT_V_OUTPUT);
+    out_opt->required = YES;
+
+    grp = G_define_standard_option(G_OPT_I_GROUP);
+    grp->required = NO;
+
+    pfile = G_define_standard_option(G_OPT_F_INPUT);
+    pfile->key = "points";
+    pfile->required = NO;
+
+    val = G_define_option();
+    val->key = "order";
+    val->type = TYPE_INTEGER;
+    val->required = YES;
+    val->description = _("Rectification polynom order (1-3)");
+    
+    sep = G_define_standard_option(G_OPT_F_SEP);
+    sep->label = _("Field separator for RMS report");
+    sep->description = _("Special characters: newline, space, comma, tab");
+
+    flag_use3d = G_define_flag();
+    flag_use3d->key = '3';
+    flag_use3d->description = _("Perform 3D transformation");
+
+    print_rms = G_define_flag();
+    print_rms->key = 'r';
+    print_rms->label = _("Print RMS errors");
+    print_rms->description = _("Print RMS errors and exit without rectifying the input map");
+
+    no_topo = G_define_flag();
+    no_topo->key = 'b';
+    no_topo->description = _("Do not build topology for output");
+
+    if (G_parser(argc, argv))
+	exit(EXIT_FAILURE);
+
+
+    if (grp->answer) {
+	G_strip(grp->answer);
+	strcpy(group, grp->answer);
+    }
+    else
+	group[0] = '\0';
+
+    points_file = pfile->answer;
+
+    if (grp->answer == NULL && points_file == NULL)
+	G_fatal_error(_("Please select a group or give an input file."));
+    else if (grp->answer != NULL && points_file != NULL)
+	G_warning(_("Points in group will be ignored, GCPs in input file are used."));
+
+    order = atoi(val->answer);
+
+    if (order < 1 || order > MAXORDER)
+	G_fatal_error(_("Invalid order (%d); please enter 1 to %d"), order,
+		      MAXORDER);
+
+    Vect_set_open_level(1);
+    Vect_open_old2(&In, in_opt->answer, "", "");
+    
+    use3d = (Vect_is_3d(&In) && flag_use3d->answer);
+    
+    if (use3d && !points_file)
+	G_fatal_error(_("A file with 3D control points is needed for 3d transformation"));
+
+    if (print_rms->answer) {
+	if (sep->answer) {
+	    rms_sep = sep->answer;
+	    if (strcmp(rms_sep, "\\t") == 0)
+		rms_sep = "\t";
+	    if (strcmp(rms_sep, "tab") == 0)
+		rms_sep = "\t";
+	    if (strcmp(rms_sep, "space") == 0)
+		rms_sep = " ";
+	    if (strcmp(rms_sep, "comma") == 0)
+		rms_sep = ",";
+	}
+	else
+	    rms_sep = "|";
+    }
+    else
+	rms_sep = NULL;
+
+    /* read the control points for the group */
+    get_control_points(group, points_file, order, use3d,
+                       print_rms->answer, rms_sep);
+    
+    if (print_rms->answer) {
+	Vect_close(&In);
+	exit(EXIT_SUCCESS);
+    }
+
+    /* get the target */
+    get_target(group);
+
+    /* Check the GRASS_OVERWRITE environment variable */
+    if ((overstr = getenv("GRASS_OVERWRITE")))  /* OK ? */
+	target_overwrite = atoi(overstr);
+
+    if (!target_overwrite) {
+	/* check if output exists in target location/mapset */
+	
+	select_target_env();
+		
+	if (G_find_vector2(out_opt->answer, G_mapset())) {
+	    G_warning(_("The vector map <%s> already exists in"), out_opt->answer);
+	    G_warning(_("target LOCATION %s, MAPSET %s:"),
+		      G_location(), G_mapset());
+	    G_fatal_error(_("Rectification cancelled."));
+	}
+	select_current_env();
+    }
+    else
+	G_debug(1, "Overwriting OK");
+
+    select_target_env();
+    Vect_open_new(&Out, out_opt->answer, Vect_is_3d(&In));
+    Vect_copy_head_data(&In, &Out);
+    Vect_hist_copy(&In, &Out);
+    Vect_hist_command(&Out);
+    select_current_env();
+    
+    Points = Vect_new_line_struct();
+    OPoints = Vect_new_line_struct();
+    Cats = Vect_new_cats_struct();
+
+    /* count lines */
+    nlines = 0;
+    while (1) {
+	type = Vect_read_next_line(&In, Points, Cats);
+	if (type == 0)
+	    continue;		/* Dead */
+
+	if (type == -1)
+	    G_fatal_error(_("Reading input vector map"));
+	if (type == -2)
+	    break;
+
+	nlines++;
+    }
+	
+    Vect_rewind(&In);
+
+    i = 0;
+    z = 0.0;
+    while ((type = Vect_read_next_line(&In, Points, Cats)) > 0) {
+	G_percent(i++, nlines, 4);
+	
+	Vect_reset_line(OPoints);
+	for (n = 0; n < Vect_get_num_line_points(Points); n++) {
+	    if (use3d) {
+		CRS_georef_3d(Points->x[n], Points->y[n], Points->z[n],
+			      &x, &y, &z, E12, N12, Z12, order);
+	    }
+	    else {
+		CRS_georef(Points->x[n], Points->y[n],
+			      &x, &y, E12, N12, order);
+		z = Points->z[n];
+	    }
+	    Vect_append_point(OPoints, x, y, z);
+	}
+	select_target_env();
+	Vect_write_line(&Out, type, OPoints, Cats);
+	select_current_env();
+    }
+    G_percent(1, 1, 1);
+    
+    select_target_env();
+    if (!no_topo->answer)
+	Vect_build(&Out);
+    /* Copy tables */
+    if (Vect_copy_tables(&In, &Out, 0))
+        G_warning(_("Failed to copy attribute table to output map"));
+    Vect_close(&Out);
+    
+    select_current_env();
+
+    Vect_close(&In);
+
+    G_done_msg(" ");
+
+    exit(EXIT_SUCCESS);
+}

vector/v.rectify/target.c

@@ -0,0 +1,44 @@
+#include <string.h>
+#include <unistd.h>
+#include <grass/glocale.h>
+#include "global.h"
+
+int get_target(char *group)
+{
+    char location[GMAPSET_MAX];
+    char mapset[GMAPSET_MAX];
+    char buf[1024];
+    int stat;
+
+    if (group && *group) {
+	if (!I_get_target(group, location, mapset)) {
+	    sprintf(buf, _("Target information for group <%s> missing"), group);
+	    goto error;
+	}
+    }
+    else {
+	sprintf(location, "%s", G_location());
+	sprintf(mapset, "%s", G_mapset());
+    }
+
+    sprintf(buf, "%s/%s", G_gisdbase(), location);
+    if (access(buf, 0) != 0) {
+	sprintf(buf, _("Target location <%s> not found"), location);
+	goto error;
+    }
+    select_target_env();
+    G__setenv("LOCATION_NAME", location);
+    stat = G__mapset_permissions(mapset);
+    if (stat > 0) {
+	G__setenv("MAPSET", mapset);
+	select_current_env();
+	return 1;
+    }
+    sprintf(buf, _("Mapset <%s> in target location <%s> - "), mapset, location);
+    strcat(buf, stat == 0 ? _("permission denied") : _("not found"));
+  error:
+    strcat(buf, _("Please run i.target for group."));
+    strcat(buf, group);
+    G_fatal_error(buf);
+    return 1;			/* never reached */
+}

vector/v.rectify/v.rectify.html

@@ -0,0 +1,115 @@
+<h2>DESCRIPTION</h2>
+
+<em>v.rectify</em> uses control points to calculate a 2D or 3D 
+transformation matrix based on a first, second, or third order 
+polynomial and then converts x,y(, z) coordinates to standard map 
+coordinates for each object in the vector map. The result is a vector 
+map with a transformed coordinate system (i.e., a different coordinate
+system than before it was rectified).
+
+<p>
+2D Ground Control Points can be identified in 
+<em><a href="i.points.html">i.points</a></em>
+or
+<em><a href="i.vpoints.html">i.vpoints</a></em>
+<p>
+3D Ground Control Points must be provided in a text file with the 
+<b>points</b> option. The 3D format is equivalent to the format for 2D 
+ground control points with an additional third coordinate:
+
+<div class="code"><pre>
+ x y z east north height status
+</pre></div>
+where x, y, z are source coordinates, east, north, height are target 
+coordinates and status (0 or 1) indicates whether a given point should 
+be used. 
+
+<p>
+If no <b>group</b> is given, the rectified vector will be written to 
+the current mapset. If a <b>group</b> is given and a target has been 
+set for this group with <em><a href="i.target.html">i.target</a></em>, 
+the rectified vector will be written to the target location and mapset.
+
+<h3>Coordinate transformation and RMSE</h3>
+<p>The desired order of transformation (1, 2, or 3) is selected with the
+<b>order</b> option.
+
+<em>v.rectify</em> will calculate the RMSE if the <b>-r</b> flag is 
+given and print out statistcs in tabular format. The last row gives a 
+summary with the first column holding the number of active points, 
+followed by average deviations for each dimension and both forward and 
+backward transformation and finally forward and backward overall RMSE.
+
+<h4>2D linear affine transformation (1st order transformation)</h4>
+
+<dl>
+	<dd> x' = a1 + b1 * x + c1 * y
+	<dd> y' = a2 + b2 * x + c2 * y
+</dl>
+
+<h4>3D linear affine transformation (1st order transformation)</h4>
+
+<dl>
+	<dd> x' = a1 + b1 * x + c1 * y + d1 * z
+	<dd> y' = a2 + b2 * x + c2 * y + d2 * z
+	<dd> z' = a3 + b3 * x + c3 * y + d3 * z
+</dl>
+
+The a,b,c,d coefficients are determined by least squares regression
+based on the control points entered.  This transformation
+applies scaling, translation and rotation. It is NOT a
+general purpose rubber-sheeting, nor is it ortho-photo
+rectification using a DEM, not second order polynomial,
+etc. It can be used if (1) you have geometrically correct
+data, and (2) the terrain or camera distortion effect can
+be ignored.
+
+<h4>Polynomial Transformation Matrix (2nd, 3d order transformation)</h4>
+
+<em>v.rectify</em> uses a first, second, or third order transformation
+matrix to calculate the registration coefficients. The number
+of control points required for a 2D transformation of the selected order
+(represented by n) is
+
+<dl>
+<dd>((n + 1) * (n + 2) / 2) 
+</dl>
+
+or 3, 6, and 10 respectively. For a 3D transformation of first, second, 
+or third order, the number of required control points is 4, 10, and 20 
+respectively. It is strongly recommended that one or more additional 
+points be identified to allow for an overly-determined transformation 
+calculation which will generate the Root Mean Square (RMS) error values 
+for each included point. The polynomial equations are determined using 
+a modified Gaussian elimination method.
+
+
+<h2>SEE ALSO</h2>
+
+The GRASS 4 <em>
+<a href="http://grass.osgeo.org/gdp/imagery/grass4_image_processing.pdf">Image
+Processing manual</a></em>
+
+<p><em>
+  <a href="m.transform.html">m.transform</a>,
+  <a href="i.rectify.html">i.rectify</a>,
+  <a href="r.proj.html">r.proj</a>,
+  <a href="v.proj.html">v.proj</a>,
+  <a href="i.group.html">i.group</a>,
+  <a href="i.points.html">i.points</a>,
+  <a href="i.vpoints.html">i.vpoints</a>,
+  <a href="i.target.html">i.target</a>
+  <br>
+  <a href="wxGUI.GCP_Manager.html">Manage Ground Control Points</a>
+</em>
+
+
+<h2>AUTHORS</h2>
+
+Markus Metz
+
+<p>
+based on <a href="i.rectify.html">i.rectify</a>
+
+
+<p><i>Last changed: $Date$</i>