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