[GRASS-SVN] r42139 - in grass/trunk/vector/lidar: . r.resamp.bspline
svn_grass at osgeo.org
svn_grass at osgeo.org
Fri May 7 04:03:11 EDT 2010
Author: mmetz
Date: 2010-05-07 04:03:03 -0400 (Fri, 07 May 2010)
New Revision: 42139
Added:
grass/trunk/vector/lidar/r.resamp.bspline/
grass/trunk/vector/lidar/r.resamp.bspline/Makefile
grass/trunk/vector/lidar/r.resamp.bspline/bspline.h
grass/trunk/vector/lidar/r.resamp.bspline/crosscorr.c
grass/trunk/vector/lidar/r.resamp.bspline/main.c
grass/trunk/vector/lidar/r.resamp.bspline/r.resamp.bspline.html
grass/trunk/vector/lidar/r.resamp.bspline/resamp.c
Log:
new module for raster resampling and NULL cell filling
Added: grass/trunk/vector/lidar/r.resamp.bspline/Makefile
===================================================================
--- grass/trunk/vector/lidar/r.resamp.bspline/Makefile (rev 0)
+++ grass/trunk/vector/lidar/r.resamp.bspline/Makefile 2010-05-07 08:03:03 UTC (rev 42139)
@@ -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
+
Added: grass/trunk/vector/lidar/r.resamp.bspline/bspline.h
===================================================================
--- grass/trunk/vector/lidar/r.resamp.bspline/bspline.h (rev 0)
+++ grass/trunk/vector/lidar/r.resamp.bspline/bspline.h 2010-05-07 08:03:03 UTC (rev 42139)
@@ -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);
Added: grass/trunk/vector/lidar/r.resamp.bspline/crosscorr.c
===================================================================
--- grass/trunk/vector/lidar/r.resamp.bspline/crosscorr.c (rev 0)
+++ grass/trunk/vector/lidar/r.resamp.bspline/crosscorr.c 2010-05-07 08:03:03 UTC (rev 42139)
@@ -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(®ion);
+
+
+ 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, ®ion, 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
Added: grass/trunk/vector/lidar/r.resamp.bspline/main.c
===================================================================
--- grass/trunk/vector/lidar/r.resamp.bspline/main.c (rev 0)
+++ grass/trunk/vector/lidar/r.resamp.bspline/main.c 2010-05-07 08:03:03 UTC (rev 42139)
@@ -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 */
Added: grass/trunk/vector/lidar/r.resamp.bspline/r.resamp.bspline.html
===================================================================
--- grass/trunk/vector/lidar/r.resamp.bspline/r.resamp.bspline.html (rev 0)
+++ grass/trunk/vector/lidar/r.resamp.bspline/r.resamp.bspline.html 2010-05-07 08:03:03 UTC (rev 42139)
@@ -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ña. CD ISBN: 978-84-690-3886-9 <br>
+
+<p><i>Last changed: $Date$</i>
Property changes on: grass/trunk/vector/lidar/r.resamp.bspline/r.resamp.bspline.html
___________________________________________________________________
Added: svn:keywords
+ Date
Added: grass/trunk/vector/lidar/r.resamp.bspline/resamp.c
===================================================================
--- grass/trunk/vector/lidar/r.resamp.bspline/resamp.c (rev 0)
+++ grass/trunk/vector/lidar/r.resamp.bspline/resamp.c 2010-05-07 08:03:03 UTC (rev 42139)
@@ -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;
+}
More information about the grass-commit
mailing list