[GRASS-SVN] r43934 - grass/trunk/imagery/i.rectify
svn_grass at osgeo.org
svn_grass at osgeo.org
Sat Oct 16 10:04:55 EDT 2010
Author: mmetz
Date: 2010-10-16 07:04:55 -0700 (Sat, 16 Oct 2010)
New Revision: 43934
Added:
grass/trunk/imagery/i.rectify/bilinear.c
grass/trunk/imagery/i.rectify/bilinear_f.c
grass/trunk/imagery/i.rectify/cubic.c
grass/trunk/imagery/i.rectify/cubic_f.c
grass/trunk/imagery/i.rectify/nearest.c
grass/trunk/imagery/i.rectify/readcell.c
Removed:
grass/trunk/imagery/i.rectify/matrix.c
grass/trunk/imagery/i.rectify/perform.c
grass/trunk/imagery/i.rectify/rowcol.h
grass/trunk/imagery/i.rectify/write.c
Modified:
grass/trunk/imagery/i.rectify/exec.c
grass/trunk/imagery/i.rectify/global.h
grass/trunk/imagery/i.rectify/i.rectify.html
grass/trunk/imagery/i.rectify/main.c
grass/trunk/imagery/i.rectify/rectify.c
grass/trunk/imagery/i.rectify/report.c
Log:
i.rectify: adjust to r.proj: transform cell center coords not cell border coords, use fast cache, offer different resampling methods
Added: grass/trunk/imagery/i.rectify/bilinear.c
===================================================================
--- grass/trunk/imagery/i.rectify/bilinear.c (rev 0)
+++ grass/trunk/imagery/i.rectify/bilinear.c 2010-10-16 14:04:55 UTC (rev 43934)
@@ -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);
+}
Added: grass/trunk/imagery/i.rectify/bilinear_f.c
===================================================================
--- grass/trunk/imagery/i.rectify/bilinear_f.c (rev 0)
+++ grass/trunk/imagery/i.rectify/bilinear_f.c 2010-10-16 14:04:55 UTC (rev 43934)
@@ -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);
+}
Added: grass/trunk/imagery/i.rectify/cubic.c
===================================================================
--- grass/trunk/imagery/i.rectify/cubic.c (rev 0)
+++ grass/trunk/imagery/i.rectify/cubic.c 2010-10-16 14:04:55 UTC (rev 43934)
@@ -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);
+}
Added: grass/trunk/imagery/i.rectify/cubic_f.c
===================================================================
--- grass/trunk/imagery/i.rectify/cubic_f.c (rev 0)
+++ grass/trunk/imagery/i.rectify/cubic_f.c 2010-10-16 14:04:55 UTC (rev 43934)
@@ -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);
+ }
+}
Modified: grass/trunk/imagery/i.rectify/exec.c
===================================================================
--- grass/trunk/imagery/i.rectify/exec.c 2010-10-15 22:57:25 UTC (rev 43933)
+++ grass/trunk/imagery/i.rectify/exec.c 2010-10-16 14:04:55 UTC (rev 43934)
@@ -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;
+ Rast_set_output_window(&target_window);
- /* 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);
- }
-
/* 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 @@
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 @@
}
else
report(name, mapset, result, (long)0, (long)0, 0);
+
+ G_free(result);
}
}
- G_done_msg("");
-
return 0;
}
Modified: grass/trunk/imagery/i.rectify/global.h
===================================================================
--- grass/trunk/imagery/i.rectify/global.h 2010-10-15 22:57:25 UTC (rev 43933)
+++ grass/trunk/imagery/i.rectify/global.h 2010-10-16 14:04:55 UTC (rev 43934)
@@ -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 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 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 @@
/* 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 *);
Modified: grass/trunk/imagery/i.rectify/i.rectify.html
===================================================================
--- grass/trunk/imagery/i.rectify/i.rectify.html 2010-10-15 22:57:25 UTC (rev 43933)
+++ grass/trunk/imagery/i.rectify/i.rectify.html 2010-10-16 14:04:55 UTC (rev 43934)
@@ -7,10 +7,10 @@
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,112 +21,56 @@
<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>
+<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 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
- >
-</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 <ESC> TO CONTINUE
- (OR<Ctrl-C> 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
- >
-</PRE>
-
-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:
+The program will calculate the RMSE and check the required number of points.
-<PRE>
-Select order of transformation --> 1st Order 2nd Order 3rd Order
-</PRE>
+<h4>Linear affine transformation (1st order transformation)</h4>
-The program will immediately recalculate the RMSE and the
-number of points required.
-
-<h3>Linear affine transformation (1st order transformation)</h3>
-
<DL>
<DD> x' = ax + by +c
<DD> y' = Ax + Bt +C
@@ -134,43 +78,65 @@
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>
+<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
-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
-
<DL>
<DD>((n + 1) * (n + 2) / 2)
</DL>
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>
+<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.
-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.
-
<P>
<!--
Note: In interactive mode it is possible to define a new file name
@@ -180,19 +146,12 @@
<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 @@
<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>
Modified: grass/trunk/imagery/i.rectify/main.c
===================================================================
--- grass/trunk/imagery/i.rectify/main.c 2010-10-15 22:57:25 UTC (rev 43933)
+++ grass/trunk/imagery/i.rectify/main.c 2010-10-16 14:04:55 UTC (rev 43934)
@@ -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 @@
char **new_name;
struct Ref ref;
+func interpolate;
+
/* georef coefficients */
double E12[10], N12[10];
@@ -60,19 +55,40 @@
*/
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 @@
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 @@
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 @@
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 @@
}
}
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,8 +271,10 @@
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 @@
{
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;
+}
Deleted: grass/trunk/imagery/i.rectify/matrix.c
===================================================================
--- grass/trunk/imagery/i.rectify/matrix.c 2010-10-15 22:57:25 UTC (rev 43933)
+++ grass/trunk/imagery/i.rectify/matrix.c 2010-10-16 14:04:55 UTC (rev 43934)
@@ -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]);
-}
Added: grass/trunk/imagery/i.rectify/nearest.c
===================================================================
--- grass/trunk/imagery/i.rectify/nearest.c (rev 0)
+++ grass/trunk/imagery/i.rectify/nearest.c 2010-10-16 14:04:55 UTC (rev 43934)
@@ -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);
+}
Deleted: grass/trunk/imagery/i.rectify/perform.c
===================================================================
--- grass/trunk/imagery/i.rectify/perform.c 2010-10-15 22:57:25 UTC (rev 43933)
+++ grass/trunk/imagery/i.rectify/perform.c 2010-10-16 14:04:55 UTC (rev 43934)
@@ -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;
-}
Added: grass/trunk/imagery/i.rectify/readcell.c
===================================================================
--- grass/trunk/imagery/i.rectify/readcell.c (rev 0)
+++ grass/trunk/imagery/i.rectify/readcell.c 2010-10-16 14:04:55 UTC (rev 43934)
@@ -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;
+}
Modified: grass/trunk/imagery/i.rectify/rectify.c
===================================================================
--- grass/trunk/imagery/i.rectify/rectify.c 2010-10-15 22:57:25 UTC (rev 43933)
+++ grass/trunk/imagery/i.rectify/rectify.c 2010-10-16 14:04:55 UTC (rev 43934)
@@ -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,75 +13,90 @@
*
*/
-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);
- win = target_window;
+ ibuffer = readcell(infd, seg_mb);
- win.west += win.ew_res / 2;
+ Rast_close(infd); /* (pmx) 17 april 2000 */
+
+ 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);
+ }
- 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;
+ /* 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 (nrows > 0) {
- if ((win.rows = nrows) > NROWS)
- win.rows = NROWS;
+ outfd = Rast_open_new(result, map_type);
+ trast = Rast_allocate_output_buf(map_type);
- compute_georef_matrix(&cellhd, &win, order);
- perform_georef(infd, rast);
- write_matrix(row, col);
+ row = 0;
+ for (n1 = target_window.north - target_window.ns_res / 2.;
+ n1 > target_window.south; n1 -= target_window.ns_res) {
- nrows -= win.rows;
- row += win.rows;
- win.north -= (win.ns_res * win.rows);
+ G_percent(row++, nrows, 2);
+
+ 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;
+
+ /* resample data point */
+ interpolate(ibuffer, tptr, map_type, &row_idx, &col_idx, &cellhd);
+
+ tptr = G_incr_void_ptr(tptr, cell_size);
}
-
- ncols -= win.cols;
- col += win.cols;
- win.west += (win.ew_res * win.cols);
- G_percent(col, col + ncols, 1);
+ 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;
cellhd.zone = target_window.zone;
@@ -97,12 +114,7 @@
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;
}
Modified: grass/trunk/imagery/i.rectify/report.c
===================================================================
--- grass/trunk/imagery/i.rectify/report.c 2010-10-15 22:57:25 UTC (rev 43933)
+++ grass/trunk/imagery/i.rectify/report.c 2010-10-16 14:04:55 UTC (rev 43934)
@@ -8,13 +8,7 @@
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;
Deleted: grass/trunk/imagery/i.rectify/rowcol.h
===================================================================
--- grass/trunk/imagery/i.rectify/rowcol.h 2010-10-15 22:57:25 UTC (rev 43933)
+++ grass/trunk/imagery/i.rectify/rowcol.h 2010-10-16 14:04:55 UTC (rev 43934)
@@ -1 +0,0 @@
-#define ROWCOL int
Deleted: grass/trunk/imagery/i.rectify/write.c
===================================================================
--- grass/trunk/imagery/i.rectify/write.c 2010-10-15 22:57:25 UTC (rev 43933)
+++ grass/trunk/imagery/i.rectify/write.c 2010-10-16 14:04:55 UTC (rev 43934)
@@ -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;
-}
More information about the grass-commit
mailing list