Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

ENH: put output STL file on z=0 instead of negative z values #12

Merged
merged 1 commit into from
Feb 14, 2022
Merged

ENH: put output STL file on z=0 instead of negative z values #12

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
128 changes: 74 additions & 54 deletions mapa/geometry.py
View file
Open in desktop
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
Open in desktop
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