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