Skip to content

W2 - Accessibility (Distance Fields)

In this workshop we will learn the foundations to quantitatively approach spatial accessibility. We will learn about distance fields, construct a euclidean distance field, and construct a manifold distance field.

0. Initialization

0.1 Importing the packages

import os
import topogenesis as tg
import pyvista as pv
import trimesh as tm
import numpy as np

#pv.set_jupyter_backend("ipyvtklink")

# convert mesh to pv_mesh
def tri_to_pv(tri_mesh):
    faces = np.pad(tri_mesh.faces, ((0, 0),(1,0)), 'constant', constant_values=3)
    pv_mesh = pv.PolyData(tri_mesh.vertices, faces)
    return pv_mesh

0.2 import meshes

envelope_path = os.path.relpath('../data/meshes/compulsory_envelope.obj')
context_path = os.path.relpath('../data/meshes/immediate_context.obj')

# load the mesh from file
envelope_mesh = tm.load(envelope_path)
context_mesh = tm.load(context_path)

# Check if the mesh is watertight
print(envelope_mesh.is_watertight)
print(context_mesh.is_watertight)

# initiating the plotter
p = pv.Plotter(notebook=True)

# adding the meshes
p.add_mesh(tri_to_pv(envelope_mesh), color='#abd8ff')
p.add_mesh(tri_to_pv(context_mesh), color='#aaaaaa')

#plotting
cpos = [(314.348192231551, 283.381630231551, 296.74455709155103),
 (65.08283250000001, -12.333572500000002, 21.07374465),
 (0.0, 0.0, 1.0)]
p.camera_position = cpos
p.window_size = 2000, 1000
p.show(use_ipyvtk=True)
p.screenshot("../data/images/streetdistance1")
print(p.camera_position)

0.3 Importing the Envelope Lattice

# loading the lattice from csv
lattice_path = os.path.relpath('../data/meshes/voxelized_envelope_highres.csv')
envelope_lattice = tg.lattice_from_csv(lattice_path)

# initiating the plotter
p = pv.Plotter()

# fast visualization of the lattice
envelope_lattice.fast_vis(p)

# adding the meshes
p.add_mesh(tri_to_pv(context_mesh), color='#aaaaaa')

#plotting
cpos = [(314.348192231551, 283.381630231551, 296.74455709155103),
 (65.08283250000001, -12.333572500000002, 21.07374465),
 (0.0, 0.0, 1.0)]
p.camera_position = cpos
p.window_size = 2000, 1000
p.show(use_ipyvtk=True)
p.screenshot("../data/images/streetdistance2")
print(p.camera_position)

0.4 Importing the Street Points

# import the streetnetwork as a point cloud
street_pc = tg.cloud_from_csv("../data/meshes/street_points.csv")

# initiating the plotter
p = pv.Plotter()

# fast visualization of the lattice
envelope_lattice.fast_vis(p)

# fast visualization of the point cloud
street_pc.fast_notebook_vis(p)

# adding the meshes
p.add_mesh(tri_to_pv(context_mesh), color='#aaaaaa')

#plotting
cpos = [(314.348192231551, 283.381630231551, 296.74455709155103),
 (65.08283250000001, -12.333572500000002, 21.07374465),
 (0.0, 0.0, 1.0)]
p.camera_position = cpos
p.window_size = 2000, 1000
p.show(use_ipyvtk=True)
p.screenshot("../data/images/streetdistance3")
print(p.camera_position)

1. Euclidean Distance Lattice

1.1 Distance Matrix

# extracting the centroid of all voxels
env_cens = envelope_lattice.centroids_threshold(-1)

# initializing the distance matrix
dist_m = []
# for each voxel ...
for voxel_cen in env_cens:
    # initializing the distance vector (per each voxel)
    dist_v = []
    # for each street point ...
    for street_point in street_pc:
        # find the difference vector
        diff = voxel_cen - street_point
        # raise the components to the power of two
        diff_p2 = diff**2
        # sum the components
        diff_p2s = diff_p2.sum()
        # compute the square root 
        dist = diff_p2s**0.5
        # add the distance to the distance vector
        dist_v.append(dist)
    # add the distance vector to the distance matrix
    dist_m.append(dist_v)
# change the distance matrix type, from list to array
dist_m = np.array(dist_m)


"""
st_m_shape = (env_cens.shape[0], street_pc.shape[0], 3)
st_m = np.broadcast_to(street_pc, st_m_shape)
ep_m_shape = (street_pc.shape[0], env_cens.shape[0], 3)
ep_m = np.broadcast_to(env_cens, ep_m_shape).transpose((1,0,2))
dist_m = np.linalg.norm(st_m - ep_m, axis=2)
"""

1.2 Distance to Closest Street Point

# find the distance to the closest street point for each voxel
min_dist = dist_m.min(axis=1)
# convert the minimum distance list to a lattice
street_eu_distance_lattice = tg.to_lattice(min_dist.reshape(envelope_lattice.shape), envelope_lattice)
# zero the value of the exterior voxels
envelope_eu_dist_lattice = street_eu_distance_lattice * envelope_lattice

# initiating the plotter
p = pv.Plotter()

l = envelope_eu_dist_lattice * envelope_lattice

# remapping
l = 250 * (l - l.min()) / l.max()

# Create the spatial reference
grid = pv.UniformGrid()

# Set the grid dimensions: shape because we want to inject our values
grid.dimensions = l.shape
# The bottom left corner of the data set
grid.origin = l.minbound
# These are the cell sizes along each axis
grid.spacing = l.unit

# Add the data values to the cell data
grid.point_arrays["Distance"] = l.flatten(order="F")  # Flatten the Lattice

# adding the meshes
p.add_mesh(tri_to_pv(context_mesh), opacity=0.1, style='wireframe')

# fast visualization of the point cloud
street_pc.fast_notebook_vis(p)

# adding the volume
opacity = np.array([0,0.6,0.6,0.6,0.6,0.6,0.6]) * 1.5
p.add_volume(grid, cmap="coolwarm", opacity=opacity, shade=True, show_scalar_bar=False)

#plotting
cpos = [(314.348192231551, 283.381630231551, 296.74455709155103),
 (65.08283250000001, -12.333572500000002, 21.07374465),
 (0.0, 0.0, 1.0)]
p.camera_position = cpos
p.window_size = 2000, 1000
p.show(use_ipyvtk=True)
p.screenshot("../data/images/streetdistance4")
print(p.camera_position)

2 Manifold Distance Lattice

2.1 Selecting the Closest Voxels

# selecting the closest voxels by setting a threshold 
street_connection_lattice = (0 < envelope_eu_dist_lattice) * (envelope_eu_dist_lattice < 14)
street_connection_lattice.shape

# initiating the plotter
p = pv.Plotter()

# fast visualization of the lattice
street_connection_lattice.fast_vis(p)

# adding the meshes
p.add_mesh(tri_to_pv(context_mesh), color='#aaaaaa')

#plotting
cpos = [(314.348192231551, 283.381630231551, 296.74455709155103),
 (65.08283250000001, -12.333572500000002, 21.07374465),
 (0.0, 0.0, 1.0)]
p.camera_position = cpos
p.window_size = 2000, 1000
p.show(use_ipyvtk=True)
p.screenshot("../data/images/streetdistance5")
print(p.camera_position)

2.2. The Stencil

# creating neighborhood definition
stencil = tg.create_stencil("von_neumann", 1, 1)
# setting the center to zero
stencil.set_index([0,0,0], 0)

print(stencil)

2.3 Initializing the Manifold Distance Lattice

# retrieve the neighbour list of each cell
neighs = street_connection_lattice.find_neighbours(stencil)

# set the maximum distance to sum of the size of the lattice in all dimensions.
max_dist = np.sum(street_connection_lattice.shape)

# initialize the street network distance lattice with all the street cells as 0, and all other cells as maximum distance possible
mn_dist_lattice = 1 - street_connection_lattice
mn_dist_lattice[mn_dist_lattice==1] = max_dist

# flatten the distance lattice for easy access
mn_dist_lattice_flat = mn_dist_lattice.flatten()

# flatten the envelope lattice
env_lat_flat = envelope_lattice.flatten()

2.4 Breadth-First Traversal

# main loop for breath-first traversal
for i in range(1, max_dist):
    # find the neighbours of the previous step
    next_step = neighs[mn_dist_lattice_flat == i - 1]
    # find the unique neighbours
    next_unq_step = np.unique(next_step.flatten())
    # check if the neighbours of the next step are inside the envelope
    validity_condition = env_lat_flat[next_unq_step]
    # select the valid neighbours
    next_valid_step = next_unq_step[validity_condition]
    # make a copy of the lattice to prevent overwriting in the memory
    mn_nex_dist_lattice_flat = np.copy(mn_dist_lattice_flat)
    # set the next step cells to the current distance
    mn_nex_dist_lattice_flat[next_valid_step] = i
    # find the minimum of the current distance and previous distances to avoid overwriting previous steps
    mn_dist_lattice_flat = np.minimum(mn_dist_lattice_flat, mn_nex_dist_lattice_flat)

    # check how many of the cells have not been traversed yet
    filled_check = mn_dist_lattice_flat * env_lat_flat == max_dist
    # if all the cells have been traversed, break the loop
    if filled_check.sum() == 0:
        print(i)
        break
mn_dist_lattice_flat=mn_dist_lattice_flat.astype(float)

dist_id=[]
for i in range(len(mn_dist_lattice_flat)):
    if mn_dist_lattice_flat[i] != 181:
        dist_id.append(i)
mn_dist_full=mn_dist_lattice_flat[mn_dist_lattice_flat!=181]
mn_dist_full_norm=1-(mn_dist_full-mn_dist_full.min())/(mn_dist_full.max()-mn_dist_full.min())


for i in range(len(dist_id)):
    mn_dist_lattice_flat[dist_id[i]]=mn_dist_full_norm[i]

# reshape and construct a lattice from the football field network distance list
mn_dist_lattice = mn_dist_lattice_flat.reshape(mn_dist_lattice.shape)

distance_lattice=mn_dist_lattice*envelope_lattice
# distance_lattice=(1-distance_lattice)*envelope_lattice

# # normalizing the lattice
# l=mn_dist_lattice*envelope_lattice
# l=(l-l.min())/(l.max()-l.min())
# l=(1-l)*envelope_lattice

# initiating the plotter
p = pv.Plotter(notebook=True)
l=distance_lattice
# Create the spatial reference
grid = pv.UniformGrid()

# Set the grid dimensions: shape because we want to inject our values
grid.dimensions = l.shape
# The bottom left corner of the data set
grid.origin = l.minbound
# These are the cell sizes along each axis
grid.spacing = l.unit

# Add the data values to the cell data
grid.point_arrays["Street Access"] = l.flatten(order="F")  # Flatten the Lattice

# adding the meshes
p.add_mesh(tri_to_pv(context_mesh), opacity=0.1, style='wireframe')

# fast visualization of the point cloud
street_pc.fast_notebook_vis(p)

# adding the volume
opacity = np.array([0,1,1,1,1,1,1])*2
p.add_volume(grid, cmap="coolwarm", clim=[0, 1.0],opacity=opacity, shade=True)

#plotting
cpos = [(314.348192231551, 283.381630231551, 296.74455709155103),
 (65.08283250000001, -12.333572500000002, 21.07374465),
 (0.0, 0.0, 1.0)]
p.camera_position = cpos
p.window_size = 2000, 1000
p.show(use_ipyvtk=True)
p.screenshot("../data/images/streetdistance6")
print(p.camera_position)

Credits

# save the street access latice to csv

csv_path = os.path.relpath('../data/fields/manifold_distance_field.csv')
distance_lattice.to_csv(csv_path)

__author__ = "Shervin Azadi"
__license__ = "MIT"
__version__ = "1.0"
__url__ = "https://github.com/shervinazadi/earthy_workshops"
__summary__ = "Earthy Design Studio"