[原][osgearth]OE地形平整代碼解讀

在FlatteningLayer文件的createHeightField函數中：使用的github在2017年1月份的代碼

if (!geoms.getComponents().empty())
    {
        osg::ref_ptr<osg::HeightField> hf = HeightFieldUtils::createReferenceHeightField(
            ex,
            257, 257,           // base tile size for elevation data
            0u,                 // no border
            true);              // initialize to HAE (0.0) heights

        // Initialize to NO DATA.
        hf->getFloatArray()->assign(hf->getNumColumns()*hf->getNumRows(), NO_DATA_VALUE);

        // Create an elevation query envelope at the LOD we are creating
        osg::ref_ptr<ElevationEnvelope> envelope = _pool->createEnvelope(workingSRS, key.getLOD());

        // Resolve the buffering widths:
        double lineWidthLocal = lineWidth()->as(workingSRS->getUnits());
        double bufferWidthLocal = bufferWidth()->as(workingSRS->getUnits());
        
        if(integrate(key, hf, &geoms, workingSRS, lineWidthLocal, bufferWidthLocal, envelope, progress) || (progress && progress->isCanceled()))
        {
            //double t_create = OE_GET_TIMER(create);
            //OE_INFO << LC << key.str() << " : t=" << t_create << "s\n";

            // If integrate made any changes, return the new heightfield.
            // (Or if the operation was canceled...return it anyway and it 
            // will be discarded).
            return hf.release();
        }
    }

 

 

　創建一個高度場，長寬都是257，邊界為0，高度引用大地水平基准面。

   用默認值初始化高度場

   在自己創建的LOD中創建一個高程查詢信

   解決緩存寬度

   整合新高程

   如果高程有任何改變，返回新的高程圖，高度場。

 

這里integrate調用的函數，調用了integratePolygons函數來創建平整的高程圖，我們看看這里具體怎么操作的

我們來看integratePolygons函數：

// Creates a heightfield that flattens an area intersecting the input polygon geometry.創建一個包含集合多邊形的高度場
    // The height of the area is found by sampling a point internal to the polygon.
    // bufferWidth = width of transition from flat area to natural terrain.
    bool integratePolygons(const TileKey& key, osg::HeightField* hf, const Geometry* geom, const SpatialReference* geomSRS,
                           double bufferWidth, ElevationEnvelope* envelope, ProgressCallback* progress)
    {
        bool wroteChanges = false;

        const GeoExtent& ex = key.getExtent();

        double col_interval = ex.width() / (double)(hf->getNumColumns()-1);
        double row_interval = ex.height() / (double)(hf->getNumRows()-1);

        POINT Pex, P, internalP;

        bool needsTransform = ex.getSRS() != geomSRS;

循環遍歷長寬間隔獲取每個頂點坐標

        for (unsigned col = 0; col < hf->getNumColumns(); ++col)
        {
            Pex.x() = ex.xMin() + (double)col * col_interval;

            for (unsigned row = 0; row < hf->getNumRows(); ++row)
            {
                // check for cancelation periodically
                //if (progress && progress->isCanceled())
                //    return false;

                Pex.y() = ex.yMin() + (double)row * row_interval;

                if (needsTransform)
                    ex.getSRS()->transform(Pex, geomSRS, P);
                else
                    P = Pex;
                
                bool done = false;
                double minD2 = bufferWidth * bufferWidth; // minimum distance(squared) to closest polygon edge

                const Polygon* bestPoly = 0L;

                ConstGeometryIterator giter(geom, false);
                while (giter.hasMore() && !done)
                {
                    const Polygon* polygon = dynamic_cast<const Polygon*>(giter.next());
                    if (polygon)
                    {
                        // Does the point P fall within the polygon?
循環檢查這里是否有點在這些幾何形狀里
                        if (polygon->contains2D(P.x(), P.y()))
                        {
                            // yes, flatten it to the polygon's centroid elevation;
                            // and we're dont with this point.
如果這點就在幾何形狀范圍里，直接跳出檢查
                            done = true;
                            bestPoly = polygon;
                            minD2 = -1.0;
                        }

                        // If not in the polygon, how far to the closest edge?
如果沒在，計算距離邊緣最近的距離的平方

                        else
                        {
                            double D2 = getDistanceSquaredToClosestEdge(P, polygon);
查看獲得值是否在緩存范圍內
                            if (D2 < minD2)
                            {
如果在范圍內，就設置好這個點在緩存內最近的位置，以便后面計算
                                minD2 = D2;
                                bestPoly = polygon;
                            }
                        }
                    }
                }

                if (bestPoly && minD2 != 0.0)
                {

　　判斷這些需要獲取的高程點，有沒有在需要關注的幾何圖形里或者緩沖區范圍內的，如果有就做以下工作，來抬高地形：

                    float h;
                    POINT internalP = getInternalPoint(bestPoly);
                    float elevInternal = envelope->getElevation(internalP.x(), internalP.y());

                    if (minD2 < 0.0)
                    {
                        h = elevInternal;
                    }
                    else
                    {
                        float elevNatural = envelope->getElevation(P.x(), P.y());
                        double blend = clamp(sqrt(minD2)/bufferWidth, 0.0, 1.0); // [0..1] 0=internal, 1=natural
                        h = smoothstep(elevInternal, elevNatural, blend);
                    }

                    hf->setHeight(col, row, h);
                    wroteChanges = true;
                }

            }
        }

        return wroteChanges;
    }

 

真正平整的函數在：integrate函數

   進入FlatteningLayer文件的integratePolygons函數

      先獲取TileKey范圍

      獲取長寬的間隔分別是多大長度

      檢查是否要做SRS轉換，這里看是需要的

     

West -180  xMin    SRS-> -20037508.343

        East  0     xMax

        South -90  xMin    SRS-> -20037508.343

        North 90    yMax

 

循環遍歷長寬間隔獲取每個頂點坐標

      POINT P是這點的世界位置點

        循環檢查這些幾何形狀

              看看點是否在這些幾何形狀里

              如果不在，計算距離邊緣最近的距離的平方

                 查看獲得值是否在平整范圍內

如果在范圍內，就設置好這個點在范圍內最近的位置，以便后面計算

              如果點就在幾何形狀范圍里，直接跳出檢查

      直接點說

就是判斷這些需要獲取的高程點，有沒有在需要關注的幾何圖形里或者緩沖區范圍內的，如果有就做以下工作，來抬高地形：

      先找到幾何圖形內部的一個頂點（多為中心，質心等）

      通過ElevationEnvelope類的getElevation函數計算這個點的高程

      判斷這個要改變高程的點與幾何圖形的位置關系：

如果在幾何圖形內就設置這個點用剛才獲取的高程點（這個方式有待商榷）

         如果在緩沖區上，獲取這個點的現實高程值。獲取當前點在緩沖區上的范圍值，做smoothstep變換。

        然后將之前得到的高程按格子塞入高程圖。

           齊活！