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 ------------------
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);
_______________________________________________ grass-dev mailing list grass-dev@grass.itc.it http://grass.itc.it/mailman/listinfo/grass-dev