[GRASS-SVN] r51546 - grass/trunk/vector/v.rectify
svn_grass at osgeo.org
svn_grass at osgeo.org
Thu Apr 26 02:54:29 EDT 2012
Author: mmetz
Date: 2012-04-25 23:54:29 -0700 (Wed, 25 Apr 2012)
New Revision: 51546
Added:
grass/trunk/vector/v.rectify/orthorot.c
grass/trunk/vector/v.rectify/svdm.c
Modified:
grass/trunk/vector/v.rectify/Makefile
grass/trunk/vector/v.rectify/cp.c
grass/trunk/vector/v.rectify/crs.h
grass/trunk/vector/v.rectify/crs3d.c
grass/trunk/vector/v.rectify/env.c
grass/trunk/vector/v.rectify/global.h
grass/trunk/vector/v.rectify/main.c
grass/trunk/vector/v.rectify/target.c
grass/trunk/vector/v.rectify/v.rectify.html
Log:
v.rectify: othogonal rotation
Modified: grass/trunk/vector/v.rectify/Makefile
===================================================================
--- grass/trunk/vector/v.rectify/Makefile 2012-04-25 12:50:56 UTC (rev 51545)
+++ grass/trunk/vector/v.rectify/Makefile 2012-04-26 06:54:29 UTC (rev 51546)
@@ -2,8 +2,9 @@
PGM = v.rectify
-LIBES = $(IMAGERYLIB) $(VECTORLIB) $(GISLIB)
-DEPENDENCIES = $(IMAGERYDEP) $(VECTORDEP) $(GISDEP)
+LIBES = $(GMATHLIB) $(IMAGERYLIB) $(VECTORLIB) $(GISLIB)
+DEPENDENCIES = $(GMATHDEP) $(IMAGERYDEP) $(VECTORDEP) $(GISDEP)
+
EXTRA_INC = $(VECT_INC)
EXTRA_CFLAGS = $(VECT_CFLAGS)
Modified: grass/trunk/vector/v.rectify/cp.c
===================================================================
--- grass/trunk/vector/v.rectify/cp.c 2012-04-25 12:50:56 UTC (rev 51545)
+++ grass/trunk/vector/v.rectify/cp.c 2012-04-26 06:54:29 UTC (rev 51546)
@@ -1,4 +1,5 @@
#include <stdlib.h>
+#include <stdio.h>
#include <string.h>
#include <math.h>
#include <unistd.h>
@@ -36,7 +37,7 @@
static void compute_rms(struct Control_Points *cp, struct Control_Points_3D *cp3,
- int order, int use3d, char *sep)
+ int order, int use3d, int orthorot, char *sep, FILE *fp)
{
int n;
int count, npoints;
@@ -52,13 +53,13 @@
/* 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",
+ fprintf(fp, "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",
+ fprintf(fp, "index%sfwd_dx%sfwd_dy%sback_dx%sback_dy%sfwd_RMS%sback_RMS",
sep, sep, sep, sep, sep, sep);
- printf("\n");
+ fprintf(fp, "\n");
if (use3d)
npoints = cp3->count;
@@ -70,7 +71,7 @@
double fx, fy, fz, fd, fd2;
double rx, ry, rz, rd, rd2;
- if (use3d) {
+ if (use3d || orthorot) {
if (cp3->status[n] <= 0)
continue;
}
@@ -83,10 +84,14 @@
/* forward: source -> target */
if (use3d) {
- CRS_georef_3d(cp3->e1[n], cp3->n1[n], cp3->z1[n],
- &e2, &n2, &z2,
- E12, N12, Z12,
- order);
+ if (orthorot)
+ CRS_georef_or(cp3->e1[n], cp3->n1[n], cp3->z1[n],
+ &e2, &n2, &z2, OR12);
+ else
+ 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]);
@@ -109,10 +114,14 @@
/* backward: target -> source */
if (use3d) {
- CRS_georef_3d(cp3->e2[n], cp3->n2[n], cp3->z2[n],
- &e1, &n1, &z1,
- E21, N21, Z21,
- order);
+ if (orthorot)
+ CRS_georef_or(cp3->e2[n], cp3->n2[n], cp3->z2[n],
+ &e1, &n1, &z1, OR21);
+ else
+ 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]);
@@ -135,17 +144,17 @@
/* print index, forward difference, backward difference,
* forward rms, backward rms */
- printf("%d", n + 1);
- printf("%s%f%s%f", sep, fx, sep, fy);
+ fprintf(fp, "%d", n + 1);
+ fprintf(fp, "%s%f%s%f", sep, fx, sep, fy);
if (use3d)
- printf("%s%f", sep, fz);
- printf("%s%f%s%f", sep, rx, sep, ry);
+ fprintf(fp, "%s%f", sep, fz);
+ fprintf(fp, "%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);
+ fprintf(fp, "%s%f", sep, rz);
+ fprintf(fp, "%s%.4f", sep, fd);
+ fprintf(fp, "%s%.4f", sep, rd);
- printf("\n");
+ fprintf(fp, "\n");
}
if (count > 0) {
@@ -162,17 +171,17 @@
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);
+ fprintf(fp, "%d", count);
+ fprintf(fp, "%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);
+ fprintf(fp, "%s%f", sep, fwd.z);
+ fprintf(fp, "%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);
+ fprintf(fp, "%s%f", sep, rev.z);
+ fprintf(fp, "%s%.4f", sep, fwd.rms);
+ fprintf(fp, "%s%.4f", sep, rev.rms);
- printf("\n");
+ fprintf(fp, "\n");
}
@@ -186,6 +195,7 @@
if (status < 0)
return 1;
+
i = (cp->count)++;
size = cp->count * sizeof(double);
cp->e1 = (double *)G_realloc(cp->e1, size);
@@ -271,7 +281,8 @@
}
-int get_control_points(char *group, char *pfile, int order, int use3d, int rms, char *sep)
+int get_control_points(char *group, char *pfile, int order, int use3d,
+ int orthorot, int rms, char *sep, FILE *fpr)
{
char msg[200];
struct Control_Points cp;
@@ -279,6 +290,12 @@
int ret = 0;
int order_pnts[2][3] = {{ 3, 6, 10 }, { 4, 10, 20 }};
+ cp.count = cp3.count = 0;
+ cp.e1 = cp.e2 = cp3.e1 = cp3.e2 = NULL;
+ cp.n1 = cp.n2 = cp3.n1 = cp3.n2 = NULL;
+ cp3.z1 = cp3.z2 = NULL;
+ cp.status = cp3.status = NULL;
+
msg[0] = '\0';
if (use3d) {
@@ -299,7 +316,11 @@
return 0;
}
- ret = CRS_compute_georef_equations_3d(&cp3, E12, N12, Z12, E21, N21, Z21, order);
+ if (orthorot)
+ ret = CRS_compute_georef_equations_or(&cp3, OR12, OR21);
+ else
+ ret = CRS_compute_georef_equations_3d(&cp3, E12, N12, Z12,
+ E21, N21, Z21, order);
}
else if (pfile) {
/* read 2D GCPs from points file */
@@ -336,7 +357,7 @@
case 0:
sprintf(&msg[strlen(msg)],
_("Not enough active control points for current order, %d are required."),
- order_pnts[use3d != 0][order]);
+ (orthorot ? 3 : order_pnts[use3d != 0][order - 1]));
break;
case -1:
strcat(msg, _("Poorly placed control points."));
@@ -355,9 +376,8 @@
G_fatal_error(msg);
if (rms) {
- compute_rms(&cp, &cp3, order, use3d, sep);
+ compute_rms(&cp, &cp3, order, use3d, orthorot, sep, fpr);
}
-
-
+
return 1;
}
Modified: grass/trunk/vector/v.rectify/crs.h
===================================================================
--- grass/trunk/vector/v.rectify/crs.h 2012-04-25 12:50:56 UTC (rev 51545)
+++ grass/trunk/vector/v.rectify/crs.h 2012-04-26 06:54:29 UTC (rev 51546)
@@ -1,4 +1,6 @@
+#ifndef CRS3D_H
+#define CRS3D_H_
#define MAXORDER 3
@@ -24,3 +26,26 @@
double *, double *, double *,
double *, double *, double *,
int);
+
+/* helmertg.c */
+int CRS_compute_georef_equations_hg(struct Control_Points_3D *,
+ double *, double *);
+int CRS_georef_hg(double, double, double,
+ double *, double *, double *,
+ double *);
+
+/* helmertq.c */
+int CRS_compute_georef_equations_hq(struct Control_Points_3D *,
+ double *, double *);
+int CRS_georef_hq(double, double, double,
+ double *, double *, double *,
+ double *);
+
+/* orthorot.c */
+int CRS_compute_georef_equations_or(struct Control_Points_3D *,
+ double *, double *);
+int CRS_georef_or(double, double, double,
+ double *, double *, double *,
+ double *);
+
+#endif
Modified: grass/trunk/vector/v.rectify/crs3d.c
===================================================================
--- grass/trunk/vector/v.rectify/crs3d.c 2012-04-25 12:50:56 UTC (rev 51545)
+++ grass/trunk/vector/v.rectify/crs3d.c 2012-04-26 06:54:29 UTC (rev 51546)
@@ -157,10 +157,8 @@
double *tempptr;
int status;
- /*
if (order < 1 || order > MAXORDER)
return MPARMERR;
- */
/* CALCULATE THE FORWARD TRANSFORMATION COEFFICIENTS */
@@ -235,6 +233,8 @@
3rd order: 10 20
*/
+ m.n = numactive + 1;
+
if (order == 1)
m.n = 4;
else if (order == 2)
@@ -477,7 +477,7 @@
/* co-linear points result in an undefined matrix, and nearly */
/* co-linear points results in a solution with rounding error */
- if (pivot == 0.0)
+ if (fabs(pivot) < GRASS_EPSILON)
return MUNSOLVABLE;
/* if row with highest pivot is not the current row, switch them */
Modified: grass/trunk/vector/v.rectify/env.c
===================================================================
--- grass/trunk/vector/v.rectify/env.c 2012-04-25 12:50:56 UTC (rev 51545)
+++ grass/trunk/vector/v.rectify/env.c 2012-04-26 06:54:29 UTC (rev 51546)
@@ -1,5 +1,8 @@
+#include <stdio.h>
#include <unistd.h>
+#include <grass/gis.h>
#include "global.h"
+
static int which_env = -1; /* 0 = cur, 1 = target */
int select_current_env(void)
Modified: grass/trunk/vector/v.rectify/global.h
===================================================================
--- grass/trunk/vector/v.rectify/global.h 2012-04-25 12:50:56 UTC (rev 51545)
+++ grass/trunk/vector/v.rectify/global.h 2012-04-26 06:54:29 UTC (rev 51546)
@@ -1,14 +1,13 @@
-#include <grass/gis.h>
-#include <grass/imagery.h>
-
/* georef coefficients */
extern double E12[20], N12[20], Z12[20];
extern double E21[20], N21[20], Z21[20];
+extern double HG12[20], HG21[20], HQ12[20], HQ21[20];
+extern double OR12[20], OR21[20];
/* cp.c */
-int get_control_points(char *, char *, int, int, int, char *);
+int get_control_points(char *, char *, int, int, int, int, char *, FILE *);
/* env.c */
int select_current_env(void);
Modified: grass/trunk/vector/v.rectify/main.c
===================================================================
--- grass/trunk/vector/v.rectify/main.c 2012-04-25 12:50:56 UTC (rev 51545)
+++ grass/trunk/vector/v.rectify/main.c 2012-04-26 06:54:29 UTC (rev 51546)
@@ -18,8 +18,11 @@
*****************************************************************************/
#include <stdlib.h>
+#include <stdio.h>
#include <string.h>
+#include <grass/gis.h>
+#include <grass/imagery.h>
#include <grass/vector.h>
#include <grass/glocale.h>
@@ -31,12 +34,14 @@
double E12[20], N12[20], Z12[20];
double E21[20], N21[20], Z21[20];
+double HG12[20], HG21[20], HQ12[20], HQ21[20];
+double OR12[20], OR21[20];
int main(int argc, char *argv[])
{
char group[INAME_LEN];
- int order; /* ADDED WITH CRS MODIFICATIONS */
+ int order, orthorot;
int n, i, nlines, type;
int target_overwrite = 0;
char *points_file, *overstr, *rms_sep;
@@ -45,15 +50,17 @@
struct line_cats *Cats;
double x, y, z;
int use3d;
+ FILE *fp;
struct Option *grp, /* imagery group */
*val, /* transformation order */
*in_opt, /* input vector name */
*out_opt, /* output vector name */
*pfile, /* text file with GCPs */
+ *rfile, /* text file to hold RMS errors */
*sep; /* field separator for RMS report */
- struct Flag *flag_use3d, *no_topo, *print_rms;
+ struct Flag *flag_use3d, *no_topo, *print_rms, *ortho;
struct GModule *module;
@@ -81,6 +88,11 @@
pfile->description = _("Name of input file with control points");
pfile->required = NO;
+ rfile = G_define_standard_option(G_OPT_F_INPUT);
+ rfile->key = "rmsfile";
+ rfile->description = _("Name of output file with RMS errors (if omitted or '-' output to stdout");
+ rfile->required = NO;
+
val = G_define_option();
val->key = "order";
val->type = TYPE_INTEGER;
@@ -97,6 +109,10 @@
flag_use3d->key = '3';
flag_use3d->description = _("Perform 3D transformation");
+ ortho = G_define_flag();
+ ortho->key = 'o';
+ ortho->description = _("Perform orthogonal 3D transformation");
+
print_rms = G_define_flag();
print_rms->key = 'r';
print_rms->label = _("Print RMS errors");
@@ -131,10 +147,16 @@
Vect_set_open_level(1);
Vect_open_old2(&In, in_opt->answer, "", "");
- use3d = (Vect_is_3d(&In) && flag_use3d->answer);
+ use3d = (Vect_is_3d(&In) &&
+ (flag_use3d->answer || ortho->answer));
+ if (!use3d && (flag_use3d->answer || ortho->answer))
+ G_fatal_error(_("3D transformation requires a 3D vector"));
+
if (use3d && !points_file)
- G_fatal_error(_("A file with 3D control points is needed for 3d transformation"));
+ G_fatal_error(_("A file with 3D control points is needed for 3D transformation"));
+
+ orthorot = ortho->answer;
if (print_rms->answer) {
if (sep->answer) {
@@ -154,9 +176,22 @@
else
rms_sep = NULL;
+ if (rfile->answer) {
+ if (strcmp(rfile->answer, "-")) {
+ fp = fopen(rfile->answer, "w");
+ if (!fp)
+ G_fatal_error(_("Unable to open file '%s' for writing"),
+ rfile->answer);
+ }
+ else
+ fp = stdout;
+ }
+ else
+ fp = stdout;
+
/* read the control points for the group */
- get_control_points(group, points_file, order, use3d,
- print_rms->answer, rms_sep);
+ get_control_points(group, points_file, order, use3d, orthorot,
+ print_rms->answer, rms_sep, fp);
if (print_rms->answer) {
Vect_close(&In);
@@ -222,8 +257,12 @@
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);
+ if (orthorot)
+ CRS_georef_or(Points->x[n], Points->y[n], Points->z[n],
+ &x, &y, &z, OR12);
+ else
+ CRS_georef_3d(Points->x[n], Points->y[n], Points->z[n],
+ &x, &y, &z, E12, N12, Z12, order);
}
else {
I_georef(Points->x[n], Points->y[n],
Added: grass/trunk/vector/v.rectify/orthorot.c
===================================================================
--- grass/trunk/vector/v.rectify/orthorot.c (rev 0)
+++ grass/trunk/vector/v.rectify/orthorot.c 2012-04-26 06:54:29 UTC (rev 51546)
@@ -0,0 +1,603 @@
+
+/***********************************************************************
+
+ orthorot.c
+
+ written by: Markus Metz
+
+ loosely based on crs3d.c
+
+ 2D/3D Transformation with orthogonal rotation
+
+************************************************************************/
+
+#include <stdio.h>
+#include <string.h>
+#include <stdlib.h>
+#include <math.h>
+#include <limits.h>
+
+#include <grass/gis.h>
+#include <grass/gmath.h>
+#include <grass/imagery.h>
+#include <grass/glocale.h>
+
+#include "crs.h"
+
+#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 *, int);
+static int calcscale(struct Control_Points_3D *, double *);
+
+void transpose_matrix(int, double **, double **);
+void matmult(int, int, double **, double **, double **);
+void copy_matrix(int, int, double **, double **);
+void init_matrix(int, int, double **);
+void scale_matrix(int, int, double, double **, double **);
+void matrix_multiply(int, int, double **, double *, double *);
+void subtract_matrix(int, int, double **, double **, double **);
+double trace(int, int, double **);
+
+void svd(int, int, double **, double **, double *, double **);
+
+/***********************************************************************
+
+ TRANSFORM A SINGLE COORDINATE PAIR.
+
+************************************************************************/
+
+int CRS_georef_or(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 OR[] /* TRANSFORMATION COEFFICIENTS */
+ )
+{
+ *e = OR[9] + OR[12] * (OR[0] * e1 + OR[1] * n1 + OR[2] * z1);
+ *n = OR[10] + OR[13] * (OR[3] * e1 + OR[4] * n1 + OR[5] * z1);
+ *z = OR[11] + OR[14] * (OR[6] * e1 + OR[7] * n1 + OR[8] * z1);
+
+ return MSUCCESS;
+}
+
+/***********************************************************************
+
+ COMPUTE THE FORWARD AND BACKWARD GEOREFFERENCING COEFFICIENTS
+ BASED ON A SET OF CONTROL POINTS
+
+ ORTHOGONAL TRANSFORMATION (ORTHOGONAL ROTATION MATRIX)
+
+ Rotation matrix
+ OR[0] OR[1] OR[2]
+ OR[3] OR[4] OR[5]
+ OR[6] OR[7] OR[8]
+
+ OR[9] x shift
+ OR[10] y shift
+ OR[11] z shift
+
+ OR[12] x or global scale
+ OR[13] y scale
+ OR[14] z scale
+
+************************************************************************/
+
+int CRS_compute_georef_equations_or(struct Control_Points_3D *cp,
+ double OR12[], double OR21[])
+{
+ double *tempptr, *OR;
+ int status, i, numactive;
+ struct Control_Points_3D cpc, /* center points */
+ cpr; /* reduced to center */
+
+ cpc.count = 1;
+ cpc.e1 = (double *)G_calloc(cpc.count, sizeof(double));
+ cpc.e2 = (double *)G_calloc(cpc.count, sizeof(double));
+ cpc.n1 = (double *)G_calloc(cpc.count, sizeof(double));
+ cpc.n2 = (double *)G_calloc(cpc.count, sizeof(double));
+ cpc.z1 = (double *)G_calloc(cpc.count, sizeof(double));
+ cpc.z2 = (double *)G_calloc(cpc.count, sizeof(double));
+ cpc.status = (int *)G_calloc(cpc.count, sizeof(int));
+
+ cpc.e1[0] = 0.0;
+ cpc.e2[0] = 0.0;
+ cpc.n1[0] = 0.0;
+ cpc.n2[0] = 0.0;
+ cpc.z1[0] = 0.0;
+ cpc.z2[0] = 0.0;
+ cpc.status[0] = 1;
+
+ /* center points */
+ for (i = numactive = 0; i < cp->count; i++) {
+ if (cp->status[i] > 0) {
+ numactive++;
+ cpc.e1[0] += cp->e1[i];
+ cpc.e2[0] += cp->e2[i];
+ cpc.n1[0] += cp->n1[i];
+ cpc.n2[0] += cp->n2[i];
+ cpc.z1[0] += cp->z1[i];
+ cpc.z2[0] += cp->z2[i];
+ }
+ }
+
+ /* this version of 3D transformation needs 3 control points */
+ if (numactive < 3)
+ return MNPTERR;
+
+ cpc.e1[0] /= numactive;
+ cpc.e2[0] /= numactive;
+ cpc.n1[0] /= numactive;
+ cpc.n2[0] /= numactive;
+ cpc.z1[0] /= numactive;
+ cpc.z2[0] /= numactive;
+
+ /* shift to center points */
+ cpr.count = numactive;
+ cpr.e1 = (double *)G_calloc(cpr.count, sizeof(double));
+ cpr.e2 = (double *)G_calloc(cpr.count, sizeof(double));
+ cpr.n1 = (double *)G_calloc(cpr.count, sizeof(double));
+ cpr.n2 = (double *)G_calloc(cpr.count, sizeof(double));
+ cpr.z1 = (double *)G_calloc(cpr.count, sizeof(double));
+ cpr.z2 = (double *)G_calloc(cpr.count, sizeof(double));
+ cpr.status = (int *)G_calloc(cpr.count, sizeof(int));
+
+ for (i = numactive = 0; i < cp->count; i++) {
+ if (cp->status[i] > 0) {
+ cpr.e1[numactive] = cp->e1[i] - cpc.e1[0];
+ cpr.e2[numactive] = cp->e2[i] - cpc.e2[0];
+ cpr.n1[numactive] = cp->n1[i] - cpc.n1[0];
+ cpr.n2[numactive] = cp->n2[i] - cpc.n2[0];
+ cpr.z1[numactive] = cp->z1[i] - cpc.z1[0];
+ cpr.z2[numactive] = cp->z2[i] - cpc.z2[0];
+ cpr.status[numactive] = 1;
+ numactive++;
+ }
+ }
+
+ /* CALCULATE THE FORWARD TRANSFORMATION COEFFICIENTS */
+
+ OR = OR12;
+ status = calccoef(&cpr, OR, 3);
+
+ calcscale(&cpr, OR);
+
+ /* CALCULATE FORWARD SHIFTS */
+
+ OR[9] = OR[10] = OR[11] = 0.0;
+
+ for (i = numactive = 0; i < cp->count; i++) {
+ if (cp->status[i] > 0) {
+ OR[9] += cp->e2[i] - OR[12] *
+ (OR[0] * cp->e1[i] +
+ OR[1] * cp->n1[i] +
+ OR[2] * cp->z1[i]);
+ OR[10] += cp->n2[i] - OR[13] *
+ (OR[3] * cp->e1[i] +
+ OR[4] * cp->n1[i] +
+ OR[5] * cp->z1[i]);
+ OR[11] += cp->z2[i] - OR[14] *
+ (OR[6] * cp->e1[i] +
+ OR[7] * cp->n1[i] +
+ OR[8] * cp->z1[i]);
+
+ numactive++;
+ }
+ }
+ OR[9] /= numactive;
+ OR[10] /= numactive;
+ OR[11] /= numactive;
+
+ /* SWITCH THE 1 AND 2 EASTING, NORTHING, AND HEIGHT ARRAYS */
+
+ tempptr = cpr.e1;
+ cpr.e1 = cpr.e2;
+ cpr.e2 = tempptr;
+ tempptr = cpr.n1;
+ cpr.n1 = cpr.n2;
+ cpr.n2 = tempptr;
+ tempptr = cpr.z1;
+ cpr.z1 = cpr.z2;
+ cpr.z2 = tempptr;
+
+ /* CALCULATE THE BACKWARD TRANSFORMATION COEFFICIENTS */
+
+ OR = OR21;
+ status = calccoef(&cpr, OR, 3);
+
+ if (status != MSUCCESS)
+ return status;
+
+ calcscale(&cpr, OR);
+
+ /* SWITCH THE 1 AND 2 EASTING, NORTHING, AND HEIGHT ARRAYS BACK */
+
+ tempptr = cpr.e1;
+ cpr.e1 = cpr.e2;
+ cpr.e2 = tempptr;
+ tempptr = cpr.n1;
+ cpr.n1 = cpr.n2;
+ cpr.n2 = tempptr;
+ tempptr = cpr.z1;
+ cpr.z1 = cpr.z2;
+ cpr.z2 = tempptr;
+
+ /* CALCULATE BACKWARD SHIFTS */
+
+ OR[9] = OR[10] = OR[11] = 0.0;
+
+ for (i = numactive = 0; i < cp->count; i++) {
+ if (cp->status[i] > 0) {
+ OR[9] += cp->e1[i] - OR[12] *
+ (OR[0] * cp->e2[i] +
+ OR[1] * cp->n2[i] +
+ OR[2] * cp->z2[i]);
+ OR[10] += cp->n1[i] - OR[13] *
+ (OR[3] * cp->e2[i] +
+ OR[4] * cp->n2[i] +
+ OR[5] * cp->z2[i]);
+ OR[11] += cp->z1[i] - OR[14] *
+ (OR[6] * cp->e2[i] +
+ OR[7] * cp->n2[i] +
+ OR[8] * cp->z2[i]);
+ numactive++;
+ }
+ }
+ OR[9] /= numactive;
+ OR[10] /= numactive;
+ OR[11] /= numactive;
+
+ OR = OR12;
+ G_debug(1, "********************************");
+ G_debug(1, "Forward transformation:");
+ G_debug(1, "Orthogonal rotation matrix:");
+ G_debug(1, "%.4f %.4f %.4f", OR[0], OR[1], OR[2]);
+ G_debug(1, "%.4f %.4f %.4f", OR[3], OR[4], OR[5]);
+ G_debug(1, "%.4f %.4f %.4f", OR[6], OR[7], OR[8]);
+ G_debug(1, "x, y, z shift: %.4f %.4f %.4f", OR[9], OR[10], OR[11]);
+ G_debug(1, "x, y, z scale: %.4f %.4f %.4f", OR[12], OR[13], OR[14]);
+
+ OR = OR21;
+ G_debug(1, "********************************");
+ G_debug(1, "Backward transformation:");
+ G_debug(1, "Orthogonal rotation matrix:");
+ G_debug(1, "%.4f %.4f %.4f", OR[0], OR[1], OR[2]);
+ G_debug(1, "%.4f %.4f %.4f", OR[3], OR[4], OR[5]);
+ G_debug(1, "%.4f %.4f %.4f", OR[6], OR[7], OR[8]);
+ G_debug(1, "x, y, z shift: %.4f %.4f %.4f", OR[9], OR[10], OR[11]);
+ G_debug(1, "x, y, z scale: %.4f %.4f %.4f", OR[12], OR[13], OR[14]);
+
+ return status;
+}
+
+/***********************************************************************
+
+ COMPUTE THE GEOREFFERENCING COEFFICIENTS
+ BASED ON A SET OF CONTROL POINTS
+
+************************************************************************/
+
+static int calccoef(struct Control_Points_3D *cp, double OR[], int ndims)
+{
+ double **src_mat = NULL;
+ double **src_mat_T = NULL;
+ double **dest_mat = NULL;
+ double **dest_mat_T = NULL;
+ double *S_vec = NULL;
+ double **R_mat = NULL;
+ double **R_mat_T = NULL;
+ double **C_mat = NULL;
+ double **C_mat_interm = NULL;
+ double **D_mat_interm = NULL;
+ double **E_mat = NULL;
+ double **P_mat = NULL;
+ double **P_mat_T = NULL;
+ double **Q_mat = NULL;
+ double *D_vec = NULL;
+ double *one_vec = NULL;
+ double trace1 = 0.0;
+ double trace2 = 0.0;
+ int numactive; /* NUMBER OF ACTIVE CONTROL POINTS */
+ int i, j;
+ int status;
+
+ /* CALCULATE THE NUMBER OF VALID CONTROL POINTS */
+
+ for (i = numactive = 0; i < cp->count; i++) {
+ if (cp->status[i] > 0)
+ numactive++;
+ }
+
+ src_mat = G_alloc_matrix(numactive, numactive);
+ dest_mat = G_alloc_matrix(numactive, numactive);
+
+ for (i = numactive = 0; i < cp->count; i++) {
+ if (cp->status[i] > 0) {
+ src_mat[numactive][0] = cp->e1[i];
+ src_mat[numactive][1] = cp->n1[i];
+ src_mat[numactive][2] = cp->z1[i];
+
+ dest_mat[numactive][0] = cp->e2[i];
+ dest_mat[numactive][1] = cp->n2[i];
+ dest_mat[numactive][2] = cp->z2[i];
+
+ numactive++;
+ }
+ }
+
+ D_vec = G_alloc_vector(ndims);
+
+ src_mat_T = G_alloc_matrix(numactive, numactive);
+ dest_mat_T = G_alloc_matrix(numactive, numactive);
+ R_mat = G_alloc_matrix(ndims, ndims);
+ R_mat_T = G_alloc_matrix(numactive, numactive);
+
+ C_mat = G_alloc_matrix(numactive, numactive);
+ C_mat_interm = G_alloc_matrix(numactive, numactive);
+ D_mat_interm = G_alloc_matrix(numactive, numactive);
+ E_mat = G_alloc_matrix(numactive, numactive);
+ P_mat = G_alloc_matrix(numactive, numactive);
+ P_mat_T = G_alloc_matrix(numactive, numactive);
+ Q_mat = G_alloc_matrix(numactive, numactive);
+
+ transpose_matrix(numactive, dest_mat, dest_mat_T);
+
+ for (i = 0; i < numactive; i++) {
+ for (j = 0; j < numactive; j++) {
+ if (i != j) {
+ E_mat[i][j] = -1.0 / (double)numactive;
+ }
+ else{
+ E_mat[i][j] = 1.0 - 1.0 / (double)numactive;
+ }
+ }
+ }
+
+ matmult(numactive, numactive, dest_mat_T, E_mat, C_mat_interm);
+ matmult(ndims, numactive, C_mat_interm, src_mat, C_mat);
+ copy_matrix(ndims, ndims, C_mat, P_mat);
+
+#if 1
+ svd(ndims, ndims, C_mat, P_mat, D_vec, Q_mat);
+ status = MSUCCESS;
+
+#else
+ /* not working... */
+ status = G_math_svduv(D_vec, C_mat, Q_mat, ndims, P_mat, ndims);
+
+ if (status == 0)
+ status = MSUCCESS;
+#endif
+
+ transpose_matrix(ndims, P_mat, P_mat_T);
+
+ /* rotation matrix */
+ matmult(ndims, ndims, Q_mat, P_mat_T, R_mat_T);
+ transpose_matrix(ndims, R_mat_T, R_mat);
+
+ /* scale */
+ matmult(numactive, ndims, C_mat, R_mat_T, C_mat_interm);
+ trace1 = trace(ndims, ndims, C_mat_interm);
+
+ transpose_matrix(numactive, src_mat, src_mat_T);
+ matmult(numactive, numactive, src_mat_T, E_mat, C_mat_interm);
+ matmult(ndims, numactive, C_mat_interm, src_mat, C_mat);
+ trace2 = trace(ndims, ndims, C_mat);
+
+ OR[14] = trace1 / trace2;
+
+ /* shifts */
+ matmult(numactive, ndims, src_mat, R_mat_T, D_mat_interm);
+ scale_matrix(numactive, ndims, OR[14], D_mat_interm, C_mat_interm);
+ subtract_matrix(numactive, ndims, dest_mat, C_mat_interm, D_mat_interm);
+ scale_matrix(numactive, ndims, 1.0 / numactive, D_mat_interm, C_mat_interm);
+ transpose_matrix(numactive, C_mat_interm, D_mat_interm);
+
+ S_vec = G_alloc_vector(ndims);
+ one_vec = G_alloc_vector(numactive);
+
+ for (i = 0; i < numactive; i++){
+ one_vec[i] = 1.0;
+ }
+
+ matrix_multiply(ndims, numactive, D_mat_interm, one_vec, S_vec);
+
+ /* matrix to vector */
+ for (i = 0; i < ndims; i++) {
+ for (j = 0; j < ndims; j++) {
+ OR[i * ndims + j] = R_mat[i][j];
+ }
+ }
+
+ G_free_matrix(src_mat);
+ G_free_matrix(src_mat_T);
+ G_free_matrix(dest_mat);
+ G_free_matrix(dest_mat_T);
+ G_free_vector(D_vec);
+ G_free_matrix(E_mat);
+ G_free_matrix(P_mat);
+ G_free_matrix(P_mat_T);
+ G_free_matrix(Q_mat);
+ G_free_matrix(R_mat);
+ G_free_matrix(R_mat_T);
+ G_free_matrix(C_mat);
+ G_free_matrix(C_mat_interm);
+ G_free_matrix(D_mat_interm);
+ G_free_vector(S_vec);
+ G_free_vector(one_vec);
+
+ return status;
+}
+
+/***********************************************************************
+
+ CALCULATE SCALE
+
+************************************************************************/
+static int calcscale(struct Control_Points_3D *cp, double OR[])
+{
+ double sumX, sumY, sumsqX, sumsqY, sumXY;
+ int numactive; /* NUMBER OF ACTIVE CONTROL POINTS */
+ int i;
+
+ /* CALCULATE SCALE */
+ sumX = sumY = sumsqX = sumsqY = sumXY = 0.0;
+
+ for (i = numactive = 0; i < cp->count; i++) {
+ if (cp->status[i] > 0) {
+ double c1, c2;
+
+ c1 = OR[0] * cp->e1[i] + OR[1] * cp->n1[i] + OR[2] * cp->z1[i];
+ c2 = cp->e2[i];
+
+ sumX += c1;
+ sumY += c2;
+ sumsqX += c1 * c1;
+ sumsqY += c2 * c2;
+ sumXY += c1 * c2;
+
+ c1 = OR[3] * cp->e1[i] + OR[4] * cp->n1[i] + OR[5] * cp->z1[i];
+ c2 = cp->n2[i];
+
+ sumX += c1;
+ sumY += c2;
+ sumsqX += c1 * c1;
+ sumsqY += c2 * c2;
+ sumXY += c1 * c2;
+
+ c1 = OR[6] * cp->e1[i] + OR[7] * cp->n1[i] + OR[8] * cp->z1[i];
+ c2 = cp->z2[i];
+
+ sumX += c1;
+ sumY += c2;
+ sumsqX += c1 * c1;
+ sumsqY += c2 * c2;
+ sumXY += c1 * c2;
+
+ numactive++;
+ }
+ }
+
+ OR[12] = (sumXY - sumX * sumY / numactive) /
+ (sumsqX - sumX * sumX / numactive);
+
+ if (fabs(OR[12] - OR[14]) > 10 * GRASS_EPSILON) {
+ G_debug(1, "Scale mismatch: %.4f %.4f", OR[12], OR[14]);
+ OR[12] = OR[14];
+ }
+
+ OR[13] = OR[14] = OR[12];
+
+ return MSUCCESS;
+}
+
+void transpose_matrix(int n, double **src_matrix, double **dest_matrix)
+{
+ int i, j;
+
+ for(i = 0; i < n; i++) {
+ for(j = 0; j < n; j++) {
+ dest_matrix[i][j] = src_matrix[j][i];
+ }
+ }
+}
+
+void matmult(int m, int n, double **mat1, double **mat2, double **mat3)
+{
+ int i, j, k;
+ double sum;
+
+ for (i = 0; i < m; i++) {
+ for (j = 0; j < n; j++) {
+ sum = 0.0;
+ for (k = 0; k < n; k++) {
+ sum += mat1[i][k] * mat2[k][j];
+ }
+ mat3[i][j] = sum;
+ }
+ }
+}
+
+void copy_matrix(int n, int m, double **src_matrix, double **dest_matrix)
+{
+ int i, j;
+
+ for (i = 0; i < n; i++) {
+ for (j = 0; j < m; j++) {
+ dest_matrix[i][j] = src_matrix[i][j];
+ }
+ }
+}
+
+double trace(int n, int m, double **matrix)
+{
+ int i;
+ double t = 0.0;
+
+ for (i = 0; i < n; i++) {
+ t += matrix[i][i];
+ }
+
+ return t;
+}
+
+void init_matrix(int n, int m, double **matrix)
+{
+ int i, j;
+
+ for (i = 0; i < n; i++) {
+ for (j = 0; j < m; j++) {
+ matrix[i][j] = 0.0;
+ }
+ }
+}
+
+void scale_matrix(int n, int m, double scal, double **src_matrix,
+ double **dest_matrix)
+{
+ int i, j;
+
+ for (i = 0; i < n; i++) {
+ for (j = 0; j < m; j++) {
+ dest_matrix[i][j] = scal * src_matrix[i][j];
+ }
+ }
+}
+
+void subtract_matrix(int n, int m, double **mat1, double **mat2,
+ double **mat3)
+{
+ int i, j;
+
+ for (i = 0; i < n; i++) {
+ for (j = 0; j < m; j++) {
+ mat3[i][j] = mat1[i][j] - mat2[i][j];
+ }
+ }
+}
+
+void matrix_multiply(int n, int m, double **mat, double *iv,
+ double *ov)
+{
+ int i, j;
+
+ for (i = 0; i < n; i++) {
+ ov[i] = 0.0;
+ for(j = 0; j < m; j++) {
+ ov[i] += mat[i][j] * iv[j];
+ }
+ }
+}
+
Property changes on: grass/trunk/vector/v.rectify/orthorot.c
___________________________________________________________________
Added: svn:mime-type
+ text/x-csrc
Added: svn:eol-style
+ native
Added: grass/trunk/vector/v.rectify/svdm.c
===================================================================
--- grass/trunk/vector/v.rectify/svdm.c (rev 0)
+++ grass/trunk/vector/v.rectify/svdm.c 2012-04-26 06:54:29 UTC (rev 51546)
@@ -0,0 +1,291 @@
+/*
+ Copyright (c) 2003, Division of Imaging Science and Biomedical Engineering,
+ University of Manchester, UK. All rights reserved.
+
+ Redistribution and use in source and binary forms, with or without modification,
+ are permitted provided that the following conditions are met:
+
+ Redistributions of source code must retain the above copyright notice, this list
+ of conditions and the following disclaimer.
+
+ Redistributions in binary form must reproduce the above copyright notice, this
+ list of conditions and the following disclaimer in the documentation and/or other
+ materials provided with the distribution.
+
+ Neither the name of the University of Manchester nor the names of its contributors
+ may be used to endorse or promote products derived from this software without
+ specific prior written permission.
+
+
+ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY
+ EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
+ OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT
+ SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+ SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+ PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+ LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+ OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ */
+
+#include <math.h>
+#include <stdlib.h>
+
+/* svd function */
+#define SIGN(u, v) ( (v)>=0.0 ? fabs(u) : -fabs(u) )
+#define MAX(x, y) ( (x) >= (y) ? (x) : (y) )
+
+static double radius(double u, double v)
+{
+ double w;
+
+ u = fabs(u);
+ v = fabs(v);
+ if (u > v) {
+ w = v / u;
+ return (u * sqrt(1. + w * w));
+ }
+ else {
+ if (v) {
+ w = u / v;
+ return (v * sqrt(1. + w * w));
+ }
+ else
+ return 0.0;
+ }
+}
+
+/*
+ Given matrix a[m][n], m>=n, using svd decomposition a = p d q' to get
+ p[m][n], diag d[n] and q[n][n].
+ */
+void svd(int m, int n, double **a, double **p, double *d, double **q)
+{
+ int flag, i, its, j, jj, k, l, nm, nm1 = n - 1, mm1 = m - 1;
+ double c, f, h, s, x, y, z;
+ double anorm = 0, g = 0, scale = 0;
+ double *r = (double *)malloc(sizeof(double) * n);
+
+ for (i = 0; i < m; i++)
+ for (j = 0; j < n; j++)
+ p[i][j] = a[i][j];
+ /*
+ for (i = m; i < n; i++)
+ p[i][j] = 0;
+ */
+
+ /* Householder reduction to bidigonal form */
+ for (i = 0; i < n; i++) {
+ l = i + 1;
+ r[i] = scale * g;
+ g = s = scale = 0.0;
+ if (i < m) {
+ for (k = i; k < m; k++)
+ scale += fabs(p[k][i]);
+ if (scale) {
+ for (k = i; k < m; k++) {
+ p[k][i] /= scale;
+ s += p[k][i] * p[k][i];
+ }
+ f = p[i][i];
+ g = -SIGN(sqrt(s), f);
+ h = f * g - s;
+ p[i][i] = f - g;
+ if (i != nm1) {
+ for (j = l; j < n; j++) {
+ for (s = 0.0, k = i; k < m; k++)
+ s += p[k][i] * p[k][j];
+ f = s / h;
+ for (k = i; k < m; k++)
+ p[k][j] += f * p[k][i];
+ }
+ }
+ for (k = i; k < m; k++)
+ p[k][i] *= scale;
+ }
+ }
+ d[i] = scale * g;
+ g = s = scale = 0.0;
+ if (i < m && i != nm1) {
+ for (k = l; k < n; k++)
+ scale += fabs(p[i][k]);
+ if (scale) {
+ for (k = l; k < n; k++) {
+ p[i][k] /= scale;
+ s += p[i][k] * p[i][k];
+ }
+ f = p[i][l];
+ g = -SIGN(sqrt(s), f);
+ h = f * g - s;
+ p[i][l] = f - g;
+ for (k = l; k < n; k++)
+ r[k] = p[i][k] / h;
+ if (i != mm1) {
+ for (j = l; j < m; j++) {
+ for (s = 0.0, k = l; k < n; k++)
+ s += p[j][k] * p[i][k];
+ for (k = l; k < n; k++)
+ p[j][k] += s * r[k];
+ }
+ }
+ for (k = l; k < n; k++)
+ p[i][k] *= scale;
+ }
+ }
+ anorm = MAX(anorm, fabs(d[i]) + fabs(r[i]));
+ }
+
+ /* Accumulation of right-hand transformations */
+ for (i = n - 1; i >= 0; i--) {
+ if (i < nm1) {
+ if (g) {
+ for (j = l; j < n; j++)
+ q[j][i] = (p[i][j] / p[i][l]) / g;
+ for (j = l; j < n; j++) {
+ for (s = 0.0, k = l; k < n; k++)
+ s += p[i][k] * q[k][j];
+ for (k = l; k < n; k++)
+ q[k][j] += s * q[k][i];
+ }
+ }
+ for (j = l; j < n; j++)
+ q[i][j] = q[j][i] = 0.0;
+ }
+ q[i][i] = 1.0;
+ g = r[i];
+ l = i;
+ }
+ /* Accumulation of left-hand transformations */
+ for (i = n - 1; i >= 0; i--) {
+ l = i + 1;
+ g = d[i];
+ if (i < nm1)
+ for (j = l; j < n; j++)
+ p[i][j] = 0.0;
+ if (g) {
+ g = 1.0 / g;
+ if (i != nm1) {
+ for (j = l; j < n; j++) {
+ for (s = 0.0, k = l; k < m; k++)
+ s += p[k][i] * p[k][j];
+ f = (s / p[i][i]) * g;
+ for (k = i; k < m; k++)
+ p[k][j] += f * p[k][i];
+ }
+ }
+ for (j = i; j < m; j++)
+ p[j][i] *= g;
+ }
+ else
+ for (j = i; j < m; j++)
+ p[j][i] = 0.0;
+ ++p[i][i];
+ }
+ /* diagonalization of the bidigonal form */
+ for (k = n - 1; k >= 0; k--) { /* loop over singlar values */
+ for (its = 0; its < 30; its++) { /* loop over allowed iterations */
+ flag = 1;
+ for (l = k; l >= 0; l--) { /* test for splitting */
+ nm = l - 1; /* note that r[l] is always
+ * zero */
+ if (fabs(r[l]) + anorm == anorm) {
+ flag = 0;
+ break;
+ }
+ if (fabs(d[nm]) + anorm == anorm)
+ break;
+ }
+ if (flag) {
+ c = 0.0; /* cancellation of r[l], if
+ * l>1 */
+ s = 1.0;
+ for (i = l; i <= k; i++) {
+ f = s * r[i];
+ if (fabs(f) + anorm != anorm) {
+ g = d[i];
+ h = radius(f, g);
+ d[i] = h;
+ h = 1.0 / h;
+ c = g * h;
+ s = (-f * h);
+ for (j = 0; j < m; j++) {
+ y = p[j][nm];
+ z = p[j][i];
+ p[j][nm] = y * c + z * s;
+ p[j][i] = z * c - y * s;
+ }
+ }
+ }
+ }
+ z = d[k];
+ if (l == k) { /* convergence */
+ if (z < 0.0) {
+ d[k] = -z;
+ for (j = 0; j < n; j++)
+ q[j][k] = (-q[j][k]);
+ }
+ break;
+ }
+ if (its == 100) {
+ /* svd: No convergence in 100 svd iterations, non_fatal error */
+ return;
+ }
+ x = d[l]; /* shift from bottom 2-by-2 minor */
+ nm = k - 1;
+ y = d[nm];
+ g = r[nm];
+ h = r[k];
+ f = ((y - z) * (y + z) + (g - h) * (g + h)) / (2.0 * h * y);
+ g = radius(f, 1.0);
+ /* next QR transformation */
+ f = ((x - z) * (x + z) + h * ((y / (f + SIGN(g, f))) - h)) / x;
+ c = s = 1.0;
+ for (j = l; j <= nm; j++) {
+ i = j + 1;
+ g = r[i];
+ y = d[i];
+ h = s * g;
+ g = c * g;
+ z = radius(f, h);
+ r[j] = z;
+ c = f / z;
+ s = h / z;
+ f = x * c + g * s;
+ g = g * c - x * s;
+ h = y * s;
+ y = y * c;
+ for (jj = 0; jj < n; jj++) {
+ x = q[jj][j];
+ z = q[jj][i];
+ q[jj][j] = x * c + z * s;
+ q[jj][i] = z * c - x * s;
+ }
+ z = radius(f, h);
+ d[j] = z; /* rotation can be arbitrary
+ * id z=0 */
+ if (z) {
+ z = 1.0 / z;
+ c = f * z;
+ s = h * z;
+ }
+ f = (c * g) + (s * y);
+ x = (c * y) - (s * g);
+ for (jj = 0; jj < m; jj++) {
+ y = p[jj][j];
+ z = p[jj][i];
+ p[jj][j] = y * c + z * s;
+ p[jj][i] = z * c - y * s;
+ }
+ }
+ r[l] = 0.0;
+ r[k] = f;
+ d[k] = x;
+ }
+ }
+ free(r);
+
+ /* dhli add: the original code does not sort the eigen value
+ should do that and change the eigen vector accordingly
+ */
+
+}
Property changes on: grass/trunk/vector/v.rectify/svdm.c
___________________________________________________________________
Added: svn:mime-type
+ text/x-csrc
Added: svn:eol-style
+ native
Modified: grass/trunk/vector/v.rectify/target.c
===================================================================
--- grass/trunk/vector/v.rectify/target.c 2012-04-25 12:50:56 UTC (rev 51545)
+++ grass/trunk/vector/v.rectify/target.c 2012-04-26 06:54:29 UTC (rev 51546)
@@ -1,5 +1,8 @@
#include <string.h>
+#include <stdio.h>
#include <unistd.h>
+#include <grass/gis.h>
+#include <grass/imagery.h>
#include <grass/glocale.h>
#include "global.h"
Modified: grass/trunk/vector/v.rectify/v.rectify.html
===================================================================
--- grass/trunk/vector/v.rectify/v.rectify.html 2012-04-25 12:50:56 UTC (rev 51545)
+++ grass/trunk/vector/v.rectify/v.rectify.html 2012-04-26 06:54:29 UTC (rev 51546)
@@ -6,6 +6,11 @@
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>
+The <em>-o</em> flag enforces orthogonal rotation (currently for 3D only)
+where the axes remain orthogonal to each other, e.g. a cube with right
+angles remains a cube with right angles after tranformation. This is not
+guaranteed even with affine (1<sup>st</sup> order) 3D transformation.
<p>
Great care should be taken with the placement of Ground Control Points.
More information about the grass-commit
mailing list