src/Scenery/hitlist.cxx

   1 // hitlist.cxx -
   2 // Height Over Terrain and Assosciated Routines for FlightGear based Scenery
   3 // Written by Norman Vine, started 2000.
   4
   5 #ifdef HAVE_CONFIG_H
   6 #  include <config.h>
   7 #endif
   8
   9 #include <float.h>
  10 #include <math.h>
  11
  12 #include <plib/sg.h>
  13 #include <plib/ssg.h>
  14
  15 #include <simgear/sg_inlines.h>
  16 #include <simgear/debug/logstream.hxx>
  17 #include <simgear/math/point3d.hxx>
  18 #include <simgear/math/sg_geodesy.hxx>
  19 #include <simgear/math/vector.hxx>
  20 #include <simgear/timing/timestamp.hxx>
  21
  22 #include <Main/globals.hxx>
  23 #include <Main/viewer.hxx>
  24 #include <Scenery/scenery.hxx>
  25
  26 #include "hitlist.hxx"
  27
  28 // Specialized version of sgMultMat4 needed because of mixed matrix
  29 // types
  30 static inline void sgMultMat4(sgdMat4 dst, sgdMat4 m1, sgMat4 m2) {
  31     for ( int j = 0 ; j < 4 ; j++ ) {
  32         dst[0][j] = m2[0][0] * m1[0][j] +
  33                     m2[0][1] * m1[1][j] +
  34                     m2[0][2] * m1[2][j] +
  35                     m2[0][3] * m1[3][j] ;
  36
  37         dst[1][j] = m2[1][0] * m1[0][j] +
  38                     m2[1][1] * m1[1][j] +
  39                     m2[1][2] * m1[2][j] +
  40                     m2[1][3] * m1[3][j] ;
  41
  42         dst[2][j] = m2[2][0] * m1[0][j] +
  43                     m2[2][1] * m1[1][j] +
  44                     m2[2][2] * m1[2][j] +
  45                     m2[2][3] * m1[3][j] ;
  46
  47         dst[3][j] = m2[3][0] * m1[0][j] +
  48                     m2[3][1] * m1[1][j] +
  49                     m2[3][2] * m1[2][j] +
  50                     m2[3][3] * m1[3][j] ;
  51     }
  52 }
  53
  54
  55 /*
  56  * Walk backwards up the tree, transforming the vertex by all the
  57  * matrices along the way.
  58  *
  59  * Upwards recursion hurts my head.
  60  */
  61 static void ssgGetEntityTransform(ssgEntity *entity, sgMat4 m ) {
  62     sgMat4 mat ;
  63
  64     // If this node has a parent - get the composite matrix for the
  65     // parent.
  66     if ( entity->getNumParents() > 0 )
  67         ssgGetEntityTransform ( entity->getParent(0), mat ) ;
  68     else
  69         sgMakeIdentMat4 ( mat ) ;
  70
  71     // If this node has a transform - then concatenate it.
  72     if ( entity -> isAKindOf ( ssgTypeTransform () ) ) {
  73         sgMat4 this_mat ;
  74         ((ssgTransform *) entity) -> getTransform ( this_mat ) ;
  75         sgPostMultMat4 ( mat, this_mat ) ;
  76     }
  77
  78     sgCopyMat4 ( m, mat ) ;
  79 }
  80
  81
  82 // return the passed entitity's bsphere's center point radius and
  83 // fully formed current model matrix for entity
  84 static inline void ssgGetCurrentBSphere( ssgEntity *entity, sgVec3 center,
  85                                   float *radius, sgMat4 m )
  86 {
  87     sgSphere *bsphere = entity->getBSphere();
  88     *radius = (double)bsphere->getRadius();
  89     sgCopyVec3( center, bsphere->getCenter() );
  90     sgMakeIdentMat4 ( m ) ;
  91     ssgGetEntityTransform( entity, m );
  92 }
  93
  94
  95 // This is same as PLib's sgdIsectInfLinePlane() and can be replaced
  96 // by it after the next PLib release
  97 static inline bool fgdIsectInfLinePlane( sgdVec3 dst,
  98                                          const sgdVec3 l_org,
  99                                          const sgdVec3 l_vec,
 100                                          const sgdVec4 plane )
 101 {
 102     SGDfloat tmp = sgdScalarProductVec3 ( l_vec, plane ) ;
 103
 104     /* Is line parallel to plane? */
 105
 106     if ( fabs ( tmp ) < FLT_EPSILON )
 107         return false ;
 108
 109     sgdScaleVec3 ( dst, l_vec, -( sgdScalarProductVec3 ( l_org, plane )
 110                                   + plane[3] ) / tmp ) ;
 111     sgdAddVec3  ( dst, l_org ) ;
 112
 113     return true ;
 114 }
 115
 116 /*
 117  * Given a point and a triangle lying on the same plane check to see
 118  * if the point is inside the triangle
 119  *
 120  * This is same as PLib's sgdPointInTriangle() and can be replaced by
 121  * it after the next PLib release
 122  */
 123 static inline bool fgdPointInTriangle( sgdVec3 point, sgdVec3 tri[3] )
 124 {
 125     sgdVec3 dif;
 126
 127     SGDfloat min, max;
 128     // punt if outside bouding cube
 129     SG_MIN_MAX3 ( min, max, tri[0][0], tri[1][0], tri[2][0] );
 130     if( (point[0] < min) || (point[0] > max) )
 131         return false;
 132     dif[0] = max - min;
 133
 134     SG_MIN_MAX3 ( min, max, tri[0][1], tri[1][1], tri[2][1] );
 135     if( (point[1] < min) || (point[1] > max) )
 136         return false;
 137     dif[1] = max - min;
 138
 139     SG_MIN_MAX3 ( min, max, tri[0][2], tri[1][2], tri[2][2] );
 140     if( (point[2] < min) || (point[2] > max) )
 141         return false;
 142     dif[2] = max - min;
 143
 144     // drop the smallest dimension so we only have to work in 2d.
 145     SGDfloat min_dim = SG_MIN3 (dif[0], dif[1], dif[2]);
 146     SGDfloat x1, y1, x2, y2, x3, y3, rx, ry;
 147     if ( fabs(min_dim-dif[0]) <= DBL_EPSILON ) {
 148         // x is the smallest dimension
 149         x1 = point[1];
 150         y1 = point[2];
 151         x2 = tri[0][1];
 152         y2 = tri[0][2];
 153         x3 = tri[1][1];
 154         y3 = tri[1][2];
 155         rx = tri[2][1];
 156         ry = tri[2][2];
 157     } else if ( fabs(min_dim-dif[1]) <= DBL_EPSILON ) {
 158         // y is the smallest dimension
 159         x1 = point[0];
 160         y1 = point[2];
 161         x2 = tri[0][0];
 162         y2 = tri[0][2];
 163         x3 = tri[1][0];
 164         y3 = tri[1][2];
 165         rx = tri[2][0];
 166         ry = tri[2][2];
 167     } else if ( fabs(min_dim-dif[2]) <= DBL_EPSILON ) {
 168         // z is the smallest dimension
 169         x1 = point[0];
 170         y1 = point[1];
 171         x2 = tri[0][0];
 172         y2 = tri[0][1];
 173         x3 = tri[1][0];
 174         y3 = tri[1][1];
 175         rx = tri[2][0];
 176         ry = tri[2][1];
 177     } else {
 178         // all dimensions are really small so lets call it close
 179         // enough and return a successful match
 180         return true;
 181     }
 182
 183     // check if intersection point is on the same side of p1 <-> p2 as p3
 184     SGDfloat tmp = (y2 - y3) / (x2 - x3);
 185     int side1 = SG_SIGN (tmp * (rx - x3) + y3 - ry);
 186     int side2 = SG_SIGN (tmp * (x1 - x3) + y3 - y1);
 187     if ( side1 != side2 ) {
 188         // printf("failed side 1 check\n");
 189         return false;
 190     }
 191
 192     // check if intersection point is on correct side of p2 <-> p3 as p1
 193     tmp = (y3 - ry) / (x3 - rx);
 194     side1 = SG_SIGN (tmp * (x2 - rx) + ry - y2);
 195     side2 = SG_SIGN (tmp * (x1 - rx) + ry - y1);
 196     if ( side1 != side2 ) {
 197         // printf("failed side 2 check\n");
 198         return false;
 199     }
 200
 201     // check if intersection point is on correct side of p1 <-> p3 as p2
 202     tmp = (y2 - ry) / (x2 - rx);
 203     side1 = SG_SIGN (tmp * (x3 - rx) + ry - y3);
 204     side2 = SG_SIGN (tmp * (x1 - rx) + ry - y1);
 205     if ( side1 != side2 ) {
 206         // printf("failed side 3  check\n");
 207         return false;
 208     }
 209
 210     return true;
 211 }
 212
 213
 214 // Check if all three vertices are the same point (or close enough)
 215 static inline int isZeroAreaTri( sgdVec3 tri[3] )
 216 {
 217     return( sgdEqualVec3(tri[0], tri[1]) ||
 218             sgdEqualVec3(tri[1], tri[2]) ||
 219             sgdEqualVec3(tri[2], tri[0]) );
 220 }
 221
 222
 223 // Constructor
 224 FGHitList::FGHitList() :
 225     last(NULL), test_dist(DBL_MAX)
 226 {
 227 }
 228
 229
 230 // Destructor
 231 FGHitList::~FGHitList() {}
 232
 233
 234 /*
 235 Find the intersection of an infinite line with a leaf the line being
 236 defined by a point and direction.
 237
 238 Variables
 239 In:
 240 ssgLeaf pointer  -- leaf
 241 qualified matrix -- m
 242 line origin      -- orig
 243 line direction   -- dir
 244 Out:
 245 result  -- intersection point
 246 normal  -- intersected tri's normal
 247
 248 Returns:
 249 true if intersection found
 250 false otherwise
 251
 252 !!! WARNING !!!
 253
 254 If you need an exhaustive list of hitpoints YOU MUST use the generic
 255 version of this function as the specialized versions will do an early
 256 out of expensive tests if the point can not be the closest one found
 257
 258 !!! WARNING !!!
 259 */
 260 int FGHitList::IntersectLeaf( ssgLeaf *leaf, sgdMat4 m,
 261                               sgdVec3 orig, sgdVec3 dir )
 262 {
 263     int num_hits = 0;
 264     int i = 0;
 265
 266     for ( ; i < leaf->getNumTriangles(); ++i ) {
 267         short i1, i2, i3;
 268         leaf->getTriangle( i, &i1, &i2, &i3 );
 269
 270         sgdVec3 tri[3];
 271         sgdSetVec3( tri[0], leaf->getVertex( i1 ) );
 272         sgdSetVec3( tri[1], leaf->getVertex( i2 ) );
 273         sgdSetVec3( tri[2], leaf->getVertex( i3 ) );
 274
 275         if( isZeroAreaTri( tri ) )
 276             continue;
 277
 278         sgdVec4 plane;
 279         sgdMakePlane( plane, tri[0], tri[1], tri[2] );
 280
 281         sgdVec3 point;
 282         if( fgdIsectInfLinePlane( point, orig, dir, plane ) ) {
 283             if( fgdPointInTriangle( point, tri ) ) {
 284                 // transform point into passed into desired coordinate frame
 285                 sgdXformPnt3( point, point, m );
 286                 sgdXformPnt4(plane,plane,m);
 287                 add(leaf,i,point,plane);
 288                 num_hits++;
 289             }
 290         }
 291     }
 292     return num_hits;
 293 }
 294
 295
 296 // Short circuit/slightly optimized version of the full IntersectLeaf()
 297 int FGHitList::IntersectLeaf( ssgLeaf *leaf, sgdMat4 m,
 298                               sgdVec3 orig, sgdVec3 dir,
 299                               GLenum primType )
 300 {
 301     double tmp_dist;
 302
 303     // number of hits but there could be more that
 304     // were not found because of short circut switch !
 305     // so you may want to use the unspecialized IntersectLeaf()
 306     int n, num_hits = 0;
 307
 308     int ntri = leaf->getNumTriangles();
 309     for ( n = 0; n < ntri; ++n )
 310     {
 311         sgdVec3 tri[3];
 312
 313         switch ( primType )
 314         {
 315         case GL_POLYGON :
 316             SG_LOG( SG_TERRAIN, SG_ALERT,
 317                     "WARNING: dubiously handled GL_POLYGON" );
 318         case GL_TRIANGLE_FAN :
 319             /* SG_LOG( SG_TERRAIN, SG_ALERT,
 320                "IntersectLeaf: GL_TRIANGLE_FAN" ); */
 321             if ( !n ) {
 322                 sgdSetVec3( tri[0], leaf->getVertex( short(0) ) );
 323                 sgdSetVec3( tri[1], leaf->getVertex( short(1) ) );
 324                 sgdSetVec3( tri[2], leaf->getVertex( short(2) ) );
 325             } else {
 326                 sgdCopyVec3( tri[1], tri[2] );
 327                 sgdSetVec3( tri[2], leaf->getVertex( short(n+2) ) );
 328             }
 329             break;
 330         case GL_TRIANGLES :
 331             /* SG_LOG( SG_TERRAIN, SG_DEBUG,
 332                "IntersectLeaf: GL_TRIANGLES" ); */
 333             sgdSetVec3( tri[0], leaf->getVertex( short(n*3) ) );
 334             sgdSetVec3( tri[1], leaf->getVertex( short(n*3+1) ) );
 335             sgdSetVec3( tri[2], leaf->getVertex( short(n*3+2) ) );
 336             break;
 337         case GL_QUAD_STRIP :
 338             SG_LOG( SG_TERRAIN, SG_ALERT,
 339                     "WARNING: dubiously handled GL_QUAD_STRIP" );
 340         case GL_TRIANGLE_STRIP :
 341             /* SG_LOG( SG_TERRAIN, SG_ALERT,
 342                "IntersectLeaf: GL_TRIANGLE_STRIP" ); */
 343             if ( !n ) {
 344                 sgdSetVec3( tri[0], leaf->getVertex( short(0) ) );
 345                 sgdSetVec3( tri[1], leaf->getVertex( short(1) ) );
 346                 sgdSetVec3( tri[2], leaf->getVertex( short(2) ) );
 347             } else {
 348                 if ( n & 1 ) {
 349                     sgdCopyVec3( tri[0], tri[2] );
 350                     sgdSetVec3( tri[2], leaf->getVertex( short(n+2) ) );
 351                 } else {
 352                     sgdCopyVec3( tri[1], tri[2] );
 353                     sgdSetVec3( tri[2], leaf->getVertex( short(n+2) ) );
 354                 }
 355             }
 356             break;
 357         case GL_QUADS :
 358             SG_LOG( SG_TERRAIN, SG_ALERT,
 359                     "WARNING: dubiously handled GL_QUADS" );
 360             sgdSetVec3( tri[0], leaf->getVertex( short(n*2) ) );
 361             sgdSetVec3( tri[1], leaf->getVertex( short(n*2+1) ) );
 362             sgdSetVec3( tri[2], leaf->getVertex( short(n*2 + 2 - (n&1)*4) ) );
 363             break;
 364         default:
 365             SG_LOG( SG_TERRAIN, SG_ALERT,
 366                     "WARNING: not-handled structure: " << primType );
 367             return IntersectLeaf( leaf, m, orig, dir);
 368         }
 369
 370         if( isZeroAreaTri( tri ) )
 371             continue;
 372
 373         sgdVec4 plane;
 374         sgdMakePlane( plane, tri[0], tri[1], tri[2] );
 375
 376         sgdVec3 point;
 377
 378         // find point of intersection of line from point org
 379         // in direction dir with triangle's plane
 380         SGDfloat tmp = sgdScalarProductVec3 ( dir, plane ) ;
 381         /* Is line parallel to plane? */
 382         if ( sgdAbs ( tmp ) < FLT_EPSILON /*DBL_EPSILON*/ )
 383             continue ;
 384
 385         // find parametric point
 386         sgdScaleVec3 ( point, dir,
 387                        -( sgdScalarProductVec3 ( orig, plane ) + plane[3] )
 388                        / tmp ) ;
 389
 390         // short circut if this point is further away then a previous hit
 391         tmp_dist = sgdDistanceSquaredVec3(point, orig );
 392         if( tmp_dist > test_dist )
 393             continue;
 394
 395         // place parametric point in world
 396         sgdAddVec3 ( point, orig ) ;
 397
 398         if( fgdPointInTriangle( point, tri ) ) {
 399             // transform point into passed coordinate frame
 400             sgdXformPnt3( point, point, m );
 401             sgdXformPnt4(plane,plane,m);
 402             add(leaf,n,point,plane);
 403             test_dist = tmp_dist;
 404             num_hits++;
 405         }
 406     }
 407     return num_hits;
 408 }
 409
 410
 411
 412 inline static bool IN_RANGE( sgdVec3 v, double radius ) {
 413     return ( sgdScalarProductVec3(v, v) < (radius*radius) );
 414 }
 415
 416
 417 void FGHitList::IntersectBranch( ssgBranch *branch, sgdMat4 m,
 418                                  sgdVec3 orig, sgdVec3 dir )
 419 {
 420     /* the lookat vector and matrix in branch's coordinate frame but
 421      * we won't determine these unless needed, This 'lazy evaluation'
 422      * is a result of profiling data */
 423     sgdVec3 orig_leaf, dir_leaf;
 424     sgdMat4 m_leaf;
 425
 426     // 'lazy evaluation' flag
 427     int first_time = 1;
 428
 429     for ( ssgEntity *kid = branch->getKid( 0 );
 430             kid != NULL;
 431             kid = branch->getNextKid() )
 432     {
 433         if ( kid->getTraversalMask() & SSGTRAV_HOT
 434              && !kid->getBSphere()->isEmpty() )
 435         {
 436             sgdVec3 center;
 437             sgdSetVec3( center,
 438                         kid->getBSphere()->getCenter()[0],
 439                         kid->getBSphere()->getCenter()[1],
 440                         kid->getBSphere()->getCenter()[2] );
 441             sgdXformPnt3( center, m ) ;
 442
 443             // sgdClosestPointToLineDistSquared( center, orig, dir )
 444             // inlined here because because of profiling results
 445             sgdVec3 u, u1, v;
 446             sgdSubVec3(u, center, orig);
 447             sgdScaleVec3( u1, dir, sgdScalarProductVec3(u,dir)
 448                           / sgdScalarProductVec3(dir,dir) );
 449             sgdSubVec3(v, u, u1);
 450
 451             // double because of possible overflow
 452             if ( IN_RANGE( v, double(kid->getBSphere()->getRadius()) ) )
 453             {
 454                 if ( kid->isAKindOf ( ssgTypeBranch() ) )
 455                 {
 456                     sgdMat4 m_new;
 457                     if ( kid->isA( ssgTypeTransform() ) )
 458                     {
 459                         sgMat4 fxform;
 460                         ((ssgTransform *)kid)->getTransform( fxform );
 461                         sgMultMat4(m_new, m, fxform);
 462                     } else {
 463                         sgdCopyMat4(m_new, m);
 464                     }
 465                     IntersectBranch( (ssgBranch *)kid, m_new, orig, dir );
 466                 }
 467                 else if ( kid->isAKindOf( ssgTypeLeaf() ) )
 468                 {
 469                     if ( first_time ) {
 470                         sgdTransposeNegateMat4( m_leaf, m );
 471                         sgdXformPnt3( orig_leaf, orig, m_leaf );
 472                         sgdXformPnt3( dir_leaf,  dir,  m_leaf );
 473                         first_time = 0;
 474                     }
 475                     // GLenum primType = ((ssgLeaf *)kid)->getPrimitiveType();
 476                     // IntersectLeaf( (ssgLeaf *)kid, m, orig_leaf, dir_leaf,
 477                     //                primType );
 478                     IntersectLeaf( (ssgLeaf *)kid, m, orig_leaf, dir_leaf );
 479                 }
 480             }  // Out of range
 481         } // branch not requested to be traversed
 482     } // end for loop
 483 }
 484
 485
 486 // a temporary hack until we get everything rewritten with sgdVec3
 487 static inline Point3D operator + (const Point3D& a, const sgdVec3 b)
 488 {
 489     return Point3D(a.x()+b[0], a.y()+b[1], a.z()+b[2]);
 490 }
 491
 492
 493 void FGHitList::Intersect( ssgBranch *scene, sgdVec3 orig, sgdVec3 dir ) {
 494     sgdMat4 m;
 495     clear();
 496     sgdMakeIdentMat4 ( m ) ;
 497     IntersectBranch( scene, m, orig, dir );
 498 }
 499
 500
 501 void FGHitList::Intersect( ssgBranch *scene, sgdMat4 m, sgdVec3 orig, sgdVec3 dir )
 502 {
 503     clear();
 504     IntersectBranch( scene, m, orig, dir );
 505 }
 506
 507
 508 // Determine scenery altitude via ssg.
 509 // returned results are in meters
 510 // static double hitlist1_time = 0.0;
 511
 512 bool fgCurrentElev( sgdVec3 abs_view_pos, double max_alt_m,
 513                     sgdVec3 scenery_center,
 514                     FGHitList *hit_list,
 515                     double *terrain_elev, double *radius, double *normal)
 516 {
 517     // SGTimeStamp start; start.stamp();
 518
 519     bool result;
 520     sgdVec3 view_pos;
 521     sgdSubVec3( view_pos, abs_view_pos, scenery_center );
 522
 523     sgdVec3 orig, dir;
 524     sgdCopyVec3(orig, view_pos );
 525     sgdCopyVec3(dir, abs_view_pos );
 526
 527     hit_list->Intersect( globals->get_scenery()->get_terrain_branch(),
 528                          orig, dir );
 529
 530     int this_hit = -1;
 531     int max_hit = -1;
 532     Point3D geoc;
 533     double hit_elev = -9999;
 534     double max_elev = -9999;
 535     Point3D sc(scenery_center[0], scenery_center[1], scenery_center[2]) ;
 536
 537     int hitcount = hit_list->num_hits();
 538     // cout << "hits = " << hitcount << endl;
 539     for ( int i = 0; i < hitcount; ++i ) {
 540         // FIXME: sgCartToGeod is slow.  Call it just once for the
 541         // "sc" point, and then handle the rest with a geodetic "up"
 542         // vector approximation.  Across one tile, this will be
 543         // acceptable.
 544         double alt = sgCartToGeod( sc + hit_list->get_point(i) ).elev();
 545         // cout << "hit " << i << " lon = " << geoc.lon() << " lat = "
 546         //      << lat_geod << " alt = " << alt << "  max alt = " << max_alt_m
 547         //      << endl;
 548         if ( alt > hit_elev && alt < max_alt_m ) {
 549             // cout << "  it's a keeper" << endl;
 550             hit_elev = alt;
 551             this_hit = i;
 552         }
 553         if ( alt > hit_elev ) {
 554             max_elev = alt;
 555             max_hit = i;
 556         }
 557     }
 558
 559     if ( this_hit < 0 ) {
 560         // no hits below us, take the max hit
 561         this_hit = max_hit;
 562         hit_elev = max_elev;
 563     }
 564
 565     if ( hit_elev > -9000 ) {
 566         *terrain_elev = hit_elev;
 567         *radius = geoc.radius();
 568         sgVec3 tmp;
 569         sgSetVec3(tmp, hit_list->get_normal(this_hit));
 570         // cout << "cur_normal: " << tmp[0] << " " << tmp[1] << " "  << tmp[2] << endl;
 571         sgdSetVec3( normal, tmp );
 572         // float *up = globals->get_current_view()->get_world_up();
 573         // cout << "world_up  : " << up[0] << " " << up[1] << " " << up[2] << endl;
 574         /* ssgState *IntersectedLeafState =
 575               ((ssgLeaf*)hit_list->get_entity(this_hit))->getState(); */
 576         result = true;
 577     } else {
 578         SG_LOG( SG_TERRAIN, SG_INFO, "no terrain intersection" );
 579         *terrain_elev = 0.0;
 580         float *up = globals->get_current_view()->get_world_up();
 581         sgdSetVec3(normal, up[0], up[1], up[2]);
 582         result = false;
 583     }
 584
 585     // SGTimeStamp finish; finish.stamp();
 586     // hitlist1_time = ( 29.0 * hitlist1_time + (finish - start) ) / 30.0;
 587     // cout << " time per call = " << hitlist1_time << endl;
 588
 589     return result;
 590 }
 591
 592
 593 // static double hitlist2_time = 0.0;
 594
 595 // Determine scenery altitude via ssg.
 596 // returned results are in meters
 597 bool fgCurrentElev( sgdVec3 abs_view_pos, double max_alt_m,
 598                     sgdVec3 scenery_center,
 599                     ssgTransform *terra_transform,
 600                     FGHitList *hit_list,
 601                     double *terrain_elev, double *radius, double *normal)
 602 {
 603     // SGTimeStamp start; start.stamp();
 604
 605     bool result;
 606     sgdVec3 view_pos;
 607     sgdSubVec3( view_pos, abs_view_pos, scenery_center );
 608
 609     sgdVec3 orig, dir;
 610     sgdCopyVec3(orig, view_pos );
 611
 612     sgdCopyVec3(dir, abs_view_pos );
 613     sgdNormalizeVec3(dir);
 614
 615     sgMat4 fxform;
 616     sgMakeIdentMat4 ( fxform ) ;
 617     ssgGetEntityTransform( terra_transform, fxform );
 618
 619     sgdMat4 xform;
 620     sgdSetMat4(xform,fxform);
 621     hit_list->Intersect( terra_transform, xform, orig, dir );
 622
 623     int this_hit = -1;
 624     int max_hit = -1;
 625     Point3D geoc;
 626     double hit_elev = -9999;
 627     double max_elev = -9999;
 628     Point3D sc(scenery_center[0], scenery_center[1], scenery_center[2]) ;
 629
 630     int hitcount = hit_list->num_hits();
 631     // cout << "hits = " << hitcount << endl;
 632     for ( int i = 0; i < hitcount; ++i ) {
 633         // FIXME: sgCartToGeod is slow.  Call it just once for the
 634         // "sc" point, and then handle the rest with a geodetic "up"
 635         // vector approximation.  Across one tile, this will be
 636         // acceptable.
 637         double alt = sgCartToGeod( sc + hit_list->get_point(i) ).elev();
 638         // cout << "hit " << i << " lon = " << geoc.lon() << " lat = "
 639         //      << lat_geod << " alt = " << alt << "  max alt = " << max_alt_m
 640         //      << endl;
 641         if ( alt > hit_elev && alt < max_alt_m ) {
 642             hit_elev = alt;
 643             this_hit = i;
 644             // cout << "  it's a keeper" << endl;
 645         }
 646         if ( alt > hit_elev ) {
 647             max_elev = alt;
 648             max_hit = i;
 649         }
 650     }
 651
 652     if ( this_hit < 0 ) {
 653         // no hits below us, take the max hit
 654         this_hit = max_hit;
 655         hit_elev = max_elev;
 656     }
 657
 658     if ( hit_elev > -9000 ) {
 659         *terrain_elev = hit_elev;
 660         *radius = geoc.radius();
 661         sgVec3 tmp;
 662         sgSetVec3(tmp, hit_list->get_normal(this_hit));
 663         // cout << "cur_normal: " << tmp[0] << " " << tmp[1] << " "  << tmp[2] << endl;
 664         sgdSetVec3( normal, tmp );
 665         // float *up = globals->get_current_view()->get_world_up();
 666        // cout << "world_up  : " << up[0] << " " << up[1] << " " << up[2] << endl;
 667         /* ssgState *IntersectedLeafState =
 668               ((ssgLeaf*)hit_list->get_entity(this_hit))->getState(); */
 669         result = true;
 670     } else {
 671         SG_LOG( SG_TERRAIN, SG_DEBUG, "DOING FULL TERRAIN INTERSECTION" );
 672         result = fgCurrentElev( abs_view_pos, max_alt_m, scenery_center,
 673                                 hit_list, terrain_elev, radius, normal);
 674     }
 675
 676     // SGTimeStamp finish; finish.stamp();
 677     // hitlist2_time = ( 29.0 * hitlist2_time + (finish - start) ) / 30.0;
 678     // cout << "time per call 2 = " << hitlist2_time << endl;
 679
 680     return result;
 681 }
 682