[GRASS-SVN] r46267 - grass/trunk/raster/r.resamp.bspline
svn_grass at osgeo.org
svn_grass at osgeo.org
Sat May 14 09:11:03 EDT 2011
Author: mmetz
Date: 2011-05-14 06:11:03 -0700 (Sat, 14 May 2011)
New Revision: 46267
Modified:
grass/trunk/raster/r.resamp.bspline/Makefile
grass/trunk/raster/r.resamp.bspline/bspline.h
grass/trunk/raster/r.resamp.bspline/crosscorr.c
grass/trunk/raster/r.resamp.bspline/main.c
grass/trunk/raster/r.resamp.bspline/resamp.c
Log:
lidarlib uses segmentlib
Modified: grass/trunk/raster/r.resamp.bspline/Makefile
===================================================================
--- grass/trunk/raster/r.resamp.bspline/Makefile 2011-05-14 13:09:50 UTC (rev 46266)
+++ grass/trunk/raster/r.resamp.bspline/Makefile 2011-05-14 13:11:03 UTC (rev 46267)
@@ -2,8 +2,8 @@
PGM = r.resamp.bspline
-LIBES = $(LIDARLIB) $(GMATHLIB) $(VECTORLIB) $(DBMILIB) $(RASTERLIB) $(GISLIB) $(MATHLIB)
-DEPENDENCIES = $(LIDARDEP) $(GMATHDEP) $(VECTORDEP) $(DBMIDEP) $(RASTERDEP) $(GISDEP)
+LIBES = $(LIDARLIB) $(GMATHLIB) $(VECTORLIB) $(DBMILIB) $(RASTERLIB) $(SEGMENTLIB) $(GISLIB) $(MATHLIB)
+DEPENDENCIES = $(LIDARDEP) $(GMATHDEP) $(VECTORDEP) $(DBMIDEP) $(RASTERDEP) $(SEGMENTDEP) $(GISDEP)
EXTRA_INC = $(VECT_INC)
EXTRA_CFLAGS = $(VECT_CFLAGS)
Modified: grass/trunk/raster/r.resamp.bspline/bspline.h
===================================================================
--- grass/trunk/raster/r.resamp.bspline/bspline.h 2011-05-14 13:09:50 UTC (rev 46266)
+++ grass/trunk/raster/r.resamp.bspline/bspline.h 2011-05-14 13:11:03 UTC (rev 46267)
@@ -40,12 +40,12 @@
};
/* resamp.c */
-struct Point *P_Read_Raster_Region_masked(char **, /**/
+struct Point *P_Read_Raster_Region_masked(SEGMENT *, /**/
struct Cell_head *, /**/
struct bound_box, /**/
struct bound_box, /**/
int *, /**/ int, /**/ double);
-double **P_Sparse_Raster_Points(double **, /**/
+int P_Sparse_Raster_Points(SEGMENT *, /**/
struct Cell_head *, /**/
struct Cell_head *, /**/
struct bound_box, /**/
@@ -66,4 +66,4 @@
/* crosscor.c */
-int cross_correlation(double **, struct Cell_head *, double, double);
+int cross_correlation(SEGMENT *, struct Cell_head *, double, double);
Modified: grass/trunk/raster/r.resamp.bspline/crosscorr.c
===================================================================
--- grass/trunk/raster/r.resamp.bspline/crosscorr.c 2011-05-14 13:09:50 UTC (rev 46266)
+++ grass/trunk/raster/r.resamp.bspline/crosscorr.c 2011-05-14 13:11:03 UTC (rev 46267)
@@ -117,7 +117,7 @@
}
/*-------------------------------------------------------------------------------------------*/
-int cross_correlation(double **matrix, struct Cell_head *src_reg, double passWE, double passNS)
+int cross_correlation(SEGMENT *in_seg, struct Cell_head *src_reg, double passWE, double passNS)
/*
matrix: matrix with raster values
passWE: spline step in West-East direction
@@ -164,7 +164,7 @@
/*Cats = Vect_new_cats_struct (); */
/* Current region is read and points recorded into observ */
- observ = P_Read_Raster_Region_Map(matrix, ®ion, src_reg, &ndata, 1024);
+ observ = P_Read_Raster_Region_Map(in_seg, ®ion, src_reg, &ndata, 1024);
G_debug(5, "CrossCorrelation: %d points read in region. ", ndata);
G_verbose_message(_("%d points read in region"),
ndata);
Modified: grass/trunk/raster/r.resamp.bspline/main.c
===================================================================
--- grass/trunk/raster/r.resamp.bspline/main.c 2011-05-14 13:09:50 UTC (rev 46266)
+++ grass/trunk/raster/r.resamp.bspline/main.c 2011-05-14 13:11:03 UTC (rev 46267)
@@ -17,12 +17,16 @@
**********************************************************************/
/* INCLUDES */
-#include <grass/config.h>
#include <stdlib.h>
#include <string.h>
+#include <sys/types.h>
+#include <sys/stat.h>
+#include <fcntl.h>
#include <math.h>
#include "bspline.h"
+#define SEGSIZE 64
+
/*--------------------------------------------------------------------*/
int main(int argc, char *argv[])
{
@@ -39,14 +43,18 @@
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 */
- char **mask_matrix; /* matrix for masking */
- /* Structs declarations */
+ SEGMENT in_seg, out_seg, mask_seg;
+ const char *in_file, *out_file, *mask_file;
+ int in_fd, out_fd, mask_fd;
+ double seg_size;
+ int seg_mb, segments_in_memory;
+ int have_mask;
+
int inrastfd, outrastfd;
- DCELL *drastbuf;
+ DCELL *drastbuf, dval;
struct History history;
struct Map_info Grid;
@@ -56,7 +64,7 @@
struct GModule *module;
struct Option *in_opt, *out_opt, *grid_opt, *stepE_opt, *stepN_opt,
- *lambda_f_opt, *method_opt, *mask_opt;
+ *lambda_f_opt, *method_opt, *mask_opt, *memory_opt;
struct Flag *null_flag, *cross_corr_flag;
struct Reg_dimens dims;
@@ -65,7 +73,7 @@
struct bound_box last_overlap_box, last_general_box;
struct Point *observ;
- struct Point *observ_masked;
+ struct Point *observ_marked;
/*----------------------------------------------------------------*/
/* Options declarations */
@@ -135,6 +143,13 @@
cross_corr_flag->description =
_("Find the best Tykhonov regularizing parameter using a \"leave-one-out\" cross validation method");
+ memory_opt = G_define_option();
+ memory_opt->key = "memory";
+ memory_opt->type = TYPE_INTEGER;
+ memory_opt->required = NO;
+ memory_opt->answer = "300";
+ memory_opt->description = _("Maximum memory to be used (in MB)");
+
/*----------------------------------------------------------------*/
/* Parsing */
G_gisinit(argv[0]);
@@ -194,6 +209,10 @@
nsplx_adj = NSPLX_MAX;
nsply_adj = NSPLY_MAX;
+ if (interp_method == P_BICUBIC) {
+ nsplx_adj = 100;
+ nsply_adj = 100;
+ }
dims.overlap = OVERLAP_SIZE * (stepN > stepE ? stepN : stepE);
P_get_edge(interp_method, &dims, stepE, stepN);
@@ -262,21 +281,41 @@
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)"));
+ seg_mb = atoi(memory_opt->answer);
+ if (seg_mb < 3)
+ G_fatal_error(_("Memory in MB must be >= 3"));
+ if (mask_opt->answer || null_flag->answer) {
+ seg_size = sizeof(double) * 2 + sizeof(char);
+ }
+ else {
+ seg_size = sizeof(double) * 2;
+ }
+ seg_size = (seg_size * SEGSIZE * SEGSIZE) / (1 << 20);
+ segments_in_memory = seg_mb / seg_size + 0.5;
+ if (segments_in_memory < 1)
+ segments_in_memory = 1;
+
+ in_file = G_tempfile();
+ in_fd = creat(in_file, 0666);
+ if (segment_format(in_fd, nrows, ncols, SEGSIZE, SEGSIZE, sizeof(double)) != 1)
+ G_fatal_error(_("Can not create temporary file"));
+ close(in_fd);
+
+ in_fd = open(in_file, 2);
+ if (segment_init(&in_seg, in_fd, segments_in_memory) != 1)
+ G_fatal_error(_("Can not initialize temporary file"));
+
+ /* read raster input */
+ inrastfd = Rast_open_old(in_opt->answer, "");
drastbuf = Rast_allocate_buf(DCELL_TYPE);
- G_message("loading input raster <%s>", in_opt->answer);
+ G_message(_("Loading input raster <%s>"), in_opt->answer);
if (1) {
int got_one = 0;
for (row = 0; row < nrows; row++) {
@@ -287,15 +326,16 @@
Rast_get_d_row_nomask(inrastfd, drastbuf, row);
for (col = 0; col < ncols; col++) {
- inrast_matrix[row][col] = drastbuf[col];
- dval = inrast_matrix[row][col];
+ dval = drastbuf[col];
if (!Rast_is_d_null_value(&dval)) {
got_one++;
}
}
+ segment_put_row(&in_seg, drastbuf, row);
+
}
if (!got_one)
- G_fatal_error("only NULL cells in input raster");
+ G_fatal_error(_("Only NULL cells in input raster"));
}
G_percent(row, nrows, 2);
G_free(drastbuf);
@@ -318,36 +358,45 @@
if (cross_corr_flag->answer) {
G_debug(1, "CrossCorrelation()");
- if (cross_correlation(inrast_matrix, &src_reg, stepE, stepN) != TRUE)
+ if (cross_correlation(&in_seg, &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);
- G_done_msg(_("Cross validation finished for se = %f and sn = %f"), stepE, stepN);
+ segment_release(&in_seg); /* release memory */
+ close(in_fd);
+ unlink(in_file);
+
exit(EXIT_SUCCESS);
}
}
/* Alloc and load masking matrix */
/* encoding: 0 = do not interpolate, 1 = interpolate */
- if (mask_opt->answer || null_flag->answer) {
+ have_mask = mask_opt->answer || null_flag->answer;
+ if (have_mask) {
int maskfd;
int null_count = 0;
- DCELL dval;
CELL cval;
CELL *maskbuf;
+ char mask_val;
G_message(_("Mark cells for interpolation"));
/* use destination window */
- mask_matrix = (char **)G_calloc(nrows, sizeof(char *));
- mask_matrix[0] = (char *)G_calloc(nrows * ncols, sizeof(char));
- for (row = 1; row < nrows; row++)
- mask_matrix[row] = mask_matrix[row - 1] + ncols;
+ mask_file = G_tempfile();
+ mask_fd = creat(mask_file, 0666);
+ if (segment_format(mask_fd, nrows, ncols, SEGSIZE, SEGSIZE, sizeof(char)) != 1)
+ G_fatal_error(_("Can not create temporary file"));
+ close(mask_fd);
+ mask_fd = open(mask_file, 2);
+ if (segment_init(&mask_seg, mask_fd, segments_in_memory) != 1)
+ G_fatal_error(_("Can not initialize temporary file"));
+
if (mask_opt->answer) {
maskfd = Rast_open_old(mask_opt->answer, "");
maskbuf = Rast_allocate_buf(CELL_TYPE);
@@ -376,21 +425,22 @@
Rast_get_d_row(inrastfd, drastbuf, row);
for (col = 0; col < ncols; col++) {
- mask_matrix[row][col] = 1;
+ mask_val = 1;
if (mask_opt->answer) {
cval = maskbuf[col];
if (Rast_is_c_null_value(&cval) || cval == 0)
- mask_matrix[row][col] = 0;
+ mask_val = 0;
}
- if (null_flag->answer && mask_matrix[row][col] == 1) {
+ if (null_flag->answer && mask_val == 1) {
dval = drastbuf[col];
if (!Rast_is_d_null_value(&dval))
- mask_matrix[row][col] = 0;
+ mask_val = 0;
else
null_count++;
}
+ segment_put(&mask_seg, &mask_val, row, col);
}
}
@@ -407,26 +457,25 @@
G_fatal_error(_("No NULL cells found in input raster."));
}
}
- else
- mask_matrix = NULL;
-
- /* Alloc raster matrix */
- if (!(outrast_matrix = G_alloc_matrix(nrows, ncols)))
- G_fatal_error(_("Cannot allocate memory for auxiliar matrix."
- "Consider changing region (resolution)"));
+ out_file = G_tempfile();
+ out_fd = creat(out_file, 0666);
+ if (segment_format(out_fd, nrows, ncols, SEGSIZE, SEGSIZE, sizeof(double)) != 1)
+ G_fatal_error(_("Can not create temporary file"));
+ close(out_fd);
+
+ out_fd = open(out_file, 2);
+ if (segment_init(&out_seg, out_fd, segments_in_memory) != 1)
+ G_fatal_error(_("Can not initialize temporary file"));
+
/* 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;
- }
- }
+ drastbuf = Rast_allocate_buf(DCELL_TYPE);
+ Rast_set_d_null_value(drastbuf, ncols);
+ for (row = 0; row < nrows; row++) {
+ G_percent(row, nrows, 2);
+ segment_put_row(&out_seg, drastbuf, row);
}
G_percent(row, nrows, 2);
@@ -548,7 +597,7 @@
dim_vect = nsplx * nsply;
observ =
- P_Read_Raster_Region_Map(inrast_matrix, &elaboration_reg,
+ P_Read_Raster_Region_Map(&in_seg, &elaboration_reg,
&src_reg, &npoints, dim_vect);
G_debug(1, "%d valid points", npoints);
@@ -565,22 +614,22 @@
G_debug(1, "Interpolation: (%d,%d): mean=%lf",
subregion_row, subregion_col, mean);
- observ_masked = NULL;
+ observ_marked = NULL;
- if (mask_matrix) {
+ if (have_mask) {
/* collect unmasked output cells */
G_debug(1, "collect unmasked output cells");
- observ_masked =
- P_Read_Raster_Region_masked(mask_matrix, &dest_reg,
+ observ_marked =
+ P_Read_Raster_Region_masked(&mask_seg, &dest_reg,
dest_box, general_box,
&npoints_marked, dim_vect, mean);
G_debug(1, "%d cells marked in general", npoints_marked);
if (npoints_marked == 0) {
- G_free(observ_masked);
- observ_masked = NULL;
+ G_free(observ_marked);
+ observ_marked = NULL;
npoints = 1; /* disable warning below */
}
}
@@ -636,30 +685,29 @@
G_free_vector(TN);
G_free_vector(Q);
- if (!observ_masked) { /* interpolate full output raster */
+ if (!observ_marked) { /* 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, NULL,
- parVect, stepN, stepE, dims.overlap, mean,
- nsplx, nsply, nrows, ncols, interp_method);
+
+ P_Regular_Points(&elaboration_reg, &dest_reg, general_box,
+ overlap_box, &out_seg,
+ parVect, stepN, stepE, dims.overlap, mean,
+ nsplx, nsply, nrows, ncols, interp_method);
}
else { /* only interpolate selected cells */
G_debug(1, "Interpolation of %d selected cells...",
npoints_marked);
- outrast_matrix =
- P_Sparse_Raster_Points(outrast_matrix,
- &elaboration_reg, &dest_reg,
- general_box, overlap_box,
- observ_masked, parVect,
- stepE, stepN,
- dims.overlap, nsplx, nsply,
- npoints_marked, interp_method, mean);
+ P_Sparse_Raster_Points(&out_seg,
+ &elaboration_reg, &dest_reg,
+ general_box, overlap_box,
+ observ_marked, parVect,
+ stepE, stepN,
+ dims.overlap, nsplx, nsply,
+ npoints_marked, interp_method, mean);
- G_free(observ_masked);
+ G_free(observ_marked);
} /* end NULL cells */
G_free_vector(parVect);
G_free_matrix(obsVect);
@@ -674,12 +722,17 @@
} /*! END WHILE; last_column = TRUE */
} /*! END WHILE; last_row = TRUE */
+ segment_release(&in_seg); /* release memory */
+ close(in_fd);
+ unlink(in_file);
+
+ if (have_mask) {
+ segment_release(&mask_seg); /* release memory */
+ close(mask_fd);
+ unlink(mask_file);
+ }
+
G_verbose_message(_("Writing output..."));
- G_free_matrix(inrast_matrix);
- if (mask_opt->answer) {
- G_free(mask_matrix[0]);
- G_free(mask_matrix);
- }
/* Writing the output raster map */
Rast_set_fp_type(DCELL_TYPE);
outrastfd = Rast_open_fp_new(out_opt->answer);
@@ -687,21 +740,29 @@
/* check */
{
int nonulls = 0;
- DCELL dval;
+ segment_flush(&out_seg);
+ drastbuf = Rast_allocate_d_buf();
+
for (row = 0; row < dest_reg.rows; row++) {
+ G_percent(row, dest_reg.rows, 2);
+ segment_get_row(&out_seg, drastbuf, row);
for (col = 0; col < dest_reg.cols; col++) {
- dval = outrast_matrix[row][col];
+ dval = drastbuf[col];
if (!Rast_is_d_null_value(&dval))
nonulls++;
}
+ Rast_put_d_row(outrastfd, drastbuf);
}
+ G_percent(1, 1, 1);
if (!nonulls)
G_warning("only NULL cells in output raster");
}
- P_Aux_to_Raster(outrast_matrix, outrastfd);
- G_free_matrix(outrast_matrix);
+ segment_release(&out_seg); /* release memory */
+ close(out_fd);
+ unlink(out_file);
+
Rast_close(outrastfd);
/* set map title */
Modified: grass/trunk/raster/r.resamp.bspline/resamp.c
===================================================================
--- grass/trunk/raster/r.resamp.bspline/resamp.c 2011-05-14 13:09:50 UTC (rev 46266)
+++ grass/trunk/raster/r.resamp.bspline/resamp.c 2011-05-14 13:11:03 UTC (rev 46267)
@@ -22,7 +22,7 @@
#include <math.h>
#include "bspline.h"
-struct Point *P_Read_Raster_Region_masked(char **mask_matrix,
+struct Point *P_Read_Raster_Region_masked(SEGMENT *mask_seg,
struct Cell_head *Original,
struct bound_box output_box,
struct bound_box General,
@@ -33,6 +33,7 @@
int pippo, npoints;
double X, Y;
struct Point *obs;
+ char mask_val;
pippo = dim_vect;
obs = (struct Point *)G_calloc(pippo, sizeof(struct Point));
@@ -77,7 +78,8 @@
for (row = startrow; row < endrow; row++) {
for (col = startcol; col < endcol; col++) {
- if (!mask_matrix[row][col])
+ segment_get(mask_seg, &mask_val, row, col);
+ if (!mask_val)
continue;
X = Rast_col_to_easting((double)(col) + 0.5, Original);
@@ -106,14 +108,14 @@
return obs;
}
-double **P_Sparse_Raster_Points(double **matrix, struct Cell_head *Elaboration,
+int P_Sparse_Raster_Points(SEGMENT *out_seg, 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;
+ double X, Y, interpolation, csi, eta, weight, dval;
int points_in_box = 0;
/* Reading points inside output region and inside general box */
@@ -125,105 +127,99 @@
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 */
+ row = (int) (floor(Rast_northing_to_row(Y, Original)) + 0.1);
+ col = (int) (floor((X - Original->west) / Original->ew_res) + 0.1);
- /* 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 (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++;
- 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);
- 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;
- interpolation += mean;
-
- if (Vect_point_in_box(X, Y, interpolation, &Overlap)) { /* (5) */
- matrix[row][col] = interpolation;
+ if (Vect_point_in_box(X, Y, interpolation, &Overlap)) { /* (5) */
+ dval = interpolation;
+ }
+ else {
+ segment_get(out_seg, &dval, row, col);
+ 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;
+ dval += 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;
+ dval = interpolation;
+ }
+ else if ((Y >= Overlap.S) && (Y <= Overlap.N)) { /* (1) */
+ weight = (General.E - X ) / overlap;
+ interpolation *= weight;
+ dval = 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;
+ dval += 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 ((Y < Overlap.S) && (Y > General.S)) { /* (2) */
+ csi = (X - General.W) / overlap;
+ eta = (Y - General.S) / overlap;
+ weight = csi * eta;
+ interpolation *= weight;
+ dval += 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;
- }
+ else if ((Y >= Overlap.S) && (Y <= Overlap.N)) { /* (2) */
+ weight = (X - General.W) / overlap;
+ interpolation *= weight;
+ dval += interpolation;
}
- } /* end not in overlap */
- } /* end if in box */
- else
- G_debug(0, "point outside General box ???");
+ }
+ else if ((X >= Overlap.W) && (X <= Overlap.E)) {
+ if ((Y > Overlap.N) && (Y < General.N)) { /* (3) */
+ weight = (General.N - Y) / overlap;
+ interpolation *= weight;
+ dval += interpolation;
+ }
+ else if ((Y < Overlap.S) && (Y > General.S)) { /* (1) */
+ weight = (Y - General.S) / overlap;
+ interpolation *= weight;
+ dval = interpolation;
+ }
+ }
+ } /* end not in overlap */
+ segment_put(out_seg, &dval, row, col);
} /* for num_points */
- return matrix;
+ return 1;
}
/* align elaboration box to source region
More information about the grass-commit
mailing list