GRASS Programmer's Manual  6.4.2(2012)
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
gsd_views.c
Go to the documentation of this file.
1 
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/gstypes.h>
30 #include "math.h"
31 
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 GL_FRONT_AND_BACK */
70  glDrawBuffer((1) ? GL_FRONT : GL_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 GL_FRONT_AND_BACK */
86  glDrawBuffer((0) ? GL_FRONT : GL_BACK);
87 
88  return (1);
89 }
90 
91 
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 
134 
135  return;
136 }
137 
143 void gsd_check_focus(geoview * gv)
144 {
145  float zmax, zmin;
146 
147  GS_get_zrange(&zmin, &zmax, 0);
148 
149  if (gv->infocus) {
150  GS_v3eq(gv->from_to[TO], gv->real_to);
151  gv->from_to[TO][Z] -= zmin;
152  GS_v3mult(gv->from_to[TO], gv->scale);
153  gv->from_to[TO][Z] *= gv->vert_exag;
154 
155  GS_v3normalize(gv->from_to[FROM], gv->from_to[TO]);
156  }
157 
158  return;
159 }
160 
167 void gsd_get_zup(geoview * gv, double *up)
168 {
169  float alpha;
170  float zup[3], fup[3];
171 
172  /* neg alpha OK since sin(-x) = -sin(x) */
173  alpha =
174  (2.0 * atan(1.0)) - acos(gv->from_to[FROM][Z] - gv->from_to[TO][Z]);
175 
176  zup[X] = gv->from_to[TO][X];
177  zup[Y] = gv->from_to[TO][Y];
178 
179  if (sin(alpha)) {
180  zup[Z] = gv->from_to[TO][Z] + 1 / sin(alpha);
181  }
182  else {
183  zup[Z] = gv->from_to[FROM][Z] + 1.0;
184  }
185 
186  GS_v3dir(gv->from_to[FROM], zup, fup);
187 
188  up[X] = fup[X];
189  up[Y] = fup[Y];
190  up[Z] = fup[Z];
191 
192  return;
193 }
194 
202 int gsd_zup_twist(geoview * gv)
203 {
204  float fr_to[2][4];
205  float look_theta, pi;
206  float alpha, beta;
207  float zup[3], yup[3], zupmag, yupmag;
208 
209  pi = 4.0 * atan(1.0);
210 
211  /* *************************************************************** */
212  /* This block of code is used to keep pos z in the up direction,
213  * correcting for SGI system default which is pos y in the up
214  * direction. Involves finding up vectors for both y up and
215  * z up, then determining angle between them. LatLon mode uses y as
216  * up direction instead of z, so no correction necessary. Next rewrite,
217  * we should use y as up for all drawing.
218  */
219  GS_v3eq(fr_to[FROM], gv->from_to[FROM]);
220  GS_v3eq(fr_to[TO], gv->from_to[TO]);
221 
222  /* neg alpha OK since sin(-x) = -sin(x) */
223  alpha = pi / 2.0 - acos(fr_to[FROM][Z] - fr_to[TO][Z]);
224 
225  zup[X] = fr_to[TO][X];
226  zup[Y] = fr_to[TO][Y];
227 
228  if (sin(alpha)) {
229  zup[Z] = fr_to[TO][Z] + 1 / sin(alpha);
230  }
231  else {
232  zup[Z] = fr_to[FROM][Z] + 1.0;
233  }
234 
235  zupmag = GS_distance(fr_to[FROM], zup);
236 
237  yup[X] = fr_to[TO][X];
238  yup[Z] = fr_to[TO][Z];
239 
240  /* neg beta OK since sin(-x) = -sin(x) */
241  beta = pi / 2.0 - acos(fr_to[TO][Y] - fr_to[FROM][Y]);
242 
243  if (sin(beta)) {
244  yup[Y] = fr_to[TO][Y] - 1 / sin(beta);
245  }
246  else {
247  yup[Y] = fr_to[FROM][Y] + 1.0;
248  }
249 
250  yupmag = GS_distance(fr_to[FROM], yup);
251 
252  look_theta = (1800.0 / pi) *
253  acos(((zup[X] - fr_to[FROM][X]) * (yup[X] - fr_to[FROM][X])
254  + (zup[Y] - fr_to[FROM][Y]) * (yup[Y] - fr_to[FROM][Y])
255  + (zup[Z] - fr_to[FROM][Z]) * (yup[Z] - fr_to[FROM][Z])) /
256  (zupmag * yupmag));
257 
258  if (fr_to[TO][X] - fr_to[FROM][X] < 0.0) {
259  look_theta = -look_theta;
260  }
261 
262  if (fr_to[TO][Z] - fr_to[FROM][Z] < 0.0) {
263  /* looking down */
264  if (fr_to[TO][Y] - fr_to[FROM][Y] < 0.0) {
265  look_theta = 1800 - look_theta;
266  }
267  }
268  else {
269  /* looking up */
270  if (fr_to[TO][Y] - fr_to[FROM][Y] > 0.0) {
271  look_theta = 1800 - look_theta;
272  }
273  }
274 
275  return ((int)(gv->twist + 1800 + look_theta));
276 }
277 
283 void gsd_do_scale(int doexag)
284 {
285  float sx, sy, sz;
286  float min, max;
287 
288  GS_get_scale(&sx, &sy, &sz, doexag);
289  gsd_scale(sx, sy, sz);
290  GS_get_zrange(&min, &max, 0);
291  gsd_translate(0.0, 0.0, -min);
292 
293  return;
294 }
295 
301 void gsd_real2model(Point3 point)
302 {
303  float sx, sy, sz;
304  float min, max, n, s, w, e;
305 
306  GS_get_region(&n, &s, &w, &e);
307  GS_get_scale(&sx, &sy, &sz, 1);
308  GS_get_zrange(&min, &max, 0);
309  point[X] = (point[X] - w) * sx;
310  point[Y] = (point[Y] - s) * sy;
311  point[Z] = (point[Z] - min) * sz;
312 
313  return;
314 }
315 
321 void gsd_model2real(Point3 point)
322 {
323  float sx, sy, sz;
324  float min, max, n, s, w, e;
325 
326  GS_get_region(&n, &s, &w, &e);
327  GS_get_scale(&sx, &sy, &sz, 1);
328  GS_get_zrange(&min, &max, 0);
329  point[X] = (sx ? point[X] / sx : 0.0) + w;
330  point[Y] = (sy ? point[Y] / sy : 0.0) + s;
331  point[Z] = (sz ? point[Z] / sz : 0.0) + min;
332 
333  return;
334 }
335 
342 void gsd_model2surf(geosurf * gs, Point3 point)
343 {
344  float min, max, sx, sy, sz;
345 
346  /* so far, only one geographic "region" allowed, so origin of
347  surface is same as origin of model space, but will need to provide
348  translations here to make up the difference, so not using gs yet */
349 
350  if (gs) {
351  /* need to undo z scaling & translate */
352  GS_get_scale(&sx, &sy, &sz, 1);
353  GS_get_zrange(&min, &max, 0);
354 
355  point[Z] = (sz ? point[Z] / sz : 0.0) + min;
356 
357  /* need to unscale x & y */
358  point[X] = (sx ? point[X] / sx : 0.0);
359  point[Y] = (sy ? point[Y] / sy : 0.0);
360  }
361 
362  return;
363 }
364 
371 void gsd_surf2real(geosurf * gs, Point3 point)
372 {
373  if (gs) {
374  point[X] += gs->ox;
375  point[Y] += gs->oy;
376  }
377 
378  return;
379 }
380 
387 void gsd_real2surf(geosurf * gs, Point3 point)
388 {
389  if (gs) {
390  point[X] -= gs->ox;
391  point[Y] -= gs->oy;
392  }
393 
394  return;
395 }