new module for raster resampling and NULL cell filling

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

vector/lidar/r.resamp.bspline/Makefile

@@ -0,0 +1,13 @@
+MODULE_TOPDIR = ../../..
+
+PGM = r.resamp.bspline
+
+LIBES = $(LIDARLIB) $(GMATHLIB) $(VECTORLIB) $(DBMILIB) $(RASTERLIB) $(GISLIB) $(MATHLIB)
+DEPENDENCIES = $(LIDARDEP) $(GMATHDEP) $(VECTORDEP) $(DBMIDEP) $(RASTERDEP) $(GISDEP)
+EXTRA_INC = $(VECT_INC)
+EXTRA_CFLAGS = $(VECT_CFLAGS)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+default: cmd
+

vector/lidar/r.resamp.bspline/bspline.h

@@ -0,0 +1,77 @@
+
+/***********************************************************************
+ *
+ * MODULE:       v.surf.bspline
+ *
+ * AUTHOR(S):    Roberto Antolin
+ *
+ * PURPOSE:      Spline Interpolation and cross correlation
+ *
+ * COPYRIGHT:    (C) 2006 by Politecnico di Milano -
+ *			     Polo Regionale di Como
+ *
+ *               This program is free software under the
+ *               GNU General Public License (>=v2).
+ *               Read the file COPYING that comes with GRASS
+ *               for details.
+ *
+ **************************************************************************/
+
+ /*INCLUDES*/
+#include <stdlib.h>
+#include <string.h>
+#include <grass/gis.h>
+#include <grass/raster.h>
+#include <grass/vector.h>
+#include <grass/dbmi.h>
+#include <grass/glocale.h>
+#include <grass/config.h>
+#include <grass/PolimiFunct.h>
+     /*STRUCTURES*/ struct Stats
+{
+    double *estima;
+    double *error;
+    int n_points;
+};
+
+struct Param
+{
+    struct Option *in, *in_ext, *out, *out_map, *dbdriver,
+	*dbdatabase, *passoE, *passoN, *lambda_f, *type;
+    struct Flag *cross_corr;
+};
+
+struct RasterPoint
+{
+    int row;
+    int col;
+};
+
+/* resamp.c */
+struct Point *P_Read_Raster_Region_Nulls(double **, /**/
+				       struct Cell_head *, /**/
+				       struct bound_box, /**/
+				       struct bound_box, /**/
+				       int *, /**/ int, /**/ double);
+double **P_Sparse_Raster_Points(double **, /**/
+			struct Cell_head *, /**/
+			struct Cell_head *, /**/
+			struct bound_box, /**/
+			struct bound_box, /**/
+			struct Point *, /**/
+			double *, /**/
+			double, /**/
+			double, /**/
+			double, /**/
+			int, /**/
+			int, /**/
+			int, /**/
+			int, /**/
+			double /**/);
+int align_elaboration_box(struct Cell_head *, struct Cell_head *, int);
+int align_interp_boxes(struct bound_box *, struct bound_box *,
+                       struct Cell_head *, struct bound_box, struct bound_box, int);
+
+				       
+/* crosscor.c */
+int cross_correlation(double **, struct Cell_head *, double, double);

vector/lidar/r.resamp.bspline/crosscorr.c

@@ -0,0 +1,368 @@
+
+/***********************************************************************
+ *
+ * MODULE:       v.surf.bspline
+ *
+ * AUTHOR(S):    Roberto Antolin
+ *
+ * PURPOSE:      Spline Interpolation and cross correlation
+ *
+ * COPYRIGHT:    (C) 2006 by Politecnico di Milano -
+ *			     Polo Regionale di Como
+ *
+ *               This program is free software under the
+ *               GNU General Public License (>=v2).
+ *               Read the file COPYING that comes with GRASS
+ *               for details.
+ *
+ **************************************************************************/
+
+ /*INCLUDES*/
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <grass/gis.h>
+#include <grass/vector.h>
+#include <grass/dbmi.h>
+#include <grass/glocale.h>
+#include <grass/config.h>
+#include <grass/PolimiFunct.h>
+#include "bspline.h"
+#define NDATA_MAX 100
+#define PARAM_LAMBDA 6
+#define PARAM_SPLINE 0
+#define SWAP(a, b) {double t=a; a=b; b=t;}
+
+/* internal functions */
+double calc_mean(double *values, int nvalues)
+{
+    int i;
+    double sum = .0;
+
+    if (nvalues == 0)
+	return .0;
+    for (i = 0; i < nvalues; i++)
+	sum += values[i];
+    return sum / nvalues;
+}
+
+
+double calc_root_mean_square(double *values, int nvalues)
+{
+    int i;
+    double rms, sum = .0;
+
+    if (nvalues == 0)
+	return .0;
+
+    for (i = 0; i < nvalues; i++)
+	sum += pow(values[i], 2) / nvalues;
+
+    rms = sqrt(sum);
+    return rms;
+
+}
+
+double calc_standard_deviation(double *values, int nvalues)
+{
+    double mean, rms, stdev;
+
+    if (nvalues == 0)
+	return .0;
+
+    rms = calc_root_mean_square(values, nvalues);
+    mean = calc_mean(values, nvalues);
+
+    stdev = sqrt(pow(rms, 2) - pow(mean, 2));
+    return stdev;
+}
+
+struct Stats alloc_Stats(int n)
+{
+    double *err, *stm;
+    struct Stats stat;
+
+    stat.n_points = n;
+    err = (double *)G_calloc(n, sizeof(double));
+    stm = (double *)G_calloc(n, sizeof(double));
+
+    stat.error = err;
+    stat.estima = stm;
+
+    return stat;
+}
+
+double find_minimum(double *values, int *l_min)
+{
+    int l;
+    double min;
+
+    min = values[0];
+
+    for (l = 0; l < PARAM_LAMBDA; l++) {
+	if (min > values[l]) {
+	    min = values[l];
+	    *l_min = l;
+	}
+    }
+    return min;
+}
+
+struct Point *swap(struct Point *point, int a, int b)
+{
+
+    SWAP(point[a].coordX, point[b].coordX);	/* Once the last value is let out, it is swap with j-value */
+    SWAP(point[a].coordY, point[b].coordY);
+    SWAP(point[a].coordZ, point[b].coordZ);
+    SWAP(point[a].cat, point[b].cat);
+
+    return point;
+}
+
+/*-------------------------------------------------------------------------------------------*/
+int cross_correlation(double **matrix, struct Cell_head *src_reg, double passWE, double passNS)
+    /*
+       matrix: matrix with raster values
+       passWE: spline step in West-East direction
+       passNS: spline step in North-South direction
+
+       RETURN:
+       TRUE on success
+       FALSE on failure
+     */
+{
+    int bilin = TRUE;		/*booleans */
+    int nsplx, nsply, nparam_spl, ndata, nnulls;
+    double *mean, *rms, *stdev;
+
+#ifdef nodef
+    double rms_min, stdev_min;
+    int lbd_min; /* lbd_min: index where minimun lambda index is found */
+#endif
+
+    /* double lambda[PARAM_LAMBDA] = { 0.0001, 0.001, 0.01, 0.1, 1.0, 10.0 }; */	/* Fixed values (by the moment) */
+    double lambda[PARAM_LAMBDA] = { 0.0001, 0.001, 0.005, 0.01, 0.02, 0.05 };	/* Fixed values (by the moment) */
+
+    double *TN, *Q, *parVect;	/* Interpolation and least-square vectors */
+    double **N, **obsVect;	/* Interpolation and least-square matrix */
+
+    struct Point *observ;
+    struct Stats stat_vect;
+
+    struct Cell_head region;
+
+    G_get_window(&region);
+
+
+    G_debug(5,
+	    "CrossCorrelation: Some tests using different lambda_i values will be done");
+
+    ndata = region.rows * region.cols;
+
+    if (ndata > NDATA_MAX)
+	G_warning(_("%d are too many cells, recommended are < 100 cells. "
+		    "The cross validation would take too much time."), ndata);
+
+    /*points = Vect_new_line_struct (); */
+    /*Cats = Vect_new_cats_struct (); */
+
+    /* Current region is read and points recorded into observ */
+    observ = P_Read_Raster_Region_Map(matrix, &region, src_reg, &ndata, &nnulls, 1024);
+    G_debug(5, "CrossCorrelation: %d points read in region. ", ndata);
+    G_verbose_message(_("%d points read in region"),
+		      ndata);
+
+    if (ndata > 50)
+	G_warning(_("Maybe it takes too long."
+		    "Consider reducing the region extends."));
+    else
+	G_debug(5, "CrossCorrelation: It shouldn't take too long.");
+
+    if (ndata > 0) {		/* If at least one point is in the region */
+	int i, j, lbd;		/* lbd: lambda index */
+	int BW;	
+	double mean_reg, *obs_mean;
+
+	mean = G_alloc_vector(PARAM_LAMBDA);	/* Alloc as much mean, rms and stdev values as the total */
+	rms = G_alloc_vector(PARAM_LAMBDA);	/* number of parameter used used for cross validation */
+	stdev = G_alloc_vector(PARAM_LAMBDA);
+
+	/* Setting number of splines as a function of WE and SN spline steps */
+	nsplx = ceil((region.east - region.west) / passWE);
+	nsply = ceil((region.north - region.south) / passNS);
+	nparam_spl = nsplx * nsply;	/* Total number of splines */
+
+	if (nparam_spl > 22900)
+	    G_fatal_error(_("Too many splines (%d x %d). "
+			    "Consider changing spline steps \"sie=\" \"sin=\"."),
+			  nsplx, nsply);
+
+	BW = P_get_BandWidth(bilin, nsply);
+	/**/
+	    /*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 */
+	Q = G_alloc_vector(ndata);	/* "a priori" var-cov matrix */
+
+	obs_mean = G_alloc_vector(ndata);
+	stat_vect = alloc_Stats(ndata);
+
+	for (lbd = 0; lbd < PARAM_LAMBDA; lbd++) {	/* For each lambda value */
+
+	    G_message(_("Begining cross validation with "
+			"lambda_i=%.4f..."), lambda[lbd]);
+	    
+	    /*
+	       How cross correlation algorithm is done:
+	       For each cicle, only the first ndata-1 "observ" elements are considered for the 
+	       interpolation. Within every interpolation mean is calculated to lowering border 
+	       errors. The point let out will be used for an estimation. The error between the 
+	       estimation and the observation is recorded for further statistics.
+	       At the end of the cicle, the last point, that is, the ndata-1 index, and the point 
+	       with j index are swapped.
+	     */
+	    for (j = 0; j < ndata; j++) {	/* Cross Correlation will use all ndata points */
+		double out_x, out_y, out_z;	/* This point is let out */
+
+		for (i = 0; i < ndata; i++) {	/* Each time, only the first ndata-1 points */
+
+		    /* Setting obsVect vector & Q matrix */
+		    Q[i] = 1;	/* Q=I */
+		    obsVect[i][0] = observ[i].coordX;
+		    obsVect[i][1] = observ[i].coordY;
+
+		    obsVect[i][2] = observ[i].coordZ;
+		    obs_mean[i] = observ[i].coordZ;
+		}		/* i index */
+
+		/* Mean calculation for every point less the last one */
+		mean_reg = calc_mean(obs_mean, ndata - 1);
+
+		for (i = 0; i < ndata; i++)
+		    obsVect[i][2] -= mean_reg;
+
+		/* This is let out */
+		out_x = observ[ndata - 1].coordX;
+		out_y = observ[ndata - 1].coordY;
+		out_z = obsVect[ndata - 1][2];
+
+		if (bilin) {	/* Bilinear interpolation */
+		    normalDefBilin(N, TN, Q, obsVect, passWE, passNS, nsplx,
+				   nsply, region.west, region.south,
+				   ndata - 1, nparam_spl, BW);
+		    nCorrectGrad(N, lambda[lbd], nsplx, nsply, passWE,
+				 passNS);
+		}
+		else {		/* Bicubic interpolation */
+		    normalDefBicubic(N, TN, Q, obsVect, passWE, passNS, nsplx,
+				     nsply, region.west, region.south,
+				     ndata - 1, nparam_spl, BW);
+		    nCorrectGrad(N, lambda[lbd], nsplx, nsply, passWE,
+				 passNS);
+		}
+
+		/* 
+		   if (bilin) interpolation (&interp, P_BILINEAR);
+		   else interpolation (&interp, P_BICUBIC);
+		 */
+		G_math_solver_cholesky_sband(N, parVect, TN, nparam_spl, BW);
+
+		/* Estimation of j-point */
+		if (bilin)
+		    stat_vect.estima[j] =
+			dataInterpolateBilin(out_x, out_y, passWE, passNS,
+					     nsplx, nsply, region.west,
+					     region.south, parVect);
+
+		else
+		    stat_vect.estima[j] =
+			dataInterpolateBilin(out_x, out_y, passWE, passNS,
+					     nsplx, nsply, region.west,
+					     region.south, parVect);
+
+		/*Difference between estimated and observated i-point */
+		stat_vect.error[j] = out_z - stat_vect.estima[j];
+		G_debug(1, "CrossCorrelation: stat_vect.error[%d]  =  %lf", j,
+			stat_vect.error[j]);
+
+		observ = swap(observ, j, ndata - 1);	/* Once the last value is let out, it is swap with j-value */
+	    }
+
+	    mean[lbd] = calc_mean(stat_vect.error, stat_vect.n_points);
+	    rms[lbd] =
+		calc_root_mean_square(stat_vect.error, stat_vect.n_points);
+	    stdev[lbd] =
+		calc_standard_deviation(stat_vect.error, stat_vect.n_points);
+
+	    G_message(_("Mean = %.5lf"), mean[lbd]);
+	    G_message(_("Root Means Square (RMS) = %.5lf"),
+		      rms[lbd]);
+	}			/* ENDFOR each lambda value */
+
+	G_free_matrix(N);
+	G_free_vector(TN);
+	G_free_vector(Q);
+	G_free_matrix(obsVect);
+	G_free_vector(parVect);
+#ifdef nodef
+	/*TODO: if the minimum lambda is wanted, the function declaration must be changed */
+	/*At this moment, consider rms only */
+	rms_min = find_minimum(rms, &lbd_min);
+	stdev_min = find_minimum(stdev, &lbd_min);
+
+	/* Writing some output */
+	G_message(_("Different number of splines and lambda_i values have "
+		    "been taken for the cross correlation"));
+	G_message(_("The minimum value for the test (rms=%lf) was "
+		    "obtained with: lambda_i = %.3f"),
+		  rms_min,
+		  lambda[lbd_min]);
+
+	*lambda_min = lambda[lbd_min];
+#endif
+
+	G_message(_("Results into a table:"));
+	G_message(_(" lambda    | mean        | rms         |"));
+	for (lbd = 0; lbd < PARAM_LAMBDA; lbd++) {
+	    G_message(_(" %-10.5f| %-12.4f| %-12.4f|"), lambda[lbd],
+		      mean[lbd], rms[lbd]);
+	}
+	
+	G_free_vector(mean);
+	G_free_vector(rms);
+    }				/* ENDIF (ndata > 0) */
+    else
+	G_warning(_("No point lies into the current region"));
+
+    G_free(observ);
+    return TRUE;
+}
+
+#ifdef nodef
+void interpolation(struct ParamInterp *interp, boolean bilin)
+{
+    if (bilin == P_BILINEAR) {	/* Bilinear interpolation */
+	normalDefBilin(interp->N, interp->TN, interp->Q, interp->obsVect,
+		       interp->passoE, interp->passoN, interp->nsplx,
+		       interp->nsply, interp->region.west,
+		       interp->region.south, interp->ndata,
+		       interp->nparam_spl, interp->BW);
+
+	nCorrectGrad(interp->N, interp->lambda[lbd], interp->nsplx,
+		     interp->nsply, interp->passoE, interp->passoN);
+    }
+    else {			/* Bicubic interpolation */
+	normalDefBicubic(interp->N, interp->TN, interp->Q, interp->obsVect,
+			 interp->passoE, interp->passoN, interp->nsplx,
+			 interp->nsply, interp->region.west,
+			 interp->region.south, interp->ndata,
+			 interp->nparam_spl, interp->BW);
+
+	nCorrectGrad(interp->N, interp->lambda[lbd], interp->nsplx,
+		     interp->nsply, interp->passoE, interp->passoN);
+    }
+    return TRUE;
+}
+#endif

vector/lidar/r.resamp.bspline/main.c

@@ -0,0 +1,647 @@
+
+/**********************************************************************
+ *
+ * MODULE:       r.resamp.bspline
+ *
+ * AUTHOR(S):    Markus Metz
+ *
+ * PURPOSE:      Spline Interpolation
+ *
+ * COPYRIGHT:    (C) 2010 GRASS development team
+ *
+ *               This program is free software under the
+ *               GNU General Public License (>=v2).
+ *               Read the file COPYING that comes with GRASS
+ *               for details.
+ *
+ **********************************************************************/
+
+/* INCLUDES */
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <grass/gis.h>
+#include <grass/raster.h>
+#include <grass/vector.h>
+#include <grass/dbmi.h>
+#include <grass/glocale.h>
+#include <grass/config.h>
+#include <grass/PolimiFunct.h>
+#include "bspline.h"
+
+/*--------------------------------------------------------------------*/
+int main(int argc, char *argv[])
+{
+    /* Variable declarations */
+    int nsply, nsplx, row, col, nrows, ncols, nsplx_adj, nsply_adj;
+    int nsubregion_col, nsubregion_row, subregion_row, subregion_col;
+    int subregion = 0, nsubregions = 0;
+    int last_row, last_column, interp_method;	/* booleans */
+    double lambda, mean;
+    double N_extension, E_extension, edgeE, edgeN;
+    double stepN, stepE;
+
+    const char *inrast, *outrast;
+    char title[64];
+
+    int dim_vect, nparameters, BW;
+    double **inrast_matrix, **outrast_matrix;	/* Matrix to store the auxiliar raster values */
+    double *TN, *Q, *parVect;	/* Interpolating and least-square vectors */
+    double **N, **obsVect;	/* Interpolation and least-square matrix */
+
+    /* Structs declarations */
+    int inrastfd, outrastfd;
+    DCELL *drastbuf;
+    struct History history;
+
+    struct Map_info Grid;
+    struct line_pnts *Points;
+    struct line_cats *Cats;
+    int cat = 1;
+
+    struct GModule *module;
+    struct Option *in_opt, *out_opt, *grid_opt, *ogrid_opt, *stepE_opt, *stepN_opt,
+		  *lambda_f_opt, *method_opt;
+    struct Flag *null_flag, *cross_corr_flag;
+
+    struct Reg_dimens dims;
+    struct Cell_head elaboration_reg, src_reg, dest_reg;
+    struct bound_box general_box, overlap_box, dest_box;
+    struct bound_box last_overlap_box, last_general_box;
+
+    struct Point *observ;
+    struct Point *observ_null;
+
+    /*----------------------------------------------------------------*/
+    /* Options declarations */
+    module = G_define_module();
+    G_add_keyword(_("raster"));
+    G_add_keyword(_("surface"));
+    G_add_keyword(_("resampling"));
+    module->description =
+	_("Bicubic or bilinear spline interpolation with Tykhonov regularization.");
+
+    in_opt = G_define_standard_option(G_OPT_R_INPUT);
+
+    out_opt = G_define_standard_option(G_OPT_R_OUTPUT);
+
+    grid_opt = G_define_standard_option(G_OPT_V_OUTPUT);
+    grid_opt->key = "grid";
+    grid_opt->required = NO;
+
+    ogrid_opt = G_define_standard_option(G_OPT_V_OUTPUT);
+    ogrid_opt->key = "ogrid";
+    ogrid_opt->required = NO;
+
+    stepE_opt = G_define_option();
+    stepE_opt->key = "se";
+    stepE_opt->type = TYPE_DOUBLE;
+    stepE_opt->required = NO;
+    stepE_opt->description =
+	_("Length of each spline step in the east-west direction. Default: ewres.");
+    stepE_opt->guisection = _("Settings");
+
+    stepN_opt = G_define_option();
+    stepN_opt->key = "sn";
+    stepN_opt->type = TYPE_DOUBLE;
+    stepN_opt->required = NO;
+    stepN_opt->description =
+	_("Length of each spline step in the north-south direction. Default: nsres.");
+    stepN_opt->guisection = _("Settings");
+
+    method_opt = G_define_option();
+    method_opt->key = "method";
+    method_opt->type = TYPE_STRING;
+    method_opt->required = NO;
+    method_opt->description = _("Spline interpolation algorithm");
+    method_opt->options = "bilinear,bicubic";
+    method_opt->answer = "bicubic";
+    method_opt->guisection = _("Settings");
+
+    lambda_f_opt = G_define_option();
+    lambda_f_opt->key = "lambda";
+    lambda_f_opt->type = TYPE_DOUBLE;
+    lambda_f_opt->required = NO;
+    lambda_f_opt->description = _("Tykhonov regularization parameter (affects smoothing)");
+    lambda_f_opt->answer = "0.001";
+    lambda_f_opt->guisection = _("Settings");
+
+    null_flag = G_define_flag();
+    null_flag->key = 'n';
+    null_flag->label = _("Only interpolate null cells in input raster map");
+    null_flag->guisection = _("Settings");
+
+    cross_corr_flag = G_define_flag();
+    cross_corr_flag->key = 'c';
+    cross_corr_flag->description =
+	_("Find the best Tykhonov regularizing parameter using a \"leave-one-out\" cross validation method");
+
+    /*----------------------------------------------------------------*/
+    /* Parsing */
+    G_gisinit(argv[0]);
+    if (G_parser(argc, argv))
+	exit(EXIT_FAILURE);
+
+    inrast = in_opt->answer;
+    outrast = out_opt->answer;
+
+    if (!strcmp(method_opt->answer, "bilinear"))
+	interp_method = P_BILINEAR;
+    else
+	interp_method = P_BICUBIC;
+
+    lambda = atof(lambda_f_opt->answer);
+
+    /* Setting regions and boxes */
+    G_debug(1, "Interpolation: Setting regions and boxes");
+    G_get_set_window(&dest_reg);
+    G_get_set_window(&elaboration_reg);
+    Vect_region_box(&dest_reg, &dest_box);
+    Vect_region_box(&elaboration_reg, &overlap_box);
+    Vect_region_box(&elaboration_reg, &general_box);
+
+    nrows = Rast_window_rows();
+    ncols = Rast_window_cols();
+
+    /* get window of input map */
+    Rast_get_cellhd(in_opt->answer, "", &src_reg);
+
+    if (stepE_opt->answer) {
+	stepE = atof(stepE_opt->answer);
+	if (stepE <= .0)
+	    G_fatal_error(_("se must be positive"));
+    }
+    else
+        stepE = src_reg.ew_res * 1.5;
+
+    if (stepN_opt->answer) {
+	stepN = atof(stepN_opt->answer);
+	if (stepN <= .0)
+	    G_fatal_error(_("sn must be positive"));
+    }
+    else
+        stepN = src_reg.ns_res * 1.5;
+
+    /*------------------------------------------------------------------
+      | Subdividing and working with tiles: 									
+      | Each original region will be divided into several subregions. 
+      | Each one will be overlaped by its neighbouring subregions. 
+      | The overlapping is calculated as a fixed OVERLAP_SIZE times
+      | the largest spline step plus 2 * orlo
+      ----------------------------------------------------------------*/
+
+    /* Fixing parameters of the elaboration region */
+    P_zero_dim(&dims);		/* Set dim struct to zero */
+
+    nsplx_adj = NSPLX_MAX;
+    nsply_adj = NSPLY_MAX;
+    if (stepN > stepE)
+	dims.overlap = OVERLAP_SIZE * stepN;
+    else
+	dims.overlap = OVERLAP_SIZE * stepE;
+    P_get_edge(interp_method, &dims, stepE, stepN);
+    P_set_dim(&dims, stepE, stepN, &nsplx_adj, &nsply_adj);
+
+    G_verbose_message(_("spline step in ew direction %g"), stepE);
+    G_verbose_message(_("spline step in ns direction %g"), stepN);
+    G_verbose_message(_("adjusted EW splines %d"), nsplx_adj);
+    G_verbose_message(_("adjusted NS splines %d"), nsply_adj);
+
+    /* calculate number of subregions */
+    edgeE = dims.ew_size - dims.overlap - 2 * dims.edge_v;
+    edgeN = dims.sn_size - dims.overlap - 2 * dims.edge_h;
+
+    N_extension = dest_reg.north - dest_reg.south;
+    E_extension = dest_reg.east - dest_reg.west;
+
+    nsubregion_col = ceil(E_extension / edgeE) + 0.5;
+    nsubregion_row = ceil(N_extension / edgeN) + 0.5;
+
+    if (nsubregion_col < 0)
+	nsubregion_col = 0;
+    if (nsubregion_row < 0)
+	nsubregion_row = 0;
+
+    nsubregions = nsubregion_row * nsubregion_col;
+
+    G_debug(1, "-------------------------------------");
+    G_debug(1, "source north %f", src_reg.north);
+    G_debug(1, "source south %f", src_reg.south);
+    G_debug(1, "source west %f", src_reg.west);
+    G_debug(1, "source east %f", src_reg.east);
+    G_debug(1, "-------------------------------------");
+
+    /* adjust source window */
+    if (1) {
+	double north = dest_reg.north + 2 * dims.edge_h;
+	double south = dest_reg.south - 2 * dims.edge_h;
+	int r0 = (int)(floor(Rast_northing_to_row(north, &src_reg)) - 0.5);
+	int r1 = (int)(floor(Rast_northing_to_row(south, &src_reg)) + 0.5);
+	double east = dest_reg.east + 2 * dims.edge_v;
+	double west = dest_reg.west - 2 * dims.edge_v;
+	/* NOTE: Rast_easting_to_col() is broken because of G_adjust_easting() */
+	/*
+	int c0 = (int)floor(Rast_easting_to_col(east, &src_reg) + 0.5);
+	int c1 = (int)floor(Rast_easting_to_col(west, &src_reg) + 0.5);
+        */
+	int c0 = (int)(floor(((east - src_reg.west) / src_reg.ew_res)) + 0.5);
+	int c1 = (int)(floor(((west - src_reg.west) / src_reg.ew_res)) - 0.5);
+
+	src_reg.north -= src_reg.ns_res * (r0);
+	src_reg.south -= src_reg.ns_res * (r1 - src_reg.rows);
+	src_reg.east += src_reg.ew_res * (c0 - src_reg.cols);
+	src_reg.west += src_reg.ew_res * (c1);
+	src_reg.rows = r1 - r0;
+	src_reg.cols = c0 - c1;
+    }
+
+    /* switch to buffered input raster window */
+    /* 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);
+    G_debug(1, "new source west %f", src_reg.west);
+    G_debug(1, "new source east %f", src_reg.east);
+    G_debug(1, "-------------------------------------");
+
+    /* read raster input */
+    inrastfd = Rast_open_old(in_opt->answer, "");
+    nrows = Rast_window_rows();
+    ncols = Rast_window_cols();
+
+    G_debug(1, "%d new rows, %d new cols", nrows, ncols);
+
+    /* Alloc raster matrix */
+    if (!(inrast_matrix = G_alloc_matrix(nrows, ncols)))
+	G_fatal_error(_("Cannot allocate memory for auxiliar matrix."
+			"Consider changing region (resolution)"));
+
+    drastbuf = Rast_allocate_buf(DCELL_TYPE);
+
+    G_message("loading input raster <%s>", in_opt->answer);
+    if (1) {
+	int got_one = 0;
+	for (row = 0; row < nrows; row++) {
+	    DCELL dval;
+	    
+	    G_percent(row, nrows, 2);
+
+	    Rast_get_d_row(inrastfd, drastbuf, row);
+
+	    for (col = 0; col < ncols; col++) {
+		inrast_matrix[row][col] = drastbuf[col];
+		dval = inrast_matrix[row][col];
+		if (!Rast_is_d_null_value(&dval)) {
+		    got_one++;
+		}
+	    }
+	}
+	if (!got_one)
+	    G_fatal_error("only NULL cells in input raster");
+    }
+    G_percent(row, nrows, 2);
+    G_free(drastbuf);
+    Rast_close(inrastfd);
+
+    /* switch back to destination = current window */
+    /* G_set_window(&dest_reg); */
+    Rast_set_window(&dest_reg);
+    nrows = Rast_window_rows();
+    ncols = Rast_window_cols();
+
+    G_debug(1, "-------------------------------------");
+    G_debug(1, "dest north %f", dest_reg.north);
+    G_debug(1, "dest south %f", dest_reg.south);
+    G_debug(1, "dest west %f", dest_reg.west);
+    G_debug(1, "dest east %f", dest_reg.east);
+    G_debug(1, "-------------------------------------");
+
+    /* cross-correlation */
+    if (cross_corr_flag->answer) {
+	G_debug(1, "CrossCorrelation()");
+
+	if (cross_correlation(inrast_matrix, &src_reg, stepE, stepN) != TRUE)
+	    G_fatal_error(_("Cross validation didn't finish correctly"));
+	else {
+	    G_debug(1, "Cross validation finished correctly");
+
+	    G_free(inrast_matrix);
+
+	    G_done_msg(_("Cross validation finished for se = %f and sn = %f"), stepE, stepN);
+	    exit(EXIT_SUCCESS);
+	}
+    }
+
+    /* 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");
+    {
+	DCELL dval;
+
+	Rast_set_d_null_value(&dval, 1);
+	for (row = 0; row < nrows; row++) {
+	    G_percent(row, nrows, 2);
+	    for (col = 0; col < ncols; col++) {
+		outrast_matrix[row][col] = dval;
+	    }
+	}
+    }
+    G_percent(row, nrows, 2);
+
+    if (grid_opt->answer) {
+	if (0 > Vect_open_new(&Grid, grid_opt->answer, WITH_Z))
+	    G_fatal_error(_("Unable to create vector map <%s>"),
+			  grid_opt->answer);
+			  
+	Points = Vect_new_line_struct();
+	Cats = Vect_new_cats_struct();
+    }
+
+    subregion_row = 0;
+    elaboration_reg.south = dest_reg.north;
+    last_row = FALSE;
+    overlap_box.S = dest_box.N;
+    general_box.S = dest_box.N;
+
+    while (last_row == FALSE) {	/* For each subregion row */
+	subregion_row++;
+	last_overlap_box.S = overlap_box.S;
+	last_general_box.S = general_box.S;
+	P_set_regions(&elaboration_reg, &general_box, &overlap_box, dims,
+		      GENERAL_ROW);
+
+	/* only works if source reg = dest reg with buffer */
+	/* messing with elaboration region is dangerous... */
+	align_elaboration_box(&elaboration_reg, &src_reg, GENERAL_ROW);
+	align_interp_boxes(&general_box, &overlap_box, &dest_reg,
+	                last_general_box, last_overlap_box, GENERAL_ROW);
+
+	if (elaboration_reg.north > dest_reg.north) {	/* First row */
+
+	    P_set_regions(&elaboration_reg, &general_box, &overlap_box, dims,
+			  FIRST_ROW);
+	    align_elaboration_box(&elaboration_reg, &src_reg, GENERAL_ROW);
+	    align_interp_boxes(&general_box, &overlap_box, &dest_reg,
+			    last_general_box, last_overlap_box, FIRST_ROW);
+	}
+
+	if (elaboration_reg.south <= dest_reg.south) {	/* Last row */
+
+	    P_set_regions(&elaboration_reg, &general_box, &overlap_box, dims,
+			  LAST_ROW);
+	    last_row = TRUE;
+	}
+
+	nsply =
+	    ceil((elaboration_reg.north -
+		  elaboration_reg.south) / stepN) + 0.5;
+	G_debug(1, "Interpolation: nsply = %d", nsply);
+
+	elaboration_reg.east = dest_reg.west;
+	last_column = FALSE;
+	subregion_col = 0;
+
+	overlap_box.E = dest_box.W;
+	general_box.E = dest_box.W;
+
+	while (last_column == FALSE) {	/* For each subregion column */
+	    int npoints = 0, npoints_null = 0, n_nulls = 0;
+
+	    subregion_col++;
+	    subregion++;
+	    if (nsubregions > 1)
+		G_message(_("subregion %d of %d"), subregion, nsubregions);
+
+	    last_overlap_box.E = overlap_box.E;
+	    last_general_box.E = general_box.E;
+
+	    P_set_regions(&elaboration_reg, &general_box, &overlap_box, dims,
+			  GENERAL_COLUMN);
+
+	    /* only works if source reg = dest reg with buffer */
+	    /* messing with elaboration region is dangerous... */
+	    align_elaboration_box(&elaboration_reg, &src_reg, GENERAL_COLUMN);
+	    align_interp_boxes(&general_box, &overlap_box, &dest_reg,
+	              last_general_box, last_overlap_box, GENERAL_COLUMN);
+
+	    if (elaboration_reg.west < dest_reg.west) {	/* First column */
+
+		P_set_regions(&elaboration_reg, &general_box, &overlap_box,
+			      dims, FIRST_COLUMN);
+		align_elaboration_box(&elaboration_reg, &src_reg, GENERAL_COLUMN);
+		align_interp_boxes(&general_box, &overlap_box, &dest_reg,
+			  last_general_box, last_overlap_box, FIRST_COLUMN);
+	    }
+
+	    if (elaboration_reg.east >= dest_reg.east) {	/* Last column */
+
+		P_set_regions(&elaboration_reg, &general_box, &overlap_box,
+			      dims, LAST_COLUMN);
+		last_column = TRUE;
+	    }
+	    nsplx =
+		ceil((elaboration_reg.east -
+		      elaboration_reg.west) / stepE) + 0.5;
+	    G_debug(1, "Interpolation: nsplx = %d", nsplx);
+
+	    if (grid_opt->answer) {
+		Vect_reset_cats(Cats);
+		Vect_cat_set(Cats, 1, cat++);
+		Vect_reset_line(Points);
+		Vect_append_point(Points, general_box.W, general_box.S, 0);
+		Vect_append_point(Points, general_box.E, general_box.S, 0);
+		Vect_append_point(Points, general_box.E, general_box.N, 0);
+		Vect_append_point(Points, general_box.W, general_box.N, 0);
+		Vect_append_point(Points, general_box.W, general_box.S, 0);
+		Vect_write_line(&Grid, GV_LINE, Points, Cats);
+		Vect_reset_line(Points);
+		Vect_append_point(Points, (general_box.E + general_box.W) / 2,
+					   (general_box.N + general_box.S) / 2, 0);
+		Vect_write_line(&Grid, GV_POINT, Points, Cats);
+	    }
+
+	    /* reading points in interpolation region */
+	    G_debug(1, "reading points from input raster...");
+	    dim_vect = nsplx * nsply;
+
+	    observ =
+		P_Read_Raster_Region_Map(inrast_matrix, &elaboration_reg,
+		                         &src_reg, &npoints, &n_nulls,
+					 dim_vect);
+
+	    G_debug(1, "%d points, %d NULL cells", npoints, n_nulls);
+
+	    G_debug(1,
+		    "Interpolation: (%d,%d): Number of points in <elaboration_box> is %d",
+		    subregion_row, subregion_col, npoints);
+
+	    /* Mean calculation for observed non-NULL points */
+	    if (npoints > 0)
+		mean = P_Mean_Calc(&elaboration_reg, observ, npoints);
+	    else
+		mean = 0;
+	    G_debug(1, "Interpolation: (%d,%d): mean=%lf",
+		    subregion_row, subregion_col, mean);
+
+	    observ_null = NULL;
+	    if (null_flag->answer) {
+		if (n_nulls) {
+		    /* read input NULL cells */
+
+		    G_debug(1, "read input NULL cells");
+
+		    observ_null =
+			P_Read_Raster_Region_Nulls(inrast_matrix, &src_reg,
+			                         dest_box, general_box,
+						 &npoints_null, 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;
+		    }
+		}
+	    }
+	    else
+		n_nulls = 1;
+
+	    if (npoints > 0 && n_nulls > 0) {	/*  */
+		int i;
+
+		nparameters = nsplx * nsply;
+		BW = P_get_BandWidth(P_BILINEAR, nsply) + 2 * (interp_method == P_BICUBIC);
+
+		/* Least Squares system */
+		N = G_alloc_matrix(nparameters, BW);	/* Normal matrix */
+		TN = G_alloc_vector(nparameters);	/* vector */
+		parVect = G_alloc_vector(nparameters);	/* Parameters vector */
+		obsVect = G_alloc_matrix(npoints, 3);	/* Observation vector */
+		Q = G_alloc_vector(npoints);	/* "a priori" var-cov matrix */
+
+		for (i = 0; i < npoints; i++) {	/* Setting obsVect vector & Q matrix */
+		    Q[i] = 1;	/* Q=I */
+		    obsVect[i][0] = observ[i].coordX;
+		    obsVect[i][1] = observ[i].coordY;
+		    obsVect[i][2] = observ[i].coordZ - mean;
+		}
+
+		G_free(observ);
+
+		/* Bilinear interpolation */
+		if (interp_method == P_BILINEAR) {
+		    G_debug(1,
+			    "Interpolation: (%d,%d): Bilinear interpolation...",
+			    subregion_row, subregion_col);
+		    normalDefBilin(N, TN, Q, obsVect, stepE, stepN, nsplx,
+				   nsply, elaboration_reg.west,
+				   elaboration_reg.south, npoints,
+				   nparameters, BW);
+		    nCorrectGrad(N, lambda, nsplx, nsply, stepE, stepN);
+		}
+		/* Bicubic interpolation */
+		else {
+		    G_debug(1,
+			    "Interpolation: (%d,%d): Bicubic interpolation...",
+			    subregion_row, subregion_col);
+		    normalDefBicubic(N, TN, Q, obsVect, stepE, stepN, nsplx,
+				     nsply, elaboration_reg.west,
+				     elaboration_reg.south, npoints,
+				     nparameters, BW);
+		    nCorrectGrad(N, lambda, nsplx, nsply, stepE, stepN);
+		}
+
+		G_math_solver_cholesky_sband(N, parVect, TN, nparameters, BW);
+
+		G_free_matrix(N);
+		G_free_vector(TN);
+		G_free_vector(Q);
+
+		if (!null_flag->answer) {	/* 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, parVect,
+					 stepN, stepE, dims.overlap, mean,
+					 nsplx, nsply, nrows, ncols, interp_method);
+		}
+		else {		/* only interpolate NULL cells */
+
+		    if (observ_null == NULL)
+			G_fatal_error("no NULL cells loaded");
+
+		    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,
+					stepE, stepN,
+					dims.overlap, nsplx, nsply,
+					npoints_null, interp_method, mean);
+
+		    G_free(observ_null);
+		}		/* end NULL cells */
+		G_free_vector(parVect);
+		G_free_matrix(obsVect);
+	    }
+	    else {
+		if (observ)
+		    G_free(observ);
+		if (npoints == 0)
+		    G_warning(_("No data within this subregion. "
+				"Consider increasing the spline step."));
+	    }
+	}			/*! END WHILE; last_column = TRUE */
+    }				/*! END WHILE; last_row = TRUE */
+
+    G_verbose_message(_("Writing output..."));
+    G_free_matrix(inrast_matrix);
+    /* Writing the output raster map */
+    Rast_set_fp_type(DCELL_TYPE);
+    outrastfd = Rast_open_fp_new(out_opt->answer);
+
+    /* check */
+    {
+	int nonulls = 0;
+	DCELL dval;
+
+	for (row = 0; row < dest_reg.rows; row++) {
+	    for (col = 0; col < dest_reg.cols; col++) {
+		dval = outrast_matrix[row][col];
+		if (!Rast_is_d_null_value(&dval))
+		    nonulls++;
+	    }
+	}
+	if (!nonulls)
+	    G_warning("only NULL cells in output raster");
+    }
+    P_Aux_to_Raster(outrast_matrix, outrastfd);
+    G_free_matrix(outrast_matrix);
+
+    Rast_close(outrastfd);
+
+    /* set map title */
+    sprintf(title, "%s interpolation with Tykhonov regularization",
+	    method_opt->answer);
+    Rast_put_cell_title(out_opt->answer, title);
+    /* write map history */
+    Rast_short_history(out_opt->answer, "raster", &history);
+    Rast_command_history(&history);
+    Rast_write_history(out_opt->answer, &history);
+
+    if (grid_opt->answer) {
+	Vect_build(&Grid);
+	Vect_close(&Grid);
+    }
+    
+    G_done_msg(" ");
+
+    exit(EXIT_SUCCESS);
+}				/*END MAIN */

vector/lidar/r.resamp.bspline/r.resamp.bspline.html

@@ -0,0 +1,123 @@
+<h2>DESCRIPTION</h2>
+<em>r.resamp.bspline</em> performs a bilinear/bicubic spline interpolation with
+Tykhonov regularization. The input is a raster surface map, e.g. elevation,
+temperature, precipitation etc. Output is a raster map. Optionally, only
+input NULL cells are interpolated, useful to fill NULL cells, an alternative
+to <a href="r.fillnulls.html">r.fillnulls</a>. Using the -n flag to only
+interpolate NULL cells will considerably speed up the module.
+<p>
+The input raster map is read at its native resolution, the output raster
+map will be produced for the current computational region set with
+<a href="g.region.html">g.region</a>. Any MASK will be respected, masked
+values will be treated as NULL cells in both the input and the output map.
+<p>
+Spline step values <b><i>se</i></b> for the east-west direction and
+<b><i>sn</i></b> for the north-south direction should not be smaller than
+the east-west and north-south resolutions of the input map. For very large
+areas with missing values (NULL cells), larger spline step values may be
+required, but most of the time the defaults (1.5 x resolution) should be fine.
+<p>
+The Tykhonov regularization parameter ("<b><i>lambda</i></b>") acts to
+smooth the interpolation. With a small <b><i>lambda</i></b>, the
+interpolated surface closely follows observation points; a larger value
+will produce a smoother interpolation. Reasonable values are 0.0001,
+0.001, 0.005, 0.01, 0.02, 0.05, 0.1. For seamless NULL cell interpolation,
+a small value is required and default is set to 0.001.
+<p>
+From a theoretical perspective, the interpolating procedure takes place in two
+parts: the first is an estimate of the linear coefficients of a spline function
+is derived from the observation points using a least squares regression; the
+second is the computation of the interpolated surface (or interpolated vector
+points). As used here, the splines are 2D piece-wise non-zero polynomial
+functions calculated within a limited, 2D area. The length of each spline step
+is defined by <b><i>se</i></b> for the east-west direction and
+<b><i>sn</i></b> for the north-south direction. For optimum performance, the
+spline step values should be no less than the east-west and north-south
+resolutions of the input map. Each non NULL cell observation is modeled as a
+linear function of the non-zero splines in the area around the observation.
+The least squares regression predicts the the coefficients of these linear functions.
+Regularization avoids the need to have one one observation and one coefficient
+for each spline (in order to avoid instability). 
+
+<p>
+A cross validation "leave-one-out" analysis is available to help to determine
+the optimal <b><i>lambda</i></b> value that produces an interpolation that
+best fits the original observation data. The more points used for
+cross-validation, the longer the time needed for computation. Empirical testing
+indicates a threshold of a maximum of 100 points is recommended. Note that cross
+validation can run very slowly if more than 100 observations are used. The
+cross-validation output reports <i>mean</i> and <i>rms</i> of the residuals from
+the true point value and the estimated from the interpolation for a fixed series
+of <b><i>lambda</i></b> values. No vector nor raster output will be created
+when cross-validation is selected. 
+
+<h2>EXAMPLES</h2>
+
+<h4>Basic interpolation</h4>
+
+<div class="code"><pre>
+r.resamp.bspline input=raster_surface output=interpolated_surface method=bicubic
+</pre></div>
+
+A bicubic spline interpolation will be done and a raster map with estimated
+(i.e., interpolated) values will be created. 
+
+<h4>Interpolation of NULL cells and patching</h4>
+
+<div class="code"><pre>
+# set region to area with NULL cells, align region to input map
+g.region n=north s=south e=east w=west align=input -p
+# interpolate NULL cells
+r.resamp.bspline -n input=input_raster output=interpolated_nulls method=bicubic
+# set region to area with NULL cells, align region to input map
+g.region rast=input -p
+# patch original map and interpolated NULLs
+r.patch input=input_raster,interpolated_nulls output=input_raster_gapfilled
+</pre></div>
+
+<h4>Estimation of <b><i>lambda</i></b> parameter with a cross validation proccess</h4>
+
+A random sample of points should be generated first with
+<a href="r.random.html">r.random</a>, and the current region should not
+include more than 100 non-NULL random cells. 
+
+<div class="code"><pre>
+r.resamp.bspline -c input=input_raster 
+</pre></div>
+
+<h2>SEE ALSO</h2>
+
+<em>
+<a href="r.fillnulls.html">r.fillnulls</a><br>
+<a href="r.resamp.rst.html">r.resamp.rst</a><br>
+<a href="r.resamp.interp.html">r.resamp.interp</a><br>
+<a href="v.surf.bspline.html">v.surf.bspline</a>
+</em>
+
+<h2>AUTHORS</h2>
+Markus Metz<br>
+<br>
+based on <a href="v.surf.bspline.html">v.surf.bspline</a> by
+<br>
+Maria Antonia Brovelli, Massimiliano Cannata, Ulisse Longoni, Mirko Reguzzoni, Roberto Antolin
+
+<h2>REFERENCES</h2>
+
+Brovelli M. A., Cannata M., and Longoni U.M., 2004, LIDAR Data
+Filtering and DTM Interpolation Within GRASS, Transactions in GIS,
+April 2004, vol. 8, iss. 2, pp. 155-174(20), Blackwell Publishing Ltd
+<p>
+Brovelli M. A. and Cannata M., 2004, Digital Terrain model
+reconstruction in urban areas from airborne laser scanning data: the
+method and an example for Pavia (Northern Italy). Computers and
+Geosciences 30, pp.325-331
+<p>
+Brovelli M. A e Longoni U.M., 2003, Software per il filtraggio di
+dati LIDAR, Rivista dell'Agenzia del Territorio, n. 3-2003, pp. 11-22
+(ISSN 1593-2192)
+<p>
+Antolin R. and Brovelli M.A., 2007, LiDAR data Filtering with GRASS GIS for the 
+Determination of Digital Terrain Models. Proceedings of Jornadas de SIG Libre,
+Girona, Espa&ntilde;a. CD ISBN: 978-84-690-3886-9 <br>
+
+<p><i>Last changed: $Date$</i>

vector/lidar/r.resamp.bspline/resamp.c

@@ -0,0 +1,401 @@
+
+/***********************************************************************
+ *
+ * MODULE:       r.resamp.bspline
+ *
+ * AUTHOR(S):    Markus Metz
+ *
+ * PURPOSE:      Spline Interpolation and cross correlation
+ *
+ * COPYRIGHT:    (C) 2010 GRASS development team
+ *
+ *               This program is free software under the
+ *               GNU General Public License (>=v2).
+ *               Read the file COPYING that comes with GRASS
+ *               for details.
+ *
+ **************************************************************************/
+
+#include "bspline.h"
+
+struct Point *P_Read_Raster_Region_Nulls(double **matrix,
+				       struct Cell_head *Original,
+				       struct bound_box output_box,
+				       struct bound_box General,
+				       int *num_points, int dim_vect,
+				       double mean)
+{
+    int col, row, startcol, endcol, startrow, endrow, nrows, ncols;
+    int pippo, npoints;
+    double X, Y, Z;
+    struct Point *obs;
+
+    pippo = dim_vect;
+    obs = (struct Point *)G_calloc(pippo, sizeof(struct Point));
+
+    /* Reading points inside input 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 */
+    if (Original->north > General.N) {
+	startrow = (double)((Original->north - General.N) / Original->ns_res - 1);
+	if (startrow < 0)
+	    startrow = 0;
+    }
+    else
+	startrow = 0;
+    if (Original->south < General.S) {
+	endrow = (double)((Original->north - General.S) / Original->ns_res + 1);
+	if (endrow > nrows)
+	    endrow = nrows;
+    }
+    else
+	endrow = nrows;
+    if (General.W > Original->west) {
+	startcol = (double)((General.W - Original->west) / Original->ew_res - 1);
+	if (startcol < 0)
+	    startcol = 0;
+    }
+    else
+	startcol = 0;
+    if (General.E < Original->east) {
+	endcol = (double)((General.E - Original->west) / Original->ew_res + 1);
+	if (endcol > ncols)
+	    endcol = ncols;
+    }
+    else
+	endcol = ncols;
+
+    for (row = startrow; row < endrow; row++) {
+	for (col = startcol; col < endcol; col++) {
+
+	    Z = matrix[row][col];
+
+	    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);
+
+		/* 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;
+
+		}
+	    }
+	}
+    }
+
+    *num_points = npoints;
+    return obs;
+}
+
+double **P_Sparse_Raster_Points(double **matrix, struct Cell_head *Elaboration,
+                struct Cell_head *Original, struct bound_box General,
+		struct bound_box Overlap, struct Point *obs,
+		double *param, double pe, double pn,
+		double overlap, int nsplx, int nsply, int num_points,
+		int bilin, double mean)
+{
+    int i, row, col;
+    double X, Y, interpolation, csi, eta, weight;
+    int points_in_box = 0;
+
+    /* Reading points inside output region and inside general box */
+    /* all points available here are inside the output box,
+     * selected by P_Read_Raster_Region_Nulls(), no check needed */
+
+    for (i = 0; i < num_points; i++) {
+
+	X = obs[i].coordX;
+	Y = obs[i].coordY;
+
+	/* Here mean is just for asking if obs point is in box */
+	if (Vect_point_in_box(X, Y, mean, &General)) { /* constraint_box */
+
+	    /* X,Y are cell center cordinates, MUST be inside General box */
+	    /* do NOT use floor, weird rounding errors, just truncate */
+	    row = (int) (floor(Rast_northing_to_row(Y, Original)) + 0.1);
+	    /* col = (int) floor(Rast_easting_to_col(X, Original)); */
+	    col = (int) (floor((X - Original->west) / Original->ew_res) + 0.1);
+
+	    if (row < 0 || row >= Original->rows) {
+		G_fatal_error("row index out of range");
+		continue;
+	    }
+
+	    if (col < 0 || col >= Original->cols) {
+		G_fatal_error("col index out of range");
+		continue;
+	    }
+	    points_in_box++;
+
+	    G_debug(3, "P_Sparse_Raster_Points: interpolate point %d...", i);
+	    if (bilin)
+		interpolation =
+		    dataInterpolateBilin(X, Y, pe, pn, nsplx,
+					 nsply, Elaboration->west,
+					 Elaboration->south, param);
+	    else
+		interpolation =
+		    dataInterpolateBicubic(X, Y, pe, pn, nsplx,
+					   nsply, Elaboration->west,
+					   Elaboration->south, param);
+
+	    interpolation += mean;
+
+	    if (Vect_point_in_box(X, Y, interpolation, &Overlap)) {	/* (5) */
+		matrix[row][col] = interpolation;
+	    }
+	    else {
+		if ((X > Overlap.E) && (X < General.E)) {
+		    if ((Y > Overlap.N) && (Y < General.N)) {	/* (3) */
+			csi = (General.E - X) / overlap;
+			eta = (General.N - Y) / overlap;
+			weight = csi * eta;
+			interpolation *= weight;
+			matrix[row][col] += interpolation;
+		    }
+		    else if ((Y < Overlap.S) && (Y > General.S)) {	/* (1) */
+			csi = (General.E - X) / overlap;
+			eta = (Y - General.S) / overlap;
+			weight = csi * eta;
+			interpolation *= weight;
+			matrix[row][col] = interpolation;
+		    }
+		    else if ((Y >= Overlap.S) && (Y <= Overlap.N)) {	/* (1) */
+			weight = (General.E - X ) / overlap;
+			interpolation *= weight;
+			matrix[row][col] = interpolation;
+		    }
+		}
+		else if ((X < Overlap.W) && (X > General.W)) {
+		    if ((Y > Overlap.N) && (Y < General.N)) {	/* (4) */
+			csi = (X - General.W) / overlap;
+			eta = (General.N - Y) / overlap;
+			weight = eta * csi;
+			interpolation *= weight;
+			matrix[row][col] += interpolation;
+		    }
+		    else if ((Y < Overlap.S) && (Y > General.S)) {	/* (2) */
+			csi = (X - General.W) / overlap;
+			eta = (Y - General.S) / overlap;
+			weight = csi * eta;
+			interpolation *= weight;
+			matrix[row][col] += interpolation;
+		    }
+		    else if ((Y >= Overlap.S) && (Y <= Overlap.N)) {	/* (2) */
+			weight = (X - General.W) / overlap;
+			interpolation *= weight;
+			matrix[row][col] += interpolation;
+		    }
+		}
+		else if ((X >= Overlap.W) && (X <= Overlap.E)) {
+		    if ((Y > Overlap.N) && (Y < General.N)) {	/* (3) */
+			weight = (General.N - Y) / overlap;
+			interpolation *= weight;
+			matrix[row][col] += interpolation;
+		    }
+		    else if ((Y < Overlap.S) && (Y > General.S)) {	/* (1) */
+			weight = (Y - General.S) / overlap;
+			interpolation *= weight;
+			matrix[row][col] = interpolation;
+		    }
+		}
+	    } /* end not in overlap */
+	}  /* end if in box */
+	else
+	    G_debug(0, "point outside General box ???");
+    }  /* for num_points */
+
+    return matrix;
+}
+
+/* align elaboration box to source region
+ * grow each side */
+int align_elaboration_box(struct Cell_head *elaboration, struct Cell_head *original, int type)
+{
+    int row, col;
+
+    switch (type) {
+    case GENERAL_ROW:		/* General case N-S direction */
+	/* northern edge */
+	row = (int)((original->north - elaboration->north) / original->ns_res);
+
+	if (row < 0)
+	    row = 0;
+
+	elaboration->north = original->north - original->ns_res * row;
+	
+	/* southern edge */
+	row = (int)((original->north - elaboration->south) / original->ns_res) + 1;
+
+	if (row > original->rows + 1)
+	    row = original->rows + 1;
+
+	elaboration->south = original->north - original->ns_res * row;
+
+	return 1;
+    
+    case GENERAL_COLUMN:	/* General case E-W direction */
+
+	/* eastern edge */
+	col = (int)((elaboration->east - original->west) / original->ew_res) + 1;
+
+	if (col > original->cols + 1)
+	    col = original->cols + 1;
+
+	elaboration->east = original->west + original->ew_res * col;
+	
+	/* western edge */
+	col = (int)((elaboration->west - original->west) / original->ew_res);
+
+	if (col < 0)
+	    col = 0;
+
+	elaboration->west = original->west + original->ew_res * col;
+
+	return 1;
+    }
+    return 0;
+}
+
+
+/* align interpolation boxes to destination region
+ * return 1 on success, 0 on failure */
+
+int align_interp_boxes(struct bound_box *general, struct bound_box *overlap,
+                       struct Cell_head *original, struct bound_box last_general,
+		       struct bound_box last_overlap, int type)
+{
+    int row, col;
+
+    switch (type) {
+    case GENERAL_ROW:		/* General case N-S direction */
+
+	/* general box */
+	/* grow north */
+	general->N = last_overlap.S;
+
+	/* shrink south */
+	row = (int)((original->north - general->S) / original->ns_res);
+
+	if (row > original->rows + 1)
+	    row = original->rows + 1;
+
+	general->S = original->north - original->ns_res * row;
+	
+	/* overlap box */
+	/* grow north */
+	overlap->N = last_general.S;
+
+	/* shrink south */
+	row = (int)((original->north - overlap->S) / original->ns_res);
+
+	if (row > original->rows + 1)
+	    row = original->rows + 1;
+
+	overlap->S = original->north - original->ns_res * row;
+
+	return 1;
+
+    case GENERAL_COLUMN:	/* General case E-W direction */
+
+	/* general box */
+	/* grow west */
+	general->W = last_overlap.E;
+
+	/* shrink east */
+	col = (int)((general->E - original->west) / original->ew_res);
+
+	if (col > original->cols + 1)
+	    col = original->cols + 1;
+
+	general->E = original->west + original->ew_res * col;
+	
+	/* overlap box */
+	/* grow west */
+	overlap->W = last_general.E;
+
+	/* shrink east */
+	col = (int)((overlap->E - original->west) / original->ew_res);
+
+	if (col > original->cols + 1)
+	    col = original->cols + 1;
+
+	overlap->E = original->west + original->ew_res * col;
+
+	return 1;
+
+    case FIRST_ROW:		/* Just started with first row */
+	general->N = original->north;
+	overlap->N = original->north;
+
+	/* shrink south */
+	row = (int)((original->north - general->S) / original->ns_res);
+
+	if (row > original->rows)
+	    row = original->rows;
+
+	general->S = original->north - original->ns_res * row;
+
+	row = (int)((original->north - overlap->S) / original->ns_res);
+
+	if (row > original->rows + 1)
+	    row = original->rows + 1;
+
+	overlap->S = original->north - original->ns_res * row;
+
+	return 1;
+
+    case LAST_ROW:		/* Reached last row */
+	general->S = original->south;
+	overlap->S = original->south;
+
+	return 1;
+
+    case FIRST_COLUMN:		/* Just started with first column */
+	general->W = original->west;
+	overlap->W = original->west;
+
+	/* shrink east */
+	col = (int)((general->E - original->west) / original->ew_res);
+
+	if (col > original->cols + 1)
+	    col = original->cols + 1;
+
+	general->E = original->west + original->ew_res * col;
+
+	col = (int)((overlap->E - original->west) / original->ew_res);
+
+	if (col > original->cols + 1)
+	    col = original->cols + 1;
+
+	overlap->E = original->west + original->ew_res * col;
+
+	return 1;
+
+    case LAST_COLUMN:		/* Reached last column */
+	general->E = original->east;
+	overlap->E = original->east;
+
+	return 1;
+    }
+
+    return 0;
+}