GRASS Programmer's Manual
6.4.2(2012)
|
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 }