GRASS Programmer's Manual  6.5.svn(2014)-r66266
 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 #if 0
92 
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 #endif
138 
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 
207 
208  return;
209 }
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 
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 
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 
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 
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 
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 
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 }
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 }
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 
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 }
int GS_v3normalize(float *v1, float *v2)
Change v2 so that v1v2 is a unit vector.
Definition: GS_util.c:322
int gs_get_datacenter(float *cen)
Get data center point.
Definition: gs.c:1232
void gsd_get_zup(geoview *gv, double *up)
Get z-up vector (z-direction)
Definition: gsd_views.c:239
#define min(x, y)
Definition: draw2.c:68
void gsd_do_scale(int doexag)
Set current scale.
Definition: gsd_views.c:355
#define Y(x)
Definition: display/draw.c:246
tuple pos
Definition: tools.py:1367
int GS_v3dir(float *v1, float *v2, float *v3)
Get a normalized direction from v1 to v2, store in v3.
Definition: GS_util.c:353
void gsd_surf2real(geosurf *gs, Point3 point)
Convert surface to real coordinates.
Definition: gsd_views.c:464
#define X(y)
Definition: display/draw.c:248
void gsd_color_func(unsigned int col)
Set current color.
Definition: gsd_prim.c:689
#define max(x, y)
Definition: draw2.c:69
int GS_get_region(float *n, float *s, float *w, float *e)
Get 2D region extent.
Definition: GS2.c:156
void gsd_translate(float dx, float dy, float dz)
Multiply the current matrix by a translation matrix.
Definition: gsd_prim.c:526
void gsd_check_focus(geoview *gv)
Check focus.
Definition: gsd_views.c:215
void gsd_update_cplanes(void)
Update cplaces.
Definition: gsd_cplane.c:86
void gsd_real2surf(geosurf *gs, Point3 point)
Convert real to surface coordinates.
Definition: gsd_views.c:480
void gsd_real2model(Point3 point)
Convert real to model coordinates.
Definition: gsd_views.c:373
int GS_get_zrange(float *min, float *max, int doexag)
Get z-extent for all loaded surfaces.
Definition: GS2.c:2687
double GS_get_aspect(void)
Get aspect value.
Definition: GS2.c:3450
void GS_v3mult(float *v1, float k)
Multiple vectors.
Definition: GS_util.c:229
void gsd_set_view(geoview *gv, geodisplay *gd)
Set view.
Definition: gsd_views.c:146
void gsd_model2surf(geosurf *gs, Point3 point)
Convert model to surface coordinates.
Definition: gsd_views.c:414
void gsd_model2real(Point3 point)
Convert model to real coordinates.
Definition: gsd_views.c:393
void GS_v3eq(float *v1, float *v2)
Copy vector values.
Definition: GS_util.c:178
void gsd_linewidth(short n)
Set width of rasterized lines.
Definition: gsd_prim.c:257
float GS_distance(float *from, float *to)
Calculate distance.
Definition: GS_util.c:141
#define FROM
Definition: y.tab.c:161
void gsd_surf2model(Point3 point)
Convert surface to model coordinates.
Definition: gsd_views.c:441
void GS_ready_draw(void)
Definition: GS2.c:2486
void gsd_scale(float xs, float ys, float zs)
Multiply the current matrix by a general scaling matrix.
Definition: gsd_prim.c:512
int n
Definition: dataquad.c:291
int gsd_get_los(float(*vect)[3], short sx, short sy)
ADD.
Definition: gsd_views.c:40
int gsd_zup_twist(geoview *gv)
ADD.
Definition: gsd_views.c:274
void GS_get_scale(float *sx, float *sy, float *sz, int doexag)
Get axis scale.
Definition: GS2.c:3238