osgEarth 2.1.1
|
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 }