[GRASS-SVN] r46177 - grass/trunk/raster/r.resamp.bspline
svn_grass at osgeo.org
svn_grass at osgeo.org
Wed May 4 03:43:26 EDT 2011
Author: mmetz
Date: 2011-05-04 00:43:26 -0700 (Wed, 04 May 2011)
New Revision: 46177
Modified:
grass/trunk/raster/r.resamp.bspline/bspline.h
grass/trunk/raster/r.resamp.bspline/crosscorr.c
grass/trunk/raster/r.resamp.bspline/main.c
grass/trunk/raster/r.resamp.bspline/resamp.c
Log:
fix bspline raster interpolation (module)
Modified: grass/trunk/raster/r.resamp.bspline/bspline.h
===================================================================
--- grass/trunk/raster/r.resamp.bspline/bspline.h 2011-05-04 07:43:01 UTC (rev 46176)
+++ grass/trunk/raster/r.resamp.bspline/bspline.h 2011-05-04 07:43:26 UTC (rev 46177)
@@ -40,12 +40,11 @@
};
/* resamp.c */
-struct Point *P_Read_Raster_Region_Nulls(double **, /**/
- char **, /**/
- struct Cell_head *, /**/
- struct bound_box, /**/
- struct bound_box, /**/
- int *, /**/ int, /**/ double);
+struct Point *P_Read_Raster_Region_masked(char **, /**/
+ struct Cell_head *, /**/
+ struct bound_box, /**/
+ struct bound_box, /**/
+ int *, /**/ int, /**/ double);
double **P_Sparse_Raster_Points(double **, /**/
struct Cell_head *, /**/
struct Cell_head *, /**/
Modified: grass/trunk/raster/r.resamp.bspline/crosscorr.c
===================================================================
--- grass/trunk/raster/r.resamp.bspline/crosscorr.c 2011-05-04 07:43:01 UTC (rev 46176)
+++ grass/trunk/raster/r.resamp.bspline/crosscorr.c 2011-05-04 07:43:26 UTC (rev 46177)
@@ -129,7 +129,7 @@
*/
{
int bilin = TRUE; /*booleans */
- int nsplx, nsply, nparam_spl, ndata, nnulls;
+ int nsplx, nsply, nparam_spl, ndata;
double *mean, *rms, *stdev;
#ifdef nodef
@@ -164,7 +164,7 @@
/*Cats = Vect_new_cats_struct (); */
/* Current region is read and points recorded into observ */
- observ = P_Read_Raster_Region_Map(matrix, NULL, ®ion, src_reg, &ndata, &nnulls, 1024);
+ observ = P_Read_Raster_Region_Map(matrix, ®ion, src_reg, &ndata, 1024);
G_debug(5, "CrossCorrelation: %d points read in region. ", ndata);
G_verbose_message(_("%d points read in region"),
ndata);
@@ -196,8 +196,8 @@
BW = P_get_BandWidth(bilin, nsply);
/**/
- /*Least Squares system */
- N = G_alloc_matrix(nparam_spl, BW); /* Normal matrix */
+ /*Least Squares system */
+ N = G_alloc_matrix(nparam_spl, BW); /* Normal matrix */
TN = G_alloc_vector(nparam_spl); /* vector */
parVect = G_alloc_vector(nparam_spl); /* Parameters vector */
obsVect = G_alloc_matrix(ndata, 3); /* Observation vector */
Modified: grass/trunk/raster/r.resamp.bspline/main.c
===================================================================
--- grass/trunk/raster/r.resamp.bspline/main.c 2011-05-04 07:43:01 UTC (rev 46176)
+++ grass/trunk/raster/r.resamp.bspline/main.c 2011-05-04 07:43:26 UTC (rev 46177)
@@ -65,7 +65,7 @@
struct bound_box last_overlap_box, last_general_box;
struct Point *observ;
- struct Point *observ_null;
+ struct Point *observ_masked;
/*----------------------------------------------------------------*/
/* Options declarations */
@@ -194,10 +194,8 @@
nsplx_adj = NSPLX_MAX;
nsply_adj = NSPLY_MAX;
- if (stepN > stepE)
- dims.overlap = OVERLAP_SIZE * stepN;
- else
- dims.overlap = OVERLAP_SIZE * stepE;
+
+ dims.overlap = OVERLAP_SIZE * (stepN > stepE ? stepN : stepE);
P_get_edge(interp_method, &dims, stepE, stepN);
P_set_dim(&dims, stepE, stepN, &nsplx_adj, &nsply_adj);
@@ -255,8 +253,8 @@
}
/* switch to buffered input raster window */
- /* G_set_window(&src_reg); */
- Rast_set_window(&src_reg);
+ G_set_window(&src_reg);
+ /* Rast_set_window(&src_reg); */
G_debug(1, "new source north %f", src_reg.north);
G_debug(1, "new source south %f", src_reg.south);
@@ -286,7 +284,7 @@
G_percent(row, nrows, 2);
- Rast_get_d_row(inrastfd, drastbuf, row);
+ Rast_get_d_row_nomask(inrastfd, drastbuf, row);
for (col = 0; col < ncols; col++) {
inrast_matrix[row][col] = drastbuf[col];
@@ -304,8 +302,8 @@
Rast_close(inrastfd);
/* switch back to destination = current window */
- /* G_set_window(&dest_reg); */
- Rast_set_window(&dest_reg);
+ G_set_window(&dest_reg);
+ /* Rast_set_window(&dest_reg); */
nrows = Rast_window_rows();
ncols = Rast_window_cols();
@@ -332,49 +330,93 @@
}
}
- /* Alloc raster matrix */
- if (!(outrast_matrix = G_alloc_matrix(nrows, ncols)))
- G_fatal_error(_("Cannot allocate memory for auxiliar matrix."
- "Consider changing region (resolution)"));
-
/* Alloc and load masking matrix */
- if (mask_opt->answer) {
+ /* encoding: 0 = do not interpolate, 1 = interpolate */
+ if (mask_opt->answer || null_flag->answer) {
int maskfd;
+ int null_count = 0;
DCELL dval;
+ CELL cval;
+ CELL *maskbuf;
- G_message(_("Load masking map"));
-
+ G_message(_("Mark cells for interpolation"));
+
+ /* use destination window */
+
mask_matrix = (char **)G_calloc(nrows, sizeof(char *));
mask_matrix[0] = (char *)G_calloc(nrows * ncols, sizeof(char));
for (row = 1; row < nrows; row++)
mask_matrix[row] = mask_matrix[row - 1] + ncols;
+ if (mask_opt->answer) {
+ maskfd = Rast_open_old(mask_opt->answer, "");
+ maskbuf = Rast_allocate_buf(CELL_TYPE);
+ }
+ else {
+ maskfd = -1;
+ maskbuf = NULL;
+ }
- maskfd = Rast_open_old(mask_opt->answer, "");
- drastbuf = Rast_allocate_buf(DCELL_TYPE);
+ if (null_flag->answer) {
+ inrastfd = Rast_open_old(in_opt->answer, "");
+ drastbuf = Rast_allocate_buf(DCELL_TYPE);
+ }
+ else {
+ inrastfd = -1;
+ drastbuf = NULL;
+ }
for (row = 0; row < nrows; row++) {
G_percent(row, nrows, 2);
- Rast_get_d_row(maskfd, drastbuf, row);
+
+ if (mask_opt->answer)
+ Rast_get_c_row(maskfd, maskbuf, row);
+
+ if (null_flag->answer)
+ Rast_get_d_row(inrastfd, drastbuf, row);
+
for (col = 0; col < ncols; col++) {
- dval = drastbuf[col];
- if (Rast_is_d_null_value(&dval) || dval == 0)
- mask_matrix[row][col] = 0;
- else
- mask_matrix[row][col] = 1;
+ mask_matrix[row][col] = 1;
+
+ if (mask_opt->answer) {
+ cval = maskbuf[col];
+ if (Rast_is_c_null_value(&cval) || cval == 0)
+ mask_matrix[row][col] = 0;
+ }
+
+ if (null_flag->answer && mask_matrix[row][col] == 1) {
+ dval = drastbuf[col];
+ if (!Rast_is_d_null_value(&dval))
+ mask_matrix[row][col] = 0;
+ else
+ null_count++;
+ }
}
}
G_percent(row, nrows, 2);
- G_free(drastbuf);
- Rast_close(maskfd);
+ if (null_flag->answer) {
+ G_free(drastbuf);
+ Rast_close(inrastfd);
+ }
+ if (mask_opt->answer) {
+ G_free(maskbuf);
+ Rast_close(maskfd);
+ }
+ if (null_flag->answer && null_count == 0 && !mask_opt->answer) {
+ G_fatal_error(_("No NULL cells found in input raster."));
+ }
}
else
mask_matrix = NULL;
+ /* Alloc raster matrix */
+ if (!(outrast_matrix = G_alloc_matrix(nrows, ncols)))
+ G_fatal_error(_("Cannot allocate memory for auxiliar matrix."
+ "Consider changing region (resolution)"));
/* initialize output */
- G_message("initializing output");
+ G_message(_("Initializing output..."));
{
DCELL dval;
@@ -445,7 +487,8 @@
general_box.E = dest_box.W;
while (last_column == FALSE) { /* For each subregion column */
- int npoints = 0, npoints_null = 0, n_nulls = 0;
+ int npoints = 0;
+ int npoints_marked = 1; /* default: all points in output */
subregion_col++;
subregion++;
@@ -505,11 +548,10 @@
dim_vect = nsplx * nsply;
observ =
- P_Read_Raster_Region_Map(inrast_matrix, mask_matrix, &elaboration_reg,
- &src_reg, &npoints, &n_nulls,
- dim_vect);
+ P_Read_Raster_Region_Map(inrast_matrix, &elaboration_reg,
+ &src_reg, &npoints, dim_vect);
- G_debug(1, "%d points, %d NULL cells", npoints, n_nulls);
+ G_debug(1, "%d valid points", npoints);
G_debug(1,
"Interpolation: (%d,%d): Number of points in <elaboration_box> is %d",
@@ -523,30 +565,27 @@
G_debug(1, "Interpolation: (%d,%d): mean=%lf",
subregion_row, subregion_col, mean);
- observ_null = NULL;
- if (null_flag->answer && n_nulls) {
- /* read input NULL cells */
+ observ_masked = NULL;
+
+ if (mask_matrix) {
+ /* collect unmasked output cells */
- G_debug(1, "read input NULL cells");
+ G_debug(1, "collect unmasked output cells");
- observ_null =
- P_Read_Raster_Region_Nulls(inrast_matrix, mask_matrix, &src_reg,
+ observ_masked =
+ P_Read_Raster_Region_masked(mask_matrix, &dest_reg,
dest_box, general_box,
- &npoints_null, dim_vect, mean);
+ &npoints_marked, dim_vect, mean);
- G_debug(1, "%d nulls in elaboration, %d nulls in general", n_nulls, npoints_null);
- if (npoints_null == 0) {
- G_free(observ_null);
- n_nulls = 0;
+ G_debug(1, "%d cells marked in general", npoints_marked);
+ if (npoints_marked == 0) {
+ G_free(observ_masked);
+ observ_masked = NULL;
+ npoints = 1; /* disable warning below */
}
}
- else if (npoints == 0 && n_nulls == 0)
- /* nothing to interpolate, disable warning below */
- npoints = 1;
- else
- n_nulls = 1;
- if (npoints > 0 && n_nulls > 0) { /* */
+ if (npoints > 0 && npoints_marked > 0) { /* */
int i;
nparameters = nsplx * nsply;
@@ -597,33 +636,30 @@
G_free_vector(TN);
G_free_vector(Q);
- if (!null_flag->answer) { /* interpolate full output raster */
+ if (!observ_masked) { /* interpolate full output raster */
G_debug(1, "Interpolation: (%d,%d): Regular_Points...",
subregion_row, subregion_col);
outrast_matrix =
P_Regular_Points(&elaboration_reg, &dest_reg, general_box,
- overlap_box, outrast_matrix, mask_matrix,
+ overlap_box, outrast_matrix, NULL,
parVect, stepN, stepE, dims.overlap, mean,
nsplx, nsply, nrows, ncols, interp_method);
}
- else { /* only interpolate NULL cells */
+ else { /* only interpolate selected cells */
- if (observ_null == NULL)
- G_fatal_error("no NULL cells loaded");
+ G_debug(1, "Interpolation of %d selected cells...",
+ npoints_marked);
- G_debug(1, "Interpolation of %d NULL cells...",
- npoints_null);
-
outrast_matrix =
P_Sparse_Raster_Points(outrast_matrix,
&elaboration_reg, &dest_reg,
general_box, overlap_box,
- observ_null, parVect,
+ observ_masked, parVect,
stepE, stepN,
dims.overlap, nsplx, nsply,
- npoints_null, interp_method, mean);
+ npoints_marked, interp_method, mean);
- G_free(observ_null);
+ G_free(observ_masked);
} /* end NULL cells */
G_free_vector(parVect);
G_free_matrix(obsVect);
Modified: grass/trunk/raster/r.resamp.bspline/resamp.c
===================================================================
--- grass/trunk/raster/r.resamp.bspline/resamp.c 2011-05-04 07:43:01 UTC (rev 46176)
+++ grass/trunk/raster/r.resamp.bspline/resamp.c 2011-05-04 07:43:26 UTC (rev 46177)
@@ -22,7 +22,7 @@
#include <math.h>
#include "bspline.h"
-struct Point *P_Read_Raster_Region_Nulls(double **matrix, char **mask_matrix,
+struct Point *P_Read_Raster_Region_masked(char **mask_matrix,
struct Cell_head *Original,
struct bound_box output_box,
struct bound_box General,
@@ -31,20 +31,20 @@
{
int col, row, startcol, endcol, startrow, endrow, nrows, ncols;
int pippo, npoints;
- double X, Y, Z;
+ double X, Y;
struct Point *obs;
pippo = dim_vect;
obs = (struct Point *)G_calloc(pippo, sizeof(struct Point));
- /* Reading points inside input box and inside General box */
+ /* Reading points inside output box and inside General box */
npoints = 0;
nrows = Original->rows;
ncols = Original->cols;
- /* original region = input raster has been adjusted to output region plus buffer
- * -> General box is somewhere inside input raster box */
+ /* original region = output region
+ * -> General box is somewhere inside output region */
if (Original->north > General.N) {
startrow = (double)((Original->north - General.N) / Original->ns_res - 1);
if (startrow < 0)
@@ -77,35 +77,27 @@
for (row = startrow; row < endrow; row++) {
for (col = startcol; col < endcol; col++) {
- if (mask_matrix) {
- if (!mask_matrix[row][col])
- continue;
- }
-
- Z = matrix[row][col];
+ if (!mask_matrix[row][col])
+ continue;
- if (Rast_is_d_null_value(&Z)) {
-
- X = Rast_col_to_easting((double)(col) + 0.5, Original);
- Y = Rast_row_to_northing((double)(row) + 0.5, Original);
+ X = Rast_col_to_easting((double)(col) + 0.5, Original);
+ Y = Rast_row_to_northing((double)(row) + 0.5, Original);
- /* Here, mean is just for asking if obs point is in box */
- if (Vect_point_in_box(X, Y, mean, &General)) { /* General */
- npoints++;
- if (npoints >= pippo) {
- pippo += dim_vect;
- obs =
- (struct Point *)G_realloc((void *)obs,
- (signed int)pippo *
- sizeof(struct Point));
- }
-
- /* Storing observation vector */
- obs[npoints - 1].coordX = X;
- obs[npoints - 1].coordY = Y;
- obs[npoints - 1].coordZ = 0;
-
+ /* Here, mean is just for asking if obs point is in box */
+ if (Vect_point_in_box(X, Y, mean, &General)) { /* General */
+ if (npoints >= pippo) {
+ pippo += dim_vect;
+ obs =
+ (struct Point *)G_realloc((void *)obs,
+ (signed int)pippo *
+ sizeof(struct Point));
}
+
+ /* Storing observation vector */
+ obs[npoints].coordX = X;
+ obs[npoints].coordY = Y;
+ obs[npoints].coordZ = 0;
+ npoints++;
}
}
}
More information about the grass-commit
mailing list