i.rectify: adjust to r.proj: transform cell center coords not cell border coords, use fast cache, offer different resampling methods

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

imagery/i.rectify/bilinear.c

@@ -0,0 +1,68 @@
+/*
+ * Name
+ *  bilinear.c -- use bilinear interpolation for given row, col
+ *
+ * Description
+ *  bilinear interpolation for the given row, column indices.
+ *  If the given row or column is outside the bounds of the input map,
+ *  the point in the output map is set to NULL.
+ *  If any of the 4 surrounding points to be used in the interpolation
+ *  is NULL it is filled with is neighbor value
+ */
+
+#include <math.h>
+#include <grass/gis.h>
+#include <grass/raster.h>
+#include "global.h"
+
+void p_bilinear(struct cache *ibuffer,	/* input buffer                  */
+		void *obufptr,	/* ptr in output buffer          */
+		int cell_type,	/* raster map type of obufptr    */
+		double *row_idx,	/* row index                     */
+		double *col_idx,	/* column index          */
+		struct Cell_head *cellhd	/* information of output map     */
+    )
+{
+    int row0,			/* lower row index for interp    */
+      row1,			/* upper row index for interp    */
+      col0,			/* lower column index for interp */
+      col1;			/* upper column index for interp */
+    DCELL t, u,			/* intermediate slope            */
+      tu,			/* t * u                         */
+      result;			/* result of interpolation       */
+    DCELL *c00, *c01, *c10, *c11;
+
+    /* cut indices to integer */
+    row0 = (int)floor(*row_idx);
+    col0 = (int)floor(*col_idx);
+    row1 = row0 + 1;
+    col1 = col0 + 1;
+
+    /* check for out of bounds - if out of bounds set NULL value and return */
+    if (row0 < 0 || row1 >= cellhd->rows || col0 < 0 || col1 >= cellhd->cols) {
+	Rast_set_null_value(obufptr, 1, cell_type);
+	return;
+    }
+
+    c00 = CPTR(ibuffer, row0, col0);
+    c01 = CPTR(ibuffer, row0, col1);
+    c10 = CPTR(ibuffer, row1, col0);
+    c11 = CPTR(ibuffer, row1, col1);
+
+    /* check for NULL values */
+    if (Rast_is_f_null_value(c00) ||
+	Rast_is_f_null_value(c01) ||
+	Rast_is_f_null_value(c10) || Rast_is_f_null_value(c11)) {
+	Rast_set_null_value(obufptr, 1, cell_type);
+	return;
+    }
+
+    /* do the interpolation  */
+    t = *col_idx - col0;
+    u = *row_idx - row0;
+    tu = t * u;
+
+    result = Rast_interp_bilinear(t, u, *c00, *c01, *c10, *c11);
+
+    Rast_set_d_value(obufptr, result, cell_type);
+}

imagery/i.rectify/bilinear_f.c

@@ -0,0 +1,49 @@
+/*
+ * Name
+ *  bilinear_f.c -- use bilinear interpolation with fallback for given row, col
+ *
+ * Description
+ *  bilinear interpolation for the given row, column indices.
+ *  If the interpolated value (but not the nearest) is null,
+ *  it falls back to nearest neighbor.
+ */
+
+#include <grass/gis.h>
+#include <grass/raster.h>
+#include <math.h>
+#include "global.h"
+
+void p_bilinear_f(struct cache *ibuffer,	/* input buffer                  */
+		void *obufptr,	/* ptr in output buffer          */
+		int cell_type,	/* raster map type of obufptr    */
+		double *row_idx,	/* row index                     */
+		double *col_idx,	/* column index          */
+	    struct Cell_head *cellhd	/* cell header of input layer    */
+    )
+{
+    /* start nearest neighbor to do some basic tests */
+    int row, col;		/* row/col of nearest neighbor   */
+    DCELL *cellp;
+
+    /* cut indices to integer */
+    row = (int)floor(*row_idx);
+    col = (int)floor(*col_idx);
+
+    /* check for out of bounds - if out of bounds set NULL value     */
+    if (row < 0 || row >= cellhd->rows || col < 0 || col >= cellhd->cols) {
+        Rast_set_null_value(obufptr, 1, cell_type);
+        return;
+    }
+
+    cellp = CPTR(ibuffer, row, col);
+    /* if nearest is null, all the other interps will be null */
+    if (Rast_is_d_null_value(cellp)) {
+        Rast_set_null_value(obufptr, 1, cell_type);
+        return;
+    }
+    
+    p_bilinear(ibuffer, obufptr, cell_type, col_idx, row_idx, cellhd);
+    /* fallback to nearest if bilinear is null */
+    if (Rast_is_d_null_value(obufptr))
+        Rast_set_d_value(obufptr, *cellp, cell_type);
+}

imagery/i.rectify/cubic.c

@@ -0,0 +1,73 @@
+/*
+ * Name
+ *  cubic.c -- use cubic convolution interpolation for given row, col
+ *
+ * Description
+ *  cubic returns the value in the buffer that is the result of cubic
+ *  convolution interpolation for the given row, column indices.
+ *  If the given row or column is outside the bounds of the input map,
+ *  the corresponding point in the output map is set to NULL.
+ *
+ *  If single surrounding points in the interpolation matrix are
+ *  NULL they where filled with their neighbor
+ */
+
+#include <grass/gis.h>
+#include <grass/raster.h>
+#include <math.h>
+#include "global.h"
+
+void p_cubic(struct cache *ibuffer,	/* input buffer                  */
+	     void *obufptr,	/* ptr in output buffer          */
+	     int cell_type,	/* raster map type of obufptr    */
+	     double *row_idx,	/* row index (decimal)           */
+	     double *col_idx,	/* column index (decimal)        */
+	     struct Cell_head *cellhd	/* information of output map     */
+    )
+{
+    int row,			/* row indices for interp        */
+      col;			/* column indices for interp     */
+    int i, j;
+    DCELL t, u,			/* intermediate slope            */
+      result,			/* result of interpolation       */
+      val[4];			/* buffer for temporary values   */
+    DCELL *cellp[4][4];
+
+    /* cut indices to integer */
+    row = (int)floor(*row_idx);
+    col = (int)floor(*col_idx);
+
+    /* check for out of bounds of map - if out of bounds set NULL value     */
+    if (row - 1 < 0 || row + 2 >= cellhd->rows ||
+	col - 1 < 0 || col + 2 >= cellhd->cols) {
+	Rast_set_null_value(obufptr, 1, cell_type);
+	return;
+    }
+
+    for (i = 0; i < 4; i++)
+	for (j = 0; j < 4; j++)
+	    cellp[i][j] = CPTR(ibuffer, row - 1 + i, col - 1 + j);
+
+    /* check for NULL value                                         */
+    for (i = 0; i < 4; i++)
+	for (j = 0; j < 4; j++) {
+	    if (Rast_is_d_null_value(cellp[i][j])) {
+		Rast_set_null_value(obufptr, 1, cell_type);
+		return;
+	    }
+	}
+
+    /* do the interpolation  */
+    t = *col_idx - col;
+    u = *row_idx - row;
+
+    for (i = 0; i < 4; i++) {
+	DCELL **tmp = cellp[i];
+
+	val[i] = Rast_interp_cubic(t, *tmp[0], *tmp[1], *tmp[2], *tmp[3]);
+    }
+
+    result = Rast_interp_cubic(u, val[0], val[1], val[2], val[3]);
+
+    Rast_set_d_value(obufptr, result, cell_type);
+}

imagery/i.rectify/cubic_f.c

@@ -0,0 +1,54 @@
+/*
+ * Name
+ *  cubic_f.c -- use cubic interpolation with fallback for given row, col
+ *
+ * Description
+ *  cubic interpolation for the given row, column indices.
+ *  If the interpolated value (but not the nearest) is null,
+ *  it falls back to bilinear.  If that interp value is null,
+ *  it falls back to nearest neighbor.
+ */
+
+#include <grass/gis.h>
+#include <grass/raster.h>
+#include <math.h>
+#include "global.h"
+
+void p_cubic_f(struct cache *ibuffer,	/* input buffer                  */
+		void *obufptr,	/* ptr in output buffer          */
+		int cell_type,	/* raster map type of obufptr    */
+		double *row_idx,	/* row index                     */
+		double *col_idx,	/* column index          */
+	    struct Cell_head *cellhd	/* cell header of input layer    */
+    )
+{
+    /* start nearest neighbor to do some basic tests */
+    int row, col;		/* row/col of nearest neighbor   */
+    DCELL *cellp;
+
+    /* cut indices to integer */
+    row = (int)floor(*row_idx);
+    col = (int)floor(*col_idx);
+
+    /* check for out of bounds - if out of bounds set NULL value     */
+    if (row < 0 || row >= cellhd->rows || col < 0 || col >= cellhd->cols) {
+        Rast_set_null_value(obufptr, 1, cell_type);
+        return;
+    }
+
+    cellp = CPTR(ibuffer, row, col);
+    /* if nearest is null, all the other interps will be null */
+    if (Rast_is_d_null_value(cellp)) {
+        Rast_set_null_value(obufptr, 1, cell_type);
+        return;
+    }
+    
+    p_cubic(ibuffer, obufptr, cell_type, col_idx, row_idx, cellhd);
+    /* fallback to bilinear if cubic is null */
+    if (Rast_is_d_null_value(obufptr)) {
+        p_bilinear(ibuffer, obufptr, cell_type, col_idx, row_idx, cellhd);
+        /* fallback to nearest if bilinear is null */
+	    if (Rast_is_d_null_value(obufptr))
+		Rast_set_d_value(obufptr, *cellp, cell_type);
+    }
+}

imagery/i.rectify/exec.c

@@ -13,38 +13,32 @@
 #include <fcntl.h>
 #include <time.h>
 #include <unistd.h>
+#include <math.h>
 
 #include <grass/raster.h>
 #include <grass/glocale.h>
 
 #include "global.h"
 
-int exec_rectify(int order, char *extension)
+int exec_rectify(int order, char *extension, char *interp_method)
 /* ADDED WITH CRS MODIFICATIONS */
 {
     char *name;
     char *mapset;
     char *result;
-    char *type;
-    int i, n;
+    char *type = "raster";
+    int n;
     struct Colors colr;
     struct Categories cats;
     struct History hist;
     int colr_ok, cats_ok;
     long start_time, rectify_time, compress_time;
 
-
-    /* allocate the output cell matrix */
-    cell_buf = (void **)G_calloc(NROWS, sizeof(void *));
-    n = NCOLS * Rast_cell_size(map_type);
-    for (i = 0; i < NROWS; i++) {
-	cell_buf[i] = (void *)G_malloc(n);
-	Rast_set_null_value(cell_buf[i], NCOLS, map_type);
-    }
+    Rast_set_output_window(&target_window);
 
     /* rectify each file */
     for (n = 0; n < ref.nfiles; n++) {
-	if ((i = ref_list[n]) < 0) {
+	if (ref_list[n] < 0) {
 	    /* continue; */
 	    name = ref.file[n].name;
 	    mapset = ref.file[n].mapset;
@@ -63,20 +57,13 @@ int exec_rectify(int order, char *extension)
 	    colr_ok = Rast_read_colors(name, mapset, &colr) > 0;
 
 	    /* Initialze History */
-	    type = "raster";
 	    Rast_short_history(name, type, &hist);
 
 	    time(&start_time);
 
-	    if (rectify(name, mapset, result, order)) {
+	    if (rectify(name, mapset, result, order, interp_method)) {
 		select_target_env();
 
-	    /***
-	     * This clobbers (with wrong values) head
-	     * written by gislib.  99% sure it should
-	     * be removed.  EGM 2002/01/03
-            Rast_put_cellhd (result,&target_window);
-	     */
 		if (cats_ok) {
 		    Rast_write_cats(result, &cats);
 		    Rast_free_cats(&cats);
@@ -101,10 +88,10 @@ int exec_rectify(int order, char *extension)
 	    }
 	    else
 		report(name, mapset, result, (long)0, (long)0, 0);
+
+	    G_free(result);
 	}
     }
 
-    G_done_msg("");
-
     return 0;
 }

imagery/i.rectify/global.h

@@ -5,25 +5,30 @@
  * (which goes on behind the scenes) may actual increase the i/o.
  */
 #include <grass/gis.h>
+#include <grass/imagery.h>
 
-#define NROWS 64
-#define NCOLS 256
+#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))
 
-#include <grass/imagery.h>
-#include "rowcol.h"
+typedef DCELL block[BDIM][BDIM];
 
-#define IDX int
+struct cache
+{
+    int fd;
+    int stride;
+    int nblocks;
+    block **grid;
+    block *blocks;
+    int *refs;
+};
 
-extern ROWCOL row_map[NROWS][NCOLS];
-extern ROWCOL col_map[NROWS][NCOLS];
-extern ROWCOL row_min[NROWS];
-extern ROWCOL row_max[NROWS];
-extern ROWCOL row_left[NROWS];
-extern ROWCOL row_right[NROWS];
-extern IDX row_idx[NROWS];
-extern int matrix_rows, matrix_cols;
+extern char *seg_mb;
+extern int seg_rows, seg_cols;
 
-extern void **cell_buf;
 extern int temp_fd;
 extern RASTER_MAP_TYPE map_type;
 extern char *temp_name;
@@ -31,6 +36,17 @@ extern int *ref_list;
 extern char **new_name;
 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];
@@ -51,20 +67,23 @@ int select_target_env(void);
 int show_env(void);
 
 /* exec.c */
-int exec_rectify(int, char *);
+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);
 
-/* matrix.c */
-int compute_georef_matrix(struct Cell_head *, struct Cell_head *, int);
+/* rectify.c */
+int rectify(char *, char *, char *, int, char *);
 
-/* perform.c */
-int perform_georef(int, void *);
+/* readcell.c */
+extern struct cache *readcell(int, const char *);
+extern block *get_block(struct cache *, int);
 
-/* rectify.c */
-int rectify(char *, char *, char *, 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(char *, char *, char *, long, long, int);
@@ -72,6 +91,19 @@ int report(char *, char *, char *, long, long, int);
 /* target.c */
 int get_target(char *);
 
-/* write.c */
-int write_matrix(int, int);
-int write_map(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.rectify/i.rectify.html

@@ -7,10 +7,10 @@ points identified in
 or
 <EM><A HREF="i.vpoints.html">i.vpoints</A></EM>
 
-to calculate a transformation matrix based on a  first,
+to calculate a transformation matrix based on a first,
 second, or third order polynomial and then converts x,y
 cell coordinates to standard map coordinates for each pixel
-in the image.  The result is a planimetric image with a
+in the image. The result is a planimetric image with a
 transformed coordinate system (i.e., a different coordinate
 system than before it was rectified).
 
@@ -21,111 +21,55 @@ or
 <EM><A HREF="i.vpoints.html">i.vpoints</A></EM>
 
 must be run before <EM>i.rectify</EM>, and both programs
-are required to rectify an image.  An image must be
+are required to rectify an image. An image must be
 rectified before it can reside in a standard coordinate
 LOCATION, and therefore be analyzed with the other map
-layers in the standard coordinate LOCATION.  Upon
+layers in the standard coordinate LOCATION. Upon
 completion of <EM>i.rectify</EM>, the rectified image is
-deposited in the target standard coordinate LOCATION.  This
+deposited in the target standard coordinate LOCATION. This
 LOCATION is selected using
 
 <EM><A HREF="i.target.html">i.target</A></EM>.
 
-<H2>Program Prompts</H2>
-
-The first prompt in the program asks for the name of
-the group containing the files to be rectified.
-
-
-<PRE>
-     Enter the group containing files to be rectified
-     Enter 'list' for a list of existing imagery groups
-     Enter 'list -f' for a verbose listing
-     Hit RETURN to cancel request
-     &gt;
-</PRE>
-
- This is the same imagery group that was selected in 
-
-<EM><A HREF="i.points.html">i.points</A></EM>
-or
-<EM><A HREF="i.vpoints.html">i.vpoints</A></EM>
-
-and the group that contains the raster maps with the marked
-points and their associated map  coordinates.  You are then
-asked to select the raster map(s) within the group to be
-rectified:
-
-
-<PRE>
-Please select the file(s) to rectify by naming an output file
-
-       spot1.1 in mapsetname      .............
-       spot1.2 in mapsetname      .............
-       spot1.3 in mapsetname      .............
-       spotclass1 in mapsetname   spotrectify1.
-    
-       spotreject1 in mapsetname  .............
-
-(enter list by any name to get a list of existing raster maps)
-
-AFTER COMPLETING ALL ANSWERS, HIT &lt;ESC&gt; TO CONTINUE
-           (OR&lt;Ctrl-C&gt; TO CANCEL)
-</PRE>
-
-More than one raster map may be rectified at a time.  Each
-cell file should be given a unique output file name.
-
-
-<P>
-
-Next, you are asked to select one of two windows regions:
-
-
-<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>More than one raster map may be rectified at a time. Each cell file
+should be given a unique output file name. The rectified image or
+rectified raster maps will be located in the target LOCATION when the
+program is completed. The original unrectified files are not modified
+or removed.
 
-The <EM>i.rectify</EM> program will only rectify that
-portion of the image or raster map that occurs within the
-chosen window region, and only that portion of the cell
-file will be relocated in the target database.  It is
+<p>
+If the <b>-c</b> flag is used, <EM>i.rectify</EM> will only rectify that
+portion of the image or raster map that occurs within the chosen window
+region in the target location, and only that portion of the cell
+file will be relocated in the target database. It is
 important therefore, to check the current mapset window in
-the target LOCATION if choice number one is selected.
+the target LOCATION if the <b>-c</b> flag is used.
 
 <P>
 
 If you are rectifying a file with plans to patch it to
 another file using the GRASS program <em>r.patch</em>,
 choose option number one, the current window in the target
-location.  This window, however, must be the default window
-for the target LOCATION.  When a file being rectified is
+location. This window, however, must be the default window
+for the target LOCATION. When a file being rectified is
 smaller than the default window in which it is being
-rectified, zeros are added to the rectified file.  Patching
-files of the same size that contain 0/non-zero data,
-eliminates the possibility of a no-data line the patched
-result.  This is because, when the images are patched, the
-zeros in the image are "covered" with non-zero pixel
-values.  When rectifying files that are going to be
+rectified, NULLs are added to the rectified file. Patching
+files of the same size that contain NULL data,
+eliminates the possibility of a no-data line in the patched
+result. This is because, when the images are patched, the
+NULLs in the image are "covered" with non-NULL pixel
+values. When rectifying files that are going to be
 patched, rectify all of the files using the same default
 window.
 
-
+<h3>Coordinate transformation</h3>
 <P>
+The desired order of transformation (1, 2, or 3) is selected with the
+<b>order</b> option.
 
-Select the order of transformation desired with the <b>order</b> option:
-
-<PRE>
-Select order of transformation --&gt;   1st Order   2nd Order  3rd Order
-</PRE>
+The program will calculate the RMSE and check the required number of points.
 
-The program will immediately recalculate the RMSE and the
-number of points required.
-
-<h3>Linear affine transformation (1st order transformation)</h3>
+<h4>Linear affine transformation (1st order transformation)</h4>
 
 <DL>
 	<DD> x' = ax + by +c
@@ -134,21 +78,19 @@ number of points required.
 
 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
+applies scaling, translation and rotation. It is NOT a
 general purpose rubber-sheeting, nor is it ortho-photo
 rectification using a DEM, not second order polynomial,
-etc.  It can be used if (1) you have geometrically correct
+etc. It can be used if (1) you have geometrically correct
 images, and (2) the terrain or camera distortion effect can
 be ignored.
 
+<H4>Polynomial Transformation Matrix (2nd, 3d order transformation)</H4>
 
-<H3>Polynomial Transformation Matrix (2nd, 3d order transformation)</H3>
-
-The ANALYZE function has been changed to support
-calculating the registration coefficients using a first,
-second, or third order transformation matrix.  The number
-of control points required for a selected order of
-transformation (represented by n) is
+<EM>i.rectify</EM> uses a first, second, or third order transformation
+matrix to calculate the registration coefficients. The number
+of control points required for a selected order of transformation
+(represented by n) is
 
 <DL>
 <DD>((n + 1) * (n + 2) / 2) 
@@ -156,20 +98,44 @@ transformation (represented by n) is
 
 or 3, 6, and 10 respectively. It is strongly recommended
 that one or more additional points be identified to allow
-for an overly- determined transformation calculation which
+for an overly-determined transformation calculation which
 will generate the Root Mean Square (RMS) error values for
-each included point.  The RMS error values for all the
+each included point. The RMS error values for all the
 included control points are immediately recalculated when
 the user selects a different transformation order from the
-menu bar.  The polynomial equations are performed using a 
+menu bar. The polynomial equations are performed using a 
 modified Gaussian elimination method.
 
-<H3>Program Execution</H3>
-
-Note:  The rectified image or rectified raster maps will be
-located in the target LOCATION when the program is
-completed.  The original unrectified files are not modified
-or removed.
+<h3>Resampling method</h3>
+<p>
+The rectified data is resampled with one of five different methods: 
+<em>nearest</em>, <em>bilinear</em>, <em>cubic</em>, <em>bilinear_f</em>
+or <em>cubic_f</em>.
+<p>
+The <em>method=nearest</em> method, which performs a nearest neighbor assignment,
+is the fastest of the three resampling methods. It is primarily used for
+categorical data such as a land use classification, since it will not change
+the values of the data cells. The <em>method=bilinear</em> method determines the new
+value of the cell based on a weighted distance average of the 4 surrounding
+cells in the input map. The <em>method=cubic</em> method determines the new value of
+the cell based on a weighted distance average of the 16 surrounding cells in
+the input map.
+<p>
+The bilinear and cubic interpolation methods are most appropriate for
+continuous data and cause some smoothing. Both options should not be used
+with categorical data, since the cell values will be altered.
+<p>
+In the bilinear and cubic methods, if any of the surrounding cells used to
+interpolate the new cell value are NULL, the resulting cell will be NULL, even if
+the nearest cell is not NULL. This will cause some thinning along NULL borders,
+such as the coasts of land areas in a DEM. The bilinear_f and cubic_f
+interpolation methods can be used if thinning along NULL edges is not desired.
+These methods "fall back" to simpler interpolation methods along NULL borders.
+That is, from cubic to bilinear to nearest.
+<p>
+If nearest neighbor assignment is used, the output map has the same raster
+format as the input map. If any of the other interpolations is used, the
+output map is written as floating point.
 
 <P>
 <!--
@@ -180,19 +146,12 @@ mode.
 
 <h2>NOTES</h2>
 
-<EM>i.rectify</EM> uses nearest neighbor resampling during
-the transformation choosing the actual pixel that has its centre nearest to
-the point location in the image. Advantage of this method is that the pixel
-brightness of the image is kept as <EM>i.rectify</EM> rearranges the
-geometry of the image pixels.
-<P>
-
 If <em>i.rectify</em> starts normally but after some time the following text is seen:
 <br><tt>
-GIS ERROR: error while writing to temp file
+ERROR: Error writing segment file
 </tt><br>
-the user may try the flag <EM>-c</EM> (or the module needs more free space
-on the hard drive).
+the user may try the <b>-c</b> flag or the module needs more free space
+on the hard drive.
 
 
 <H2>SEE ALSO</H2>
@@ -210,8 +169,9 @@ Processing manual</A></EM>
 <A HREF="i.points.html">i.points</A>,
 <A HREF="i.vpoints.html">i.vpoints</A>,
 <A HREF="i.target.html">i.target</A>
-</EM><br>
-<em><a href="gm_georect.html">gis.m: GEORECTIFY TOOL</a></em>
+<br>
+<a href="wxGUI.GCP_Manager.html">Manage Ground Control Points</a>,
+<a href="gm_georect.html">gis.m: GEORECTIFY TOOL</a></em>
 
 
 <H2>AUTHORS</H2>

imagery/i.rectify/main.c

@@ -10,7 +10,8 @@
  *               Markus Neteler <neteler itc.it>, 
  *               Bernhard Reiter <bernhard intevation.de>, 
  *               Glynn Clements <glynn gclements.plus.com>, 
- *               Hamish Bowman <hamish_b yahoo.com>
+ *               Hamish Bowman <hamish_b yahoo.com>,
+ *               Markus Metz
  * PURPOSE:      calculate a transformation matrix and then convert x,y cell 
  *               coordinates to standard map coordinates for each pixel in the 
  *               image (control points can come from i.points or i.vpoints)
@@ -31,17 +32,9 @@
 #include "global.h"
 #include "crs.h"
 
+char *seg_mb;
+int seg_rows, seg_cols;
 
-ROWCOL row_map[NROWS][NCOLS];
-ROWCOL col_map[NROWS][NCOLS];
-ROWCOL row_min[NROWS];
-ROWCOL row_max[NROWS];
-ROWCOL row_left[NROWS];
-ROWCOL row_right[NROWS];
-IDX row_idx[NROWS];
-int matrix_rows, matrix_cols;
-
-void **cell_buf;
 int temp_fd;
 RASTER_MAP_TYPE map_type;
 char *temp_name;
@@ -49,6 +42,8 @@ int *ref_list;
 char **new_name;
 struct Ref ref;
 
+func interpolate;
+
 /* georef coefficients */
 
 double E12[10], N12[10];
@@ -60,19 +55,40 @@ double E21[10], N21[10];
  */
 struct Cell_head target_window;
 
-#define NFILES 15
+#define NFILES 15   /* ??? */
 
 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 group[INAME_LEN], extension[INAME_LEN];
     char result[NFILES][15];
     int order;			/* ADDED WITH CRS MODIFICATIONS */
+    char *ipolname;		/* name of interpolation method */
+    int method;
     int n, i, m, k = 0;
     int got_file = 0;
 
-    struct Option *grp, *val, *ifile, *ext, *tres;
+    struct Option *grp,         /* imagery group */
+     *val,                      /* transformation order */
+     *ifile,			/* input files */
+     *ext,			/* extension */
+     *tres,			/* target resolution */
+     *mem,			/* amount of memory for cache */
+     *interpol;			/* interpolation method:
+				   nearest neighbor, bilinear, cubic */
     struct Flag *c, *a;
     struct GModule *module;
 
@@ -115,6 +131,24 @@ int main(int argc, char *argv[])
     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");
+
     c = G_define_flag();
     c->key = 'c';
     c->description =
@@ -124,15 +158,30 @@ int main(int argc, char *argv[])
     a->key = 'a';
     a->description = _("Rectify all raster maps in group");
 
-
     if (G_parser(argc, argv))
 	exit(EXIT_FAILURE);
 
+    /* get the method */
+    for (method = 0; (ipolname = menu[method].name); method++)
+	if (strcmp(ipolname, interpol->answer) == 0)
+	    break;
+
+    if (!ipolname)
+	G_fatal_error(_("<%s=%s> unknown %s"),
+		      interpol->key, interpol->answer, interpol->key);
+    interpolate = menu[method].method;
+
     G_strip(grp->answer);
     strcpy(group, grp->answer);
     strcpy(extension, ext->answer);
     order = atoi(val->answer);
 
+    seg_mb = NULL;
+    if (mem->answer) {
+	if (atoi(mem->answer) > 0)
+	    seg_mb = mem->answer;
+    }
+
     if (!ifile->answers)
 	a->answer = 1;		/* force all */
 
@@ -148,6 +197,9 @@ int main(int argc, char *argv[])
 	G_fatal_error(_("Invalid order (%d); please enter 1 to %d"), order,
 		      MAXORDER);
 
+    if (!(seg_mb > 0))
+	G_fatal_error(_("Amount of memory to use in MB must be > 0"));
+
     /* determine the number of files in this group */
     if (I_get_group_ref(group, &ref) <= 0) {
 	G_warning(_("Location: %s"), G_location());
@@ -173,10 +225,17 @@ int main(int argc, char *argv[])
 	}
     }
     else {
+	char xname[GNAME_MAX], xmapset[GMAPSET_MAX], *name;
 	for (m = 0; m < k; m++) {
 	    got_file = 0;
+	    if (G_name_is_fully_qualified(ifile->answers[m], xname, xmapset))
+		name = xname;
+	    else
+		name = ifile->answers[m];
+
 	    for (n = 0; n < ref.nfiles; n++) {
-		if (strcmp(ifile->answers[m], ref.file[n].name) == 0) {
+		ref_list[n] = 1;
+		if (strcmp(name, ref.file[n].name) == 0) {
 		    got_file = 1;
 		    ref_list[n] = -1;
 		    break;
@@ -212,7 +271,9 @@ int main(int argc, char *argv[])
     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);
+    exec_rectify(order, extension, interpol->answer);
+
+    G_done_msg(" ");
 
     exit(EXIT_SUCCESS);
 }
@@ -222,12 +283,33 @@ void err_exit(char *file, char *grp)
 {
     int n;
 
-    fprintf(stderr,
-	    _("Input raster map <%s> does not exist in group <%s>.\n Try:\n"),
+    G_warning(_("Input raster map <%s> does not exist in group <%s>."),
 	    file, grp);
+    G_message(_("Try:"));
 
     for (n = 0; n < ref.nfiles; n++)
-	fprintf(stderr, "%s\n", ref.file[n].name);
+	G_message("%s", ref.file[n].name);
 
     G_fatal_error(_("Exit!"));
 }
+
+static char *make_ipol_list(void)
+{
+    int size = 0;
+    int i;
+    char *buf;
+
+    for (i = 0; menu[i].name; i++)
+	size += strlen(menu[i].name) + 1;
+
+    buf = G_malloc(size);
+    *buf = '\0';
+
+    for (i = 0; menu[i].name; i++) {
+	if (i)
+	    strcat(buf, ",");
+	strcat(buf, menu[i].name);
+    }
+
+    return buf;
+}

imagery/i.rectify/matrix.c

@@ -1,100 +0,0 @@
-
-/**************************************************************
- * compute_georef_matrix (win1, win2)
- *
- */
-#include <stdlib.h>
-#include "global.h"
-#include "crs.h"		/* CRS HEADER FILE */
-static int cmp(const void *, const void *);
-
-int compute_georef_matrix(struct Cell_head *win1,
-			  struct Cell_head *win2, int order)
-{
-    ROWCOL *rmap, *cmap, rr, cc;
-    int nrow1, nrow2;
-    int ncol1, ncol2;
-    double n1, w1, ns_res1, ew_res1;
-    double n2, e2, ns_res2, ew_res2;
-    double nx, ex;
-
-    /*    double NX, EX; DELETED WITH CRS MODIFICATIONS */
-    int row, col;
-    int min, max;
-
-    ns_res1 = win1->ns_res;
-    ew_res1 = win1->ew_res;
-    nrow1 = win1->rows;
-    ncol1 = win1->cols;
-    n1 = win1->north;
-    w1 = win1->west;
-
-    ns_res2 = win2->ns_res;
-    ew_res2 = win2->ew_res;
-    nrow2 = win2->rows;
-    ncol2 = win2->cols;
-    n2 = win2->north;
-    e2 = win2->west;
-    matrix_rows = nrow2;
-    matrix_cols = ncol2;
-
-    /* georef equation is
-     * ex = E21a + E21b * e2 + E21c * n2
-     * nx = N21a + N21b * e2 + N21c * n2
-     *
-     * compute common code (for northing) outside east loop
-     */
-    for (n2 = win2->north, row = 0; row < nrow2; row++, n2 -= ns_res2) {
-	rmap = row_map[row];
-	cmap = col_map[row];
-	min = max = -1;
-
-	/* compute common code */
-	/*  DELETED WITH CRS MODIFICATIONS
-	   EX = E21a + E21c * n2;
-	   NX = N21a + N21c * n2;
-	 */
-
-	for (e2 = win2->west, col = 0; col < ncol2; col++, e2 += ew_res2) {
-	    /* georef e2,n2 */
-
-	    CRS_georef(e2, n2, &ex, &nx, E21, N21, order);	/* BACKWARDS TRANSFORMATION */
-
-	    /* commented out by CRS
-	       ex = EX + E21b * e2;
-	       nx = NX + N21b * e2;
-	     */
-
-	    rr = (n1 - nx) / ns_res1;
-	    if (rr < 0 || rr >= nrow1)
-		rr = -1;
-	    else if (min < 0)
-		min = max = rr;
-	    else if (rr < min)
-		min = rr;
-	    else if (rr > max)
-		max = rr;
-	    *rmap++ = rr;
-
-	    cc = (ex - w1) / ew_res1;
-	    if (cc < 0 || cc >= ncol1)
-		cc = -1;
-	    *cmap++ = cc;
-	}
-	row_min[row] = min;
-	row_max[row] = max;
-	row_left[row] = 0;
-	row_right[row] = matrix_cols - 1;
-	row_idx[row] = row;
-    }
-    qsort(row_idx, nrow2, sizeof(IDX), cmp);
-
-    return 0;
-}
-
-static int cmp(const void *aa, const void *bb)
-{
-    const IDX *a = aa, *b = bb;
-
-    return (int)(row_min[*a] - row_min[*b]);
-}

imagery/i.rectify/nearest.c

@@ -0,0 +1,41 @@
+
+/*
+ *      nearest.c - returns the nearest neighbor to a given
+ *                  x,y position
+ */
+
+#include <math.h>
+#include <grass/gis.h>
+#include <grass/raster.h>
+#include "global.h"
+
+void p_nearest(struct cache *ibuffer,	/* input buffer                  */
+	       void *obufptr,	/* ptr in output buffer          */
+	       int cell_type,	/* raster map type of obufptr    */
+	       double *row_idx,	/* row index in input matrix     */
+	       double *col_idx,	/* column index in input matrix  */
+	       struct Cell_head *cellhd	/* cell header of input layer    */
+    )
+{
+    int row, col;		/* row/col of nearest neighbor   */
+    DCELL *cellp;
+
+    /* cut indices to integer */
+    row = (int)floor(*row_idx);
+    col = (int)floor(*col_idx);
+
+    /* check for out of bounds - if out of bounds set NULL value     */
+    if (row < 0 || row >= cellhd->rows || col < 0 || col >= cellhd->cols) {
+	Rast_set_null_value(obufptr, 1, cell_type);
+	return;
+    }
+
+    cellp = CPTR(ibuffer, row, col);
+
+    if (Rast_is_d_null_value(cellp)) {
+	Rast_set_null_value(obufptr, 1, cell_type);
+	return;
+    }
+
+    Rast_set_d_value(obufptr, *cellp, cell_type);
+}

imagery/i.rectify/perform.c

@@ -1,99 +0,0 @@
-
-#include <grass/raster.h>
-
-#include "global.h"
-
-static ROWCOL *rmap, *cmap, left, right;
-static int do_cell(int, void *, void *);
-
-int rast_size;
-
-int perform_georef(int infd, void *rast)
-{
-    int row;
-    int curidx;
-    int idx;
-    int i;
-
-    rast_size = Rast_cell_size(map_type);
-
-    for (row = 0; row < matrix_rows; row++)
-	Rast_set_null_value(cell_buf[row], matrix_cols, map_type);
-
-    curidx = 0;
-    while (1) {
-	/* find first row */
-	while (curidx < matrix_rows) {
-	    idx = row_idx[curidx];
-	    row = row_min[idx];
-	    if (row >= 0)
-		break;
-	    curidx++;
-	}
-	/*
-	   fprintf (stderr, " curidx %d\n", curidx);
-	 */
-	if (curidx >= matrix_rows)
-	    break;
-	/*
-	   fprintf (stderr, "read row %d\n", row);
-	 */
-
-	Rast_get_row_nomask(infd, G_incr_void_ptr(rast, rast_size), row, map_type);
-
-	for (i = curidx; i < matrix_rows; i++) {
-	    idx = row_idx[i];
-	    if (row != row_min[idx])
-		break;
-	    /*
-	       fprintf (stderr, "  process matrix row %d\n", idx);
-	     */
-	    rmap = row_map[idx];
-	    cmap = col_map[idx];
-	    left = row_left[idx];
-	    right = row_right[idx];
-	    do_cell(row, G_incr_void_ptr(rast, rast_size), cell_buf[idx]);
-
-	    row_min[idx]++;
-	    if (row_min[idx] > row_max[idx])
-		row_min[idx] = -1;
-	    row_left[idx] = left;
-	    row_right[idx] = right;
-	}
-    }
-
-    return 0;
-}
-
-static int do_cell(int row, void *in, void *out)
-{
-    int col;
-    void *inptr, *outptr;
-
-    for (; left <= right; left++) {
-	inptr = G_incr_void_ptr(in, cmap[left] * rast_size);
-	outptr = G_incr_void_ptr(out, left * rast_size);
-	if (rmap[left] < 0)
-	    continue;
-	if (rmap[left] != row)
-	    break;
-	Rast_raster_cpy(outptr, inptr, 1, map_type);
-    }
-    for (; left <= right; right--) {
-	inptr = G_incr_void_ptr(in, cmap[right] * rast_size);
-	outptr = G_incr_void_ptr(out, right * rast_size);
-	if (rmap[right] < 0)
-	    continue;
-	if (rmap[right] != row)
-	    break;
-	Rast_raster_cpy(outptr, inptr, 1, map_type);
-    }
-    for (col = left; col <= right; col++) {
-	inptr = G_incr_void_ptr(in, cmap[col] * rast_size);
-	outptr = G_incr_void_ptr(out, col * rast_size);
-	if (rmap[col] == row)
-	    Rast_raster_cpy(outptr, inptr, 1, map_type);
-    }
-
-    return 0;
-}

imagery/i.rectify/readcell.c

@@ -0,0 +1,132 @@
+/*
+ * readcell.c - reads an entire cell layer into a buffer
+ *
+ */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <sys/types.h>
+#include <sys/stat.h>
+#include <fcntl.h>
+#include <unistd.h>
+#include <grass/gis.h>
+#include <grass/raster.h>
+#include <grass/glocale.h>
+#include "global.h"
+
+struct cache *readcell(int fdi, const char *size)
+{
+    DCELL *tmpbuf;
+    struct cache *c;
+    int nrows;
+    int ncols;
+    int row;
+    char *filename;
+    int nx, ny;
+    int nblocks;
+    int i;
+
+    nrows = Rast_input_window_rows();
+    ncols = Rast_input_window_cols();
+
+    ny = (nrows + BDIM - 1) / BDIM;
+    nx = (ncols + BDIM - 1) / BDIM;
+
+    if (size)
+	nblocks = atoi(size) * ((1 << 20) / sizeof(block));
+    else
+	nblocks = (nx + ny) * 2;	/* guess */
+
+    if (nblocks > nx * ny)
+	nblocks = nx * ny;
+
+    c = G_malloc(sizeof(struct cache));
+    c->stride = nx;
+    c->nblocks = nblocks;
+    c->grid = (block **) G_calloc(nx * ny, sizeof(block *));
+    c->blocks = (block *) G_malloc(nblocks * sizeof(block));
+    c->refs = (int *)G_malloc(nblocks * sizeof(int));
+
+    if (nblocks < nx * ny) {
+	/* Temporary file must be created in output location */
+	select_target_env();
+	filename = G_tempfile();
+	select_current_env();
+	c->fd = open(filename, O_RDWR | O_CREAT | O_EXCL, 0600);
+	if (c->fd < 0)
+	    G_fatal_error(_("Unable to open temporary file"));
+	remove(filename);
+    }
+    else
+	c->fd = -1;
+
+    G_important_message(_("Allocating memory and reading input map..."));
+    G_percent(0, nrows, 5);
+
+    for (i = 0; i < c->nblocks; i++)
+	c->refs[i] = -1;
+
+    tmpbuf = (DCELL *) G_malloc(nx * sizeof(block));
+
+    for (row = 0; row < nrows; row += BDIM) {
+	int x, y;
+
+	for (y = 0; y < BDIM; y++) {
+	    G_percent(row + y, nrows, 5);
+
+	    if (row + y >= nrows)
+		break;
+
+	    Rast_get_d_row(fdi, &tmpbuf[y * nx * BDIM], row + y);
+	}
+
+	for (x = 0; x < nx; x++)
+	    for (y = 0; y < BDIM; y++)
+		if (c->fd >= 0) {
+		    if (write
+			(c->fd, &tmpbuf[(y * nx + x) * BDIM],
+			 BDIM * sizeof(DCELL)) < 0)
+			G_fatal_error(_("Error writing segment file"));
+		}
+		else
+		    memcpy(&c->blocks[BKIDX(c, HI(row), x)][LO(y)][0],
+			   &tmpbuf[(y * nx + x) * BDIM],
+			   BDIM * sizeof(DCELL));
+    }
+
+    G_free(tmpbuf);
+
+    if (c->fd < 0)
+	for (i = 0; i < c->nblocks; i++) {
+	    c->grid[i] = &c->blocks[i];
+	    c->refs[i] = i;
+	}
+
+    return c;
+}
+
+block *get_block(struct cache * c, int idx)
+{
+    int replace = rand() % c->nblocks;
+    block *p = &c->blocks[replace];
+    int ref = c->refs[replace];
+    off_t offset = (off_t) idx * sizeof(DCELL) << L2BSIZE;
+
+    if (c->fd < 0)
+	G_fatal_error(_("Internal error: cache miss on fully-cached map"));
+
+    if (ref >= 0)
+	c->grid[ref] = NULL;
+
+    c->grid[idx] = p;
+    c->refs[replace] = idx;
+
+    if (lseek(c->fd, offset, SEEK_SET) < 0)
+	G_fatal_error(_("Error seeking on segment file"));
+
+    if (read(c->fd, p, sizeof(block)) < 0)
+	G_fatal_error(_("Error writing segment file"));
+
+    return p;
+}

imagery/i.rectify/rectify.c

@@ -1,9 +1,11 @@
 #include <unistd.h>
+#include <string.h>
 
 #include <grass/raster.h>
 #include <grass/glocale.h>
 
 #include "global.h"
+#include "crs.h"		/* CRS HEADER FILE */
 
 /* Modified to support Grass 5.0 fp format 11 april 2000
  *
@@ -11,74 +13,89 @@
  *
  */
 
-int rectify(char *name, char *mapset, char *result, int order)
+int rectify(char *name, char *mapset, char *result, int order, char *interp_method)
 {
-    struct Cell_head cellhd, win;
+    struct Cell_head cellhd;
     int ncols, nrows;
-    int row, col;
-    int infd;
-    void *rast;
+    int row;
+    double row_idx, col_idx;
+    int infd, cell_size, outfd;
+    void *trast, *tptr;
+    double n1, e1, nx, ex;
+    struct cache *ibuffer;
 
     select_current_env();
     Rast_get_cellhd(name, mapset, &cellhd);
 
-    /* open the 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
-     *
-     * also tell open that raster map will have the same format
-     * (ie number of bytes per cell) as the file being rectified
-     */
-
-    select_target_env();
-    Rast_set_cell_format(cellhd.format);
-    select_current_env();
-
     /* open the file to be rectified
      * set window to cellhd first to be able to read file exactly
      */
     Rast_set_input_window(&cellhd);
     infd = Rast_open_old(name, mapset);
     map_type = Rast_get_map_type(infd);
-    rast = (void *)G_calloc(cellhd.cols + 1, Rast_cell_size(map_type));
-    Rast_set_null_value(rast, cellhd.cols + 1, map_type);
+    cell_size = Rast_cell_size(map_type);
+
+    ibuffer = readcell(infd, seg_mb);
 
-    win = target_window;
+    Rast_close(infd);		/* (pmx) 17 april 2000 */
 
-    win.west += win.ew_res / 2;
+    G_message(_("Rectify <%s@%s> (location <%s>)"),
+	      name, mapset, G_location());
+    select_target_env();
+    G_message(_("into  <%s@%s> (location <%s>) ..."),
+	      result, G_mapset(), G_location());
+
+    nrows = target_window.rows;
     ncols = target_window.cols;
-    col = 0;
 
-    temp_fd = 0;
+    if (strcmp(interp_method, "nearest") != 0) {
+	map_type = DCELL_TYPE;
+	cell_size = Rast_cell_size(map_type);
+    }
+
+    /* open the result file into target window
+     * this open must be first since we change the window later
+     * raster maps open for writing are not affected by window changes
+     * but those open for reading are
+     */
 
-    while (ncols > 0) {
-	if ((win.cols = ncols) > NCOLS)
-	    win.cols = NCOLS;
-	win.north = target_window.north - win.ns_res / 2;
-	nrows = target_window.rows;
-	row = 0;
+    outfd = Rast_open_new(result, map_type);
+    trast = Rast_allocate_output_buf(map_type);
 
-	while (nrows > 0) {
-	    if ((win.rows = nrows) > NROWS)
-		win.rows = NROWS;
+    row = 0;
+    for (n1 = target_window.north - target_window.ns_res / 2.;
+         n1 > target_window.south; n1 -= target_window.ns_res) {
 
-	    compute_georef_matrix(&cellhd, &win, order);
-	    perform_georef(infd, rast);
-	    write_matrix(row, col);
+	G_percent(row++, nrows, 2);
 
-	    nrows -= win.rows;
-	    row += win.rows;
-	    win.north -= (win.ns_res * win.rows);
-	}
+	Rast_set_null_value(trast, ncols, map_type);
+	tptr = trast;
+	for (e1 = target_window.west + target_window.ew_res / 2.;
+	     e1 < target_window.east; e1 += target_window.ew_res) {
+
+	    /* backwards transformation of target cell center */
+	    CRS_georef(e1, n1, &ex, &nx, E21, N21, order);
+
+	    /* convert to row/column indices of source raster */
+	    row_idx = (cellhd.north - nx) / cellhd.ns_res;
+	    col_idx = (ex - cellhd.west) / cellhd.ew_res;
 
-	ncols -= win.cols;
-	col += win.cols;
-	win.west += (win.ew_res * win.cols);
-	G_percent(col, col + ncols, 1);
+	    /* resample data point */
+	    interpolate(ibuffer, tptr, map_type, &row_idx, &col_idx, &cellhd);
+
+	    tptr = G_incr_void_ptr(tptr, cell_size);
+	}
+	Rast_put_row(outfd, trast, map_type);
     }
+    G_percent(1, 1, 1);
 
-    select_target_env();
+    Rast_close(outfd);		/* (pmx) 17 april 2000 */
+    G_free(trast);
+
+    close(ibuffer->fd);
+    G_free(ibuffer);
+
+    Rast_get_cellhd(result, G_mapset(), &cellhd);
 
     if (cellhd.proj == 0) {	/* x,y imagery */
 	cellhd.proj = target_window.proj;
@@ -97,12 +114,7 @@ int rectify(char *name, char *mapset, char *result, int order)
 		  name, mapset);
     }
 
-    target_window.compressed = cellhd.compressed;
-    Rast_close(infd);		/* (pmx) 17 april 2000 */
-    write_map(result);
     select_current_env();
 
-    G_free(rast);
-
     return 1;
 }

imagery/i.rectify/report.c

@@ -8,13 +8,7 @@ int report(char *name, char *mapset, char *result,
     long seconds;
     long ncells;
 
-    select_current_env();
-    G_message(_("Rectify <%s@%s> (location <%s>)"),
-	      name, mapset, G_location());
-    select_target_env();
-    G_message(_("into  <%s@%s> (location <%s>) ... %s"),
-	      result, G_mapset(), G_location(), ok ? _("complete") : _("failed"));
-    select_current_env();
+    G_message("%s", ok ? _("complete") : _("failed"));
 
     if (!ok)
 	return 1;

imagery/i.rectify/rowcol.h

@@ -1 +0,0 @@
-#define ROWCOL int

imagery/i.rectify/write.c

@@ -1,64 +0,0 @@
-#include <sys/types.h>
-#include <sys/stat.h>
-#include <fcntl.h>
-#include <unistd.h>
-
-#include <grass/raster.h>
-#include <grass/glocale.h>
-
-#include "global.h"
-
-int write_matrix(int row, int col)
-{
-    int n;
-    off_t offset;
-
-    select_target_env();
-    if (!temp_fd) {
-	temp_name = G_tempfile();
-	temp_fd = creat(temp_name, 0660);
-    }
-    for (n = 0; n < matrix_rows; n++) {
-	offset =
-	    ((off_t) row++ * target_window.cols +
-	     col) * Rast_cell_size(map_type);
-	lseek(temp_fd, offset, SEEK_SET);
-
-	if (write(temp_fd, cell_buf[n], Rast_cell_size(map_type) * matrix_cols)
-	    != Rast_cell_size(map_type) * matrix_cols) {
-	    unlink(temp_name);
-	    G_fatal_error(_("Error while writing to temp file"));
-	}
-	/*Rast_put_c_row_random (outfd, cell_buf[n], row++, col, matrix_cols); */
-    }
-    select_current_env();
-
-    return 0;
-}
-
-int write_map(char *name)
-{
-    int fd, row;
-    void *rast;
-
-    Rast_set_output_window(&target_window);
-
-    /* working with split windows, can not use Rast_allocate_buf(map_type); */
-    rast = G_malloc(target_window.cols * Rast_cell_size(map_type));
-    close(temp_fd);
-    temp_fd = open(temp_name, 0);
-    fd = Rast_open_new(name, map_type);
-
-    for (row = 0; row < target_window.rows; row++) {
-	if (read(temp_fd, rast, target_window.cols * Rast_cell_size(map_type))
-	    != target_window.cols * Rast_cell_size(map_type))
-	    G_fatal_error(_("Error writing row %d"), row);
-	Rast_put_row(fd, rast, map_type);
-    }
-    close(temp_fd);
-    unlink(temp_name);
-    Rast_close(fd);
-    G_free(rast);
-
-    return 0;
-}