1 /*!
2 \file lib/ogsf/gsd_views.c
3
4 \brief OGSF library - manipulating views (lower level functions)
5
6 GRASS OpenGL gsurf OGSF Library
7
8 (C) 1999-2008 by the GRASS Development Team
9
10 This program is free software under the
11 GNU General Public License (>=v2).
12 Read the file COPYING that comes with GRASS
13 for details.
14
15 \author Bill Brown USACERL (January 1993)
16 \author Doxygenized by Martin Landa <landa.martin gmail.com> (May 2008)
17 */
18
19 #include <grass/config.h>
20
21 #if defined(OPENGL_X11) || defined(OPENGL_WINDOWS)
22 #include <GL/gl.h>
23 #include <GL/glu.h>
24 #elif defined(OPENGL_AQUA)
25 #include <OpenGL/gl.h>
26 #include <OpenGL/glu.h>
27 #endif
28
29 #include <grass/ogsf.h>
30 #include "math.h"
31
32 /*!
33 \brief ADD
34
35 \param vect
36 \param sx, sy screen coordinates
37
38 \return 1
39 */
gsd_get_los(float (* vect)[3],short sx,short sy)40 int gsd_get_los(float (*vect)[3], short sx, short sy)
41 {
42 double fx, fy, fz, tx, ty, tz;
43 GLdouble modelMatrix[16], projMatrix[16];
44 GLint viewport[4];
45
46 GS_ready_draw();
47 glPushMatrix();
48
49 gsd_do_scale(1);
50 glGetDoublev(GL_MODELVIEW_MATRIX, modelMatrix);
51 glGetDoublev(GL_PROJECTION_MATRIX, projMatrix);
52 glGetIntegerv(GL_VIEWPORT, viewport);
53 glPopMatrix();
54
55 /* OGLXXX XXX I think this is backwards gluProject(XXX); */
56 /* WAS: mapw(Vobj, sx, sy, &fx, &fy, &fz, &tx, &ty, &tz); */
57 gluUnProject((GLdouble) sx, (GLdouble) sy, 0.0, modelMatrix,
58 projMatrix, viewport, &fx, &fy, &fz);
59 gluUnProject((GLdouble) sx, (GLdouble) sy, 1.0, modelMatrix,
60 projMatrix, viewport, &tx, &ty, &tz);
61 vect[FROM][X] = fx;
62 vect[FROM][Y] = fy;
63 vect[FROM][Z] = fz;
64 vect[TO][X] = tx;
65 vect[TO][Y] = ty;
66 vect[TO][Z] = tz;
67
68 /* DEBUG - should just be a dot */
69 /* OGLXXX frontbuffer: other possibilities include GSD_BOTH */
70 GS_set_draw((1) ? GSD_FRONT : GSD_BACK);
71 glPushMatrix();
72 gsd_do_scale(1);
73 gsd_linewidth(3);
74 gsd_color_func(0x8888FF);
75
76 /* OGLXXX for multiple, independent line segments: use GL_LINES */
77 glBegin(GL_LINE_STRIP);
78 glVertex3fv(vect[FROM]);
79 glVertex3fv(vect[TO]);
80 glEnd();
81
82 gsd_linewidth(1);
83 glPopMatrix();
84
85 /* OGLXXX frontbuffer: other possibilities include GSD_BOTH */
86 GS_set_draw((0) ? GSD_FRONT : GSD_BACK);
87
88 return (1);
89 }
90
91 #if 0
92 /*!
93 \brief Set view
94
95 Establishes viewing & projection matrices
96
97 \param gv view (geoview)
98 \param dp display (geodisplay)
99 */
100 void gsd_set_view(geoview * gv, geodisplay * gd)
101 {
102 double up[3];
103 GLint mm;
104
105 /* will expand when need to check for in focus, ortho, etc. */
106
107 gsd_check_focus(gv);
108 gsd_get_zup(gv, up);
109
110 gd->aspect = GS_get_aspect();
111
112 glGetIntegerv(GL_MATRIX_MODE, &mm);
113 glMatrixMode(GL_PROJECTION);
114 glLoadIdentity();
115 gluPerspective((double).1 * (gv->fov), (double)gd->aspect,
116 (double)gd->nearclip, (double)gd->farclip);
117
118 glMatrixMode(mm);
119
120 glLoadIdentity();
121
122 /* update twist parm */
123 glRotatef((float)(gv->twist / 10.), 0.0, 0.0, 1.0);
124
125 /* OGLXXX lookat: replace UPx with vector */
126 gluLookAt((double)gv->from_to[FROM][X], (double)gv->from_to[FROM][Y],
127 (double)gv->from_to[FROM][Z], (double)gv->from_to[TO][X],
128 (double)gv->from_to[TO][Y], (double)gv->from_to[TO][Z],
129 (double)up[X], (double)up[Y], (double)up[Z]);
130
131 /* have to redefine clipping planes when view changes */
132
133 gsd_update_cplanes();
134
135 return;
136 }
137 #endif
138 /*!
139 \brief Set view
140
141 Establishes viewing & projection matrices
142
143 \param gv view (geoview)
144 \param dp display (geodisplay)
145 */
gsd_set_view(geoview * gv,geodisplay * gd)146 void gsd_set_view(geoview * gv, geodisplay * gd)
147 {
148 double up[3];
149 float pos[3];
150 int i;
151 GLdouble modelMatrix[16];
152 GLint mm;
153
154 /* will expand when need to check for in focus, ortho, etc. */
155
156 gsd_check_focus(gv);
157 gsd_get_zup(gv, up);
158
159 gd->aspect = GS_get_aspect();
160
161 glGetIntegerv(GL_MATRIX_MODE, &mm);
162 glMatrixMode(GL_PROJECTION);
163 glLoadIdentity();
164 gluPerspective((double).1 * (gv->fov), (double)gd->aspect,
165 (double)gd->nearclip, (double)gd->farclip);
166
167 glMatrixMode(mm);
168
169 glLoadIdentity();
170
171 /* update twist parm */
172 glRotatef((float)(gv->twist / 10.), 0.0, 0.0, 1.0);
173
174 /* OGLXXX lookat: replace UPx with vector */
175 gluLookAt((double)gv->from_to[FROM][X], (double)gv->from_to[FROM][Y],
176 (double)gv->from_to[FROM][Z], (double)gv->from_to[TO][X],
177 (double)gv->from_to[TO][Y], (double)gv->from_to[TO][Z],
178 (double)up[X], (double)up[Y], (double)up[Z]);
179
180 /* rotate to get rotation matrix and then save it*/
181 if (gv->rotate.do_rot) {
182
183 glPushMatrix();
184 glLoadMatrixd(gv->rotate.rotMatrix);
185
186 glRotated(gv->rotate.rot_angle, gv->rotate.rot_axes[0],
187 gv->rotate.rot_axes[1], gv->rotate.rot_axes[2]);
188 glGetDoublev(GL_MODELVIEW_MATRIX, modelMatrix);
189
190 for (i = 0; i < 16; i++) {
191 gv->rotate.rotMatrix[i] = modelMatrix[i];
192 }
193
194 glPopMatrix();
195 }
196
197 gs_get_datacenter(pos);
198 gsd_surf2model(pos);
199 /* translate rotation center to view center, rotate and translate back */
200 glTranslatef(pos[0], pos[1], pos[2]);
201 glMultMatrixd(gv->rotate.rotMatrix);
202 glTranslatef(-pos[0], -pos[1], -pos[2]);
203
204 /* have to redefine clipping planes when view changes */
205
206 gsd_update_cplanes();
207
208 return;
209 }
210 /*!
211 \brief Check focus
212
213 \param gv view (geoview)
214 */
gsd_check_focus(geoview * gv)215 void gsd_check_focus(geoview * gv)
216 {
217 float zmax, zmin;
218
219 GS_get_zrange(&zmin, &zmax, 0);
220
221 if (gv->infocus) {
222 GS_v3eq(gv->from_to[TO], gv->real_to);
223 gv->from_to[TO][Z] -= zmin;
224 GS_v3mult(gv->from_to[TO], gv->scale);
225 gv->from_to[TO][Z] *= gv->vert_exag;
226
227 GS_v3normalize(gv->from_to[FROM], gv->from_to[TO]);
228 }
229
230 return;
231 }
232
233 /*!
234 \brief Get z-up vector (z-direction)
235
236 \param gv view (geoview)
237 \param up up vector
238 */
gsd_get_zup(geoview * gv,double * up)239 void gsd_get_zup(geoview * gv, double *up)
240 {
241 float alpha;
242 float zup[3], fup[3];
243
244 /* neg alpha OK since sin(-x) = -sin(x) */
245 alpha =
246 (2.0 * atan(1.0)) - acos(gv->from_to[FROM][Z] - gv->from_to[TO][Z]);
247
248 zup[X] = gv->from_to[TO][X];
249 zup[Y] = gv->from_to[TO][Y];
250
251 if (sin(alpha)) {
252 zup[Z] = gv->from_to[TO][Z] + 1 / sin(alpha);
253 }
254 else {
255 zup[Z] = gv->from_to[FROM][Z] + 1.0;
256 }
257
258 GS_v3dir(gv->from_to[FROM], zup, fup);
259
260 up[X] = fup[X];
261 up[Y] = fup[Y];
262 up[Z] = fup[Z];
263
264 return;
265 }
266
267 /*!
268 \brief ADD
269
270 \param gv view (geoview)
271
272 \return ?
273 */
gsd_zup_twist(geoview * gv)274 int gsd_zup_twist(geoview * gv)
275 {
276 float fr_to[2][4];
277 float look_theta, pi;
278 float alpha, beta;
279 float zup[3], yup[3], zupmag, yupmag;
280
281 pi = 4.0 * atan(1.0);
282
283 /* *************************************************************** */
284 /* This block of code is used to keep pos z in the up direction,
285 * correcting for SGI system default which is pos y in the up
286 * direction. Involves finding up vectors for both y up and
287 * z up, then determining angle between them. LatLon mode uses y as
288 * up direction instead of z, so no correction necessary. Next rewrite,
289 * we should use y as up for all drawing.
290 */
291 GS_v3eq(fr_to[FROM], gv->from_to[FROM]);
292 GS_v3eq(fr_to[TO], gv->from_to[TO]);
293
294 /* neg alpha OK since sin(-x) = -sin(x) */
295 alpha = pi / 2.0 - acos(fr_to[FROM][Z] - fr_to[TO][Z]);
296
297 zup[X] = fr_to[TO][X];
298 zup[Y] = fr_to[TO][Y];
299
300 if (sin(alpha)) {
301 zup[Z] = fr_to[TO][Z] + 1 / sin(alpha);
302 }
303 else {
304 zup[Z] = fr_to[FROM][Z] + 1.0;
305 }
306
307 zupmag = GS_distance(fr_to[FROM], zup);
308
309 yup[X] = fr_to[TO][X];
310 yup[Z] = fr_to[TO][Z];
311
312 /* neg beta OK since sin(-x) = -sin(x) */
313 beta = pi / 2.0 - acos(fr_to[TO][Y] - fr_to[FROM][Y]);
314
315 if (sin(beta)) {
316 yup[Y] = fr_to[TO][Y] - 1 / sin(beta);
317 }
318 else {
319 yup[Y] = fr_to[FROM][Y] + 1.0;
320 }
321
322 yupmag = GS_distance(fr_to[FROM], yup);
323
324 look_theta = (1800.0 / pi) *
325 acos(((zup[X] - fr_to[FROM][X]) * (yup[X] - fr_to[FROM][X])
326 + (zup[Y] - fr_to[FROM][Y]) * (yup[Y] - fr_to[FROM][Y])
327 + (zup[Z] - fr_to[FROM][Z]) * (yup[Z] - fr_to[FROM][Z])) /
328 (zupmag * yupmag));
329
330 if (fr_to[TO][X] - fr_to[FROM][X] < 0.0) {
331 look_theta = -look_theta;
332 }
333
334 if (fr_to[TO][Z] - fr_to[FROM][Z] < 0.0) {
335 /* looking down */
336 if (fr_to[TO][Y] - fr_to[FROM][Y] < 0.0) {
337 look_theta = 1800 - look_theta;
338 }
339 }
340 else {
341 /* looking up */
342 if (fr_to[TO][Y] - fr_to[FROM][Y] > 0.0) {
343 look_theta = 1800 - look_theta;
344 }
345 }
346
347 return ((int)(gv->twist + 1800 + look_theta));
348 }
349
350 /*!
351 \brief Set current scale
352
353 \param doexag use z-exaggeration
354 */
gsd_do_scale(int doexag)355 void gsd_do_scale(int doexag)
356 {
357 float sx, sy, sz;
358 float min, max;
359
360 GS_get_scale(&sx, &sy, &sz, doexag);
361 gsd_scale(sx, sy, sz);
362 GS_get_zrange(&min, &max, 0);
363 gsd_translate(0.0, 0.0, -min);
364
365 return;
366 }
367
368 /*!
369 \brief Convert real to model coordinates
370
371 \param point[in,out] 3d point (Point3)
372 */
gsd_real2model(Point3 point)373 void gsd_real2model(Point3 point)
374 {
375 float sx, sy, sz;
376 float min, max, n, s, w, e;
377
378 GS_get_region(&n, &s, &w, &e);
379 GS_get_scale(&sx, &sy, &sz, 1);
380 GS_get_zrange(&min, &max, 0);
381 point[X] = (point[X] - w) * sx;
382 point[Y] = (point[Y] - s) * sy;
383 point[Z] = (point[Z] - min) * sz;
384
385 return;
386 }
387
388 /*!
389 \brief Convert model to real coordinates
390
391 \param point[in,out] 3d point (x,y,z)
392 */
gsd_model2real(Point3 point)393 void gsd_model2real(Point3 point)
394 {
395 float sx, sy, sz;
396 float min, max, n, s, w, e;
397
398 GS_get_region(&n, &s, &w, &e);
399 GS_get_scale(&sx, &sy, &sz, 1);
400 GS_get_zrange(&min, &max, 0);
401 point[X] = (sx ? point[X] / sx : 0.0) + w;
402 point[Y] = (sy ? point[Y] / sy : 0.0) + s;
403 point[Z] = (sz ? point[Z] / sz : 0.0) + min;
404
405 return;
406 }
407
408 /*!
409 \brief Convert model to surface coordinates
410
411 \param gs surface (geosurf)
412 \param point 3d point (Point3)
413 */
gsd_model2surf(geosurf * gs,Point3 point)414 void gsd_model2surf(geosurf * gs, Point3 point)
415 {
416 float min, max, sx, sy, sz;
417
418 /* so far, only one geographic "region" allowed, so origin of
419 surface is same as origin of model space, but will need to provide
420 translations here to make up the difference, so not using gs yet */
421
422 if (gs) {
423 /* need to undo z scaling & translate */
424 GS_get_scale(&sx, &sy, &sz, 1);
425 GS_get_zrange(&min, &max, 0);
426
427 point[Z] = (sz ? point[Z] / sz : 0.0) + min;
428
429 /* need to unscale x & y */
430 point[X] = (sx ? point[X] / sx : 0.0);
431 point[Y] = (sy ? point[Y] / sy : 0.0);
432 }
433
434 return;
435 }
436 /*!
437 \brief Convert surface to model coordinates
438
439 \param point 3d point (Point3)
440 */
gsd_surf2model(Point3 point)441 void gsd_surf2model(Point3 point)
442 {
443 float min, max, sx, sy, sz;
444
445 /* need to undo z scaling & translate */
446 GS_get_scale(&sx, &sy, &sz, 1);
447 GS_get_zrange(&min, &max, 0);
448
449 point[Z] = (sz ? (point[Z] - min) * sz : 0.0);
450
451 /* need to unscale x & y */
452 point[X] = (sx ? point[X] * sx : 0.0);
453 point[Y] = (sy ? point[Y] * sy : 0.0);
454
455
456 return;
457 }
458 /*!
459 \brief Convert surface to real coordinates
460
461 \param gs surface (geosurf)
462 \param[in,out] point 3d point (Point3)
463 */
gsd_surf2real(geosurf * gs,Point3 point)464 void gsd_surf2real(geosurf * gs, Point3 point)
465 {
466 if (gs) {
467 point[X] += gs->ox;
468 point[Y] += gs->oy;
469 }
470
471 return;
472 }
473
474 /*!
475 \brief Convert real to surface coordinates
476
477 \param gs surface (geosurf)
478 \param[in,out] point 3d point (Point3)
479 */
gsd_real2surf(geosurf * gs,Point3 point)480 void gsd_real2surf(geosurf * gs, Point3 point)
481 {
482 if (gs) {
483 point[X] -= gs->ox;
484 point[Y] -= gs->oy;
485 }
486
487 return;
488 }
489