PyMartini

一个用于快速RTIN地形网格生成的Cython端口 Martini，比Node中的Martini快2-3倍。唯一的依赖项是Numpy。

大峡谷的线框渲染图。网格是用 PyMartini 创建的，使用 quantized-mesh-encoder 编码，使用 dem-tiler 按需提供，并使用 deck.gl 渲染。

安装

使用pip

pip install pymartini

或使用Conda

conda install -c conda-forge pymartini

使用

示例

API是模仿Martini设计的。

from pymartini import Martini

# set up mesh generator for a certain 2^k+1 grid size
# Usually either 257 or 513
martini = Martini(257)

# generate RTIN hierarchy from terrain data (an array of size^2 length)
tile = martini.create_tile(terrain)

# get a mesh (vertices and triangles indices) for a 10m error
vertices, triangles = tile.get_mesh(10)

API

《Martini》类和 `create_tile` 和 `get_mesh` 方法是直接从JS Martini库移植的。

此外，我还包括两个辅助函数：`decode_ele` 将Mapbox Terrain RGB或Terrarium PNG数组解码为高程；以及 `rescale_positions`，它为每个顶点添加高程，并可选择将每个顶点的XY坐标线性缩放到一个新的边界框。

Martini

一个类，用于实例化 `create_tile` 和 `get_mesh` 步骤所需的常量。如以下基准测试所示，实例化 `Martini` 类是三个函数中最慢的。如果您计划创建许多相同大小的网格，请创建一个 `Martini` 类，并从中创建许多瓦片。

参数

• grid_size (int, 默认 257)：在生成网格时使用的网格大小。必须是 2^k+1。如果您的源高程图是 256x256 像素，请使用 grid_size=257 并填充边界像素。

返回

返回一个 Martini 实例，您可以在其上调用 create_tile。

Martini.create_tile

从地形数据生成 RTIN 层级。这比创建 Martini 实例要快，但比创建具有给定最大误差的网格要慢。如果您需要为同一瓦片创建具有不同误差的多个网格，则应重用 Tile 实例。

参数

• terrain (numpy ndarray)：一个表示输入高度图的 dtype float32 数组。数组可以是展平的，形状为 (2^k+1 * 2^k+1)，或者是一个形状为 (2^k+1, 2^k+1) 的二维数组。注意，对于二维数组，pymartini 预期索引顺序为 (列，行)，因此您可能需要先转置数组。目前，如果输入数组的 dtype 不是 np.float32，则会生成错误。

返回

返回一个您可以调用 get_mesh 的 Tile 实例。

Tile.get_mesh

获取具有给定最大误差的网格。

参数

• max_error (float，默认 0)：输出网格中每个三角形的最大垂直误差。例如，如果输入高度图的单位是米，使用 max_error=5 意味着网格会不断细化，直到每个三角形在 5 米范围内逼近高度图表面。

返回

返回一个包含 (vertices，triangles) 的元组。

每个都是一个平坦的 numpy 数组。Vertices 表示每个顶点的交错 2D 坐标，例如 [x0, y0, x1, y1, ...]。如果您需要 3D 坐标，可以使用下面描述的 rescale_positions 辅助函数。

triangles 表示 vertices 数组内的 索引。因此 [0, 1, 3, ...] 将使用 vertices 数组中的第一个、第二个和第四个顶点作为一个单独的三角形。

decode_ele

一个辅助函数，用于将 PNG 地形瓦片解码为高程。

参数

• png (np.ndarray)：一个三维 ndarray，包含以红色、绿色和蓝色通道编码的高程。必须具有形状 (tile_size, tile_size, >=3) 或 (>=3, tile_size, tile_size)，其中 tile_size 通常为 256 或 512。
• encoding (str)：'mapbox' 或 'terrarium' 中的一个，这是两种主要的 RGB 高程值编码。
• backfill (bool，默认 True)：是否创建一个大小为 (tile_size + 1, tile_size + 1) 的数组，填充底部和右侧边缘。这是由于 Martini 需要大小为 2^n + 1 的网格。

返回

• (np.ndarray) 包含解码后的高程值的数组。如果 backfill 为 True，则返回形状为 (tile_size + 1, tile_size + 1)，否则返回形状为 (tile_size, tile_size)，其中 tile_size 是输入数组的形状。

示例

from imageio import imread
from pymartini import decode_ele

path = './test/data/fuji.png'
fuji = imread(path)
terrain = decode_ele(fuji, 'mapbox')

rescale_positions

一个辅助函数，用于缩放 vertices 输出并添加高程。输出是一个形式为 [[x1, y1, z1], [x2, y2, z2], ...] 的 numpy ndarray。

参数

• vertices：(np.array) 从 Martini 输出的顶点
• terrain：(np.ndarray) 2d 高程数组，作为 decode_ele 输出的高程。期望具有形状 (grid_size, grid_size)。**terrain 预期是传递给 Martini.create_tile 的完全相同的数组**。如果您使用不同的或转置的数组，网格看起来会很奇怪。请参阅 #15。如果您需要转置数组，请在传递给 Martini.create_tile 之前进行操作。
• bounds：(List[float]，默认 None) 线性缩放位置值到这个范围，期望为 [minx, miny, maxx, maxy]。如果不提供，则不执行缩放
• flip_y：(bool，默认 False) 翻转 y 坐标。当原始数据源是 PNG 时可能很有用，因为 PNG 的原点位于右上角。

示例

from imageio import imread
from pymartini import decode_ele, Martini, rescale_positions

path = './test/data/terrarium.png'
png = imread(path)
terrain = decode_ele(png, 'mapbox')
martini = Martini(png.shape[0] + 1)
tile = martini.create_tile(terrain)
vertices, triangles = tile.get_mesh(10)

# Use mercantile to find the bounds in WGS84 of this tile
import mercantile
bounds = mercantile.bounds(mercantile.Tile(385, 803, 11))

# Rescale positions to WGS84
rescaled = rescale_positions(
    vertices,
    terrain,
    bounds=bounds,
    flip_y=True
    column_row=True
)

马提尼还是德尔塔林？

两种流行的地形网格生成算法是："马提尼"算法，可在JavaScript martini 库和这个Python pymartini 库中找到，以及 "德尔塔林"算法，可在C++ hmm 库、Python pydelatin 库以及JavaScript delatin 库中找到。

应该使用哪一个？

对于大多数用途，优先使用pydelatin而不是pymartini。一个很好的分析来自一个马提尼问题

Martini

• 仅适用于2^n+1 x 2^n+1的正方形网格。
• 生成网格的层次结构（单次运行后选择任意细节）
• 优化了网格速度而非质量。

德尔塔林

• 适用于任意栅格网格。
• 为特定细节生成单个网格。
• 优化了质量（给定误差下尽可能少的三角形）。

正确性

pymartini通过了原始马提尼JS库中包含的（唯一）测试用例。我还编写了一些额外的兼容性测试，以比较pymartini和马提尼的输出。我在第二步结束时发现了浮点值的一些小差异。

第二步，martini.create_tile(terrain)，计算每个可能三角形的最大误差并将它们累加。因此，小的浮点误差似乎在误差累加到较大的三角形中时被放大。这些误差似乎在JS输出的1e-5之内。

当使用512px瓦片与256px瓦片相比，这些差异更大，这加强了我的假设，即差异与Python和JavaScript之间的小的低级浮点或位运算差异有关。

如果您想更深入地了解这一点，请查看martini.pyx中的Tile.update()以及相应的马提尼代码。

类型检查

截至pymartini 0.4.0，提供了类型，可用于与mypy之类的检查器一起使用。如果您想获得最大好处，请确保启用Numpy的mypy插件。

基准测试

准备步骤在Python中比在Node中快3倍；生成网格在Python中比在Node中快2倍。

Python

git clone https://github.com/kylebarron/pymartini
cd pymartini
pip install '.[test]'
python bench.py

init tileset: 14.860ms
create tile: 5.862ms
mesh (max_error=30): 1.010ms
vertices: 9700.0, triangles: 19078.0
mesh 0: 18.350ms
mesh 1: 17.581ms
mesh 2: 15.245ms
mesh 3: 13.853ms
mesh 4: 11.284ms
mesh 5: 12.360ms
mesh 6: 8.293ms
mesh 7: 8.342ms
mesh 8: 7.166ms
mesh 9: 5.678ms
mesh 10: 5.886ms
mesh 11: 5.092ms
mesh 12: 3.732ms
mesh 13: 3.420ms
mesh 14: 3.524ms
mesh 15: 3.101ms
mesh 16: 2.892ms
mesh 17: 2.358ms
mesh 18: 2.250ms
mesh 19: 2.293ms
mesh 20: 2.281ms
20 meshes total: 155.559ms

JS (Node)

git clone https://github.com/mapbox/martini
cd martini
npm install
node -r esm bench.js

init tileset: 54.293ms
create tile: 17.307ms
mesh: 6.230ms
vertices: 9704, triangles: 19086
mesh 0: 43.181ms
mesh 1: 33.102ms
mesh 2: 30.735ms
mesh 3: 25.935ms
mesh 4: 20.643ms
mesh 5: 17.511ms
mesh 6: 15.066ms
mesh 7: 13.334ms
mesh 8: 11.180ms
mesh 9: 9.651ms
mesh 10: 9.240ms
mesh 11: 10.996ms
mesh 12: 7.520ms
mesh 13: 6.617ms
mesh 14: 5.860ms
mesh 15: 5.693ms
mesh 16: 4.907ms
mesh 17: 4.469ms
mesh 18: 4.267ms
mesh 19: 4.267ms
mesh 20: 3.619ms
20 meshes total: 290.256ms

许可协议

这个库是从Mapbox的马提尼移植过来的，该库采用ISC许可。我的添加采用MIT许可。

ISC许可

版权所有（c）2019，Mapbox

允许免费或付费使用、复制、修改和/或分发此软件的目的是任何目的，前提是上述版权声明和本许可声明出现在所有副本中。

软件按“原样”提供，作者放弃与此软件相关的所有保证，包括所有暗示的适销性和适用性的保证。在任何情况下，作者均不对任何特殊、直接、间接或后果性损害或任何因使用或性能此软件而产生的任何损害承担任何责任，无论此类损害是由于合同、疏忽或其他侵权行为，无论此类损害是由于使用或性能此软件而产生。