C++ GeoExtent::width方法代码示例

osg::BoundingSphered
FeatureModelGraph::getBoundInWorldCoords(const GeoExtent& extent,
                                         const MapFrame*  mapf ) const
{
    osg::Vec3d center, corner;
    GeoExtent workingExtent;

    if ( !extent.isValid() )
    {
        return osg::BoundingSphered();
    }

    if ( extent.getSRS()->isEquivalentTo( _usableMapExtent.getSRS() ) )
    {
        workingExtent = extent;
    }
    else
    {
        workingExtent = extent.transform( _usableMapExtent.getSRS() ); // safe.
    }

    workingExtent.getCentroid( center.x(), center.y() );
    
    double centerZ = 0.0;    
    if ( mapf )
    {
        // Use an appropriate resolution for this extents width
        double resolution = workingExtent.width();
        ElevationQuery query( *mapf );
        GeoPoint p( mapf->getProfile()->getSRS(), center, ALTMODE_ABSOLUTE );
        query.getElevation( p, center.z(), resolution );
        centerZ = center.z();
    }    

    corner.x() = workingExtent.xMin();
    corner.y() = workingExtent.yMin();
    corner.z() = 0;

    if ( _session->getMapInfo().isGeocentric() )
    {
        const SpatialReference* ecefSRS = workingExtent.getSRS()->getECEF();
        workingExtent.getSRS()->transform( center, ecefSRS, center );
        workingExtent.getSRS()->transform( corner, ecefSRS, corner );
        //workingExtent.getSRS()->transformToECEF( center, center );
        //workingExtent.getSRS()->transformToECEF( corner, corner );
    }

    if (workingExtent.getSRS()->isGeographic() &&
        ( workingExtent.width() >= 90 || workingExtent.height() >= 90 ) )
    {
        return osg::BoundingSphered( osg::Vec3d(0,0,0), 2*center.length() );
    }

    return osg::BoundingSphered( center, (center-corner).length() );
}

osg::BoundingSphere SimplePager::getBounds(const TileKey& key) const
{
    int samples = 6;

    GeoExtent extent = key.getExtent();

    double xSample = extent.width() / (double)samples;
    double ySample = extent.height() / (double)samples;

    osg::BoundingSphere bs;
    for (int c = 0; c < samples+1; c++)
    {
        double x = extent.xMin() + (double)c * xSample;
        for (int r = 0; r < samples+1; r++)
        {
            double y = extent.yMin() + (double)r * ySample;
            osg::Vec3d world;

            GeoPoint samplePoint(extent.getSRS(), x, y, 0, ALTMODE_ABSOLUTE);

            GeoPoint wgs84 = samplePoint.transform(osgEarth::SpatialReference::create("epsg:4326"));
            wgs84.toWorld(world);
            bs.expandBy(world);
        }
    }
    return bs;
}

void
FeatureNode::updateClusterCulling()
{
    // install a cluster culler.
    if ( getMapNode()->isGeocentric() && _clusterCulling && !_clusterCullingCallback)
    {
        const GeoExtent& ccExtent = _extent;
        if ( ccExtent.isValid() )
        {
            // if the extent is more than 90 degrees, bail
            GeoExtent geodeticExtent = ccExtent.transform( ccExtent.getSRS()->getGeographicSRS() );
            if ( geodeticExtent.width() < 90.0 && geodeticExtent.height() < 90.0 )
            {
                // get the geocentric tile center:
                osg::Vec3d tileCenter;
                ccExtent.getCentroid( tileCenter.x(), tileCenter.y() );

                osg::Vec3d centerECEF;
                ccExtent.getSRS()->transform( tileCenter, getMapNode()->getMapSRS()->getECEF(), centerECEF );
                _clusterCullingCallback = ClusterCullingFactory::create2( this, centerECEF );
                if ( _clusterCullingCallback )
                    this->addCullCallback( _clusterCullingCallback );
            }
        }
    }
    else if (!_clusterCulling && _clusterCullingCallback)
    {
        this->removeCullCallback( _clusterCullingCallback );
        _clusterCullingCallback = 0;
    }
}

GeoLocator::GeoLocator( const osgTerrain::Locator& prototype, const GeoExtent& dataExtent, const GeoExtent& displayExtent ) :
osgTerrain::Locator( prototype ),
_dataExtent( dataExtent ),
_inverseCalculated(false)
{
    // assume they are the same SRS
    _x0 = osg::clampBetween( (displayExtent.xMin()-dataExtent.xMin())/dataExtent.width(), 0.0, 1.0 );
    _x1 = osg::clampBetween( (displayExtent.xMax()-dataExtent.xMin())/dataExtent.width(), 0.0, 1.0 );
    _y0 = osg::clampBetween( (displayExtent.yMin()-dataExtent.yMin())/dataExtent.height(), 0.0, 1.0 );
    _y1 = osg::clampBetween( (displayExtent.yMax()-dataExtent.yMin())/dataExtent.height(), 0.0, 1.0 );
}

bool
GeoLocator::createScaleBiasMatrix(const GeoExtent& window, osg::Matrixf& out) const
{
    float scalex = window.width() / _dataExtent.width();
    float scaley = window.height() / _dataExtent.height();
    float biasx  = (window.xMin()-_dataExtent.xMin()) / _dataExtent.width();
    float biasy  = (window.yMin()-_dataExtent.yMin()) / _dataExtent.height();

    out(0,0) = scalex;
    out(1,1) = scaley;
    out(3,0) = biasx;
    out(3,1) = biasy;

    return true;
}

osg::HeightField*
HeightFieldUtils::createReferenceHeightField(const GeoExtent& ex,
                                             unsigned         numCols,
                                             unsigned         numRows,
                                             bool             expressAsHAE)
{
    osg::HeightField* hf = new osg::HeightField();
    hf->allocate( numCols, numRows );
    hf->setOrigin( osg::Vec3d( ex.xMin(), ex.yMin(), 0.0 ) );
    hf->setXInterval( (ex.xMax() - ex.xMin())/(double)(numCols-1) );
    hf->setYInterval( (ex.yMax() - ex.yMin())/(double)(numRows-1) );

    const VerticalDatum* vdatum = ex.isValid() ? ex.getSRS()->getVerticalDatum() : 0L;

    if ( vdatum && expressAsHAE )
    {
        // need the lat/long extent for geoid queries:
        GeoExtent geodeticExtent = ex.getSRS()->isGeographic() ? ex : ex.transform( ex.getSRS()->getGeographicSRS() );
        double latMin = geodeticExtent.yMin();
        double lonMin = geodeticExtent.xMin();
        double lonInterval = geodeticExtent.width() / (double)(numCols-1);
        double latInterval = geodeticExtent.height() / (double)(numRows-1);

        for( unsigned r=0; r<numRows; ++r )
        {            
            double lat = latMin + latInterval*(double)r;
            for( unsigned c=0; c<numCols; ++c )
            {
                double lon = lonMin + lonInterval*(double)c;
                double offset = vdatum->msl2hae(lat, lon, 0.0);
                hf->setHeight( c, r, offset );
            }
        }
    }
    else
    {
        for(unsigned int i=0; i<hf->getHeightList().size(); i++ )
        {
            hf->getHeightList()[i] = 0.0;
        }
    }

    hf->setBorderWidth( 0 );
    return hf;    
}

void
HeightFieldUtils::resolveInvalidHeights(osg::HeightField* grid,
                                        const GeoExtent&  ex,
                                        float             invalidValue,
                                        const Geoid*      geoid)
{
    if ( geoid )
    {
        // need the lat/long extent for geoid queries:
        unsigned numRows = grid->getNumRows();
        unsigned numCols = grid->getNumColumns();
        GeoExtent geodeticExtent = ex.getSRS()->isGeographic() ? ex : ex.transform( ex.getSRS()->getGeographicSRS() );
        double latMin = geodeticExtent.yMin();
        double lonMin = geodeticExtent.xMin();
        double lonInterval = geodeticExtent.width() / (double)(numCols-1);
        double latInterval = geodeticExtent.height() / (double)(numRows-1);

        for( unsigned r=0; r<numRows; ++r )
        {
            double lat = latMin + latInterval*(double)r;
            for( unsigned c=0; c<numCols; ++c )
            {
                double lon = lonMin + lonInterval*(double)c;
                if ( grid->getHeight(c, r) == invalidValue )
                {
                    grid->setHeight( c, r, geoid->getHeight(lat, lon) );
                }
            }
        }
    }
    else
    {
        for(unsigned int i=0; i<grid->getHeightList().size(); i++ )
        {
            if ( grid->getHeightList()[i] == invalidValue )
            {
                grid->getHeightList()[i] = 0.0;
            }
        }
    }
}

GeoHeightField
GeoHeightField::createSubSample( const GeoExtent& destEx, ElevationInterpolation interpolation) const
{
    double div = destEx.width()/_extent.width();
    if ( div >= 1.0f )
        return GeoHeightField::INVALID;

    int w = _heightField->getNumColumns();
    int h = _heightField->getNumRows();
    //double dx = _heightField->getXInterval() * div;
    //double dy = _heightField->getYInterval() * div;
    double xInterval = _extent.width() / (double)(_heightField->getNumColumns()-1);
    double yInterval = _extent.height() / (double)(_heightField->getNumRows()-1);
    double dx = xInterval * div;
    double dy = yInterval * div;

    osg::HeightField* dest = new osg::HeightField();
    dest->allocate( w, h );
    dest->setXInterval( dx );
    dest->setYInterval( dy );

    // copy over the skirt height, adjusting it for tile size.
    dest->setSkirtHeight( _heightField->getSkirtHeight() * div );

    double x, y;
    int col, row;

    for( x = destEx.xMin(), col=0; col < w; x += dx, col++ )
    {
        for( y = destEx.yMin(), row=0; row < h; y += dy, row++ )
        {
            float height = HeightFieldUtils::getHeightAtLocation( _heightField.get(), x, y, _extent.xMin(), _extent.yMin(), xInterval, yInterval, interpolation);
            dest->setHeight( col, row, height );
        }
    }

    osg::Vec3d orig( destEx.xMin(), destEx.yMin(), _heightField->getOrigin().z() );
    dest->setOrigin( orig );

    return GeoHeightField( dest, destEx, _vsrs.get() );
}

osg::BoundingSphered
FeatureModelGraph::getBoundInWorldCoords(const GeoExtent& extent,
                                         const MapFrame*  mapf ) const
{
    osg::Vec3d center, corner;
    //double z = 0.0;
    GeoExtent workingExtent;

    if ( extent.getSRS()->isEquivalentTo( _usableMapExtent.getSRS() ) )
    {
        workingExtent = extent;
    }
    else
    {
        workingExtent = extent.transform( _usableMapExtent.getSRS() ); // safe.
    }

    workingExtent.getCentroid( center.x(), center.y() );
    
    double centerZ = 0.0;    
    if ( mapf )
    {
        // Use an appropriate resolution for this extents width
        double resolution = workingExtent.width();             
        ElevationQuery query( *mapf );
        query.getElevation( GeoPoint(mapf->getProfile()->getSRS(),center), center.z(), resolution );
        centerZ = center.z();
    }    

    corner.x() = workingExtent.xMin();
    corner.y() = workingExtent.yMin();
    corner.z() = 0;

    if ( _session->getMapInfo().isGeocentric() )
    {
        workingExtent.getSRS()->transformToECEF( center, center );
        workingExtent.getSRS()->transformToECEF( corner, corner );
    }

    return osg::BoundingSphered( center, (center-corner).length() );
}

void
FeaturesToNodeFilter::computeLocalizers( const FilterContext& context )
{
    if ( context.isGeoreferenced() )
    {
        if ( context.getSession()->getMapInfo().isGeocentric() )
        {
            const SpatialReference* geogSRS = context.profile()->getSRS()->getGeographicSRS();
            GeoExtent geodExtent = context.extent()->transform( geogSRS );
            if ( geodExtent.width() < 180.0 )
            {
                osg::Vec3d centroid, centroidECEF;
                geodExtent.getCentroid( centroid.x(), centroid.y() );
                geogSRS->transform( centroid, geogSRS->getECEF(), centroidECEF );
                geogSRS->getECEF()->createLocalToWorld( centroidECEF, _local2world );
                _world2local.invert( _local2world );
            }
        }

        else // projected
        {
            if ( context.extent().isSet() )
            {
                osg::Vec3d centroid;
                context.extent()->getCentroid(centroid.x(), centroid.y());

                context.extent()->getSRS()->transform(
                    centroid,
                    context.getSession()->getMapInfo().getProfile()->getSRS(),
                    centroid );

                _world2local.makeTranslate( -centroid );
                _local2world.invert( _world2local );
            }
        }
    }
}

GeoGraph::GeoGraph(const GeoExtent& extent, float maxRange, unsigned maxObjects,
                   unsigned splitDim, float splitRangeFactor,
                   unsigned rootWidth, unsigned rootHeight ) :
GeoCell( extent, maxRange, maxObjects, splitDim, splitRangeFactor, 0 )
{
    _rootWidth = osg::maximum( rootWidth, (unsigned)2 );
    _rootHeight = osg::maximum( rootHeight, (unsigned)2 );

    if ( _depth == 0 )
    {
        double xinterval = extent.width() / (double)_rootWidth;
        double yinterval = extent.height() / (double)_rootHeight;

        for( unsigned y=0; y<_rootHeight; ++y )
        {
            for( unsigned x=0; x<_rootWidth; ++x )
            {
                GeoExtent cellExtent(
                    _extent.getSRS(),
                    _extent.xMin() + xinterval*(double)x,
                    _extent.yMin() + yinterval*(double)y,
                    _extent.xMin() + xinterval*(double)(x+1),
                    _extent.yMin() + yinterval*(double)(y+1) );

                GeoCell * child = new GeoCell(
                    cellExtent,
                    _maxRange,
                    _maxObjects,
                    _splitDim,
                    _splitRangeFactor,
                    1 );

                this->addChild( child, 0, maxRange ); //FLT_MAX );
            }
        }                    
    }
}

void
FeaturesToNodeFilter::computeLocalizers( const FilterContext& context, const osgEarth::GeoExtent &extent, osg::Matrixd &out_w2l, osg::Matrixd &out_l2w )
{
    if ( context.isGeoreferenced() )
    {
        if ( context.getSession()->getMapInfo().isGeocentric() )
        {
            const SpatialReference* geogSRS = context.profile()->getSRS()->getGeographicSRS();
            GeoExtent geodExtent = extent.transform( geogSRS );
            if ( geodExtent.width() < 180.0 )
            {
                osg::Vec3d centroid, centroidECEF;
                geodExtent.getCentroid( centroid.x(), centroid.y() );
                geogSRS->transform( centroid, geogSRS->getECEF(), centroidECEF );
                geogSRS->getECEF()->createLocalToWorld( centroidECEF, out_l2w );
                out_w2l.invert( out_l2w );
            }
        }

        else // projected
        {
            if ( extent.isValid() )
            {
                osg::Vec3d centroid;
                extent.getCentroid(centroid.x(), centroid.y());

                extent.getSRS()->transform(
                    centroid,
                    context.getSession()->getMapInfo().getProfile()->getSRS(),
                    centroid );

                out_w2l.makeTranslate( -centroid );
                out_l2w.invert( out_w2l );
            }
        }
    }
}

void
Profile::addIntersectingTiles(const GeoExtent& key_ext, unsigned localLOD, std::vector<TileKey>& out_intersectingKeys) const
{
    // assume a non-crossing extent here.
    if ( key_ext.crossesAntimeridian() )
    {
        OE_WARN << "Profile::addIntersectingTiles cannot process date-line cross" << std::endl;
        return;
    }

    int tileMinX, tileMaxX;
    int tileMinY, tileMaxY;

    // Special path for mercator (does NOT work for cube, e.g.)
    if ( key_ext.getSRS()->isMercator() )
    {
        int precision = 5;
        double eps = 0.001;

        double keyWidth = round(key_ext.width(), precision);
        int destLOD = 0;
        double w, h;
        getTileDimensions(0, w, h);
        for(; (round(w,precision) - keyWidth) > eps; w*=0.5, h*=0.5, destLOD++ );

        double destTileWidth, destTileHeight;
        getTileDimensions( destLOD, destTileWidth, destTileHeight );
        destTileWidth = round(destTileWidth, precision);
        destTileHeight = round(destTileHeight, precision);

        tileMinX = quantize( ((key_ext.xMin() - _extent.xMin()) / destTileWidth), eps );
        tileMaxX = (int)((key_ext.xMax() - _extent.xMin()) / destTileWidth);

        tileMinY = quantize( ((_extent.yMax() - key_ext.yMax()) / destTileHeight), eps );
        tileMaxY = (int) ((_extent.yMax() - key_ext.yMin()) / destTileHeight);
    }

    else
    {
        double destTileWidth, destTileHeight;
        getTileDimensions(localLOD, destTileWidth, destTileHeight);

        //OE_DEBUG << std::fixed << "  Source Tile: " << key.getLevelOfDetail() << " (" << keyWidth << ", " << keyHeight << ")" << std::endl;
        //OE_DEBUG << std::fixed << "  Dest Size: " << destLOD << " (" << destTileWidth << ", " << destTileHeight << ")" << std::endl;

        tileMinX = (int)((key_ext.xMin() - _extent.xMin()) / destTileWidth);
        tileMaxX = (int)((key_ext.xMax() - _extent.xMin()) / destTileWidth);

        tileMinY = (int)((_extent.yMax() - key_ext.yMax()) / destTileHeight);
        tileMaxY = (int)((_extent.yMax() - key_ext.yMin()) / destTileHeight);
    }

    unsigned int numWide, numHigh;
    getNumTiles(localLOD, numWide, numHigh);

    // bail out if the tiles are out of bounds.
    if ( tileMinX >= (int)numWide || tileMinY >= (int)numHigh ||
            tileMaxX < 0 || tileMaxY < 0 )
    {
        return;
    }

    tileMinX = osg::clampBetween(tileMinX, 0, (int)numWide-1);
    tileMaxX = osg::clampBetween(tileMaxX, 0, (int)numWide-1);
    tileMinY = osg::clampBetween(tileMinY, 0, (int)numHigh-1);
    tileMaxY = osg::clampBetween(tileMaxY, 0, (int)numHigh-1);

    OE_DEBUG << std::fixed << "  Dest Tiles: " << tileMinX << "," << tileMinY << " => " << tileMaxX << "," << tileMaxY << std::endl;

    for (int i = tileMinX; i <= tileMaxX; ++i)
    {
        for (int j = tileMinY; j <= tileMaxY; ++j)
        {
            //TODO: does not support multi-face destination keys.
            out_intersectingKeys.push_back( TileKey(localLOD, i, j, this) );
        }
    }
}

void
    TFSPackager::package( FeatureSource* features, const std::string& destination, const std::string& layername, const std::string& description )
{   
    if (!_destSRSString.empty())
    {
        _srs = SpatialReference::create( _destSRSString );
    }

    //Get the destination SRS from the feature source if it's not already set
    if (!_srs.valid())
    {
        _srs = features->getFeatureProfile()->getSRS();
    }
	
    //Get the extent of the dataset, or use the custom extent value
    GeoExtent srsExtent = _customExtent;
    if (!srsExtent.isValid())
        srsExtent = features->getFeatureProfile()->getExtent();

    //Transform to lat/lon extents
    GeoExtent extent = srsExtent.transform( _srs.get() );

    osg::ref_ptr< const osgEarth::Profile > profile = osgEarth::Profile::create(extent.getSRS(), extent.xMin(), extent.yMin(), extent.xMax(), extent.yMax(), 1, 1);


    TileKey rootKey = TileKey(0, 0, 0, profile );    


    osg::ref_ptr< FeatureTile > root = new FeatureTile( rootKey );
    //Loop through all the features and try to insert them into the quadtree
    osg::ref_ptr< FeatureCursor > cursor = features->createFeatureCursor( _query );
    int added = 0;
    int failed = 0;
    int skipped = 0;
    int highestLevel = 0;

    while (cursor.valid() && cursor->hasMore())
    {        
        osg::ref_ptr< Feature > feature = cursor->nextFeature();

        //Reproject the feature to the dest SRS if it's not already
        if (!feature->getSRS()->isEquivalentTo( _srs ) )
        {
            feature->transform( _srs );
        }

        if (feature->getGeometry() && feature->getGeometry()->getBounds().valid() && feature->getGeometry()->isValid())
        {

            AddFeatureVisitor v(feature.get(), _maxFeatures, _firstLevel, _maxLevel, _method);
            root->accept( &v );
            if (!v._added)
            {
                OE_NOTICE << "Failed to add feature " << feature->getFID() << std::endl;
                failed++;
            }
            else
            {                
                if (highestLevel < v._levelAdded)
                {
                    highestLevel = v._levelAdded;
                }
                added++;
                OE_DEBUG << "Added " << added << std::endl;
            }   
        }
        else
        {
            OE_NOTICE << "Skipping feature " << feature->getFID() << " with null or invalid geometry" << std::endl;
            skipped++;
        }
    }   
    OE_NOTICE << "Added=" << added << " Skipped=" << skipped << " Failed=" << failed << std::endl;

#if 1
    // Print the width of tiles at each level
    for (int i = 0; i <= highestLevel; ++i)
    {
        TileKey tileKey(i, 0, 0, profile);
        GeoExtent tileExtent = tileKey.getExtent();
        OE_NOTICE << "Level " << i << " tile size: " << tileExtent.width() << std::endl;
    }
#endif

    WriteFeaturesVisitor write(features, destination, _method, _srs);
    root->accept( &write );

    //Write out the meta doc
    TFSLayer layer;
    layer.setTitle( layername );
    layer.setAbstract( description );
    layer.setFirstLevel( _firstLevel );
    layer.setMaxLevel( highestLevel );
    layer.setExtent( profile->getExtent() );
    layer.setSRS( _srs.get() );
    TFSReaderWriter::write( layer, osgDB::concatPaths( destination, "tfs.xml"));

}

/**
 * Prepares a geometry into a grid if it is too big geospatially to have a sensible local tangent plane
 * We will also tile the geometry if it just has too many points to speed up the tesselator.
 */
void prepareForTesselation(Geometry* geometry, const SpatialReference* featureSRS, double targetTileSizeDeg, unsigned int maxPointsPerTile, GeometryCollection& out)
{
    // Clear the output list.
    GeometryCollection tiles;
    
    unsigned int count = geometry->size();

    unsigned int tx = 1;
    unsigned int ty = 1;

    // Tile the geometry if it's geospatial size is too large to have a sensible local tangent plane.
    GeoExtent featureExtentDeg = GeoExtent(featureSRS, geometry->getBounds()).transform(SpatialReference::create("wgs84"));

    // Tile based on the extent
    if ( featureExtentDeg.width() > targetTileSizeDeg  || featureExtentDeg.height() > targetTileSizeDeg)
    {
        // Determine the tile size based on the extent.
        tx = ceil( featureExtentDeg.width() / targetTileSizeDeg );
        ty = ceil (featureExtentDeg.height() / targetTileSizeDeg );        
    }
    else if (count > maxPointsPerTile)
    {
        // Determine the size based on the number of points.
        unsigned numTiles = ((double)count / (double)maxPointsPerTile) + 1u;
        tx = ceil(sqrt((double)numTiles));
        ty = tx;        
    }

    if (tx == 1 && ty == 1)
    {
        // The geometry doesn't need modified so just add it to the list.
        tiles.push_back( geometry );
    }
    else
    {
        tileGeometry( geometry, featureSRS, tx, ty, tiles );
    }

    out.clear();

#if 1
    // Just copy the output tiles to the output.
    std::copy(tiles.begin(), tiles.end(), std::back_inserter(out));
#else
    // Calling this code will recursively subdivide the cells based on the number of points they have.
    // This works but it will produces a non-regular grid which doesn't render well in geocentric
    // due to the curvature of the earth so we disable it for now.
    //
    // Reduce the size of the tiles if needed.
    for (unsigned int i = 0; i < tiles.size(); i++)
    {
        if (tiles[i]->size() > maxPointsPerTile)
        {
            GeometryCollection tmp;
            downsizeGeometry(tiles[i].get(), featureSRS, maxPointsPerTile, tmp);
            std::copy(tmp.begin(), tmp.end(), std::back_inserter(out));
        }
        else
        {
            out.push_back( tiles[i].get() );
        }
    }
#endif
}