高度图生成算法？

Posted 2023-02-18

【中文标题】高度图生成算法？

【英文标题】：Heightmap generation algorithm?

【发布时间】：2010-10-08 03:13:57

【问题描述】：

我在互联网上四处寻找，找不到解决这个特定问题的完美算法：

我们的客户有一组点和重量数据以及每个点，如下图所示：

weighted points http://chakrit.net/files/***/so_heightmap_points.png

其中，我们有一个 GIS 程序，可以从这些点及其权重值生成“高度图”或某种地形数据，但由于我们有近一千个数据点并且这些数据会随着时间而变化，所以我们希望创建我们自己的工具来自动生成这些高度图。

到目前为止，我已尝试使用Sqrt((x1 - x2) ^ 2 + (y1 - y2) ^ 2) 计算每个像素到最近数据点的距离的权重，并将权重和距离因子应用于数据点的颜色以生成该特定像素的渐变颜色：

heightmap result http://chakrit.net/files/***/so_heightmap_result.png

您可以看到，某些数据点的配置仍然存在问题，当数据点很多时，算法有时会生成一个相当多边形的图像。理想的结果应该看起来更像一个省略号，而不是一个多边形。

这是来自***关于梯度上升的文章的一个示例图片，它展示了我想要的结果：

mountains http://chakrit.net/files/***/so_gradient_descent.png

梯度上升算法不是我感兴趣的。我对什么感兴趣；是首先计算该图片中的原始函数的算法，提供具有权重的数据点。

我没有上过任何拓扑数学课程，但我可以做一些微积分。我想我可能遗漏了一些东西，并且不知道应该在那个 Google 搜索框中输入什么内容。

我需要一些指示。

谢谢！

【问题讨论】：

您希望生成的地形有多平滑？你有负位移和正位移吗？我假设重量不是指直接占据的“小山”而是最大点。

@ShuggyCoUk 1. 尽可能流畅，没有太多复杂性。 2. 所有值都是正数。 3. 是的，它可以解释为山的最大点并且仍然是正确的。

@chakrit 图片链接已损坏，请尽可能修复它们

【参考方案1】：

我知道这是一个相当老的问题，但我在尝试解决类似问题时偶然发现了它。

有一个名为Surfit 的开源项目正是实现了这种类型的功能。

【讨论】：

哦，它永远不会过时，因为即使在我离开公司后，也必须有人维护它，我只需要把那个人联系起来，其他人也可能会觉得它有用:)

【参考方案2】：

我不久前在 Winamp AVS 中实现了类似的东西。它使用“metaballs”类型的方法来计算每个数据点的平方反比距离（以避免速度的 sqrt），将其封顶（例如到 1.0），并为 2D 网格上的每个点计算这些距离的总和。这将给出一个平滑变化的颜色/高度图。

如果您想查看代码，它位于我的J10 AVS pack 的“Glowy”预设中。

编辑：只是看着它，我添加了一些其他爵士乐以使其看起来更漂亮，最重要的部分是：

d1=s/(sqr(px1-rx)+sqr(py1-ry));
d2=s/(sqr(px2-rx)+sqr(py2-ry));
d3=s/(sqr(px3-rx)+sqr(py3-ry));
d4=s/(sqr(px4-rx)+sqr(py4-ry));
d5=s/(sqr(px5-rx)+sqr(py5-ry));
d6=s/(sqr(px6-rx)+sqr(py6-ry));
d=d1+d2+d3+d4+d5+d6;

这需要 6 个点的总和。对红色、绿色和蓝色输出值所做的一切都是为了让它看起来更漂亮。 6 分不算多，但请记住，当它是新机器时，我试图在 400MHz 机器上的 320x200 网格上实时运行它（它以 ~20fps 的速度运行）。 :)

用 red = d 替换 red =、green = 和 blue = ... 行；等等......看看我的意思。所有的美感都消失了，你只剩下数据点周围平滑变化的斑点的灰度图像。

另一个编辑：我忘了说“s”是所有点的共享权重，为每个点更改它会赋予每个点单独的权重，例如d1 = 2/(...) 和 d2 = 1/(...) 将使 d1 在其中心的高度是 d2 的两倍。您可能还想用 d1 = 2/max(..., 1.0) 之类的东西来限制底部的表达式，以平滑点的顶部，这样它们就不会在中间的无穷大处达到峰值。 :)

抱歉，答案很混乱...我认为发布代码示例就足够了，但经过检查，我的代码令人困惑且难以阅读。 :(

【参考方案3】：

您已询问有关不规则数据的二维插值算法的信息，这是一个相当复杂的领域。既然您说您拥有 ArcGIS，我强烈建议您使用其内置的 features 在 ArcGIS 中插入 features 以进行自动计算。我相信这将比编写自己的插值算法容易。我已经对 ArcGIS 进行了一些自动化操作，它相当简单。

如果您确实编写了自己的插值代码（我建议您不要这样做），首先要选择合适的算法，因为有几种算法各有优缺点。这是从优秀插值工具Surfer 的帮助中摘录的一些建议（顺便说一句，它也可以很容易地自动化）。还有更多算法，这些只是我尝试过的。

克里金法

是更灵活的方法之一，可用于对几乎任何类型的数据集进行网格化。对于大多数数据集，使用默认线性变差函数的克里金法非常有效。一般来说，我们最常推荐这种方法。克里金法是默认的网格化方法，因为它可以为大多数数据集生成良好的地图。对于较大的数据集，克里金法可能会相当慢。克里金法可以推断超出数据 Z 范围的网格值。

反距离加权

速度很快，但倾向于在数据点周围生成同心轮廓的“靶心”模式。幂的反距离不会外推超出数据范围的 Z 值。一个简单的逆距离加权算法很容易实现，但速度较慢。

线性插值三角剖分

速度很快。当您使用小型数据集时，带有线性插值的三角剖分会在数据点之间生成不同的三角形面。使用线性插值进行三角剖分不会将 Z 值外推到数据范围之外。

Shephard 方法

类似于 Inverse Distance to a Power，但不会产生“靶心”模式，尤其是在使用平滑因子时。

Shepard 方法

可以推断超出数据 Z 范围的值。

要实现算法：您可以尝试谷歌搜索，或点击其他答案中的链接。有一些包含插值的开源 GIS 包，所以如果你喜欢通过 C++ 挖坑，也许你可以从中提取算法。或者大卫·沃森的this book 显然被认为是经典之作，尽管阅读起来很棘手，并且示例代码是意大利面条基本！但是，据我所知，这是最好的。如果 Stack Overflow 上的其他人知道得更好，请纠正我，我也不敢相信。

【讨论】：

实际上，一位操作 ArcGIS 的同事向我提出了这个问题。自动化可能是一个不错的选择，我会尝试一下。谢谢！

顺便说一句，如果工作流不满足您的需求，您可以使用 ArcGIS 宏或编写插件 DLL 等。

【参考方案4】：

Kriging 是执行此操作的重量级方法之一，尤其是在 GIS 领域。它有几个很好的数学属性 - 缺点是它可能会很慢，具体取决于您的variogram。

如果你想要更简单的东西，有很多插值例程可以很好地处理这个问题。如果您能获得一份Numerical Recipes 的副本，第 3 章将专门解释许多插值变体，并包含代码示例及其功能属性的描述。

【讨论】：

据我所知，C 2nd Edition 中的数字配方仅包含一个用于 2D 线性插值的例程。有点限制，您可能还想考虑克里金法或反距离加权或其他答案中建议的其他方法之一。

【参考方案5】：

您正在寻找 Blender 调用“metaballs”（Wikipedia article with links、example）的东西。这样想：

您的对象是从地面伸出的圆锥体。它们都是抛物线，重量表明它们从地面伸出多远。或者，使它们都具有相同的高度并相应地调整抛物线的“平面度”，因此较大的重量会使圆锥体非常宽，而较低的重量会使圆锥体锋利。甚至在某种程度上两者兼而有之。

我建议你实现这个，看看它的样子。

接下来，您需要在结果上挂一块布或橡皮布。布料会拉伸一定量，通常会由于重力而下垂。锥体保持它。

只要你靠近一个圆锥的中心，Z坐标就是圆锥表面上的位置。当您离开锥体中心时，重力开始下降，其他锥体的影响增加。

【讨论】：

他确实希望制作一个二维隐式曲面。可视化的有用方法，但是如何计算呢？ :-)

***文章包含公式和链接。

【参考方案6】：

表面插值似乎是一个困难的数学问题。 另一种更便宜的方法是：

For each pixel:For each point:pixel.addWeight(weight(point, pixel))

def addWeight(w):totalweight += wnumberofweights += 1weight = totalweight / numberofweights

权重函数示例：

def weight(point, pixel):return point.weight * 1/(1 + sqrt((point.x - pixel.x)^2 + (point.y - pixel.y)^2))

这是一种蛮力的方法，但很简单。

【讨论】：

有趣...会尝试一下并回复您。

【参考方案7】：

是计算的算法 那张图片中的原始功能 第一名，提供数据点 带重量。

这是可能的。如果您从单点开始，您将始终以圆形结束，但如果您对数据点进行加权并考虑到这一点，您可以将圆形挤压成椭圆形，如图所示。

最终得到多边形的原因是您在计算中使用了离散函数 - 首先找到最接近的颜色，然后确定颜色。

您应该研究梯度算法，该算法根据将点包围在三角形中的三个数据点的距离和权重为该点分配颜色。

梯度算法

这取决于您要显示的内容。一个简单的算法是：

对于每个像素：

找出构成围绕该像素的最小三角形的三个点

将此点设置为受权重和到每个数据点的距离影响的颜色（HSV 颜色系统）：

pixel.color = datapoint[1].weight * distance(pixel, datapoint[1]) * datapoint[1].color + datapoint[2].weight * distance(pixel, datapoint[2]) * datapoint[2].color + datapoint[3].weight * distance(pixel, datapoint[3]) * datapoint[3].color

我在这里使用 +，但您需要确定适合您的应用程序的“平均”算法。

-亚当

【讨论】：

嗯，这取决于你想要的结果。理想情况下，您会为每个像素考虑宇宙中的每个数据点，但这是处理密集型的，可能不是您真正想要的。但是，它可能是您需要的（例如磁场）

请注意在输出地图中为每个 m 个像素迭代所有 n 个点。这是 O(n*m)，对于 1000x1000 的图像和 1000 个数据点，这是十亿次操作。这不会扩展。使用 Shepherd 算法平铺平面或类似的东西。

【参考方案8】：

您正在寻找的是表面插值。

有些产品可以做到这一点（这里是one）

然后可以以所需的分辨率询问生成的函数/样条线/其他数学结构以提供高度图。

你的插值函数

Sqrt((x1 - x2) ^ 2 + (y1 - y2) ^ 2) 

类似于Inverse Distance Weighted 方法，只是您应用了任意过滤器并丢弃了许多其他数据点。

这些技术中的大多数都依赖于合理数量的样本和支撑这些值的“地形”行为。

我建议使用权重作为高度样本，并在第二个链接中尝试简单的 Shepard 方法（不要过滤任何像素开始），方法是获取样本点对插值点的整体高度值的贡献比例您可以以这些比率混合样本的颜色，也可以为点着色。使用强度（粗略地说是简单 RGB 空间中的灰度）来显示高度或添加黑色轮廓线，就像示例图像一样。

【讨论】：

这应该是我要找的...而且我已经在使用 ArcGIS.. 谢谢！

另请注意，您可以以较低的分辨率进行渲染并进行简单（快速）的双线性插值以生成最终高度图，但是如果您提供低分辨率高度图，您的工具可能会自行执行此操作无论如何。

【参考方案9】：

这个问题并不像表面上看起来那么简单。你的问题是两个区域的边界两边需要有相同的高度，也就是说，给定像素的高度是由不止一个最近邻决定的。

如果我理解正确，您至少需要两种算法（和第三条行话）。

要正确执行此操作，您需要将飞机分解为Voronoi tesselation。

您可能想要使用kd-tree 来帮助您找到最近的邻居。不是采用 O(n^2)，而是将其降低到 O(n log(n))（额外的好处是您的 Voronoi 区域生成阶段在开发中将足够快，可以在高度计算阶段工作）。

现在您有一个将每个点索引到其最近邻居 i 的二维地图，您需要遍历地图上的每个 x,y 点并计算其高度。

要对给定的点 x,y 执行此操作，首先获取其最近的邻居 i 并将其粘贴到列表中，然后收集 Voronoi 图上的所有连续区域。一个简单的方法是使用flood fill 找到该区域中的所有点，然后查看边界并收集其他身份。

使用所有最近邻居的列表，您现在可以正确插值！ （有关插值方案，请参见其他答案）。

【讨论】：

+1 kd-tree 是公平的，但 Voronoi 是什么？ ...看来我毕竟需要复杂的数学...

只要看一下 Voronoi 图文章中的图片，很明显，生成最近邻地图会使图片与文章中显示的图片相似。

“这个问题并不像表面上看起来那么简单”好极了:)

@ShuggyCoUk：感谢分享 Shepherd 的方法链接。几个月来我一直在为这个问题感到困惑，知道什么对 Google 来说是一个很大的帮助！

对于 KDtree + Python 中的反距离加权，请参阅***.com/questions/3104781/…