add i.ortho.elev

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

imagery/i.ortho.photo/Makefile

@@ -2,20 +2,14 @@
 MODULE_TOPDIR = ../..
 
 SUBDIRS1 = \
-	menu \
-	photo.2image \
-	photo.2target \
-	photo.camera \
-	photo.elev \
-	photo.init \
-	photo.rectify \
-	photo.target
+	i.ortho.elev \
+	i.photo.rectify
 
-SUBDIRS = libes $(SUBDIRS1)
+SUBDIRS = lib $(SUBDIRS1)
 
 include $(MODULE_TOPDIR)/include/Make/Dir.make
 
-$(SUBDIRS1): libes
+$(SUBDIRS1): lib
 
 default: parsubdirs
 

imagery/i.ortho.photo/i.ortho.elev/Makefile

@@ -0,0 +1,12 @@
+MODULE_TOPDIR = ../../..
+
+PGM = i.ortho.elev
+
+LIBES     = $(IMAGERYLIB) $(GISLIB) $(IORTHOLIB) $(GMATHLIB)
+DEPENDENCIES= $(IMAGERYDEP) $(IORTHODEP) $(GISDEP) $(GMATHDEP)
+
+EXTRA_CFLAGS = -I../lib
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+default: cmd

imagery/i.ortho.photo/i.ortho.elev/main.c

@@ -0,0 +1,189 @@
+
+/****************************************************************************
+ *
+ * MODULE:       i.photo.elev
+ * AUTHOR(S):    Mike Baba (original contributor)
+ *               Markus Neteler <neteler itc.it>, 
+ *               Roberto Flor <flor itc.it>, 
+ *               Bernhard Reiter <bernhard intevation.de>, 
+ *               Glynn Clements <glynn gclements.plus.com>
+ *               Hamish Bowman
+ *               Markus Metz
+ *
+ * PURPOSE:      Select the elevation model
+ * COPYRIGHT:    (C) 1999-2012 by the GRASS Development Team
+ *
+ *               This program is free software under the GNU General Public
+ *               License (>=v2). Read the file COPYING that comes with GRASS
+ *               for details.
+ *
+ *****************************************************************************/
+
+/* read the target for the block and cast it into the alternate GRASS env */
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <unistd.h>
+#include <string.h>
+#include <grass/gis.h>
+#include <grass/imagery.h>
+#include <grass/glocale.h>
+#include "orthophoto.h"
+
+static int which_env;
+
+int select_current_env(void);
+int select_target_env(void);
+
+int main(int argc, char *argv[])
+{
+
+    struct GModule *module;
+    struct Option *group_opt, *elev_opt;
+    struct Flag *list_flag;
+
+    char location[GMAPSET_MAX];
+    char mapset[GMAPSET_MAX];
+    char group[GNAME_MAX];
+
+    char *elev_layer;
+    char *mapset_elev;
+    char *tl;
+    char *math_exp;
+    char *units;
+    char *nd;
+
+    char buf[100];
+    int stat;
+    int overwrite;
+
+    G_gisinit(argv[0]);
+
+    module = G_define_module();
+    G_add_keyword(_("imagery"));
+    G_add_keyword(_("orthorectify"));
+    module->description =
+	_("Select or modify the target elevation model.");
+
+    group_opt = G_define_standard_option(G_OPT_I_GROUP);
+    group_opt->description =
+	_("Name of imagery group for ortho-rectification");
+
+    elev_opt = G_define_standard_option(G_OPT_R_ELEV);
+    elev_opt->required = NO;
+    elev_opt->description =
+	_("Name of elevation map to use for ortho-rectification");
+
+    list_flag = G_define_flag();
+    list_flag->key = 'l';
+    list_flag->description =
+	_("List available raster maps in target mapset and exit");
+
+    if (G_parser(argc, argv))
+	exit(EXIT_FAILURE);
+
+    elev_layer = (char *)G_malloc(GNAME_MAX * sizeof(char));
+    mapset_elev = (char *)G_malloc(GMAPSET_MAX * sizeof(char));
+    tl = (char *)G_malloc(80 * sizeof(char));
+    math_exp = (char *)G_malloc(80 * sizeof(char));
+    units = (char *)G_malloc(80 * sizeof(char));
+    nd = (char *)G_malloc(80 * sizeof(char));
+
+    *elev_layer = 0;
+    *mapset_elev = 0;
+    *tl = 0;
+    *math_exp = 0;
+    *units = 0;
+    *nd = 0;
+
+    strcpy(group, group_opt->answer);
+
+    G_suppress_warnings(1);
+    if (!I_get_target(group, location, mapset)) {
+	G_fatal_error(_("Please select a target for group [%s] first"), group);
+    }
+
+    G_suppress_warnings(0);
+    sprintf(buf, "%s/%s", G_gisdbase(), location);
+    if (access(buf, 0) != 0) {
+	G_fatal_error(_("Target location [%s] not found\n"), location);
+    }
+
+    G__create_alt_env();
+    G__setenv("LOCATION_NAME", location);
+
+    stat = G__mapset_permissions(mapset);
+    if (stat > 0) {
+	G__setenv("MAPSET", mapset);
+	G__create_alt_search_path();
+	G__switch_env();
+	G__switch_search_path();
+	which_env = 0;
+
+	select_target_env();
+	
+	if (list_flag->answer) {
+	    /* list raster maps in target location */
+
+	    /* select current location */
+	    select_current_env();
+	}
+	else {
+	    if (!elev_opt->answer) {
+		/* select current location */
+		select_current_env();
+		G_fatal_error(_("Elevation map name is missing. Please set '%s' option"),
+		    elev_opt->key);
+	    }
+	    /* elevation map exists in target ? */
+	    
+	    if (G_find_raster2(elev_opt->answer, mapset) == NULL) {
+		/* select current location */
+		select_current_env();
+		G_fatal_error(_("Raster map <%s> not found"), elev_opt->answer);
+	    }
+
+	    I_get_group_elev(group, elev_layer, mapset_elev, tl, math_exp, units, nd);
+	    overwrite = G__getenv("OVERWRITE") != NULL;
+	    
+	    if (!overwrite && *elev_layer) {
+		G_warning(_("Elevation for group <%s> is already set to <%s>"),
+		          group, elev_layer);
+	    }
+
+	    I_put_group_elev(group, elev_opt->answer, mapset_elev, tl,
+			     math_exp, units, nd);
+	    /* select current location */
+	    select_current_env();
+	}
+
+	exit(EXIT_SUCCESS);
+    }
+
+    G_fatal_error(_("Mapset [%s] in target location [%s] - %s "),
+                  mapset, location,
+		  stat == 0 ? _("permission denied\n") : _("not found\n"));
+}
+
+
+int select_current_env(void)
+{
+    if (which_env != 0) {
+	G__switch_env();
+	G__switch_search_path();
+	which_env = 0;
+    }
+
+    return 0;
+}
+
+int select_target_env(void)
+{
+    if (which_env != 1) {
+	G__switch_env();
+	G__switch_search_path();
+	which_env = 1;
+    }
+
+    return 0;
+}

imagery/i.ortho.photo/i.photo.rectify/Makefile

@@ -0,0 +1,15 @@
+MODULE_TOPDIR = ../../..
+
+PGM = i.photo.rectify
+
+# for DEBUG3 see README!
+#EXTRA_CFLAGS = -I../libes -DDEBUG3
+EXTRA_CFLAGS = -I../libes
+
+LIBES     = $(IMAGERYLIB) $(GISLIB) $(IORTHOLIB) $(VASKLIB) $(CURSES) $(GMATHLIB)
+DEPENDENCIES= $(IMAGERYDEP) $(IORTHODEP) $(GISDEP) $(VASKDEP) $(GMATHDEP)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+default: cmd
+

imagery/i.ortho.photo/i.photo.rectify/README

@@ -0,0 +1,13 @@
+1/2002: Markus Neteler
+
+ - have updated write routines to GRASS 5 functions (taken from i.rectify2)
+ - have changed elev reading from CELL to raster_type to support
+   also FP elevations
+
+12/2010: Markus Metz
+
+ - sync'ed to i.rectify as far as possible
+ - the module is now non-interactive and stand-alone: 
+   typing i.photo.rectify on the commandline pops up a gui if GRASS is started with a gui
+ - i.ortho.photo has a new interactive part to set the options for i.photo.rectify 
+

imagery/i.ortho.photo/i.photo.rectify/angle.c

@@ -0,0 +1,219 @@
+/* calculate camera angle to local surface
+ * 90 degrees: orthogonal to local surface
+ * 0 degress: parallel to local surface
+ * < 0 degrees: not visible by camera
+ * 
+ * Earth curvature is not considered, assuming that the extends of the 
+ * imagery to be orthorectified are rather small
+ * Shadowing effects by ridges and peaks are not considered */
+
+#include <stdlib.h>
+#include <string.h>
+#include <unistd.h>
+#include <math.h>
+#include "global.h"
+
+int camera_angle(char *name)
+{
+    int row, col, nrows, ncols;
+    double XC = group.XC;
+    double YC = group.YC;
+    double ZC = group.ZC;
+    double c_angle, c_angle_min, c_alt, c_az, slope, aspect;
+    double radians_to_degrees = 180.0 / M_PI;
+    /* double degrees_to_radians = M_PI / 180.0; */
+    DCELL e1, e2, e3, e4, e5, e6, e7, e8, e9;
+    double factor, V, H, dx, dy, dz, key;
+    double north, south, east, west, ns_med;
+    FCELL *fbuf0, *fbuf1, *fbuf2, *tmpbuf, *outbuf;
+    int elevfd, outfd;
+    struct Cell_head cellhd;
+    struct Colors colr;
+    FCELL clr_min, clr_max;
+    struct History hist;
+    char *type;
+
+    G_message(_("Calculating camera angle to local surface..."));
+    
+    select_target_env();
+    
+    /* align target window to elevation map, otherwise we get artefacts
+     * like in r.slope.aspect -a */
+     
+    if (G_get_cellhd(elev_name, elev_mapset, &cellhd) < 0)
+	return 0;
+
+    G_align_window(&target_window, &cellhd);
+    G_set_window(&target_window);
+    
+    elevfd = G_open_cell_old(elev_name, elev_mapset);
+    if (elevfd < 0) {
+	G_fatal_error(_("Could not open elevation raster"));
+	return 1;
+    }
+
+    nrows = target_window.rows;
+    ncols = target_window.cols;
+    
+    outfd = G_open_raster_new(name, FCELL_TYPE);
+    fbuf0 = G_allocate_raster_buf(FCELL_TYPE);
+    fbuf1 = G_allocate_raster_buf(FCELL_TYPE);
+    fbuf2 = G_allocate_raster_buf(FCELL_TYPE);
+    outbuf = G_allocate_raster_buf(FCELL_TYPE);
+    
+    /* give warning if location units are different from meters and zfactor=1 */
+    factor = G_database_units_to_meters_factor();
+    if (factor != 1.0)
+	G_warning(_("Converting units to meters, factor=%.6f"), factor);
+
+    G_begin_distance_calculations();
+    north = G_row_to_northing(0.5, &target_window);
+    ns_med = G_row_to_northing(1.5, &target_window);
+    south = G_row_to_northing(2.5, &target_window);
+    east = G_col_to_easting(2.5, &target_window);
+    west = G_col_to_easting(0.5, &target_window);
+    V = G_distance(east, north, east, south) * 4;
+    H = G_distance(east, ns_med, west, ns_med) * 4;
+    
+    c_angle_min = 90;
+    G_get_raster_row(elevfd, fbuf1, 0, FCELL_TYPE);
+    G_get_raster_row(elevfd, fbuf2, 1, FCELL_TYPE);
+
+    for (row = 0; row < nrows; row++) {
+	G_percent(row, nrows, 2);
+	
+	G_set_null_value(outbuf, ncols, FCELL_TYPE);
+
+	/* first and last row */
+	if (row == 0 || row == nrows - 1) {
+	    G_put_raster_row(outfd, outbuf, FCELL_TYPE);
+	    continue;
+	}
+	
+	tmpbuf = fbuf0;
+	fbuf0 = fbuf1;
+	fbuf1 = fbuf2;
+	fbuf2 = tmpbuf;
+	
+	G_get_raster_row(elevfd, fbuf2, row + 1, FCELL_TYPE);
+
+	north = G_row_to_northing(row + 0.5, &target_window);
+
+	for (col = 1; col < ncols - 1; col++) {
+	    
+	    e1 = fbuf0[col - 1];
+	    if (G_is_d_null_value(&e1))
+		continue;
+	    e2 = fbuf0[col];
+	    if (G_is_d_null_value(&e2))
+		continue;
+	    e3 = fbuf0[col + 1];
+	    if (G_is_d_null_value(&e3))
+		continue;
+	    e4 = fbuf1[col - 1];
+	    if (G_is_d_null_value(&e4))
+		continue;
+	    e5 = fbuf1[col];
+	    if (G_is_d_null_value(&e5))
+		continue;
+	    e6 = fbuf1[col + 1];
+	    if (G_is_d_null_value(&e6))
+		continue;
+	    e7 = fbuf2[col - 1];
+	    if (G_is_d_null_value(&e7))
+		continue;
+	    e8 = fbuf2[col];
+	    if (G_is_d_null_value(&e8))
+		continue;
+	    e9 = fbuf2[col + 1];
+	    if (G_is_d_null_value(&e9))
+		continue;
+	    
+	    dx = ((e1 + e4 + e4 + e7) - (e3 + e6 + e6 + e9)) / H;
+	    dy = ((e7 + e8 + e8 + e9) - (e1 + e2 + e2 + e3)) / V;
+	    
+	    /* compute topographic parameters */
+	    key = dx * dx + dy * dy;
+	    /* slope in radians */
+	    slope = atan(sqrt(key));
+
+	    /* aspect in radians */
+	    if (key == 0.)
+		aspect = 0.;
+	    else if (dx == 0) {
+		if (dy > 0)
+		    aspect = M_PI / 2;
+		else
+		    aspect = 1.5 * M_PI;
+	    }
+	    else {
+		aspect = atan2(dy, dx);
+		if (aspect <= 0.)
+		    aspect = 2 * M_PI + aspect;
+	    }
+	    
+	    /* camera altitude angle in radians */
+	    east = G_col_to_easting(col + 0.5, &target_window);
+	    dx = east - XC;
+	    dy = north - YC;
+	    dz = ZC - e5;
+	    c_alt = atan(sqrt(dx * dx + dy * dy) / dz);
+
+	    /* camera azimuth angle in radians */
+	    c_az = atan(dy / dx);
+	    if (east < XC && north != YC)
+		c_az += M_PI;
+	    else if (north < YC && east > XC)
+		c_az += 2 * M_PI;
+		
+	    /* camera angle to real ground */
+	    /* orthogonal to ground: 90 degrees */
+	    /* parallel to ground: 0 degrees */
+	    c_angle = asin(cos(c_alt) * cos(slope) - sin(c_alt) * sin(slope) * cos(c_az - aspect));
+	    
+	    outbuf[col] = c_angle * radians_to_degrees;
+	    if (c_angle_min > outbuf[col])
+		c_angle_min = outbuf[col];
+	}
+	G_put_raster_row(outfd, outbuf, FCELL_TYPE);
+    }
+    G_percent(row, nrows, 2);
+
+    G_close_cell(elevfd);
+    G_close_cell(outfd);
+    G_free(fbuf0);
+    G_free(fbuf1);
+    G_free(fbuf2);
+    G_free(outbuf);
+
+    type = "raster";
+    G_short_history(name, type, &hist);
+    G_command_history(&hist);
+    G_write_history(name, &hist);
+    
+    G_init_colors(&colr);
+    if (c_angle_min < 0) {
+	clr_min = (FCELL)((int)(c_angle_min / 10 - 1)) * 10;
+	clr_max = 0;
+	G_add_f_raster_color_rule(&clr_min, 0, 0, 0, &clr_max, 0,
+				  0, 0, &colr);
+    }
+    clr_min = 0;
+    clr_max = 10;
+    G_add_f_raster_color_rule(&clr_min, 0, 0, 0, &clr_max, 255,
+			      0, 0, &colr);
+    clr_min = 10;
+    clr_max = 40;
+    G_add_f_raster_color_rule(&clr_min, 255, 0, 0, &clr_max, 255,
+			      255, 0, &colr);
+    clr_min = 40;
+    clr_max = 90;
+    G_add_f_raster_color_rule(&clr_min, 255, 255, 0, &clr_max, 0,
+			      255, 0, &colr);
+
+    G_write_colors(name, G_mapset(), &colr);
+
+    select_current_env();
+
+    return 1;
+}

imagery/i.ortho.photo/i.photo.rectify/aver_z.c

@@ -0,0 +1,31 @@
+#include "global.h"
+
+int get_aver_elev(struct Ortho_Control_Points *cpz, double *aver_z)
+{
+    double meanz;
+    double *cp = cpz->z2;
+    int i, n;
+
+    /*  Need 1 control points */
+    if (cpz->count <= 0) {
+	return (-1);
+    }
+
+    /* set average elevation from mean values of CONZ points */
+    meanz = 0;
+    n = 0;
+    for (i = 0; i < cpz->count; i++) {
+	if (cpz->status[i] <= 0)
+	    continue;
+
+	n++;
+	meanz += *(cp++);
+	G_debug(3, "In ortho meanz = %f", meanz);
+    }
+
+    *aver_z = meanz / n;
+
+    G_debug(1, "In ortho aver_z = %f", *aver_z);
+
+    return 0;
+}

imagery/i.ortho.photo/i.photo.rectify/bilinear.c

@@ -0,0 +1,58 @@
+/*
+ * 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 "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;
+    double 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) {
+	G_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 (G_is_d_null_value(cellp)) {
+		G_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 = G_interp_bilinear(t, u, c[0][0], c[0][1], c[1][0], c[1][1]);
+
+    G_set_raster_value_d(obufptr, result, cell_type);
+}

imagery/i.ortho.photo/i.photo.rectify/bilinear_f.c

@@ -0,0 +1,48 @@
+/*
+ * 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 <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) {
+        G_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 (G_is_d_null_value(cellp)) {
+        G_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 (G_is_d_null_value(obufptr))
+        G_set_raster_value_d(obufptr, cell, cell_type);
+}

imagery/i.ortho.photo/i.photo.rectify/cp.c

@@ -0,0 +1,66 @@
+#include <stdlib.h>
+#include <string.h>
+#include "global.h"
+
+int get_conz_points(void)
+{
+    char msg[200];
+    /* struct Ortho_Control_Points cpz; */
+
+    if (!I_get_con_points(group.name, &group.control_points))
+	exit(0);
+
+    sprintf(msg, _("Control Z Point file for group [%s] in [%s] \n \n"),
+	    group.name, G_mapset());
+
+    G_verbose_message(_("Computing equations..."));
+
+    Compute_ortho_equation();
+
+    switch (group.con_equation_stat) {
+    case -1:
+	strcat(msg, _("Poorly placed Control Points!\n"));
+	strcat(msg, _("Can not generate the transformation equation.\n"));
+	strcat(msg, _("Run OPTION 7 of i.ortho.photo again!\n"));
+	break;
+    case 0:
+	strcat(msg, _("No active Control Points!\n"));
+	strcat(msg, _("Can not generate the transformation equation.\n"));
+	strcat(msg, _("Run OPTION 7 of i.ortho.photo!\n"));
+	break;
+    default:
+	return 1;
+    }
+    G_fatal_error(msg);
+}
+
+int get_ref_points(void)
+{
+    char msg[200];
+    /* struct Ortho_Photo_Points cp; */
+
+    if (!I_get_ref_points(group.name, &group.photo_points))
+	exit(0);
+
+    sprintf(msg, _("Reference Point file for group [%s] in [%s] \n \n"),
+	    group.name, G_mapset());
+
+    Compute_ref_equation();
+    switch (group.ref_equation_stat) {
+    case -1:
+	strcat(msg, _("Poorly placed Reference Points!\n"));
+	strcat(msg, _("Can not generate the transformation equation.\n"));
+	strcat(msg, _("Run OPTION 5 of i.ortho.photo again!\n"));
+	break;
+
+    case 0:
+	strcat(msg, _("No active Reference Points!\n"));
+	strcat(msg, _("Can not generate the transformation equation.\n"));
+	strcat(msg, _("Run OPTION 5 of i.ortho.photo!\n"));
+	break;
+    default:
+	return 1;
+    }
+    G_fatal_error(msg);
+    /* exit(1);   shouldn't get here */
+}

imagery/i.ortho.photo/i.photo.rectify/cubic.c

@@ -0,0 +1,66 @@
+/*
+ * 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 <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;
+    double 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) {
+	G_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 (G_is_d_null_value(cellp)) {
+		G_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] = G_interp_cubic(t, cell[i][0], cell[i][1], cell[i][2], cell[i][3]);
+    }
+
+    result = G_interp_cubic(u, val[0], val[1], val[2], val[3]);
+
+    G_set_raster_value_d(obufptr, result, cell_type);
+}

imagery/i.ortho.photo/i.photo.rectify/cubic_f.c

@@ -0,0 +1,53 @@
+/*
+ * 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 <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) {
+        G_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 (G_is_d_null_value(cellp)) {
+        G_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 (G_is_d_null_value(obufptr)) {
+        p_bilinear(ibuffer, obufptr, cell_type, row_idx, col_idx, cellhd);
+        /* fallback to nearest if bilinear is null */
+	if (G_is_d_null_value(obufptr))
+	    G_set_raster_value_d(obufptr, cell, cell_type);
+    }
+}

imagery/i.ortho.photo/i.photo.rectify/defs.h

@@ -0,0 +1,35 @@
+/* cache for raster data, code taken from r.proj */
+
+#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];   /* FCELL sufficient ? */
+
+struct cache
+{
+    int fd;
+    int stride;
+    int nblocks;
+    block **grid;
+    block *blocks;
+    int *refs;
+};
+
+typedef void (*func) (struct cache *, void *, int, double *, double *, struct Cell_head *);
+
+#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))])
+
+struct menu
+{
+    func method;		/* routine to interpolate new value      */
+    char *name;			/* method name                           */
+    char *text;			/* menu display - full description       */
+};
+

imagery/i.ortho.photo/i.photo.rectify/description.html

@@ -0,0 +1,110 @@
+<H2>DESCRIPTION</H2>
+
+<EM>i.photo.rectify</EM> rectifies an image by using the image to photo
+coordinate transformation matrix created by <A HREF="i.photo.2image.html">i.photo.2image</A>
+and the rectification parameters created by <A HREF="i.photo.2target.html">i.photo.2target</A>.
+Rectification is the process by which the geometry of an image is made
+planimetric.  This is accomplished by mapping an image from one coordinate
+system to another. In <EM>i.photo.rectify</EM> the parameters computed by
+<A HREF="i.photo.2image.html">i.photo.2image</A> and
+<A HREF="i.photo.2target.html">i.photo.2target</A> are used in equations to
+convert x,y image coordinates to standard map coordinates for each pixel in
+the image.  The result is an image with a standard map coordinate system,
+compensated for relief distortions and photographic tilt. Upon completion of
+the program the rectified image is deposited in a previously targeted GRASS
+LOCATION.
+<P>
+Images can be resampled with various different interpolation methods: 
+nearest neighbor assignment, bilinear and bicubic interpolation. The 
+bilinear and bicubic interpolation methods are also available with a
+fallback option. These methods "fall back" to simpler interpolation 
+methods along NULL borders. That is, from bicubic to bilinear to nearest.
+<P>
+The process may take an hour or more depending on the size of the image,
+the speed of the computer, the number files, and the size and resolution
+of the selected window.
+<P>
+The rectified image will be located in the target LOCATION when the program
+is completed. The original unrectified files are not modified or removed.
+<P>
+The optional <em>angle</em> output holds the camera angle in degrees to
+the local surface, considering local slope and aspect. A value of 90 
+degrees indicates that the camera angle was orthogonal to the local 
+surface, a value of 0 degrees indicates that the camera angle was 
+parallel to the local surface and negative values indicate that the 
+surface was invisible to the camera. As a rule of thumb, values below 30
+degrees indicate problem areas where the orthorectified output will 
+appear blurred. Because terrain shadowing effects are not considered, 
+areas with high camera angles may also appear blurred if they are located
+(viewed from the camera position) behind mountain ridges or peaks.
+<P>
+<EM>i.photo.rectify</EM> can be run directly, specifying options in the 
+command line or the GUI, or it can be invoked as OPTION 8 through 
+<A HREF="i.ortho.photo.html">i.ortho.photo</A>. If invoked though 
+<A HREF="i.ortho.photo.html">i.ortho.photo</A>, an interactive terminal 
+is used to determine the options.
+
+<H4>Interactive mode</H4>
+<P> You are first asked if all images within the imagery group should 
+be rectified. If this option is not chosen, you are asked to specify for 
+each image within the imagery group whether it should be rectified or not.
+<P>
+More than one file may be rectified at a time. Each file
+should have a unique output file name. The next prompt asks you for an 
+extension to be appended to the rectified images.
+<P>
+The next prompt will ask you whether a camera angle map should be 
+produced and if yes, what should be its name.
+<P>
+After that you are asked if overwriting existing maps in the target 
+location and mapset should be allowed.
+<P>
+The next prompt asks you to select one of two windows:
+<P>
+<PRE>
+      Please select one of the following options
+      1.   Use the current window in the target location
+      2.   Determine the smallest window which covers the image
+      &gt;
+</PRE>
+<P>
+If you choose option 2, you can also specify a desired target resolution.
+<P>
+<EM>i.photo.rectify</EM> will only rectify that portion of the
+image that occurs within the chosen window.  Only that portion will be
+relocated in the target database. It is therefore important to check the
+current window in the target LOCATION if choice number one is selected.
+<P>
+Next you are asked to select an interpolation method.
+<PRE>
+      Please select one of the following interpolation methods
+      1. nearest neighbor
+      2. bilinear
+      3. bicubic
+      4. bilinear with fallback
+      5. bicubic with fallback
+      &gt;
+</PRE>
+<P>
+The last prompt will ask you about the amount of memory to be used by
+<EM>i.photo.rectify</EM>.
+
+<H2>SEE ALSO</H2>
+
+<EM>
+<A HREF="i.ortho.photo.html">i.ortho.photo</A><br>
+<A HREF="i.photo.camera.html">i.photo.camera</A><br>
+<A HREF="i.photo.2image.html">i.photo.2image</A><br>
+<A HREF="i.photo.2target.html">i.photo.2target</A><br>
+<A HREF="i.photo.init.html">i.photo.init</A><br>
+<A HREF="i.rectify.html">i.rectify</A><br>
+</EM>
+
+
+<H2>AUTHOR</H2>
+Mike Baba,  DBA Systems, Inc.<br>
+Updated rectification and elevation map to FP 1/2002 Markus Neteler<br>
+Bugfixes and enhancements 12/2010 Markus Metz
+
+<p>
+<i>Last changed: $Date$</i>

imagery/i.ortho.photo/i.photo.rectify/env.c

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

imagery/i.ortho.photo/i.photo.rectify/equ.c

@@ -0,0 +1,60 @@
+#include "global.h"
+
+/* compute target to photo equation */
+int Compute_ortho_equation(void)
+{
+
+    /* struct Ortho_Control_Points  temp_points; */
+    double e0, e1, e2, n0, n1, n2, z1, z2;
+    int status, i;
+    struct Ortho_Control_Points temp_points;
+
+    /* 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);
+	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);
+
+    return 0;
+}
+
+/* compute photo to image equation */
+int Compute_ref_equation(void)
+{
+    group.ref_equation_stat = I_compute_ref_equations(&group.photo_points,
+						      group.E12, group.N12,
+						      group.E21, group.N21);
+
+    return 0;
+}

imagery/i.ortho.photo/i.photo.rectify/exec.c

@@ -0,0 +1,121 @@
+/*
+   exec.c --
+
+   Loop through all files to be rectified and do the retification.
+   Handles things like support files.
+ */
+
+#include <string.h>
+#include <unistd.h>
+#include <time.h>
+#include <sys/types.h>
+#include <sys/stat.h>
+#include <fcntl.h>
+#include "global.h"
+
+int exec_rectify(char *extension, char *interp_method, char *angle_map)
+{
+    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;
+    double aver_z;
+    int elevfd;
+    struct cache *ebuffer;
+
+    G_debug(1, "Open elevation raster: ");
+
+    /* open elevation raster */
+    select_target_env();
+    G_set_window(&target_window);
+    G_debug(1, "target window: rs=%d cs=%d n=%f s=%f w=%f e=%f\n",
+	    target_window.rows, target_window.cols, target_window.north,
+	    target_window.south, target_window.west, target_window.east);
+
+    elevfd = G_open_cell_old(elev_name, elev_mapset);
+    if (elevfd < 0) {
+	G_fatal_error(_("Could not open elevation raster"));
+	return 1;
+    }
+    ebuffer = readcell(elevfd, seg_mb_elev, 1);
+    select_target_env();
+    G_close_cell(elevfd);
+
+    /* get an average elevation of the control points */
+    /* this is used only if target cells have no elevation */
+    get_aver_elev(&group.control_points, &aver_z);
+
+    G_message("-----------------------------------------------");
+
+    /* rectify each file */
+    for (n = 0; n < group.group_ref.nfiles; n++) {
+	if (!ref_list[n])
+	    continue;
+
+	name = group.group_ref.file[n].name;
+	mapset = group.group_ref.file[n].mapset;
+	result =
+	    G_malloc(strlen(group.group_ref.file[n].name) + strlen(extension) + 1);
+	strcpy(result, group.group_ref.file[n].name);
+	strcat(result, extension);
+
+	G_debug(2, "ORTHO RECTIFYING:");
+	G_debug(2, "NAME %s", name);
+	G_debug(2, "MAPSET %s", mapset);
+	G_debug(2, "RESULT %s", result);
+	G_debug(2, "select_current_env...");
+
+	select_current_env();
+
+	cats_ok = G_read_cats(name, mapset, &cats) >= 0;
+	colr_ok = G_read_colors(name, mapset, &colr) > 0;
+
+	/* Initialze History */
+	if (G_read_history(name, mapset, &hist) < 0)
+	    G_short_history(result, type, &hist);
+	G_debug(2, "reading was fine...");
+
+	time(&start_time);
+
+	G_debug(2, "Starting the rectification...");
+
+	if (rectify(name, mapset, ebuffer, aver_z, result, interp_method)) {
+	    G_debug(2, "Done. Writing results...");
+	    select_target_env();
+	    if (cats_ok) {
+		G_write_cats(result, &cats);
+		G_free_cats(&cats);
+	    }
+	    if (colr_ok) {
+		G_write_colors(result, G_mapset(), &colr);
+		G_free_colors(&colr);
+	    }
+	    /* Write out History */
+	    G_command_history(&hist);
+	    G_write_history(result, &hist);
+
+	    select_current_env();
+	    time(&rectify_time);
+	    report(rectify_time - start_time, 1);
+	}
+	else
+	    report((long)0, 0);
+
+	G_free(result);
+    }
+    
+    close(ebuffer->fd);
+    release_cache(ebuffer);
+
+    if (angle_map) {
+	camera_angle(angle_map);
+    }
+    
+    return 0;
+}

imagery/i.ortho.photo/i.photo.rectify/get_wind.c

@@ -0,0 +1,216 @@
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include "global.h"
+
+int get_ref_window(struct Cell_head *cellhd)
+{
+    int i;
+    int count;
+    struct Cell_head win;
+
+    /* from all the files in the group, get max extends and min resolutions */
+    count = 0;
+    for (i = 0; i < group.group_ref.nfiles; i++) {
+	if (!ref_list[i])
+	    continue;
+
+	if (count++ == 0) {
+	    G_get_cellhd(group.group_ref.file[i].name,
+			 group.group_ref.file[i].mapset, cellhd);
+	}
+	else {
+	    G_get_cellhd(group.group_ref.file[i].name,
+			 group.group_ref.file[i].mapset, &win);
+	    /* max extends */
+	    if (cellhd->north < win.north)
+		cellhd->north = win.north;
+	    if (cellhd->south > win.south)
+		cellhd->south = win.south;
+	    if (cellhd->west > win.west)
+		cellhd->west = win.west;
+	    if (cellhd->east < win.east)
+		cellhd->east = win.east;
+	    /* min resolution */
+	    if (cellhd->ns_res > win.ns_res)
+		cellhd->ns_res = win.ns_res;
+	    if (cellhd->ew_res > win.ew_res)
+		cellhd->ew_res = win.ew_res;
+	}
+    }
+
+    /* if the north-south is not multiple of the resolution,
+     *    round the south downward
+     */
+    cellhd->rows = (cellhd->north - cellhd->south) /cellhd->ns_res + 0.5;
+    cellhd->south = cellhd->north - cellhd->rows * cellhd->ns_res;
+
+    /* do the same for the west */
+    cellhd->cols = (cellhd->east - cellhd->west) / cellhd->ew_res + 0.5;
+    cellhd->west = cellhd->east - cellhd->cols * cellhd->ew_res;
+
+    return 1;
+}
+
+int georef_window(struct Cell_head *w1, struct Cell_head *w2, double res)
+{
+    double n, e, z1, ad;
+    double n0, e0;
+    double aver_z;
+    struct _corner {
+        double n, e;
+    } nw, ne, se, sw;
+
+    /* get an average elevation from the active control points */
+    /* for mountainous areas this is a very rough approximation
+     * which would become more accurate only if actual elevation 
+     * values are used */
+    get_aver_elev(&group.control_points, &aver_z);
+    G_debug(1, "Aver elev = %f", aver_z);
+
+    /* compute ortho ref of all corners */
+
+    I_georef(w1->west, w1->north, &e0, &n0, group.E12, group.N12);
+    I_inverse_ortho_ref(e0, n0, aver_z, &e, &n, &z1, &group.camera_ref,
+			group.XC, group.YC, group.ZC, group.MI);
+
+    G_debug(1, "NORTH WEST CORNER");
+    G_debug(1, "group.E12 = %f %f %f,", group.E12[0], group.E12[1],
+	    group.E12[2]);
+    G_debug(1, "group.N12 = %f %f %f,", group.N12[0], group.N12[1],
+	    group.N12[2]);
+    G_debug(1, "image  x = %f y = %f, photo x = %f y = %f", w1->west,
+	    w1->north, e0, n0);
+    G_debug(1, "target x = %f y = %f", e, n);
+
+    w2->north = w2->south = n;
+    w2->west = w2->east = e;
+    nw.n = n;
+    nw.e = e;
+
+    I_georef(w1->east, w1->north, &e0, &n0, group.E12, group.N12);
+    I_inverse_ortho_ref(e0, n0, aver_z, &e, &n, &z1, &group.camera_ref,
+			group.XC, group.YC, group.ZC, group.MI);
+
+    G_debug(1, "NORTH EAST CORNER");
+    G_debug(1, "image  x = %f y = %f, photo x = %f y = %f", w1->east,
+	    w1->north, e0, n0);
+    G_debug(1, "target x = %f y = %f", e, n);
+
+    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;
+
+    I_georef(w1->west, w1->south, &e0, &n0, group.E12, group.N12);
+    I_inverse_ortho_ref(e0, n0, aver_z, &e, &n, &z1, &group.camera_ref,
+			group.XC, group.YC, group.ZC, group.MI);
+
+    G_debug(1, "SOUTH WEST CORNER");
+    G_debug(1, "image  x = %f y = %f, photo x = %f y = %f", w1->west,
+	    w1->south, e0, n0);
+    G_debug(1, "target x = %f y = %f", e, n);
+
+    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;
+
+    I_georef(w1->east, w1->south, &e0, &n0, group.E12, group.N12);
+    I_inverse_ortho_ref(e0, n0, aver_z, &e, &n, &z1, &group.camera_ref,
+			group.XC, group.YC, group.ZC, group.MI);
+
+    G_debug(1, "SOUTH EAST CORNER");
+    G_debug(1, "image  x = %f y = %f, photo x = %f y = %f", w1->east,
+	    w1->south, e0, n0);
+    G_debug(1, "target x = %f y = %f", e, n);
+
+    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_debug(1, "FINAL");
+    G_debug(1, "east = %f \n west = %f \n north = %f \n south = %f",
+	    w2->east, w2->west, w2->north, w2->south);
+    G_debug(1, "RESOLUTION");
+    G_debug(1, "EW = %f", w2->ew_res);
+    G_debug(1, "NS = %f", w2->ns_res);
+
+    return 0;
+}

imagery/i.ortho.photo/i.photo.rectify/global.h

@@ -0,0 +1,18 @@
+#include <grass/gis.h>
+#include <grass/imagery.h>
+#include <grass/ortholib.h>
+#include <grass/glocale.h>
+#include "orthophoto.h"
+#include "defs.h"
+
+extern func interpolate;	/* interpolation routine */
+
+extern int seg_mb_img, seg_mb_elev;
+extern char *extension;
+extern int *ref_list;
+extern struct Ortho_Image_Group group;
+extern char *elev_name;
+extern char *elev_mapset;
+extern struct Cell_head target_window;
+
+#include "local_proto.h"

imagery/i.ortho.photo/i.photo.rectify/local_proto.h

@@ -0,0 +1,59 @@
+/* angle.c */
+int camera_angle(char *);
+
+/* aver_z.c */
+int get_aver_elev(struct Ortho_Control_Points *, double *);
+
+/* cp.c */
+int get_conz_points(void);
+int get_ref_points(void);
+
+/* elev_data.c */
+int elev_data(char *, int);
+
+/* env.c */
+int select_current_env(void);
+int select_target_env(void);
+int show_env(void);
+
+/* equ.c */
+int Compute_ortho_equation(void);
+int Compute_ref_equation(void);
+
+/* exec.c */
+int exec_rectify(char *, char *, char *);
+
+/* get_wind.c */
+int get_ref_window(struct Cell_head *);
+int georef_window(struct Cell_head *, struct Cell_head *, double);
+
+/* rectify.c */
+int rectify(char *, char *, struct cache *, double, char *, char *);
+
+/* readcell.c */
+struct cache *readcell(int, int, int);
+block *get_block(struct cache *, int);
+void release_cache(struct cache *);
+
+/* 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 *);

imagery/i.ortho.photo/i.photo.rectify/main.c

@@ -0,0 +1,412 @@
+
+/****************************************************************************
+ *
+ * MODULE:       i.photo.rectify
+ * AUTHOR(S):    Mike Baba,  DBA Systems, Inc. (original contributor)
+ *               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:      Rectifies an image by using the image to photo coordinate 
+ *               and photo to target transformation matrices
+ * COPYRIGHT:    (C) 1999-2010 by the GRASS Development Team
+ *
+ *               This program is free software under the GNU General Public
+ *               License (>=v2). Read the file COPYING that comes with GRASS
+ *               for details.
+ *
+ *****************************************************************************/
+
+#include <stdlib.h>
+#include <string.h>
+#include "global.h"
+
+int seg_mb_img, seg_mb_elev;
+
+int *ref_list;
+struct Ortho_Image_Group group;
+
+char *elev_name;
+char *elev_mapset;
+
+func interpolate;
+
+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_bilinear_f, "bilinear_f", "bilinear with fallback"},
+    {p_cubic_f, "cubic_f", "cubic convolution with fallback"},
+    {NULL, NULL, NULL}
+};
+
+static char *make_ipol_list(void);
+
+int main(int argc, char *argv[])
+{
+    char extension[INAME_LEN];
+    char *ipolname;		/* name of interpolation method */
+    int method;
+    char *seg_mb;
+    int i, m, k = 0;
+    int got_file = 0, target_overwrite = 0;
+    char *overstr;
+
+    char *camera;
+    int n, nfiles;
+    char tl[100];
+    char math_exp[100];
+    char units[100];
+    char nd[100];
+
+    struct Cell_head cellhd, elevhd;
+
+    struct Option *grp,         /* imagery group */
+     *ifile,			/* input files */
+     *ext,			/* extension */
+     *tres,			/* target resolution */
+     *mem,			/* amount of memory for cache */
+     *interpol,			/* interpolation method:
+				   nearest neighbor, bilinear, cubic */
+     *angle;			/* camera angle relative to ground surface */
+
+    struct Flag *c, *a;
+    struct GModule *module;
+
+    G_gisinit(argv[0]);
+
+    module = G_define_module();
+    module->keywords = _("imagery, orthorectify");
+    module->description =
+	_("Orthorectifies an image by using the image to photo coordinate transformation matrix.");
+
+    grp = G_define_standard_option(G_OPT_I_GROUP);
+
+    ifile = G_define_standard_option(G_OPT_R_INPUTS);
+    ifile->required = NO;
+
+    ext = G_define_option();
+    ext->key = "extension";
+    ext->type = TYPE_STRING;
+    ext->required = YES;
+    ext->multiple = NO;
+    ext->description = _("Output raster map(s) suffix");
+
+    tres = G_define_option();
+    tres->key = "resolution";
+    tres->type = TYPE_DOUBLE;
+    tres->required = NO;
+    tres->description = _("Target resolution (ignored if -c flag used)");
+
+    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");
+
+    ipolname = make_ipol_list();
+
+    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");
+    
+    angle = G_define_standard_option(G_OPT_R_OUTPUT);
+    angle->key = "angle";
+    angle->required = NO;
+    angle->description = _("Raster map with camera angle relative to ground surface");
+
+    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.name, grp->answer);
+    strcpy(extension, ext->answer);
+
+    seg_mb = NULL;
+    if (mem->answer) {
+	if (atoi(mem->answer) > 0)
+	    seg_mb = mem->answer;
+    }
+
+    if (!ifile->answers)
+	a->answer = 1;		/* force all */
+
+    /* Find out how many files on command line */
+    if (!a->answer) {
+	for (k = 0; ifile->answers[k]; k++);
+    }
+
+    camera = (char *)G_malloc(GNAME_MAX * sizeof(char));
+    elev_name = (char *)G_malloc(GNAME_MAX * sizeof(char));
+    elev_mapset = (char *)G_malloc(GMAPSET_MAX * sizeof(char));
+
+    /* find group */
+    if (!I_find_group(group.name))
+	G_fatal_error(_("Group <%s> not found"), group.name);
+
+    /* get the group ref */
+    if (!I_get_group_ref(group.name, (struct Ref *)&group.group_ref))
+	G_fatal_error(_("Could not read REF file for group <%s>"), group.name);
+    nfiles = group.group_ref.nfiles;
+    if (nfiles <= 0) {
+	G_important_message(_("Group <%s> contains no raster maps; run i.group"),
+			    grp->answer);
+	exit(EXIT_SUCCESS);
+    }
+
+    ref_list = (int *)G_malloc(nfiles * sizeof(int));
+
+    if (a->answer) {
+	for (n = 0; n < group.group_ref.nfiles; n++) {
+	    ref_list[n] = 1;
+	}
+    }
+    else {
+	char xname[GNAME_MAX], xmapset[GMAPSET_MAX], *name, *mapset;
+
+	for (n = 0; n < group.group_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 < group.group_ref.nfiles; n++) {
+		if (mapset) {
+		    if (strcmp(name, group.group_ref.file[n].name) == 0 &&
+		        strcmp(mapset, group.group_ref.file[n].mapset) == 0) {
+			got_file = 1;
+			ref_list[n] = 1;
+			break;
+		    }
+		}
+		else {
+		    if (strcmp(name, group.group_ref.file[n].name) == 0) {
+			got_file = 1;
+			ref_list[n] = 1;
+			break;
+		    }
+		}
+	    }
+	    if (got_file == 0)
+		err_exit(ifile->answers[m], group.name);
+	}
+    }
+
+    /** look for camera info for this block **/
+    if (!I_get_group_camera(group.name, camera))
+	G_fatal_error(_("No camera reference file selected for group <%s>"),
+		      group.name);
+
+    if (!I_get_cam_info(camera, &group.camera_ref))
+	G_fatal_error(_("Bad format in camera file for group <%s>"),
+		      group.name);
+
+    /* get initial camera exposure station, if any */
+    if (I_find_initial(group.name)) {
+	if (!I_get_init_info(group.name, &group.camera_exp))
+	    G_warning(_("Bad format in initial exposure station file for group <%s>"),
+		      group.name);
+    }
+
+    /* read the reference points for the group, compute image-to-photo trans. */
+    get_ref_points();
+
+    /* read the control points for the group, convert to photo coords. */
+    get_conz_points();
+
+    /* get the target */
+    get_target(group.name);
+    
+    /* Check the GRASS_OVERWRITE environment variable */
+    if ((overstr = getenv("GRASS_OVERWRITE")))  /* OK ? */
+	target_overwrite = atoi(overstr);
+
+    if (!target_overwrite) {
+	/* check if output exists in target location/mapset */
+	char result[GNAME_MAX];
+	
+	select_target_env();
+	for (i = 0; i < group.group_ref.nfiles; i++) {
+	    if (!ref_list[i])
+		continue;
+
+	    strcpy(result, group.group_ref.file[i].name);
+	    strcat(result, extension);
+	    
+	    if (G_legal_filename(result) < 0)
+		G_fatal_error(_("Extension <%s> is illegal"), extension);
+		
+	    if (G_find_cell(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 (angle->answer) {
+	    if (G_find_cell(angle->answer, 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>", angle->answer);
+		G_fatal_error(_("Orthorectification cancelled."));
+	    }
+	}
+	
+	select_current_env();
+    }
+    else
+	G_debug(1, "Overwriting OK");
+
+    /* do not use current region in target location */
+    if (!c->answer) {
+	double res = -1;
+	
+	if (tres->answer) {
+	    if (!((res = atof(tres->answer)) > 0))
+		G_warning(_("Target resolution must be > 0, ignored"));
+	}
+	/* get reference window from imagery group */
+	get_ref_window(&cellhd);
+	georef_window(&cellhd, &target_window, res);
+    }
+
+    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);
+
+    G_debug(1, "Looking for elevation file in group: <%s>", group.name);
+
+    /* get the block elevation layer raster map in target location */
+    if (!I_get_group_elev(group.name, elev_name, elev_mapset, tl,
+			 math_exp, units, nd))
+	G_fatal_error(_("No target elevation model selected for group <%s>"),
+		      group.name);
+
+    G_debug(1, "Block elevation: <%s> in <%s>", elev_name, elev_mapset);
+
+    /* get the elevation layer header in target location */
+    select_target_env();
+    G_get_cellhd(elev_name, elev_mapset, &elevhd);
+    select_current_env();
+    
+    /* determine memory for elevation and imagery */
+    seg_mb_img = seg_mb_elev = -1;
+    if (seg_mb) {
+	int max_rows, max_cols;
+	int nx, ny;
+	double max_mb_img, max_mb_elev;
+	int seg_mb_total;
+
+	max_rows = max_cols = 0;
+	for (i = 0; i < group.group_ref.nfiles; i++) {
+	    if (ref_list[i]) {
+		G_get_cellhd(group.group_ref.file[i].name,
+			     group.group_ref.file[i].mapset, &cellhd);
+		if (max_rows < cellhd.rows)
+		    max_rows = cellhd.rows;
+		if (max_cols < cellhd.cols)
+		    max_cols = cellhd.cols;
+	    }
+	}
+
+	ny = (max_rows + BDIM - 1) / BDIM;
+	nx = (max_cols + BDIM - 1) / BDIM;
+
+	max_mb_img = ((double)nx * ny * sizeof(block)) / (1<<20);
+
+	ny = (target_window.rows + BDIM - 1) / BDIM;
+	nx = (target_window.cols + BDIM - 1) / BDIM;
+
+	max_mb_elev = ((double)nx * ny * sizeof(block)) / (1<<20);
+
+	if ((seg_mb_total = atoi(seg_mb)) > 0) {
+	    seg_mb_elev = max_mb_elev * (seg_mb_total / (max_mb_img + max_mb_elev)) + 0.5;
+	    seg_mb_img = max_mb_img * (seg_mb_total / (max_mb_img + max_mb_elev)) + 0.5;
+	}
+    }
+
+    /* go do it */
+    exec_rectify(extension, interpol->answer, angle->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 < group.group_ref.nfiles; n++)
+	G_message("%s@%s", group.group_ref.file[n].name, group.group_ref.file[n].mapset);
+
+    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;
+}

imagery/i.ortho.photo/i.photo.rectify/nearest.c

@@ -0,0 +1,39 @@
+/*
+ *      nearest.c - returns the nearest neighbor to a given
+ *                  x,y position
+ */
+
+#include <math.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) {
+	G_set_null_value(obufptr, 1, cell_type);
+	return;
+    }
+
+    cellp = CPTR(ibuffer, row, col);
+
+    if (G_is_d_null_value(cellp)) {
+	G_set_null_value(obufptr, 1, cell_type);
+	return;
+    }
+
+    G_set_raster_value_d(obufptr, *cellp, cell_type);
+}

imagery/i.ortho.photo/i.photo.rectify/readcell.c

@@ -0,0 +1,153 @@
+/*
+ * 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 "global.h"
+
+struct cache *readcell(int fdi, int size, int target_env)
+{
+    DCELL *tmpbuf;
+    struct cache *c;
+    int nrows;
+    int ncols;
+    int row;
+    char *filename;
+    int nx, ny;
+    int nblocks;
+    int i;
+
+    if (target_env)
+	select_target_env();
+    else
+	select_current_env();
+
+    nrows = G_window_rows();
+    ncols = G_window_cols();
+
+    /* Temporary file must be created in the same location/mapset 
+     * where the module was called */
+    if (target_env)
+	select_current_env();
+
+    ny = (nrows + BDIM - 1) / BDIM;
+    nx = (ncols + BDIM - 1) / BDIM;
+
+    if (size > 0)
+	nblocks = 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) {
+	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_debug(1, "%d of %d blocks in memory", nblocks, nx * ny);
+
+    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));
+
+    if (target_env)
+	select_target_env();
+    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;
+
+	    G_get_d_raster_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;
+	}
+
+    if (target_env)
+	select_current_env();
+
+    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 reading segment file"));
+
+    return p;
+}
+
+void release_cache(struct cache *c)
+{
+    G_free(c->refs);
+    G_free(c->blocks);
+    G_free(c->grid);
+    
+    G_free(c);
+    c = NULL;
+}

imagery/i.ortho.photo/i.photo.rectify/rectify.c

@@ -0,0 +1,144 @@
+/* rectification code */
+
+/* 1/2002: updated to GRASS 5 write routines and 
+   CELL/FP elevation - Markus Neteler
+ */
+
+#include <stdlib.h>
+#include <string.h>
+#include <unistd.h>
+#include "global.h"
+
+int rectify(char *name, char *mapset, struct cache *ebuffer,
+            double aver_z, char *result, char *interp_method)
+{
+    struct Cell_head cellhd;
+    int ncols, nrows;
+    int row, col;
+    double row_idx, col_idx;
+    int infd, outfd;
+    RASTER_MAP_TYPE map_type;
+    int cell_size;
+    void *trast, *tptr;
+    double n1, e1, z1;
+    double nx, ex, nx1, ex1, zx1;
+    struct cache *ibuffer;
+
+    select_current_env();
+    if (G_get_cellhd(name, mapset, &cellhd) < 0)
+	return 0;
+
+    /* open the file to be rectified
+     * set window to cellhd first to be able to read file exactly
+     */
+    G_set_window(&cellhd);
+    infd = G_open_cell_old(name, mapset);
+    if (infd < 0) {
+	return 0;
+    }
+    map_type = G_get_raster_map_type(infd);
+    cell_size = G_raster_size(map_type);
+
+    ibuffer = readcell(infd, seg_mb_img, 0);
+
+    G_close_cell(infd);		/* (pmx) 17 april 2000 */
+
+    G_message(_("Rectify <%s@%s> (location <%s>)"),
+	      name, mapset, G_location());
+    select_target_env();
+    G_set_window(&target_window);
+    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 = G_raster_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 = G_open_raster_new(result, map_type);
+    trast = G_allocate_raster_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);
+
+	G_set_null_value(trast, ncols, map_type);
+	tptr = trast;
+	for (col = 0; col < ncols; col++) {
+	    DCELL *zp = CPTR(ebuffer, row, col);
+
+	    e1 = target_window.west + (col + 0.5) * target_window.ew_res;
+	    
+	    /* if target cell has no elevation, set to aver_z */
+	    if (G_is_d_null_value(zp)) {
+		G_warning(_("No elevation available at row = %d, col = %d"), row, col);
+		z1 = aver_z;
+	    }
+	    else
+		z1 = *zp;
+
+	    /* target coordinates e1, n1 to photo coordinates ex1, nx1 */
+	    I_ortho_ref(e1, n1, z1, &ex1, &nx1, &zx1, &group.camera_ref,
+			group.XC, group.YC, group.ZC, group.M);
+
+	    G_debug(5, "\t\tAfter ortho ref (photo cords): ex = %f \t nx =  %f",
+		    ex1, nx1);
+
+	    /* photo coordinates ex1, nx1 to image coordinates ex, nx */
+	    I_georef(ex1, nx1, &ex, &nx, group.E21, group.N21);
+
+	    G_debug(5, "\t\tAfter geo ref: ex = %f \t nx =  %f", ex, nx);
+
+	    /* 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);
+	}
+	G_put_raster_row(outfd, trast, map_type);
+    }
+    G_percent(1, 1, 1);
+
+    G_close_cell(outfd);
+    G_free(trast);
+
+    close(ibuffer->fd);
+    release_cache(ibuffer);
+
+    if (G_get_cellhd(result, G_mapset(), &cellhd) < 0)
+	return 0;
+
+    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;
+}

imagery/i.ortho.photo/i.photo.rectify/report.c

@@ -0,0 +1,33 @@
+#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;
+}

imagery/i.ortho.photo/i.photo.rectify/target.c

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