高度图生成算法？

Question

我在互联网上搜索了很久，但并没有找到一个完美的算法解决这个问题:

我们的客户有一组带权重数据的点，如下图所示：

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

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

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

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

你可以看到，在某些数据点的配置中仍然存在问题，当有很多数据点时，该算法有时会产生一个相当多边形的图像。理想的结果应该更像一个椭圆，而不是一个多边形。

这里有一张来自维基百科关于梯度上升的文章的图片，展示了我想要的结果：

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

梯度上升算法并不是我的兴趣所在。 我感兴趣的是，提供带权重数据点后计算那个图片中原始函数的算法。

我没有学过拓扑数学课程，但我可以做一些微积分。 我认为我可能遗漏了什么，现在非常迷茫，不知道应该在Google搜索框中输入什么。

我需要一些指导。

谢谢！

- chakrit

---

你希望生成的地形有多光滑？ 你是否具有负位移和正位移。我假设重量指的不是“山丘”直接占据的区域，而是它的最高点。 - ShuggyCoUk

1. 尽可能平稳，不要过于复杂。
2. 所有的值都是正数。
3. 是的，它可以被解释为山顶的最高点，并且仍然是正确的。

- chakrit

@chakrit，图片链接已经失效，请尽可能修复它们。 - Aditya

Answer 1

你需要的是表面插值。

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

然后，得到的函数/样条线/其他数学结构可以在所需的分辨率下进行查询，以提供高度图。

你的插值函数。

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

类似于反距离加权方法，但您正在应用任意过滤器并丢弃许多其他数据点。

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

我建议使用权重作为高度样本，并在第二个链接中尝试简单的Shepard方法（一开始不要过滤任何像素），通过在插值点处将样本点对整体高度值的贡献比例混合，可以按这些比例混合样本的颜色来着色该点。使用强度（简单RGB空间中的灰度）来显示高度或添加黑色等高线，就像示例图像所做的那样。

- ShuggyCoUk

---

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

Answer 2

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

如果我理解正确，您需要至少两种算法（和第三个术语）。

要正确执行此操作，您需要将平面分成 Voronoi tessellation。

您可能希望使用 kd-tree来帮助您找到最近的邻居。这将使其从O(n^2)降至O(n log(n))（额外好处是您的Voronoi区域生成阶段足够快，以便在高度计算阶段中进行开发）。

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

为了对给定的点x,y进行操作，首先要获取其最近邻居i并将其放入列表中，然后收集Voronoi图上所有连续的区域。一种简单的方法是使用flood fill查找该区域中的所有点，然后在边界周围查找并收集其他标识。

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

- Jared Updike

---

@Jared 啊...如果我移除梯度计算，只是将每个像素的颜色设置为最近点的颜色，那么我会得到一个沃罗诺伊图，对吧？ - chakrit

@chakrit：确切地说，这只是一个直观概念的花哨名称。实际上，我会使用一个二维整数数组来跟踪最近邻居而不是颜色本身，但从概念上讲，它是相同的，而你的颜色使得可视化这个东西成为可能。 - Jared Updike

有关 Python 中的 KDtree + 反距离加权，请参见以下链接： https://dev59.com/n3A75IYBdhLWcg3wy8Qi - denis

Answer 3

您询问了关于二维不规则数据插值算法的信息，这是一个相当复杂的领域。既然您说您拥有ArcGIS，我强烈建议您使用其内置的功能进行自动计算，以在ArcGIS中自动插值。我确信这比编写自己的插值算法要简单得多。我已经对ArcGIS进行了一些自动化处理，它非常直观易懂。

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

• Kriging是更加灵活的方法之一，适用于几乎任何类型的数据集。对于大多数数据集来说，使用默认的线性变异函数进行Kriging是非常有效的。总的来说，我们通常会推荐这种方法。Kriging是默认的插值方法，因为它可以为大多数数据集生成良好的地图。对于较大的数据集，Kriging可能会比较慢。Kriging可以将网格值外推到数据的Z范围之外。
• 反距离加权法速度快，但有生成“靶心”式同心轮廓的趋势。幂反距离法不会将Z值外推到数据范围之外。简单的反距离加权算法易于实现，但速度较慢。
• 三角剖分与线性插值速度快。当您使用小数据集时，三角剖分与线性插值在数据点之间生成明显的三角形面。三角剖分与线性插值不会将Z值外推到数据范围之外。
• Shephard's Method类似于幂反距离法，但不会产生“靶心”式模式，特别是在使用平滑因子时。 Shepard's Method可以将值外推到数据的Z范围之外。

要实现算法：您可以尝试使用Google搜索，或者按照其他答案中的链接进行。有一些开源GIS软件包包括插值，因此如果您喜欢通过C++挖掘，可能可以从中提取算法。或者David Watson的这本书据说是经典之作，尽管阅读起来很棘手，示例代码也像意大利面条一样混乱！但是据我所知，它是最好的选择。如果Stack Overflow上还有其他更好的人，请纠正我，因为我也不敢相信这一点。

- MarkJ

Answer 4

Kriging 是其中一种重量级的方法，特别是在GIS领域中。它具有几个很好的数学特性 - 缺点是它可能会因为你的变异函数而变慢。

如果您想要更简单的东西，有许多插值例程可以很好地处理这个问题。如果您可以获得数值计算方法, 第3章专门解释了许多插值变体，并包括代码示例和它们的功能属性描述。

- jasedit

---

据我所记，C语言数值计算方法（Numerical Recipes in C）第二版仅包含2D线性插值例程。这有点受限制，您可能还想考虑Kriging或反距离加权或其他答案中建议的方法之一。 - MarkJ

Answer 5

这是一种计算给定数据点权重下的原始函数的算法。如果你从单个点开始，最终会得到圆形，但如果你考虑数据点的权重并将其纳入考虑，就可以将圆形压缩成椭圆形状，如图所示。

之所以会出现多边形，是因为你在计算中使用了离散函数 - 首先找到最接近的颜色，然后确定颜色。

相反，你应该研究渐变算法，该算法基于三个数据点的距离和权重来为一个点分配颜色，并将其包围在三角形中。

渐变算法取决于你要显示什么内容。一个简单的算法如下：

对于每个像素： - 找到围绕该像素的最小三角形的三个点 - 将该点设置为受两个数据点的权重和距离影响的颜色（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 这里使用了+，但你需要确定适合你的应用程序的“平均”算法。

- Adam Davis

---

嗯，这取决于你想要的结果是什么。理想情况下，你会考虑到每个像素宇宙中的每个数据点，但这需要大量的处理，并可能不是你真正想要的。然而，这可能是你需要的（例如磁场）。 - Adam Davis

注意不要为输出地图中的每个像素迭代所有n个点。这是O(n*m)的，对于一个1000x1000的图像和1000个数据点，这将需要十亿次操作。这样做无法扩展。使用Shepherd算法来切割平面或类似的方法。 - Jared Updike

Answer 6

我以前在Winamp AVS中实现了类似的东西。它使用“元球”类型的方法计算每个数据点的倒二次距离（为了提高速度而避免sqrt），将其限制（例如为1.0），并对2D网格上的每个点计算这些距离之和。这将提供一个平滑变化的颜色/高度图。

如果您想查看代码，可以在我的 J10 AVS包中找到“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个点并不多，但请记住，当它是新的时候，我试图让它在320x200网格上实时运行在一台400MHz的机器上（它以约20fps的速度运行）。 :)

用red = d;等替换red =，green =和blue = ...行，以了解我的意思。所有美观性都消失了，你只剩下一个平滑变化的灰度图像，围绕数据点。

另一个编辑：我忘了说“s”是所有点的共享权重，为每个点更改它会给每个点单独的权重，例如d1 = 2 / (...)和d2 = 1 / (...)会使d1在其中心处的高度比d2高两倍。您还可以使用类似于d1 = 2 / max（...，1.0）的表达式将表达式底部的峰值平滑掉，以便它们不会在中间无限制地达到峰值。 :)

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

- jheriko

Answer 7

表面插值似乎是一个困难且数学问题。

另一种更便宜的方法是：

对于每个像素：
对于每个点：
像素.addWeight(weight(point, pixel))

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

示例权重函数：

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

这是一种相当暴力的方法，但很简单。

- rjprins

Answer 8

你正在寻找Blender所称的 "metaballs" (带有链接的维基百科文章, 例子)。从这个角度来看：

你的物体是从地面上伸出的锥形。它们都是抛物线，重量决定了它们从地面上伸出多远。或者，让它们都有相同的高度，并相应地调整抛物线的“平坦度”，大重量使锥形体很宽，而小重量则使它尖锐，甚至在一定程度上两者都可以。

我建议你实现这个并看看它的效果。

接下来，你需要将一块布或橡胶板悬挂在结果上。布会被拉伸一定量，通常由于重力而下垂。锥形体使其保持上升。

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

- Aaron Digulla