在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 "; // 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函数:
1 // Creates a heightfield that flattens an area intersecting the input polygon geometry.创建一个包含集合多边形的高度场 2 // The height of the area is found by sampling a point internal to the polygon. 3 // bufferWidth = width of transition from flat area to natural terrain. 4 bool integratePolygons(const TileKey& key, osg::HeightField* hf, const Geometry* geom, const SpatialReference* geomSRS, 5 double bufferWidth, ElevationEnvelope* envelope, ProgressCallback* progress) 6 { 7 bool wroteChanges = false; 8 9 const GeoExtent& ex = key.getExtent(); 10 11 double col_interval = ex.width() / (double)(hf->getNumColumns()-1); 12 double row_interval = ex.height() / (double)(hf->getNumRows()-1); 13 14 POINT Pex, P, internalP; 15 16 bool needsTransform = ex.getSRS() != geomSRS; 17
循环遍历长宽间隔获取每个顶点坐标
18 for (unsigned col = 0; col < hf->getNumColumns(); ++col) 19 { 20 Pex.x() = ex.xMin() + (double)col * col_interval; 21 22 for (unsigned row = 0; row < hf->getNumRows(); ++row) 23 { 24 // check for cancelation periodically 25 //if (progress && progress->isCanceled()) 26 // return false; 27 28 Pex.y() = ex.yMin() + (double)row * row_interval; 29 30 if (needsTransform) 31 ex.getSRS()->transform(Pex, geomSRS, P); 32 else 33 P = Pex; 34 35 bool done = false; 36 double minD2 = bufferWidth * bufferWidth; // minimum distance(squared) to closest polygon edge 37 38 const Polygon* bestPoly = 0L; 39 40 ConstGeometryIterator giter(geom, false); 41 while (giter.hasMore() && !done) 42 { 43 const Polygon* polygon = dynamic_cast<const Polygon*>(giter.next()); 44 if (polygon) 45 { 46 // Does the point P fall within the polygon?
循环检查这里是否有点在这些几何形状里 47 if (polygon->contains2D(P.x(), P.y())) 48 { 49 // yes, flatten it to the polygon's centroid elevation; 50 // and we're dont with this point.
如果这点就在几何形状范围里,直接跳出检查 51 done = true; 52 bestPoly = polygon; 53 minD2 = -1.0; 54 } 55 56 // If not in the polygon, how far to the closest edge?
如果没在,计算距离边缘最近的距离的平方
57 else 58 { 59 double D2 = getDistanceSquaredToClosestEdge(P, polygon);
查看获得值是否在缓存范围内 60 if (D2 < minD2) 61 {
如果在范围内,就设置好这个点在缓存内最近的位置,以便后面计算 62 minD2 = D2; 63 bestPoly = polygon; 64 } 65 } 66 } 67 } 68 69 if (bestPoly && minD2 != 0.0) 70 {
判断这些需要获取的高程点,有没有在需要关注的几何图形里或者缓冲区范围内的,如果有就做以下工作,来抬高地形:
71 float h; 72 POINT internalP = getInternalPoint(bestPoly); 73 float elevInternal = envelope->getElevation(internalP.x(), internalP.y()); 74 75 if (minD2 < 0.0) 76 { 77 h = elevInternal; 78 } 79 else 80 { 81 float elevNatural = envelope->getElevation(P.x(), P.y()); 82 double blend = clamp(sqrt(minD2)/bufferWidth, 0.0, 1.0); // [0..1] 0=internal, 1=natural 83 h = smoothstep(elevInternal, elevNatural, blend); 84 } 85 86 hf->setHeight(col, row, h); 87 wroteChanges = true; 88 } 89 90 } 91 } 92 93 return wroteChanges; 94 }
真正平整的函数在: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变换。
然后将之前得到的高程按格子塞入高程图。
齐活!