[osgearth]OE地形平整代码解读

最新推荐文章于 2022-11-29 11:13:16 发布

pplander

最新推荐文章于 2022-11-29 11:13:16 发布

阅读量410

点赞数

分类专栏： osgEarth 文章标签： opengl

osgEarth 专栏收录该内容

1 篇文章 0 订阅

订阅专栏

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

if (!geoms.getComponents().empty())
{
osg::ref_ptr<osg::HeightField> hf = HeightFieldUtils::createReferenceHeightField(
ex,
, , // 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()-);
double row_interval = ex.height() / (double)(hf->getNumRows()-);
POINT Pex, P, internalP;
bool needsTransform = ex.getSRS() != geomSRS;

循环遍历长宽间隔获取每个顶点坐标

for (unsigned col = ; col < hf->getNumColumns(); ++col)
{
Pex.x() = ex.xMin() + (double)col * col_interval;
for (unsigned row = ; 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;
}