[GRASS-SVN] r49825 - in grass/trunk/vector: . v.rectify
svn_grass at osgeo.org
svn_grass at osgeo.org
Mon Dec 19 04:53:35 EST 2011
Author: mmetz
Date: 2011-12-19 01:53:35 -0800 (Mon, 19 Dec 2011)
New Revision: 49825
Added:
grass/trunk/vector/v.rectify/
grass/trunk/vector/v.rectify/crs3d.c
grass/trunk/vector/v.rectify/v.rectify.html
Removed:
grass/trunk/vector/v.rectify/README
grass/trunk/vector/v.rectify/bilinear.c
grass/trunk/vector/v.rectify/bilinear_f.c
grass/trunk/vector/v.rectify/cubic.c
grass/trunk/vector/v.rectify/cubic_f.c
grass/trunk/vector/v.rectify/exec.c
grass/trunk/vector/v.rectify/get_wind.c
grass/trunk/vector/v.rectify/i.rectify.html
grass/trunk/vector/v.rectify/lanczos.c
grass/trunk/vector/v.rectify/nearest.c
grass/trunk/vector/v.rectify/readcell.c
grass/trunk/vector/v.rectify/rectify.c
grass/trunk/vector/v.rectify/report.c
Modified:
grass/trunk/vector/Makefile
grass/trunk/vector/v.rectify/Makefile
grass/trunk/vector/v.rectify/cp.c
grass/trunk/vector/v.rectify/crs.c
grass/trunk/vector/v.rectify/crs.h
grass/trunk/vector/v.rectify/global.h
grass/trunk/vector/v.rectify/main.c
grass/trunk/vector/v.rectify/target.c
Log:
add new module v.rectify
Modified: grass/trunk/vector/Makefile
===================================================================
--- grass/trunk/vector/Makefile 2011-12-19 06:24:44 UTC (rev 49824)
+++ grass/trunk/vector/Makefile 2011-12-19 09:53:35 UTC (rev 49825)
@@ -70,6 +70,7 @@
v.qcount \
v.random \
v.reclass \
+ v.rectify \
v.sample \
v.segment \
v.select \
Modified: grass/trunk/vector/v.rectify/Makefile
===================================================================
--- grass/trunk/imagery/i.rectify/Makefile 2011-12-16 13:42:01 UTC (rev 49787)
+++ grass/trunk/vector/v.rectify/Makefile 2011-12-19 09:53:35 UTC (rev 49825)
@@ -1,9 +1,9 @@
MODULE_TOPDIR = ../..
-PGM = i.rectify
+PGM = v.rectify
-LIBES = $(IMAGERYLIB) $(RASTERLIB) $(GISLIB)
-DEPENDENCIES = $(IMAGERYDEP) $(RASTERDEP) $(GISDEP)
+LIBES = $(IMAGERYLIB) $(VECTORLIB) $(GISLIB)
+DEPENDENCIES = $(IMAGERYDEP) $(VECTORDEP) $(GISDEP)
include $(MODULE_TOPDIR)/include/Make/Module.make
Deleted: grass/trunk/vector/v.rectify/README
===================================================================
--- grass/trunk/imagery/i.rectify/README 2011-12-16 13:42:01 UTC (rev 49787)
+++ grass/trunk/vector/v.rectify/README 2011-12-19 09:53:35 UTC (rev 49825)
@@ -1,3 +0,0 @@
-As i.rectify functionality is integrated,
-i.rectify2 is renamed to i.rectify
-
Deleted: grass/trunk/vector/v.rectify/bilinear.c
===================================================================
--- grass/trunk/imagery/i.rectify/bilinear.c 2011-12-16 13:42:01 UTC (rev 49787)
+++ grass/trunk/vector/v.rectify/bilinear.c 2011-12-19 09:53:35 UTC (rev 49825)
@@ -1,60 +0,0 @@
-/*
- * Name
- * bilinear.c -- use bilinear interpolation for given row, col
- *
- * Description
- * bilinear interpolation for the given row, column indices.
- * If the given row or column is outside the bounds of the input map,
- * the point in the output map is set to NULL.
- * If any of the 4 surrounding points to be used in the interpolation
- * is NULL it is filled with is neighbor value
- */
-
-#include <math.h>
-#include <grass/gis.h>
-#include <grass/raster.h>
-#include "global.h"
-
-void p_bilinear(struct cache *ibuffer, /* input buffer */
- void *obufptr, /* ptr in output buffer */
- int cell_type, /* raster map type of obufptr */
- double *row_idx, /* row index */
- double *col_idx, /* column index */
- struct Cell_head *cellhd /* information of output map */
- )
-{
- int row; /* row indices for interp */
- int col; /* column indices for interp */
- int i, j;
- DCELL t, u; /* intermediate slope */
- DCELL result; /* result of interpolation */
- DCELL c[2][2];
-
- /* cut indices to integer */
- row = (int)floor(*row_idx - 0.5);
- col = (int)floor(*col_idx - 0.5);
-
- /* check for out of bounds - if out of bounds set NULL value and return */
- if (row < 0 || row + 1 >= cellhd->rows || col < 0 || col + 1 >= cellhd->cols) {
- Rast_set_null_value(obufptr, 1, cell_type);
- return;
- }
-
- for (i = 0; i < 2; i++)
- for (j = 0; j < 2; j++) {
- const DCELL *cellp = CPTR(ibuffer, row + i, col + j);
- if (Rast_is_d_null_value(cellp)) {
- Rast_set_null_value(obufptr, 1, cell_type);
- return;
- }
- c[i][j] = *cellp;
- }
-
- /* do the interpolation */
- t = *col_idx - 0.5 - col;
- u = *row_idx - 0.5 - row;
-
- result = Rast_interp_bilinear(t, u, c[0][0], c[0][1], c[1][0], c[1][1]);
-
- Rast_set_d_value(obufptr, result, cell_type);
-}
Deleted: grass/trunk/vector/v.rectify/bilinear_f.c
===================================================================
--- grass/trunk/imagery/i.rectify/bilinear_f.c 2011-12-16 13:42:01 UTC (rev 49787)
+++ grass/trunk/vector/v.rectify/bilinear_f.c 2011-12-19 09:53:35 UTC (rev 49825)
@@ -1,50 +0,0 @@
-/*
- * Name
- * bilinear_f.c -- use bilinear interpolation with fallback for given row, col
- *
- * Description
- * bilinear interpolation for the given row, column indices.
- * If the interpolated value (but not the nearest) is null,
- * it falls back to nearest neighbor.
- */
-
-#include <grass/gis.h>
-#include <grass/raster.h>
-#include <math.h>
-#include "global.h"
-
-void p_bilinear_f(struct cache *ibuffer, /* input buffer */
- void *obufptr, /* ptr in output buffer */
- int cell_type, /* raster map type of obufptr */
- double *row_idx, /* row index */
- double *col_idx, /* column index */
- struct Cell_head *cellhd /* cell header of input layer */
- )
-{
- /* start nearest neighbor to do some basic tests */
- int row, col; /* row/col of nearest neighbor */
- DCELL *cellp, cell;
-
- /* cut indices to integer */
- row = (int)floor(*row_idx);
- col = (int)floor(*col_idx);
-
- /* check for out of bounds - if out of bounds set NULL value */
- if (row < 0 || row >= cellhd->rows || col < 0 || col >= cellhd->cols) {
- Rast_set_null_value(obufptr, 1, cell_type);
- return;
- }
-
- cellp = CPTR(ibuffer, row, col);
- /* if nearest is null, all the other interps will be null */
- if (Rast_is_d_null_value(cellp)) {
- Rast_set_null_value(obufptr, 1, cell_type);
- return;
- }
- cell = *cellp;
-
- p_bilinear(ibuffer, obufptr, cell_type, row_idx, col_idx, cellhd);
- /* fallback to nearest if bilinear is null */
- if (Rast_is_d_null_value(obufptr))
- Rast_set_d_value(obufptr, cell, cell_type);
-}
Modified: grass/trunk/vector/v.rectify/cp.c
===================================================================
--- grass/trunk/imagery/i.rectify/cp.c 2011-12-16 13:42:01 UTC (rev 49787)
+++ grass/trunk/vector/v.rectify/cp.c 2011-12-19 09:53:35 UTC (rev 49825)
@@ -1,25 +1,340 @@
#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" /* CRS HEADER FILE */
-int get_control_points(char *group, int order /* THIS HAS BEEN ADDED WITH THE CRS MODIFICATIONS */
- )
+#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 (!I_get_control_points(group, &cp))
- exit(0);
+ 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");
- sprintf(msg, _("Control Point file for group <%s@%s> - "),
- group, G_mapset());
+ 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;
+ }
- switch (CRS_compute_georef_equations(&cp, E12, N12, E21, N21, order)) {
+ 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 + 1) * (order + 2) / 2);
+ order_pnts[use3d != 0][order]);
break;
case -1:
strcat(msg, _("Poorly placed control points."));
@@ -32,15 +347,15 @@
strcat(msg, _("Invalid order"));
break;
default:
- /* COMMENTED OUT WHEN SUPPORT FOR 3rd ORDER WAS ADDED BY 'CRS'
- E12a = E12[0]; E12b = E12[1]; E12c = E12[2];
- N12a = N12[0]; N12b = N12[1]; N12c = N12[2];
- E21a = E21[0]; E21b = E21[1]; E21c = E21[2];
- N21a = N21[0]; N21b = N21[1]; N21c = N21[2];
- */
- return 1;
+ break;
}
- G_fatal_error(msg);
+ if (ret != 1)
+ G_fatal_error(msg);
+
+ if (rms) {
+ compute_rms(&cp, &cp3, order, use3d, sep);
+ }
+
- return 0; /* G_fatal_error() calls exit() */
+ return 1;
}
Modified: grass/trunk/vector/v.rectify/crs.c
===================================================================
--- grass/trunk/imagery/i.rectify/crs.c 2011-12-16 13:42:01 UTC (rev 49787)
+++ grass/trunk/vector/v.rectify/crs.c 2011-12-19 09:53:35 UTC (rev 49825)
@@ -1,8 +1,6 @@
-/***************************************************************************/
+/***********************************************************************
-/***************************************************************************/
-/*
CRS.C - Center for Remote Sensing rectification routines
Written By: Brian J. Buckley
@@ -20,45 +18,18 @@
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>
-/*
- #define MSDOS 1
- */
-
-/*
- #define BDEBUG
- */
-#ifdef BDEBUG
-FILE *fp;
-#endif
-
-#ifdef MSDOS
-
-#include "stdlib.h"
-#include "dummy.h"
-
-typedef double DOUBLE;
-
-#else
-
+#include <grass/gis.h>
#include <grass/imagery.h>
-typedef double DOUBLE;
-
-#endif
-
#include "crs.h"
/* STRUCTURE FOR USE INTERNALLY WITH THESE FUNCTIONS. THESE FUNCTIONS EXPECT
@@ -67,19 +38,13 @@
struct MATRIX
{
int n; /* SIZE OF THIS MATRIX (N x N) */
- DOUBLE *v;
+ 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 */
@@ -87,148 +52,45 @@
#define MPARMERR -3 /* PARAMETER ERROR */
#define MINTERR -4 /* INTERNAL ERROR */
-/***************************************************************************/
-/*
- FUNCTION PROTOTYPES FOR STATIC (INTERNAL) FUNCTIONS
- */
+/***********************************************************************
-/***************************************************************************/
+ FUNCTION PROTOTYPES FOR STATIC (INTERNAL) FUNCTIONS
-#ifdef MSDOS
+************************************************************************/
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);
+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);
-#ifdef BDEBUG
-static int checkgeoref(struct Control_Points *, double *, double *, int, int,
- FILE * fp);
-#endif
+/***********************************************************************
-#else
+ TRANSFORM A SINGLE COORDINATE PAIR.
-static int calccoef();
-static int calcls();
-static int exactdet();
-static int solvemat();
-static DOUBLE term();
+************************************************************************/
-#ifdef BDEBUG
-static int checkgeoref();
-#endif
-
-#endif
-
-/***************************************************************************/
-/*
- USE THIS TRANSFORMATION FUNCTION IF YOU WANT TO DO ARRAYS.
- */
-
-/***************************************************************************/
-
-#ifdef BETTERGEOREF
-
-extern void CRS_georef(int, DOUBLE *, DOUBLE *, DOUBLE *, DOUBLE *, int);
-
-void CRS_georef2(int order, /* ORDER OF TRANSFORMATION TO BE PERFORMED, MUST MATCH THE
- ORDER USED TO CALCULATE THE COEFFICIENTS */
- double E[], /* EASTING COEFFICIENTS */
- double N[], /* NORTHING COEFFICIENTS */
- double e[], /* EASTINGS TO BE TRANSFORMED */
- double n[], /* NORTHINGS TO BE TRANSFORMED */
- int numpts /* NUMBER OF POINTS TO BE TRANSFORMED */
- )
-{
- DOUBLE e3, e2n, e2, en2, en, e1, n3, n2, n1;
- int i;
-
- switch (order) {
- case 3:
-
- for (i = 0; i < numpts; i++) {
- e1 = e[i];
- n1 = n[i];
- e2 = e1 * e1;
- en = e1 * n1;
- n2 = n1 * n1;
- e3 = e1 * e2;
- e2n = e2 * n1;
- en2 = e1 * n2;
- n3 = n1 * n2;
-
- e[i] = 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[i] = 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;
-
- case 2:
-
- for (i = 0; i < numpts; i++) {
- e1 = e[i];
- n1 = n[i];
- e2 = e1 * e1;
- n2 = n1 * n1;
- en = e1 * n1;
-
- e[i] = E[0] + E[1] * e1 + E[2] * n1 +
- E[3] * e2 + E[4] * en + E[5] * n2;
- n[i] = N[0] + N[1] * e1 + N[2] * n1 +
- N[3] * e2 + N[4] * en + N[5] * n2;
- }
- break;
-
- case 1:
-
- for (i = 0; i < numpts; i++) {
- e1 = e[i];
- n1 = n[i];
- e[i] = E[0] + E[1] * e1 + E[2] * n1;
- n[i] = N[0] + N[1] * e1 + N[2] * n1;
- }
- break;
- }
-}
-
-#endif
-
-/***************************************************************************/
-/*
- 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 */
+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;
+ 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;
@@ -238,7 +100,6 @@
break;
case 3:
-
e2 = e1 * e1;
en = e1 * n1;
n2 = n1 * n1;
@@ -258,48 +119,36 @@
break;
default:
-
- return (MPARMERR);
- break;
+ return MPARMERR;
}
- return (MSUCCESS);
+ return MSUCCESS;
}
-/***************************************************************************/
-/*
- COMPUTE THE GEOREFFERENCING COEFFICIENTS BASED ON A SET OF CONTROL POINTS
- */
+/***********************************************************************
-/***************************************************************************/
+ 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)
+************************************************************************/
+
+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);
+ return MPARMERR;
-#ifdef BDEBUG
- fp = fopen("error.dat", "w");
- if (fp == NULL)
- return (-1);
-#endif
-
/* CALCULATE THE FORWARD TRANSFORMATION COEFFICIENTS */
status = calccoef(cp, E12, N12, order);
+
if (status != MSUCCESS)
- return (status);
+ return status;
-#ifdef BDEBUG
- checkgeoref(cp, E12, N12, order, 1, fp);
-#endif
-
/* SWITCH THE 1 AND 2 EASTING AND NORTHING ARRAYS */
tempptr = cp->e1;
@@ -313,11 +162,6 @@
status = calccoef(cp, E21, N21, order);
-#ifdef BDEBUG
- checkgeoref(cp, E21, N21, order, 0, fp);
- fclose(fp);
-#endif
-
/* SWITCH THE 1 AND 2 EASTING AND NORTHING ARRAYS BACK */
tempptr = cp->e1;
@@ -327,22 +171,22 @@
cp->n1 = cp->n2;
cp->n2 = tempptr;
- return (status);
+ return status;
}
-/***************************************************************************/
-/*
- COMPUTE THE GEOREFFERENCING COEFFICIENTS BASED ON A SET OF CONTROL POINTS
- */
+/***********************************************************************
-/***************************************************************************/
+ COMPUTE THE GEOREFFERENCING COEFFICIENTS
+ BASED ON A SET OF CONTROL POINTS
-static int
-calccoef(struct Control_Points *cp, double E[], double N[], int order)
+************************************************************************/
+
+static int calccoef(struct Control_Points *cp, double E[], double N[],
+ int order)
{
struct MATRIX m;
- DOUBLE *a;
- DOUBLE *b;
+ double *a;
+ double *b;
int numactive; /* NUMBER OF ACTIVE CONTROL POINTS */
int status, i;
@@ -359,47 +203,36 @@
m.n = ((order + 1) * (order + 2)) / 2;
if (numactive < m.n)
- return (MNPTERR);
+ return MNPTERR;
/* INITIALIZE MATRIX */
- m.v = (DOUBLE *) G_calloc(m.n * m.n, sizeof(DOUBLE));
- if (m.v == NULL) {
- return (MMEMERR);
- }
- a = (DOUBLE *) G_calloc(m.n, sizeof(DOUBLE));
- if (a == NULL) {
- G_free((char *)m.v);
- return (MMEMERR);
- }
- b = (DOUBLE *) G_calloc(m.n, sizeof(DOUBLE));
- if (b == NULL) {
- G_free((char *)m.v);
- G_free((char *)a);
- return (MMEMERR);
- }
+ 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((char *)m.v);
- G_free((char *)a);
- G_free((char *)b);
+ G_free(m.v);
+ G_free(a);
+ G_free(b);
- return (status);
+ return status;
}
-/***************************************************************************/
-/*
- CALCULATE THE TRANSFORMATION COEFFICIENTS WITH EXACTLY THE MINIMUM
- NUMBER OF CONTROL POINTS REQUIRED FOR THIS TRANSFORMATION.
- */
+/***********************************************************************
-/***************************************************************************/
+ 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 */
+************************************************************************/
+
+static int exactdet(struct Control_Points *cp, struct MATRIX *m,
+ double a[], double b[],
+ double E[], /* EASTING COEFFICIENTS */
double N[] /* NORTHING COEFFICIENTS */
)
{
@@ -410,56 +243,35 @@
if (cp->status[pntnow] > 0) {
/* POPULATE MATRIX M */
-#ifdef BDEBUG
- fprintf(fp, "%2d ", pntnow + 1);
-#endif
-
- for (j = 1; j <= m->n; j++) {
+ for (j = 1; j <= m->n; j++)
M(currow, j) = term(j, cp->e1[pntnow], cp->n1[pntnow]);
-#ifdef BDEBUG
- fprintf(fp, "%+14.7le ", M(currow, j));
- if (j == 5)
- fprintf(fp, "\n ");
-#endif
- }
-#ifdef BDEBUG
- fprintf(fp, "\n");
-#endif
/* POPULATE MATRIX A AND B */
a[currow - 1] = cp->e2[pntnow];
b[currow - 1] = cp->n2[pntnow];
-#ifdef BDEBUG
- fprintf(fp, " %+14.7le ", a[currow - 1]);
- fprintf(fp, "%+14.7le\n", b[currow - 1]);
-#endif
currow++;
}
-#ifdef BDEBUG
- else {
- fprintf(fp, "%2d UNUSED\n", pntnow + 1);
- }
-#endif
}
if (currow - 1 != m->n)
- return (MINTERR);
+ return MINTERR;
- return (solvemat(m, a, b, E, N));
+ 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.
- */
+/***********************************************************************
-/***************************************************************************/
+ 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 */
+************************************************************************/
+
+static int calcls(struct Control_Points *cp, struct MATRIX *m,
+ double a[], double b[],
+ double E[], /* EASTING COEFFICIENTS */
double N[] /* NORTHING COEFFICIENTS */
)
{
@@ -492,83 +304,81 @@
}
if (numactive <= m->n)
- return (MINTERR);
+ return MINTERR;
/* TRANSPOSE VALUES IN UPPER HALF OF M TO OTHER HALF */
- for (i = 2; i <= m->n; i++) {
+ 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));
+ 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
- */
+ 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)
+************************************************************************/
+
+static double term(int term, double e, double n)
{
switch (term) {
case 1:
- return ((DOUBLE) 1.0);
+ return 1.0;
case 2:
- return ((DOUBLE) e);
+ return e;
case 3:
- return ((DOUBLE) n);
+ return n;
case 4:
- return ((DOUBLE) (e * e));
+ return e * e;
case 5:
- return ((DOUBLE) (e * n));
+ return e * n;
case 6:
- return ((DOUBLE) (n * n));
+ return n * n;
case 7:
- return ((DOUBLE) (e * e * e));
+ return e * e * e;
case 8:
- return ((DOUBLE) (e * e * n));
+ return e * e * n;
case 9:
- return ((DOUBLE) (e * n * n));
+ return e * n * n;
case 10:
- return ((DOUBLE) (n * n * n));
+ return n * n * n;
}
- return ((DOUBLE) 0.0);
+
+ 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 |
+ SOLVE FOR THE 'E' AND 'N' COEFFICIENTS BY USING A SOMEWHAT MODIFIED
+ GAUSSIAN ELIMINATION METHOD.
- and
+ | 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 | | 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[])
+************************************************************************/
+
+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 */
+ double factor, temp;
+ double pivot; /* ACTUAL VALUE OF THE LARGEST PIVOT CANDIDATE */
for (i = 1; i <= m->n; i++) {
j = i;
@@ -590,7 +400,7 @@
/* co-linear points results in a solution with rounding error */
if (pivot == 0.0)
- return (MUNSOLVABLE);
+ return MUNSOLVABLE;
/* if row with highest pivot is not the current row, switch them */
@@ -632,80 +442,5 @@
N[i - 1] = b[i - 1] / M(i, i);
}
- return (MSUCCESS);
+ return MSUCCESS;
}
-
-/***************************************************************************/
-/*
- */
-
-/***************************************************************************/
-
-#ifdef BDEBUG
-
-static int checkgeoref(struct Control_Points *cp,
- double E[], double N[], int order, int forward,
- FILE * fp)
-{
- DOUBLE xrms, yrms, dx, dy, dx2, dy2, totaldist, dist;
- double tempx, tempy;
- int i, n, numactive;
-
- n = ((order + 1) * (order + 2)) / 2;
-
- if (forward)
- fprintf(fp, "FORWARD:\n");
- else
- fprintf(fp, "BACKWARD:\n");
-
- fprintf(fp, "%d order\n", order);
- for (i = 0; i < n; i++)
- fprintf(fp, "%+.17E %+.17E\n", E[i], N[i]);
-
- xrms = yrms = dx2 = dy2 = totaldist = 0.0;
- numactive = 0;
- for (i = 0; i < cp->count; i++) {
- fprintf(fp, "\nCONTROL POINT: %d\n", i + 1);
-
- fprintf(fp, "%20s: %+.20lE %+.20lE\n", "ORIGINAL POINT",
- cp->e1[i], cp->n1[i]);
- fprintf(fp, "%20s: %+.20lE %+.20lE\n", "DESIRED POINT",
- cp->e2[i], cp->n2[i]);
-
- if (cp->status[i] > 0) {
- numactive++;
- CRS_georef(cp->e1[i], cp->n1[i], &tempx, &tempy, E, N, order);
-
- fprintf(fp, "%20s: %+.20lE %+.20lE\n", "CALCULATED POINT", tempx,
- tempy);
- dx = tempx - cp->e2[i];
- dy = tempy - cp->n2[i];
- fprintf(fp, "%20s: %+.20lE %+.20lE\n", "RESIDUAL ERROR", dx, dy);
- dx2 = dx * dx;
- dy2 = dy * dy;
- dist = sqrt(dx2 + dy2);
- fprintf(fp, "%20s: %+.20lE\n", "DISTANCE (RMS) ERROR", dist);
-
- xrms += dx2;
- yrms += dy2;
-
- totaldist += dist;
- }
- else
- fprintf(fp, "NOT USED\n");
- }
- xrms = sqrt(xrms / (DOUBLE) numactive);
- yrms = sqrt(yrms / (DOUBLE) numactive);
-
- fprintf(fp, "\n%20s: %+.20lE %+.20lE\n", "RMS ERROR", xrms, yrms);
-
- fprintf(fp, "\n%20s: %+.20lE\n", "TOTAL RMS ERROR",
- sqrt(xrms * xrms + yrms * yrms));
-
- fprintf(fp, "\n%20s: %+.20lE\n", "AVG. DISTANCE ERROR",
- totaldist / numactive);
-
- return (0);
-}
-
-#endif
Modified: grass/trunk/vector/v.rectify/crs.h
===================================================================
--- grass/trunk/imagery/i.rectify/crs.h 2011-12-16 13:42:01 UTC (rev 49787)
+++ grass/trunk/vector/v.rectify/crs.h 2011-12-19 09:53:35 UTC (rev 49825)
@@ -24,24 +24,31 @@
#define MAXORDER 3
-#ifdef MSDOS
+struct Control_Points_3D
+{
+ int count;
+ double *e1;
+ double *n1;
+ double *z1;
+ double *e2;
+ double *n2;
+ double *z2;
+ int *status;
+};
-extern int CRS_compute_georef_equations(struct Control_Points *, double *,
+
+/* crs.c */
+int CRS_compute_georef_equations(struct Control_Points *, double *,
double *, double *, double *, int);
-extern int CRS_georef(double, double, double *, double *, double *, double *,
+int CRS_georef(double, double, double *, double *, double *, double *,
int);
-#else
-
-#ifdef NO_PROTO
-extern int CRS_compute_georef_equations();
-extern int CRS_georef();
-#else
-/* crs.c */
-void CRS_georef2(int, double[], double[], double[], double[], int);
-int CRS_georef(double, double, double *, double *, double[], double[], int);
-int CRS_compute_georef_equations(struct Control_Points *,
- double[], double[], double[], double[], int);
-#endif
-
-#endif
+/* 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);
Copied: grass/trunk/vector/v.rectify/crs3d.c (from rev 49787, grass/trunk/imagery/i.rectify/crs.c)
===================================================================
--- grass/trunk/vector/v.rectify/crs3d.c (rev 0)
+++ grass/trunk/vector/v.rectify/crs3d.c 2011-12-19 09:53:35 UTC (rev 49825)
@@ -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;
+}
Deleted: grass/trunk/vector/v.rectify/cubic.c
===================================================================
--- grass/trunk/imagery/i.rectify/cubic.c 2011-12-16 13:42:01 UTC (rev 49787)
+++ grass/trunk/vector/v.rectify/cubic.c 2011-12-19 09:53:35 UTC (rev 49825)
@@ -1,68 +0,0 @@
-/*
- * Name
- * cubic.c -- use cubic convolution interpolation for given row, col
- *
- * Description
- * cubic returns the value in the buffer that is the result of cubic
- * convolution interpolation for the given row, column indices.
- * If the given row or column is outside the bounds of the input map,
- * the corresponding point in the output map is set to NULL.
- *
- * If single surrounding points in the interpolation matrix are
- * NULL they where filled with their neighbor
- */
-
-#include <grass/gis.h>
-#include <grass/raster.h>
-#include <math.h>
-#include "global.h"
-
-void p_cubic(struct cache *ibuffer, /* input buffer */
- void *obufptr, /* ptr in output buffer */
- int cell_type, /* raster map type of obufptr */
- double *row_idx, /* row index (decimal) */
- double *col_idx, /* column index (decimal) */
- struct Cell_head *cellhd /* information of output map */
- )
-{
- int row; /* row indices for interp */
- int col; /* column indices for interp */
- int i, j;
- DCELL t, u; /* intermediate slope */
- DCELL result; /* result of interpolation */
- DCELL val[4]; /* buffer for temporary values */
- DCELL cell[4][4];
-
- /* cut indices to integer */
- row = (int)floor(*row_idx - 0.5);
- col = (int)floor(*col_idx - 0.5);
-
- /* check for out of bounds of map - if out of bounds set NULL value */
- if (row - 1 < 0 || row + 2 >= cellhd->rows ||
- col - 1 < 0 || col + 2 >= cellhd->cols) {
- Rast_set_null_value(obufptr, 1, cell_type);
- return;
- }
-
- for (i = 0; i < 4; i++)
- for (j = 0; j < 4; j++) {
- const DCELL *cellp = CPTR(ibuffer, row - 1 + i, col - 1 + j);
- if (Rast_is_d_null_value(cellp)) {
- Rast_set_null_value(obufptr, 1, cell_type);
- return;
- }
- cell[i][j] = *cellp;
- }
-
- /* do the interpolation */
- t = *col_idx - 0.5 - col;
- u = *row_idx - 0.5 - row;
-
- for (i = 0; i < 4; i++) {
- val[i] = Rast_interp_cubic(t, cell[i][0], cell[i][1], cell[i][2], cell[i][3]);
- }
-
- result = Rast_interp_cubic(u, val[0], val[1], val[2], val[3]);
-
- Rast_set_d_value(obufptr, result, cell_type);
-}
Deleted: grass/trunk/vector/v.rectify/cubic_f.c
===================================================================
--- grass/trunk/imagery/i.rectify/cubic_f.c 2011-12-16 13:42:01 UTC (rev 49787)
+++ grass/trunk/vector/v.rectify/cubic_f.c 2011-12-19 09:53:35 UTC (rev 49825)
@@ -1,55 +0,0 @@
-/*
- * Name
- * cubic_f.c -- use cubic interpolation with fallback for given row, col
- *
- * Description
- * cubic interpolation for the given row, column indices.
- * If the interpolated value (but not the nearest) is null,
- * it falls back to bilinear. If that interp value is null,
- * it falls back to nearest neighbor.
- */
-
-#include <grass/gis.h>
-#include <grass/raster.h>
-#include <math.h>
-#include "global.h"
-
-void p_cubic_f(struct cache *ibuffer, /* input buffer */
- void *obufptr, /* ptr in output buffer */
- int cell_type, /* raster map type of obufptr */
- double *row_idx, /* row index */
- double *col_idx, /* column index */
- struct Cell_head *cellhd /* cell header of input layer */
- )
-{
- /* start nearest neighbor to do some basic tests */
- int row, col; /* row/col of nearest neighbor */
- DCELL *cellp, cell;
-
- /* cut indices to integer */
- row = (int)floor(*row_idx);
- col = (int)floor(*col_idx);
-
- /* check for out of bounds - if out of bounds set NULL value */
- if (row < 0 || row >= cellhd->rows || col < 0 || col >= cellhd->cols) {
- Rast_set_null_value(obufptr, 1, cell_type);
- return;
- }
-
- cellp = CPTR(ibuffer, row, col);
- /* if nearest is null, all the other interps will be null */
- if (Rast_is_d_null_value(cellp)) {
- Rast_set_null_value(obufptr, 1, cell_type);
- return;
- }
- cell = *cellp;
-
- p_cubic(ibuffer, obufptr, cell_type, row_idx, col_idx, cellhd);
- /* fallback to bilinear if cubic is null */
- if (Rast_is_d_null_value(obufptr)) {
- p_bilinear(ibuffer, obufptr, cell_type, row_idx, col_idx, cellhd);
- /* fallback to nearest if bilinear is null */
- if (Rast_is_d_null_value(obufptr))
- Rast_set_d_value(obufptr, cell, cell_type);
- }
-}
Deleted: grass/trunk/vector/v.rectify/exec.c
===================================================================
--- grass/trunk/imagery/i.rectify/exec.c 2011-12-16 13:42:01 UTC (rev 49787)
+++ grass/trunk/vector/v.rectify/exec.c 2011-12-19 09:53:35 UTC (rev 49825)
@@ -1,91 +0,0 @@
-/* CMD mode from Bob Covill 2001
- *
- * small fixes: MN
- *
- * Bug left: extension overwrites input name 1/2002
- */
-
-#include <stdlib.h>
-#include <stdio.h>
-#include <string.h>
-#include <sys/types.h>
-#include <sys/stat.h>
-#include <fcntl.h>
-#include <time.h>
-#include <unistd.h>
-#include <math.h>
-
-#include <grass/raster.h>
-#include <grass/glocale.h>
-
-#include "global.h"
-
-int exec_rectify(int order, char *extension, char *interp_method)
-/* ADDED WITH CRS MODIFICATIONS */
-{
- char *name;
- char *mapset;
- char *result;
- char *type = "raster";
- int n;
- struct Colors colr;
- struct Categories cats;
- struct History hist;
- int colr_ok, cats_ok;
- long start_time, rectify_time;
-
- Rast_set_output_window(&target_window);
- G_message("-----------------------------------------------");
-
- /* rectify each file */
- for (n = 0; n < ref.nfiles; n++) {
- if (ref_list[n]) {
- name = ref.file[n].name;
- mapset = ref.file[n].mapset;
-
- /* generate out name, add extension to output */
- result =
- G_malloc(strlen(ref.file[n].name) + strlen(extension) + 1);
- strcpy(result, ref.file[n].name);
- strcat(result, extension);
-
- select_current_env();
-
- cats_ok = Rast_read_cats(name, mapset, &cats) >= 0;
- colr_ok = Rast_read_colors(name, mapset, &colr) > 0;
-
- /* Initialze History */
- if (Rast_read_history(name, mapset, &hist) < 0)
- Rast_short_history(result, type, &hist);
-
- time(&start_time);
-
- if (rectify(name, mapset, result, order, interp_method)) {
- select_target_env();
-
- if (cats_ok) {
- Rast_write_cats(result, &cats);
- Rast_free_cats(&cats);
- }
- if (colr_ok) {
- Rast_write_colors(result, G_mapset(), &colr);
- Rast_free_colors(&colr);
- }
-
- /* Write out History */
- Rast_command_history(&hist);
- Rast_write_history(result, &hist);
-
- select_current_env();
- time(&rectify_time);
- report(rectify_time - start_time, 1);
- }
- else
- report((long)0, 0);
-
- G_free(result);
- }
- }
-
- return 0;
-}
Deleted: grass/trunk/vector/v.rectify/get_wind.c
===================================================================
--- grass/trunk/imagery/i.rectify/get_wind.c 2011-12-16 13:42:01 UTC (rev 49787)
+++ grass/trunk/vector/v.rectify/get_wind.c 2011-12-19 09:53:35 UTC (rev 49825)
@@ -1,118 +0,0 @@
-#include <math.h>
-#include <grass/glocale.h>
-#include "global.h"
-#include "crs.h" /* CRS HEADER FILE */
-
-int georef_window(struct Cell_head *w1, struct Cell_head *w2, int order, double res)
-{
- double n, e, ad;
- struct _corner {
- double n, e;
- } nw, ne, se, sw;
-
- /* extends */
- CRS_georef(w1->west, w1->north, &e, &n, E12, N12, order);
- w2->north = w2->south = n;
- w2->west = w2->east = e;
- nw.n = n;
- nw.e = e;
-
- CRS_georef(w1->east, w1->north, &e, &n, E12, N12, order);
- ne.n = n;
- ne.e = e;
- if (n > w2->north)
- w2->north = n;
- if (n < w2->south)
- w2->south = n;
- if (e > w2->east)
- w2->east = e;
- if (e < w2->west)
- w2->west = e;
-
- CRS_georef(w1->west, w1->south, &e, &n, E12, N12, order);
- sw.n = n;
- sw.e = e;
- if (n > w2->north)
- w2->north = n;
- if (n < w2->south)
- w2->south = n;
- if (e > w2->east)
- w2->east = e;
- if (e < w2->west)
- w2->west = e;
-
- CRS_georef(w1->east, w1->south, &e, &n, E12, N12, order);
- se.n = n;
- se.e = e;
- if (n > w2->north)
- w2->north = n;
- if (n < w2->south)
- w2->south = n;
- if (e > w2->east)
- w2->east = e;
- if (e < w2->west)
- w2->west = e;
-
- /* resolution */
- if (res > 0)
- w2->ew_res = w2->ns_res = res;
- else {
- /* this results in ugly res values, and ns_res != ew_res */
- /* and is no good for rotation */
- /*
- w2->ns_res = (w2->north - w2->south) / w1->rows;
- w2->ew_res = (w2->east - w2->west) / w1->cols;
- */
-
- /* alternative: account for rotation and order > 1 */
-
- /* N-S extends along western and eastern edge */
- w2->ns_res = (sqrt((nw.n - sw.n) * (nw.n - sw.n) +
- (nw.e - sw.e) * (nw.e - sw.e)) +
- sqrt((ne.n - se.n) * (ne.n - se.n) +
- (ne.e - se.e) * (ne.e - se.e))) / (2.0 * w1->rows);
-
- /* E-W extends along northern and southern edge */
- w2->ew_res = (sqrt((nw.n - ne.n) * (nw.n - ne.n) +
- (nw.e - ne.e) * (nw.e - ne.e)) +
- sqrt((sw.n - se.n) * (sw.n - se.n) +
- (sw.e - se.e) * (sw.e - se.e))) / (2.0 * w1->cols);
-
- /* make ew_res = ns_res */
- w2->ns_res = (w2->ns_res + w2->ew_res) / 2.0;
- w2->ew_res = w2->ns_res;
-
- /* nice round values */
- if (w2->ns_res > 1) {
- if (w2->ns_res < 10) {
- /* round to first decimal */
- w2->ns_res = (int)(w2->ns_res * 10 + 0.5) / 10.0;
- w2->ew_res = w2->ns_res;
- }
- else {
- /* round to integer */
- w2->ns_res = (int)(w2->ns_res + 0.5);
- w2->ew_res = w2->ns_res;
- }
- }
- }
-
- /* adjust extends */
- ad = w2->north > 0 ? 0.5 : -0.5;
- w2->north = (int) (ceil(w2->north / w2->ns_res) + ad) * w2->ns_res;
- ad = w2->south > 0 ? 0.5 : -0.5;
- w2->south = (int) (floor(w2->south / w2->ns_res) + ad) * w2->ns_res;
- ad = w2->east > 0 ? 0.5 : -0.5;
- w2->east = (int) (ceil(w2->east / w2->ew_res) + ad) * w2->ew_res;
- ad = w2->west > 0 ? 0.5 : -0.5;
- w2->west = (int) (floor(w2->west / w2->ew_res) + ad) * w2->ew_res;
-
- w2->rows = (w2->north - w2->south + w2->ns_res / 2.0) / w2->ns_res;
- w2->cols = (w2->east - w2->west + w2->ew_res / 2.0) / w2->ew_res;
-
- G_message(_("Region N=%f S=%f E=%f W=%f"), w2->north, w2->south,
- w2->east, w2->west);
- G_message(_("Resolution EW=%f NS=%f"), w2->ew_res, w2->ns_res);
-
- return 0;
-}
Modified: grass/trunk/vector/v.rectify/global.h
===================================================================
--- grass/trunk/imagery/i.rectify/global.h 2011-12-16 13:42:01 UTC (rev 49787)
+++ grass/trunk/vector/v.rectify/global.h 2011-12-19 09:53:35 UTC (rev 49825)
@@ -1,110 +1,19 @@
-/* These next defines determine the size of the sub-window that will
- * be held in memory. Larger values will require
- * more memory (but less i/o) If you increase these values, keep in
- * mind that although you think the i/o will decrease, system paging
- * (which goes on behind the scenes) may actual increase the i/o.
- */
#include <grass/gis.h>
#include <grass/imagery.h>
-#define L2BDIM 6
-#define BDIM (1<<(L2BDIM))
-#define L2BSIZE (2*(L2BDIM))
-#define BSIZE (1<<(L2BSIZE))
-#define HI(i) ((i)>>(L2BDIM))
-#define LO(i) ((i)&((BDIM)-1))
-typedef DCELL block[BDIM][BDIM];
-
-struct cache
-{
- int fd;
- int stride;
- int nblocks;
- block **grid;
- block *blocks;
- int *refs;
-};
-
-extern char *seg_mb;
-
-extern RASTER_MAP_TYPE map_type;
-extern int *ref_list;
-extern struct Ref ref;
-
-typedef void (*func) (struct cache *, void *, int, double *, double *, struct Cell_head *);
-
-extern func interpolate; /* interpolation routine */
-
-struct menu
-{
- func method; /* routine to interpolate new value */
- char *name; /* method name */
- char *text; /* menu display - full description */
-};
-
/* georef coefficients */
-extern double E12[10], N12[10];
-extern double E21[10], N21[10];
+extern double E12[20], N12[20], Z12[20];
+extern double E21[20], N21[20], Z21[20];
-/* DELETED WITH CRS MODIFICATIONS
- extern double E12a, E12b, E12c, N12a, N12b, N12c;
- extern double E21a, E21b, E21c, N21a, N21b, N21c;
- */
-extern struct Cell_head target_window;
-
/* cp.c */
-int get_control_points(char *, int);
+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);
-/* exec.c */
-int exec_rectify(int, char *, char *);
-
-/* get_wind.c */
-int get_target_window(int);
-int georef_window(struct Cell_head *, struct Cell_head *, int, double);
-
-/* rectify.c */
-int rectify(char *, char *, char *, int, char *);
-
-/* readcell.c */
-extern struct cache *readcell(int, const char *);
-extern block *get_block(struct cache *, int);
-
-#define BKIDX(c,y,x) ((y) * (c)->stride + (x))
-#define BKPTR(c,y,x) ((c)->grid[BKIDX((c),(y),(x))])
-#define BLOCK(c,y,x) (BKPTR((c),(y),(x)) ? BKPTR((c),(y),(x)) : get_block((c),BKIDX((c),(y),(x))))
-#define CPTR(c,y,x) (&(*BLOCK((c),HI((y)),HI((x))))[LO((y))][LO((x))])
-
-/* report.c */
-int report(long, int);
-
/* target.c */
int get_target(char *);
-
-/* declare resampling methods */
-/* bilinear.c */
-extern void p_bilinear(struct cache *, void *, int, double *, double *,
- struct Cell_head *);
-/* cubic.c */
-extern void p_cubic(struct cache *, void *, int, double *, double *,
- struct Cell_head *);
-/* nearest.c */
-extern void p_nearest(struct cache *, void *, int, double *, double *,
- struct Cell_head *);
-/* bilinear_f.c */
-extern void p_bilinear_f(struct cache *, void *, int, double *, double *,
- struct Cell_head *);
-/* cubic_f.c */
-extern void p_cubic_f(struct cache *, void *, int, double *, double *,
- struct Cell_head *);
-/* lanczos.c */
-extern void p_lanczos(struct cache *, void *, int, double *, double *,
- struct Cell_head *);
-extern void p_lanczos_f(struct cache *, void *, int, double *, double *,
- struct Cell_head *);
Deleted: grass/trunk/vector/v.rectify/i.rectify.html
===================================================================
--- grass/trunk/imagery/i.rectify/i.rectify.html 2011-12-16 13:42:01 UTC (rev 49787)
+++ grass/trunk/vector/v.rectify/i.rectify.html 2011-12-19 09:53:35 UTC (rev 49825)
@@ -1,180 +0,0 @@
-<h2>DESCRIPTION</h2>
-
-<em>i.rectify</em> uses the control
-points identified in
-<em><a href="i.points.html">i.points</a></em>
-or
-<em><a href="i.vpoints.html">i.vpoints</a></em>
-
-to calculate a transformation matrix based on a first,
-second, or third order polynomial and then converts x,y
-cell coordinates to standard map coordinates for each pixel
-in the image. The result is a planimetric image with a
-transformed coordinate system (i.e., a different coordinate
-system than before it was rectified).
-
-<p>
-<em><a href="i.points.html">i.points</a></em>
-or
-<em><a href="i.vpoints.html">i.vpoints</a></em>
-
-must be run before <em>i.rectify</em>, and both programs
-are required to rectify an image. An image must be
-rectified before it can reside in a standard coordinate
-LOCATION, and therefore be analyzed with the other map
-layers in the standard coordinate LOCATION. Upon
-completion of <em>i.rectify</em>, the rectified image is
-deposited in the target standard coordinate LOCATION. This
-LOCATION is selected using
-
-<em><a href="i.target.html">i.target</a></em>.
-
-<p>More than one raster map may be rectified at a time. Each cell file
-should be given a unique output file name. The rectified image or
-rectified raster maps will be located in the target LOCATION when the
-program is completed. The original unrectified files are not modified
-or removed.
-
-<p>If the <b>-c</b> flag is used, <em>i.rectify</em> will only rectify that
-portion of the image or raster map that occurs within the chosen window
-region in the target location, and only that portion of the cell
-file will be relocated in the target database. It is
-important therefore, to check the current mapset window in
-the target LOCATION if the <b>-c</b> flag is used.
-
-<p>
-If you are rectifying a file with plans to patch it to
-another file using the GRASS program <em>r.patch</em>,
-choose option number one, the current window in the target
-location. This window, however, must be the default window
-for the target LOCATION. When a file being rectified is
-smaller than the default window in which it is being
-rectified, NULLs are added to the rectified file. Patching
-files of the same size that contain NULL data,
-eliminates the possibility of a no-data line in the patched
-result. This is because, when the images are patched, the
-NULLs in the image are "covered" with non-NULL pixel
-values. When rectifying files that are going to be
-patched, rectify all of the files using the same default
-window.
-
-<h3>Coordinate transformation</h3>
-<p>The desired order of transformation (1, 2, or 3) is selected with the
-<b>order</b> option.
-
-The program will calculate the RMSE and check the required number of points.
-
-<h4>Linear affine transformation (1st order transformation)</h4>
-
-<dl>
- <dd> x' = ax + by +c
- <dd> y' = Ax + Bt +C
-</dl>
-
-The a,b,c,A,B,C 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
-images, and (2) the terrain or camera distortion effect can
-be ignored.
-
-<h4>Polynomial Transformation Matrix (2nd, 3d order transformation)</h4>
-
-<em>i.rectify</em> uses a first, second, or third order transformation
-matrix to calculate the registration coefficients. The number
-of control points required for a selected order of transformation
-(represented by n) is
-
-<dl>
-<dd>((n + 1) * (n + 2) / 2)
-</dl>
-
-or 3, 6, and 10 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 RMS error values for all the
-included control points are immediately recalculated when
-the user selects a different transformation order from the
-menu bar. The polynomial equations are performed using a
-modified Gaussian elimination method.
-
-<h3>Resampling method</h3>
-<p>The rectified data is resampled with one of seven different methods:
-<em>nearest</em>, <em>bilinear</em>, <em>cubic</em>, <em>lanczos</em>,
-<em>bilinear_f</em>, <em>cubic_f</em>, or <em>lanczos_f</em>.
-<p>The <em>method=nearest</em> method, which performs a nearest neighbor assignment,
-is the fastest of the resampling methods. It is primarily used for
-categorical data such as a land use classification, since it will not change
-the values of the data cells. The <em>method=bilinear</em> method determines the new
-value of the cell based on a weighted distance average of the 4 surrounding
-cells in the input map. The <em>method=cubic</em> method determines the new value of
-the cell based on a weighted distance average of the 16 surrounding cells in
-the input map. The <em>method=lanczos</em> method determines the new value of
-the cell based on a weighted distance average of the 25 surrounding cells in
-the input map.
-<p>The bilinear, cubic and lanczos interpolation methods are most appropriate for
-continuous data and cause some smoothing. These options should not be used
-with categorical data, since the cell values will be altered.
-<p>In the bilinear, cubic and lanczos methods, if any of the surrounding cells used to
-interpolate the new cell value are NULL, the resulting cell will be NULL, even if
-the nearest cell is not NULL. This will cause some thinning along NULL borders,
-such as the coasts of land areas in a DEM. The bilinear_f, cubic_f and lanczos_f
-interpolation methods can be used if thinning along NULL edges is not desired.
-These methods "fall back" to simpler interpolation methods along NULL borders.
-That is, from lanczos to cubic to bilinear to nearest.
-<p>If nearest neighbor assignment is used, the output map has the same raster
-format as the input map. If any of the other interpolations is used, the
-output map is written as floating point.
-
-<p><!--
-Note: In interactive mode it is possible to define a new file name
-for the target images. This is (currently) not provided in command line
-mode.
--->
-
-<h2>NOTES</h2>
-
-If <em>i.rectify</em> starts normally but after some time the following text is seen:
-<br><tt>
-ERROR: Error writing segment file
-</tt><br>
-the user may try the <b>-c</b> flag or the module needs more free space
-on the hard drive.
-
-
-<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="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>
-
-William R. Enslin,
-Michigan State University,
-Center for Remote Sensing
-
-<p>Modified for GRASS 5.0 by:<br>
-Luca Palmeri (palmeri at ux1.unipd.it)<br>
-Bill Hughes<br>
-Pierre de Mouveaux (pmx at audiovu.com)
-<br>
-CMD mode by Bob Covill
-
-<p><i>Last changed: $Date$</i>
Deleted: grass/trunk/vector/v.rectify/lanczos.c
===================================================================
--- grass/trunk/imagery/i.rectify/lanczos.c 2011-12-16 13:42:01 UTC (rev 49787)
+++ grass/trunk/vector/v.rectify/lanczos.c 2011-12-19 09:53:35 UTC (rev 49825)
@@ -1,109 +0,0 @@
-/*
- * Name
- * lanczos.c -- use lanczos interpolation for given row, col
- *
- * Description
- * lanczos returns the value in the buffer that is the result of
- * lanczos interpolation for the given row, column indices.
- * If the given row or column is outside the bounds of the input map,
- * the corresponding point in the output map is set to NULL.
- *
- * If single surrounding points in the interpolation matrix are
- * NULL they where filled with their neighbor
- */
-
-#include <grass/gis.h>
-#include <grass/raster.h>
-#include <math.h>
-#include "global.h"
-
-void p_lanczos(struct cache *ibuffer, /* input buffer */
- void *obufptr, /* ptr in output buffer */
- int cell_type, /* raster map type of obufptr */
- double *row_idx, /* row index (decimal) */
- double *col_idx, /* column index (decimal) */
- struct Cell_head *cellhd /* information of output map */
- )
-{
- int row; /* row indices for interp */
- int col; /* column indices for interp */
- int i, j, k;
- DCELL t, u; /* intermediate slope */
- DCELL result; /* result of interpolation */
- DCELL cell[25];
-
- /* cut indices to integer */
- row = (int)floor(*row_idx);
- col = (int)floor(*col_idx);
-
- /* check for out of bounds of map - if out of bounds set NULL value */
- if (row - 2 < 0 || row + 2 >= cellhd->rows ||
- col - 2 < 0 || col + 2 >= cellhd->cols) {
- Rast_set_null_value(obufptr, 1, cell_type);
- return;
- }
-
- k = 0;
- for (i = 0; i < 5; i++) {
- for (j = 0; j < 5; j++) {
- const DCELL *cellp = CPTR(ibuffer, row - 2 + i, col - 2 + j);
- if (Rast_is_d_null_value(cellp)) {
- Rast_set_null_value(obufptr, 1, cell_type);
- return;
- }
- cell[k++] = *cellp;
- }
- }
-
- /* do the interpolation */
- t = *col_idx - 0.5 - col;
- u = *row_idx - 0.5 - row;
-
- result = Rast_interp_lanczos(t, u, cell);
-
- Rast_set_d_value(obufptr, result, cell_type);
-}
-
-void p_lanczos_f(struct cache *ibuffer, /* input buffer */
- void *obufptr, /* ptr in output buffer */
- int cell_type, /* raster map type of obufptr */
- double *row_idx, /* row index (decimal) */
- double *col_idx, /* column index (decimal) */
- struct Cell_head *cellhd /* information of output map */
- )
-{
- int row; /* row indices for interp */
- int col; /* column indices for interp */
- DCELL *cellp, cell;
-
- /* cut indices to integer */
- row = (int)floor(*row_idx);
- col = (int)floor(*col_idx);
-
- /* check for out of bounds - if out of bounds set NULL value */
- if (row < 0 || row >= cellhd->rows || col < 0 || col >= cellhd->cols) {
- Rast_set_null_value(obufptr, 1, cell_type);
- return;
- }
-
- cellp = CPTR(ibuffer, row, col);
- /* if nearest is null, all the other interps will be null */
- if (Rast_is_d_null_value(cellp)) {
- Rast_set_null_value(obufptr, 1, cell_type);
- return;
- }
- cell = *cellp;
-
- p_lanczos(ibuffer, obufptr, cell_type, row_idx, col_idx, cellhd);
- /* fallback to bicubic if lanczos is null */
- if (Rast_is_d_null_value(obufptr)) {
- p_cubic(ibuffer, obufptr, cell_type, row_idx, col_idx, cellhd);
- /* fallback to bilinear if cubic is null */
- if (Rast_is_d_null_value(obufptr)) {
- p_bilinear(ibuffer, obufptr, cell_type, row_idx, col_idx, cellhd);
- /* fallback to nearest if bilinear is null */
- if (Rast_is_d_null_value(obufptr))
- Rast_set_d_value(obufptr, cell, cell_type);
- }
- }
-}
Modified: grass/trunk/vector/v.rectify/main.c
===================================================================
--- grass/trunk/imagery/i.rectify/main.c 2011-12-16 13:42:01 UTC (rev 49787)
+++ grass/trunk/vector/v.rectify/main.c 2011-12-19 09:53:35 UTC (rev 49825)
@@ -1,20 +1,14 @@
/****************************************************************************
*
- * MODULE: i.rectify
- * AUTHOR(S): William R. Enslin, Michigan State U. (original contributor)
- * Luca Palmeri <palmeri ux1.unipd.it>
- * Bill Hughes,
- * Pierre de Mouveaux <pmx audiovu.com>,
- * Bob Covill (original CMD version),
- * Markus Neteler <neteler itc.it>,
- * Bernhard Reiter <bernhard intevation.de>,
- * Glynn Clements <glynn gclements.plus.com>,
- * Hamish Bowman <hamish_b yahoo.com>,
- * Markus Metz
- * PURPOSE: calculate a transformation matrix and then convert x,y cell
- * coordinates to standard map coordinates for each pixel in the
- * image (control points can come from i.points or i.vpoints)
+ * 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
@@ -26,232 +20,148 @@
#include <stdlib.h>
#include <string.h>
-#include <grass/raster.h>
+#include <grass/vector.h>
#include <grass/glocale.h>
#include "global.h"
#include "crs.h"
-char *seg_mb;
-RASTER_MAP_TYPE map_type;
-int *ref_list;
-struct Ref ref;
-
-func interpolate;
-
/* georef coefficients */
-double E12[10], N12[10];
-double E21[10], N21[10];
+double E12[20], N12[20], Z12[20];
+double E21[20], N21[20], Z21[20];
-/* DELETED WITH CRS MODIFICATIONS
- double E12a, E12b, E12c, N12a, N12b, N12c;
- double E21a, E21b, E21c, N21a, N21b, N21c;
- */
-struct Cell_head target_window;
-void err_exit(char *, char *);
-
-/* modify this table to add new methods */
-struct menu menu[] = {
- {p_nearest, "nearest", "nearest neighbor"},
- {p_bilinear, "bilinear", "bilinear"},
- {p_cubic, "cubic", "cubic convolution"},
- {p_lanczos, "lanczos", "lanczos filter"},
- {p_bilinear_f, "bilinear_f", "bilinear with fallback"},
- {p_cubic_f, "cubic_f", "cubic convolution with fallback"},
- {p_lanczos_f, "lanczos_f", "lanczos filter with fallback"},
- {NULL, NULL, NULL}
-};
-
-static char *make_ipol_list(void);
-
int main(int argc, char *argv[])
{
- char group[INAME_LEN], extension[INAME_LEN];
+ char group[INAME_LEN];
int order; /* ADDED WITH CRS MODIFICATIONS */
- char *ipolname; /* name of interpolation method */
- int method;
- int n, i, m, k = 0;
- int got_file = 0, target_overwrite = 0;
- char *overstr;
- struct Cell_head cellhd;
+ 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 */
- *ifile, /* input files */
- *ext, /* extension */
- *tres, /* target resolution */
- *mem, /* amount of memory for cache */
- *interpol; /* interpolation method:
- nearest neighbor, bilinear, cubic */
- struct Flag *c, *a;
+ *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(_("imagery"));
+ G_add_keyword(_("vector"));
G_add_keyword(_("rectify"));
module->description =
- _("Rectifies an image by computing a coordinate "
- "transformation for each pixel in the image based on the "
+ _("Rectifies a vector by computing a coordinate "
+ "transformation for each object in the vector based on the "
"control points.");
- grp = G_define_standard_option(G_OPT_I_GROUP);
+ in_opt = G_define_standard_option(G_OPT_V_INPUT);
+ in_opt->required = YES;
- ifile = G_define_standard_option(G_OPT_R_INPUTS);
- ifile->required = NO;
+ out_opt = G_define_standard_option(G_OPT_V_OUTPUT);
+ out_opt->required = YES;
- ext = G_define_option();
- ext->key = "extension";
- ext->type = TYPE_STRING;
- ext->required = YES;
- ext->multiple = NO;
- ext->description = _("Output raster map(s) suffix");
+ 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");
- tres = G_define_option();
- tres->key = "res";
- tres->type = TYPE_DOUBLE;
- tres->required = NO;
- tres->description = _("Target resolution (ignored if -c flag used)");
+ flag_use3d = G_define_flag();
+ flag_use3d->key = '3';
+ flag_use3d->description = _("Perform 3D transformation");
- mem = G_define_option();
- mem->key = "memory";
- mem->type = TYPE_DOUBLE;
- mem->key_desc = "memory in MB";
- mem->required = NO;
- mem->answer = "300";
- mem->description = _("Amount of memory to use in MB");
+ 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");
- ipolname = make_ipol_list();
+ no_topo = G_define_flag();
+ no_topo->key = 'b';
+ no_topo->description = _("Do not build topology for output");
- interpol = G_define_option();
- interpol->key = "method";
- interpol->type = TYPE_STRING;
- interpol->required = NO;
- interpol->answer = "nearest";
- interpol->options = ipolname;
- interpol->description = _("Interpolation method to use");
-
- c = G_define_flag();
- c->key = 'c';
- c->description =
- _("Use current region settings in target location (def.=calculate smallest area)");
-
- a = G_define_flag();
- a->key = 'a';
- a->description = _("Rectify all raster maps in group");
-
if (G_parser(argc, argv))
exit(EXIT_FAILURE);
- /* get the method */
- for (method = 0; (ipolname = menu[method].name); method++)
- if (strcmp(ipolname, interpol->answer) == 0)
- break;
- if (!ipolname)
- G_fatal_error(_("<%s=%s> unknown %s"),
- interpol->key, interpol->answer, interpol->key);
- interpolate = menu[method].method;
-
- G_strip(grp->answer);
- strcpy(group, grp->answer);
- strcpy(extension, ext->answer);
- order = atoi(val->answer);
-
- seg_mb = NULL;
- if (mem->answer) {
- if (atoi(mem->answer) > 0)
- seg_mb = mem->answer;
+ if (grp->answer) {
+ G_strip(grp->answer);
+ strcpy(group, grp->answer);
}
+ else
+ group[0] = '\0';
- if (!ifile->answers)
- a->answer = 1; /* force all */
+ points_file = pfile->answer;
- /* Find out how many files on command line */
- if (!a->answer) {
- for (m = 0; ifile->answers[m]; m++) {
- k = m;
- }
- k++;
- }
+ 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);
- /* determine the number of files in this group */
- if (I_get_group_ref(group, &ref) <= 0) {
- G_warning(_("Location: %s"), G_location());
- G_warning(_("Mapset: %s"), G_mapset());
- G_fatal_error(_("Group <%s> does not exist"), grp->answer);
- }
+ 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 (ref.nfiles <= 0) {
- G_important_message(_("Group <%s> contains no raster maps; run i.group"),
- grp->answer);
- exit(EXIT_SUCCESS);
- }
-
- ref_list = (int *)G_malloc(ref.nfiles * sizeof(int));
-
- if (a->answer) {
- for (n = 0; n < ref.nfiles; n++) {
- ref_list[n] = 1;
+ 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 {
- char xname[GNAME_MAX], xmapset[GMAPSET_MAX], *name, *mapset;
+ else
+ rms_sep = NULL;
- for (n = 0; n < ref.nfiles; n++)
- ref_list[n] = 0;
-
- for (m = 0; m < k; m++) {
- got_file = 0;
- if (G_name_is_fully_qualified(ifile->answers[m], xname, xmapset)) {
- name = xname;
- mapset = xmapset;
- }
- else {
- name = ifile->answers[m];
- mapset = NULL;
- }
-
- got_file = 0;
- for (n = 0; n < ref.nfiles; n++) {
- if (mapset) {
- if (strcmp(name, ref.file[n].name) == 0 &&
- strcmp(mapset, ref.file[n].mapset) == 0) {
- got_file = 1;
- ref_list[n] = 1;
- break;
- }
- }
- else {
- if (strcmp(name, ref.file[n].name) == 0) {
- got_file = 1;
- ref_list[n] = 1;
- break;
- }
- }
- }
- if (got_file == 0)
- err_exit(ifile->answers[m], group);
- }
+ /* 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);
}
- /* read the control points for the group */
- get_control_points(group, order);
-
/* get the target */
get_target(group);
@@ -261,92 +171,85 @@
if (!target_overwrite) {
/* check if output exists in target location/mapset */
- char result[GNAME_MAX];
select_target_env();
- for (i = 0; i < ref.nfiles; i++) {
- if (!ref_list[i])
- continue;
-
- strcpy(result, ref.file[i].name);
- strcat(result, extension);
-
- if (G_legal_filename(result) < 0)
- G_fatal_error(_("Extension <%s> is illegal"), extension);
- if (G_find_raster2(result, G_mapset())) {
- G_warning(_("The following raster map already exists in"));
- G_warning(_("target LOCATION %s, MAPSET %s:"),
- G_location(), G_mapset());
- G_warning("<%s>", result);
- G_fatal_error(_("Orthorectification cancelled."));
- }
+ 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");
- /* do not use current region in target location */
- if (!c->answer) {
- double res = -1;
+ 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++;
+ }
- if (tres->answer) {
- if (!((res = atof(tres->answer)) > 0))
- G_warning(_("Target resolution must be > 0, ignored"));
- }
- /* Calculate smallest region */
- if (a->answer)
- Rast_get_cellhd(ref.file[0].name, ref.file[0].mapset, &cellhd);
- else
- Rast_get_cellhd(ifile->answers[0], ref.file[0].mapset, &cellhd);
+ Vect_rewind(&In);
- georef_window(&cellhd, &target_window, order, res);
+ 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();
- G_verbose_message(_("Using region: N=%f S=%f, E=%f W=%f"), target_window.north,
- target_window.south, target_window.east, target_window.west);
+ Vect_close(&In);
- exec_rectify(order, extension, interpol->answer);
-
G_done_msg(" ");
exit(EXIT_SUCCESS);
}
-
-
-void err_exit(char *file, char *grp)
-{
- int n;
-
- G_warning(_("Input raster map <%s> does not exist in group <%s>."),
- file, grp);
- G_message(_("Try:"));
-
- for (n = 0; n < ref.nfiles; n++)
- G_message("%s", ref.file[n].name);
-
- G_fatal_error(_("Exit!"));
-}
-
-static char *make_ipol_list(void)
-{
- int size = 0;
- int i;
- char *buf;
-
- for (i = 0; menu[i].name; i++)
- size += strlen(menu[i].name) + 1;
-
- buf = G_malloc(size);
- *buf = '\0';
-
- for (i = 0; menu[i].name; i++) {
- if (i)
- strcat(buf, ",");
- strcat(buf, menu[i].name);
- }
-
- return buf;
-}
Deleted: grass/trunk/vector/v.rectify/nearest.c
===================================================================
--- grass/trunk/imagery/i.rectify/nearest.c 2011-12-16 13:42:01 UTC (rev 49787)
+++ grass/trunk/vector/v.rectify/nearest.c 2011-12-19 09:53:35 UTC (rev 49825)
@@ -1,42 +0,0 @@
-
-/*
- * nearest.c - returns the nearest neighbor to a given
- * x,y position
- */
-
-#include <math.h>
-#include <grass/gis.h>
-#include <grass/raster.h>
-#include "global.h"
-
-void p_nearest(struct cache *ibuffer, /* input buffer */
- void *obufptr, /* ptr in output buffer */
- int cell_type, /* raster map type of obufptr */
- double *row_idx, /* row index in input matrix */
- double *col_idx, /* column index in input matrix */
- struct Cell_head *cellhd /* cell header of input layer */
- )
-{
- int row, col; /* row/col of nearest neighbor */
- DCELL *cellp;
-
- /* cut indices to integer and get nearest cell */
- /* the row_idx, col_idx correction for bilinear/bicubic does not apply here */
- row = (int)floor(*row_idx);
- col = (int)floor(*col_idx);
-
- /* check for out of bounds - if out of bounds set NULL value */
- if (row < 0 || row >= cellhd->rows || col < 0 || col >= cellhd->cols) {
- Rast_set_null_value(obufptr, 1, cell_type);
- return;
- }
-
- cellp = CPTR(ibuffer, row, col);
-
- if (Rast_is_d_null_value(cellp)) {
- Rast_set_null_value(obufptr, 1, cell_type);
- return;
- }
-
- Rast_set_d_value(obufptr, *cellp, cell_type);
-}
Deleted: grass/trunk/vector/v.rectify/readcell.c
===================================================================
--- grass/trunk/imagery/i.rectify/readcell.c 2011-12-16 13:42:01 UTC (rev 49787)
+++ grass/trunk/vector/v.rectify/readcell.c 2011-12-19 09:53:35 UTC (rev 49825)
@@ -1,130 +0,0 @@
-/*
- * readcell.c - reads an entire cell layer into a buffer
- *
- */
-
-#include <stdio.h>
-#include <stdlib.h>
-#include <string.h>
-#include <sys/types.h>
-#include <sys/stat.h>
-#include <fcntl.h>
-#include <unistd.h>
-#include <grass/gis.h>
-#include <grass/raster.h>
-#include <grass/glocale.h>
-#include "global.h"
-
-struct cache *readcell(int fdi, const char *size)
-{
- DCELL *tmpbuf;
- struct cache *c;
- int nrows;
- int ncols;
- int row;
- char *filename;
- int nx, ny;
- int nblocks;
- int i;
-
- nrows = Rast_input_window_rows();
- ncols = Rast_input_window_cols();
-
- ny = (nrows + BDIM - 1) / BDIM;
- nx = (ncols + BDIM - 1) / BDIM;
-
- if (size)
- nblocks = atoi(size) * ((1 << 20) / sizeof(block));
- else
- nblocks = (nx + ny) * 2; /* guess */
-
- if (nblocks > nx * ny)
- nblocks = nx * ny;
-
- c = G_malloc(sizeof(struct cache));
- c->stride = nx;
- c->nblocks = nblocks;
- c->grid = (block **) G_calloc(nx * ny, sizeof(block *));
- c->blocks = (block *) G_malloc(nblocks * sizeof(block));
- c->refs = (int *)G_malloc(nblocks * sizeof(int));
-
- if (nblocks < nx * ny) {
- /* Temporary file must be created in input location */
- filename = G_tempfile();
- c->fd = open(filename, O_RDWR | O_CREAT | O_EXCL, 0600);
- if (c->fd < 0)
- G_fatal_error(_("Unable to open temporary file"));
- remove(filename);
- }
- else
- c->fd = -1;
-
- G_important_message(_("Allocating memory and reading input map..."));
- G_percent(0, nrows, 5);
-
- for (i = 0; i < c->nblocks; i++)
- c->refs[i] = -1;
-
- tmpbuf = (DCELL *) G_malloc(nx * sizeof(block));
-
- for (row = 0; row < nrows; row += BDIM) {
- int x, y;
-
- for (y = 0; y < BDIM; y++) {
- G_percent(row + y, nrows, 5);
-
- if (row + y >= nrows)
- break;
-
- Rast_get_d_row(fdi, &tmpbuf[y * nx * BDIM], row + y);
- }
-
- for (x = 0; x < nx; x++)
- for (y = 0; y < BDIM; y++)
- if (c->fd >= 0) {
- if (write
- (c->fd, &tmpbuf[(y * nx + x) * BDIM],
- BDIM * sizeof(DCELL)) < 0)
- G_fatal_error(_("Error writing segment file"));
- }
- else
- memcpy(&c->blocks[BKIDX(c, HI(row), x)][LO(y)][0],
- &tmpbuf[(y * nx + x) * BDIM],
- BDIM * sizeof(DCELL));
- }
-
- G_free(tmpbuf);
-
- if (c->fd < 0)
- for (i = 0; i < c->nblocks; i++) {
- c->grid[i] = &c->blocks[i];
- c->refs[i] = i;
- }
-
- return c;
-}
-
-block *get_block(struct cache * c, int idx)
-{
- int replace = rand() % c->nblocks;
- block *p = &c->blocks[replace];
- int ref = c->refs[replace];
- off_t offset = (off_t) idx * sizeof(DCELL) << L2BSIZE;
-
- if (c->fd < 0)
- G_fatal_error(_("Internal error: cache miss on fully-cached map"));
-
- if (ref >= 0)
- c->grid[ref] = NULL;
-
- c->grid[idx] = p;
- c->refs[replace] = idx;
-
- if (lseek(c->fd, offset, SEEK_SET) < 0)
- G_fatal_error(_("Error seeking on segment file"));
-
- if (read(c->fd, p, sizeof(block)) < 0)
- G_fatal_error(_("Error writing segment file"));
-
- return p;
-}
Deleted: grass/trunk/vector/v.rectify/rectify.c
===================================================================
--- grass/trunk/imagery/i.rectify/rectify.c 2011-12-16 13:42:01 UTC (rev 49787)
+++ grass/trunk/vector/v.rectify/rectify.c 2011-12-19 09:53:35 UTC (rev 49825)
@@ -1,119 +0,0 @@
-#include <unistd.h>
-#include <string.h>
-
-#include <grass/raster.h>
-#include <grass/glocale.h>
-
-#include "global.h"
-#include "crs.h" /* CRS HEADER FILE */
-
-/* Modified to support Grass 5.0 fp format 11 april 2000
- *
- * Pierre de Mouveaux - pmx at audiovu.com
- *
- */
-
-int rectify(char *name, char *mapset, char *result, int order, char *interp_method)
-{
- struct Cell_head cellhd;
- int ncols, nrows;
- int row, col;
- double row_idx, col_idx;
- int infd, cell_size, outfd;
- void *trast, *tptr;
- double n1, e1, nx, ex;
- struct cache *ibuffer;
-
- select_current_env();
- Rast_get_cellhd(name, mapset, &cellhd);
-
- /* open the file to be rectified
- * set window to cellhd first to be able to read file exactly
- */
- Rast_set_input_window(&cellhd);
- infd = Rast_open_old(name, mapset);
- map_type = Rast_get_map_type(infd);
- cell_size = Rast_cell_size(map_type);
-
- ibuffer = readcell(infd, seg_mb);
-
- Rast_close(infd); /* (pmx) 17 april 2000 */
-
- G_message(_("Rectify <%s@%s> (location <%s>)"),
- name, mapset, G_location());
- select_target_env();
- G_message(_("into <%s@%s> (location <%s>) ..."),
- result, G_mapset(), G_location());
-
- nrows = target_window.rows;
- ncols = target_window.cols;
-
- if (strcmp(interp_method, "nearest") != 0) {
- map_type = DCELL_TYPE;
- cell_size = Rast_cell_size(map_type);
- }
-
- /* open the result file into target window
- * this open must be first since we change the window later
- * raster maps open for writing are not affected by window changes
- * but those open for reading are
- */
-
- outfd = Rast_open_new(result, map_type);
- trast = Rast_allocate_output_buf(map_type);
-
- for (row = 0; row < nrows; row++) {
- n1 = target_window.north - (row + 0.5) * target_window.ns_res;
-
- G_percent(row, nrows, 2);
-
- Rast_set_null_value(trast, ncols, map_type);
- tptr = trast;
- for (col = 0; col < ncols; col++) {
- e1 = target_window.west + (col + 0.5) * target_window.ew_res;
-
- /* backwards transformation of target cell center */
- CRS_georef(e1, n1, &ex, &nx, E21, N21, order);
-
- /* convert to row/column indices of source raster */
- row_idx = (cellhd.north - nx) / cellhd.ns_res;
- col_idx = (ex - cellhd.west) / cellhd.ew_res;
-
- /* resample data point */
- interpolate(ibuffer, tptr, map_type, &row_idx, &col_idx, &cellhd);
-
- tptr = G_incr_void_ptr(tptr, cell_size);
- }
- Rast_put_row(outfd, trast, map_type);
- }
- G_percent(1, 1, 1);
-
- Rast_close(outfd); /* (pmx) 17 april 2000 */
- G_free(trast);
-
- close(ibuffer->fd);
- G_free(ibuffer);
-
- Rast_get_cellhd(result, G_mapset(), &cellhd);
-
- if (cellhd.proj == 0) { /* x,y imagery */
- cellhd.proj = target_window.proj;
- cellhd.zone = target_window.zone;
- }
-
- if (target_window.proj != cellhd.proj) {
- cellhd.proj = target_window.proj;
- G_warning(_("Raster map <%s@%s>: projection don't match current settings"),
- name, mapset);
- }
-
- if (target_window.zone != cellhd.zone) {
- cellhd.zone = target_window.zone;
- G_warning(_("Raster map <%s@%s>: zone don't match current settings"),
- name, mapset);
- }
-
- select_current_env();
-
- return 1;
-}
Deleted: grass/trunk/vector/v.rectify/report.c
===================================================================
--- grass/trunk/imagery/i.rectify/report.c 2011-12-16 13:42:01 UTC (rev 49787)
+++ grass/trunk/vector/v.rectify/report.c 2011-12-19 09:53:35 UTC (rev 49825)
@@ -1,33 +0,0 @@
-#include <grass/glocale.h>
-#include "global.h"
-
-int report(long rectify, int ok)
-{
- int minutes, hours;
- long seconds;
- long ncells;
-
- G_message("%s", ok ? _("complete") : _("failed"));
-
- if (!ok)
- return 1;
-
- seconds = rectify;
- minutes = seconds / 60;
- hours = minutes / 60;
- minutes -= hours * 60;
- ncells = target_window.rows * target_window.cols;
- G_verbose_message(_("%d rows, %d cols (%ld cells) completed in"),
- target_window.rows, target_window.cols, ncells);
- if (hours)
- G_verbose_message(_("%d:%02d:%02ld hours"), hours, minutes, seconds % 60);
- else
- G_verbose_message(_("%d:%02ld minutes"), minutes, seconds % 60);
- if (seconds)
- G_verbose_message(_("%.1f cells per minute"),
- (60.0 * ncells) / ((double)seconds));
-
- G_message("-----------------------------------------------");
-
- return 1;
-}
Modified: grass/trunk/vector/v.rectify/target.c
===================================================================
--- grass/trunk/imagery/i.rectify/target.c 2011-12-16 13:42:01 UTC (rev 49787)
+++ grass/trunk/vector/v.rectify/target.c 2011-12-19 09:53:35 UTC (rev 49825)
@@ -10,10 +10,16 @@
char buf[1024];
int stat;
- if (!I_get_target(group, location, mapset)) {
- sprintf(buf, _("Target information for group <%s> missing"), group);
- goto error;
+ 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) {
@@ -25,7 +31,6 @@
stat = G__mapset_permissions(mapset);
if (stat > 0) {
G__setenv("MAPSET", mapset);
- G_get_window(&target_window);
select_current_env();
return 1;
}
Copied: grass/trunk/vector/v.rectify/v.rectify.html (from rev 49787, grass/trunk/imagery/i.rectify/i.rectify.html)
===================================================================
--- grass/trunk/vector/v.rectify/v.rectify.html (rev 0)
+++ grass/trunk/vector/v.rectify/v.rectify.html 2011-12-19 09:53:35 UTC (rev 49825)
@@ -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>
More information about the grass-commit
mailing list