kylebarron/quantized-mesh-encoder

A fast Python Quantized Mesh encoder

quantized-mesh-encoder

A fast Python Quantized Mesh encoder. Encodes a mesh with 100k coordinates and 180k triangles in 20ms. Example viewer.

The Grand Canyon and Walhalla Plateau. The mesh is created using pydelatin or pymartini, encoded using quantized-mesh-encoder, served on-demand using dem-tiler, and rendered with deck.gl.

Overview

Quantized Mesh is a format to encode terrain meshes for efficient client-side terrain rendering. Such files are supported in Cesium and deck.gl.

This library is designed to support performant server-side on-demand terrain mesh generation.

Install

With pip:

pip install quantized-mesh-encoder

or with Conda:

conda install -c conda-forge quantized-mesh-encoder

Using

API

quantized_mesh_encoder.encode

Arguments:

• f: a writable file-like object in which to write encoded bytes
• positions: (array[float]): either a 1D Numpy array or a 2D Numpy array of shape (-1, 3) containing 3D positions.
• indices (array[int]): either a 1D Numpy array or a 2D Numpy array of shape (-1, 3) indicating triples of coordinates from positions to make triangles. For example, if the first three values of indices are 0, 1, 2, then that defines a triangle formed by the first 9 values in positions, three for the first vertex (index 0), three for the second vertex, and three for the third vertex.

Keyword arguments:

• bounds (List[float], optional): a list of bounds, [minx, miny, maxx, maxy]. By default, inferred as the minimum and maximum values of positions.
• sphere_method (str, optional): As part of the header information when encoding Quantized Mesh, it's necessary to compute a bounding sphere, which contains all positions of the mesh. sphere_method designates the algorithm to use for creating the bounding sphere. Must be one of 'bounding_box', 'naive', 'ritter' or None. Default is None.
  • 'bounding_box': Finds the bounding box of all positions, then defines the center of the sphere as the center of the bounding box, and defines the radius as the distance back to the corner. This method produces the largest bounding sphere, but is the fastest: roughly 70 µs on my computer.
  • 'naive': Finds the bounding box of all positions, then defines the center of the sphere as the center of the bounding box. It then checks the distance to every other point and defines the radius as the maximum of these distances. This method will produce a slightly smaller bounding sphere than the bounding_box method when points are not in the 3D corners. This is the next fastest at roughly 160 µs on my computer.
  • 'ritter': Implements the Ritter Method for bounding spheres. It first finds the center of the longest span, then checks every point for containment, enlarging the sphere if necessary. This can produce smaller bounding spheres than the naive method, but it does not always, so often both are run, see next option. This is the slowest method, at roughly 300 µs on my computer.
  • None: Runs both the naive and the ritter methods, then returns the smaller of the two. Since this runs both algorithms, it takes around 500 µs on my computer
• ellipsoid (quantized_mesh_encoder.Ellipsoid, optional): ellipsoid defined by its semi-major a and semi-minor b axes. Default: WGS84 ellipsoid.
• extensions: list of extensions to encode in quantized mesh object. These must be Extension instances. See Quantized Mesh Extensions.

quantized_mesh_encoder.Ellipsoid

Ellipsoid used for mesh calculations.

Arguments:

• a (float): semi-major axis
• b (float): semi-minor axis

quantized_mesh_encoder.WGS84

Default WGS84 ellipsoid. Has a semi-major axis a of 6378137.0 meters and semi-minor axis b of 6356752.3142451793 meters.

Quantized Mesh Extensions

There are a variety of extensions to the Quantized Mesh spec.

quantized_mesh_encoder.VertexNormalsExtension

Implements the Terrain Lighting extension. Per-vertex normals will be generated from your mesh data.

Keyword Arguments:

• indices: mesh indices
• positions: mesh positions
• ellipsoid: instance of Ellipsoid class, default: WGS84 ellipsoid

quantized_mesh_encoder.WaterMaskExtension

Implements the Water Mask extension.

Keyword Arguments:

• data (Union[np.ndarray, np.uint8, int]): Data for water mask.

quantized_mesh_encoder.MetadataExtension

Implements the Metadata extension.

• data (Union[Dict, bytes]): Metadata data to encode. If a dictionary, json.dumps will be called to create bytes in UTF-8 encoding.

Examples

Write to file

from quantized_mesh_encoder import encode
with open('output.terrain', 'wb') as f:
    encode(f, positions, indices)

Quantized mesh files are usually saved gzipped. An easy way to create a gzipped file is to use gzip.open:

import gzip
from quantized_mesh_encoder import encode
with gzip.open('output.terrain', 'wb') as f:
    encode(f, positions, indices)

Write to buffer

It's also pretty simple to write to an in-memory buffer instead of a file

from io import BytesIO
from quantized_mesh_encoder import encode
with BytesIO() as bio:
    encode(bio, positions, indices)

Or to gzip the in-memory buffer:

import gzip
from io import BytesIO
with BytesIO() as bio:
    with gzip.open(bio, 'wb') as gzipf:
        encode(gzipf, positions, indices)

Alternate Ellipsoid

By default, the WGS84 ellipsoid is used for all calculations. An alternate ellipsoid may be useful for non-Earth planetary bodies.

from quantized_mesh_encoder import encode, Ellipsoid

# From https://ui.adsabs.harvard.edu/abs/2010EM%26P..106....1A/abstract
mars_ellipsoid = Ellipsoid(3_395_428, 3_377_678)

with open('output.terrain', 'wb') as f:
    encode(f, positions, indices, ellipsoid=mars_ellipsoid)

Quantized Mesh Extensions

from quantized_mesh_encoder import encode, VertexNormalsExtension, MetadataExtension

vertex_normals = VertexNormalsExtension(positions=positions, indices=indices)
metadata = MetadataExtension(data={'hello': 'world'})

with open('output.terrain', 'wb') as f:
    encode(f, positions, indices, extensions=(vertex_normals, metadata))

Generating the mesh

To encode a mesh into a quantized mesh file, you first need a mesh! This project was designed to be used with pydelatin or pymartini, fast elevation heightmap to terrain mesh generators.

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

png = imread(png_path)
terrain = decode_ele(png, 'terrarium')
terrain = terrain.T
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
bounds = mercantile.bounds(mercantile.Tile(x, y, z))

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

with BytesIO() as f:
    quantized_mesh_encoder.encode(f, rescaled, triangles)
    f.seek(0)
    return ("OK", "application/vnd.quantized-mesh", f.read())

You can also look at the source of _mesh() in dem-tiler for a working reference.

License

Much of this code is ported or derived from quantized-mesh-tile in some way. quantized-mesh-tile is also released under the MIT license.

quantized-mesh-encoder's Issues

Docs about using gzip as file object

Quantized mesh tiles are by gzip encoded per the spec. Show how easy it is for the file object to be a gzip file descriptor

to a buffer:

gzip.GzipFile(fileobj=io.BytesIO(), mode='w').writable()

https://stackoverflow.com/a/32794964

to a file:

with gzip.open("path/to/file.gz", 'wb') as f:
    f.write(bytes)

write required number of bytes

https://github.com/kylebarron/qmesh-py/blob/87bddb4e75db040d8134062d01f1c4498f4174f1/qmesh/encode.py#L156-L161

---

This issue was generated by todo based on a TODO comment in 87bddb4. It's been assigned to @kylebarron because they committed the code.

Probably wrong encoding on big-endian computers

I often dump the bytes of a numpy array to a buffer for fast writing to file

quantized-mesh-encoder/quantized_mesh_encoder/encode.py

Line 144 in efd52a2

f.write(u_zz.tobytes())

This probably uses system endianness. Since quantized mesh files must be little-endian, presumably files written on big-endian computers would be incorrect.

Look at the Numpy byte-swapping docs for a fix

Add tests for correctness

e.g. decode existing tile with quantized-mesh-tile, encode with mine and check?

Finish ritter method for bounding sphere

Currently I only have the "naive" method of finding an axis-aligned bounding box, then finding the sphere from that box.

It would be good to finish the ritter method, then to take the minimum of the two of them.

References:

• https://help.agi.com/AGIComponents/html/BlogBoundingSphere.htm#!
• http://geomalgorithms.com/a08-_containers.html
• http://www.realtimerendering.com/resources/GraphicsGems/gems/BoundSphere.c
• https://en.wikipedia.org/wiki/Bounding_sphere#Ritter's_bounding_sphere

Example with Trimesh

https://trimsh.org/

- change to onViewportLoad to align with Tile3DLayer

quantized-mesh-encoder/site/src/tile-layer/tile-layer.js

Lines 12 to 17 in 204bcfc

	// TODO - change to onViewportLoad to align with Tile3DLayer
	onViewportLoad: {type: 'function', optional: true, value: null, compare: false},
	onTileLoad: {type: 'function', value: tile => {}, compare: false},
	// eslint-disable-next-line
	onTileError: {type: 'function', value: err => console.error(err), compare: false},
	extent: {type: 'array', optional: true, value: null, compare: true},

---

This issue was generated by todo based on a TODO comment in 204bcfc. It's been assigned to @kylebarron because they committed the code.

this could probably be better without >>31

https://github.com/kylebarron/qmesh-py/blob/5fc5bb2a4524dd8bc8d4a83f03bec17c52b128f2/qmesh/util.py#L56-L61

---

This issue was generated by todo based on a TODO comment in 5fc5bb2. It's been assigned to @kylebarron because they committed the code.

Check triangleIndices

This is decoded with quantized-mesh-decoder in js:

The triangleIndices looks a little suspect. Especially since there are two zeros in a row.

whether to use 16 or 32 bits

https://github.com/kylebarron/qmesh-py/blob/87bddb4e75db040d8134062d01f1c4498f4174f1/qmesh/encode.py#L145-L150

---

This issue was generated by todo based on a TODO comment in 87bddb4. It's been assigned to @kylebarron because they committed the code.

Move CI testing to github actions

Normals pointing the wrong way

14/3090/6430.terrain?url=terrarium viewed in quantized-mesh-viewer

Example demo is not working properly

Hi @kylebarron,
I can't visualize proper your example demo: https://kylebarron.dev/quantized-mesh-encoder/

This is what I see:

Looks like it's related to CORS

Cross-Origin Request Blocked: The Same Origin Policy disallows reading the remote resource at 
https://us-east-1-lambda.kylebarron.dev/dem/mesh/13/1507/2820.terrain?url=terrarium&mesh_max_error=9.41&mesh_algorithm=pydelatin&flip_y=true. 
(Reason: CORS header ‘Access-Control-Allow-Origin’ missing). Status code: 530.

Vertex normals

Optionally compute vertex normals extension. There's two parts:

Normal Generation

Based on the quantized-mesh-tile implementation, essentially for every vertex you look at all the triangles that touch it. For each triangle, compute its normal, and then average the normals of the touching triangles to find each per-vertex normal.

I think that to make it performant, you'd want to first compute all the normals of the triangles. And then you can use that array of normals somehow as a lookup table, so that you can find all the triangles per vertex at once, and then average within the axis...

Normal Encoding

Quantizing and encoding each x, y, z 96-bit floating point unit vector into an x,y 16-bit representation.

References

• Spec
• Quantized mesh tile generation
• Encoding technique paper

Check winding order of input triangles

It looks like (py)martini creates a mesh with vertices (at least sometimes) in clockwise order, while the Quantized Mesh spec requires counter-clockwise winding order.

e.g. here are the first three 2D vertices, and it looks like they define a clockwise triangle

So I should add an option, probably True by default, that checks the input geometries and flips the order of indices if the triangle is in clockwise winding order.

It looks like it shouldn't actually be too slow to check, though it might be best in Cython... not sure if I can vectorize it easily

• Stackoverflow answer, a fast method for checking winding order
• explanatory blog post
• wikipedia

This might solve #9 , #11

obj to quantized-mesh

how to obj import to lod tiles without using dem?

Refactor example to use TerrainLayer subclass

It would look much cleaner to just subclass TerrainLayer, because the only thing you'd have to overwrite is loadTerrain

Allow indices to be a 2d array

e.g. pydelatin currently returns a 2d array of shape (-1, 3) (which is sensible) but encode_indices expects a 1d array.

You should also cast to np.uint32, the data type expected by the cython encode_indices implementation.

You should also use .astype(np.uint32, casting='safe') so that an error would be produced if the data can't be transformed effectively.

Water mask encoding (not generation)

I don't want to add code to generate a water mask, since you can't assume elevations < 0 are ocean. But you could add an argument to encode to allow a user to pass a pregenerated water mask.

allow 6 input elements, for min/max elevation too?

https://github.com/kylebarron/qmesh-py/blob/5fc5bb2a4524dd8bc8d4a83f03bec17c52b128f2/qmesh/encode.py#L80-L85

---

This issue was generated by todo based on a TODO comment in 5fc5bb2. It's been assigned to @kylebarron because they committed the code.

loop through points a second time, testing for containment in sphere

quantized-mesh-encoder/quantized_mesh_encoder/bounding_sphere.py

Lines 45 to 47 in 0e6b16e

	# TODO: loop through points a second time, testing for containment in sphere
	return center, radius

---

This issue was generated by todo based on a TODO comment in 0e6b16e. It's been assigned to @kylebarron because they committed the code.

Comparing deck.gl rendering vs maptiler

url:
https://api.maptiler.com/tiles/terrain-quantized-mesh/13/3090/5737.terrain?key=...

bounds:
-112.08251953125,
36.0791015625,
-112.1044921875,
36.057128906249986

x, y, level = 3090, 5737, 13

import numpy as np
self =GeographicTilingScheme()
class GeographicTilingScheme:
    def __init__(self):
        self._numberOfLevelZeroTilesX = 2
        self._numberOfLevelZeroTilesY = 1
        self._rectangle = (-np.pi, -np.pi/2, np.pi, np.pi/2)

    def getNumberOfXTilesAtLevel(self, level):
        """
        Gets the total number of tiles in the X direction at a specified level-of-detail

         @param {Number} level The level-of-detail.
         @returns {Number} The number of tiles in the X direction at the given level.

        """
        return self._numberOfLevelZeroTilesX << level

    def getNumberOfYTilesAtLevel(self, level):
        """
        Gets the total number of tiles in the Y direction at a specified level-of-detail.

         @param {Number} level The level-of-detail.
         @returns {Number} The number of tiles in the Y direction at the given level.

        """
        return self._numberOfLevelZeroTilesY << level

    def tileXYToRectangle(self, x, y, level):
        rectangle = self._rectangle

        xTiles = self.getNumberOfXTilesAtLevel(level)
        yTiles = self.getNumberOfYTilesAtLevel(level)

        width = rectangle[2] - rectangle[0]
        xTileWidth = width / xTiles
        west = x * xTileWidth + rectangle[0]
        east = (x + 1) * xTileWidth + rectangle[0]

        height = rectangle[3] - rectangle[1]
        yTileHeight = height / yTiles
        north = rectangle[3] - y * yTileHeight
        south = rectangle[3] - (y + 1) * yTileHeight

        return np.array([east, south, west, north])

    def tileXYToNativeRectangle(self, x, y, level):
        """
        Converts tile x, y coordinates and level to a rectangle expressed in the native coordinates
        of the tiling scheme.

        Args:
            - x: The integer x coordinate of the tile.
            - y: The integer y coordinate of the tile.
            - level: level The tile level-of-detail.  Zero is the least detailed.
        """
        return np.degrees(self.tileXYToRectangle(x, y, level))

from shapely.geometry import box
bounds = list(np.degrees(self.tileXYToRectangle(x, y, level)))
bounds[1] = np.abs(bounds[1])
bounds[3] = np.abs(bounds[3])

print(bounds)
g = box(bounds[0], bounds[3], bounds[2], bounds[1])
from keplergl_cli import Visualize
g.exterior.coords
Visualize(g)
type(g)
g.exterior.coords
box(*np.degrees(self.tileXYToRectangle(x, y, level)))