C++ GeoPoint类代码示例

void
PolyhedralLineOfSightNode::updateSamples()
{
    if ( _geode->getNumDrawables() == 0 )
        rebuildGeometry();

    osg::Geometry* geom  = _geode->getDrawable(0)->asGeometry();
    osg::Vec3Array* verts = dynamic_cast<osg::Vec3Array*>( geom->getVertexArray() );

    double distance = _distance.as(Units::METERS);

    // get the world coords and a l2w transform for the origin:
    osg::Vec3d  originWorld;
    osg::Matrix local2world, world2local;

    Terrain* t = getMapNode()->getTerrain();
    GeoPoint origin = getPosition();
    origin.makeAbsolute( t );
    origin.toWorld( originWorld, t );
    origin.createLocalToWorld( local2world );
    world2local.invert( local2world );

    // set up an intersector:
    osgUtil::LineSegmentIntersector* lsi = new osgUtil::LineSegmentIntersector( originWorld, originWorld );
    osgUtil::IntersectionVisitor iv( lsi );

    // intersect the verts (skip the origin point) with the map node.
    for( osg::Vec3Array::iterator v = verts->begin()+1; v != verts->end(); ++v )
    {
        osg::Vec3d unit = *v;
        unit.normalize();
        unit *= distance;

        osg::Vec3d world = unit * local2world;

        if ( osg::equivalent(unit.length(), 0.0) )
        {
            OE_WARN << "problem." << std::endl;
        }

        lsi->reset();
        lsi->setStart( originWorld );
        lsi->setEnd( world );

        OE_DEBUG << LC << "Ray: " <<
            originWorld.x() << "," << originWorld.y() << "," << originWorld.z() << " => "
            << world.x() << "," << world.y() << "," << world.z()
            << std::endl;


        getMapNode()->getTerrain()->accept( iv );

        osgUtil::LineSegmentIntersector::Intersections& hits = lsi->getIntersections();
        if ( !hits.empty() )
        {
            const osgUtil::LineSegmentIntersector::Intersection& hit = *hits.begin();
            osg::Vec3d newV = hit.getWorldIntersectPoint() * world2local;
            if ( newV.length() > 1.0 )
                *v = newV;
            else
                *v = unit;
        }
        else
        {
            *v = unit;
        }

        //OE_NOTICE << "Radial point = " << v->x() << ", " << v->y() << ", " << v->z() << std::endl;
    }

    verts->dirty();
    geom->dirtyBound();
}

void
AutoClipPlaneCullCallback::operator()( osg::Node* node, osg::NodeVisitor* nv )
{
    if ( _active )
    {
        osgUtil::CullVisitor* cv = Culling::asCullVisitor(nv);
        if ( cv )
        {
            osgEarth::Map* map = _mapNode.valid() ? _mapNode->getMap() : 0;
            osg::Camera* cam = cv->getCurrentCamera();
            osg::ref_ptr<osg::CullSettings::ClampProjectionMatrixCallback>& clamper = _clampers.get(cam);
            if ( !clamper.valid() )
            {
                clamper = new CustomProjClamper();
                cam->setClampProjectionMatrixCallback( clamper.get() );
                OE_INFO << LC << "Installed custom projeciton matrix clamper" << std::endl;
            }
            else
            {
                CustomProjClamper* c = static_cast<CustomProjClamper*>(clamper.get());

                osg::Vec3d eye, center, up;
                cam->getViewMatrixAsLookAt( eye, center, up );

                // clamp the far clipping plane to the approximate horizon distance
                if ( _autoFarPlaneClamping )
                {
                    double d = eye.length();
                    c->_maxFar = sqrt( d*d - _rp2 );
                }
                else
                {
                    c->_maxFar = DBL_MAX;
                }

                // get the height-above-ellipsoid. If we need to be more accurate, we can use 
                // ElevationQuery in the future..
                //osg::Vec3d loc;
                GeoPoint loc;
                if ( map )
                {
                    loc.fromWorld( map->getSRS(), eye );
                    //map->worldPointToMapPoint( eye, loc );
                }
                else
                {
                    static osg::EllipsoidModel em;
                    osg::Vec3d t;
                    em.convertXYZToLatLongHeight( eye.x(), eye.y(), eye.z(), loc.y(), loc.x(), loc.z() );
                }
                
                //double hae = loc.z();
                double hae = loc.z();
                if (_mapNode.valid())
                {
                    double height = 0.0;
                    _mapNode->getTerrain()->getHeight(loc.getSRS(), loc.x(), loc.y(), &height);
                    //OE_NOTICE << "got height " << height << std::endl;
                    hae -= height;
                    //OE_NOTICE << "HAE=" << hae <<  std::endl;
                }

                // ramp a new near/far ratio based on the HAE.
                c->_nearFarRatio = Utils::remap( hae, 0.0, _haeThreshold, _minNearFarRatio, _maxNearFarRatio );
            }

#if 0
            {
                double n, f, a, v;
                cv->getProjectionMatrix()->getPerspective(v, a, n, f);
                OE_INFO << std::setprecision(16) << "near = " << n << ", far = " << f << ", ratio = " << n/f << std::endl;
            }
#endif
        }
    }
    traverse( node, nv );
}

void TerrainProfileGraph::drawHoverCursor(const QPointF& position)
{
  if (_hoverLine)
  {
    _scene->removeItem(_hoverLine);
    delete _hoverLine;
    _hoverLine = 0L;
  }

  if (_graphField.width() < 2 || _graphField.height() < 2)
    return;

  double xPos = position.x() < _graphField.x() ? _graphField.x() : (position.x() > _graphField.x() + _graphField.width() ? _graphField.x() + _graphField.width() : position.x());

  QLineF vLine(xPos, _graphField.y(), xPos, _graphField.y() + _graphField.height());

  QPointF* intersect = new QPointF;
  bool foundIntersect = false;
  for (int i=0; i < _graphLines.count(); i++)
  {
    if (vLine.intersect(_graphLines[i], intersect) == QLineF::BoundedIntersection)
    {
      foundIntersect = true;
      break;
    }
  }

  if (foundIntersect)
  {
    // Draw the upper line segment.  Also serves as the parent item.
    _hoverLine = new QGraphicsLineItem(xPos, _graphField.y(), xPos, intersect->y() - 3);
    _hoverLine->setPen(_hoverPen);
    _hoverLine->setZValue(OVERLAY_Z);
    _scene->addItem(_hoverLine);

    // Draw the box around the intersect point
    QGraphicsRectItem* hoverBox = new QGraphicsRectItem(xPos - 3, intersect->y() - 3, 6, 6);
    hoverBox->setPen(_hoverPen);
    hoverBox->setBrush(Qt::NoBrush);
    hoverBox->setZValue(OVERLAY_Z);
    hoverBox->setParentItem(_hoverLine);

    // Draw the lower line segment
    QGraphicsLineItem* lowerLine = new QGraphicsLineItem(xPos, intersect->y() + 3, xPos, _graphField.y() + _graphField.height() + 5);
    lowerLine->setPen(_hoverPen);
    lowerLine->setZValue(OVERLAY_Z);
    lowerLine->setParentItem(_hoverLine);

    // Draw the text and background
    double y = (1.0 - ((intersect->y() - _graphField.y()) / _graphField.height())) * (_graphMaxY - _graphMinY) + _graphMinY;
    int textOffset = 10;

    QGraphicsSimpleTextItem* hoverText = new QGraphicsSimpleTextItem(QString::number(y) + tr("m"));
    hoverText->setBrush(QBrush(_axesColor));
    hoverText->setFont(_graphFont);
    hoverText->setZValue(OVERLAY_Z);

    if (intersect->x() + textOffset + hoverText->boundingRect().width() < _graphField.x() + _graphField.width())
      hoverText->setPos(intersect->x() + textOffset, intersect->y() - hoverText->boundingRect().height());
    else
      hoverText->setPos(intersect->x() - textOffset - hoverText->boundingRect().width(), intersect->y() - hoverText->boundingRect().height());

    QGraphicsRectItem* hoverTextBackground = new QGraphicsRectItem(hoverText->x() - 3, hoverText->y() - 1, 
                                                                   hoverText->boundingRect().width() + 6,
                                                                   hoverText->boundingRect().height() + 1);
    hoverTextBackground->setPen(_axesPen);
    hoverTextBackground->setBrush(QBrush(_graphColor));
    hoverTextBackground->setZValue(OVERLAY_Z);
    hoverTextBackground->setParentItem(_hoverLine);

    hoverText->setParentItem(_hoverLine);

    // Update callback
    if (_positionCallback.valid())
    {
      double distanceFactor = ((xPos - _graphField.x()) / (double)_graphField.width());

      osg::Vec3d worldStart, worldEnd;
      _calculator->getStart(ALTMODE_ABSOLUTE).toWorld(worldStart);
      _calculator->getEnd(ALTMODE_ABSOLUTE).toWorld(worldEnd);

      double worldX = (worldEnd.x() - worldStart.x()) * distanceFactor + worldStart.x();
      double worldY = (worldEnd.y() - worldStart.y()) * distanceFactor + worldStart.y();
      double worldZ = (worldEnd.z() - worldStart.z()) * distanceFactor + worldStart.z();

      GeoPoint mapPos;
      mapPos.fromWorld(_calculator->getStart().getSRS(), osg::Vec3d(worldX, worldY, worldZ));

      _positionCallback->updatePosition(mapPos.y(), mapPos.x(), hoverText->text().toStdString());
    }
  }
  else
  {
    // No intersect found so just draw the full line at xPos
    _hoverLine = new QGraphicsLineItem(xPos, _graphField.y(), xPos, _graphField.y() + _graphField.height() + 5);
    _hoverLine->setPen(_hoverPen);
    _hoverLine->setZValue(OVERLAY_Z);
    _scene->addItem(_hoverLine);
  }

//.........这里部分代码省略.........

    int Ardb::GeoSearch(const DBID& db, const Slice& key, const GeoSearchOptions& options, ValueDataDeque& results)
    {
        uint64 start_time = get_current_epoch_micros();
        GeoHashBitsSet ress;
        double x = options.x, y = options.y;
        if (options.coord_type == GEO_WGS84_TYPE)
        {
            x = GeoHashHelper::GetMercatorX(options.x);
            y = GeoHashHelper::GetMercatorY(options.y);
        }
        if (options.by_member)
        {
            ValueData score, attr;
            int err = ZGetNodeValue(db, key, options.member, score, attr);
            if (0 != err)
            {
                return err;
            }
            Buffer attr_content(const_cast<char*>(attr.bytes_value.data()), 0, attr.bytes_value.size());
            GeoPoint point;
            point.Decode(attr_content);
            x = point.x;
            y = point.y;
        }
        ZSetCacheElementSet subset;
        if (options.in_members)
        {
            StringSet::const_iterator sit = options.submembers.begin();
            while (sit != options.submembers.end())
            {
                ZSetCaheElement ele;
                ValueData score, attr;
                if (0 == ZGetNodeValue(db, key, *sit, score, attr))
                {
                    ele.score = score.NumberValue();
                    Buffer buf1, buf2;
                    ValueData vv;
                    vv.SetValue(*sit, true);
                    vv.Encode(buf1);
                    attr.Encode(buf2);
                    ele.value.assign(buf1.GetRawReadBuffer(), buf1.ReadableBytes());
                    ele.attr.assign(buf2.GetRawReadBuffer(), buf2.ReadableBytes());
                    subset.insert(ele);
                }
                sit++;
            }
        }
        GeoHashHelper::GetAreasByRadius(GEO_MERCATOR_TYPE, y, x, options.radius, ress);
        /*
         * Merge areas if possible to avoid disk search
         */
        std::vector<ZRangeSpec> range_array;
        GeoHashBitsSet::iterator rit = ress.begin();
        GeoHashBitsSet::iterator next_it = ress.begin();
        next_it++;
        while (rit != ress.end())
        {
            GeoHashBits& hash = *rit;
            GeoHashBits next = hash;
            next.bits++;
            while (next_it != ress.end() && next.bits == next_it->bits)
            {
                next.bits++;
                next_it++;
                rit++;
            }
            ZRangeSpec range;
            range.contain_min = true;
            range.contain_max = false;
            range.min.SetIntValue(GeoHashHelper::Allign52Bits(hash));
            range.max.SetIntValue(GeoHashHelper::Allign52Bits(next));
            range_array.push_back(range);
            rit++;
            next_it++;
        }
        DEBUG_LOG("After areas merging, reduce searching area size from %u to %u", ress.size(), range_array.size());

        GeoPointArray points;
        std::vector<ZRangeSpec>::iterator hit = range_array.begin();
        Iterator* iter = NULL;
        while (hit != range_array.end())
        {
            ZSetQueryOptions z_options;
            z_options.withscores = false;
            z_options.withattr = true;
            ValueDataArray values;
            if (options.in_members)
            {
                ZSetCache::GetRangeInZSetCache(subset, *hit, z_options.withscores, z_options.withattr,
                        ZSetValueStoreCallback, &values);
            }
            else
            {
                ZRangeByScoreRange(db, key, *hit, iter, z_options, true, ZSetValueStoreCallback, &values);
            }
            ValueDataArray::iterator vit = values.begin();
            while (vit != values.end())
            {
                GeoPoint point;
                vit->ToString(point.value);
//.........这里部分代码省略.........

void
RadialLineOfSightNode::compute_fill(osg::Node* node, bool backgroundThread)
{
    if ( !getMapNode() )
        return;

    GeoPoint centerMap;
    _center.transform( getMapNode()->getMapSRS(), centerMap );
    centerMap.toWorld( _centerWorld, getMapNode()->getTerrain() );

    bool isProjected = getMapNode()->getMapSRS()->isProjected();

    osg::Vec3d up = isProjected ? osg::Vec3d(0,0,1) : osg::Vec3d(_centerWorld);
    up.normalize();

    //Get the "side" vector
    osg::Vec3d side = isProjected ? osg::Vec3d(1,0,0) : up ^ osg::Vec3d(0,0,1);

    //Get the number of spokes
    double delta = osg::PI * 2.0 / (double)_numSpokes;
    
    osg::Geometry* geometry = new osg::Geometry;
    geometry->setUseVertexBufferObjects(true);

    osg::Vec3Array* verts = new osg::Vec3Array();
    verts->reserve(_numSpokes * 2);
    geometry->setVertexArray( verts );

    osg::Vec4Array* colors = new osg::Vec4Array();
    colors->reserve( _numSpokes * 2 );

    geometry->setColorArray( colors );
    geometry->setColorBinding(osg::Geometry::BIND_PER_VERTEX);

    osg::ref_ptr<osgUtil::IntersectorGroup> ivGroup = new osgUtil::IntersectorGroup();

    for (unsigned int i = 0; i < (unsigned int)_numSpokes; i++)
    {
        double angle = delta * (double)i;
        osg::Quat quat(angle, up );
        osg::Vec3d spoke = quat * (side * _radius);
        osg::Vec3d end = _centerWorld + spoke;        
        osg::ref_ptr<DPLineSegmentIntersector> dplsi = new DPLineSegmentIntersector( _centerWorld, end );
        ivGroup->addIntersector( dplsi.get() );
    }

    osgUtil::IntersectionVisitor iv;
    iv.setIntersector( ivGroup.get() );

    node->accept( iv );

    for (unsigned int i = 0; i < (unsigned int)_numSpokes; i++)
    {
        //Get the current hit
        DPLineSegmentIntersector* los = dynamic_cast<DPLineSegmentIntersector*>(ivGroup->getIntersectors()[i].get());
        DPLineSegmentIntersector::Intersections& hits = los->getIntersections();

        osg::Vec3d currEnd = los->getEnd();
        bool currHasLOS = hits.empty();
        osg::Vec3d currHit = currHasLOS ? osg::Vec3d() : hits.begin()->getWorldIntersectPoint();

        //Get the next hit
        unsigned int nextIndex = i + 1;
        if (nextIndex == _numSpokes) nextIndex = 0;
        DPLineSegmentIntersector* losNext = static_cast<DPLineSegmentIntersector*>(ivGroup->getIntersectors()[nextIndex].get());
        DPLineSegmentIntersector::Intersections& hitsNext = losNext->getIntersections();

        osg::Vec3d nextEnd = losNext->getEnd();
        bool nextHasLOS = hitsNext.empty();
        osg::Vec3d nextHit = nextHasLOS ? osg::Vec3d() : hitsNext.begin()->getWorldIntersectPoint();
        
        if (currHasLOS && nextHasLOS)
        {
            //Both rays have LOS            
            verts->push_back( _centerWorld - _centerWorld );
            colors->push_back( _goodColor );
            
            verts->push_back( nextEnd - _centerWorld );
            colors->push_back( _goodColor );
            
            verts->push_back( currEnd - _centerWorld );
            colors->push_back( _goodColor );
        }        
        else if (!currHasLOS && !nextHasLOS)
        {
            //Both rays do NOT have LOS

            //Draw the "good triangle"            
            verts->push_back( _centerWorld - _centerWorld );
            colors->push_back( _goodColor );
            
            verts->push_back( nextHit - _centerWorld );
            colors->push_back( _goodColor );
            
            verts->push_back( currHit - _centerWorld );                       
            colors->push_back( _goodColor );

            //Draw the two bad triangles
            verts->push_back( currHit - _centerWorld );
            colors->push_back( _badColor );
//.........这里部分代码省略.........

PyObject* xcsoar_Airspaces_addPolygon(Pyxcsoar_Airspaces *self, PyObject *args) {
  PyObject *py_points = nullptr,
           *py_name = nullptr,
           *py_as_class = nullptr,
           *py_base_ref = nullptr,
           *py_top_ref = nullptr;
  double base_alt, top_alt;

  if (!PyArg_ParseTuple(args, "OOOdOdO", &py_points, &py_name, &py_as_class,
                                         &base_alt, &py_base_ref,
                                         &top_alt, &py_top_ref)) {
    return nullptr;
  }

  /* Parse points */
  std::vector<GeoPoint> points;

  if (!PySequence_Check(py_points)) {
    PyErr_SetString(PyExc_ValueError, "First argument is no sequence");
    return nullptr;
  }

  Py_ssize_t num_items = PySequence_Fast_GET_SIZE(py_points);

  for (Py_ssize_t i = 0; i < num_items; ++i) {
    PyObject *py_location = PySequence_Fast_GET_ITEM(py_points, i);

    GeoPoint location = Python::ReadLonLat(py_location);

    if (!location.IsValid()) {
      if (PyErr_Occurred() == nullptr)
        PyErr_SetString(PyExc_RuntimeError, "Unknown error while parsing location");

      return nullptr;
    }

    points.push_back(location);
  }

  if (points.size() < 3) {
    PyErr_SetString(PyExc_ValueError, "Polygon has not enough points");
    return nullptr;
  }

  /* Parse airspace name */
  tstring name;

  if (!Python::PyStringToString(py_name, name)) {
    PyErr_SetString(PyExc_ValueError, "Can't parse airspace name.");
    return nullptr;
  }

  /* Parse airspace class */
  tstring as_class;
  AirspaceClass type = AirspaceClass::OTHER;

  if (!Python::PyStringToString(py_as_class, as_class)) {
    PyErr_SetString(PyExc_ValueError, "Can't parse airspace class.");
    return nullptr;
  }

  for (unsigned i = 0; i < ARRAY_SIZE(airspace_class_strings); i++) {
    if (as_class.compare(airspace_class_strings[i].string) == 0)
      type = airspace_class_strings[i].type;
  }

  /* Parse airspace base and top */
  tstring base_ref, top_ref;
  AirspaceAltitude base, top;

  if (!Python::PyStringToString(py_base_ref, base_ref)) {
    PyErr_SetString(PyExc_ValueError, "Can't parse airspace base reference.");
    return nullptr;
  }

  if (!Python::PyStringToString(py_top_ref, top_ref)) {
    PyErr_SetString(PyExc_ValueError, "Can't parse airspace top reference.");
    return nullptr;
  }

  if (base_ref.compare("MSL") == 0) {
    base.reference = AltitudeReference::MSL;
    base.altitude = base_alt;
  } else if (base_ref.compare("FL") == 0) {
    base.reference = AltitudeReference::STD;
    base.flight_level = base_alt;
  } else if (base_ref.compare("AGL") == 0) {
    base.reference = AltitudeReference::AGL;
    base.altitude_above_terrain = base_alt;
  } else {
    PyErr_SetString(PyExc_ValueError, "Can't parse airspace base.");
    return nullptr;
  }

  if (top_ref.compare("MSL") == 0) {
    top.reference = AltitudeReference::MSL;
    top.altitude = top_alt;
  } else if (top_ref.compare("FL") == 0) {
    top.reference = AltitudeReference::STD;
    top.flight_level = top_alt;
//.........这里部分代码省略.........

GeoVector::GeoVector(const GeoPoint &source, const GeoPoint &target)
{
  *this = source.DistanceBearing(target);
}

osg::Geometry*
GeometryPool::createGeometry(const TileKey& tileKey,
                             const MapInfo& mapInfo,
                             MaskGenerator* maskSet) const
{    
    // Establish a local reference frame for the tile:
    osg::Vec3d centerWorld;
    GeoPoint centroid;
    tileKey.getExtent().getCentroid( centroid );
    centroid.toWorld( centerWorld );

    osg::Matrix world2local, local2world;
    centroid.createWorldToLocal( world2local );
    local2world.invert( world2local );

    // Attempt to calculate the number of verts in the surface geometry.
    bool createSkirt = _options.heightFieldSkirtRatio() > 0.0f;

    unsigned numVertsInSurface    = (_tileSize*_tileSize);
    unsigned numVertsInSkirt      = createSkirt ? _tileSize*4u - 4u : 0;
    unsigned numVerts             = numVertsInSurface + numVertsInSkirt;    
    unsigned numIndiciesInSurface = (_tileSize-1) * (_tileSize-1) * 6;
    unsigned numIncidesInSkirt    = getNumSkirtElements();
    
    // TODO: reconsider this ... 
    GLenum mode = (_options.gpuTessellation() == true) ? GL_PATCHES : GL_TRIANGLES;

    // Pre-allocate enough space for all triangles.
    osg::DrawElements* primSet = new osg::DrawElementsUShort(mode);

    primSet->reserveElements(numIndiciesInSurface + numIncidesInSkirt);

    osg::BoundingSphere tileBound;

    // the geometry:
    osg::Geometry* geom = new osg::Geometry();
    geom->setUseVertexBufferObjects(true);
    geom->setUseDisplayList(false);

    geom->addPrimitiveSet( primSet );

    // the vertex locations:
    osg::Vec3Array* verts = new osg::Vec3Array();
    verts->reserve( numVerts );
    geom->setVertexArray( verts );

    // the surface normals (i.e. extrusion vectors)
    osg::Vec3Array* normals = new osg::Vec3Array();
    normals->reserve( numVerts );
    geom->setNormalArray( normals );
    geom->setNormalBinding( geom->BIND_PER_VERTEX );

#if 0
    // colors
    osg::Vec4Array* colors = new osg::Vec4Array();
    colors->push_back(osg::Vec4f(1,1,1,1));
    geom->setColorArray(colors);
    geom->setColorBinding(osg::Geometry::BIND_OVERALL);
#endif

    osg::Vec3Array* neighbors = 0L;
    if ( _options.morphTerrain() == true )
    {
        // neighbor positions (for morphing)
        neighbors = new osg::Vec3Array();
        neighbors->reserve( numVerts );
        geom->setTexCoordArray( 1, neighbors );
    }

    // tex coord is [0..1] across the tile. The 3rd dimension tracks whether the
    // vert is masked: 0=yes, 1=no
#ifdef SHARE_TEX_COORDS
    bool populateTexCoords = false;
    if ( !_sharedTexCoords.valid() )
    {
        _sharedTexCoords = new osg::Vec3Array();
        _sharedTexCoords->reserve( numVerts );
        populateTexCoords = true;
    }    
    osg::Vec3Array* texCoords = _sharedTexCoords.get();
#else
    bool populateTexCoords = true;
    osg::Vec3Array* texCoords = new osg::Vec3Array();
    texCoords->reserve( numVerts );
#endif

    geom->setTexCoordArray( 0, texCoords );
    
    float delta = 1.0/(_tileSize-1);
    osg::Vec3d tdelta(delta,0,0);
    tdelta.normalize();
    osg::Vec3d vZero(0,0,0);

    osg::ref_ptr<GeoLocator> locator = GeoLocator::createForKey( tileKey, mapInfo );

    for(unsigned row=0; row<_tileSize; ++row)
    {
        float ny = (float)row/(float)(_tileSize-1);
        for(unsigned col=0; col<_tileSize; ++col)
        {
//.........这里部分代码省略.........

bool
LocalizedNode::updateTransforms( const GeoPoint& p, osg::Node* patch )
{
    if ( p.isValid() )
    {
        GeoPoint absPos(p);
        if ( !makeAbsolute(absPos, patch) )
            return false;

        OE_DEBUG << LC << "Update transforms for position: " << absPos.x() << ", " << absPos.y() << ", " << absPos.z()
            << std::endl;

        osg::Matrixd local2world;
        absPos.createLocalToWorld( local2world );
        
        // apply the local offsets
        local2world.preMult( osg::Matrix::translate(_localOffset) );

        if ( _autoTransform )
        {
            static_cast<osg::AutoTransform*>(_xform.get())->setPosition( local2world.getTrans() );
            static_cast<osg::AutoTransform*>(_xform.get())->setScale( _scale );
            static_cast<osg::AutoTransform*>(_xform.get())->setRotation( _localRotation );
        }
        else
        {
            static_cast<osg::MatrixTransform*>(_xform.get())->setMatrix( 
                osg::Matrix::scale(_scale) * 
                osg::Matrix::rotate(_localRotation) *
                local2world  );
        }

        
        CullNodeByHorizon* culler = dynamic_cast<CullNodeByHorizon*>(_xform->getCullCallback());
        if ( culler )
            culler->_world = local2world.getTrans();
    }
    else
    {
        osg::Vec3d absPos = p.vec3d() + _localOffset;

        if ( _autoTransform )
        {
            static_cast<osg::AutoTransform*>(_xform.get())->setPosition( absPos );
            static_cast<osg::AutoTransform*>(_xform.get())->setScale( _scale );
            static_cast<osg::AutoTransform*>(_xform.get())->setRotation( _localRotation );
        }
        else
        {
            static_cast<osg::MatrixTransform*>(_xform.get())->setMatrix(
                osg::Matrix::scale(_scale) * 
                osg::Matrix::rotate(_localRotation) *
                osg::Matrix::translate(absPos) );
        }
    }
    

    dirtyBound();

    return true;
}

bool
ElevationQuery::getElevationImpl(const GeoPoint& point,
                                 float&          out_elevation,
                                 double          desiredResolution,
                                 double*         out_actualResolution)
{
    // assertion.
    if ( !point.isAbsolute() )
    {
        OE_WARN << LC << "Assertion failure; input must be absolute" << std::endl;
        return false;
    }

    osg::Timer_t begin = osg::Timer::instance()->tick();

    // first try the terrain patches.
    if ( _patchLayers.size() > 0 )
    {
        osgUtil::IntersectionVisitor iv;

        if ( _ivrc.valid() )
            iv.setReadCallback(_ivrc.get());

        for(std::vector<ModelLayer*>::iterator i = _patchLayers.begin(); i != _patchLayers.end(); ++i)
        {
            // find the scene graph for this layer:
            osg::Node* node = (*i)->getSceneGraph( _mapf.getUID() );
            if ( node )
            {
                // configure for intersection:
                osg::Vec3d surface;
                point.toWorld( surface );

                // trivial bounds check:
                if ( node->getBound().contains(surface) )
                {
                    osg::Vec3d nvector;
                    point.createWorldUpVector(nvector);

                    osg::Vec3d start( surface + nvector*5e5 );
                    osg::Vec3d end  ( surface - nvector*5e5 );

                    // first time through, set up the intersector on demand
                    if ( !_patchLayersLSI.valid() )
                    {
                        _patchLayersLSI = new DPLineSegmentIntersector(start, end);
                        _patchLayersLSI->setIntersectionLimit( _patchLayersLSI->LIMIT_NEAREST );
                    }
                    else
                    {
                        _patchLayersLSI->reset();
                        _patchLayersLSI->setStart( start );
                        _patchLayersLSI->setEnd  ( end );
                    }

                    // try it.
                    iv.setIntersector( _patchLayersLSI.get() );
                    node->accept( iv );

                    // check for a result!!
                    if ( _patchLayersLSI->containsIntersections() )
                    {
                        osg::Vec3d isect = _patchLayersLSI->getIntersections().begin()->getWorldIntersectPoint();

                        // transform back to input SRS:
                        GeoPoint output;
                        output.fromWorld( point.getSRS(), isect );
                        out_elevation = (float)output.z();
                        if ( out_actualResolution )
                            *out_actualResolution = 0.0;

                        return true;
                    }
                }
                else
                {
                    //OE_INFO << LC << "Trivial rejection (bounds check)" << std::endl;
                }
            }
        }
    }

    if ( _mapf.elevationLayers().empty() )
    {
        // this means there are no heightfields.
        out_elevation = 0.0;
        return true;
    }

    // tile size (resolution of elevation tiles)
    unsigned tileSize = 257; // yes?

    // default LOD:
    unsigned lod = 23u;

    // attempt to map the requested resolution to an LOD:
    if (desiredResolution > 0.0)
    {
        int level = _mapf.getProfile()->getLevelOfDetailForHorizResolution(desiredResolution, tileSize);
        if ( level > 0 )
//.........这里部分代码省略.........

bool
ElevationQuery::getElevationImpl(const GeoPoint& point,
                                 double&         out_elevation,
                                 double          desiredResolution,
                                 double*         out_actualResolution)
{
    osg::Timer_t start = osg::Timer::instance()->tick();

    if ( _maxDataLevel == 0 || _tileSize == 0 )
    {
        // this means there are no heightfields.
        out_elevation = 0.0;
        return true;
    }

    
    //This is the max resolution that we actually have data at this point
    unsigned int bestAvailLevel = getMaxLevel( point.x(), point.y(), point.getSRS(), _mapf.getProfile());

    if (desiredResolution > 0.0)
    {
        unsigned int desiredLevel = _mapf.getProfile()->getLevelOfDetailForHorizResolution( desiredResolution, _tileSize );
        if (desiredLevel < bestAvailLevel) bestAvailLevel = desiredLevel;
    }

    OE_DEBUG << "Best available data level " << point.x() << ", " << point.y() << " = "  << bestAvailLevel << std::endl;

    // transform the input coords to map coords:
    GeoPoint mapPoint = point;
    if ( point.isValid() && !point.getSRS()->isEquivalentTo( _mapf.getProfile()->getSRS() ) )
    {
        mapPoint = point.transform(_mapf.getProfile()->getSRS());
        if ( !mapPoint.isValid() )
        {
            OE_WARN << LC << "Fail: coord transform failed" << std::endl;
            return false;
        }
    }

    osg::ref_ptr<osg::HeightField> tile;

    // get the tilekey corresponding to the tile we need:
    TileKey key = _mapf.getProfile()->createTileKey( mapPoint.x(), mapPoint.y(), bestAvailLevel );
    if ( !key.valid() )
    {
        OE_WARN << LC << "Fail: coords fall outside map" << std::endl;
        return false;
    }

    // Check the tile cache. Note that the TileSource already likely has a MemCache
    // attached to it. We employ a secondary cache here for a couple reasons. One, this
    // cache will store not only the heightfield, but also the tesselated tile in the event
    // that we're using GEOMETRIC mode. Second, since the call the getHeightField can 
    // fallback on a lower resolution, this cache will hold the final resolution heightfield
    // instead of trying to fetch the higher resolution one each item.

    TileCache::Record record = _tileCache.get( key );
    if ( record.valid() )
        tile = record.value().get();

    // if we didn't find it, build it.
    if ( !tile.valid() )
    {
        // generate the heightfield corresponding to the tile key, automatically falling back
        // on lower resolution if necessary:
        _mapf.getHeightField( key, true, tile, 0L );

        // bail out if we could not make a heightfield a all.
        if ( !tile.valid() )
        {
            OE_WARN << LC << "Unable to create heightfield for key " << key.str() << std::endl;
            return false;
        }

        _tileCache.insert(key, tile.get());
    }

    OE_DEBUG << LC << "LRU Cache, hit ratio = " << _tileCache.getStats()._hitRatio << std::endl;

    // see what the actual resolution of the heightfield is.
    if ( out_actualResolution )
        *out_actualResolution = (double)tile->getXInterval();

    bool result = true;

    const GeoExtent& extent = key.getExtent();
    double xInterval = extent.width()  / (double)(tile->getNumColumns()-1);
    double yInterval = extent.height() / (double)(tile->getNumRows()-1);
    
    out_elevation = (double) HeightFieldUtils::getHeightAtLocation( 
        tile.get(), 
        mapPoint.x(), mapPoint.y(), 
        extent.xMin(), extent.yMin(), 
        xInterval, yInterval, _mapf.getMapInfo().getElevationInterpolation() );

    osg::Timer_t end = osg::Timer::instance()->tick();
    _queries++;
    _totalTime += osg::Timer::instance()->delta_s( start, end );

    return result;
//.........这里部分代码省略.........

void
RasterMap::ScanLine(const GeoPoint &start, const GeoPoint &end,
                    short *buffer, unsigned size, bool interpolate) const
{
  assert(buffer != NULL);
  assert(size > 0);

  const short invalid = RasterBuffer::TERRAIN_INVALID;

  const fixed total_distance = start.DistanceS(end);
  if (!positive(total_distance)) {
    std::fill_n(buffer, size, invalid);
    return;
  }

  /* clip the line to the map bounds */

  GeoPoint clipped_start = start, clipped_end = end;
  const GeoClip clip(GetBounds());
  if (!clip.ClipLine(clipped_start, clipped_end)) {
    std::fill_n(buffer, size, invalid);
    return;
  }

  fixed clipped_start_distance =
    std::max(clipped_start.DistanceS(start), fixed(0));
  fixed clipped_end_distance =
    std::max(clipped_end.DistanceS(start), fixed(0));

  /* calculate the offsets of the clipped range within the buffer */

  unsigned clipped_start_offset =
    (unsigned)(size * clipped_start_distance / total_distance);
  unsigned clipped_end_offset =
    uround(size * clipped_end_distance / total_distance);
  if (clipped_end_offset > size)
    clipped_end_offset = size;
  if (clipped_start_offset + 2 > clipped_end_offset) {
    std::fill_n(buffer, size, invalid);
    return;
  }

  assert(clipped_start_offset < size);
  assert(clipped_end_offset <= size);

  /* fill the two regions which are outside the map  */

  std::fill(buffer, buffer + clipped_start_offset, invalid);
  std::fill(buffer + clipped_end_offset, buffer + size, invalid);

  /* now scan the middle part which is within the map */

  const unsigned max_x = raster_tile_cache.GetFineWidth();
  const unsigned max_y = raster_tile_cache.GetFineHeight();

  RasterLocation raster_start = projection.ProjectFine(clipped_start);
  if (raster_start.x >= max_x)
    raster_start.x = max_x - 1;
  if (raster_start.y >= max_y)
    raster_start.y = max_y - 1;

  RasterLocation raster_end = projection.ProjectFine(clipped_end);
  if (raster_end.x >= max_x)
    raster_end.x = max_x - 1;
  if (raster_end.y >= max_y)
    raster_end.y = max_y - 1;

  raster_tile_cache.ScanLine(raster_start, raster_end,
                             buffer + clipped_start_offset,
                             clipped_end_offset - clipped_start_offset,
                             interpolate);
}

bool
OLCTriangle::FindClosingPairs(unsigned old_size)
{
  if (predict) {
    return closing_pairs.insert(ClosingPair(0, n_points-1));
  }

  struct TracePointNodeAccessor {
    gcc_pure
    int GetX(const TracePointNode &node) const {
      return node.point->GetFlatLocation().longitude;
    }

    gcc_pure
    int GetY(const TracePointNode &node) const {
      return node.point->GetFlatLocation().latitude;
    }
  };

  QuadTree<TracePointNode, TracePointNodeAccessor> search_point_tree;

  for (unsigned i = old_size; i < n_points; ++i) {
    TracePointNode node;
    node.point = &GetPoint(i);
    node.index = i;

    search_point_tree.insert(node);
  }

  search_point_tree.Optimise();

  bool new_pair = false;

  for (unsigned i = old_size; i < n_points; ++i) {
    TracePointNode point;
    point.point = &GetPoint(i);
    point.index = i;

    const unsigned max_range =
      trace_master.ProjectRange(GetPoint(i).GetLocation(), max_distance);

    const GeoPoint start = GetPoint(i).GetLocation();
    const int min_altitude = GetMinimumFinishAltitude(GetPoint(i));
    const int max_altitude = GetMaximumStartAltitude(GetPoint(i));

    unsigned last = 0, first = i;

    const auto visitor = [this, i, start,
                          min_altitude, max_altitude,
                          &first, &last]
      (const TracePointNode &node) {
      const SearchPoint dest = GetPoint(node.index);

      if (node.index + 2 < i &&
          GetPoint(node.index).GetIntegerAltitude() <= max_altitude &&
          start.Distance(dest.GetLocation()) <= max_distance) {
        // point i is last point
        first = std::min(node.index, first);
        last = i;
      } else if (node.index > i + 2 &&
                 GetPoint(node.index).GetIntegerAltitude() >= min_altitude &&
                 start.Distance(dest.GetLocation()) <= max_distance) {
        // point i is first point
        first = i;
        last = std::max(node.index, last);
      }
    };

    search_point_tree.VisitWithinRange(point, max_range, visitor);

    if (last != 0 && closing_pairs.insert(ClosingPair(first, last)))
      new_pair = true;
  }

  return new_pair;
}

void
ModelSplatter::operator()(const TileKey& key, osg::Node* node)
{
    TerrainTileNode* tile = osgEarth::findTopMostNodeOfType<TerrainTileNode>(node);
    if ( !tile )
        return;

    if ( key.getLOD() >= _minLOD && _model.valid() )
    {
        // make sure the correct model is loaded
        establish();

        // elevation texture and matrix are required
        osg::Texture* elevationTex = tile->getElevationTexture();
        if ( !elevationTex )
        {
            //OE_WARN << LC << "No elevation texture for key " << key.str() << "\n";
            return;
        }

        osg::RefMatrix* elevationTexMat = tile->getElevationTextureMatrix();
        if ( !elevationTexMat )
        {
            //OE_WARN << LC << "No elevation texture matrix for key " << key.str() << "\n";
            return;
        }
        
        tile->addChild( _model.get() );

        osg::StateSet* ss = tile->getOrCreateStateSet();

        // first, a rotation vector to make trees point up.
        GeoPoint p;
        key.getExtent().getCentroid(p);
        osg::Vec3d up;
        p.createWorldUpVector(up);
        osg::Quat q;
        q.makeRotate(osg::Vec3d(0,0,1), up);
        osg::Matrixd zup = osg::Matrixd::rotate(q);

        // matrices to resolve the weird terrain localization into a usable LTP.
        osg::Matrix tile2world = tile->getMatrix();
        osg::Matrix world2ltp;
        p.createWorldToLocal(world2ltp);
        osg::Matrix local2ltp = tile2world * world2ltp;
        osg::Matrix ltp2local;
        ltp2local.invert(local2ltp);

        // after inverting the matrix, combine the ZUP (optimization)
        local2ltp.preMult( zup );

        ss->addUniform( new osg::Uniform("oe_trees_local2ltp", osg::Matrixf(local2ltp)) );
        ss->addUniform( new osg::Uniform("oe_trees_ltp2local", osg::Matrixf(ltp2local)) );

        // calculate the scatter area:
        float h = key.getExtent().height() * 111320.0f;
        float w = key.getExtent().width() * 111320.0f * cos(fabs(osg::DegreesToRadians(p.y())));
        ss->addUniform( new osg::Uniform("oe_trees_span", osg::Vec2f(w,h)) );
        
        ss->setTextureAttributeAndModes(2, tile->getElevationTexture(), 1);
        ss->addUniform( new osg::Uniform("oe_terrain_tex_matrix", osg::Matrixf(*elevationTexMat)) );        
    }
}

 /**
  * distance from this to the reference
  */
 fixed DistanceTo(const GeoPoint &ref) const {
   return reference.Distance(ref);
 }