Heightmaps

This notebook introduces surfaces parameterized by heightmaps.

# !pip install git+https://github.com/tum-pbs/PhiFlow@3.0
from phi.flow import *

Heightmaps can be used in 2D or 3D. They encode a grid surface made up of lines or planes, respectively. The grid vertices are equally spaced and the displacement (height) of each point is specified via a 1D or 2D tensor.

The following code creates a heightmap in 2D from a list of height values.

height = wrap([.1, .02, 0, 0, 1, .95, .8, .5, 0], spatial('x'))
bounds = Box(x=2, y=1)
heightmap = geom.Heightmap(height, bounds, max_dist=.1)
plot(heightmap)

The heightmap surface separates the inside from the outside. Whether the inside is below or above the height values can be set in the Heightmap constructor.

is_inside = CenteredGrid(lambda loc: heightmap.lies_inside(loc), x=100, y=100, bounds=Box(x=2, y=2))
plot(is_inside)

/opt/hostedtoolcache/Python/3.12.7/x64/lib/python3.12/site-packages/phi/vis/_matplotlib/_matplotlib_plots.py:167: UserWarning: This figure includes Axes that are not compatible with tight_layout, so results might be incorrect.
  plt.tight_layout()  # because subplot titles can be added after figure creation

While the inside check is always exact, the distance from the surface is an approximation. For small distances, the face directly underneath a point will usually be closest but for farther points this is not always true. To account for this, each face links to one other face that will be queried for larger distances. Finding these secondary faces is performed during heightmap construction where the parameter max_dist influences which other faces are the most important to link.

Above, we set max_dist=0.1. In the following cell, we visualize the distance from the surface.

distance = CenteredGrid(lambda loc: heightmap.approximate_signed_distance(loc), x=100, y=100, bounds=Box(x=2, y=2))
plot(distance, heightmap, overlay='args')

With max_dist=0.1, the values further than 0.1 away from the cliff get inaccurate distances. Let's increase `max_dist´ to see the difference.

heightmap = geom.Heightmap(height, bounds, max_dist=1)
distance = CenteredGrid(lambda loc: heightmap.approximate_signed_distance(loc), x=100, y=100, bounds=Box(x=2, y=2))
plot(distance, heightmap, overlay='args')

Now larger distances are covered as well. This runs with the same number of computations but trades off accuracy closer to the surface.

We can additionally query the direction to the closest surface point as well as the corresponding normal vector.

Let's plot these vectors for random points inside (orange) and outside (red) of the heightmap.

points = Box(x=3, y=(-1, 1)).sample_uniform(instance(points=200))
sgn_dist, delta, normal, _, face_index = heightmap.approximate_closest_surface(points)
plot({'closest': [heightmap, PointCloud(points, .1 * math.vec_normalize(delta))],
      'normal': [heightmap, PointCloud(points, .2 * normal)]}, overlay='list', color=[0, heightmap.lies_inside(points)*2+1])

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Cell In[6], line 3
      1 points = Box(x=3, y=(-1, 1)).sample_uniform(instance(points=200))
      2 sgn_dist, delta, normal, _, face_index = heightmap.approximate_closest_surface(points)
----> 3 plot({'closest': [heightmap, PointCloud(points, .1 * math.vec_normalize(delta))],
      4       'normal': [heightmap, PointCloud(points, .2 * normal)]}, overlay='list', color=[0, heightmap.lies_inside(points)*2+1])

File /opt/hostedtoolcache/Python/3.12.7/x64/lib/python3.12/site-packages/phi/vis/_vis.py:322, in plot(lib, row_dims, col_dims, animate, overlay, title, size, same_scale, log_dims, show_color_bar, color, alpha, err, frame_time, repeat, plt_params, max_subfigures, *fields)
    320     for i, f in enumerate(fields):
    321         idx = indices[pos][i]
--> 322         plots.plot(f, figure, axes[pos], subplots[pos], min_val, max_val, show_color_bar, color[pos][i], alpha[idx], err[idx])
    323 plots.finalize(figure)
    324 LAST_FIGURE[0] = figure

File /opt/hostedtoolcache/Python/3.12.7/x64/lib/python3.12/site-packages/phi/vis/_vis_base.py:382, in PlottingLibrary.plot(self, data, figure, subplot, space, *args, **kwargs)
    380 for recipe in self.recipes:
    381     if recipe.can_plot(data, space):
--> 382         recipe.plot(data, figure, subplot, space, *args, **kwargs)
    383         return
    384 raise NotImplementedError(f"No {self.name} recipe found for {data}. Recipes: {self.recipes}")

File /opt/hostedtoolcache/Python/3.12.7/x64/lib/python3.12/site-packages/phi/vis/_matplotlib/_matplotlib_plots.py:458, in VectorCloud2D.plot(self, data, figure, subplot, space, min_val, max_val, show_color_bar, color, alpha, err)
    456 u, v = reshaped_numpy(c_data.values.vector[dims], [vector, c_data.shape.without('vector')])
    457 color_i = color[idx]
--> 458 if (color[idx] == 'cmap').all:
    459     col = _next_line_color(subplot, kind='collections')  # ToDo
    460 elif color[idx].shape:

AttributeError: 'bool' object has no attribute 'all'

Vertical Heightmaps

The grid orientation is determined by which spatial dimension is not part of the grid. E.g. to create a vertical heightmap, the displacement values should be listed along the y direction.

height = wrap([.1, .02, 0, 0, 1, .95, .8, .5, 0], spatial('y'))
bounds = Box(x=1, y=2)
plot(geom.Heightmap(height, bounds, max_dist=.1))

2D Heightmaps

Adding a third dimension z yields a 2D heightmap in 3D space. Here, we create height values from a function, but the rest stays the same.

bounds = Box(x=2, y=2, z=1)
height = CenteredGrid(lambda pos: math.exp(-math.vec_squared(pos-1) * 3), 0, bounds['x,y'], x=10, y=10).values
heightmap = geom.Heightmap(height, bounds, max_dist=.1)
show(heightmap)

---------------------------------------------------------------------------
AssertionError                            Traceback (most recent call last)
Cell In[8], line 4
      2 height = CenteredGrid(lambda pos: math.exp(-math.vec_squared(pos-1) * 3), 0, bounds['x,y'], x=10, y=10).values
      3 heightmap = geom.Heightmap(height, bounds, max_dist=.1)
----> 4 show(heightmap)

File /opt/hostedtoolcache/Python/3.12.7/x64/lib/python3.12/site-packages/phi/vis/_vis.py:58, in show(lib, row_dims, col_dims, animate, overlay, title, size, same_scale, log_dims, show_color_bar, color, alpha, err, frame_time, repeat, plt_params, max_subfigures, *fields)
     56 kwargs = locals()
     57 del kwargs['fields']
---> 58 fig = plot(*fields, **kwargs)
     59 plots = get_plots_by_figure(fig)
     60 if isinstance(fig, Tensor):

File /opt/hostedtoolcache/Python/3.12.7/x64/lib/python3.12/site-packages/phi/vis/_vis.py:322, in plot(lib, row_dims, col_dims, animate, overlay, title, size, same_scale, log_dims, show_color_bar, color, alpha, err, frame_time, repeat, plt_params, max_subfigures, *fields)
    320     for i, f in enumerate(fields):
    321         idx = indices[pos][i]
--> 322         plots.plot(f, figure, axes[pos], subplots[pos], min_val, max_val, show_color_bar, color[pos][i], alpha[idx], err[idx])
    323 plots.finalize(figure)
    324 LAST_FIGURE[0] = figure

File /opt/hostedtoolcache/Python/3.12.7/x64/lib/python3.12/site-packages/phi/vis/_vis_base.py:382, in PlottingLibrary.plot(self, data, figure, subplot, space, *args, **kwargs)
    380 for recipe in self.recipes:
    381     if recipe.can_plot(data, space):
--> 382         recipe.plot(data, figure, subplot, space, *args, **kwargs)
    383         return
    384 raise NotImplementedError(f"No {self.name} recipe found for {data}. Recipes: {self.recipes}")

File /opt/hostedtoolcache/Python/3.12.7/x64/lib/python3.12/site-packages/phi/vis/_dash/_plotly_plots.py:494, in Scatter3D.plot(self, data, figure, subplot, space, min_val, max_val, show_color_bar, color, alpha, err)
    492     color_i = plotly_color(color[idx].native())
    493 if spatial(data.geometry):
--> 494     assert spatial(data.geometry).rank == 1
    495     xyz = math.reshaped_numpy(data[idx].points.vector[dims], [vector, data.shape.non_channel.non_spatial, spatial])
    496     xyz_padded = [[i.tolist() + [None] for i in c] for c in xyz]

AssertionError: 

Again, we can sample directions and distances from the surface.

points = bounds.sample_uniform(instance(points=200))
sgn_dist, delta, normal, _, face_index = heightmap.approximate_closest_surface(points)
show({'closest': [heightmap, PointCloud(points, .1 * math.vec_normalize(delta))],
      'normal': [heightmap, PointCloud(points, .2 * normal)]}, overlay='list', color=[0, heightmap.lies_inside(points)*2+1])

---------------------------------------------------------------------------
AssertionError                            Traceback (most recent call last)
Cell In[9], line 3
      1 points = bounds.sample_uniform(instance(points=200))
      2 sgn_dist, delta, normal, _, face_index = heightmap.approximate_closest_surface(points)
----> 3 show({'closest': [heightmap, PointCloud(points, .1 * math.vec_normalize(delta))],
      4       'normal': [heightmap, PointCloud(points, .2 * normal)]}, overlay='list', color=[0, heightmap.lies_inside(points)*2+1])

File /opt/hostedtoolcache/Python/3.12.7/x64/lib/python3.12/site-packages/phi/vis/_vis.py:58, in show(lib, row_dims, col_dims, animate, overlay, title, size, same_scale, log_dims, show_color_bar, color, alpha, err, frame_time, repeat, plt_params, max_subfigures, *fields)
     56 kwargs = locals()
     57 del kwargs['fields']
---> 58 fig = plot(*fields, **kwargs)
     59 plots = get_plots_by_figure(fig)
     60 if isinstance(fig, Tensor):

File /opt/hostedtoolcache/Python/3.12.7/x64/lib/python3.12/site-packages/phi/vis/_vis.py:322, in plot(lib, row_dims, col_dims, animate, overlay, title, size, same_scale, log_dims, show_color_bar, color, alpha, err, frame_time, repeat, plt_params, max_subfigures, *fields)
    320     for i, f in enumerate(fields):
    321         idx = indices[pos][i]
--> 322         plots.plot(f, figure, axes[pos], subplots[pos], min_val, max_val, show_color_bar, color[pos][i], alpha[idx], err[idx])
    323 plots.finalize(figure)
    324 LAST_FIGURE[0] = figure

File /opt/hostedtoolcache/Python/3.12.7/x64/lib/python3.12/site-packages/phi/vis/_vis_base.py:382, in PlottingLibrary.plot(self, data, figure, subplot, space, *args, **kwargs)
    380 for recipe in self.recipes:
    381     if recipe.can_plot(data, space):
--> 382         recipe.plot(data, figure, subplot, space, *args, **kwargs)
    383         return
    384 raise NotImplementedError(f"No {self.name} recipe found for {data}. Recipes: {self.recipes}")

File /opt/hostedtoolcache/Python/3.12.7/x64/lib/python3.12/site-packages/phi/vis/_dash/_plotly_plots.py:494, in Scatter3D.plot(self, data, figure, subplot, space, min_val, max_val, show_color_bar, color, alpha, err)
    492     color_i = plotly_color(color[idx].native())
    493 if spatial(data.geometry):
--> 494     assert spatial(data.geometry).rank == 1
    495     xyz = math.reshaped_numpy(data[idx].points.vector[dims], [vector, data.shape.non_channel.non_spatial, spatial])
    496     xyz_padded = [[i.tolist() + [None] for i in c] for c in xyz]

AssertionError: 

Combining Heightmaps

Like all Geometry objects, heightmaps can be merged using the union function. If we only pass Heightmap objects with the same grid resolution, Φ_Flow will automatically vectorize computations so that the latency does not increase.

We can use this to create complex obstacles with detailed surfaces on all sides or to create cavities.

height = (1 - math.linspace(-1, 1, spatial(x=10)) ** 2) ** .5
upper_heightmap = geom.Heightmap(height, Box(x=(-1, 1), y=None), max_dist=.1, fill_below=False, extrapolation=0)
lower_heightmap = geom.Heightmap(-height, Box(x=(-1, 1), y=None), max_dist=.1, fill_below=True, extrapolation=0)
heightmap = union(lower_heightmap, upper_heightmap)

points = Box(x=(-2, 2), y=(-2, 2)).sample_uniform(instance(points=200))
sgn_dist, delta, normal, _, face_index = heightmap.approximate_closest_surface(points)
show({'closest': [heightmap, PointCloud(points, .1 * math.vec_normalize(delta))],
      'normal': [heightmap, PointCloud(points, .2 * normal)]}, overlay='list', color=[0, heightmap.lies_inside(points)*2+1])

---------------------------------------------------------------------------
IncompatibleShapes                        Traceback (most recent call last)
Cell In[10], line 8
      6 points = Box(x=(-2, 2), y=(-2, 2)).sample_uniform(instance(points=200))
      7 sgn_dist, delta, normal, _, face_index = heightmap.approximate_closest_surface(points)
----> 8 show({'closest': [heightmap, PointCloud(points, .1 * math.vec_normalize(delta))],
      9       'normal': [heightmap, PointCloud(points, .2 * normal)]}, overlay='list', color=[0, heightmap.lies_inside(points)*2+1])

File /opt/hostedtoolcache/Python/3.12.7/x64/lib/python3.12/site-packages/phi/vis/_vis.py:58, in show(lib, row_dims, col_dims, animate, overlay, title, size, same_scale, log_dims, show_color_bar, color, alpha, err, frame_time, repeat, plt_params, max_subfigures, *fields)
     56 kwargs = locals()
     57 del kwargs['fields']
---> 58 fig = plot(*fields, **kwargs)
     59 plots = get_plots_by_figure(fig)
     60 if isinstance(fig, Tensor):

File /opt/hostedtoolcache/Python/3.12.7/x64/lib/python3.12/site-packages/phi/vis/_vis.py:252, in plot(lib, row_dims, col_dims, animate, overlay, title, size, same_scale, log_dims, show_color_bar, color, alpha, err, frame_time, repeat, plt_params, max_subfigures, *fields)
    250 ncols = uniform_bound(col_dims).volume
    251 assert nrows * ncols <= max_subfigures, f"Too many subfigures ({nrows * ncols}) for max_subfigures={max_subfigures}. If you want to plot this many subfigures, increase the limit."
--> 252 positioning, indices = layout_sub_figures(data, row_dims, col_dims, animate, overlay, 0, 0)
    253 # --- Process arguments ---
    254 plots = default_plots(positioning) if lib is None else get_plots(lib)

File /opt/hostedtoolcache/Python/3.12.7/x64/lib/python3.12/site-packages/phi/vis/_vis.py:375, in layout_sub_figures(data, row_dims, col_dims, animate, overlay, offset_row, offset_col, positioning, indices, base_index)
    373     offset += shape(e).only(row_dims).volume
    374 elif dim0.only(col_dims):
--> 375     layout_sub_figures(e, row_dims, col_dims, animate, overlay, offset_row, offset_col + offset, positioning, indices, index)
    376     offset += shape(e).only(col_dims).volume
    377 else:

File /opt/hostedtoolcache/Python/3.12.7/x64/lib/python3.12/site-packages/phi/vis/_vis.py:375, in layout_sub_figures(data, row_dims, col_dims, animate, overlay, offset_row, offset_col, positioning, indices, base_index)
    373     offset += shape(e).only(row_dims).volume
    374 elif dim0.only(col_dims):
--> 375     layout_sub_figures(e, row_dims, col_dims, animate, overlay, offset_row, offset_col + offset, positioning, indices, index)
    376     offset += shape(e).only(col_dims).volume
    377 else:

File /opt/hostedtoolcache/Python/3.12.7/x64/lib/python3.12/site-packages/phi/vis/_vis.py:362, in layout_sub_figures(data, row_dims, col_dims, animate, overlay, offset_row, offset_col, positioning, indices, base_index)
    359 if dim0.only(overlay):
    360     for overlay_index in dim0.only(overlay).meshgrid(names=True):  # overlay these fields
    361         # ToDo expand constants along rows/cols
--> 362         layout_sub_figures(data[overlay_index], row_dims, col_dims, animate, overlay, offset_row, offset_col, positioning, indices, {**base_index, **overlay_index})
    363     return positioning, indices
    364 elif dim0.only(animate):

File /opt/hostedtoolcache/Python/3.12.7/x64/lib/python3.12/site-packages/phi/vis/_vis.py:382, in layout_sub_figures(data, row_dims, col_dims, animate, overlay, offset_row, offset_col, positioning, indices, base_index)
    380 # --- data must be a plottable object ---
    381 data = to_field(data)
--> 382 overlay = data.shape.only(overlay)
    383 animate = data.shape.only(animate).without(overlay)
    384 row_shape = data.shape.only(row_dims).without(animate).without(overlay)

File /opt/hostedtoolcache/Python/3.12.7/x64/lib/python3.12/site-packages/phi/field/_field.py:211, in Field.shape(self)
    209     return batch(self._geometry) & self.resolution & non_dual(self._values).without(self.resolution) & self._geometry.shape['vector']
    210 set_shape = self._geometry.sets[self.sampled_at]
--> 211 return batch(self._geometry) & (channel(self._geometry) - 'vector') & set_shape & self._values

File /opt/hostedtoolcache/Python/3.12.7/x64/lib/python3.12/site-packages/phiml/math/_shape.py:735, in Shape.__and__(self, other)
    733 if other is dual:
    734     return concat_shapes(self, self.primal.as_dual())
--> 735 return merge_shapes(self, other)

File /opt/hostedtoolcache/Python/3.12.7/x64/lib/python3.12/site-packages/phiml/math/_shape.py:1926, in merge_shapes(order, allow_varying_sizes, *objs)
   1924 else:
   1925     if not sizes_match:
-> 1926         raise IncompatibleShapes(f"Cannot merge shapes {shapes} because dimension '{dim.name}' exists with different sizes.", *shapes)
   1927     names1 = type_group.get_item_names(dim)
   1928     names2 = sh.get_item_names(dim)

IncompatibleShapes: Cannot merge shapes [(unionⁱ=2, xˢ=9), (unionⁱ=1, xˢ=9)] because dimension 'union' exists with different sizes.

Here, we set extrapolation=0 to extend the heightmap beyond its bounds. Consequently, points with x < -1 or x > 1 count as inside the mesh and are colored red correspondingly.