绘制地形图
我一直在致力于二维连续数据的可视化项目。 您可以使用它来研究二维地图上的海拔数据或温度模式。 从本质上讲，它实际上是一种将 3 维扁平化为二维加颜色的方法。 在我的特定研究领域，我实际上并没有使用地理海拔数据，但这是一个很好的比喻，所以我将在整篇文章中坚持使用它。
无论如何，此时，我已经有了一个令我非常满意的“连续颜色”渲染器：
>渐变是标准色轮，其中红色像素表示高值坐标，紫色像素表示低值。
底层数据结构使用一些非常聪明的（如果我自己这么说的话）插值算法来实现任意深度缩放地图的细节。
此时，我想绘制一些地形轮廓线（使用二次贝塞尔曲线），但我还没有找到任何好的文献来描述用于查找这些曲线的有效算法。
为了让您了解我的想法，这里有一个穷人的实现（渲染器在遇到与轮廓线相交的像素时只使用黑色 RGB 值）：
不过，这种方法存在几个问题：
图形中具有较陡斜率的区域会导致较薄（且经常破碎）拓扑线。 理想情况下，所有地形线都应该是连续的。
图形中具有较平坦坡度的区域会导致更宽的地形线（通常是整个黑色区域，尤其是在渲染区域的外周）。
因此，我正在寻找一种矢量绘制方法来获得那些漂亮、完美的 1 像素厚曲线。 该算法的基本结构必须包括以下步骤：
在我想要绘制地形线的每个离散高程处，找到一组坐标，其中该坐标处的高程非常接近（给定任意 epsilon 值）到所需的高程。
消除冗余点。 例如，如果三个点位于一条完美的直线上，那么中心点就是多余的，因为可以在不改变曲线形状的情况下消除它。 同样，对于贝塞尔曲线，通常可以通过调整相邻控制点的位置来消除某些锚点。
将剩余的点组装成一个序列，使得两点之间的每条线段都近似于高程中性轨迹，并且没有两条线段交叉路径。 每个点序列必须创建一个封闭的多边形，或者必须与渲染区域的边界框相交。
对于每个顶点，找到一对控制点，使得所得曲线相对于步骤 #2 中消除的冗余点表现出最小误差。
确保当前渲染比例下可见的所有地形特征均由适当的地形线表示。 例如，如果数据包含高度较高但直径极小的尖峰，则仍应绘制地形线。 仅当垂直要素直径小于图像的整体渲染粒度时，才应忽略垂直要素。
但即使在这些限制下，我仍然可以想到几种不同的启发式方法来查找线条：
查找渲染边界框内的高点。 从那个高点，沿着几条不同的轨迹下山。 每当遍历线穿过高程阈值时，都会将该点添加到特定于高程的存储桶中。 当遍历路径达到局部最小值时，改变路线并上坡。
沿着渲染区域的矩形边界框执行高分辨率遍历。 在每个高程阈值处（以及在坡度反转方向的拐点处），将这些点添加到特定于高程的存储桶中。 边界遍历完成后，从这些桶中的边界点开始向内追踪。
扫描整个渲染区域，以稀疏的规则间隔进行高程测量。 对于每个测量，使用它与海拔阈值的接近程度作为决定是否对其邻居进行插值测量的机制。 使用此技术可以更好地保证整个渲染区域的覆盖，但很难将结果点组装成合理的顺序来构建路径。
所以，这些是我的一些想法......
在深入实施之前，我想看看 StackOverflow 上的其他人是否有解决此类问题的经验，并且可以为准确而高效的实施提供指导。
编辑：
我对 ellisbben 提出的“渐变”建议特别感兴趣。 我的核心数据结构（忽略一些优化插值快捷方式）可以表示为一组二维高斯函数的总和，这是完全可微的。
我想我需要一个数据结构来表示三维斜率，以及一个用于计算任意点的斜率向量的函数。 在我的脑海中，我不知道该怎么做（虽然看起来应该很容易），但如果你有一个解释数学的链接，我将非常感激！
更新：
感谢 ellisbben 和 Azim 的出色贡献，我现在可以计算场中任意点的轮廓角。 绘制真正的地形线很快就会完成！
以下是更新后的渲染图，包含和不包含我一直在使用的基于贫民窟栅格的拓扑渲染器。 每个图像包含一千个随机样本点，用红点表示。 该点的轮廓角由白线表示。 在某些情况下，在给定点处无法测量斜率（基于插值的粒度），因此出现红点时没有相应的轮廓角线。
享受！
（注意：这些渲染使用与之前的渲染不同的表面形貌——因为我在制作原型时在每次迭代中随机生成数据结构——但核心渲染方法是相同的，所以我确保您明白了。）
这是一个有趣的事实：在这些渲染的右侧，您会在完美的水平和垂直角度看到一堆奇怪的轮廓线。 这些是插值过程的产物，插值过程使用插值器网格来减少执行核心渲染操作所需的计算数量（约 500%）。 所有这些奇怪的轮廓线都出现在两个插值器网格单元之间的边界上。
幸运的是，这些文物实际上并不重要。 尽管在斜率计算期间可以检测到伪影，但最终渲染器不会注意到它们，因为它以不同的位深度运行。
再次更新：
啊啊啊啊啊，作为我入睡前的最后一次放纵，这里有另一对渲染图，一个是老式的“连续颜色”风格，另一个有 20,000 个渐变样本。 在这组渲染中，我消除了点样本的红点，因为它不必要地使图像混乱。
在这里，由于插值器集合的网格结构，您可以真正看到我之前提到的那些插值工件。 我应该强调的是，这些伪影在最终轮廓渲染中将完全不可见（因为任何两个相邻插值器单元之间的幅度差异小于渲染图像的位深度）。
胃口好！！
I've been working on a visualization project for 2-dimensional continuous data. It's the kind of thing you could use to study elevation data or temperature patterns on a 2D map. At its core, it's really a way of flattening 3-dimensions into two-dimensions-plus-color. In my particular field of study, I'm not actually working with geographical elevation data, but it's a good metaphor, so I'll stick with it throughout this post.
Anyhow, at this point, I have a "continuous color" renderer that I'm very pleased with:
The gradient is the standard color-wheel, where red pixels indicate coordinates with high values, and violet pixels indicate low values.
The underlying data structure uses some very clever (if I do say so myself) interpolation algorithms to enable arbitrarily deep zooming into the details of the map.
At this point, I want to draw some topographical contour lines (using quadratic bezier curves), but I haven't been able to find any good literature describing efficient algorithms for finding those curves.
To give you an idea for what I'm thinking about, here's a poor-man's implementation (where the renderer just uses a black RGB value whenever it encounters a pixel that intersects a contour line):
There are several problems with this approach, though:
Areas of the graph with a steeper slope result in thinner (and often broken) topo lines. Ideally, all topo lines should be continuous.
Areas of the graph with a flatter slope result in wider topo lines (and often entire regions of blackness, especially at the outer perimeter of the rendering region).
So I'm looking at a vector-drawing approach for getting those nice, perfect 1-pixel-thick curves. The basic structure of the algorithm will have to include these steps:
At each discrete elevation where I want to draw a topo line, find a set of coordinates where the elevation at that coordinate is extremely close (given an arbitrary epsilon value) to the desired elevation.
Eliminate redundant points. For example, if three points are in a perfectly-straight line, then the center point is redundant, since it can be eliminated without changing the shape of the curve. Likewise, with bezier curves, it is often possible to eliminate cetain anchor points by adjusting the position of adjacent control points.
Assemble the remaining points into a sequence, such that each segment between two points approximates an elevation-neutral trajectory, and such that no two line segments ever cross paths. Each point-sequence must either create a closed polygon, or must intersect the bounding box of the rendering region.
For each vertex, find a pair of control points such that the resultant curve exhibits a minimum error, with respect to the redundant points eliminated in step #2.
Ensure that all features of the topography visible at the current rendering scale are represented by appropriate topo lines. For example, if the data contains a spike with high altitude, but with extremely small diameter, the topo lines should still be drawn. Vertical features should only be ignored if their feature diameter is smaller than the overall rendering granularity of the image.
But even under those constraints, I can still think of several different heuristics for finding the lines:
Find the high-point within the rendering bounding-box. From that high point, travel downhill along several different trajectories. Any time the traversal line crossest an elevation threshold, add that point to an elevation-specific bucket. When the traversal path reaches a local minimum, change course and travel uphill.
Perform a high-resolution traversal along the rectangular bounding-box of the rendering region. At each elevation threshold (and at inflection points, wherever the slope reverses direction), add those points to an elevation-specific bucket. After finishing the boundary traversal, start tracing inward from the boundary points in those buckets.
Scan the entire rendering region, taking an elevation measurement at a sparse regular interval. For each measurement, use it's proximity to an elevation threshold as a mechanism to decide whether or not to take an interpolated measurement of its neighbors. Using this technique would provide better guarantees of coverage across the whole rendering region, but it'd be difficult to assemble the resultant points into a sensible order for constructing paths.
So, those are some of my thoughts...
Before diving deep into an implementation, I wanted to see whether anyone else on StackOverflow has experience with this sort of problem and could provide pointers for an accurate and efficient implementation.
Edit:
I'm especially interested in the "Gradient" suggestion made by ellisbben. And my core data structure (ignoring some of the optimizing interpolation shortcuts) can be represented as the summation of a set of 2D gaussian functions, which is totally differentiable.
I suppose I'll need a data structure to represent a three-dimensional slope, and a function for calculating that slope vector for at arbitrary point. Off the top of my head, I don't know how to do that (though it seems like it ought to be easy), but if you have a link explaining the math, I'd be much obliged!
UPDATE:
Thanks to the excellent contributions by ellisbben and Azim, I can now calculate the contour angle for any arbitrary point in the field. Drawing the real topo lines will follow shortly!
Here are updated renderings, with and without the ghetto raster-based topo-renderer that I've been using. Each image includes a thousand random sample points, represented by red dots. The angle-of-contour at that point is represented by a white line. In certain cases, no slope could be measured at the given point (based on the granularity of interpolation), so the red dot occurs without a corresponding angle-of-contour line.
Enjoy!
(NOTE: These renderings use a different surface topography than the previous renderings -- since I randomly generate the data structures on each iteration, while I'm prototyping -- but the core rendering method is the same, so I'm sure you get the idea.)
Here's a fun fact: over on the right-hand-side of these renderings, you'll see a bunch of weird contour lines at perfect horizontal and vertical angles. These are artifacts of the interpolation process, which uses a grid of interpolators to reduce the number of computations (by about 500%) necessary to perform the core rendering operations. All of those weird contour lines occur on the boundary between two interpolator grid cells.
Luckily, those artifacts don't actually matter. Although the artifacts are detectable during slope calculation, the final renderer won't notice them, since it operates at a different bit depth.
UPDATE AGAIN:
Aaaaaaaand, as one final indulgence before I go to sleep, here's another pair of renderings, one in the old-school "continuous color" style, and one with 20,000 gradient samples. In this set of renderings, I've eliminated the red dot for point-samples, since it unnecessarily clutters the image.
Here, you can really see those interpolation artifacts that I referred to earlier, thanks to the grid-structure of the interpolator collection. I should emphasize that those artifacts will be completely invisible on the final contour rendering (since the difference in magnitude between any two adjacent interpolator cells is less than the bit depth of the rendered image).
Bon appetit!!
如果你对这篇内容有疑问，欢迎到本站社区发帖提问 参与讨论，获取更多帮助，或者扫码二维码加入 Web 技术交流群。
绑定邮箱获取回复消息
由于您还没有绑定你的真实邮箱，如果其他用户或者作者回复了您的评论，将不能在第一时间通知您！
发布评论
评论（8）
梯度是一个可能对您有帮助的数学运算符。
如果您可以将插值转换为可微函数，则高度的梯度将始终指向最陡上升的方向。 所有等高的曲线都垂直于该点评估的高度梯度。
您从最高点开始的想法是明智的，但如果存在多个局部最大值，则可能会错过功能。
我建议
The gradient is a mathematical operator that may help you.
If you can turn your interpolation into a differentiable function, the gradient of the height will always point in the direction of steepest ascent. All curves of equal height are perpendicular to the gradient of height evaluated at that point.
Your idea about starting from the highest point is sensible, but might miss features if there is more than one local maximum.
I'd suggest
回应您对 @erickson 的评论并回答有关计算函数梯度的问题。 您可以按如下方式进行数值微分，而不是计算 300 项函数的导数。
给定图像中的一个点 [x,y]，您可以计算梯度（最速下降方向），
其中 dx 和 dy 可以是网格中的间距。 等高线将垂直于梯度延伸。 因此，为了获得轮廓方向 c，我们可以将 g=[v,w] 乘以矩阵 A=[0 -1, 1 0] 给出
In response to your comment to @erickson and to answer the point about calculating the gradient of your function. Instead of calculating the derivatives of your 300 term function you could do a numeric differentiation as follows.
Given a point [x,y] in your image you could calculate the gradient (direction of steepest decent)
where dx and dy could be the spacing in your grid. The contour line will run perpendicular to the gradient. So, to get the contour direction, c, we can multiply g=[v,w] by matrix, A=[0 -1, 1 0] giving
或者，有 行进方块 算法似乎适合您的问题，尽管您可能想要如果使用粗网格，请平滑结果。
您要绘制的地形曲线是二维标量场的等值面。 对于 3 维等值面，有行进立方体算法。
Alternately, there is the marching squares algorithm which seems appropriate to your problem, although you may want to smooth the results if you use a coarse grid.
The topo curves you want to draw are isosurfaces of a scalar field over 2 dimensions. For isosurfaces in 3 dimensions, there is the marching cubes algorithm.
我自己也想要这样的东西，但还没有找到基于矢量的解决方案。
不过，基于栅格的解决方案并没有那么糟糕，特别是如果您的数据是基于栅格的。 如果您的数据也是基于矢量的（换句话说，您有一个表面的 3D 模型），您应该能够进行一些实际的数学运算来找到不同高度处水平面的相交曲线。
对于基于光栅的方法，我查看每对相邻像素。 如果一个位于等高线上方，一个位于等高线下方，则显然它们之间有一条等高线。 我用来消除轮廓线锯齿的技巧是将轮廓线颜色混合到两个像素中，与它们与理想化轮廓线的接近程度成正比。
也许一些例子会有所帮助。 假设当前像素的“标高”为 12 英尺，相邻像素的标高为 8 英尺，等高线每隔 10 英尺。那么，中间有一条等高线； 使用不透明度为 50% 的轮廓线颜色绘制当前像素。 另一个像素位于 11 英尺处，并且有一个邻居位于 6 英尺处。 以 80% 的不透明度为当前像素着色。
不幸的是，我手边没有代码，而且可能还有更多代码（我依稀记得也查看过对角邻居，并通过
sqrt(2) / 2进行调整）。 我希望这足以给你要点。I've wanted something like this myself, but haven't found a vector-based solution.
A raster-based solution isn't that bad, though, especially if your data is raster-based. If your data is vector-based too (in other words, you have a 3D model of your surface), you should be able to do some real math to find the intersection curves with horizontal planes at varying elevations.
For a raster-based approach, I look at each pair of neighboring pixels. If one is above a contour level, and one is below, obviously a contour line runs between them. The trick I used to anti-alias the contour line is to mix the contour line color into both pixels, proportional to their closeness to the idealized contour line.
Maybe some examples will help. Suppose that the current pixel is at an "elevation" of 12 ft, a neighbor is at an elevation of 8 ft, and contour lines are every 10 ft. Then, there is a contour line half way between; paint the current pixel with the contour line color at 50% opacity. Another pixel is at 11 feet and has a neighbor at 6 feet. Color the current pixel at 80% opacity.
Unfortunately, I don't have the code handy, and there might have been a bit more to it (I vaguely recall looking at diagonal neighbors too, and adjusting by
sqrt(2) / 2). I hope this enough to give you the gist.我发现您想要做的事情在 MATLAB 中使用轮廓函数很容易完成。 诸如对轮廓进行低密度近似之类的事情可能可以通过对轮廓进行一些相当简单的后处理来完成。
幸运的是，GNU Octave（MATLAB 的克隆版本）具有各种等高线绘图函数的实现。 您可以查看该代码以获得几乎可以肯定在数学上合理的算法和实现。 或者，您可以将处理卸载到 Octave。 查看与其他语言交互页面，看看是否会更容易。
披露：我没有太多使用 Octave，也没有实际测试过它的等高线绘图。 然而，根据我使用 MATLAB 的经验，我可以说，只要您将数据输入 MATLAB，只需几行代码，它就能满足您所需的几乎所有内容。
另外，恭喜你制作了一个非常梵高式的坡地图。
It occurred to me that what you're trying to do would be pretty easy to do in MATLAB, using the contour function. Doing things like making low-density approximations to your contours can probably be done with some fairly simple post-processing of the contours.
Fortunately, GNU Octave, a MATLAB clone, has implementations of the various contour plotting functions. You could look at that code for an algorithm and implementation that's almost certainly mathematically sound. Or, you might just be able to offload the processing to Octave. Check out the page on interfacing with other languages to see if that would be easier.
Disclosure: I haven't used Octave very much, and I haven't actually tested it's contour plotting. However, from my experience with MATLAB, I can say that it will give you almost everything you're asking for in just a few lines of code, provided you get your data into MATLAB.
Also, congratulations on making a very VanGough-esque slopefield plot.
将您渲染的内容与真实世界的地形图进行比较 - 它们对我来说看起来一模一样！ 我不会改变任何事情...
compare what you have rendered with a real-world topo map - they look identical to me! i wouldn't change a thing...
将数据写为 HGT 文件（非常简单的数字高程数据USGS 使用的格式）并使用免费开源的 gdal_contour 工具创建等高线。 这对于陆地地图来说非常有效，限制是数据点是带符号的 16 位数字，这非常适合地球的高度范围（以米为单位），但可能不足以满足您的数据，我认为这不是一个实际地形图 - 尽管您确实提到了地形图。
Write the data out as an HGT file (very simple digital elevation data format used by USGS) and use the free and open-source gdal_contour tool to create contours. That works very well for terrestrial maps, the constraint being that the data points are signed 16-bit numbers, which fits the earthly range of heights in metres very well, but may not be enough for your data, which I assume not to be a map of actual terrain - although you do mention terrain maps.
我推荐 CONREC 方法：
如果线条太锯齿，请使用较小的网格。 如果线条足够平滑并且算法花费的时间太长，请使用更大的网格。
I recommend the CONREC approach:
If the lines are too jagged, use a smaller grid. If the lines are smooth enough and the algorithm is taking too long, use a larger grid.