GRASS Programmer's Manual
6.4.2(2012)
|
00001 00019 #include <grass/config.h> 00020 00021 #if defined(OPENGL_X11) || defined(OPENGL_WINDOWS) 00022 #include <GL/gl.h> 00023 #include <GL/glu.h> 00024 #elif defined(OPENGL_AQUA) 00025 #include <OpenGL/gl.h> 00026 #include <OpenGL/glu.h> 00027 #endif 00028 00029 #include <grass/gstypes.h> 00030 #include "math.h" 00031 00040 int gsd_get_los(float (*vect)[3], short sx, short sy) 00041 { 00042 double fx, fy, fz, tx, ty, tz; 00043 GLdouble modelMatrix[16], projMatrix[16]; 00044 GLint viewport[4]; 00045 00046 GS_ready_draw(); 00047 glPushMatrix(); 00048 00049 gsd_do_scale(1); 00050 glGetDoublev(GL_MODELVIEW_MATRIX, modelMatrix); 00051 glGetDoublev(GL_PROJECTION_MATRIX, projMatrix); 00052 glGetIntegerv(GL_VIEWPORT, viewport); 00053 glPopMatrix(); 00054 00055 /* OGLXXX XXX I think this is backwards gluProject(XXX); */ 00056 /* WAS: mapw(Vobj, sx, sy, &fx, &fy, &fz, &tx, &ty, &tz); */ 00057 gluUnProject((GLdouble) sx, (GLdouble) sy, 0.0, modelMatrix, 00058 projMatrix, viewport, &fx, &fy, &fz); 00059 gluUnProject((GLdouble) sx, (GLdouble) sy, 1.0, modelMatrix, 00060 projMatrix, viewport, &tx, &ty, &tz); 00061 vect[FROM][X] = fx; 00062 vect[FROM][Y] = fy; 00063 vect[FROM][Z] = fz; 00064 vect[TO][X] = tx; 00065 vect[TO][Y] = ty; 00066 vect[TO][Z] = tz; 00067 00068 /* DEBUG - should just be a dot */ 00069 /* OGLXXX frontbuffer: other possibilities include GL_FRONT_AND_BACK */ 00070 glDrawBuffer((1) ? GL_FRONT : GL_BACK); 00071 glPushMatrix(); 00072 gsd_do_scale(1); 00073 gsd_linewidth(3); 00074 gsd_color_func(0x8888FF); 00075 00076 /* OGLXXX for multiple, independent line segments: use GL_LINES */ 00077 glBegin(GL_LINE_STRIP); 00078 glVertex3fv(vect[FROM]); 00079 glVertex3fv(vect[TO]); 00080 glEnd(); 00081 00082 gsd_linewidth(1); 00083 glPopMatrix(); 00084 00085 /* OGLXXX frontbuffer: other possibilities include GL_FRONT_AND_BACK */ 00086 glDrawBuffer((0) ? GL_FRONT : GL_BACK); 00087 00088 return (1); 00089 } 00090 00091 00100 void gsd_set_view(geoview * gv, geodisplay * gd) 00101 { 00102 double up[3]; 00103 GLint mm; 00104 00105 /* will expand when need to check for in focus, ortho, etc. */ 00106 00107 gsd_check_focus(gv); 00108 gsd_get_zup(gv, up); 00109 00110 gd->aspect = GS_get_aspect(); 00111 00112 glGetIntegerv(GL_MATRIX_MODE, &mm); 00113 glMatrixMode(GL_PROJECTION); 00114 glLoadIdentity(); 00115 gluPerspective((double).1 * (gv->fov), (double)gd->aspect, 00116 (double)gd->nearclip, (double)gd->farclip); 00117 00118 glMatrixMode(mm); 00119 00120 glLoadIdentity(); 00121 00122 /* update twist parm */ 00123 glRotatef((float)(gv->twist / 10.), 0.0, 0.0, 1.0); 00124 00125 /* OGLXXX lookat: replace UPx with vector */ 00126 gluLookAt((double)gv->from_to[FROM][X], (double)gv->from_to[FROM][Y], 00127 (double)gv->from_to[FROM][Z], (double)gv->from_to[TO][X], 00128 (double)gv->from_to[TO][Y], (double)gv->from_to[TO][Z], 00129 (double)up[X], (double)up[Y], (double)up[Z]); 00130 00131 /* have to redefine clipping planes when view changes */ 00132 00133 gsd_update_cplanes(); 00134 00135 return; 00136 } 00137 00143 void gsd_check_focus(geoview * gv) 00144 { 00145 float zmax, zmin; 00146 00147 GS_get_zrange(&zmin, &zmax, 0); 00148 00149 if (gv->infocus) { 00150 GS_v3eq(gv->from_to[TO], gv->real_to); 00151 gv->from_to[TO][Z] -= zmin; 00152 GS_v3mult(gv->from_to[TO], gv->scale); 00153 gv->from_to[TO][Z] *= gv->vert_exag; 00154 00155 GS_v3normalize(gv->from_to[FROM], gv->from_to[TO]); 00156 } 00157 00158 return; 00159 } 00160 00167 void gsd_get_zup(geoview * gv, double *up) 00168 { 00169 float alpha; 00170 float zup[3], fup[3]; 00171 00172 /* neg alpha OK since sin(-x) = -sin(x) */ 00173 alpha = 00174 (2.0 * atan(1.0)) - acos(gv->from_to[FROM][Z] - gv->from_to[TO][Z]); 00175 00176 zup[X] = gv->from_to[TO][X]; 00177 zup[Y] = gv->from_to[TO][Y]; 00178 00179 if (sin(alpha)) { 00180 zup[Z] = gv->from_to[TO][Z] + 1 / sin(alpha); 00181 } 00182 else { 00183 zup[Z] = gv->from_to[FROM][Z] + 1.0; 00184 } 00185 00186 GS_v3dir(gv->from_to[FROM], zup, fup); 00187 00188 up[X] = fup[X]; 00189 up[Y] = fup[Y]; 00190 up[Z] = fup[Z]; 00191 00192 return; 00193 } 00194 00202 int gsd_zup_twist(geoview * gv) 00203 { 00204 float fr_to[2][4]; 00205 float look_theta, pi; 00206 float alpha, beta; 00207 float zup[3], yup[3], zupmag, yupmag; 00208 00209 pi = 4.0 * atan(1.0); 00210 00211 /* *************************************************************** */ 00212 /* This block of code is used to keep pos z in the up direction, 00213 * correcting for SGI system default which is pos y in the up 00214 * direction. Involves finding up vectors for both y up and 00215 * z up, then determining angle between them. LatLon mode uses y as 00216 * up direction instead of z, so no correction necessary. Next rewrite, 00217 * we should use y as up for all drawing. 00218 */ 00219 GS_v3eq(fr_to[FROM], gv->from_to[FROM]); 00220 GS_v3eq(fr_to[TO], gv->from_to[TO]); 00221 00222 /* neg alpha OK since sin(-x) = -sin(x) */ 00223 alpha = pi / 2.0 - acos(fr_to[FROM][Z] - fr_to[TO][Z]); 00224 00225 zup[X] = fr_to[TO][X]; 00226 zup[Y] = fr_to[TO][Y]; 00227 00228 if (sin(alpha)) { 00229 zup[Z] = fr_to[TO][Z] + 1 / sin(alpha); 00230 } 00231 else { 00232 zup[Z] = fr_to[FROM][Z] + 1.0; 00233 } 00234 00235 zupmag = GS_distance(fr_to[FROM], zup); 00236 00237 yup[X] = fr_to[TO][X]; 00238 yup[Z] = fr_to[TO][Z]; 00239 00240 /* neg beta OK since sin(-x) = -sin(x) */ 00241 beta = pi / 2.0 - acos(fr_to[TO][Y] - fr_to[FROM][Y]); 00242 00243 if (sin(beta)) { 00244 yup[Y] = fr_to[TO][Y] - 1 / sin(beta); 00245 } 00246 else { 00247 yup[Y] = fr_to[FROM][Y] + 1.0; 00248 } 00249 00250 yupmag = GS_distance(fr_to[FROM], yup); 00251 00252 look_theta = (1800.0 / pi) * 00253 acos(((zup[X] - fr_to[FROM][X]) * (yup[X] - fr_to[FROM][X]) 00254 + (zup[Y] - fr_to[FROM][Y]) * (yup[Y] - fr_to[FROM][Y]) 00255 + (zup[Z] - fr_to[FROM][Z]) * (yup[Z] - fr_to[FROM][Z])) / 00256 (zupmag * yupmag)); 00257 00258 if (fr_to[TO][X] - fr_to[FROM][X] < 0.0) { 00259 look_theta = -look_theta; 00260 } 00261 00262 if (fr_to[TO][Z] - fr_to[FROM][Z] < 0.0) { 00263 /* looking down */ 00264 if (fr_to[TO][Y] - fr_to[FROM][Y] < 0.0) { 00265 look_theta = 1800 - look_theta; 00266 } 00267 } 00268 else { 00269 /* looking up */ 00270 if (fr_to[TO][Y] - fr_to[FROM][Y] > 0.0) { 00271 look_theta = 1800 - look_theta; 00272 } 00273 } 00274 00275 return ((int)(gv->twist + 1800 + look_theta)); 00276 } 00277 00283 void gsd_do_scale(int doexag) 00284 { 00285 float sx, sy, sz; 00286 float min, max; 00287 00288 GS_get_scale(&sx, &sy, &sz, doexag); 00289 gsd_scale(sx, sy, sz); 00290 GS_get_zrange(&min, &max, 0); 00291 gsd_translate(0.0, 0.0, -min); 00292 00293 return; 00294 } 00295 00301 void gsd_real2model(Point3 point) 00302 { 00303 float sx, sy, sz; 00304 float min, max, n, s, w, e; 00305 00306 GS_get_region(&n, &s, &w, &e); 00307 GS_get_scale(&sx, &sy, &sz, 1); 00308 GS_get_zrange(&min, &max, 0); 00309 point[X] = (point[X] - w) * sx; 00310 point[Y] = (point[Y] - s) * sy; 00311 point[Z] = (point[Z] - min) * sz; 00312 00313 return; 00314 } 00315 00321 void gsd_model2real(Point3 point) 00322 { 00323 float sx, sy, sz; 00324 float min, max, n, s, w, e; 00325 00326 GS_get_region(&n, &s, &w, &e); 00327 GS_get_scale(&sx, &sy, &sz, 1); 00328 GS_get_zrange(&min, &max, 0); 00329 point[X] = (sx ? point[X] / sx : 0.0) + w; 00330 point[Y] = (sy ? point[Y] / sy : 0.0) + s; 00331 point[Z] = (sz ? point[Z] / sz : 0.0) + min; 00332 00333 return; 00334 } 00335 00342 void gsd_model2surf(geosurf * gs, Point3 point) 00343 { 00344 float min, max, sx, sy, sz; 00345 00346 /* so far, only one geographic "region" allowed, so origin of 00347 surface is same as origin of model space, but will need to provide 00348 translations here to make up the difference, so not using gs yet */ 00349 00350 if (gs) { 00351 /* need to undo z scaling & translate */ 00352 GS_get_scale(&sx, &sy, &sz, 1); 00353 GS_get_zrange(&min, &max, 0); 00354 00355 point[Z] = (sz ? point[Z] / sz : 0.0) + min; 00356 00357 /* need to unscale x & y */ 00358 point[X] = (sx ? point[X] / sx : 0.0); 00359 point[Y] = (sy ? point[Y] / sy : 0.0); 00360 } 00361 00362 return; 00363 } 00364 00371 void gsd_surf2real(geosurf * gs, Point3 point) 00372 { 00373 if (gs) { 00374 point[X] += gs->ox; 00375 point[Y] += gs->oy; 00376 } 00377 00378 return; 00379 } 00380 00387 void gsd_real2surf(geosurf * gs, Point3 point) 00388 { 00389 if (gs) { 00390 point[X] -= gs->ox; 00391 point[Y] -= gs->oy; 00392 } 00393 00394 return; 00395 }