Skip to content

Commit

Permalink
ENH: put output STL file on z=0 instead of negative z values (#12)
Browse files Browse the repository at this point in the history
  • Loading branch information
@fgebhart
fgebhart authored Feb 14, 2022
1 parent e34aefe commit ada21a1
Show file tree
Hide file tree
Showing 2 changed files with 75 additions and 54 deletions.
128 changes: 74 additions & 54 deletions mapa/geometry.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -39,7 +39,14 @@ def _create_raster(array: npt.ArrayLike, max_x: int, max_y: int) -> np.ndarray:
@timing
@nb.njit(fastmath=True, cache=True)
def _compute_triangles_of_3d_surface(
raster: npt.ArrayLike, array: npt.ArrayLike, max_x: int, max_y: int, x_scale: float, y_scale: float, z_scale: float
raster: npt.ArrayLike,
array: npt.ArrayLike,
max_x: int,
max_y: int,
x_scale: float,
y_scale: float,
z_scale: float,
z_offset: float,
) -> np.ndarray:
triangles = np.full((max_x, max_y, 4, 3, 3), -1.0)
for ix in range(0, max_x):
Expand All @@ -51,64 +58,64 @@ def _compute_triangles_of_3d_surface(
# first vertex
triangles[ix, iy, 0, 0, 0] = (ix + 1 / 2) * x_scale
triangles[ix, iy, 0, 0, 1] = (iy + 1 / 2) * y_scale
triangles[ix, iy, 0, 0, 2] = (array[ix, iy]) * z_scale
triangles[ix, iy, 0, 0, 2] = (array[ix, iy]) * z_scale + z_offset
# second vertex
triangles[ix, iy, 0, 1, 0] = ix * x_scale
triangles[ix, iy, 0, 1, 1] = iy * y_scale
triangles[ix, iy, 0, 1, 2] = (raster[ix, iy]) * z_scale
triangles[ix, iy, 0, 1, 2] = (raster[ix, iy]) * z_scale + z_offset
# third vertex
triangles[ix, iy, 0, 2, 0] = (ix + 1) * x_scale
triangles[ix, iy, 0, 2, 1] = iy * y_scale
triangles[ix, iy, 0, 2, 2] = (raster[ix + 1, iy]) * z_scale
triangles[ix, iy, 0, 2, 2] = (raster[ix + 1, iy]) * z_scale + z_offset

# left triangle
# first vertex
triangles[ix, iy, 1, 0, 0] = ix * x_scale
triangles[ix, iy, 1, 0, 1] = (iy + 1) * y_scale
triangles[ix, iy, 1, 0, 2] = (raster[ix, iy + 1]) * z_scale
triangles[ix, iy, 1, 0, 2] = (raster[ix, iy + 1]) * z_scale + z_offset
# second vertex
triangles[ix, iy, 1, 1, 0] = ix * x_scale
triangles[ix, iy, 1, 1, 1] = iy * y_scale
triangles[ix, iy, 1, 1, 2] = (raster[ix, iy]) * z_scale
triangles[ix, iy, 1, 1, 2] = (raster[ix, iy]) * z_scale + z_offset
# third vertex
triangles[ix, iy, 1, 2, 0] = (ix + 1 / 2) * x_scale
triangles[ix, iy, 1, 2, 1] = (iy + 1 / 2) * y_scale
triangles[ix, iy, 1, 2, 2] = (array[ix, iy]) * z_scale
triangles[ix, iy, 1, 2, 2] = (array[ix, iy]) * z_scale + z_offset

# bottom triangle
# first vertex
triangles[ix, iy, 2, 0, 0] = (ix + 1) * x_scale
triangles[ix, iy, 2, 0, 1] = (iy + 1) * y_scale
triangles[ix, iy, 2, 0, 2] = (raster[ix + 1, iy + 1]) * z_scale
triangles[ix, iy, 2, 0, 2] = (raster[ix + 1, iy + 1]) * z_scale + z_offset
# second vertex
triangles[ix, iy, 2, 1, 0] = ix * x_scale
triangles[ix, iy, 2, 1, 1] = (iy + 1) * y_scale
triangles[ix, iy, 2, 1, 2] = (raster[ix, iy + 1]) * z_scale
triangles[ix, iy, 2, 1, 2] = (raster[ix, iy + 1]) * z_scale + z_offset
# third vertex
triangles[ix, iy, 2, 2, 0] = (ix + 1 / 2) * x_scale
triangles[ix, iy, 2, 2, 1] = (iy + 1 / 2) * y_scale
triangles[ix, iy, 2, 2, 2] = (array[ix, iy]) * z_scale
triangles[ix, iy, 2, 2, 2] = (array[ix, iy]) * z_scale + z_offset

# right triangle
# first vertex
triangles[ix, iy, 3, 0, 0] = (ix + 1 / 2) * x_scale
triangles[ix, iy, 3, 0, 1] = (iy + 1 / 2) * y_scale
triangles[ix, iy, 3, 0, 2] = (array[ix, iy]) * z_scale
triangles[ix, iy, 3, 0, 2] = (array[ix, iy]) * z_scale + z_offset
# second vertex
triangles[ix, iy, 3, 1, 0] = (ix + 1) * x_scale
triangles[ix, iy, 3, 1, 1] = iy * y_scale
triangles[ix, iy, 3, 1, 2] = (raster[ix + 1, iy]) * z_scale
triangles[ix, iy, 3, 1, 2] = (raster[ix + 1, iy]) * z_scale + z_offset
# third vertex
triangles[ix, iy, 3, 2, 0] = (ix + 1) * x_scale
triangles[ix, iy, 3, 2, 1] = (iy + 1) * y_scale
triangles[ix, iy, 3, 2, 2] = (raster[ix + 1, iy + 1]) * z_scale
triangles[ix, iy, 3, 2, 2] = (raster[ix + 1, iy + 1]) * z_scale + z_offset

return triangles.reshape((max_x * max_y * 4, 3, 3))


@timing
def _compute_triangles_of_body_side(
raster: npt.ArrayLike, max_x: int, max_y: int, x_scale: float, y_scale: float, z_scale: float, ground: float
raster: npt.ArrayLike, max_x: int, max_y: int, x_scale: float, y_scale: float, z_scale: float, z_offset: float
) -> np.ndarray:
# loop over raster and build triangles when in first and last col and row
triangles = []
Expand All @@ -121,84 +128,92 @@ def _compute_triangles_of_body_side(
if ix == 0: # first row
triangles.append( # triangle with two points at top of mesh
[
[0, iy * y_scale, raster[ix][iy] * z_scale], # first point in col
[0, (iy + 1) * y_scale, raster[ix][iy + 1] * z_scale], # second point in col
[0, iy * y_scale, ground], # first point on ground
[0, iy * y_scale, raster[ix][iy] * z_scale + z_offset], # first point in col
[0, (iy + 1) * y_scale, raster[ix][iy + 1] * z_scale + z_offset], # second point in col
[0, iy * y_scale, 0], # first point on ground
]
)
triangles.append( # triangle with two points at ground
[
[0, iy * y_scale, ground],
[0, (iy + 1) * y_scale, raster[ix][iy + 1] * z_scale],
[0, (iy + 1) * y_scale, ground],
[0, iy * y_scale, 0],
[0, (iy + 1) * y_scale, raster[ix][iy + 1] * z_scale + z_offset],
[0, (iy + 1) * y_scale, 0],
]
)
if ix == max_x - 1: # last row
triangles.append( # two points at top 3d mesh
[
[max_x * x_scale, (iy + 1) * y_scale, raster[ix][iy + 1] * z_scale], # second point in col
[max_x * x_scale, iy * y_scale, raster[ix][iy] * z_scale], # first point in col
[max_x * x_scale, iy * y_scale, ground], # first point on ground
[
max_x * x_scale,
(iy + 1) * y_scale,
raster[ix][iy + 1] * z_scale + z_offset,
], # second point in col
[max_x * x_scale, iy * y_scale, raster[ix][iy] * z_scale + z_offset], # first point in col
[max_x * x_scale, iy * y_scale, 0], # first point on ground
]
)
triangles.append( # two points at ground
[
[max_x * x_scale, (iy + 1) * y_scale, raster[ix][iy + 1] * z_scale],
[max_x * x_scale, iy * y_scale, ground],
[max_x * x_scale, (iy + 1) * y_scale, ground],
[max_x * x_scale, (iy + 1) * y_scale, raster[ix][iy + 1] * z_scale + z_offset],
[max_x * x_scale, iy * y_scale, 0],
[max_x * x_scale, (iy + 1) * y_scale, 0],
]
)
if iy == 0: # first col
# two points at top 3d mesh
triangles.append(
[
[(ix + 1) * x_scale, 0, raster[ix + 1][iy] * z_scale], # second point in col
[ix * x_scale, 0, raster[ix][iy] * z_scale], # first point in col
[ix * x_scale, 0, ground], # first point on ground
[(ix + 1) * x_scale, 0, raster[ix + 1][iy] * z_scale + z_offset], # second point in col
[ix * x_scale, 0, raster[ix][iy] * z_scale + z_offset], # first point in col
[ix * x_scale, 0, 0], # first point on ground
]
)
# two points at ground
triangles.append(
[
[(ix + 1) * x_scale, 0, raster[ix + 1][iy] * z_scale],
[ix * x_scale, 0, ground],
[(ix + 1) * x_scale, 0, ground],
[(ix + 1) * x_scale, 0, raster[ix + 1][iy] * z_scale + z_offset],
[ix * x_scale, 0, 0],
[(ix + 1) * x_scale, 0, 0],
]
)
if iy == max_y - 1: # last col
# two points at top 3d mesh
triangles.append(
[
[ix * x_scale, max_y * y_scale, raster[ix][iy] * z_scale], # first point in col
[(ix + 1) * x_scale, max_y * y_scale, raster[ix + 1][iy] * z_scale], # second point in col
[ix * x_scale, max_y * y_scale, ground], # first point on ground
[ix * x_scale, max_y * y_scale, raster[ix][iy] * z_scale + z_offset], # first point in col
[
(ix + 1) * x_scale,
max_y * y_scale,
raster[ix + 1][iy] * z_scale + z_offset,
], # second point in col
[ix * x_scale, max_y * y_scale, 0], # first point on ground
]
)
# two points at ground
triangles.append(
[
[ix * x_scale, max_y * y_scale, ground],
[(ix + 1) * x_scale, max_y * y_scale, raster[ix + 1][iy] * z_scale],
[(ix + 1) * x_scale, max_y * y_scale, ground],
[ix * x_scale, max_y * y_scale, 0],
[(ix + 1) * x_scale, max_y * y_scale, raster[ix + 1][iy] * z_scale + z_offset],
[(ix + 1) * x_scale, max_y * y_scale, 0],
]
)
return triangles


def _compute_triangles_of_bottom(max_x: int, max_y: int, x_scale: float, y_scale: float, ground: float) -> np.ndarray:
def _compute_triangles_of_bottom(max_x: int, max_y: int, x_scale: float, y_scale: float) -> np.ndarray:
triangles = []
triangles.append(
[
[0, 0, ground],
[0, max_y * y_scale, ground],
[max_x * x_scale, 0, ground],
[0, 0, 0],
[0, max_y * y_scale, 0],
[max_x * x_scale, 0, 0],
]
)
triangles.append(
[
[0, max_y * y_scale, ground],
[max_x * x_scale, max_y * y_scale, ground],
[max_x * x_scale, 0, ground],
[0, max_y * y_scale, 0],
[max_x * x_scale, max_y * y_scale, 0],
[max_x * x_scale, 0, 0],
]
)
return triangles
Expand All @@ -219,21 +234,28 @@ def compute_all_triangles(
click.echo(f"{'🗺 creating base raster for tiff...':<50s}", nl=False)
raster = _create_raster(array, max_x, max_y)

# determine height from ground
# determine z_offset
if z_offset is None:
# using the natural height, i.e. islands will have an z_offset of ~0 and mountains will have a larger z_offset
ground = -raster.min()
ground = 0 if ground > 0 else ground
z_offset = raster.min()
z_offset = 0 if z_offset > 0 else z_offset
else:
# use given input offset as height to ground
ground = -z_offset
z_offset = z_offset
if z_offset < 0:
click.echo("☝️ Warning: Be careful using negative z_offsets, as it might break your 3D model.")

# compute triangles
click.echo(f"{'⛰ computing triangles of 3d surface...':<50s}", nl=False)
dem_triangles = _compute_triangles_of_3d_surface(
raster=raster, array=array, max_x=max_x, max_y=max_y, x_scale=xy_scale, y_scale=xy_scale, z_scale=z_scale
raster=raster,
array=array,
max_x=max_x,
max_y=max_y,
x_scale=xy_scale,
y_scale=xy_scale,
z_scale=z_scale,
z_offset=z_offset,
)
click.echo(f"{'📐 computing triangles of body sides...':<50s}", nl=False)
side_triangles = _compute_triangles_of_body_side(
Expand All @@ -243,11 +265,9 @@ def compute_all_triangles(
x_scale=xy_scale,
y_scale=xy_scale,
z_scale=z_scale,
ground=ground,
)
bottom_triangles = _compute_triangles_of_bottom(
max_x=max_x, max_y=max_y, x_scale=xy_scale, y_scale=xy_scale, ground=ground
z_offset=z_offset,
)
bottom_triangles = _compute_triangles_of_bottom(max_x=max_x, max_y=max_y, x_scale=xy_scale, y_scale=xy_scale)
return np.vstack((dem_triangles, side_triangles, bottom_triangles))


Expand Down
1 change: 1 addition & 0 deletions tests/test_geometry.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -52,6 +52,7 @@ def test_compute_triangles_of_3d_surface() -> None:
x_scale=x_scale,
y_scale=y_scale,
z_scale=z_scale,
z_offset=0.0,
)
expected = np.array(
[
Expand Down

0 comments on commit ada21a1

Please sign in to comment.