GRASS Programmer's Manual  6.4.2(2012)
gs_norms.c
Go to the documentation of this file.
00001 
00019 #include <math.h>
00020 
00021 #include <grass/gis.h>
00022 #include <grass/gstypes.h>
00023 
00024 #include "gsget.h"
00025 #include "rowcol.h"
00026 
00027 #define NTOP 0x00001000
00028 #define NBOT 0x00000100
00029 #define NLFT 0x00000010
00030 #define NRGT 0x00000001
00031 
00032 #define NALL 0x00001111
00033 
00034 #define NTL  0x00001010
00035 #define NTR  0x00001001
00036 #define NBL  0x00000110
00037 #define NBR  0x00000101
00038 
00042 #define SET_NORM(i) \
00043        dz1 = z1 - z2; \
00044        dz2 = z3 - z4; \
00045        temp[0] = (float) -dz1 * y_res_z2; \
00046        temp[1] = (float) dz2 * x_res_z2; \
00047        temp[2] = c_z2; \
00048        normalizer = sqrt(temp[0] * temp[0] + temp[1] * temp[1] + c_z2_sq); \
00049        if (!normalizer) normalizer= 1.0; \
00050        temp[X] /= normalizer; \
00051        temp[Y] /= normalizer; \
00052        temp[Z] /= normalizer; \
00053        PNORM(i,temp);
00054 
00055 static long slice;
00056 static float x_res_z2, y_res_z2;
00057 static float c_z2, c_z2_sq;
00058 static typbuff *elbuf;
00059 static unsigned long *norm;
00060 
00061 /*
00062    #define USE_GL_NORMALIZE
00063  */
00064 
00072 void init_vars(geosurf * gs)
00073 {
00074     /* optimized - these are static - global to this file */
00075     norm = gs->norms;
00076     elbuf = gs_get_att_typbuff(gs, ATT_TOPO, 0);
00077 
00078 #ifdef USE_GL_NORMALIZE
00079     c_z2 =
00080         2.0 * gs->xres * gs->yres * gs->x_mod * gs->y_mod / GS_global_exag();
00081     c_z2_sq = c_z2 * c_z2;
00082     x_res_z2 = 2.0 * gs->xres * gs->z_exag * gs->x_mod;
00083     y_res_z2 = 2.0 * gs->yres * gs->z_exag * gs->y_mod;
00084 #else
00085 
00086     {
00087         float sx, sy, sz;
00088 
00089         GS_get_scale(&sx, &sy, &sz, 1);
00090 
00091         c_z2 = 2.0 * gs->xres * gs->yres * gs->x_mod * gs->y_mod;
00092         c_z2_sq = c_z2 * c_z2;
00093         x_res_z2 = 2.0 * gs->xres * gs->z_exag * gs->x_mod;
00094         y_res_z2 = 2.0 * gs->yres * gs->z_exag * gs->y_mod;
00095     }
00096 #endif
00097 
00098     slice = gs->y_mod * gs->cols;
00099 
00100     return;
00101 }
00102 
00124 int gs_calc_normals(geosurf * gs)
00125 {
00126     int row, col;
00127     int xcnt, ycnt;
00128     int xmod, ymod;
00129 
00130     if (!gs->norm_needupdate || !gs->norms) {
00131         return (0);
00132     }
00133 
00134     gs->norm_needupdate = 0;
00135     gs_update_curmask(gs);
00136 
00137     xmod = gs->x_mod;
00138     ymod = gs->y_mod;
00139 
00140     xcnt = VCOLS(gs);
00141     ycnt = VROWS(gs);
00142 
00143     init_vars(gs);
00144 
00145     G_debug(5, "gs_calc_normals(): id=%d", gs->gsurf_id);
00146 
00147     /* first row - just use single cell */
00148     /* first col - use bottom & right neighbors */
00149     calc_norm(gs, 0, 0, NBR);
00150 
00151     for (col = 1; col < xcnt; col++) {
00152         /* turn off top neighbor for first row */
00153         calc_norm(gs, 0, col * xmod, ~NTOP);
00154     }
00155 
00156     /* use bottom & left neighbors for last col */
00157     calc_norm(gs, 0, col * xmod, NBL);
00158 
00159     /* now use four neighboring points for rows 1 - (n-1) */
00160     for (row = 1; row < ycnt; row++) {
00161         if (!(row % 100))
00162             G_debug(5, "gs_calc_normals(): row=%d", row);
00163 
00164         /* turn off left neighbor for first col */
00165         calc_norm(gs, row * ymod, 0, ~NLFT);
00166 
00167         /* use all 4 neighbors until last col */
00168         for (col = 1; col < xcnt; col++) {
00169             calc_norm(gs, row * ymod, col * xmod, NALL);
00170         }
00171 
00172         /* turn off right neighbor for last col */
00173         calc_norm(gs, row * ymod, col * xmod, ~NRGT);
00174     }
00175 
00176     /* last row */
00177     /* use top & right neighbors for first col */
00178     calc_norm(gs, row * ymod, 0, NTR);
00179 
00180     for (col = 1; col < xcnt; col++) {
00181         /* turn off bottom neighbor for last row */
00182         calc_norm(gs, row * ymod, col * xmod, ~NBOT);
00183     }
00184 
00185     /* use top & left neighbors for last column */
00186     calc_norm(gs, row * ymod, col * xmod, NTL);
00187 
00188     return (1);
00189 }
00190 
00206 int calc_norm(geosurf * gs, int drow, int dcol, unsigned int neighbors)
00207 {
00208     long noffset;
00209     float temp[3], normalizer, dz1, dz2, z0, z1, z2, z3, z4;
00210 
00211     if (gs->curmask) {
00212         /* need to check masked neighbors */
00213         /* NOTE: this should automatically eliminate nullvals */
00214         if (neighbors & NTOP) {
00215             if (BM_get(gs->curmask, dcol, drow - gs->y_mod)) {
00216                 /* masked */
00217                 neighbors &= ~NTOP;
00218             }
00219         }
00220 
00221         if (neighbors & NBOT) {
00222             if (BM_get(gs->curmask, dcol, drow + gs->y_mod)) {
00223                 /* masked */
00224                 neighbors &= ~NBOT;
00225             }
00226         }
00227 
00228         if (neighbors & NLFT) {
00229             if (BM_get(gs->curmask, dcol - gs->x_mod, drow)) {
00230                 /* masked */
00231                 neighbors &= ~NLFT;
00232             }
00233         }
00234 
00235         if (neighbors & NRGT) {
00236             if (BM_get(gs->curmask, dcol + gs->x_mod, drow)) {
00237                 /* masked */
00238                 neighbors &= ~NRGT;
00239             }
00240         }
00241     }
00242 
00243     if (!neighbors) {
00244         /* none */
00245         return (0);
00246     }
00247 
00248     noffset = DRC2OFF(gs, drow, dcol);
00249 
00250     if (!GET_MAPATT(elbuf, noffset, z0)) {
00251         return (0);
00252     }
00253 
00254     z1 = z2 = z3 = z4 = z0;
00255 
00256     /* we know these aren't null now, maybe use faster GET_MAPATT? */
00257     if (neighbors & NRGT) {
00258         GET_MAPATT(elbuf, noffset + gs->x_mod, z1);
00259         if (!(neighbors & NLFT)) {
00260             z2 = z0 + (z0 - z1);
00261         }
00262     }
00263 
00264     if (neighbors & NLFT) {
00265         GET_MAPATT(elbuf, noffset - gs->x_mod, z2);
00266 
00267         if (!(neighbors & NRGT)) {
00268             z1 = z0 + (z0 - z2);
00269         }
00270     }
00271 
00272     if (neighbors & NTOP) {
00273         GET_MAPATT(elbuf, noffset - slice, z4);
00274 
00275         if (!(neighbors & NBOT)) {
00276             z3 = z0 + (z0 - z4);
00277         }
00278     }
00279 
00280     if (neighbors & NBOT) {
00281         GET_MAPATT(elbuf, noffset + slice, z3);
00282 
00283         if (!(neighbors & NTOP)) {
00284             z4 = z0 + (z0 - z3);
00285         }
00286     }
00287 
00288     SET_NORM(norm[noffset]);
00289 
00290     return (1);
00291 }
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines