sample.c

Go to the documentation of this file.
00001 /*
00002  * Copyright (C) 1994. James Darrell McCauley.  (darrell@mccauley-usa.com)
00003  *                                              http://mccauley-usa.com/
00004  *
00005  * This program is free software under the GPL (>=v2)
00006  * Read the file GPL.TXT coming with GRASS for details.
00007  *
00008  * 1/2006: moved to libgis from v.sample/v.drape for clone removal
00009  ***************************************************
00010  */
00011 
00012 #include <string.h>
00013 #include <unistd.h>
00014 #include <math.h>
00015 #include <grass/gis.h>
00016 #include <grass/glocale.h>
00017 
00018 
00019 /* prototypes */
00020 static double raster_sample_nearest (int fd,
00021                   struct Cell_head *window, 
00022                   struct Categories *cats,
00023                   double north, double east, int usedesc);
00024 static double raster_sample_bilinear (int fd,
00025                   struct Cell_head *window, 
00026                   struct Categories *cats,
00027                   double north, double east, int usedesc);
00028 static double raster_sample_cubic (int fd,
00029                   struct Cell_head *window, 
00030                   struct Categories *cats,
00031                   double north, double east, int usedesc);
00032 static double scancatlabel (char *);
00033 static void raster_row_error(struct Cell_head *window, double north, double east);
00034 
00035 
00054 double G_get_raster_sample (int fd,
00055                 struct Cell_head *window,
00056                 struct Categories *cats,
00057                 double north, double east,
00058                 int usedesc, INTERP_TYPE itype)
00059 {
00060     double retval;
00061 
00062     switch (itype)
00063     {
00064     case NEAREST:
00065         retval = raster_sample_nearest(fd, window, cats, north, east, usedesc);
00066     break;
00067     case BILINEAR:
00068         retval = raster_sample_bilinear(fd, window, cats, north, east, usedesc);
00069     break;
00070     case CUBIC:
00071         retval = raster_sample_cubic(fd, window, cats, north, east, usedesc);
00072     break;
00073     default:
00074         G_fatal_error(_("unknown interpolation type"));
00075     }
00076 
00077     return retval;
00078 }
00079 
00080 
00081 static double raster_sample_nearest (int fd,
00082                   struct Cell_head *window,
00083                   struct Categories *cats,
00084                   double north, double east, int usedesc)
00085 {
00086     int row, col;
00087     double predicted;
00088     DCELL *maprow = G_allocate_d_raster_buf();
00089 
00090     /* convert northing and easting to row and col, resp */
00091     row = (int) G_northing_to_row (north, window);
00092     col = (int) G_easting_to_col (east, window);
00093 
00094     if (G_get_d_raster_row (fd, maprow, row) < 0)
00095         raster_row_error(window, north, east);
00096 
00097     if (G_is_d_null_value(&(maprow[col])))
00098     {
00099         predicted = 0.0;
00100     }
00101     else if (usedesc)
00102     {
00103         char *buf = G_get_cat(maprow[col], cats);
00104 
00105         G_squeeze(buf);
00106         predicted = scancatlabel(buf);
00107     }
00108     else
00109     {
00110         predicted = maprow[col];
00111     }
00112 
00113     G_free(maprow);
00114 
00115     return predicted;
00116 }
00117 
00118 
00119 static double raster_sample_bilinear (int fd,
00120                   struct Cell_head *window,
00121                   struct Categories *cats,
00122                   double north, double east, int usedesc)
00123 {
00124     int i, row, col;
00125     double grid[2][2], tmp1, tmp2;
00126     DCELL *arow = G_allocate_d_raster_buf ();
00127     DCELL *brow = G_allocate_d_raster_buf ();
00128 
00129     /* convert northing and easting to row and col, resp */
00130     row = (int) G_northing_to_row (north, window);
00131     col = (int) G_easting_to_col (east, window);
00132 
00133     if (G_get_d_raster_row (fd, arow, row) < 0)
00134         raster_row_error(window, north, east);
00135 
00136     /*-
00137      * we need 2x2 pixels to do the interpolation. First we decide if we
00138      * need the previous or next map row
00139      */
00140     if (row == 0)
00141     {
00142         /* arow is at top, must get row below */
00143         if (G_get_d_raster_row (fd, brow, row + 1) < 0)
00144             raster_row_error(window, north, east);
00145     }
00146     else if (row+1 == G_window_rows ())
00147     {
00148         /* amaprow is at bottom, get the one above it */
00149         /* brow = arow; */
00150         for(i = 0; i < G_window_cols(); ++i)
00151             brow[i] = arow[i];
00152 
00153         row--;
00154 
00155         if (G_get_d_raster_row (fd, arow, row) < 0)
00156             raster_row_error(window, north, east);
00157     }
00158     else if (north - G_row_to_northing ((double) row + 0.5, window) > 0)
00159     {
00160         /* north is above a horizontal centerline going through arow */
00161         /* brow = arow; */
00162         for(i = 0; i < G_window_cols(); ++i)
00163             brow[i] = arow[i];
00164 
00165         row--;
00166 
00167         if (G_get_d_raster_row (fd, arow, row) < 0)
00168             raster_row_error(window, north, east);
00169     }
00170     else
00171     {
00172         /* north is below a horizontal centerline going through arow */
00173         if (G_get_d_raster_row (fd, brow, row+1) < 0)
00174             raster_row_error(window, north, east);
00175     }
00176 
00177     /*-
00178      * next, we decide if we need the column to the right or left of the
00179      * current column using a procedure similar to above
00180      */
00181     if (col+1 == G_window_cols())
00182         col--;
00183     else if (east - G_col_to_easting((double) col + 0.5, window) < 0)
00184         col--;
00185 
00186     /*-
00187      * now were ready to do bilinear interpolation over
00188      * arow[col], arow[col+1],
00189      * brow[col], brow[col+1]
00190      */
00191  
00192     if (usedesc)
00193     {
00194         char *buf;
00195 
00196         G_squeeze(buf = G_get_cat ((int)arow[col], cats));
00197         grid[0][0] = scancatlabel (buf);
00198         G_squeeze(buf = G_get_cat ((int)arow[col+1], cats));
00199         grid[0][1] = scancatlabel (buf);
00200         G_squeeze(buf = G_get_cat ((int)brow[col], cats));
00201         grid[1][0] = scancatlabel (buf);
00202         G_squeeze(buf = G_get_cat ((int)brow[col+1], cats));
00203         grid[1][1] = scancatlabel (buf);
00204     }
00205     else
00206     {
00207         grid[0][0] = (double) arow[col];
00208         grid[0][1] = (double) arow[col + 1];
00209         grid[1][0] = (double) brow[col];
00210         grid[1][1] = (double) brow[col + 1];
00211     }
00212 
00213     /* Treat NULL's as zero */
00214     if (G_is_d_null_value(&(arow[col])))
00215         grid[0][0] = 0.0;
00216     if (G_is_d_null_value(&(arow[col+1])))
00217         grid[0][1] = 0.0;
00218     if (G_is_d_null_value(&(brow[col])))
00219         grid[1][0] = 0.0;
00220     if (G_is_d_null_value(&(brow[col+1])))
00221         grid[1][1] = 0.0;
00222 
00223     east = fabs(G_col_to_easting((double)col, window) - east);
00224     while (east > window->ew_res)
00225         east -= window->ew_res;
00226 
00227     north = fabs(G_row_to_northing((double)row, window) - north);
00228     while (north > window->ns_res)
00229         north -= window->ns_res;
00230 
00231     /* we do two linear interpolations along the rows */
00232     tmp1 = east * grid[0][1] + (window->ew_res - east) * grid[0][0];
00233     tmp1 /= window->ew_res;
00234     tmp2 = east * grid[1][1] + (window->ew_res - east) * grid[1][0];
00235     tmp2 /= window->ew_res;
00236     G_debug(3, "DIAG: r=%d c=%d t1=%g t2=%g\te=%g n=%g",
00237             row, col, tmp1, tmp2, east, north);
00238     G_debug(3, "DIAG: %g %g\n      %g %g",
00239             grid[0][0], grid[0][1], grid[1][0], grid[1][1]);
00240 
00241     /*-
00242      * Now we interpolate along a line parallel to columns
00243      * and passing through easting
00244      */
00245     G_free(arow);
00246     G_free(brow);
00247 
00248     return (north * tmp2 + (window->ns_res - north) * tmp1) / window->ns_res;
00249 }
00250 
00251 
00252 static double raster_sample_cubic (int fd,
00253                   struct Cell_head *window, 
00254                   struct Categories *cats,
00255                   double north, double east, int usedesc)
00256 {
00257     int i, row, col;
00258     double grid[4][4], tmp[4];
00259     DCELL *arow = G_allocate_d_raster_buf();
00260     DCELL *brow = G_allocate_d_raster_buf();
00261     DCELL *crow = G_allocate_d_raster_buf();
00262     DCELL *drow = G_allocate_d_raster_buf();
00263 
00264     /* convert northing and easting to row and col, resp */
00265     row = G_northing_to_row(north, window);
00266     col = G_easting_to_col(east, window);
00267 
00268     if (G_get_d_raster_row(fd, arow, row) < 0)
00269         raster_row_error(window, north, east);
00270 
00271     /* we need 4x4 pixels to do the interpolation. */
00272     if (row == 0)
00273     {
00274         /* row containing sample is at top, must get three rows below */
00275         if (G_get_d_raster_row(fd, brow, row + 1) < 0)
00276             raster_row_error(window, north, east);
00277         if (G_get_d_raster_row(fd, crow, row + 2) < 0)
00278             raster_row_error(window, north, east);
00279         if (G_get_d_raster_row(fd, drow, row + 3) < 0)
00280             raster_row_error(window, north, east);
00281     }
00282     else if (row == 1)
00283     {
00284         /* must get row above and tow rows below */
00285         for (i = 0; i < G_window_cols(); ++i)
00286             brow[i] = arow[i];
00287 
00288         if (G_get_d_raster_row(fd, arow, row - 1) < 0)
00289             raster_row_error(window, north, east);
00290         if (G_get_d_raster_row(fd, crow, row + 1) < 0)
00291             raster_row_error(window, north, east);
00292         if (G_get_d_raster_row(fd, drow, row + 2) < 0)
00293             raster_row_error(window, north, east);
00294 
00295         row--;
00296     }
00297     else if (row + 1 == G_window_rows())
00298     {
00299         /* arow is at bottom, get the three above it */
00300         for (i = 0; i < G_window_cols(); ++i)
00301             drow[i] = arow[i];
00302 
00303         if (G_get_d_raster_row(fd, arow, row - 3) < 0)
00304             raster_row_error(window, north, east);
00305         if (G_get_d_raster_row(fd, brow, row - 2) < 0)
00306             raster_row_error(window, north, east);
00307         if (G_get_d_raster_row(fd, crow, row - 1) < 0)
00308             raster_row_error(window, north, east);
00309 
00310         row -= 3;
00311     }
00312     else if (row + 2 == G_window_rows())
00313     {
00314         /* arow is next to bottom, get the one below and two above it */
00315         for (i = 0; i < G_window_cols(); ++i)
00316             crow[i] = arow[i];
00317 
00318         if (G_get_d_raster_row(fd, arow, row - 2) < 0)
00319             raster_row_error(window, north, east);
00320         if (G_get_d_raster_row(fd, brow, row - 1) < 0)
00321             raster_row_error(window, north, east);
00322         if (G_get_d_raster_row(fd, drow, row + 1) < 0)
00323             raster_row_error(window, north, east);
00324 
00325         row -= 2;
00326     }
00327     else if (north - G_row_to_northing((double)row + 0.5, window) > 0)
00328     {
00329         /*
00330          * north is above a horizontal centerline going through arow. need two
00331          * above and one below
00332          */
00333         for (i = 0; i < G_window_cols(); ++i)
00334             crow[i] = arow[i];
00335 
00336         if (G_get_d_raster_row(fd, arow, row - 2) < 0)
00337             raster_row_error(window, north, east);
00338         if (G_get_d_raster_row(fd, brow, row - 1) < 0)
00339             raster_row_error(window, north, east);
00340         if (G_get_d_raster_row(fd, drow, row + 1) < 0)
00341             raster_row_error(window, north, east);
00342 
00343         row -= 2;
00344     }
00345     else
00346     {
00347         /*
00348          * north is below a horizontal centerline going through arow need one
00349          * above and two below
00350          */
00351         for (i = 0; i < G_window_cols(); ++i)
00352             brow[i] = arow[i];
00353 
00354         if (G_get_d_raster_row(fd, arow, row - 1) < 0)
00355             raster_row_error(window, north, east);
00356         if (G_get_d_raster_row(fd, crow, row + 1) < 0)
00357             raster_row_error(window, north, east);
00358         if (G_get_d_raster_row(fd, drow, row + 2) < 0)
00359             raster_row_error(window, north, east);
00360 
00361         row--;
00362     }
00363 
00364     /*
00365      * Next, we decide if we need columns to the right and/or left of the
00366      * current column using a procedure similar to above
00367      */
00368     if (col == 0 || col == 1)
00369         col = 0;
00370     else if (col + 1 == G_window_cols())
00371         col -= 3;
00372     else if (col + 2 == G_window_cols())
00373         col -= 2;
00374     else if (east - G_col_to_easting((double)col + 0.5, window) < 0)
00375         /* east is left of center */
00376         col -= 2;
00377     else
00378         col--;
00379 
00380     /*
00381      * now were ready to do cubic interpolation over
00382      * arow[col], arow[col+1], arow[col+2], arow[col+3],
00383      * brow[col], brow[col+1], brow[col+2], brow[col+3],
00384      * crow[col], crow[col+1], crow[col+2], crow[col+3],
00385      * drow[col], drow[col+1], drow[col+2], drow[col+3],
00386      */
00387 
00388     if (usedesc)
00389     {
00390         char *buf;
00391 
00392         for (i = 0; i < 4; ++i) {
00393             G_squeeze(buf = G_get_cat(arow[col + i], cats));
00394             grid[0][i] = scancatlabel(buf);
00395             G_squeeze(buf = G_get_cat(brow[col + i], cats));
00396             grid[1][i] = scancatlabel(buf);
00397             G_squeeze(buf = G_get_cat(crow[col + i], cats));
00398             grid[2][i] = scancatlabel(buf);
00399             G_squeeze(buf = G_get_cat(drow[col + i], cats));
00400             grid[3][i] = scancatlabel(buf);
00401         }
00402     }
00403     else
00404     {
00405         for (i = 0; i < 4; ++i) {
00406             grid[0][i] = (double)arow[col + i];
00407             grid[1][i] = (double)brow[col + i];
00408             grid[2][i] = (double)crow[col + i];
00409             grid[3][i] = (double)drow[col + i];
00410         }
00411     }
00412 
00413     /* Treat NULL cells as 0.0 */
00414     for (i = 0; i < 4; i++) {
00415         if (G_is_d_null_value(&(arow[col + i])))
00416             grid[0][i] = 0.0;
00417         if (G_is_d_null_value(&(brow[col + i])))
00418             grid[1][i] = 0.0;
00419         if (G_is_d_null_value(&(crow[col + i])))
00420             grid[2][i] = 0.0;
00421         if (G_is_d_null_value(&(drow[col + i])))
00422             grid[3][i] = 0.0;
00423     }
00424 
00425     /* this needs work here */
00426     east = fabs(G_col_to_easting((double)col + 1, window) - east);
00427     while (east > window->ew_res)
00428         east -= window->ew_res;
00429     east /= window->ew_res;
00430 
00431     north = fabs(G_row_to_northing((double)row + 1, window) - north);
00432     while (north > window->ns_res)
00433         north -= window->ns_res;
00434     north /= window->ns_res;
00435 
00436     /* we do four cubic convolutions along the rows */
00437     for (i = 0; i < 4; ++i)
00438         tmp[i] = east * (east * (east * (grid[i][3] - grid[i][2] + grid[i][1] - grid[i][0])
00439                     + (grid[i][2] - grid[i][3] - 2 * grid[i][1] + 2 * grid[i][0]))
00440                     + (grid[i][2] - grid[i][0])) + grid[i][1];
00441 
00442 #ifdef DEBUG
00443     for (i = 0; i < 4; ++i) {
00444         for (j = 0; j < 4; ++j)
00445             fprintf(stderr, "%g ", grid[i][j]);
00446         fprintf(stderr, "\n");
00447     }
00448 
00449     G_message("DIAG: (%d,%d) 1=%3.2g 2=%3.2g 3=%3.2g 4=%3.2g\te=%g n=%g",
00450               row, col, tmp[0], tmp[1], tmp[2], tmp[3], east, north);
00451 #endif
00452 
00453     G_free(arow);
00454     G_free(brow);
00455     G_free(crow);
00456     G_free(drow);
00457 
00458     /* user horner's method again for the final interpolation */
00459     return (north * (north * (north * (tmp[3] - tmp[2] + tmp[1] - tmp[0])
00460                            + (tmp[2] - tmp[3] - 2 * tmp[1] + 2 * tmp[0]))
00461                            + (tmp[2] - tmp[0])) + tmp[1]);
00462 }
00463 
00464 
00465 static double scancatlabel (char *str)
00466 {
00467     double val;
00468 
00469     if (strcmp(str,"no data") != 0)
00470         sscanf(str, "%lf", &val);
00471     else
00472     {
00473         G_warning(_("\"no data\" label found; setting to zero"));
00474         val = 0.0;
00475     }
00476 
00477     return val;
00478 }
00479 
00480 
00481 static void raster_row_error(struct Cell_head *window, double north, double east)
00482 {
00483     G_debug(3, "DIAG: \tRegion is: n=%g s=%g e=%g w=%g",
00484          window->north, window->south, window->east, window->west);
00485     G_debug(3, "      \tData point is north=%g east=%g", north, east);
00486 
00487     G_fatal_error(_("problem reading raster cell file"));
00488 }

Generated on Wed Dec 19 14:59:06 2007 for GRASS by  doxygen 1.5.4