[GRASS-SVN] r46222 -
grass/branches/releasebranch_6_4/imagery/i.rectify
svn_grass at osgeo.org
svn_grass at osgeo.org
Mon May 9 11:33:35 EDT 2011
Author: mmetz
Date: 2011-05-09 08:33:35 -0700 (Mon, 09 May 2011)
New Revision: 46222
Added:
grass/branches/releasebranch_6_4/imagery/i.rectify/bilinear.c
grass/branches/releasebranch_6_4/imagery/i.rectify/bilinear_f.c
grass/branches/releasebranch_6_4/imagery/i.rectify/cubic.c
grass/branches/releasebranch_6_4/imagery/i.rectify/cubic_f.c
grass/branches/releasebranch_6_4/imagery/i.rectify/nearest.c
grass/branches/releasebranch_6_4/imagery/i.rectify/readcell.c
Removed:
grass/branches/releasebranch_6_4/imagery/i.rectify/matrix.c
grass/branches/releasebranch_6_4/imagery/i.rectify/perform.c
grass/branches/releasebranch_6_4/imagery/i.rectify/rowcol.h
grass/branches/releasebranch_6_4/imagery/i.rectify/write.c
Modified:
grass/branches/releasebranch_6_4/imagery/i.rectify/description.html
grass/branches/releasebranch_6_4/imagery/i.rectify/exec.c
grass/branches/releasebranch_6_4/imagery/i.rectify/get_wind.c
grass/branches/releasebranch_6_4/imagery/i.rectify/global.h
grass/branches/releasebranch_6_4/imagery/i.rectify/main.c
grass/branches/releasebranch_6_4/imagery/i.rectify/rectify.c
grass/branches/releasebranch_6_4/imagery/i.rectify/report.c
Log:
backport i.rectify
Added: grass/branches/releasebranch_6_4/imagery/i.rectify/bilinear.c
===================================================================
--- grass/branches/releasebranch_6_4/imagery/i.rectify/bilinear.c (rev 0)
+++ grass/branches/releasebranch_6_4/imagery/i.rectify/bilinear.c 2011-05-09 15:33:35 UTC (rev 46222)
@@ -0,0 +1,59 @@
+/*
+ * 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 "global.h"
+
+void p_bilinear(struct cache *ibuffer, /* input buffer */
+ void *obufptr, /* ptr in output buffer */
+ int cell_type, /* raster map type of obufptr */
+ double *row_idx, /* row index */
+ double *col_idx, /* column index */
+ struct Cell_head *cellhd /* information of output map */
+ )
+{
+ int row; /* row indices for interp */
+ int col; /* column indices for interp */
+ int i, j;
+ DCELL t, u; /* intermediate slope */
+ DCELL result; /* result of interpolation */
+ DCELL c[2][2];
+
+ /* cut indices to integer */
+ row = (int)floor(*row_idx - 0.5);
+ col = (int)floor(*col_idx - 0.5);
+
+ /* check for out of bounds - if out of bounds set NULL value and return */
+ if (row < 0 || row + 1 >= cellhd->rows || col < 0 || col + 1 >= cellhd->cols) {
+ G_set_null_value(obufptr, 1, cell_type);
+ return;
+ }
+
+ for (i = 0; i < 2; i++)
+ for (j = 0; j < 2; j++) {
+ const DCELL *cellp = CPTR(ibuffer, row + i, col + j);
+ if (G_is_d_null_value(cellp)) {
+ G_set_null_value(obufptr, 1, cell_type);
+ return;
+ }
+ c[i][j] = *cellp;
+ }
+
+ /* do the interpolation */
+ t = *col_idx - 0.5 - col;
+ u = *row_idx - 0.5 - row;
+
+ result = G_interp_bilinear(t, u, c[0][0], c[0][1], c[1][0], c[1][1]);
+
+ G_set_raster_value_d(obufptr, result, cell_type);
+}
Added: grass/branches/releasebranch_6_4/imagery/i.rectify/bilinear_f.c
===================================================================
--- grass/branches/releasebranch_6_4/imagery/i.rectify/bilinear_f.c (rev 0)
+++ grass/branches/releasebranch_6_4/imagery/i.rectify/bilinear_f.c 2011-05-09 15:33:35 UTC (rev 46222)
@@ -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 <math.h>
+#include "global.h"
+
+void p_bilinear_f(struct cache *ibuffer, /* input buffer */
+ void *obufptr, /* ptr in output buffer */
+ int cell_type, /* raster map type of obufptr */
+ double *row_idx, /* row index */
+ double *col_idx, /* column index */
+ struct Cell_head *cellhd /* cell header of input layer */
+ )
+{
+ /* start nearest neighbor to do some basic tests */
+ int row, col; /* row/col of nearest neighbor */
+ DCELL *cellp, cell;
+
+ /* cut indices to integer */
+ row = (int)floor(*row_idx);
+ col = (int)floor(*col_idx);
+
+ /* check for out of bounds - if out of bounds set NULL value */
+ if (row < 0 || row >= cellhd->rows || col < 0 || col >= cellhd->cols) {
+ G_set_null_value(obufptr, 1, cell_type);
+ return;
+ }
+
+ cellp = CPTR(ibuffer, row, col);
+ /* if nearest is null, all the other interps will be null */
+ if (G_is_d_null_value(cellp)) {
+ G_set_null_value(obufptr, 1, cell_type);
+ return;
+ }
+ cell = *cellp;
+
+ p_bilinear(ibuffer, obufptr, cell_type, row_idx, col_idx, cellhd);
+ /* fallback to nearest if bilinear is null */
+ if (G_is_d_null_value(obufptr))
+ G_set_raster_value_d(obufptr, cell, cell_type);
+}
Added: grass/branches/releasebranch_6_4/imagery/i.rectify/cubic.c
===================================================================
--- grass/branches/releasebranch_6_4/imagery/i.rectify/cubic.c (rev 0)
+++ grass/branches/releasebranch_6_4/imagery/i.rectify/cubic.c 2011-05-09 15:33:35 UTC (rev 46222)
@@ -0,0 +1,67 @@
+/*
+ * 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 <math.h>
+#include "global.h"
+
+void p_cubic(struct cache *ibuffer, /* input buffer */
+ void *obufptr, /* ptr in output buffer */
+ int cell_type, /* raster map type of obufptr */
+ double *row_idx, /* row index (decimal) */
+ double *col_idx, /* column index (decimal) */
+ struct Cell_head *cellhd /* information of output map */
+ )
+{
+ int row; /* row indices for interp */
+ int col; /* column indices for interp */
+ int i, j;
+ DCELL t, u; /* intermediate slope */
+ DCELL result; /* result of interpolation */
+ DCELL val[4]; /* buffer for temporary values */
+ DCELL cell[4][4];
+
+ /* cut indices to integer */
+ row = (int)floor(*row_idx - 0.5);
+ col = (int)floor(*col_idx - 0.5);
+
+ /* check for out of bounds of map - if out of bounds set NULL value */
+ if (row - 1 < 0 || row + 2 >= cellhd->rows ||
+ col - 1 < 0 || col + 2 >= cellhd->cols) {
+ G_set_null_value(obufptr, 1, cell_type);
+ return;
+ }
+
+ for (i = 0; i < 4; i++)
+ for (j = 0; j < 4; j++) {
+ const DCELL *cellp = CPTR(ibuffer, row - 1 + i, col - 1 + j);
+ if (G_is_d_null_value(cellp)) {
+ G_set_null_value(obufptr, 1, cell_type);
+ return;
+ }
+ cell[i][j] = *cellp;
+ }
+
+ /* do the interpolation */
+ t = *col_idx - 0.5 - col;
+ u = *row_idx - 0.5 - row;
+
+ for (i = 0; i < 4; i++) {
+ val[i] = G_interp_cubic(t, cell[i][0], cell[i][1], cell[i][2], cell[i][3]);
+ }
+
+ result = G_interp_cubic(u, val[0], val[1], val[2], val[3]);
+
+ G_set_raster_value_d(obufptr, result, cell_type);
+}
Added: grass/branches/releasebranch_6_4/imagery/i.rectify/cubic_f.c
===================================================================
--- grass/branches/releasebranch_6_4/imagery/i.rectify/cubic_f.c (rev 0)
+++ grass/branches/releasebranch_6_4/imagery/i.rectify/cubic_f.c 2011-05-09 15:33:35 UTC (rev 46222)
@@ -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 <math.h>
+#include "global.h"
+
+void p_cubic_f(struct cache *ibuffer, /* input buffer */
+ void *obufptr, /* ptr in output buffer */
+ int cell_type, /* raster map type of obufptr */
+ double *row_idx, /* row index */
+ double *col_idx, /* column index */
+ struct Cell_head *cellhd /* cell header of input layer */
+ )
+{
+ /* start nearest neighbor to do some basic tests */
+ int row, col; /* row/col of nearest neighbor */
+ DCELL *cellp, cell;
+
+ /* cut indices to integer */
+ row = (int)floor(*row_idx);
+ col = (int)floor(*col_idx);
+
+ /* check for out of bounds - if out of bounds set NULL value */
+ if (row < 0 || row >= cellhd->rows || col < 0 || col >= cellhd->cols) {
+ G_set_null_value(obufptr, 1, cell_type);
+ return;
+ }
+
+ cellp = CPTR(ibuffer, row, col);
+ /* if nearest is null, all the other interps will be null */
+ if (G_is_d_null_value(cellp)) {
+ G_set_null_value(obufptr, 1, cell_type);
+ return;
+ }
+ cell = *cellp;
+
+ p_cubic(ibuffer, obufptr, cell_type, row_idx, col_idx, cellhd);
+ /* fallback to bilinear if cubic is null */
+ if (G_is_d_null_value(obufptr)) {
+ p_bilinear(ibuffer, obufptr, cell_type, row_idx, col_idx, cellhd);
+ /* fallback to nearest if bilinear is null */
+ if (G_is_d_null_value(obufptr))
+ G_set_raster_value_d(obufptr, cell, cell_type);
+ }
+}
Modified: grass/branches/releasebranch_6_4/imagery/i.rectify/description.html
===================================================================
--- grass/branches/releasebranch_6_4/imagery/i.rectify/description.html 2011-05-09 14:14:19 UTC (rev 46221)
+++ grass/branches/releasebranch_6_4/imagery/i.rectify/description.html 2011-05-09 15:33:35 UTC (rev 46222)
@@ -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,12 +21,12 @@
<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>.
@@ -164,6 +164,39 @@
menu bar. The polynomial equations are performed using a
modified Gaussian elimination method.
+<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. The <em>method=lanczos</em> method determines the new value of
+the cell based on a weighted distance average of the 25 surrounding cells in
+the input map.
+<p>
+The bilinear, cubic and lanczos interpolation methods are most appropriate for
+continuous data and cause some smoothing. These options should not be used
+with categorical data, since the cell values will be altered.
+<p>
+In the bilinear, cubic and lanczos 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, cubic_f and lanczos_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 lanczos to 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.
+
<H3>Program Execution</H3>
Note: The rectified image or rectified raster maps will be
Modified: grass/branches/releasebranch_6_4/imagery/i.rectify/exec.c
===================================================================
--- grass/branches/releasebranch_6_4/imagery/i.rectify/exec.c 2011-05-09 14:14:19 UTC (rev 46221)
+++ grass/branches/releasebranch_6_4/imagery/i.rectify/exec.c 2011-05-09 15:33:35 UTC (rev 46222)
@@ -16,33 +16,25 @@
#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;
+ long start_time, rectify_time;
+ G_message("-----------------------------------------------");
- /* allocate the output cell matrix */
- cell_buf = (void **)G_calloc(NROWS, sizeof(void *));
- n = NCOLS * G_raster_size(map_type);
- for (i = 0; i < NROWS; i++) {
- cell_buf[i] = (void *)G_malloc(n);
- G_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) {
- /* continue; */
+ if (ref_list[n]) {
name = ref.file[n].name;
mapset = ref.file[n].mapset;
@@ -51,8 +43,6 @@
G_malloc(strlen(ref.file[n].name) + strlen(extension) + 1);
strcpy(result, ref.file[n].name);
strcat(result, extension);
- G_message(_("Rectified input raster map <%s> will be saved as <%s>"),
- name, result);
select_current_env();
@@ -60,20 +50,14 @@
colr_ok = G_read_colors(name, mapset, &colr) > 0;
/* Initialze History */
- type = "raster";
- G_short_history(name, type, &hist);
+ if (G_read_history(name, mapset, &hist) < 0)
+ G_short_history(result, 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
- G_put_cellhd (result,&target_window);
- */
if (cats_ok) {
G_write_cats(result, &cats);
G_free_cats(&cats);
@@ -83,26 +67,20 @@
G_free_colors(&colr);
}
- /* Write out History Structure History */
- sprintf(hist.title, "%s", result);
- sprintf(hist.datsrc_1, "%s", name);
- sprintf(hist.edhist[0], "Created from: i.rectify");
- sprintf(hist.edhist[1], "Transformation order = %d", order);
- hist.edlinecnt = 2;
+ /* Write out History */
+ G_command_history(&hist);
G_write_history(result, &hist);
select_current_env();
time(&rectify_time);
- compress_time = rectify_time;
- report(name, mapset, result, rectify_time - start_time,
- compress_time - rectify_time, 1);
+ report(rectify_time - start_time, 1);
}
else
- report(name, mapset, result, (long)0, (long)0, 0);
+ report((long)0, 0);
+
+ G_free(result);
}
}
- G_done_msg("");
-
return 0;
}
Modified: grass/branches/releasebranch_6_4/imagery/i.rectify/get_wind.c
===================================================================
--- grass/branches/releasebranch_6_4/imagery/i.rectify/get_wind.c 2011-05-09 14:14:19 UTC (rev 46221)
+++ grass/branches/releasebranch_6_4/imagery/i.rectify/get_wind.c 2011-05-09 15:33:35 UTC (rev 46222)
@@ -1,3 +1,4 @@
+#include <math.h>
#include <grass/glocale.h>
#include "global.h"
#include "crs.h" /* CRS HEADER FILE */
@@ -2,11 +3,19 @@
-int georef_window(struct Cell_head *w1, struct Cell_head *w2, int order)
+int georef_window(struct Cell_head *w1, struct Cell_head *w2, int order, double res)
{
- double n, e;
+ double n, e, ad;
+ struct _corner {
+ double n, e;
+ } nw, ne, se, sw;
+ /* extends */
CRS_georef(w1->west, w1->north, &e, &n, E12, N12, order);
w2->north = w2->south = n;
w2->west = w2->east = e;
+ nw.n = n;
+ nw.e = e;
CRS_georef(w1->east, w1->north, &e, &n, E12, N12, order);
+ ne.n = n;
+ ne.e = e;
if (n > w2->north)
@@ -21,6 +30,8 @@
w2->west = e;
CRS_georef(w1->west, w1->south, &e, &n, E12, N12, order);
+ sw.n = n;
+ sw.e = e;
if (n > w2->north)
w2->north = n;
if (n < w2->south)
@@ -31,6 +42,8 @@
w2->west = e;
CRS_georef(w1->east, w1->south, &e, &n, E12, N12, order);
+ se.n = n;
+ se.e = e;
if (n > w2->north)
w2->north = n;
if (n < w2->south)
@@ -40,9 +53,63 @@
if (e < w2->west)
w2->west = e;
- w2->ns_res = (w2->north - w2->south) / w1->rows;
- w2->ew_res = (w2->east - w2->west) / w1->cols;
+ /* resolution */
+ if (res > 0)
+ w2->ew_res = w2->ns_res = res;
+ else {
+ /* this results in ugly res values, and ns_res != ew_res */
+ /* and is no good for rotation */
+ /*
+ w2->ns_res = (w2->north - w2->south) / w1->rows;
+ w2->ew_res = (w2->east - w2->west) / w1->cols;
+ */
+ /* alternative: account for rotation and order > 1 */
+
+ /* N-S extends along western and eastern edge */
+ w2->ns_res = (sqrt((nw.n - sw.n) * (nw.n - sw.n) +
+ (nw.e - sw.e) * (nw.e - sw.e)) +
+ sqrt((ne.n - se.n) * (ne.n - se.n) +
+ (ne.e - se.e) * (ne.e - se.e))) / (2.0 * w1->rows);
+
+ /* E-W extends along northern and southern edge */
+ w2->ew_res = (sqrt((nw.n - ne.n) * (nw.n - ne.n) +
+ (nw.e - ne.e) * (nw.e - ne.e)) +
+ sqrt((sw.n - se.n) * (sw.n - se.n) +
+ (sw.e - se.e) * (sw.e - se.e))) / (2.0 * w1->cols);
+
+ /* make ew_res = ns_res */
+ w2->ns_res = (w2->ns_res + w2->ew_res) / 2.0;
+ w2->ew_res = w2->ns_res;
+
+ /* nice round values */
+ if (w2->ns_res > 1) {
+ if (w2->ns_res < 10) {
+ /* round to first decimal */
+ w2->ns_res = (int)(w2->ns_res * 10 + 0.5) / 10.0;
+ w2->ew_res = w2->ns_res;
+ }
+ else {
+ /* round to integer */
+ w2->ns_res = (int)(w2->ns_res + 0.5);
+ w2->ew_res = w2->ns_res;
+ }
+ }
+ }
+
+ /* adjust extends */
+ ad = w2->north > 0 ? 0.5 : -0.5;
+ w2->north = (int) (ceil(w2->north / w2->ns_res) + ad) * w2->ns_res;
+ ad = w2->south > 0 ? 0.5 : -0.5;
+ w2->south = (int) (floor(w2->south / w2->ns_res) + ad) * w2->ns_res;
+ ad = w2->east > 0 ? 0.5 : -0.5;
+ w2->east = (int) (ceil(w2->east / w2->ew_res) + ad) * w2->ew_res;
+ ad = w2->west > 0 ? 0.5 : -0.5;
+ w2->west = (int) (floor(w2->west / w2->ew_res) + ad) * w2->ew_res;
+
+ w2->rows = (w2->north - w2->south + w2->ns_res / 2.0) / w2->ns_res;
+ w2->cols = (w2->east - w2->west + w2->ew_res / 2.0) / w2->ew_res;
+
G_message(_("Region N=%f S=%f E=%f W=%f"), w2->north, w2->south,
w2->east, w2->west);
G_message(_("Resolution EW=%f NS=%f"), w2->ew_res, w2->ns_res);
Modified: grass/branches/releasebranch_6_4/imagery/i.rectify/global.h
===================================================================
--- grass/branches/releasebranch_6_4/imagery/i.rectify/global.h 2011-05-09 14:14:19 UTC (rev 46221)
+++ grass/branches/releasebranch_6_4/imagery/i.rectify/global.h 2011-05-09 15:33:35 UTC (rev 46222)
@@ -5,32 +5,43 @@
* (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 void **cell_buf;
-extern int temp_fd;
+extern char *seg_mb;
extern RASTER_MAP_TYPE map_type;
-extern char *temp_name;
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,27 +62,48 @@
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);
+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);
+int report(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 *);
+/* lanczos.c */
+extern void p_lanczos(struct cache *, void *, int, double *, double *,
+ struct Cell_head *);
+extern void p_lanczos_f(struct cache *, void *, int, double *, double *,
+ struct Cell_head *);
Modified: grass/branches/releasebranch_6_4/imagery/i.rectify/main.c
===================================================================
--- grass/branches/releasebranch_6_4/imagery/i.rectify/main.c 2011-05-09 14:14:19 UTC (rev 46221)
+++ grass/branches/releasebranch_6_4/imagery/i.rectify/main.c 2011-05-09 15:33:35 UTC (rev 46222)
@@ -10,11 +10,12 @@
* Markus Neteler <neteler itc.it>,
* Bernhard Reiter <bernhard intevation.de>,
* Glynn Clements <glynn gclements.plus.com>,
- * Hamish Bowman <hamish_nospam 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)
- * COPYRIGHT: (C) 2002-2006 by the GRASS Development Team
+ * COPYRIGHT: (C) 2002-2011 by the GRASS Development Team
*
* This program is free software under the GNU General Public
* License (>=v2). Read the file COPYING that comes with GRASS
@@ -24,28 +25,18 @@
#include <stdlib.h>
#include <string.h>
+#include "global.h"
#include <grass/glocale.h>
-#include "global.h"
#include "crs.h"
+char *seg_mb;
-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;
int *ref_list;
-char **new_name;
struct Ref ref;
+func interpolate;
+
/* georef coefficients */
double E12[10], N12[10];
@@ -57,27 +48,42 @@
*/
struct Cell_head target_window;
-#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 */
- int n, i, m, k;
- int got_file = 0;
+ char *ipolname; /* name of interpolation method */
+ int method;
+ int n, i, m, k = 0;
+ int got_file = 0, target_overwrite = 0;
+ char *overstr;
+ struct Cell_head cellhd;
- struct Option *grp, *val, *ifile, *ext;
+ 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;
- struct Cell_head cellhd;
-
- setbuf(stdout, NULL);
- setbuf(stderr, NULL);
-
G_gisinit(argv[0]);
module = G_define_module();
@@ -105,6 +111,30 @@
val->required = YES;
val->description = _("Rectification polynom order (1-3)");
+ tres = G_define_option();
+ tres->key = "res";
+ tres->type = TYPE_DOUBLE;
+ tres->required = NO;
+ tres->description = _("Target resolution (ignored if -c flag used)");
+
+ mem = G_define_option();
+ mem->key = "memory";
+ mem->type = TYPE_DOUBLE;
+ mem->key_desc = "memory in MB";
+ mem->required = NO;
+ mem->answer = "300";
+ mem->description = _("Amount of memory to use in MB");
+
+ ipolname = make_ipol_list();
+
+ interpol = G_define_option();
+ interpol->key = "method";
+ interpol->type = TYPE_STRING;
+ interpol->required = NO;
+ interpol->answer = "nearest";
+ interpol->options = ipolname;
+ interpol->description = _("Interpolation method to use");
+
c = G_define_flag();
c->key = 'c';
c->description =
@@ -114,19 +144,34 @@
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 */
- /* Find out how files on command line */
+ /* Find out how many files on command line */
if (!a->answer) {
for (m = 0; ifile->answers[m]; m++) {
k = m;
@@ -148,26 +193,47 @@
exit(EXIT_SUCCESS);
}
- for (i = 0; i < NFILES; i++)
- result[i][0] = 0;
-
ref_list = (int *)G_malloc(ref.nfiles * sizeof(int));
- new_name = (char **)G_malloc(ref.nfiles * sizeof(char *));
if (a->answer) {
for (n = 0; n < ref.nfiles; n++) {
- ref_list[n] = -1;
+ ref_list[n] = 1;
}
}
else {
+ char xname[GNAME_MAX], xmapset[GMAPSET_MAX], *name, *mapset;
+
+ for (n = 0; n < ref.nfiles; n++)
+ ref_list[n] = 0;
+
for (m = 0; m < k; m++) {
got_file = 0;
+ if (G__name_is_fully_qualified(ifile->answers[m], xname, xmapset)) {
+ name = xname;
+ mapset = xmapset;
+ }
+ else {
+ name = ifile->answers[m];
+ mapset = NULL;
+ }
+
+ got_file = 0;
for (n = 0; n < ref.nfiles; n++) {
- if (strcmp(ifile->answers[m], ref.file[n].name) == 0) {
- got_file = 1;
- ref_list[n] = -1;
- break;
+ if (mapset) {
+ if (strcmp(name, ref.file[n].name) == 0 &&
+ strcmp(mapset, ref.file[n].mapset) == 0) {
+ got_file = 1;
+ ref_list[n] = 1;
+ break;
+ }
}
+ else {
+ if (strcmp(name, ref.file[n].name) == 0) {
+ got_file = 1;
+ ref_list[n] = 1;
+ break;
+ }
+ }
}
if (got_file == 0)
err_exit(ifile->answers[m], group);
@@ -180,12 +246,47 @@
/* get the target */
get_target(group);
+ /* Check the GRASS_OVERWRITE environment variable */
+ if ((overstr = getenv("GRASS_OVERWRITE"))) /* OK ? */
+ target_overwrite = atoi(overstr);
- if (c->answer) {
- /* Use current Region */
- G_get_window(&target_window);
+ if (!target_overwrite) {
+ /* check if output exists in target location/mapset */
+ char result[GNAME_MAX];
+
+ select_target_env();
+ for (i = 0; i < ref.nfiles; i++) {
+ if (!ref_list[i])
+ continue;
+
+ strcpy(result, ref.file[i].name);
+ strcat(result, extension);
+
+ if (G_legal_filename(result) < 0)
+ G_fatal_error(_("Extension <%s> is illegal"), extension);
+
+ if (G_find_cell(result, G_mapset())) {
+ G_warning(_("The following raster map already exists in"));
+ G_warning(_("target LOCATION %s, MAPSET %s:"),
+ G_location(), G_mapset());
+ G_warning("<%s>", result);
+ G_fatal_error(_("Orthorectification cancelled."));
+ }
+ }
+
+ select_current_env();
}
- else {
+ else
+ G_debug(1, "Overwriting OK");
+
+ /* do not use current region in target location */
+ if (!c->answer) {
+ double res = -1;
+
+ if (tres->answer) {
+ if (!((res = atof(tres->answer)) > 0))
+ G_warning(_("Target resolution must be > 0, ignored"));
+ }
/* Calculate smallest region */
if (a->answer) {
if (G_get_cellhd(ref.file[0].name, ref.file[0].mapset, &cellhd) <
@@ -199,14 +300,16 @@
G_fatal_error(_("Unable to read header of raster map <%s>"),
ifile->answers[0]);
}
- georef_window(&cellhd, &target_window, order);
+ georef_window(&cellhd, &target_window, order, res);
}
G_verbose_message(_("Using region: N=%f S=%f, E=%f W=%f"), target_window.north,
target_window.south, target_window.east, target_window.west);
- exec_rectify(order, extension);
+ exec_rectify(order, extension, interpol->answer);
+ G_done_msg(" ");
+
exit(EXIT_SUCCESS);
}
@@ -215,12 +318,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/branches/releasebranch_6_4/imagery/i.rectify/matrix.c
===================================================================
--- grass/branches/releasebranch_6_4/imagery/i.rectify/matrix.c 2011-05-09 14:14:19 UTC (rev 46221)
+++ grass/branches/releasebranch_6_4/imagery/i.rectify/matrix.c 2011-05-09 15:33:35 UTC (rev 46222)
@@ -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/branches/releasebranch_6_4/imagery/i.rectify/nearest.c
===================================================================
--- grass/branches/releasebranch_6_4/imagery/i.rectify/nearest.c (rev 0)
+++ grass/branches/releasebranch_6_4/imagery/i.rectify/nearest.c 2011-05-09 15:33:35 UTC (rev 46222)
@@ -0,0 +1,41 @@
+
+/*
+ * nearest.c - returns the nearest neighbor to a given
+ * x,y position
+ */
+
+#include <math.h>
+#include <grass/gis.h>
+#include "global.h"
+
+void p_nearest(struct cache *ibuffer, /* input buffer */
+ void *obufptr, /* ptr in output buffer */
+ int cell_type, /* raster map type of obufptr */
+ double *row_idx, /* row index in input matrix */
+ double *col_idx, /* column index in input matrix */
+ struct Cell_head *cellhd /* cell header of input layer */
+ )
+{
+ int row, col; /* row/col of nearest neighbor */
+ DCELL *cellp;
+
+ /* cut indices to integer and get nearest cell */
+ /* the row_idx, col_idx correction for bilinear/bicubic does not apply here */
+ row = (int)floor(*row_idx);
+ col = (int)floor(*col_idx);
+
+ /* check for out of bounds - if out of bounds set NULL value */
+ if (row < 0 || row >= cellhd->rows || col < 0 || col >= cellhd->cols) {
+ G_set_null_value(obufptr, 1, cell_type);
+ return;
+ }
+
+ cellp = CPTR(ibuffer, row, col);
+
+ if (G_is_d_null_value(cellp)) {
+ G_set_null_value(obufptr, 1, cell_type);
+ return;
+ }
+
+ G_set_raster_value_d(obufptr, *cellp, cell_type);
+}
Deleted: grass/branches/releasebranch_6_4/imagery/i.rectify/perform.c
===================================================================
--- grass/branches/releasebranch_6_4/imagery/i.rectify/perform.c 2011-05-09 14:14:19 UTC (rev 46221)
+++ grass/branches/releasebranch_6_4/imagery/i.rectify/perform.c 2011-05-09 15:33:35 UTC (rev 46222)
@@ -1,98 +0,0 @@
-#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 = G_raster_size(map_type);
-
- for (row = 0; row < matrix_rows; row++)
- G_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);
- */
-
- if (G_get_raster_row_nomask
- (infd, G_incr_void_ptr(rast, rast_size), row, map_type) < 0)
- return 0;
-
- 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;
- G_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;
- G_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)
- G_raster_cpy(outptr, inptr, 1, map_type);
- }
-
- return 0;
-}
Added: grass/branches/releasebranch_6_4/imagery/i.rectify/readcell.c
===================================================================
--- grass/branches/releasebranch_6_4/imagery/i.rectify/readcell.c (rev 0)
+++ grass/branches/releasebranch_6_4/imagery/i.rectify/readcell.c 2011-05-09 15:33:35 UTC (rev 46222)
@@ -0,0 +1,129 @@
+/*
+ * 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/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 = G_window_rows();
+ ncols = G_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 input location */
+ filename = G_tempfile();
+ c->fd = open(filename, O_RDWR | O_CREAT | O_EXCL, 0600);
+ if (c->fd < 0)
+ G_fatal_error(_("Unable to open temporary file"));
+ remove(filename);
+ }
+ else
+ c->fd = -1;
+
+ G_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;
+
+ G_get_d_raster_row(fdi, &tmpbuf[y * nx * BDIM], row + y);
+ }
+
+ for (x = 0; x < nx; x++)
+ for (y = 0; y < BDIM; y++)
+ if (c->fd >= 0) {
+ if (write
+ (c->fd, &tmpbuf[(y * nx + x) * BDIM],
+ BDIM * sizeof(DCELL)) < 0)
+ G_fatal_error(_("Error writing segment file"));
+ }
+ else
+ memcpy(&c->blocks[BKIDX(c, HI(row), x)][LO(y)][0],
+ &tmpbuf[(y * nx + x) * BDIM],
+ BDIM * sizeof(DCELL));
+ }
+
+ G_free(tmpbuf);
+
+ if (c->fd < 0)
+ for (i = 0; i < c->nblocks; i++) {
+ c->grid[i] = &c->blocks[i];
+ c->refs[i] = i;
+ }
+
+ 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/branches/releasebranch_6_4/imagery/i.rectify/rectify.c
===================================================================
--- grass/branches/releasebranch_6_4/imagery/i.rectify/rectify.c 2011-05-09 14:14:19 UTC (rev 46221)
+++ grass/branches/releasebranch_6_4/imagery/i.rectify/rectify.c 2011-05-09 15:33:35 UTC (rev 46222)
@@ -1,6 +1,8 @@
#include <unistd.h>
+#include <string.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
*
@@ -8,82 +10,95 @@
*
*/
-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;
- char buf[2048] = "";
+ double row_idx, col_idx;
+ int infd, cell_size, outfd;
+ void *trast, *tptr;
+ double n1, e1, nx, ex;
+ struct cache *ibuffer;
select_current_env();
if (G_get_cellhd(name, mapset, &cellhd) < 0)
return 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
- *
- * 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();
- G_set_window(&target_window);
- G_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
*/
G_set_window(&cellhd);
infd = G_open_cell_old(name, mapset);
if (infd < 0) {
- close(infd);
return 0;
}
map_type = G_get_raster_map_type(infd);
- rast = (void *)G_calloc(G_window_cols() + 1, G_raster_size(map_type));
- G_set_null_value(rast, G_window_cols() + 1, map_type);
+ cell_size = G_raster_size(map_type);
- G_copy(&win, &target_window, sizeof(win));
+ ibuffer = readcell(infd, seg_mb);
- win.west += win.ew_res / 2;
+ G_close_cell(infd); /* (pmx) 17 april 2000 */
+
+ G_message(_("Rectify <%s@%s> (location <%s>)"),
+ name, mapset, G_location());
+ select_target_env();
+ G_set_window(&target_window);
+ G_message(_("into <%s@%s> (location <%s>) ..."),
+ result, G_mapset(), G_location());
+
+ nrows = target_window.rows;
ncols = target_window.cols;
- col = 0;
- temp_fd = 0;
+ if (strcmp(interp_method, "nearest") != 0) {
+ map_type = DCELL_TYPE;
+ cell_size = G_raster_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 = G_open_raster_new(result, map_type);
+ trast = G_allocate_raster_buf(map_type);
- compute_georef_matrix(&cellhd, &win, order);
- perform_georef(infd, rast);
- write_matrix(row, col);
+ for (row = 0; row < nrows; row++) {
+ n1 = target_window.north - (row + 0.5) * target_window.ns_res;
- nrows -= win.rows;
- row += win.rows;
- win.north -= (win.ns_res * win.rows);
+ G_percent(row, nrows, 2);
+
+ G_set_null_value(trast, ncols, map_type);
+ tptr = trast;
+ for (col = 0; col < ncols; col++) {
+ e1 = target_window.west + (col + 0.5) * 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);
+ G_put_raster_row(outfd, trast, map_type);
}
+ G_percent(1, 1, 1);
- select_target_env();
+ G_close_cell(outfd); /* (pmx) 17 april 2000 */
+ G_free(trast);
+ close(ibuffer->fd);
+ G_free(ibuffer);
+
+ if (G_get_cellhd(result, G_mapset(), &cellhd) < 0)
+ return 0;
+
if (cellhd.proj == 0) { /* x,y imagery */
cellhd.proj = target_window.proj;
cellhd.zone = target_window.zone;
@@ -101,12 +116,7 @@
name, mapset);
}
- target_window.compressed = cellhd.compressed;
- G_close_cell(infd); /* (pmx) 17 april 2000 */
- write_map(result);
select_current_env();
- G_free(rast);
-
return 1;
}
Modified: grass/branches/releasebranch_6_4/imagery/i.rectify/report.c
===================================================================
--- grass/branches/releasebranch_6_4/imagery/i.rectify/report.c 2011-05-09 14:14:19 UTC (rev 46221)
+++ grass/branches/releasebranch_6_4/imagery/i.rectify/report.c 2011-05-09 15:33:35 UTC (rev 46222)
@@ -1,20 +1,13 @@
#include <grass/glocale.h>
#include "global.h"
-int report(char *name, char *mapset, char *result,
- long rectify, long compress, int ok)
+int report(long rectify, int ok)
{
int minutes, hours;
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;
@@ -27,30 +20,14 @@
G_verbose_message(_("%d rows, %d cols (%ld cells) completed in"),
target_window.rows, target_window.cols, ncells);
if (hours)
- G_verbose_message("%d:%02d:%02ld", hours, minutes, seconds % 60);
+ G_verbose_message(_("%d:%02d:%02ld hours"), hours, minutes, seconds % 60);
else
- G_verbose_message("%d:%02ld", minutes, seconds % 60);
+ G_verbose_message(_("%d:%02ld minutes"), minutes, seconds % 60);
if (seconds)
G_verbose_message(_("%.1f cells per minute"),
(60.0 * ncells) / ((double)seconds));
- seconds = compress;
-
- if (seconds <= 0) {
- G_message("-----------------------------------------------");
- return 1;
- }
-
- minutes = seconds / 60;
- hours = minutes / 60;
- minutes -= hours * 60;
- G_verbose_message(_("data compression required an additional"));
- if (hours)
- G_verbose_message("%d:%02d:%02ld", hours, minutes, seconds % 60);
- else
- G_verbose_message("%d:%02ld\n", minutes, seconds % 60);
-
G_message("-----------------------------------------------");
- return 0;
+ return 1;
}
Deleted: grass/branches/releasebranch_6_4/imagery/i.rectify/rowcol.h
===================================================================
--- grass/branches/releasebranch_6_4/imagery/i.rectify/rowcol.h 2011-05-09 14:14:19 UTC (rev 46221)
+++ grass/branches/releasebranch_6_4/imagery/i.rectify/rowcol.h 2011-05-09 15:33:35 UTC (rev 46222)
@@ -1 +0,0 @@
-#define ROWCOL short
Deleted: grass/branches/releasebranch_6_4/imagery/i.rectify/write.c
===================================================================
--- grass/branches/releasebranch_6_4/imagery/i.rectify/write.c 2011-05-09 14:14:19 UTC (rev 46221)
+++ grass/branches/releasebranch_6_4/imagery/i.rectify/write.c 2011-05-09 15:33:35 UTC (rev 46222)
@@ -1,66 +0,0 @@
-#include <sys/types.h>
-#include <sys/stat.h>
-#include <fcntl.h>
-#include <unistd.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) * G_raster_size(map_type);
- lseek(temp_fd, offset, SEEK_SET);
-
- if (write(temp_fd, cell_buf[n], G_raster_size(map_type) * matrix_cols)
- != G_raster_size(map_type) * matrix_cols) {
- unlink(temp_name);
- G_fatal_error(_("Error while writing to temp file"));
- }
- /*G_put_map_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;
-
- G_set_window(&target_window);
-
- rast = G_allocate_raster_buf(map_type);
- close(temp_fd);
- temp_fd = open(temp_name, 0);
- fd = G_open_raster_new(name, map_type);
-
- if (fd <= 0)
- G_fatal_error(_("Unable to create raster map <%s>"), name);
-
- for (row = 0; row < target_window.rows; row++) {
- if (read(temp_fd, rast, target_window.cols * G_raster_size(map_type))
- != target_window.cols * G_raster_size(map_type))
- G_fatal_error(_("Error writing row %d"), row);
- if (G_put_raster_row(fd, rast, map_type) < 0) {
- G_fatal_error(_("Failed writing raster map <%s> row %d"),
- name, row);
- unlink(temp_name);
- }
- }
- close(temp_fd);
- unlink(temp_name);
- G_close_cell(fd);
-
- return 0;
-}
More information about the grass-commit
mailing list