[GRASS-dev] Update of r.surf.idw to FP
Markus Neteler
neteler at itc.it
Mon Oct 15 15:51:42 EDT 2007
On Wed, Sep 05, 2007 at 02:54:57PM +0200, Markus Neteler wrote:
> Glynn Clements wrote on 03/20/2007 07:24 PM:
> > Markus Neteler wrote:
> >
> >> I have tried to update r.surf.idw to floating point:
> >>
> >> but I am not quite sure if this is efficient like this.
> >> Currently it doesn't compile due to some "extern" magic
> >> which I don't quite understand.
> >> I have been inspired by r.example.
> >>
> >
> > The problem is that you have renamed some (but not all) occurrences of
> > "cell" to "outrast", although the variable is still called "cell".
> >
> >
> (given that r.surf.idw should stay for more time)
>
> I have updated the patch (find attached) to the current CVS version,
> still crashing.
> Advice/fix most welcome...
Still struggling with the FP patch (find attached).
Markus
------------------
ITC -> dall'1 marzo 2007 Fondazione Bruno Kessler
ITC -> since 1 March 2007 Fondazione Bruno Kessler
------------------
-------------- next part --------------
Index: raster/r.surf.idw/main.c
===================================================================
RCS file: /grassrepository/grass6/raster/r.surf.idw/main.c,v
retrieving revision 2.13
diff -u -r2.13 main.c
--- raster/r.surf.idw/main.c 20 Jul 2007 22:33:57 -0000 2.13
+++ raster/r.surf.idw/main.c 15 Oct 2007 19:49:45 -0000
@@ -60,6 +60,8 @@
struct Flag *e;
} flag;
int n, fd, maskfd;
+ unsigned char *outrast;
+ RASTER_MAP_TYPE data_type;
/* Initialize the GIS calls */
G_gisinit(argv[0]);
@@ -119,8 +121,6 @@
/* initialize function pointers */
lookup_and_function_ptrs(nrows, ncols);
- /* allocate buffers for row i/o */
- cell = G_allocate_cell_buf();
if ((maskfd = G_maskfd()) >= 0 || error_flag) { /* apply mask to output */
if (error_flag) /* use input as mask when -e option chosen */
maskfd = G_open_cell_old(input, layer_mapset);
@@ -134,20 +134,24 @@
if (fd < 0)
G_fatal_error(_("Unable to open raster map <%s>"), input);
+ /* allocate buffers for row i/o */
+ data_type = G_get_raster_map_type(fd);
+ outrast = G_allocate_raster_buf(data_type);
+
/* Store input data in array-indexed doubly-linked lists and close input file */
- rowlist = row_lists(nrows, ncols, &datarows, &n, fd, cell);
+ rowlist = row_lists(nrows, ncols, &datarows, &n, fd, outrast, data_type);
G_close_cell(fd);
if (npoints > n)
npoints = n;
/* open cell layer for writing output */
- fd = G_open_cell_new(output);
+ fd = G_open_raster_new(output, data_type);
if (fd < 0)
G_fatal_error(_("Unable to create raster map <%s>"), output);
/* call the interpolation function */
- interpolate(rowlist, nrows, ncols, datarows, npoints, fd, maskfd);
+ interpolate(rowlist, nrows, ncols, datarows, npoints, fd, maskfd, data_type);
/* free allocated memory */
free_row_lists(rowlist, nrows);
@@ -209,10 +213,9 @@
int
interpolate(MELEMENT rowlist[], SHORT nrows, SHORT ncols, SHORT datarows,
- int npoints, int out_fd, int maskfd)
+ int npoints, int out_fd, int maskfd, RASTER_MAP_TYPE data_type)
{
- extern CELL *cell;
-
+ extern unsigned char *outrast; /* output buffer */
MELEMENT *Rptr;
EW *search, *ewptr, *current_row, /* start row for north/south search */
*lastrow; /* last element in search array */
@@ -236,7 +239,7 @@
G_percent(row, nrows, 2);
/* if mask occurs, read current row of the mask */
- if (mask && G_get_map_row(maskfd, mask, row) < 0)
+ if (mask && G_get_raster_row(maskfd, mask, row, data_type) < 0)
G_fatal_error(_("Cannot read row"));
/* prepare search array for next row of interpolations */
@@ -248,14 +251,25 @@
/* if (row != 279 && col != 209) continue; */
/* don't interpolate outside of the mask */
- if (mask && mask[col] == 0) {
- cell[col] = 0;
- continue;
- }
+ if (mask && mask[col] == 0) {
+ switch (data_type) {
+ case CELL_TYPE:
+ ((CELL *) outrast)[col] = 0;
+ break;
+ case FCELL_TYPE:
+ ((FCELL *) outrast)[col] = 0.;
+ break;
+ case DCELL_TYPE:
+ ((DCELL *) outrast)[col] = 0.;
+ break;
+ }
+ continue;
+ }
/* make a list of npoints neighboring data pts */
nbr_head->next = NULL;
- if (make_neighbors_list(search, lastrow, current_row, row, col, nbr_head, npoints)) { /* otherwise, known data value assigned */
+ if (make_neighbors_list(search, lastrow, current_row, row, col, nbr_head, npoints, data_type)) {
+ /* otherwise, known data value assigned */
/* calculate value to be set for the cell from the data values
* of npoints closest neighboring points */
@@ -268,15 +282,36 @@
Nptr = Nptr->next;
} while (Nptr); /* to end of list */
- cell[col] = (CELL) (sum1 / sum2 + .5);
- /* fprintf (stdout,"%d,%d = %d\n", col, row, cell[col]); */
-
- if (error_flag) /* output interpolation error for this cell */
- cell[col] -= mask[col];
- }
+ switch (data_type) {
+ case CELL_TYPE:
+ ((CELL *) outrast)[col] = (CELL) (sum1 / sum2 + .5);
+ break;
+ case FCELL_TYPE:
+ ((FCELL *) outrast)[col] = sum1 / sum2 + .5;
+ break;
+ case DCELL_TYPE:
+ ((DCELL *) outrast)[col] = sum1 / sum2 + .5;
+ break;
+ }
+
+ /* fprintf (stdout,"%d,%d = %d\n", col, row, cell[col]); */
+
+ if (error_flag) /* output interpolation error for this cell */
+ switch (data_type) {
+ case CELL_TYPE:
+ ((CELL *) outrast)[col] -= mask[col];
+ break;
+ case FCELL_TYPE:
+ ((FCELL *) outrast)[col] -= mask[col];
+ break;
+ case DCELL_TYPE:
+ ((DCELL *) outrast)[col] -= mask[col];
+ break;
+ }
+ }
} /* end of loop over columns */
- G_put_raster_row(out_fd, cell, CELL_TYPE);
+ G_put_raster_row(out_fd, outrast, data_type);
/* advance current row pointer if necessary */
if (current_row->start->y == row && current_row != lastrow)
@@ -299,10 +334,10 @@
/* to be interpolated using data value of its neighbors */
int make_neighbors_list(EW * firstrow, EW * lastrow, EW * curr_row, SHORT row, SHORT col, NEIGHBOR * head, /* head points to dummy plus npoints neighbors */
- int npoints)
+ int npoints, RASTER_MAP_TYPE data_type)
{
- extern CELL *cell;
+ extern unsigned char *outrast; /* output buffer */
SHORT neighbors = 0, /* number of neighbors in current list */
nsearch = 1, ssearch = 1; /* expand search north and south */
EW *north, *south;
@@ -320,9 +355,19 @@
else
north->east = north->east->next;
}
- else { /* no interpolation required */
- cell[col] = north->east->value;
- return (0);
+ else { /* no interpolation required */
+ switch (data_type) {
+ case CELL_TYPE:
+ ((CELL *) outrast)[col] = north->east->value;
+ break;
+ case FCELL_TYPE:
+ ((FCELL *) outrast)[col] = north->east->value;
+ break;
+ case DCELL_TYPE:
+ ((DCELL *) outrast)[col] = north->east->value;
+ break;
+ }
+ return (0);
}
}
/* initialize south search routine */
@@ -664,7 +709,8 @@
SHORT * datarows, /* number of rows with non-zero input data */
int *npts, /* number of data points available */
int fd, /* file descriptor, input */
- CELL * cell /* array of data for a single row */
+ unsigned char *cell, /* array of data for a single row */
+ RASTER_MAP_TYPE data_type
)
{
int row, col; /* row and column indices */
@@ -684,7 +730,7 @@
for (row = 0, Rptr = rowlist; row < rows; row++) {
G_percent(row, rows, 1);
- if (G_get_map_row_nomask(fd, cell, row) < 0)
+ if (G_get_raster_row_nomask(fd, cell, row, data_type) < 0)
G_fatal_error(_("Cannot read row"));
for (col = 0; col < cols; col++) {
@@ -744,6 +790,4 @@
*nextcol = (double)i *i * ew2;
return 0;
}
-
-
Index: raster/r.surf.idw/main.h
===================================================================
RCS file: /grassrepository/grass6/raster/r.surf.idw/main.h,v
retrieving revision 2.0
diff -u -r2.0 main.h
--- raster/r.surf.idw/main.h 9 Nov 2004 14:02:05 -0000 2.0
+++ raster/r.surf.idw/main.h 15 Oct 2007 19:49:45 -0000
@@ -67,6 +67,7 @@
double offset_distance_LL(SHORT);
double (*check_offset) (SHORT); /* function pointer */
+unsigned char *outrast;
#endif
/* dist.c */
@@ -83,8 +84,8 @@
int LL_lookup_tables(SHORT, SHORT);
/* main.c */
int lookup_and_function_ptrs(SHORT, SHORT);
-int interpolate(MELEMENT [], SHORT, SHORT, SHORT, int, int, int);
-int make_neighbors_list(EW *, EW *, EW *, SHORT, SHORT, NEIGHBOR *, int);
+int interpolate(MELEMENT [], SHORT, SHORT, SHORT, int, int, int, RASTER_MAP_TYPE);
+int make_neighbors_list(EW *, EW *, EW *, SHORT, SHORT, NEIGHBOR *, int, RASTER_MAP_TYPE);
int search(EW **, NEIGHBOR *, SHORT, SHORT, int, SHORT *, EW *, SHORT);
int exhaust(EW **, NEIGHBOR *, SHORT, SHORT);
EW *next_row(EW *, EW *, SHORT *, SHORT);
@@ -93,5 +94,5 @@
int replace_neighbor(MELEMENT **, NEIGHBOR *, double);
int sort_neighbors(NEIGHBOR *, double);
int free_row_lists(MELEMENT *, SHORT);
-MELEMENT *row_lists(SHORT, SHORT, SHORT *, int *, int, CELL *);
+MELEMENT *row_lists(SHORT, SHORT, SHORT *, int *, int, unsigned char *, RASTER_MAP_TYPE);
int lookup_tables(SHORT, SHORT);
More information about the grass-dev
mailing list