[GRASS-SVN] r55088 - in grass/trunk/imagery/i.ortho.photo: . i.ortho.transform
svn_grass at osgeo.org
svn_grass at osgeo.org
Sun Feb 17 12:50:13 PST 2013
Author: mmetz
Date: 2013-02-17 12:50:13 -0800 (Sun, 17 Feb 2013)
New Revision: 55088
Added:
grass/trunk/imagery/i.ortho.photo/i.ortho.transform/
grass/trunk/imagery/i.ortho.photo/i.ortho.transform/i.ortho.transform.html
Removed:
grass/trunk/imagery/i.ortho.photo/i.ortho.transform/m.transform.html
Modified:
grass/trunk/imagery/i.ortho.photo/i.ortho.transform/Makefile
grass/trunk/imagery/i.ortho.photo/i.ortho.transform/main.c
Log:
add i.ortho.transform
Modified: grass/trunk/imagery/i.ortho.photo/i.ortho.transform/Makefile
===================================================================
--- grass/trunk/misc/m.transform/Makefile 2013-02-10 00:49:30 UTC (rev 54992)
+++ grass/trunk/imagery/i.ortho.photo/i.ortho.transform/Makefile 2013-02-17 20:50:13 UTC (rev 55088)
@@ -1,10 +1,12 @@
-MODULE_TOPDIR = ../..
+MODULE_TOPDIR = ../../..
-PGM = m.transform
+PGM = i.ortho.transform
-LIBES = $(IMAGERYLIB) $(GISLIB) $(MATHLIB)
-DEPENDENCIES = $(IMAGERYDEP) $(GISDEP)
+EXTRA_CFLAGS = -I../lib
+LIBES = $(IMAGERYLIB) $(IORTHOLIB) $(GISLIB) $(MATHLIB)
+DEPENDENCIES = $(IMAGERYDEP) $(IORTHODEP) $(GISDEP)
+
include $(MODULE_TOPDIR)/include/Make/Module.make
default: cmd
Copied: grass/trunk/imagery/i.ortho.photo/i.ortho.transform/i.ortho.transform.html (from rev 54992, grass/trunk/misc/m.transform/m.transform.html)
===================================================================
--- grass/trunk/imagery/i.ortho.photo/i.ortho.transform/i.ortho.transform.html (rev 0)
+++ grass/trunk/imagery/i.ortho.photo/i.ortho.transform/i.ortho.transform.html 2013-02-17 20:50:13 UTC (rev 55088)
@@ -0,0 +1,63 @@
+<h2>DESCRIPTION</h2>
+
+<em>m.transform</em> is an utility to compute transformation
+based upon GCPs and output error measurements.
+
+
+<h2>NOTES</h2>
+
+For coordinates given with the <b>coords</b> file option or fed from
+<tt>stdin</tt>, the input format is "x y" with one coordinate pair per
+line.
+
+<p>The transformations are:
+<p>order=1:
+<div class="code"><pre>
+ e = [E0 E1][1].[1]
+ [E2 0][e] [n]
+
+ n = [N0 N1][1].[1]
+ [N2 0][e] [n]
+</pre></div>
+
+order=2:
+<div class="code"><pre>
+ e = [E0 E1 E3][1 ] [1 ]
+ [E2 E4 0][e ].[n ]
+ [E5 0 0][e²] [n²]
+
+ n = [N0 N1 N3][1 ] [1 ]
+ [N2 N4 0][e ].[n ]
+ [N5 0 0][e²] [n²]
+</pre></div>
+
+order=3:
+<div class="code"><pre>
+ e = [E0 E1 E3 E6][1 ] [1 ]
+ [E2 E4 E7 0][e ].[n ]
+ [E5 E8 0 0][e²] [n²]
+ [E9 0 0 0][e³] [n³]
+
+ n = [N0 N1 N3 N6][1 ] [1 ]
+ [N2 N4 N7 0][e ].[n ]
+ [N5 N8 0 0][e²] [n²]
+ [N9 0 0 0][e³] [n³]
+</pre></div>
+
+["." = dot-product, (AE).N = N'EA.]
+<p>In other words, order=1 and order=2 are equivalent to order=3 with
+the higher coefficients equal to zero.
+
+
+<h2>SEE ALSO</h2>
+
+<em><a href="i.rectify.html">i.rectify</a></em>
+
+
+<h2>AUTHORS</h2>
+
+Brian J. Buckley<br>
+Glynn Clements<br>
+Hamish Bowman
+
+<p><i>Last changed: $Date$</i>
Deleted: grass/trunk/imagery/i.ortho.photo/i.ortho.transform/m.transform.html
===================================================================
--- grass/trunk/misc/m.transform/m.transform.html 2013-02-10 00:49:30 UTC (rev 54992)
+++ grass/trunk/imagery/i.ortho.photo/i.ortho.transform/m.transform.html 2013-02-17 20:50:13 UTC (rev 55088)
@@ -1,63 +0,0 @@
-<h2>DESCRIPTION</h2>
-
-<em>m.transform</em> is an utility to compute transformation
-based upon GCPs and output error measurements.
-
-
-<h2>NOTES</h2>
-
-For coordinates given with the <b>coords</b> file option or fed from
-<tt>stdin</tt>, the input format is "x y" with one coordinate pair per
-line.
-
-<p>The transformations are:
-<p>order=1:
-<div class="code"><pre>
- e = [E0 E1][1].[1]
- [E2 0][e] [n]
-
- n = [N0 N1][1].[1]
- [N2 0][e] [n]
-</pre></div>
-
-order=2:
-<div class="code"><pre>
- e = [E0 E1 E3][1 ] [1 ]
- [E2 E4 0][e ].[n ]
- [E5 0 0][e²] [n²]
-
- n = [N0 N1 N3][1 ] [1 ]
- [N2 N4 0][e ].[n ]
- [N5 0 0][e²] [n²]
-</pre></div>
-
-order=3:
-<div class="code"><pre>
- e = [E0 E1 E3 E6][1 ] [1 ]
- [E2 E4 E7 0][e ].[n ]
- [E5 E8 0 0][e²] [n²]
- [E9 0 0 0][e³] [n³]
-
- n = [N0 N1 N3 N6][1 ] [1 ]
- [N2 N4 N7 0][e ].[n ]
- [N5 N8 0 0][e²] [n²]
- [N9 0 0 0][e³] [n³]
-</pre></div>
-
-["." = dot-product, (AE).N = N'EA.]
-<p>In other words, order=1 and order=2 are equivalent to order=3 with
-the higher coefficients equal to zero.
-
-
-<h2>SEE ALSO</h2>
-
-<em><a href="i.rectify.html">i.rectify</a></em>
-
-
-<h2>AUTHORS</h2>
-
-Brian J. Buckley<br>
-Glynn Clements<br>
-Hamish Bowman
-
-<p><i>Last changed: $Date$</i>
Modified: grass/trunk/imagery/i.ortho.photo/i.ortho.transform/main.c
===================================================================
--- grass/trunk/misc/m.transform/main.c 2013-02-10 00:49:30 UTC (rev 54992)
+++ grass/trunk/imagery/i.ortho.photo/i.ortho.transform/main.c 2013-02-17 20:50:13 UTC (rev 55088)
@@ -1,13 +1,14 @@
/****************************************************************************
*
- * MODULE: m.transform (nee g.transform)
+ * MODULE: i.ortho.transform (cloned from m.tranform nee g.transform)
* AUTHOR(S): Brian J. Buckley
* Glynn Clements
* Hamish Bowman
+ * Markus Metz
* PURPOSE: Utility to compute transformation based upon GCPs and
* output error measurements
- * COPYRIGHT: (C) 2006-2010 by the GRASS Development Team
+ * COPYRIGHT: (C) 2006-2013 by the GRASS Development Team
*
* This program is free software under the GNU General Public
* License (>=v2). Read the file COPYING that comes with GRASS
@@ -24,6 +25,7 @@
#include <grass/glocale.h>
#include <grass/imagery.h>
#include <grass/glocale.h>
+#include "orthophoto.h"
struct Max
{
@@ -37,8 +39,6 @@
double sum2, rms;
};
-static char *name;
-static int order;
static int summary;
static int forward;
static char **columns;
@@ -48,12 +48,10 @@
static int need_rd;
static char *coord_file;
-static double E12[10], N12[10], E21[10], N21[10];
+struct Ortho_Image_Group group;
-static struct Control_Points points;
+static struct Ortho_Control_Points *points;
-static int equation_stat;
-
static int count;
static struct Stats fwd, rev;
@@ -82,37 +80,88 @@
static void compute_transformation(void)
{
- static const int order_pnts[3] = { 3, 6, 10 };
- int n, i;
+ int n, i, status;
+ double e0, e1, e2, n0, n1, n2, z1, z2;
+ struct Ortho_Control_Points temp_points;
- equation_stat =
- I_compute_georef_equations(&points, E12, N12, E21, N21, order);
+ /* compute photo <-> image equations */
+ group.ref_equation_stat = I_compute_ref_equations(&group.photo_points,
+ group.E12, group.N12,
+ group.E21, group.N21);
- if (equation_stat == 0)
- G_fatal_error(_("Not enough points, %d are required"),
- order_pnts[order - 1]);
+ if (group.ref_equation_stat <= 0)
+ G_fatal_error(_("Error conducting transform (%d)"),
+ group.ref_equation_stat);
- if (equation_stat <= 0)
- G_fatal_error(_("Error conducting transform (%d)"), equation_stat);
+ /* compute target <-> photo equations */
+ /* alloc and fill temp control points */
+ temp_points.count = 0;
+ temp_points.status = NULL;
+ temp_points.e1 = NULL;
+ temp_points.n1 = NULL;
+ temp_points.z1 = NULL;
+ temp_points.e2 = NULL;
+ temp_points.n2 = NULL;
+ temp_points.z2 = NULL;
+
+ /* e0, n0, equal photo coordinates not image coords */
+ for (i = 0; i < group.control_points.count; i++) {
+ status = group.control_points.status[i];
+ e1 = group.control_points.e1[i];
+ n1 = group.control_points.n1[i];
+ z1 = group.control_points.z1[i];
+ e2 = group.control_points.e2[i];
+ n2 = group.control_points.n2[i];
+ z2 = group.control_points.z2[i];
+
+ /* image to photo transformation */
+ I_georef(e1, n1, &e0, &n0, group.E12, group.N12, 1);
+ I_new_con_point(&temp_points, e0, n0, z1, e2, n2, z2, status);
+ }
+
+
+ group.con_equation_stat = I_compute_ortho_equations(&temp_points,
+ &group.camera_ref,
+ &group.camera_exp,
+ &group.XC, &group.YC,
+ &group.ZC,
+ &group.omega,
+ &group.phi,
+ &group.kappa,
+ &group.M,
+ &group.MI);
+
+ if (group.con_equation_stat <= 0)
+ G_fatal_error(_("Error conducting transform (%d)"),
+ group.con_equation_stat);
+
count = 0;
- for (n = 0; n < points.count; n++) {
+ for (n = 0; n < points->count; n++) {
double e1, n1, e2, n2;
double fx, fy, fd, fd2;
double rx, ry, rd, rd2;
- if (points.status[n] <= 0)
+ if (points->status[n] <= 0)
continue;
count++;
if (need_fwd) {
- I_georef(points.e1[n], points.n1[n], &e2, &n2, E12, N12, order);
+ /* image -> photo -> target */
- fx = fabs(e2 - points.e2[n]);
- fy = fabs(n2 - points.n2[n]);
+ /* image coordinates ex, nx to photo coordinates ex1, nx1 */
+ I_georef(points->e1[n], points->n1[n], &e1, &n1, group.E12, group.N12, 1);
+ /* photo coordinates ex1, nx1 to target coordinates e1, n1 */
+ I_inverse_ortho_ref(e1, n1, points->z1[n], &e2, &n2, &z2,
+ &group.camera_ref,
+ group.XC, group.YC, group.ZC, group.MI);
+
+ fx = fabs(e2 - points->e2[n]);
+ fy = fabs(n2 - points->n2[n]);
+
if (need_fd)
diagonal(&fd, &fd2, fx, fy);
@@ -121,11 +170,19 @@
}
if (need_rev) {
- I_georef(points.e2[n], points.n2[n], &e1, &n1, E21, N21, order);
+ /* target -> photo -> image */
- rx = fabs(e1 - points.e1[n]);
- ry = fabs(n1 - points.n1[n]);
+ /* target coordinates e1, n1 to photo coordinates ex1, nx1 */
+ I_ortho_ref(points->e2[n], points->n2[n], points->z2[n],
+ &e2, &n2, &z2, &group.camera_ref,
+ group.XC, group.YC, group.ZC, group.M);
+ /* photo coordinates ex1, nx1 to image coordinates ex, nx */
+ I_georef(e2, n2, &e1, &n1, group.E21, group.N21, 1);
+
+ rx = fabs(e1 - points->e1[n]);
+ ry = fabs(n1 - points->n1[n]);
+
if (need_rd)
diagonal(&rd, &rd2, rx, ry);
@@ -148,9 +205,9 @@
if (strcmp("idx", col) == 0)
printf(" %d", n);
if (strcmp("src", col) == 0)
- printf(" %f %f", points.e1[n], points.n1[n]);
+ printf(" %f %f", points->e1[n], points->n1[n]);
if (strcmp("dst", col) == 0)
- printf(" %f %f", points.e2[n], points.n2[n]);
+ printf(" %f %f", points->e2[n], points->n2[n]);
if (strcmp("fwd", col) == 0)
printf(" %f %f", e2, n2);
if (strcmp("rev", col) == 0)
@@ -190,13 +247,15 @@
static void analyze(void)
{
- if (equation_stat == -1)
- G_warning(_("Poorly placed control points"));
- else if (equation_stat == -2)
+ if (group.ref_equation_stat == -1)
+ G_warning(_("Poorly placed image to photo control points"));
+ else if (group.con_equation_stat == -1)
+ G_warning(_("Poorly placed image to target control points"));
+ else if (group.ref_equation_stat == -2 || group.con_equation_stat == -2)
G_fatal_error(_("Insufficient memory"));
- else if (equation_stat < 0)
+ else if (group.ref_equation_stat < 0 || group.con_equation_stat < 0)
G_fatal_error(_("Parameter error"));
- else if (equation_stat == 0)
+ else if (group.ref_equation_stat == 0 || group.con_equation_stat == 0)
G_fatal_error(_("No active control points"));
else if (summary) {
printf("Number of active points: %d\n", count);
@@ -241,30 +300,43 @@
static void dump_cooefs(void)
{
int i;
- static const int order_pnts[3] = { 3, 6, 10 };
- for (i = 0; i < order_pnts[order - 1]; i++)
- fprintf(stdout, "E%d=%.15g\n", i, forward ? E12[i] : E21[i]);
+ for (i = 0; i < 3; i++)
+ fprintf(stdout, "E%d=%.15g\n", i, forward ? group.E12[i] : group.E21[i]);
- for (i = 0; i < order_pnts[order - 1]; i++)
- fprintf(stdout, "N%d=%.15g\n", i, forward ? N12[i] : N21[i]);
+ for (i = 0; i < 3; i++)
+ fprintf(stdout, "N%d=%.15g\n", i, forward ? group.N12[i] : group.N21[i]);
+
+ /* print ortho transformation matrix ? */
}
-static void xform_value(double east, double north)
+static void xform_value(double east, double north, double height)
{
- double xe, xn;
+ double e1, n1, z1, xe, xn, xz;
- if(forward)
- I_georef(east, north, &xe, &xn, E12, N12, order);
- else
- I_georef(east, north, &xe, &xn, E21, N21, order);
+ if (forward) {
+ /* image -> photo -> target */
+ I_georef(east, north, &e1, &n1, group.E12, group.N12, 1);
+ z1 = height;
+ I_inverse_ortho_ref(e1, n1, z1, &xe, &xn, &xz, &group.camera_ref,
+ group.XC, group.YC, group.ZC, group.MI);
+ xz = z1;
+ }
+ else {
+ /* target -> photo -> image */
+ I_ortho_ref(east, north, height, &e1, &n1, &z1, &group.camera_ref,
+ group.XC, group.YC, group.ZC, group.M);
- fprintf(stdout, "%.15g %.15g\n", xe, xn);
+ I_georef(e1, n1, &xe, &xn, group.E21, group.N21, 1);
+ xz = 0.;
+ }
+
+ fprintf(stdout, "%.15g %.15g %.15g\n", xe, xn, xz);
}
static void do_pt_xforms(void)
{
- double easting, northing;
+ double easting, northing, height;
int ret;
FILE *fp;
@@ -277,7 +349,7 @@
}
for (;;) {
- char buf[64];
+ char buf[1024];
if (!G_getl2(buf, sizeof(buf), fp))
break;
@@ -290,11 +362,11 @@
? G_scan_northing(,,-1) */
/* ? muliple delims with sscanf(buf, "%[ ,|\t]", &dummy) ? */
- ret = sscanf(buf, "%lf %lf", &easting, &northing);
- if (ret != 2)
+ ret = sscanf(buf, "%lf %lf %lf", &easting, &northing, &height);
+ if (ret != 3)
G_fatal_error(_("Invalid coordinates: [%s]"), buf);
- xform_value(easting, northing);
+ xform_value(easting, northing, height);
}
if (fp != stdin)
@@ -304,7 +376,7 @@
int main(int argc, char **argv)
{
- struct Option *grp, *val, *fmt, *xfm_pts;
+ struct Option *grp, *fmt, *xfm_pts;
struct Flag *sum, *rev_flag, *dump_flag;
struct GModule *module;
char *desc;
@@ -313,7 +385,8 @@
/* Get Args */
module = G_define_module();
- G_add_keyword(_("general"));
+ G_add_keyword(_("imagery"));
+ G_add_keyword(_("orthorectify"));
G_add_keyword(_("transformation"));
G_add_keyword(_("GCP"));
module->description =
@@ -321,13 +394,6 @@
grp = G_define_standard_option(G_OPT_I_GROUP);
- val = G_define_option();
- val->key = "order";
- val->type = TYPE_INTEGER;
- val->required = YES;
- val->options = "1-3";
- val->description = _("Rectification polynomial order");
-
fmt = G_define_option();
fmt->key = "format";
fmt->type = TYPE_STRING;
@@ -374,27 +440,35 @@
exit(EXIT_FAILURE);
- name = grp->answer;
- order = atoi(val->answer);
+ G_strip(grp->answer);
+ strcpy(group.name, grp->answer);
+
summary = !!sum->answer;
columns = fmt->answers;
forward = !rev_flag->answer;
coord_file = xfm_pts->answer;
- I_get_control_points(name, &points);
+ if (!I_get_ref_points(group.name, &group.photo_points)) {
+ G_fatal_error(_("Can not read reference points for group <%s>"),
+ group.name);
+ }
+ if (!I_get_con_points(group.name, &group.control_points)) {
+ G_fatal_error(_("Can not read control points for group <%s>"),
+ group.name);
+ }
+
+ points = &group.control_points;
parse_format();
compute_transformation();
- I_put_control_points(name, &points);
-
analyze();
- if(dump_flag->answer)
+ if (dump_flag->answer)
dump_cooefs();
- if(coord_file)
+ if (coord_file)
do_pt_xforms();
return 0;
More information about the grass-commit
mailing list