src/Scenery/tilemgr.cxx

   1 // tilemgr.cxx -- routines to handle dynamic management of scenery tiles
   2 //
   3 // Written by Curtis Olson, started January 1998.
   4 //
   5 // Copyright (C) 1997  Curtis L. Olson  - curt@infoplane.com
   6 //
   7 // This program is free software; you can redistribute it and/or
   8 // modify it under the terms of the GNU General Public License as
   9 // published by the Free Software Foundation; either version 2 of the
  10 // License, or (at your option) any later version.
  11 //
  12 // This program is distributed in the hope that it will be useful, but
  13 // WITHOUT ANY WARRANTY; without even the implied warranty of
  14 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  15 // General Public License for more details.
  16 //
  17 // You should have received a copy of the GNU General Public License
  18 // along with this program; if not, write to the Free Software
  19 // Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  20 //
  21 // $Id$
  22
  23
  24 #ifdef HAVE_CONFIG_H
  25 #  include <config.h>
  26 #endif
  27
  28 #ifdef HAVE_WINDOWS_H
  29 #  include <windows.h>
  30 #endif
  31
  32 #include <GL/glut.h>
  33 #include <XGL/xgl.h>
  34
  35 #include <Aircraft/aircraft.hxx>
  36
  37 #include <Debug/logstream.hxx>
  38 // #include <Bucket/bucketutils.hxx>
  39 #include <Include/fg_constants.h>
  40 #include <Main/options.hxx>
  41 #include <Main/views.hxx>
  42 #include <Math/fg_geodesy.hxx>
  43 #include <Math/mat3.h>
  44 #include <Math/point3d.hxx>
  45 #include <Math/polar3d.hxx>
  46 #include <Math/vector.hxx>
  47 #include <Objects/materialmgr.hxx>
  48 #include <Objects/obj.hxx>
  49
  50 #ifndef FG_OLD_WEATHER
  51 #  include <WeatherCM/FGLocalWeatherDatabase.h>
  52 #else
  53 #  include <Weather/weather.hxx>
  54 #endif
  55
  56 #include "scenery.hxx"
  57 #include "tilecache.hxx"
  58 #include "tileentry.hxx"
  59 #include "tilemgr.hxx"
  60
  61
  62 // to test clipping speedup in fgTileMgrRender()
  63 #if defined ( USE_FAST_FOV_CLIP )
  64 // #define TEST_FOV_CLIP
  65 // #define TEST_ELEV
  66 #endif
  67
  68
  69 extern ssgRoot *scene;
  70
  71
  72 // the tile manager
  73 FGTileMgr global_tile_mgr;
  74
  75
  76 // Constructor
  77 FGTileMgr::FGTileMgr ( void ):
  78     state( Start )
  79 {
  80 }
  81
  82
  83 // Destructor
  84 FGTileMgr::~FGTileMgr ( void ) {
  85 }
  86
  87
  88 // Initialize the Tile Manager subsystem
  89 int FGTileMgr::init( void ) {
  90     FG_LOG( FG_TERRAIN, FG_INFO, "Initializing Tile Manager subsystem." );
  91
  92     // load default material library
  93     if ( ! material_mgr.loaded() ) {
  94         material_mgr.load_lib();
  95     }
  96
  97     global_tile_cache.init();
  98
  99     state = Inited;
 100
 101     return 1;
 102 }
 103
 104
 105 // schedule a tile for loading
 106 static void disable_tile( int cache_index ) {
 107     // see if tile already exists in the cache
 108     // cout << "DISABLING CACHE ENTRY = " << cache_index << endl;
 109     FGTileEntry *t = global_tile_cache.get_tile( cache_index );
 110     t->ssg_disable();
 111 }
 112
 113
 114 // schedule a tile for loading
 115 int FGTileMgr::sched_tile( const FGBucket& b ) {
 116     // see if tile already exists in the cache
 117     int cache_index = global_tile_cache.exists( b );
 118
 119     if ( cache_index >= 0 ) {
 120         // tile exists in cache, reenable it.
 121         // cout << "REENABLING DISABLED TILE" << endl;
 122         FGTileEntry *t = global_tile_cache.get_tile( cache_index );
 123         t->select_ptr->select( 1 );
 124         t->mark_loaded();
 125     } else {
 126         // find the next available cache entry and mark it as
 127         // scheduled
 128         cache_index = global_tile_cache.next_avail();
 129         FGTileEntry *t = global_tile_cache.get_tile( cache_index );
 130         t->mark_scheduled_for_use();
 131
 132         // register a load request
 133         FGLoadRec request;
 134         request.b = b;
 135         request.cache_index = cache_index;
 136         load_queue.push_back( request );
 137     }
 138
 139     return cache_index;
 140 }
 141
 142
 143 // load a tile
 144 void FGTileMgr::load_tile( const FGBucket& b, int cache_index) {
 145
 146     FG_LOG( FG_TERRAIN, FG_DEBUG, "Loading tile " << b );
 147
 148     global_tile_cache.fill_in(cache_index, b);
 149
 150     FG_LOG( FG_TERRAIN, FG_DEBUG, "Loaded for cache index: " << cache_index );
 151 }
 152
 153
 154 // Calculate shortest distance from point to line
 155 static double point_line_dist_squared( const Point3D& tc, const Point3D& vp,
 156                                        MAT3vec d )
 157 {
 158     MAT3vec p, p0;
 159
 160     p[0] = tc.x(); p[1] = tc.y(); p[2] = tc.z();
 161     p0[0] = vp.x(); p0[1] = vp.y(); p0[2] = vp.z();
 162
 163     return fgPointLineSquared(p, p0, d);
 164 }
 165
 166
 167 // Determine scenery altitude.  Normally this just happens when we
 168 // render the scene, but we'd also like to be able to do this
 169 // explicitely.  lat & lon are in radians.  abs_view_pos in meters.
 170 // Returns result in meters.
 171 double
 172 FGTileMgr::current_elev_new( const FGBucket& p ) {
 173     FGTileEntry *t;
 174     fgFRAGMENT *frag_ptr;
 175     Point3D abs_view_pos = current_view.get_abs_view_pos();
 176     Point3D earth_center(0.0);
 177     Point3D result;
 178     MAT3vec local_up;
 179     double dist, lat_geod, alt, sea_level_r;
 180     int index;
 181
 182     local_up[0] = abs_view_pos.x();
 183     local_up[1] = abs_view_pos.y();
 184     local_up[2] = abs_view_pos.z();
 185
 186     // Find current translation offset
 187     // fgBucketFind(lon * RAD_TO_DEG, lat * RAD_TO_DEG, &p);
 188     index = global_tile_cache.exists(p);
 189     if ( index < 0 ) {
 190         FG_LOG( FG_TERRAIN, FG_WARN, "Tile not found" );
 191         return 0.0;
 192     }
 193
 194     t = global_tile_cache.get_tile(index);
 195
 196     scenery.next_center = t->center;
 197
 198     FG_LOG( FG_TERRAIN, FG_DEBUG,
 199             "Current bucket = " << p << "  Index = " << p.gen_index_str() );
 200     FG_LOG( FG_TERRAIN, FG_DEBUG,
 201             "abs_view_pos = " << abs_view_pos );
 202
 203     // calculate tile offset
 204     // x = (t->offset.x = t->center.x - scenery.center.x);
 205     // y = (t->offset.y = t->center.y - scenery.center.y);
 206     // z = (t->offset.z = t->center.z - scenery.center.z);
 207
 208     // calc current terrain elevation calculate distance from
 209     // vertical tangent line at current position to center of
 210     // tile.
 211
 212     /* printf("distance squared = %.2f, bounding radius = %.2f\n",
 213        point_line_dist_squared(&(t->offset), &(v->view_pos),
 214        v->local_up), t->bounding_radius); */
 215
 216     dist = point_line_dist_squared( t->center, abs_view_pos, local_up );
 217     if ( dist < FG_SQUARE(t->bounding_radius) ) {
 218
 219         // traverse fragment list for tile
 220         FGTileEntry::FragmentIterator current = t->begin();
 221         FGTileEntry::FragmentIterator last = t->end();
 222
 223         for ( ; current != last; ++current ) {
 224             frag_ptr = &(*current);
 225             /* printf("distance squared = %.2f, bounding radius = %.2f\n",
 226                point_line_dist_squared( &(frag_ptr->center),
 227                &abs_view_pos), local_up),
 228                frag_ptr->bounding_radius); */
 229
 230             dist = point_line_dist_squared( frag_ptr->center,
 231                                             abs_view_pos,
 232                                             local_up);
 233             if ( dist <= FG_SQUARE(frag_ptr->bounding_radius) ) {
 234                 if ( frag_ptr->intersect( abs_view_pos,
 235                                           earth_center, 0, result ) ) {
 236                     FG_LOG( FG_TERRAIN, FG_DEBUG, "intersection point " <<
 237                             result );
 238                     // compute geocentric coordinates of tile center
 239                     Point3D pp = fgCartToPolar3d(result);
 240                     FG_LOG( FG_TERRAIN, FG_DEBUG, "  polar form = " << pp );
 241                     // convert to geodetic coordinates
 242                     fgGeocToGeod(pp.lat(), pp.radius(), &lat_geod,
 243                                  &alt, &sea_level_r);
 244
 245                     // printf("alt = %.2f\n", alt);
 246                     // exit since we found an intersection
 247                     if ( alt > -9999.0 ) {
 248                         // printf("returning alt\n");
 249                         return alt;
 250                     } else {
 251                         // printf("returning 0\n");
 252                         return 0.0;
 253                     }
 254                 }
 255             }
 256         }
 257     }
 258
 259     FG_LOG( FG_TERRAIN, FG_INFO, "(new) no terrain intersection found" );
 260
 261     return 0.0;
 262 }
 263
 264
 265 // Determine scenery altitude.  Normally this just happens when we
 266 // render the scene, but we'd also like to be able to do this
 267 // explicitely.  lat & lon are in radians.  abs_view_pos in meters.
 268 // Returns result in meters.
 269 double
 270 FGTileMgr::current_elev( double lon, double lat, const Point3D& abs_view_pos ) {
 271     FGTileCache *c;
 272     FGTileEntry *t;
 273     fgFRAGMENT *frag_ptr;
 274     Point3D earth_center(0.0);
 275     Point3D result;
 276     MAT3vec local_up;
 277     double dist, lat_geod, alt, sea_level_r;
 278     int index;
 279
 280     c = &global_tile_cache;
 281
 282     local_up[0] = abs_view_pos.x();
 283     local_up[1] = abs_view_pos.y();
 284     local_up[2] = abs_view_pos.z();
 285
 286     FG_LOG( FG_TERRAIN, FG_DEBUG, "Absolute view pos = " << abs_view_pos );
 287
 288     // Find current translation offset
 289     FGBucket p( lon * RAD_TO_DEG, lat * RAD_TO_DEG );
 290     index = c->exists(p);
 291     if ( index < 0 ) {
 292         FG_LOG( FG_TERRAIN, FG_WARN, "Tile not found" );
 293         return 0.0;
 294     }
 295
 296     t = c->get_tile(index);
 297
 298     scenery.next_center = t->center;
 299
 300     FG_LOG( FG_TERRAIN, FG_DEBUG,
 301             "Pos = (" << lon * RAD_TO_DEG << ", " << lat * RAD_TO_DEG
 302             << ")  Current bucket = " << p
 303             << "  Index = " << p.gen_index_str() );
 304
 305     FG_LOG( FG_TERRAIN, FG_DEBUG, "Tile center " << t->center
 306             << "  bounding radius = " << t->bounding_radius );
 307
 308     // calculate tile offset
 309     // x = (t->offset.x = t->center.x - scenery.center.x);
 310     // y = (t->offset.y = t->center.y - scenery.center.y);
 311     // z = (t->offset.z = t->center.z - scenery.center.z);
 312
 313     // calc current terrain elevation calculate distance from
 314     // vertical tangent line at current position to center of
 315     // tile.
 316
 317     /* printf("distance squared = %.2f, bounding radius = %.2f\n",
 318        point_line_dist_squared(&(t->offset), &(v->view_pos),
 319        v->local_up), t->bounding_radius); */
 320
 321     dist = point_line_dist_squared( t->center, abs_view_pos, local_up );
 322     FG_LOG( FG_TERRAIN, FG_DEBUG, "(gross check) dist squared = " << dist );
 323
 324     if ( dist < FG_SQUARE(t->bounding_radius) ) {
 325
 326         // traverse fragment list for tile
 327         FGTileEntry::FragmentIterator current = t->begin();
 328         FGTileEntry::FragmentIterator last = t->end();
 329
 330         for ( ; current != last; ++current ) {
 331             frag_ptr = &(*current);
 332             /* printf("distance squared = %.2f, bounding radius = %.2f\n",
 333                point_line_dist_squared( &(frag_ptr->center),
 334                &abs_view_pos), local_up),
 335                frag_ptr->bounding_radius); */
 336
 337             dist = point_line_dist_squared( frag_ptr->center,
 338                                             abs_view_pos,
 339                                             local_up);
 340             if ( dist <= FG_SQUARE(frag_ptr->bounding_radius) ) {
 341                 if ( frag_ptr->intersect( abs_view_pos,
 342                                           earth_center, 0, result ) ) {
 343                     FG_LOG( FG_TERRAIN, FG_DEBUG, "intersection point " <<
 344                             result );
 345                     // compute geocentric coordinates of tile center
 346                     Point3D pp = fgCartToPolar3d(result);
 347                     FG_LOG( FG_TERRAIN, FG_DEBUG, "  polar form = " << pp );
 348                     // convert to geodetic coordinates
 349                     fgGeocToGeod(pp.lat(), pp.radius(), &lat_geod,
 350                                  &alt, &sea_level_r);
 351
 352                     // printf("alt = %.2f\n", alt);
 353                     // exit since we found an intersection
 354                     if ( alt > -9999.0 ) {
 355                         // printf("returning alt\n");
 356                         return alt;
 357                     } else {
 358                         // printf("returning 0\n");
 359                         return 0.0;
 360                     }
 361                 }
 362             }
 363         }
 364     }
 365
 366     FG_LOG( FG_TERRAIN, FG_INFO, "(old) no terrain intersection found" );
 367
 368     return 0.0;
 369 }
 370
 371
 372 inline int fg_sign( const double x ) {
 373     return x < 0 ? -1 : 1;
 374 }
 375
 376 inline double fg_min( const double a, const double b ) {
 377     return b < a ? b : a;
 378 }
 379
 380 inline double fg_max( const double a, const double b ) {
 381     return a < b ? b : a;
 382 }
 383
 384 // return the minimum of the three values
 385 inline double fg_min3( const double a, const double b, const double c ) {
 386     return a > b ? fg_min(b, c) : fg_min(a, c);
 387 }
 388
 389 // return the maximum of the three values
 390 inline double fg_max3 (const double a, const double b, const double c ) {
 391     return a < b ? fg_max(b, c) : fg_max(a, c);
 392 }
 393
 394
 395 // check for an instersection with the individual triangles of a leaf
 396 static bool my_ssg_instersect_leaf( string s, ssgLeaf *leaf, sgMat4 m,
 397                                     const sgVec3 p, const sgVec3 dir,
 398                                     sgVec3 result )
 399 {
 400     sgVec3 v1, v2, n;
 401     sgVec3 p1, p2, p3;
 402     double x, y, z;  // temporary holding spot for result
 403     double a, b, c, d;
 404     double x0, y0, z0, x1, y1, z1, a1, b1, c1;
 405     double t1, t2, t3;
 406     double xmin, xmax, ymin, ymax, zmin, zmax;
 407     double dx, dy, dz, min_dim, x2, y2, x3, y3, rx, ry;
 408     float *tmp;
 409     int side1, side2;
 410     short i1, i2, i3;
 411
 412     // cout << s << "Intersecting" << endl;
 413
 414     // traverse the triangle list for this leaf
 415     for ( int i = 0; i < leaf->getNumTriangles(); ++i ) {
 416         // cout << s << "testing triangle = " << i << endl;
 417
 418         leaf->getTriangle( i, &i1, &i2, &i3 );
 419
 420         // get triangle vertex coordinates
 421
 422         tmp = leaf->getVertex( i1 );
 423         // cout << s << "orig point 1 = " << tmp[0] << " " << tmp[1]
 424         //      << " " << tmp[2] << endl;
 425         sgXformPnt3( p1, tmp, m ) ;
 426
 427         tmp = leaf->getVertex( i2 );
 428         // cout << s << "orig point 2 = " << tmp[0] << " " << tmp[1]
 429         //      << " " << tmp[2] << endl;
 430         sgXformPnt3( p2, tmp, m ) ;
 431
 432         tmp = leaf->getVertex( i3 );
 433         // cout << s << "orig point 3 = " << tmp[0] << " " << tmp[1]
 434         //      << " " << tmp[2] << endl;
 435         sgXformPnt3( p3, tmp, m ) ;
 436
 437         // cout << s << "point 1 = " << p1[0] << " " << p1[1] << " " << p1[2]
 438         //      << endl;
 439         // cout << s << "point 2 = " << p2[0] << " " << p2[1] << " " << p2[2]
 440         //      << endl;
 441         // cout << s << "point 3 = " << p3[0] << " " << p3[1] << " " << p3[2]
 442         //      << endl;
 443
 444         // calculate two edge vectors, and the face normal
 445         sgSubVec3(v1, p2, p1);
 446         sgSubVec3(v2, p3, p1);
 447         sgVectorProductVec3(n, v1, v2);
 448
 449         // calculate the plane coefficients for the plane defined by
 450         // this face.  If n is the normal vector, n = (a, b, c) and p1
 451         // is a point on the plane, p1 = (x0, y0, z0), then the
 452         // equation of the line is a(x-x0) + b(y-y0) + c(z-z0) = 0
 453         a = n[0];
 454         b = n[1];
 455         c = n[2];
 456         d = a * p1[0] + b * p1[1] + c * p1[2];
 457         // printf("a, b, c, d = %.2f %.2f %.2f %.2f\n", a, b, c, d);
 458
 459         // printf("p1(d) = %.2f\n", a * p1[0] + b * p1[1] + c * p1[2]);
 460         // printf("p2(d) = %.2f\n", a * p2[0] + b * p2[1] + c * p2[2]);
 461         // printf("p3(d) = %.2f\n", a * p3[0] + b * p3[1] + c * p3[2]);
 462
 463         // calculate the line coefficients for the specified line
 464         x0 = p[0];  x1 = p[0] + dir[0];
 465         y0 = p[1];  y1 = p[1] + dir[1];
 466         z0 = p[2];  z1 = p[2] + dir[2];
 467
 468         if ( fabs(x1 - x0) > FG_EPSILON ) {
 469             a1 = 1.0 / (x1 - x0);
 470         } else {
 471             // we got a big divide by zero problem here
 472             a1 = 0.0;
 473         }
 474         b1 = y1 - y0;
 475         c1 = z1 - z0;
 476
 477         // intersect the specified line with this plane
 478         t1 = b * b1 * a1;
 479         t2 = c * c1 * a1;
 480
 481         // printf("a = %.2f  t1 = %.2f  t2 = %.2f\n", a, t1, t2);
 482
 483         if ( fabs(a + t1 + t2) > FG_EPSILON ) {
 484             x = (t1*x0 - b*y0 + t2*x0 - c*z0 + d) / (a + t1 + t2);
 485             t3 = a1 * (x - x0);
 486             y = b1 * t3 + y0;
 487             z = c1 * t3 + z0;
 488             // printf("result(d) = %.2f\n", a * x + b * y + c * z);
 489         } else {
 490             // no intersection point
 491             continue;
 492         }
 493
 494 #if 0
 495         if ( side_flag ) {
 496             // check to see if end0 and end1 are on opposite sides of
 497             // plane
 498             if ( (x - x0) > FG_EPSILON ) {
 499                 t1 = x;
 500                 t2 = x0;
 501                 t3 = x1;
 502             } else if ( (y - y0) > FG_EPSILON ) {
 503                 t1 = y;
 504                 t2 = y0;
 505                 t3 = y1;
 506             } else if ( (z - z0) > FG_EPSILON ) {
 507                 t1 = z;
 508                 t2 = z0;
 509                 t3 = z1;
 510             } else {
 511                 // everything is too close together to tell the difference
 512                 // so the current intersection point should work as good
 513                 // as any
 514                 sgSetVec3( result, x, y, z );
 515                 return true;
 516             }
 517             side1 = fg_sign (t1 - t2);
 518             side2 = fg_sign (t1 - t3);
 519             if ( side1 == side2 ) {
 520                 // same side, punt
 521                 continue;
 522             }
 523         }
 524 #endif
 525
 526         // check to see if intersection point is in the bounding
 527         // cube of the face
 528 #ifdef XTRA_DEBUG_STUFF
 529         xmin = fg_min3 (p1[0], p2[0], p3[0]);
 530         xmax = fg_max3 (p1[0], p2[0], p3[0]);
 531         ymin = fg_min3 (p1[1], p2[1], p3[1]);
 532         ymax = fg_max3 (p1[1], p2[1], p3[1]);
 533         zmin = fg_min3 (p1[2], p2[2], p3[2]);
 534         zmax = fg_max3 (p1[2], p2[2], p3[2]);
 535         printf("bounding cube = %.2f,%.2f,%.2f  %.2f,%.2f,%.2f\n",
 536                xmin, ymin, zmin, xmax, ymax, zmax);
 537 #endif
 538         // punt if outside bouding cube
 539         if ( x < (xmin = fg_min3 (p1[0], p2[0], p3[0])) ) {
 540             continue;
 541         } else if ( x > (xmax = fg_max3 (p1[0], p2[0], p3[0])) ) {
 542             continue;
 543         } else if ( y < (ymin = fg_min3 (p1[1], p2[1], p3[1])) ) {
 544             continue;
 545         } else if ( y > (ymax = fg_max3 (p1[1], p2[1], p3[1])) ) {
 546             continue;
 547         } else if ( z < (zmin = fg_min3 (p1[2], p2[2], p3[2])) ) {
 548             continue;
 549         } else if ( z > (zmax = fg_max3 (p1[2], p2[2], p3[2])) ) {
 550             continue;
 551         }
 552
 553         // (finally) check to see if the intersection point is
 554         // actually inside this face
 555
 556         //first, drop the smallest dimension so we only have to work
 557         //in 2d.
 558         dx = xmax - xmin;
 559         dy = ymax - ymin;
 560         dz = zmax - zmin;
 561         min_dim = fg_min3 (dx, dy, dz);
 562         if ( fabs(min_dim - dx) <= FG_EPSILON ) {
 563             // x is the smallest dimension
 564             x1 = p1[1];
 565             y1 = p1[2];
 566             x2 = p2[1];
 567             y2 = p2[2];
 568             x3 = p3[1];
 569             y3 = p3[2];
 570             rx = y;
 571             ry = z;
 572         } else if ( fabs(min_dim - dy) <= FG_EPSILON ) {
 573             // y is the smallest dimension
 574             x1 = p1[0];
 575             y1 = p1[2];
 576             x2 = p2[0];
 577             y2 = p2[2];
 578             x3 = p3[0];
 579             y3 = p3[2];
 580             rx = x;
 581             ry = z;
 582         } else if ( fabs(min_dim - dz) <= FG_EPSILON ) {
 583             // z is the smallest dimension
 584             x1 = p1[0];
 585             y1 = p1[1];
 586             x2 = p2[0];
 587             y2 = p2[1];
 588             x3 = p3[0];
 589             y3 = p3[1];
 590             rx = x;
 591             ry = y;
 592         } else {
 593             // all dimensions are really small so lets call it close
 594             // enough and return a successful match
 595             sgSetVec3( result, x, y, z );
 596             return true;
 597         }
 598
 599         // check if intersection point is on the same side of p1 <-> p2 as p3
 600         t1 = (y1 - y2) / (x1 - x2);
 601         side1 = fg_sign (t1 * ((x3) - x2) + y2 - (y3));
 602         side2 = fg_sign (t1 * ((rx) - x2) + y2 - (ry));
 603         if ( side1 != side2 ) {
 604             // printf("failed side 1 check\n");
 605             continue;
 606         }
 607
 608         // check if intersection point is on correct side of p2 <-> p3 as p1
 609         t1 = (y2 - y3) / (x2 - x3);
 610         side1 = fg_sign (t1 * ((x1) - x3) + y3 - (y1));
 611         side2 = fg_sign (t1 * ((rx) - x3) + y3 - (ry));
 612         if ( side1 != side2 ) {
 613             // printf("failed side 2 check\n");
 614             continue;
 615         }
 616
 617         // check if intersection point is on correct side of p1 <-> p3 as p2
 618         t1 = (y1 - y3) / (x1 - x3);
 619         side1 = fg_sign (t1 * ((x2) - x3) + y3 - (y2));
 620         side2 = fg_sign (t1 * ((rx) - x3) + y3 - (ry));
 621         if ( side1 != side2 ) {
 622             // printf("failed side 3  check\n");
 623             continue;
 624         }
 625
 626         // printf( "intersection point = %.2f %.2f %.2f\n", x, y, z);
 627         sgSetVec3( result, x, y, z );
 628         return true;
 629     }
 630
 631     // printf("\n");
 632
 633     return false;
 634 }
 635
 636
 637 void FGTileMgr::my_ssg_los( string s, ssgBranch *branch, sgMat4 m,
 638                             const sgVec3 p, const sgVec3 dir )
 639 {
 640     sgSphere *bsphere;
 641     for ( ssgEntity *kid = branch->getKid( 0 );
 642           kid != NULL;
 643           kid = branch->getNextKid() )
 644     {
 645         if ( kid->getTraversalMask() & SSGTRAV_HOT ) {
 646             bsphere = kid->getBSphere();
 647             sgVec3 center;
 648             sgCopyVec3( center, bsphere->getCenter() );
 649             sgXformPnt3( center, m ) ;
 650             // cout << s << "entity bounding sphere:" << endl;
 651             // cout << s << "center = " << center[0] << " "
 652             //      << center[1] << " " << center[2] << endl;
 653             // cout << s << "radius = " << bsphere->getRadius() << endl;
 654             double radius_sqd = bsphere->getRadius() * bsphere->getRadius();
 655             if ( sgPointLineDistSquared( center, p, dir ) < radius_sqd ) {
 656                 // possible intersections
 657                 if ( kid->isAKindOf ( ssgTypeBranch() ) ) {
 658                     sgMat4 m_new;
 659                     sgCopyMat4(m_new, m);
 660                     if ( kid->isA( ssgTypeTransform() ) ) {
 661                         sgMat4 xform;
 662                         ((ssgTransform *)kid)->getTransform( xform );
 663                         sgPreMultMat4( m_new, xform );
 664                     }
 665                     my_ssg_los( s + " ", (ssgBranch *)kid, m_new, p, dir );
 666                 } else if ( kid->isAKindOf ( ssgTypeLeaf() ) ) {
 667                     sgVec3 result;
 668                     if ( my_ssg_instersect_leaf( s, (ssgLeaf *)kid, m, p, dir,
 669                                                  result ) )
 670                     {
 671                         cout << "sgLOS hit: " << result[0] << ","
 672                              << result[1] << "," << result[2] << endl;
 673                         hit_pts[hitcount] = result;
 674                         hitcount++;
 675                     }
 676                 }
 677             } else {
 678                 // end of the line for this branch
 679             }
 680         } else {
 681             // branch requested not to be traversed
 682         }
 683     }
 684 }
 685
 686
 687 // Determine scenery altitude via ssg.  Normally this just happens
 688 // when we render the scene, but we'd also like to be able to do this
 689 // explicitely.  lat & lon are in radians.  view_pos in current world
 690 // coordinate translated near (0,0,0) (in meters.)  Returns result in
 691 // meters.
 692 double
 693 FGTileMgr::current_elev_ssg( const Point3D& abs_view_pos,
 694                              const Point3D& view_pos )
 695 {
 696     hitcount = 0;
 697
 698     sgMat4 m;
 699     sgMakeIdentMat4 ( m ) ;
 700
 701     sgVec3 sgavp, sgvp;
 702     sgSetVec3(sgavp, abs_view_pos.x(), abs_view_pos.y(), abs_view_pos.z() );
 703     sgSetVec3(sgvp, view_pos.x(), view_pos.y(), view_pos.z() );
 704
 705     cout << "starting ssg_los, abs view pos = " << abs_view_pos[0] << " "
 706          << abs_view_pos[1] << " " << abs_view_pos[2] << endl;
 707     cout << "starting ssg_los, view pos = " << view_pos[0] << " "
 708          << view_pos[1] << " " << view_pos[2] << endl;
 709     my_ssg_los( "", scene, m, sgvp, sgavp );
 710
 711     double result = -9999;
 712
 713     for ( int i = 0; i < hitcount; ++i ) {
 714         Point3D rel_cart( hit_pts[i][0], hit_pts[i][1], hit_pts[i][2] );
 715         Point3D abs_cart = rel_cart + scenery.center;
 716         Point3D pp = fgCartToPolar3d( abs_cart );
 717         FG_LOG( FG_TERRAIN, FG_INFO, "  polar form = " << pp );
 718         // convert to geodetic coordinates
 719         double lat_geod, alt, sea_level_r;
 720         fgGeocToGeod(pp.lat(), pp.radius(), &lat_geod,
 721                      &alt, &sea_level_r);
 722
 723         // printf("alt = %.2f\n", alt);
 724         // exit since we found an intersection
 725         if ( alt > result && alt < 10000 ) {
 726             // printf("returning alt\n");
 727             result = alt;
 728         }
 729     }
 730
 731     if ( result > -9000 ) {
 732         return result;
 733     } else {
 734         FG_LOG( FG_TERRAIN, FG_INFO, "no terrain intersection" );
 735         return 0.0;
 736     }
 737 }
 738
 739
 740 // given the current lon/lat, fill in the array of local chunks.  If
 741 // the chunk isn't already in the cache, then read it from disk.
 742 int FGTileMgr::update( void ) {
 743     FGTileCache *c;
 744     FGInterface *f;
 745     FGBucket p2;
 746     static FGBucket p_last(false);
 747     static double last_lon = -1000.0;  // in degrees
 748     static double last_lat = -1000.0;  // in degrees
 749     int tile_diameter;
 750     int i, j, dw, dh;
 751
 752     c = &global_tile_cache;
 753     f = current_aircraft.fdm_state;
 754
 755     tile_diameter = current_options.get_tile_diameter();
 756
 757     FGBucket p1( f->get_Longitude() * RAD_TO_DEG,
 758                  f->get_Latitude() * RAD_TO_DEG );
 759     dw = tile_diameter / 2;
 760     dh = tile_diameter / 2;
 761
 762     if ( (p1 == p_last) && (state == Running) ) {
 763         // same bucket as last time
 764         FG_LOG( FG_TERRAIN, FG_DEBUG, "Same bucket as last time" );
 765     } else if ( (state == Start) || (state == Inited) ) {
 766         state = Running;
 767
 768         // First time through or we have teleported, initialize the
 769         // system and load all relavant tiles
 770
 771         FG_LOG( FG_TERRAIN, FG_INFO, "Updating Tile list for " << p1 );
 772         FG_LOG( FG_TERRAIN, FG_INFO, "  First time through ... " );
 773         FG_LOG( FG_TERRAIN, FG_INFO, "  Updating Tile list for " << p1 );
 774         FG_LOG( FG_TERRAIN, FG_INFO, "  Loading "
 775                 << tile_diameter * tile_diameter << " tiles" );
 776
 777         // wipe/initialize tile cache
 778         c->init();
 779         p_last.make_bad();
 780
 781         // build the local area list and schedule tiles for loading
 782
 783         // start with the center tile and work out in concentric
 784         // "rings"
 785
 786         p2 = fgBucketOffset( f->get_Longitude() * RAD_TO_DEG,
 787                              f->get_Latitude() * RAD_TO_DEG,
 788                              0, 0 );
 789         sched_tile( p2 );
 790
 791         // prime scenery center calculations
 792         Point3D geod_center( p2.get_center_lon() * RAD_TO_DEG,
 793                              p2.get_center_lat() * RAD_TO_DEG, 0.0 );
 794         scenery.center = scenery.next_center = fgGeodToCart( geod_center );
 795
 796         Point3D geod_view_center( p2.get_center_lon() * RAD_TO_DEG,
 797                                   p2.get_center_lat() * RAD_TO_DEG,
 798                                   cur_fdm_state->get_Altitude()*FEET_TO_METER +
 799                                   3 );
 800         current_view.abs_view_pos = fgGeodToCart( geod_view_center );
 801         current_view.view_pos = current_view.abs_view_pos - scenery.center;
 802
 803         for ( i = 3; i <= tile_diameter; i = i + 2 ) {
 804             int span = i / 2;
 805
 806             // bottom row
 807             for ( j = -span; j <= span; ++j ) {
 808                 p2 = fgBucketOffset( f->get_Longitude() * RAD_TO_DEG,
 809                                      f->get_Latitude() * RAD_TO_DEG,
 810                                      j, -span );
 811                 sched_tile( p2 );
 812             }
 813
 814             // top row
 815             for ( j = -span; j <= span; ++j ) {
 816                 p2 = fgBucketOffset( f->get_Longitude() * RAD_TO_DEG,
 817                                      f->get_Latitude() * RAD_TO_DEG,
 818                                      j, span );
 819                 sched_tile( p2 );
 820             }
 821
 822             // middle rows
 823             for ( j = -span + 1; j <= span - 1; ++j ) {
 824                 p2 = fgBucketOffset( f->get_Longitude() * RAD_TO_DEG,
 825                                      f->get_Latitude() * RAD_TO_DEG,
 826                                      -span, j );
 827                 sched_tile( p2 );
 828                 p2 = fgBucketOffset( f->get_Longitude() * RAD_TO_DEG,
 829                                      f->get_Latitude() * RAD_TO_DEG,
 830                                      span, j );
 831                 sched_tile( p2 );
 832             }
 833
 834         }
 835
 836         /* for ( j = 0; j < tile_diameter; j++ ) {
 837             for ( i = 0; i < tile_diameter; i++ ) {
 838                 // fgBucketOffset(&p1, &p2, i - dw, j - dh);
 839                 p2 = fgBucketOffset( f->get_Longitude() * RAD_TO_DEG,
 840                                      f->get_Latitude() * RAD_TO_DEG,
 841                                      i - dw, j -dh );
 842                 sched_tile( p2 );
 843             }
 844         } */
 845
 846         // Now force a load of the center tile and inner ring so we
 847         // have something to see in our first frame.
 848         for ( i = 0; i < 9; ++i ) {
 849             if ( load_queue.size() ) {
 850                 FG_LOG( FG_TERRAIN, FG_DEBUG,
 851                         "Load queue not empty, loading a tile" );
 852
 853                 FGLoadRec pending = load_queue.front();
 854                 load_queue.pop_front();
 855                 load_tile( pending.b, pending.cache_index );
 856             }
 857         }
 858
 859     } else {
 860         // We've moved to a new bucket, we need to scroll our
 861         // structures, and load in the new tiles
 862
 863 #if 0
 864         // make sure load queue is flushed before doing shift
 865         while ( load_queue.size() ) {
 866             FG_LOG( FG_TERRAIN, FG_DEBUG,
 867                     "Load queue not empty, flushing queue before tile shift." );
 868
 869             FGLoadRec pending = load_queue.front();
 870             load_queue.pop_front();
 871             load_tile( pending.b, pending.index );
 872         }
 873 #endif
 874
 875         // CURRENTLY THIS ASSUMES WE CAN ONLY MOVE TO ADJACENT TILES.
 876         // AT ULTRA HIGH SPEEDS THIS ASSUMPTION MAY NOT BE VALID IF
 877         // THE AIRCRAFT CAN SKIP A TILE IN A SINGLE ITERATION.
 878
 879         FG_LOG( FG_TERRAIN, FG_INFO, "Updating Tile list for " << p1 );
 880
 881         if ( (p1.get_lon() > p_last.get_lon()) ||
 882              ( (p1.get_lon() == p_last.get_lon()) &&
 883                (p1.get_x() > p_last.get_x()) ) ) {
 884             FG_LOG( FG_TERRAIN, FG_INFO,
 885                     "  (East) Loading " << tile_diameter << " tiles" );
 886             for ( j = 0; j < tile_diameter; j++ ) {
 887                 // scrolling East
 888                 // schedule new column
 889                 p2 = fgBucketOffset( last_lon, last_lat, dw + 1, j - dh );
 890                 sched_tile( p2 );
 891             }
 892         } else if ( (p1.get_lon() < p_last.get_lon()) ||
 893                     ( (p1.get_lon() == p_last.get_lon()) &&
 894                       (p1.get_x() < p_last.get_x()) ) ) {
 895             FG_LOG( FG_TERRAIN, FG_INFO,
 896                     "  (West) Loading " << tile_diameter << " tiles" );
 897             for ( j = 0; j < tile_diameter; j++ ) {
 898                 // scrolling West
 899                 // schedule new column
 900                 p2 = fgBucketOffset( last_lon, last_lat, -dw - 1, j - dh );
 901                 sched_tile( p2 );
 902             }
 903         }
 904
 905         if ( (p1.get_lat() > p_last.get_lat()) ||
 906              ( (p1.get_lat() == p_last.get_lat()) &&
 907                (p1.get_y() > p_last.get_y()) ) ) {
 908             FG_LOG( FG_TERRAIN, FG_INFO,
 909                     "  (North) Loading " << tile_diameter << " tiles" );
 910             for ( i = 0; i < tile_diameter; i++ ) {
 911                 // scrolling North
 912                 // schedule new row
 913                 p2 = fgBucketOffset( last_lon, last_lat, i - dw, dh + 1);
 914                 sched_tile( p2 );
 915             }
 916         } else if ( (p1.get_lat() < p_last.get_lat()) ||
 917                     ( (p1.get_lat() == p_last.get_lat()) &&
 918                       (p1.get_y() < p_last.get_y()) ) ) {
 919             FG_LOG( FG_TERRAIN, FG_INFO,
 920                     "  (South) Loading " << tile_diameter << " tiles" );
 921             for ( i = 0; i < tile_diameter; i++ ) {
 922                 // scrolling South
 923                 // schedule new row
 924                 p2 = fgBucketOffset( last_lon, last_lat, i - dw, -dh - 1);
 925                 sched_tile( p2 );
 926             }
 927         }
 928     }
 929
 930     if ( load_queue.size() ) {
 931         FG_LOG( FG_TERRAIN, FG_DEBUG, "Load queue not empty, loading a tile" );
 932
 933         FGLoadRec pending = load_queue.front();
 934         load_queue.pop_front();
 935         load_tile( pending.b, pending.cache_index );
 936     }
 937
 938     // find our current elevation (feed in the current bucket to save work)
 939     Point3D geod_pos = Point3D( f->get_Longitude(), f->get_Latitude(), 0.0);
 940     Point3D tmp_abs_view_pos = fgGeodToCart(geod_pos);
 941
 942     cout << "current elevation (old) == "
 943          << current_elev( f->get_Longitude(), f->get_Latitude(),
 944                           tmp_abs_view_pos )
 945          << endl;
 946     scenery.cur_elev = current_elev_ssg( current_view.abs_view_pos,
 947                                          current_view.view_pos );
 948     cout << "current elevation (ssg) == " << scenery.cur_elev << endl;
 949
 950     p_last = p1;
 951     last_lon = f->get_Longitude() * RAD_TO_DEG;
 952     last_lat = f->get_Latitude() * RAD_TO_DEG;
 953
 954     return 1;
 955 }
 956
 957
 958 // NEW
 959
 960 // inrange() IS THIS POINT WITHIN POSSIBLE VIEWING RANGE ?
 961 //      calculate distance from vertical tangent line at
 962 //      current position to center of object.
 963 //      this is equivalent to
 964 //      dist = point_line_dist_squared( &(t->center), &(v->abs_view_pos),
 965 //                                      v->local_up );
 966 //      if ( dist < FG_SQUARE(t->bounding_radius) ) {
 967 //
 968 // the compiler should inline this for us
 969
 970 static int
 971 inrange( const double radius, const Point3D& center, const Point3D& vp,
 972          const MAT3vec up)
 973 {
 974     MAT3vec u, u1, v;
 975     //  double tmp;
 976
 977     // u = p - p0
 978     u[0] = center.x() - vp.x();
 979     u[1] = center.y() - vp.y();
 980     u[2] = center.z() - vp.z();
 981
 982     // calculate the projection, u1, of u along d.
 983     // u1 = ( dot_prod(u, d) / dot_prod(d, d) ) * d;
 984
 985     MAT3_SCALE_VEC(u1, up,
 986                    (MAT3_DOT_PRODUCT(u, up) / MAT3_DOT_PRODUCT(up, up)) );
 987
 988     // v = u - u1 = vector from closest point on line, p1, to the
 989     // original point, p.
 990     MAT3_SUB_VEC(v, u, u1);
 991
 992     return( FG_SQUARE(radius) >= MAT3_DOT_PRODUCT(v, v));
 993 }
 994
 995
 996 // NEW for legibility
 997
 998 // update this tile's geometry for current view
 999 // The Compiler should inline this
1000 static void
1001 update_tile_geometry( FGTileEntry *t, GLdouble *MODEL_VIEW)
1002 {
1003     GLfloat *m;
1004     double x, y, z;
1005
1006     // calculate tile offset
1007     t->offset = t->center - scenery.center;
1008
1009     x = t->offset.x();
1010     y = t->offset.y();
1011     z = t->offset.z();
1012
1013     m = t->model_view;
1014
1015     // Calculate the model_view transformation matrix for this tile
1016     FG_MEM_COPY( m, MODEL_VIEW, 16*sizeof(GLdouble) );
1017
1018     // This is equivalent to doing a glTranslatef(x, y, z);
1019     m[12] += (m[0]*x + m[4]*y + m[8] *z);
1020     m[13] += (m[1]*x + m[5]*y + m[9] *z);
1021     m[14] += (m[2]*x + m[6]*y + m[10]*z);
1022     // m[15] += (m[3]*x + m[7]*y + m[11]*z);
1023     // m[3] m7[] m[11] are 0.0 see LookAt() in views.cxx
1024     // so m[15] is unchanged
1025 }
1026
1027
1028 // Prepare the ssg nodes ... for each tile, set it's proper
1029 // transform and update it's range selector based on current
1030 // visibilty
1031 void FGTileMgr::prep_ssg_nodes( void ) {
1032     FGTileEntry *t;
1033
1034     float ranges[2];
1035     ranges[0] = 0.0f;
1036
1037     // traverse the potentially viewable tile list and update range
1038     // selector and transform
1039     for ( int i = 0; i < (int)global_tile_cache.get_size(); i++ ) {
1040         t = global_tile_cache.get_tile( i );
1041
1042         if ( t->is_loaded() ) {
1043             // set range selector (LOD trick) to be distance to center
1044             // of tile + bounding radius
1045 #ifndef FG_OLD_WEATHER
1046             ranges[1] = WeatherDatabase->getWeatherVisibility()
1047                 + t->bounding_radius;
1048 #else
1049             ranges[1] = current_weather.get_visibility()+t->bounding_radius;
1050 #endif
1051             t->range_ptr->setRanges( ranges, 2 );
1052
1053             // calculate tile offset
1054             t->SetOffset( scenery.center );
1055
1056             // calculate ssg transform
1057             sgCoord sgcoord;
1058             sgSetCoord( &sgcoord,
1059                         t->offset.x(), t->offset.y(), t->offset.z(),
1060                         0.0, 0.0, 0.0 );
1061             t->transform_ptr->setTransform( &sgcoord );
1062         }
1063     }
1064 }