sync i.rectify and i.ortho.rectify

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

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

@@ -4,8 +4,8 @@ PGM = i.ortho.rectify
 
 EXTRA_CFLAGS = -I../lib
 
-LIBES     = $(IMAGERYLIB) $(RASTERLIB) $(GISLIB) $(IORTHOLIB) $(GMATHLIB)
-DEPENDENCIES= $(IMAGERYDEP) $(RASTERDEP) $(IORTHODEP) $(GISDEP) $(GMATHDEP)
+LIBES = $(IMAGERYLIB) $(RASTERLIB) $(GISLIB) $(IORTHOLIB) $(GMATHLIB)
+DEPENDENCIES = $(IMAGERYDEP) $(RASTERDEP) $(IORTHODEP) $(GISDEP) $(GMATHDEP)
 
 include $(MODULE_TOPDIR)/include/Make/Module.make
 

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

@@ -11,3 +11,4 @@
    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 
 
+keep in sync with i.rectify

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

@@ -13,12 +13,12 @@
 #include <math.h>
 #include "global.h"
 
-int camera_angle(char *name)
+int camera_angle(struct Ortho_Image_Group * group, char *name)
 {
     int row, col, nrows, ncols;
-    double XC = group.XC;
-    double YC = group.YC;
-    double ZC = group.ZC;
+    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; */

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

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

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

@@ -1,5 +1,12 @@
 /* cache for raster data, code taken from r.proj */
 
+/* These next defines determine the size of the sub-window that will
+ * be held in memory.  Larger values will require
+ * more memory (but less i/o) If you increase these values, keep  in
+ * mind that although you think the i/o will decrease, system paging
+ * (which goes on behind the scenes) may actual increase the i/o.
+ */
+
 #define L2BDIM 6
 #define BDIM (1<<(L2BDIM))
 #define L2BSIZE (2*(L2BDIM))

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

@@ -1,10 +1,8 @@
 #include "global.h"
 
 /* compute target to photo equation */
-int Compute_ortho_equation(void)
+int Compute_ortho_equation(struct Ortho_Image_Group *group)
 {
-
-    /* struct Ortho_Control_Points  temp_points; */
     double e0, e1, e2, n0, n1, n2, z1, z2;
     int status, i;
     struct Ortho_Control_Points temp_points;
@@ -20,41 +18,41 @@ int Compute_ortho_equation(void)
     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];
+    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_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);
+    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)
+int Compute_ref_equation(struct Ortho_Image_Group *group)
 {
-    group.ref_equation_stat = I_compute_ref_equations(&group.photo_points,
-						      group.E12, group.N12,
-						      group.E21, group.N21);
+    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.ortho.rectify/exec.c

@@ -1,19 +1,23 @@
 /*
-   exec.c --
-
-   Loop through all files to be rectified and do the retification.
-   Handles things like support files.
+ * exec.c --
+ *
+ * Loop through all files to be rectified and do the retification.
+ * Handles things like support files.
  */
 
+#include <stdlib.h>
+#include <stdio.h>
 #include <string.h>
-#include <unistd.h>
-#include <time.h>
 #include <sys/types.h>
 #include <sys/stat.h>
 #include <fcntl.h>
+#include <time.h>
+#include <unistd.h>
+#include <math.h>
 #include "global.h"
 
-int exec_rectify(char *extension, char *interp_method, char *angle_map)
+int exec_rectify(struct Ortho_Image_Group *group, int *ref_list,
+                 char *extension, char *interp_method, char *angle_map)
 {
     char *name;
     char *mapset;
@@ -49,20 +53,20 @@ int exec_rectify(char *extension, char *interp_method, char *angle_map)
 
     /* 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);
+    get_aver_elev(&group->control_points, &aver_z);
 
     G_message("-----------------------------------------------");
 
     /* rectify each file */
-    for (n = 0; n < group.group_ref.nfiles; n++) {
+    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;
+	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);
+	    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:");
@@ -85,7 +89,7 @@ int exec_rectify(char *extension, char *interp_method, char *angle_map)
 
 	G_debug(2, "Starting the rectification...");
 
-	if (rectify(name, mapset, ebuffer, aver_z, result, interp_method)) {
+	if (rectify(group, name, mapset, ebuffer, aver_z, result, interp_method)) {
 	    G_debug(2, "Done. Writing results...");
 	    select_target_env();
 	    if (cats_ok) {
@@ -114,7 +118,7 @@ int exec_rectify(char *extension, char *interp_method, char *angle_map)
     release_cache(ebuffer);
 
     if (angle_map) {
-	camera_angle(angle_map);
+	camera_angle(group, angle_map);
     }
     
     return 0;

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

@@ -3,7 +3,7 @@
 #include <math.h>
 #include "global.h"
 
-int get_ref_window(struct Cell_head *cellhd)
+int get_ref_window(struct Ref *ref, int *ref_list, struct Cell_head *cellhd)
 {
     int i;
     int count;
@@ -11,17 +11,17 @@ int get_ref_window(struct Cell_head *cellhd)
 
     /* from all the files in the group, get max extends and min resolutions */
     count = 0;
-    for (i = 0; i < group.group_ref.nfiles; i++) {
+    for (i = 0; i < ref->nfiles; i++) {
 	if (!ref_list[i])
 	    continue;
 
 	if (count++ == 0) {
-	    Rast_get_cellhd(group.group_ref.file[i].name,
-			 group.group_ref.file[i].mapset, cellhd);
+	    Rast_get_cellhd(ref->file[i].name,
+			    ref->file[i].mapset, cellhd);
 	}
 	else {
-	    Rast_get_cellhd(group.group_ref.file[i].name,
-			 group.group_ref.file[i].mapset, &win);
+	    Rast_get_cellhd(ref->file[i].name,
+			    ref->file[i].mapset, &win);
 	    /* max extends */
 	    if (cellhd->north < win.north)
 		cellhd->north = win.north;
@@ -42,7 +42,7 @@ int get_ref_window(struct Cell_head *cellhd)
     /* 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->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 */
@@ -52,7 +52,8 @@ int get_ref_window(struct Cell_head *cellhd)
     return 1;
 }
 
-int georef_window(struct Cell_head *w1, struct Cell_head *w2, double res)
+int georef_window(struct Ortho_Image_Group *group, struct Cell_head *w1,
+                  struct Cell_head *w2, double res)
 {
     double n, e, z1, ad;
     double n0, e0;
@@ -65,20 +66,20 @@ int georef_window(struct Cell_head *w1, struct Cell_head *w2, double res)
     /* 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);
+    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, 1);
-    I_inverse_ortho_ref(e0, n0, aver_z, &e, &n, &z1, &group.camera_ref,
-			group.XC, group.YC, group.ZC, group.MI);
+    I_georef(w1->west, w1->north, &e0, &n0, group->E12, group->N12, 1);
+    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, "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);
@@ -88,9 +89,9 @@ int georef_window(struct Cell_head *w1, struct Cell_head *w2, double res)
     nw.n = n;
     nw.e = e;
 
-    I_georef(w1->east, w1->north, &e0, &n0, group.E12, group.N12, 1);
-    I_inverse_ortho_ref(e0, n0, aver_z, &e, &n, &z1, &group.camera_ref,
-			group.XC, group.YC, group.ZC, group.MI);
+    I_georef(w1->east, w1->north, &e0, &n0, group->E12, group->N12, 1);
+    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,
@@ -108,9 +109,9 @@ int georef_window(struct Cell_head *w1, struct Cell_head *w2, double res)
     if (e < w2->west)
 	w2->west = e;
 
-    I_georef(w1->west, w1->south, &e0, &n0, group.E12, group.N12, 1);
-    I_inverse_ortho_ref(e0, n0, aver_z, &e, &n, &z1, &group.camera_ref,
-			group.XC, group.YC, group.ZC, group.MI);
+    I_georef(w1->west, w1->south, &e0, &n0, group->E12, group->N12, 1);
+    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,
@@ -128,9 +129,9 @@ int georef_window(struct Cell_head *w1, struct Cell_head *w2, double res)
     if (e < w2->west)
 	w2->west = e;
 
-    I_georef(w1->east, w1->south, &e0, &n0, group.E12, group.N12, 1);
-    I_inverse_ortho_ref(e0, n0, aver_z, &e, &n, &z1, &group.camera_ref,
-			group.XC, group.YC, group.ZC, group.MI);
+    I_georef(w1->east, w1->south, &e0, &n0, group->E12, group->N12, 1);
+    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,
@@ -192,7 +193,7 @@ int georef_window(struct Cell_head *w1, struct Cell_head *w2, double res)
 	}
     }
 
-    /* adjust extends */
+    /* adjust extents */
     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;

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

@@ -8,9 +8,6 @@
 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;

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

@@ -1,12 +1,12 @@
 /* angle.c */
-int camera_angle(char *);
+int camera_angle(struct Ortho_Image_Group *, 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);
+int get_conz_points(struct Ortho_Image_Group *);
+int get_ref_points(struct Ortho_Image_Group *);
 
 /* elev_data.c */
 int elev_data(char *, int);
@@ -17,18 +17,18 @@ int select_target_env(void);
 int show_env(void);
 
 /* equ.c */
-int Compute_ortho_equation(void);
-int Compute_ref_equation(void);
+int Compute_ortho_equation(struct Ortho_Image_Group *);
+int Compute_ref_equation(struct Ortho_Image_Group *);
 
 /* exec.c */
-int exec_rectify(char *, char *, char *);
+int exec_rectify(struct Ortho_Image_Group *, int *, char *, char *, char *);
 
 /* get_wind.c */
-int get_ref_window(struct Cell_head *);
-int georef_window(struct Cell_head *, struct Cell_head *, double);
+int get_ref_window(struct Ref *, int *, struct Cell_head *);
+int georef_window(struct Ortho_Image_Group *, struct Cell_head *, struct Cell_head *, double);
 
 /* rectify.c */
-int rectify(char *, char *, struct cache *, double, char *, char *);
+int rectify(struct Ortho_Image_Group *, char *, char *, struct cache *, double, char *, char *);
 
 /* readcell.c */
 struct cache *readcell(int, int, int);
@@ -36,7 +36,7 @@ block *get_block(struct cache *, int);
 void release_cache(struct cache *);
 
 /* report.c */
-int report(int, int);
+int report(time_t, int);
 
 /* target.c */
 int get_target(char *);

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

@@ -25,9 +25,6 @@
 
 int seg_mb_img, seg_mb_elev;
 
-int *ref_list;
-struct Ortho_Image_Group group;
-
 char *elev_name;
 char *elev_mapset;
 
@@ -35,7 +32,7 @@ func interpolate;
 
 struct Cell_head target_window;
 
-void err_exit(char *, char *);
+void err_exit(struct Ref *, char *, char *);
 
 /* modify this table to add new methods */
 struct menu menu[] = {
@@ -61,8 +58,10 @@ int main(int argc, char *argv[])
     int got_file = 0, target_overwrite = 0;
     char *overstr;
 
+    struct Ortho_Image_Group group;
+    int *ref_list;
+    int n;
     char *camera;
-    int n, nfiles;
     char tl[100];
     char math_exp[100];
     char units[100];
@@ -75,9 +74,9 @@ int main(int argc, char *argv[])
      *ext,			/* extension */
      *tres,			/* target resolution */
      *mem,			/* amount of memory for cache */
-     *interpol,			/* interpolation method:
+     *angle,			/* camera angle relative to ground surface */
+     *interpol;			/* interpolation method:
 				   nearest neighbor, bilinear, cubic */
-     *angle;			/* camera angle relative to ground surface */
 
     struct Flag *c, *a;
     struct GModule *module;
@@ -176,20 +175,25 @@ int main(int argc, char *argv[])
     elev_mapset = (char *)G_malloc(GMAPSET_MAX * sizeof(char));
 
     /* find group */
-    if (!I_find_group(group.name))
+    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) {
+    /* determine the number of files in this group */
+    if (!I_get_group_ref(group.name, &group.group_ref)) {
+	G_warning(_("Location: %s"), G_location());
+	G_warning(_("Mapset: %s"), G_mapset());
+	G_fatal_error(_("Could not read REF file for group <%s>"),
+	              group.name);
+    }
+
+    if (group.group_ref.nfiles <= 0) {
 	G_important_message(_("Group <%s> contains no raster maps; run i.group"),
 			    grp->answer);
 	exit(EXIT_SUCCESS);
     }
 
-    ref_list = (int *)G_malloc(nfiles * sizeof(int));
+    ref_list = (int *)G_malloc(group.group_ref.nfiles * sizeof(int));
 
     if (a->answer) {
 	for (n = 0; n < group.group_ref.nfiles; n++) {
@@ -232,7 +236,7 @@ int main(int argc, char *argv[])
 		}
 	    }
 	    if (got_file == 0)
-		err_exit(ifile->answers[m], group.name);
+		err_exit(&group.group_ref, ifile->answers[m], group.name);
 	}
     }
 
@@ -253,14 +257,14 @@ int main(int argc, char *argv[])
     }
 
     /* read the reference points for the group, compute image-to-photo trans. */
-    get_ref_points();
+    get_ref_points(&group);
 
     /* read the control points for the group, convert to photo coords. */
-    get_conz_points();
+    get_conz_points(&group);
 
     /* get the target */
     get_target(group.name);
-    
+
     /* Check the GRASS_OVERWRITE environment variable */
     if ((overstr = getenv("GRASS_OVERWRITE")))  /* OK ? */
 	target_overwrite = atoi(overstr);
@@ -312,8 +316,8 @@ int main(int argc, char *argv[])
 		G_warning(_("Target resolution must be > 0, ignored"));
 	}
 	/* get reference window from imagery group */
-	get_ref_window(&cellhd);
-	georef_window(&cellhd, &target_window, res);
+	get_ref_window(&group.group_ref, ref_list, &cellhd);
+	georef_window(&group, &cellhd, &target_window, res);
     }
 
     G_verbose_message(_("Using region: N=%f S=%f, E=%f W=%f"), target_window.north,
@@ -371,7 +375,7 @@ int main(int argc, char *argv[])
     }
 
     /* go do it */
-    exec_rectify(extension, interpol->answer, angle->answer);
+    exec_rectify(&group, ref_list, extension, interpol->answer, angle->answer);
 
     G_done_msg(" ");
 
@@ -379,7 +383,7 @@ int main(int argc, char *argv[])
 }
 
 
-void err_exit(char *file, char *grp)
+void err_exit(struct Ref *ref, char *file, char *grp)
 {
     int n;
 
@@ -387,8 +391,8 @@ void err_exit(char *file, char *grp)
 	    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);
+    for (n = 0; n < ref->nfiles; n++)
+	G_message("%s@%s", ref->file[n].name, ref->file[n].mapset);
 
     G_fatal_error(_("Exit!"));
 }

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

@@ -123,14 +123,14 @@ block *get_block(struct cache * c, int idx)
 {
     int replace = G_lrand48() % c->nblocks;
     block *p = &c->blocks[replace];
-    int ref = c->refs[replace];
+    int cref = 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;
+    if (cref >= 0)
+	c->grid[cref] = NULL;
 
     c->grid[idx] = p;
     c->refs[replace] = idx;

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

@@ -9,8 +9,9 @@
 #include <unistd.h>
 #include "global.h"
 
-int rectify(char *name, char *mapset, struct cache *ebuffer,
-            double aver_z, char *result, char *interp_method)
+int rectify(struct Ortho_Image_Group *group, char *name, char *mapset,
+            struct cache *ebuffer, double aver_z, char *result,
+	    char *interp_method)
 {
     struct Cell_head cellhd;
     int ncols, nrows;
@@ -84,14 +85,14 @@ int rectify(char *name, char *mapset, struct cache *ebuffer,
 		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);
+	    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, 1);
+	    I_georef(ex1, nx1, &ex, &nx, group->E21, group->N21, 1);
 
 	    G_debug(5, "\t\tAfter geo ref: ex = %f \t nx =  %f", ex, nx);
 

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

@@ -1,10 +1,10 @@
 #include <grass/glocale.h>
 #include "global.h"
 
-int report(int rectify, int ok)
+int report(time_t rectify, int ok)
 {
     int minutes, hours;
-    int seconds;
+    time_t seconds;
     long ncells;
 
     G_message("%s", ok ? _("complete") : _("failed"));
@@ -20,9 +20,9 @@ int report(int rectify, int ok)
     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:%02d hours"), hours, minutes, seconds % 60);
+	G_verbose_message(_("%d:%02d:%02d hours"), hours, minutes, (int)(seconds % 60));
     else
-	G_verbose_message(_("%d:%02d minutes"), minutes, seconds % 60);
+	G_verbose_message(_("%d:%02d minutes"), minutes, (int)(seconds % 60));
     if (seconds)
 	G_verbose_message(_("%.1f cells per minute"),
 			  (60.0 * ncells) / ((double)seconds));

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

@@ -31,5 +31,9 @@ int get_target(char *group)
     sprintf(buf, _("Mapset <%s> in target location <%s> - "), mapset, location);
     strcat(buf, stat == 0 ? _("permission denied") : _("not found"));
   error:
+    strcat(buf, "\n");
+    strcat(buf, _("Please run i.target for group "));
+    strcat(buf, group);
     G_fatal_error("%s", buf);
+    return 1;			/* never reached */
 }

imagery/i.rectify/README

@@ -1,3 +1,4 @@
 As i.rectify functionality is integrated,
 i.rectify2 is renamed to i.rectify
 
+keep in sync with i.ortho.rectify

imagery/i.rectify/bilinear.c

@@ -11,22 +11,20 @@
  */
 
 #include <math.h>
-#include <grass/gis.h>
-#include <grass/raster.h>
 #include "global.h"
 
-void p_bilinear(struct cache *ibuffer,	/* input buffer                  */
-		void *obufptr,	/* ptr in output buffer          */
-		int cell_type,	/* raster map type of obufptr    */
-		double *row_idx,	/* row index                     */
-		double *col_idx,	/* column index          */
-		struct Cell_head *cellhd	/* information of output map     */
+void p_bilinear(struct cache *ibuffer,	  /* input buffer                */
+		void *obufptr,		  /* ptr in output buffer        */
+		int cell_type,		  /* raster map type of obufptr  */
+		double *row_idx,	  /* row index                   */
+		double *col_idx,	  /* column index                */
+		struct Cell_head *cellhd  /* information of output map   */
     )
 {
     int row;			/* row indices for interp        */
     int col;			/* column indices for interp     */
     int i, j;
-    DCELL t, u;			/* intermediate slope            */
+    double t, u;		/* intermediate slope            */
     DCELL result;		/* result of interpolation       */
     DCELL c[2][2];
 

imagery/i.rectify/bilinear_f.c

@@ -8,17 +8,15 @@
  *  it falls back to nearest neighbor.
  */
 
-#include <grass/gis.h>
-#include <grass/raster.h>
 #include <math.h>
 #include "global.h"
 
-void p_bilinear_f(struct cache *ibuffer,	/* input buffer                  */
-		void *obufptr,	/* ptr in output buffer          */
-		int cell_type,	/* raster map type of obufptr    */
-		double *row_idx,	/* row index                     */
-		double *col_idx,	/* column index          */
-	    struct Cell_head *cellhd	/* cell header of input layer    */
+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 */

imagery/i.rectify/cp.c

@@ -1,47 +1,47 @@
 #include <stdlib.h>
 #include <string.h>
-#include <grass/glocale.h>
 #include "global.h"
 
-int get_control_points(char *group, int order	/* THIS HAS BEEN ADDED WITH THE CRS MODIFICATIONS */
-    )
+int get_control_points(struct Image_Group *group, int order)
 {
     char msg[200];
 
-    if (!I_get_control_points(group, &cp))
+    if (!I_get_control_points(group->name, &group->control_points))
 	exit(0);
 
     sprintf(msg, _("Control Point file for group <%s@%s> - "),
-	    group, G_mapset());
+	    group->name, G_mapset());
+
+    G_verbose_message(_("Computing equations..."));
 
     if (order == 0) {
-	switch (I_compute_georef_equations_tps(&cp, &E12_t, &N12_t, &E21_t, &N21_t)) {
+	switch (I_compute_georef_equations_tps(&group->control_points,
+	                                       &group->E12_t,
+					       &group->N12_t,
+					       &group->E21_t,
+					       &group->N21_t)) {
 	case 0:
-	    sprintf(&msg[strlen(msg)],
-		    _("Not enough active control points for thin plate spline."));
+	    strcat(msg, _("Not enough active control points for thin plate spline."));
 	    break;
 	case -1:
 	    strcat(msg, _("Poorly placed control points."));
 	    strcat(msg, _(" Can not generate the transformation equation."));
 	    break;
 	case -2:
-	    strcat(msg, _("Not enough memory to solve for transformation equation"));
+	    strcat(msg, _("Not enough memory to solve transformation equations."));
 	    break;
 	case -3:
-	    strcat(msg, _("Invalid order"));
+	    strcat(msg, _("Invalid order."));
 	    break;
 	default:
-	    /* COMMENTED OUT WHEN SUPPORT FOR 3rd ORDER WAS ADDED BY 'CRS'
-	       E12a = E12[0]; E12b = E12[1]; E12c = E12[2];
-	       N12a = N12[0]; N12b = N12[1]; N12c = N12[2];
-	       E21a = E21[0]; E21b = E21[1]; E21c = E21[2];
-	       N21a = N21[0]; N21b = N21[1]; N21c = N21[2];
-	     */
 	    return 1;
 	}
     }
     else {
-	switch (I_compute_georef_equations(&cp, E12, N12, E21, N21, order)) {
+	switch (I_compute_georef_equations(&group->control_points,
+	                                   group->E12, group->N12,
+					   group->E21, group->N21,
+					   order)) {
 	case 0:
 	    sprintf(&msg[strlen(msg)],
 		    _("Not enough active control points for current order, %d are required."),
@@ -52,22 +52,16 @@ int get_control_points(char *group, int order	/* THIS HAS BEEN ADDED WITH THE CR
 	    strcat(msg, _(" Can not generate the transformation equation."));
 	    break;
 	case -2:
-	    strcat(msg, _("Not enough memory to solve for transformation equation"));
+	    strcat(msg, _("Not enough memory to solve transformation equations."));
 	    break;
 	case -3:
 	    strcat(msg, _("Invalid order"));
 	    break;
 	default:
-	    /* COMMENTED OUT WHEN SUPPORT FOR 3rd ORDER WAS ADDED BY 'CRS'
-	       E12a = E12[0]; E12b = E12[1]; E12c = E12[2];
-	       N12a = N12[0]; N12b = N12[1]; N12c = N12[2];
-	       E21a = E21[0]; E21b = E21[1]; E21c = E21[2];
-	       N21a = N21[0]; N21b = N21[1]; N21c = N21[2];
-	     */
 	    return 1;
 	}
     }
     G_fatal_error("%s", msg);
 
-    return 0; /* G_fatal_error() calls exit() */
+    return 0;
 }

imagery/i.rectify/cubic.c

@@ -12,23 +12,21 @@
  *  NULL they where filled with their neighbor
  */
 
-#include <grass/gis.h>
-#include <grass/raster.h>
 #include <math.h>
 #include "global.h"
 
-void p_cubic(struct cache *ibuffer,	/* input buffer                  */
-	     void *obufptr,	/* ptr in output buffer          */
-	     int cell_type,	/* raster map type of obufptr    */
-	     double *row_idx,	/* row index (decimal)           */
-	     double *col_idx,	/* column index (decimal)        */
-	     struct Cell_head *cellhd	/* information of output map     */
+void p_cubic(struct cache *ibuffer,    /* input buffer                */
+	     void *obufptr,	       /* ptr in output buffer        */
+	     int cell_type,	       /* raster map type of obufptr  */
+	     double *row_idx,	       /* row index (decimal)         */
+	     double *col_idx,	       /* column index (decimal)      */
+	     struct Cell_head *cellhd  /* information of output map   */
     )
 {
     int row;			/* row indices for interp        */
     int col;			/* column indices for interp     */
     int i, j;
-    DCELL t, u;			/* intermediate slope            */
+    double t, u;		/* intermediate slope            */
     DCELL result;		/* result of interpolation       */
     DCELL val[4];		/* buffer for temporary values   */
     DCELL cell[4][4];

imagery/i.rectify/cubic_f.c

@@ -9,17 +9,15 @@
  *  it falls back to nearest neighbor.
  */
 
-#include <grass/gis.h>
-#include <grass/raster.h>
 #include <math.h>
 #include "global.h"
 
-void p_cubic_f(struct cache *ibuffer,	/* input buffer                  */
-		void *obufptr,	/* ptr in output buffer          */
-		int cell_type,	/* raster map type of obufptr    */
-		double *row_idx,	/* row index                     */
-		double *col_idx,	/* column index          */
-	    struct Cell_head *cellhd	/* cell header of input layer    */
+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 */

imagery/i.rectify/defs.h

@@ -0,0 +1,42 @@
+/* cache for raster data, code taken from r.proj */
+
+/* These next defines determine the size of the sub-window that will
+ * be held in memory.  Larger values will require
+ * more memory (but less i/o) If you increase these values, keep  in
+ * mind that although you think the i/o will decrease, system paging
+ * (which goes on behind the scenes) may actual increase the i/o.
+ */
+
+#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.rectify/env.c

@@ -1,5 +1,6 @@
 #include <unistd.h>
 #include "global.h"
+
 static int which_env = -1;	/* 0 = cur, 1 = target */
 
 int select_current_env(void)

imagery/i.rectify/exec.c

@@ -1,8 +1,10 @@
 /* CMD mode from Bob Covill 2001 
+ * exec.c --
  *
- * small fixes: MN
+ * Loop through all files to be rectified and do the retification.
+ * Handles things like support files.
  *
- * Bug left: extension overwrites input name 1/2002
+ * small fixes: MN
  */
 
 #include <stdlib.h>
@@ -14,14 +16,10 @@
 #include <time.h>
 #include <unistd.h>
 #include <math.h>
-
-#include <grass/raster.h>
-#include <grass/glocale.h>
-
 #include "global.h"
 
-int exec_rectify(int order, char *extension, char *interp_method)
-/* ADDED WITH CRS MODIFICATIONS */
+int exec_rectify(struct Image_Group *group, int *ref_list, char *extension,
+                 char *interp_method, int order)
 {
     char *name;
     char *mapset;
@@ -38,53 +36,53 @@ int exec_rectify(int order, char *extension, char *interp_method)
     G_message("-----------------------------------------------");
 
     /* rectify each file */
-    for (n = 0; n < ref.nfiles; n++) {
-	if (ref_list[n]) {
-	    name = ref.file[n].name;
-	    mapset = ref.file[n].mapset;
+    for (n = 0; n < group->ref.nfiles; n++) {
+	if (!ref_list[n])
+	    continue;
 
-	    /* generate out name, add extension to output */
-	    result =
-		G_malloc(strlen(ref.file[n].name) + strlen(extension) + 1);
-	    strcpy(result, ref.file[n].name);
-	    strcat(result, extension);
+	name = group->ref.file[n].name;
+	mapset = group->ref.file[n].mapset;
 
-	    select_current_env();
+	/* generate out name, add extension to output */
+	result =
+	    G_malloc(strlen(group->ref.file[n].name) + strlen(extension) + 1);
+	strcpy(result, group->ref.file[n].name);
+	strcat(result, extension);
 
-	    cats_ok = Rast_read_cats(name, mapset, &cats) >= 0;
-	    colr_ok = Rast_read_colors(name, mapset, &colr) > 0;
+	select_current_env();
 
-	    /* Initialze History */
-	    if (Rast_read_history(name, mapset, &hist) < 0)
-		Rast_short_history(result, type, &hist);
+	cats_ok = Rast_read_cats(name, mapset, &cats) >= 0;
+	colr_ok = Rast_read_colors(name, mapset, &colr) > 0;
 
-	    time(&start_time);
+	/* Initialze History */
+	if (Rast_read_history(name, mapset, &hist) < 0)
+	    Rast_short_history(result, type, &hist);
 
-	    if (rectify(name, mapset, result, order, interp_method)) {
-		select_target_env();
+	time(&start_time);
 
-		if (cats_ok) {
-		    Rast_write_cats(result, &cats);
-		    Rast_free_cats(&cats);
-		}
-		if (colr_ok) {
-		    Rast_write_colors(result, G_mapset(), &colr);
-		    Rast_free_colors(&colr);
-		}
+	if (rectify(group, name, mapset, result, order, interp_method)) {
+	    select_target_env();
 
-		/* Write out History */
-		Rast_command_history(&hist);
-		Rast_write_history(result, &hist);
-
-		select_current_env();
-		time(&rectify_time);
-		report(rectify_time - start_time, 1);
+	    if (cats_ok) {
+		Rast_write_cats(result, &cats);
+		Rast_free_cats(&cats);
+	    }
+	    if (colr_ok) {
+		Rast_write_colors(result, G_mapset(), &colr);
+		Rast_free_colors(&colr);
 	    }
-	    else
-		report((time_t)0, 0);
+	    /* Write out History */
+	    Rast_command_history(&hist);
+	    Rast_write_history(result, &hist);
 
-	    G_free(result);
+	    select_current_env();
+	    time(&rectify_time);
+	    report(rectify_time - start_time, 1);
 	}
+	else
+	    report(0, 0);
+
+	G_free(result);
     }
 
     return 0;

imagery/i.rectify/get_wind.c

@@ -1,28 +1,82 @@
+#include <stdlib.h>
+#include <string.h>
 #include <math.h>
-#include <grass/glocale.h>
 #include "global.h"
 
-int georef_window(struct Cell_head *w1, struct Cell_head *w2, int order, double res)
+int get_ref_window(struct Ref *ref, int *ref_list, 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 < ref->nfiles; i++) {
+	if (!ref_list[i])
+	    continue;
+
+	if (count++ == 0) {
+	    Rast_get_cellhd(ref->file[i].name,
+			    ref->file[i].mapset, cellhd);
+	}
+	else {
+	    Rast_get_cellhd(ref->file[i].name,
+			    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 Image_Group *group, struct Cell_head *w1, struct Cell_head *w2, int order, double res)
 {
     double n, e, ad;
     struct _corner {
         double n, e;
     } nw, ne, se, sw;
 
-    /* extends */
+    /* extents */
+
+    /* compute geo ref of all corners */
     if (order == 0)
-	I_georef_tps(w1->west, w1->north, &e, &n, E12_t, N12_t, &cp, 1);
+	I_georef_tps(w1->west, w1->north, &e, &n, group->E12_t,
+	             group->N12_t, &group->control_points, 1);
     else
-	I_georef(w1->west, w1->north, &e, &n, E12, N12, order);
+	I_georef(w1->west, w1->north, &e, &n, group->E12, group->N12, order);
     w2->north = w2->south = n;
     w2->west = w2->east = e;
     nw.n = n;
     nw.e = e;
 
     if (order == 0)
-	I_georef_tps(w1->east, w1->north, &e, &n, E12_t, N12_t, &cp, 1);
+	I_georef_tps(w1->east, w1->north, &e, &n, group->E12_t,
+	             group->N12_t, &group->control_points, 1);
     else
-	I_georef(w1->east, w1->north, &e, &n, E12, N12, order);
+	I_georef(w1->east, w1->north, &e, &n, group->E12, group->N12, order);
     ne.n = n;
     ne.e = e;
     if (n > w2->north)
@@ -35,9 +89,10 @@ int georef_window(struct Cell_head *w1, struct Cell_head *w2, int order, double
 	w2->west = e;
 
     if (order == 0)
-	I_georef_tps(w1->west, w1->south, &e, &n, E12_t, N12_t, &cp, 1);
+	I_georef_tps(w1->west, w1->south, &e, &n, group->E12_t,
+	             group->N12_t, &group->control_points, 1);
     else
-	I_georef(w1->west, w1->south, &e, &n, E12, N12, order);
+	I_georef(w1->west, w1->south, &e, &n, group->E12, group->N12, order);
     sw.n = n;
     sw.e = e;
     if (n > w2->north)
@@ -50,9 +105,10 @@ int georef_window(struct Cell_head *w1, struct Cell_head *w2, int order, double
 	w2->west = e;
 
     if (order == 0)
-	I_georef_tps(w1->east, w1->south, &e, &n, E12_t, N12_t, &cp, 1);
+	I_georef_tps(w1->east, w1->south, &e, &n, group->E12_t,
+	             group->N12_t, &group->control_points, 1);
     else
-	I_georef(w1->east, w1->south, &e, &n, E12, N12, order);
+	I_georef(w1->east, w1->south, &e, &n, group->E12, group->N12, order);
     se.n = n;
     se.e = e;
     if (n > w2->north)
@@ -108,7 +164,7 @@ int georef_window(struct Cell_head *w1, struct Cell_head *w2, int order, double
 	}
     }
 
-    /* adjust extends */
+    /* adjust extents */
     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;

imagery/i.rectify/global.h

@@ -1,116 +1,24 @@
-/* These next defines determine the size of the sub-window that will
- * be held in memory.  Larger values will require
- * more memory (but less i/o) If you increase these values, keep  in
- * mind that although you think the i/o will decrease, system paging
- * (which goes on behind the scenes) may actual increase the i/o.
- */
-#include <time.h>
 #include <grass/gis.h>
 #include <grass/imagery.h>
+#include <grass/glocale.h>
+#include "defs.h"
 
-#define L2BDIM 6
-#define BDIM (1<<(L2BDIM))
-#define L2BSIZE (2*(L2BDIM))
-#define BSIZE (1<<(L2BSIZE))
-#define HI(i) ((i)>>(L2BDIM))
-#define LO(i) ((i)&((BDIM)-1))
+extern func interpolate;	/* interpolation routine */
 
-typedef DCELL block[BDIM][BDIM];
+extern int seg_mb_img;
 
-struct cache
+struct Image_Group
 {
-    int fd;
-    int stride;
-    int nblocks;
-    block **grid;
-    block *blocks;
-    int *refs;
+    char name[GNAME_MAX];
+    struct Ref ref;
+    struct Control_Points control_points;
+    int equation_stat;
+    /* georef coefficients */
+    double E12[10], N12[10], E21[10], N21[10];
+    /* TPS coefficients */
+    double *E12_t, *N12_t, *E21_t, *N21_t;
 };
 
-extern char *seg_mb;
-
-extern RASTER_MAP_TYPE map_type;
-extern int *ref_list;
-extern struct Ref ref;
-
-typedef void (*func) (struct cache *, void *, int, double *, double *, struct Cell_head *);
-
-extern func interpolate;		/* interpolation routine        */
-
-struct menu
-{
-    func method;		/* routine to interpolate new value      */
-    char *name;			/* method name                           */
-    char *text;			/* menu display - full description       */
-};
-
-/* georef coefficients */
-
-extern double E12[10], N12[10];
-extern double E21[10], N21[10];
-
-extern double *E12_t, *N12_t;
-extern double *E21_t, *N21_t;
-
-extern struct Control_Points cp;
-
-/* DELETED WITH CRS MODIFICATIONS
-   extern double E12a, E12b, E12c, N12a, N12b, N12c;
-   extern double E21a, E21b, E21c, N21a, N21b, N21c;
- */
 extern struct Cell_head target_window;
 
-/* cp.c */
-int get_control_points(char *, int);
-
-/* env.c */
-int select_current_env(void);
-int select_target_env(void);
-int show_env(void);
-
-/* exec.c */
-int exec_rectify(int, char *, char *);
-
-/* get_wind.c */
-int get_target_window(int);
-int georef_window(struct Cell_head *, struct Cell_head *, int, double);
-
-/* rectify.c */
-int rectify(char *, char *, char *, int, char *);
-
-/* readcell.c */
-extern struct cache *readcell(int, const char *);
-extern block *get_block(struct cache *, int);
-
-#define BKIDX(c,y,x) ((y) * (c)->stride + (x))
-#define BKPTR(c,y,x) ((c)->grid[BKIDX((c),(y),(x))])
-#define BLOCK(c,y,x) (BKPTR((c),(y),(x)) ? BKPTR((c),(y),(x)) : get_block((c),BKIDX((c),(y),(x))))
-#define CPTR(c,y,x) (&(*BLOCK((c),HI((y)),HI((x))))[LO((y))][LO((x))])
-
-/* report.c */
-int report(time_t, int);
-
-/* target.c */
-int get_target(char *);
-
-/* declare resampling methods */
-/* bilinear.c */
-extern void p_bilinear(struct cache *, void *, int, double *, double *,
-		       struct Cell_head *);
-/* cubic.c */
-extern void p_cubic(struct cache *, void *, int, double *, double *,
-		    struct Cell_head *);
-/* nearest.c */
-extern void p_nearest(struct cache *, void *, int, double *, double *,
-		      struct Cell_head *);
-/* bilinear_f.c */
-extern void p_bilinear_f(struct cache *, void *, int, double *, double *,
-		       struct Cell_head *);
-/* cubic_f.c */
-extern void p_cubic_f(struct cache *, void *, int, double *, double *,
-		    struct Cell_head *);
-/* lanczos.c */
-extern void p_lanczos(struct cache *, void *, int, double *, double *,
-		    struct Cell_head *);
-extern void p_lanczos_f(struct cache *, void *, int, double *, double *,
-		    struct Cell_head *);
+#include "local_proto.h"

imagery/i.rectify/i.rectify.html

@@ -66,14 +66,14 @@ The program will calculate the RMSE and check the required number of points.
 <h4>Linear affine transformation (1st order transformation)</h4>
 
 <dl>
-	<dd> x' = ax + by +c
-	<dd> y' = Ax + Bt +C
+	<dd> x' = ax + by + c
+	<dd> y' = Ax + By + C
 </dl>
 
 The a,b,c,A,B,C are determined by least squares regression
 based on the control points entered.  This transformation
 applies scaling, translation and rotation. It is NOT a
-general purpose rubber-sheeting, nor is it ortho-photo
+general purpose rubber-sheeting like TPS, nor is it ortho-photo
 rectification using a DEM, not second order polynomial,
 etc. It can be used if (1) you have geometrically correct
 images, and (2) the terrain or camera distortion effect can

imagery/i.rectify/local_proto.h

@@ -0,0 +1,51 @@
+
+/* cp.c */
+int get_control_points(struct Image_Group *, int);
+
+/* env.c */
+int select_current_env(void);
+int select_target_env(void);
+int show_env(void);
+
+/* exec.c */
+int exec_rectify(struct Image_Group *, int *, char *, char *, int);
+
+/* get_wind.c */
+int get_ref_window(struct Ref *, int *, struct Cell_head *);
+int georef_window(struct Image_Group *, struct Cell_head *, struct Cell_head *, int, double);
+
+/* rectify.c */
+int rectify(struct Image_Group *, char *, char *, char *, int, char *);
+
+/* readcell.c */
+struct cache *readcell(int, int);
+block *get_block(struct cache *, int);
+void release_cache(struct cache *);
+
+/* report.c */
+int report(time_t, int);
+
+/* target.c */
+int get_target(char *);
+
+/* declare resampling methods */
+/* bilinear.c */
+extern void p_bilinear(struct cache *, void *, int, double *, double *,
+		       struct Cell_head *);
+/* cubic.c */
+extern void p_cubic(struct cache *, void *, int, double *, double *,
+		    struct Cell_head *);
+/* nearest.c */
+extern void p_nearest(struct cache *, void *, int, double *, double *,
+		      struct Cell_head *);
+/* bilinear_f.c */
+extern void p_bilinear_f(struct cache *, void *, int, double *, double *,
+		       struct Cell_head *);
+/* cubic_f.c */
+extern void p_cubic_f(struct cache *, void *, int, double *, double *,
+		    struct Cell_head *);
+/* lanczos.c */
+extern void p_lanczos(struct cache *, void *, int, double *, double *,
+		    struct Cell_head *);
+extern void p_lanczos_f(struct cache *, void *, int, double *, double *,
+		    struct Cell_head *);

imagery/i.rectify/main.c

@@ -25,37 +25,15 @@
 
 #include <stdlib.h>
 #include <string.h>
-
-#include <grass/raster.h>
-#include <grass/glocale.h>
-
 #include "global.h"
 
-char *seg_mb;
-
-RASTER_MAP_TYPE map_type;
-int *ref_list;
-struct Ref ref;
+int seg_mb_img;
 
 func interpolate;
 
-/* georef coefficients */
-
-double E12[10], N12[10];
-double E21[10], N21[10];
-
-double *E12_t, *N12_t;
-double *E21_t, *N21_t;
-
-struct Control_Points cp;
-
-/* DELETED WITH CRS MODIFICATIONS
-   double E12a, E12b, E12c, N12a, N12b, N12c;
-   double E21a, E21b, E21c, N21a, N21b, N21c;
- */
 struct Cell_head target_window;
 
-void err_exit(char *, char *);
+void err_exit(struct Ref *, char *, char *);
 
 /* modify this table to add new methods */
 struct menu menu[] = {
@@ -63,7 +41,7 @@ struct menu menu[] = {
     {p_bilinear, "linear", "linear interpolation"},
     {p_cubic, "cubic", "cubic convolution"},
     {p_lanczos, "lanczos", "lanczos filter"},
-    {p_bilinear_f, "linear_f", "linear  interpolation with fallback"},
+    {p_bilinear_f, "linear_f", "linear interpolation with fallback"},
     {p_cubic_f, "cubic_f", "cubic convolution with fallback"},
     {p_lanczos_f, "lanczos_f", "lanczos filter with fallback"},
     {NULL, NULL, NULL}
@@ -73,13 +51,19 @@ static char *make_ipol_list(void);
 
 int main(int argc, char *argv[])
 {
-    char group[INAME_LEN], extension[INAME_LEN];
+    char extension[INAME_LEN];
     int order;			/* ADDED WITH CRS MODIFICATIONS */
     char *ipolname;		/* name of interpolation method */
     int method;
-    int n, i, m, k = 0;
+    char *seg_mb;
+    int i, m, k = 0;
     int got_file = 0, target_overwrite = 0;
     char *overstr;
+
+    struct Image_Group group;
+    int *ref_list;
+    int n;
+
     struct Cell_head cellhd;
 
     struct Option *grp,         /* imagery group */
@@ -174,7 +158,7 @@ int main(int argc, char *argv[])
     interpolate = menu[method].method;
 
     G_strip(grp->answer);
-    strcpy(group, grp->answer);
+    strcpy(group.name, grp->answer);
     strcpy(extension, ext->answer);
     order = atoi(val->answer);
 
@@ -189,10 +173,7 @@ int main(int argc, char *argv[])
 
     /* Find out how many files on command line */
     if (!a->answer) {
-	for (m = 0; ifile->answers[m]; m++) {
-	    k = m;
-	}
-	k++;
+	for (k = 0; ifile->answers[k]; k++);
     }
 
     if (!t->answer && (order < 1 || order > 3))  /* MAXORDER in lib/imagery/georef.c */
@@ -201,30 +182,36 @@ int main(int argc, char *argv[])
     if (t->answer)
 	order = 0;
 
+    /* find group */
+    if (!I_find_group(group.name)) {
+	G_fatal_error(_("Group <%s> not found"), group.name);
+    }
+
     /* determine the number of files in this group */
-    if (I_get_group_ref(group, &ref) <= 0) {
+    if (!I_get_group_ref(group.name, &group.ref)) {
 	G_warning(_("Location: %s"), G_location());
 	G_warning(_("Mapset: %s"), G_mapset());
-	G_fatal_error(_("Group <%s> does not exist"), grp->answer);
+	G_fatal_error(_("Could not read REF file for group <%s>"),
+	              group.name);
     }
 
-    if (ref.nfiles <= 0) {
+    if (group.ref.nfiles <= 0) {
 	G_important_message(_("Group <%s> contains no raster maps; run i.group"),
 			    grp->answer);
 	exit(EXIT_SUCCESS);
     }
 
-    ref_list = (int *)G_malloc(ref.nfiles * sizeof(int));
+    ref_list = (int *)G_malloc(group.ref.nfiles * sizeof(int));
 
     if (a->answer) {
-	for (n = 0; n < ref.nfiles; n++) {
+	for (n = 0; n < group.ref.nfiles; n++) {
 	    ref_list[n] = 1;
 	}
     }
     else {
 	char xname[GNAME_MAX], xmapset[GMAPSET_MAX], *name, *mapset;
 
-	for (n = 0; n < ref.nfiles; n++)
+	for (n = 0; n < group.ref.nfiles; n++)
 		ref_list[n] = 0;
 
 	for (m = 0; m < k; m++) {
@@ -239,17 +226,17 @@ int main(int argc, char *argv[])
 	    }
 
 	    got_file = 0;
-	    for (n = 0; n < ref.nfiles; n++) {
+	    for (n = 0; n < group.ref.nfiles; n++) {
 		if (mapset) {
-		    if (strcmp(name, ref.file[n].name) == 0 &&
-		        strcmp(mapset, ref.file[n].mapset) == 0) {
+		    if (strcmp(name, group.ref.file[n].name) == 0 &&
+		        strcmp(mapset, group.ref.file[n].mapset) == 0) {
 			got_file = 1;
 			ref_list[n] = 1;
 			break;
 		    }
 		}
 		else {
-		    if (strcmp(name, ref.file[n].name) == 0) {
+		    if (strcmp(name, group.ref.file[n].name) == 0) {
 			got_file = 1;
 			ref_list[n] = 1;
 			break;
@@ -257,15 +244,15 @@ int main(int argc, char *argv[])
 		}
 	    }
 	    if (got_file == 0)
-		err_exit(ifile->answers[m], group);
+		err_exit(&group.ref, ifile->answers[m], group.name);
 	}
     }
 
     /* read the control points for the group */
-    get_control_points(group, order);
+    get_control_points(&group, order);
 
     /* get the target */
-    get_target(group);
+    get_target(group.name);
 
     /* Check the GRASS_OVERWRITE environment variable */
     if ((overstr = getenv("GRASS_OVERWRITE")))  /* OK ? */
@@ -276,11 +263,11 @@ int main(int argc, char *argv[])
 	char result[GNAME_MAX];
 	
 	select_target_env();
-	for (i = 0; i < ref.nfiles; i++) {
+	for (i = 0; i < group.ref.nfiles; i++) {
 	    if (!ref_list[i])
 		continue;
 
-	    strcpy(result, ref.file[i].name);
+	    strcpy(result, group.ref.file[i].name);
 	    strcat(result, extension);
 	    
 	    if (G_legal_filename(result) < 0)
@@ -308,19 +295,22 @@ int main(int argc, char *argv[])
 	    if (!((res = atof(tres->answer)) > 0))
 		G_warning(_("Target resolution must be > 0, ignored"));
 	}
-	/* Calculate smallest region */
-	if (a->answer)
-	    Rast_get_cellhd(ref.file[0].name, ref.file[0].mapset, &cellhd);
-	else
-	    Rast_get_cellhd(ifile->answers[0], ref.file[0].mapset, &cellhd);
-
-	georef_window(&cellhd, &target_window, order, res);
+	/* get reference window from imagery group */
+	get_ref_window(&group.ref, ref_list, &cellhd);
+	georef_window(&group, &cellhd, &target_window, order, 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);
 
-    exec_rectify(order, extension, interpol->answer);
+    /* determine memory for elevation and imagery */
+    seg_mb_img = -1;
+    if (seg_mb) {
+	seg_mb_img = atoi(seg_mb);
+    }
+
+    /* go do it */
+    exec_rectify(&group, ref_list, extension, interpol->answer, order);
 
     G_done_msg(" ");
 
@@ -328,7 +318,7 @@ int main(int argc, char *argv[])
 }
 
 
-void err_exit(char *file, char *grp)
+void err_exit(struct Ref *ref, char *file, char *grp)
 {
     int n;
 
@@ -336,8 +326,8 @@ void err_exit(char *file, char *grp)
 	    file, grp);
     G_message(_("Try:"));
 
-    for (n = 0; n < ref.nfiles; n++)
-	G_message("%s", ref.file[n].name);
+    for (n = 0; n < ref->nfiles; n++)
+	G_message("%s@%s", ref->file[n].name, ref->file[n].mapset);
 
     G_fatal_error(_("Exit!"));
 }

imagery/i.rectify/nearest.c

@@ -1,20 +1,17 @@
-
 /*
  *      nearest.c - returns the nearest neighbor to a given
  *                  x,y position
  */
 
 #include <math.h>
-#include <grass/gis.h>
-#include <grass/raster.h>
 #include "global.h"
 
-void p_nearest(struct cache *ibuffer,	/* input buffer                  */
-	       void *obufptr,	/* ptr in output buffer          */
-	       int cell_type,	/* raster map type of obufptr    */
-	       double *row_idx,	/* row index in input matrix     */
-	       double *col_idx,	/* column index in input matrix  */
-	       struct Cell_head *cellhd	/* cell header of input layer    */
+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   */

imagery/i.rectify/readcell.c

@@ -10,12 +10,9 @@
 #include <sys/stat.h>
 #include <fcntl.h>
 #include <unistd.h>
-#include <grass/gis.h>
-#include <grass/raster.h>
-#include <grass/glocale.h>
 #include "global.h"
 
-struct cache *readcell(int fdi, const char *size)
+struct cache *readcell(int fdi, int size)
 {
     DCELL *tmpbuf;
     struct cache *c;
@@ -35,8 +32,8 @@ struct cache *readcell(int fdi, const char *size)
     ny = (nrows + BDIM - 1) / BDIM;
     nx = (ncols + BDIM - 1) / BDIM;
 
-    if (size)
-	nblocks = atoi(size) * ((1 << 20) / sizeof(block));
+    if (size > 0)
+	nblocks = size * ((1 << 20) / sizeof(block));
     else
 	nblocks = (nx + ny) * 2;	/* guess */
 
@@ -60,6 +57,8 @@ struct cache *readcell(int fdi, const char *size)
     }
     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);
@@ -126,7 +125,17 @@ block *get_block(struct cache * c, int idx)
 	G_fatal_error(_("Error seeking on segment file"));
 
     if (read(c->fd, p, sizeof(block)) < 0)
-	G_fatal_error(_("Error writing segment file"));
+	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.rectify/rectify.c

@@ -1,10 +1,4 @@
-#include <unistd.h>
-#include <string.h>
-
-#include <grass/raster.h>
-#include <grass/glocale.h>
-
-#include "global.h"
+/* rectification code */
 
 /* Modified to support Grass 5.0 fp format 11 april 2000
  *
@@ -12,13 +6,21 @@
  *
  */
 
-int rectify(char *name, char *mapset, char *result, int order, char *interp_method)
+#include <stdlib.h>
+#include <string.h>
+#include <unistd.h>
+#include "global.h"
+
+int rectify(struct Image_Group *group, char *name, char *mapset,
+            char *result, int order, char *interp_method)
 {
     struct Cell_head cellhd;
     int ncols, nrows;
     int row, col;
     double row_idx, col_idx;
-    int infd, cell_size, outfd;
+    int infd, outfd;
+    RASTER_MAP_TYPE map_type;
+    int cell_size;
     void *trast, *tptr;
     double n1, e1, nx, ex;
     struct cache *ibuffer;
@@ -34,7 +36,7 @@ int rectify(char *name, char *mapset, char *result, int order, char *interp_meth
     map_type = Rast_get_map_type(infd);
     cell_size = Rast_cell_size(map_type);
 
-    ibuffer = readcell(infd, seg_mb);
+    ibuffer = readcell(infd, seg_mb_img);
 
     Rast_close(infd);		/* (pmx) 17 april 2000 */
 
@@ -73,9 +75,9 @@ int rectify(char *name, char *mapset, char *result, int order, char *interp_meth
 
 	    /* backwards transformation of target cell center */
 	    if (order == 0)
-		I_georef_tps(e1, n1, &ex, &nx, E21_t, N21_t, &cp, 0);
+		I_georef_tps(e1, n1, &ex, &nx, group->E21_t, group->N21_t, &group->control_points, 0);
 	    else
-		I_georef(e1, n1, &ex, &nx, E21, N21, order);
+		I_georef(e1, n1, &ex, &nx, group->E21, group->N21, order);
 
 	    /* convert to row/column indices of source raster */
 	    row_idx = (cellhd.north - nx) / cellhd.ns_res;
@@ -94,7 +96,7 @@ int rectify(char *name, char *mapset, char *result, int order, char *interp_meth
     G_free(trast);
 
     close(ibuffer->fd);
-    G_free(ibuffer);
+    release_cache(ibuffer);
 
     Rast_get_cellhd(result, G_mapset(), &cellhd);
 

imagery/i.rectify/report.c

@@ -1,4 +1,3 @@
-#include <time.h>
 #include <grass/glocale.h>
 #include "global.h"
 

imagery/i.rectify/target.c

@@ -1,6 +1,5 @@
-#include <string.h>
 #include <unistd.h>
-#include <grass/glocale.h>
+#include <string.h>
 #include "global.h"
 
 int get_target(char *group)