osgEarth 2.1.1

/home/cube/sources/osgearth/src/osgEarth/ElevationQuery.cpp

Go to the documentation of this file.
00001 #include <osgEarth/ElevationQuery>
00002 #include <osgEarth/Locators>
00003 #include <osgTerrain/TerrainTile>
00004 #include <osgTerrain/GeometryTechnique>
00005 #include <osgUtil/IntersectionVisitor>
00006 #include <osgUtil/LineSegmentIntersector>
00007 
00008 #define LC "[ElevationQuery] "
00009 
00010 using namespace osgEarth;
00011 using namespace OpenThreads;
00012 
00013 ElevationQuery::ElevationQuery( const Map* map ) :
00014 _mapf( map, Map::ELEVATION_LAYERS )
00015 {
00016     postCTOR();
00017 }
00018 
00019 ElevationQuery::ElevationQuery( const MapFrame& mapFrame ) :
00020 _mapf( mapFrame )
00021 {
00022     postCTOR();
00023 }
00024 
00025 void
00026 ElevationQuery::postCTOR()
00027 {
00028     _tileSize         = 0;
00029     _maxDataLevel     = 0;
00030     _technique        = TECHNIQUE_PARAMETRIC;
00031     _interpolation    = INTERP_BILINEAR;
00032     _maxLevelOverride = -1;
00033 
00034     // Limit the size of the cache we'll use to cache heightfields. This is an
00035     // LRU cache.
00036     _tileCache.setMaxSize( 50 );
00037 }
00038 
00039 void
00040 ElevationQuery::sync()
00041 {
00042     if ( _mapf.sync() || _tileSize == 0 || _maxDataLevel == 0 )
00043     {
00044         _tileSize = 0;
00045         _maxDataLevel = 0;
00046 
00047         for( ElevationLayerVector::const_iterator i = _mapf.elevationLayers().begin(); i != _mapf.elevationLayers().end(); ++i )
00048         {
00049             // we need the maximum tile size
00050             int layerTileSize = i->get()->getTileSize();
00051             if ( layerTileSize > _tileSize )
00052                 _tileSize = layerTileSize;
00053 
00054             // we also need the maximum available data level.
00055             unsigned int layerMaxDataLevel = i->get()->getMaxDataLevel();
00056             if ( layerMaxDataLevel > _maxDataLevel )
00057                 _maxDataLevel = layerMaxDataLevel;
00058         }
00059     }
00060 }
00061 
00062 ElevationQuery::Technique
00063 ElevationQuery::getTechnique() const
00064 {
00065     return _technique;
00066 }
00067 
00068 void
00069 ElevationQuery::setTechnique( ElevationQuery::Technique technique )
00070 {
00071     _technique = technique;
00072 }
00073 
00074 void
00075 ElevationQuery::setMaxTilesToCache( int value )
00076 {
00077     _tileCache.setMaxSize( value );
00078 }
00079 
00080 int
00081 ElevationQuery::getMaxTilesToCache() const
00082 {
00083     return _tileCache.getMaxSize();
00084 }
00085 
00086 void
00087 ElevationQuery::setInterpolation( ElevationInterpolation interp)
00088 {
00089     _interpolation = interp;
00090 }
00091 
00092 ElevationInterpolation
00093 ElevationQuery::getElevationInterpolation() const
00094 {
00095     return _interpolation;
00096 }
00097 
00098 bool
00099 ElevationQuery::getElevation(const osg::Vec3d&       point,
00100                              const SpatialReference* pointSRS,
00101                              double&                 out_elevation,
00102                              double                  desiredResolution,
00103                              double*                 out_actualResolution)
00104 {
00105     sync();
00106     return getElevationImpl( point, pointSRS, out_elevation, desiredResolution, out_actualResolution );
00107 }
00108 
00109 bool
00110 ElevationQuery::getElevations(std::vector<osg::Vec3d>& points,
00111                               const SpatialReference*  pointsSRS,
00112                               bool                     ignoreZ,
00113                               double                   desiredResolution )
00114 {
00115     sync();
00116     for( osg::Vec3dArray::iterator i = points.begin(); i != points.end(); ++i )
00117     {
00118         double elevation;
00119         double z = (*i).z();
00120         if ( getElevationImpl( *i, pointsSRS, elevation, desiredResolution ) )
00121         {
00122             (*i).z() = ignoreZ ? elevation : elevation + z;
00123         }
00124     }
00125     return true;
00126 }
00127 
00128 bool
00129 ElevationQuery::getElevations(const std::vector<osg::Vec3d>& points,
00130                               const SpatialReference*        pointsSRS,
00131                               std::vector<double>&           out_elevations,
00132                               double                         desiredResolution )
00133 {
00134     sync();
00135     for( osg::Vec3dArray::const_iterator i = points.begin(); i != points.end(); ++i )
00136     {
00137         double elevation;
00138         if ( getElevationImpl( *i, pointsSRS, elevation, desiredResolution ) )
00139         {
00140             out_elevations.push_back( elevation );
00141         }
00142     }
00143     return true;
00144 }
00145 
00146 bool
00147 ElevationQuery::getElevationImpl(const osg::Vec3d&       point,
00148                                  const SpatialReference* pointSRS,
00149                                  double&                 out_elevation,
00150                                  double                  desiredResolution,
00151                                  double*                 out_actualResolution)
00152 {
00153     if ( _maxDataLevel == 0 || _tileSize == 0 )
00154     {
00155         // this means there are no heightfields.
00156         out_elevation = 0.0;
00157         return true;
00158     }
00159    
00160     // this is the ideal LOD for the requested resolution:
00161     unsigned int idealLevel = desiredResolution > 0.0
00162         ? _mapf.getProfile()->getLevelOfDetailForHorizResolution( desiredResolution, _tileSize )
00163         : _maxDataLevel;        
00164 
00165     // based on the heightfields available, this is the best we can theorically do:
00166     unsigned int bestAvailLevel = osg::minimum( idealLevel, _maxDataLevel );
00167     if (_maxLevelOverride >= 0)
00168     {
00169         bestAvailLevel = osg::minimum(bestAvailLevel, (unsigned int)_maxLevelOverride);
00170     }
00171     
00172     // transform the input coords to map coords:
00173     osg::Vec3d mapPoint = point;
00174     if ( pointSRS && !pointSRS->isEquivalentTo( _mapf.getProfile()->getSRS() ) )
00175     {
00176         if ( !pointSRS->transform2D( point.x(), point.y(), _mapf.getProfile()->getSRS(), mapPoint.x(), mapPoint.y() ) )
00177         {
00178             OE_WARN << LC << "Fail: coord transform failed" << std::endl;
00179             return false;
00180         }
00181     }
00182 
00183     osg::ref_ptr<osg::HeightField> hf;
00184     osg::ref_ptr<osgTerrain::TerrainTile> tile;
00185 
00186     // get the tilekey corresponding to the tile we need:
00187     TileKey key = _mapf.getProfile()->createTileKey( mapPoint.x(), mapPoint.y(), bestAvailLevel );
00188     if ( !key.valid() )
00189     {
00190         OE_WARN << LC << "Fail: coords fall outside map" << std::endl;
00191         return false;
00192     }
00193 
00194     // Check the tile cache. Note that the TileSource already likely has a MemCache
00195     // attached to it. We employ a secondary cache here for a couple reasons. One, this
00196     // cache will store not only the heightfield, but also the tesselated tile in the event
00197     // that we're using GEOMETRIC mode. Second, since the call the getHeightField can 
00198     // fallback on a lower resolution, this cache will hold the final resolution heightfield
00199     // instead of trying to fetch the higher resolution one each tiem.
00200 
00201     TileCache::Record record = _tileCache.get( key );
00202     if ( record.valid() )
00203         tile = record.value().get();
00204          
00205     // if we found it, make sure it has a heightfield in it:
00206     if ( tile.valid() )
00207     {
00208         osgTerrain::HeightFieldLayer* layer = dynamic_cast<osgTerrain::HeightFieldLayer*>(tile->getElevationLayer());
00209         if ( layer )
00210             hf = layer->getHeightField();
00211 
00212         if ( !hf.valid() )
00213             tile = 0L;
00214     }
00215 
00216     // if we didn't find it (or it didn't have heightfield data), build it.
00217     if ( !tile.valid() )
00218     {
00219         // generate the heightfield corresponding to the tile key, automatically falling back
00220         // on lower resolution if necessary:
00221         _mapf.getHeightField( key, true, hf, 0L, _interpolation );
00222 
00223         // bail out if we could not make a heightfield a all.
00224         if ( !hf.valid() )
00225         {
00226             OE_WARN << LC << "Unable to create heightfield for key " << key.str() << std::endl;
00227             return false;
00228         }
00229 
00230         // All this stuff is requires for GEOMETRIC mode. An optimization would be to
00231         // defer this so that PARAMETRIC mode doesn't waste time
00232         GeoLocator* locator = GeoLocator::createForKey( key, _mapf.getMapInfo() );
00233 
00234         tile = new osgTerrain::TerrainTile();
00235 
00236         osgTerrain::HeightFieldLayer* layer = new osgTerrain::HeightFieldLayer( hf.get() );
00237         layer->setLocator( locator );
00238 
00239         tile->setElevationLayer( layer );
00240         tile->setRequiresNormals( false );
00241         tile->setTerrainTechnique( new osgTerrain::GeometryTechnique );
00242 
00243         // store it in the local tile cache.
00244         _tileCache.insert( key, tile.get() );
00245     }
00246 
00247     OE_DEBUG << LC << "LRU Cache, hit ratio = " << _tileCache.getStats()._hitRatio << std::endl;
00248 
00249     // see what the actual resolution of the heightfield is.
00250     if ( out_actualResolution )
00251         *out_actualResolution = (double)hf->getXInterval();
00252 
00253     // finally it's time to get a height value:
00254     if ( _technique == TECHNIQUE_PARAMETRIC )
00255     {
00256         const GeoExtent& extent = key.getExtent();
00257         double xInterval = extent.width()  / (double)(hf->getNumColumns()-1);
00258         double yInterval = extent.height() / (double)(hf->getNumRows()-1);
00259         out_elevation = (double) HeightFieldUtils::getHeightAtLocation( 
00260             hf.get(), mapPoint.x(), mapPoint.y(), extent.xMin(), extent.yMin(), xInterval, yInterval );
00261         return true;
00262     }
00263     else // ( _technique == TECHNIQUE_GEOMETRIC )
00264     {
00265         osg::Vec3d start, end, zero;
00266 
00267         if ( _mapf.getMapInfo().isGeocentric() )
00268         {
00269             const SpatialReference* mapSRS = _mapf.getProfile()->getSRS();
00270 
00271             mapSRS->transformToECEF( osg::Vec3d(mapPoint.y(), mapPoint.x(),  50000.0), start );
00272             mapSRS->transformToECEF( osg::Vec3d(mapPoint.y(), mapPoint.x(), -50000.0), end );
00273             mapSRS->transformToECEF( osg::Vec3d(mapPoint.y(), mapPoint.x(),      0.0), zero );
00274         }
00275         else // PROJECTED
00276         {
00277             start.set( mapPoint.x(), mapPoint.y(),  50000.0 );
00278             end.set  ( mapPoint.x(), mapPoint.y(), -50000.0 );
00279             zero.set ( mapPoint.x(), mapPoint.y(),      0.0 );
00280         }
00281 
00282         osgUtil::LineSegmentIntersector* i = new osgUtil::LineSegmentIntersector( start, end );
00283         osgUtil::IntersectionVisitor iv;
00284         iv.setIntersector( i );
00285 
00286         tile->accept( iv );
00287 
00288         osgUtil::LineSegmentIntersector::Intersections& results = i->getIntersections();
00289         if ( !results.empty() )
00290         {
00291             const osgUtil::LineSegmentIntersector::Intersection& result = *results.begin();
00292             osg::Vec3d isectPoint = result.getWorldIntersectPoint();
00293             out_elevation = (isectPoint-end).length2() > (zero-end).length2()
00294                 ? (isectPoint-zero).length()
00295                 : -(isectPoint-zero).length();
00296             return true;            
00297         }
00298 
00299         OE_DEBUG << LC << "No intersection" << std::endl;
00300         return false;
00301     }
00302 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines