import copy
import errno
import json
import logging
import math
from numbers import Integral
import os
from typing import List, Union
from typeguard import typechecked
import numpy as np
import open3d as o3d
import pyproj
import scipy.spatial.transform
from matplotlib.colors import Colormap
from py3dtiles.tileset.content.b3dm import B3dm, B3dmHeader
from py3dtiles.export import BoundingBox, Node, tile_extent, Feature, BatchTable, TriangleSoup, GlTF, divide
from py3dtiles.tileset.bounding_volume_box import BoundingVolumeBox
from tqdm import tqdm
from functools import cached_property, lru_cache, partial
from.utils import _make_colors_dict
BYTELENGTH = 28
def _Y_to_Z(cube):
tmp = cube.get_axis_aligned_bounding_box()
tmp = o3d.geometry.OrientedBoundingBox.create_from_axis_aligned_bounding_box(
tmp)
R = tmp.get_rotation_matrix_from_xyz((np.pi / 2, 0, 0))
tmp.rotate(R, center=(0, 0, 0))
return [list(tmp.get_min_bound()), list(tmp.get_max_bound())]
def _rot_ecef(origin_wsg):
rot1 = scipy.spatial.transform.Rotation.from_euler(
'x', -(90 - origin_wsg[1]), degrees=True).as_matrix()
rot3 = scipy.spatial.transform.Rotation.from_euler(
'z', -(90 + origin_wsg[0]), degrees=True).as_matrix()
return rot1.dot(rot3)
def _json_to_array(jsond):
json_str = json.dumps(jsond).replace(" ", "")
return np.fromstring(json_str, dtype=np.uint8)
[docs]class VoxelTiles:
"""
3d B3dm tiling from numpy array
"""
def __init__(self, xyz: np.ndarray, dxdydz: np.ndarray, values: dict, result: int, cmap: Union[str, Colormap, type(None)] = None):
"""
xyz: array-like, shape (n, 3)
Voxel points positions.
dxdydz: array-like, shape (n, 3)
Voxels size
values: dict
dictionary holding voxel data / metadata
result: int
key of the column you want to visualize
cmap : list of str or `~matplotlib.colors.Colormap`, default: 'gist_rainbow'
A `.Colormap` instance or registered colormap name. *cmap* is only
used if *c* is an array of floats.
"""
self.xyz = xyz
self.dxdydz = dxdydz
self.values = values
# pick result to color
self.result = result
self.cmap = cmap
self.tiles_directory = None
self.origin = None
# @property
# def origin(self) -> np.ndarray:
# """
# Origin (south-west corner) of the xyz dataset.
# """
# return np.mean(self.xyz, axis=0)
@cached_property
def voxel_models(self) -> np.ndarray:
"""
Numpy.array of object containing n*2 open3d.TriangleMesh voxel and values.
"""
models = []
for pos, delta, val in zip(self.xyz, self.dxdydz, np.array(self.values['label_values']).T):
cube = o3d.geometry.TriangleMesh.create_box(*delta)
cube.translate(pos)
models.append([cube, val])
models = np.array(models, dtype=object)
self.origin = np.array([mesh.vertices for mesh in np.array(models)[:,0]]).reshape(-1,3).mean(axis=0)
for cube in models[:,0]:
cube.translate(-self.origin)
rotation = cube.get_rotation_matrix_from_xyz((-np.pi / 2, 0, 0))
cube.rotate(rotation, center=(0, 0, 0))
return models
@property
def enu_ecf_transform(self)-> np.ndarray:
project_wsg = pyproj.Transformer.from_crs('28992', '4326', always_xy=True)
project_cesium = pyproj.Transformer.from_crs('28992', '4978', always_xy=True)
origin_cesium = np.array(project_cesium.transform(*self.origin))
origin_wsg = np.array(project_wsg.transform(*self.origin))
rot1 = scipy.spatial.transform.Rotation.from_euler('x', -(90 - origin_wsg[1]),
degrees=True).as_matrix()
rot3 = scipy.spatial.transform.Rotation.from_euler('z', -(90 + origin_wsg[0]),
degrees=True).as_matrix()
rotMatrix = rot1.dot(rot3)
transform = np.column_stack([np.vstack([rotMatrix, origin_cesium]), [0, 0, 0, 1]])
return transform
[docs] def write_tiles(self, directory: str = '.'):
"""
Generate Cesium b3dm tiles (https://github.com/CesiumGS/3d-tiles/tree/main/specification) from voxels
Parameters
----------
directory: str
Output directory path
"""
if self.tiles_directory:
directory = self.tiles_directory
cubes = self.voxel_models[:, 0]
data = np.vstack(self.voxel_models[:, 1])
cubes_coords = []
logging.info("Loading voxels")
for c in tqdm(cubes):
coord = [list(np.array(c.vertices, dtype=np.float32)[idx])
for idx in np.array(c.triangles)]
cubes_coords.append(coord)
geoms = []
ids = []
logging.info("Creating binary triangles")
for i, coord in tqdm(enumerate(cubes_coords)):
triangles_array = []
triangles_array.append(np.array(coord))
ts = TriangleSoup()
ts.triangles = triangles_array
# Create a bounding volume for the TriangleSoup
min_bound = np.array(coord).reshape(-1,3).min(axis=0)
max_bound = np.array(coord).reshape(-1,3).max(axis=0)
bounding_volume = BoundingVolumeBox()
bounding_volume.set_from_mins_maxs(np.concatenate([min_bound,max_bound]))
ts.bounding_volume = bounding_volume
geoms.append(ts)
ids.append(i)
positions = [ts.get_position_array() for ts in geoms]
normals = [ts.get_normal_array() for ts in geoms]
bboxes = [ts.get_bbox() for ts in geoms]
transform = self.enu_ecf_transform.flatten()
logging.info("Creating tileset...")
max_tile_size = 2000
if 100 <= len(cubes) < 1000:
features_per_tile = 99
else:
features_per_tile = 10000
z_up_bboxes = [_Y_to_Z(c) for c in cubes]
z_up_bboxes_flat = np.vstack(z_up_bboxes)
x_min = z_up_bboxes_flat[:, 0].min()
y_min = z_up_bboxes_flat[:, 1].min()
x_max = z_up_bboxes_flat[:, 0].max()
y_max = z_up_bboxes_flat[:, 1].max()
extent = BoundingBox([x_min, y_min], [x_max, y_max])
extent_x = x_max - x_min
extent_y = y_max - y_min
# Create quadtree
tree = Node()
for i in range(0, int(math.ceil(extent_x / max_tile_size))):
for j in range(0, int(math.ceil(extent_y / max_tile_size))):
tile = tile_extent(extent, max_tile_size, i, j)
geos = []
for idx, box in enumerate(z_up_bboxes):
bbox = BoundingBox(box[0], box[1])
if tile.inside(bbox.center()):
geos.append(Feature(idx, bbox))
if len(geos) == 0:
continue
if len(geos) > features_per_tile:
node = Node(geos[0:features_per_tile])
tree.add(node)
divide(tile, geos[features_per_tile:len(geos)], i * 2,
j * 2, max_tile_size / 2., features_per_tile, node)
else:
node = Node(geos)
tree.add(node)
# Export b3dm & tileset
logging.info("Creating tiles...")
nodes = tree.all_nodes()
identity = np.identity(4).flatten('F')
try:
os.makedirs(os.path.realpath(os.path.join(directory, "tiles")))
except OSError as e:
if e.errno != errno.EEXIST:
raise
for node in nodes:
if len(node.features) != 0:
bin_arrays = []
gids = []
gdata = []
for feature in node.features:
pos = feature.index
bin_arrays.append({
'position': positions[pos],
'normal': normals[pos],
'bbox': np.array((bboxes[pos])).tolist(),
})
if ids is not None:
gids.append(ids[pos])
if data is not None:
gdata.append(data[pos])
gdata = np.vstack(gdata)
gltf = GlTF.from_binary_arrays(bin_arrays, identity)
bt = None
if ids is not None:
bt = BatchTable()
bt.add_property_as_json("id", gids)
for n, label in enumerate(self.values['label_names']):
bt.add_property_as_json(label, gdata[:, n].tolist())
b3dm = B3dm.from_gltf(gltf, bt)
# make batch_length json
batch_length = {'BATCH_LENGTH': len(
b3dm.body.batch_table.header.data['id'])}
ft_array = _json_to_array(batch_length)
# new header
glTF_arr = B3dm.from_gltf(b3dm.body.gltf, b3dm.body.batch_table).to_array()
bth_arr = b3dm.body.batch_table.to_array()
b3dm.header.tile_byte_length = BYTELENGTH + \
len(ft_array) + len(glTF_arr) + len(bth_arr)
b3dm.header.ft_json_byte_length = len(ft_array)
# new header
new_header_arr = b3dm.header.to_array()
new_body_array = np.concatenate(
[ft_array, b3dm.body.to_array()])
b3dm = np.concatenate((new_header_arr, new_body_array))
filename = os.path.realpath(os.path.join(
directory, "tiles/{0}.b3dm".format(node.id)))
with open(filename, 'wb') as f:
f.write(b3dm)
filename = os.path.realpath(os.path.join(directory, "tileset.json"))
tileset = tree.to_tileset(transform)
with open(filename, 'w') as f:
f.write(json.dumps(tileset))
figFileName = os.path.realpath(
os.path.join(directory, "colorbar.json"))
colors_dict = _make_colors_dict(self.cmap, self.values['label_names'][self.result], min(
self.values['label_values'][self.result]), max(self.values['label_values'][self.result]))
colorbar = {}
colorbar["label"] = self.values['label_names'][self.result]
colorbar["stops"] = colors_dict
with open(figFileName, "w") as outfile:
json.dump(colorbar, outfile)