Street safety
A lot of this stuff hinges on having the right initial data. We cant run the agents at this moment. I Added the stuff here that was already in workshop 4, but like I said, it will not run.
Personally, I believe that we have to store the eventual data of every script in a csv file. This script needs a program of requirements.
Load Libraries
import os
import topogenesis as tg
import pyvista as pv
import trimesh as tm
import pandas as pd
import numpy as np
import functions
np.random.seed(0)
# 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
Define 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)
Define the environment
# loading the lattice from csv
lattice_path = os.path.relpath('../data/meshes/voxelized_envelope_highres.csv')
avail_lattice = tg.lattice_from_csv(lattice_path)
init_avail_lattice = tg.to_lattice(np.copy(avail_lattice), avail_lattice)
load program (insert the program of requirements here)
#loading the relative relations excell
# relative_rels = pd.read_excel('../data/relationships/relative_relations.xlsx', index_col=0)
# relative_rels_norm = relative_rels.div(relative_rels.sum(axis=1), axis=0)
# relative_rels_norm = relative_rels_norm.fillna(0)
#loading the relative preferences excell
relative_prefs = pd.read_excel('../data/relationships/relative_preferences.xlsx')
relative_prefs.iloc[:, 1:] = relative_prefs.iloc[:, 1:].div(relative_prefs.iloc[:, 1:].sum(axis=1), axis=0)
relative_prefs.iloc[:, 1:] = relative_prefs.iloc[:, 1:].fillna(0)
relative_prefs_norm = relative_prefs
relative_prefs_norm = relative_prefs_norm.head(3)
program_prefs = relative_prefs_norm.drop(["space_name", "facade_connection","penthouse_factor","street_sight", "noise_approval", "sunlight_access", "sky_view_factor", "football_field"], 1)
program_prefs
Load the value fields
# loading the lattice from csv
fields = {}
for f in program_prefs.columns:
lattice_path = os.path.relpath('../data/fields/' + f + '.csv')
fields[f] = tg.lattice_from_csv(lattice_path)
Initialize the agents
# initialize the occupation lattice
occ_lattice = avail_lattice * 0 - 1
# Finding the index of the available voxels in avail_lattice
avail_flat = avail_lattice.flatten()
avail_index = np.array(np.where(avail_lattice == 1)).T
# Randomly choosing three available voxels
agn_num = len(program_prefs)
np.random.seed()
select_id = np.random.choice(len(avail_index), agn_num)
agn_origins = avail_index[select_id]
# adding the origins to the agents locations
agn_locs = []
# for each agent origin ...
for a_id, a_origin in enumerate(agn_origins):
# add the origin to the list of agent locations
agn_locs.append([a_origin])
# set the origin in availablity lattice as 0 (UNavailable)
avail_lattice[tuple(a_origin)] = 0
# set the origin in occupation lattice as the agent id (a_id)
occ_lattice[tuple(a_origin)] = a_id
printing the ammount of agents
print(a_id)
show the agents
p = pv.Plotter(notebook=True)
# Set the grid dimensions: shape + 1 because we want to inject our values on the CELL data
grid = pv.UniformGrid()
grid.dimensions = np.array(occ_lattice.shape) + 1
# The bottom left corner of the data set
grid.origin = occ_lattice.minbound - occ_lattice.unit * 0.5
# These are the cell sizes along each axis
grid.spacing = occ_lattice.unit
# adding the boundingbox wireframe
p.add_mesh(grid.outline(), color="grey", label="Domain")
# adding axes
p.add_axes()
p.show_bounds(grid="back", location="back", color="#777777")
# Add the data values to the cell data
grid.cell_arrays["Agents"] = occ_lattice.flatten(order="F").astype(int) # Flatten the array!
# filtering the voxels
threshed = grid.threshold([-0.1, agn_num - 0.9])
# adding the voxels
p.add_mesh(threshed, show_edges=True, opacity=1.0, show_scalar_bar=False)
# adding the availability lattice
init_avail_lattice.fast_vis(p)
p.show(use_ipyvtk=True)
Run the simulation
# make a deep copy of occupation lattice
cur_occ_lattice = tg.to_lattice(np.copy(occ_lattice), occ_lattice)
# initialzing the list of frames
frames = [cur_occ_lattice]
# setting the time variable to 0
t = 0
n_frames = 30
# main feedback loop of the simulation (for each time step ...)
while t<n_frames:
# for each agent ...
for a_id, a_prefs in relative_prefs_norm.iterrows():
# retrieve the list of the locations of the current agent
a_locs = agn_locs[a_id]
# initialize the list of free neighbours
free_neighs = []
# for each location of the agent
for loc in a_locs:
# retrieve the list of neighbours of the agent based on the stencil
neighs = avail_lattice.find_neighbours_masked(stencil, loc = loc)
# for each neighbour ...
for n in neighs:
# compute 3D index of neighbour
neigh_3d_id = np.unravel_index(n, avail_lattice.shape)
# if the neighbour is available...
if avail_lattice[neigh_3d_id]:
# add the neighbour to the list of free neighbours
free_neighs.append(neigh_3d_id)
# check if found any free neighbour
if len(free_neighs)>0:
# convert free neighbours to a numpy array
fns = np.array(free_neighs)
# find the value of neighbours
# init the agent value array
a_eval = np.ones(len(fns))
# for each field...
for f in program_prefs.columns:
# find the raw value of free neighbours...
vals = fields[f][fns[:,0], fns[:,1], fns[:,2]]
# raise the the raw value to the power of preference weight of the agent
a_weighted_vals = vals ** a_prefs[f]
# multiply them to the previous weighted values
a_eval *= a_weighted_vals
#post-processing
functions.squareness(square_weight=0.9, free_neighs=free_neighs, a_eval=a_eval)
# select the neighbour with highest evaluation
selected_int = np.argmax(a_eval)
# find 3D integer index of selected neighbour
selected_neigh_3d_id = free_neighs[selected_int]
# find the location of the newly selected neighbour
selected_neigh_loc = np.array(selected_neigh_3d_id).flatten()
# add the newly selected neighbour location to agent locations
agn_locs[a_id].append(selected_neigh_loc)
# set the newly selected neighbour as UNavailable (0) in the availability lattice
avail_lattice[selected_neigh_3d_id] = 0
# set the newly selected neighbour as OCCUPIED by current agent
# (-1 means not-occupied so a_id)
occ_lattice[selected_neigh_3d_id] = a_id
# constructing the new lattice
new_occ_lattice = tg.to_lattice(np.copy(occ_lattice), occ_lattice)
# adding the new lattice to the list of frames
frames.append(new_occ_lattice)
# adding one to the time counter
t += 1
Visualize the result
p = pv.Plotter(notebook=True)
base_lattice = frames[0]
# Set the grid dimensions: shape + 1 because we want to inject our values on the CELL data
grid = pv.UniformGrid()
grid.dimensions = np.array(base_lattice.shape) + 1
# The bottom left corner of the data set
grid.origin = base_lattice.minbound - base_lattice.unit * 0.5
# These are the cell sizes along each axis
grid.spacing = base_lattice.unit
# adding the boundingbox wireframe
p.add_mesh(grid.outline(), color="grey", label="Domain")
# adding the availability lattice
init_avail_lattice.fast_vis(p)
# adding axes
p.add_axes()
p.show_bounds(grid="back", location="back", color="#aaaaaa")
def create_mesh(value):
f = int(value)
lattice = frames[f]
# Add the data values to the cell data
grid.cell_arrays["Agents"] = lattice.flatten(order="F").astype(int) # Flatten the array!
# filtering the voxels
threshed = grid.threshold([-0.1, agn_num - 0.9])
# adding the voxels
p.add_mesh(threshed, name='sphere', show_edges=True, opacity=1.0, show_scalar_bar=False)
return
p.add_slider_widget(create_mesh, [0, n_frames], title='Time', value=0, event_type="always", style="classic")
p.show(use_ipyvtk=True)
# inputting the street grid
street_grid = os.path.relpath("../data/meshes/grid_points_2022.csv")
street_g = tg.cloud_from_csv(street_grid)
# street_grid_full = os.path.relpath("../data/meshes/grid_points_2022_full.csv")
# street_g = tg.cloud_from_csv(street_grid_full)
# inputting the surrounding buildings meshes
context_path = os.path.relpath('../data/meshes/immediate_context.obj')
context_mesh = tm.load(context_path)
Collecting Agent voxels and creating a mesh out of them
# duplicating the ocuupied lattice from the agent growth process
occupied_lattice = tg.to_lattice(np.copy(avail_lattice), avail_lattice)
# setting the values of the occupied lattice to 0
occupied_lattice[occupied_lattice==1]=0
# retrieving the agent locations
for i in range(len(agn_locs)):
for loc in agn_locs[i]:
occupied_lattice[tuple(loc)]=1
# initiating the plotter
p = pv.Plotter()
# fast visualization of the lattice
occupied_lattice.fast_vis(p)
# plotting
p.show(use_ipyvtk=True)
# create a transformation matrix for each agent location
def transformation_matrix (arr):
return [[1,0,0,int(arr[0])],[0,1,0,int(arr[1])],[0,0,1,int(arr[2])],[0,0,0,1]]
# collecting center points of agent voxels
agn_cens = occupied_lattice.centroids_threshold(0)
# creating a box mesh and transforming it to each centroid of the agent voxels with the transformation matrix above
meshes = [tm.creation.box(extents=avail_lattice.unit,transform=transformation_matrix(agn_cens[i])) for i in range(len(agn_cens))]
agent_mesh = tm.util.concatenate(meshes)
# initiating the plotter
p = pv.Plotter(notebook=True)
# env_cens.fast_vis(p)
# adding the meshes
p.add_mesh(tri_to_pv(agent_mesh), color='#abd8ff')
# plotting
p.show()
Collecting facade voxels and their centroids
# creating neighborhood definition
stencil = tg.create_stencil("von_neumann", 1, 1)
# setting the center to zero
stencil.set_index([0,0,0], 0)
stencil.set_index([0,0,1], 0)
stencil.set_index([0,0,-1], 0)
stencil.function = tg.sfunc.sum
print(stencil)
# collecting neighbours count for occupied lattice
occupied_lat_sten=occupied_lattice.apply_stencil(stencil)
facade_lattice=avail_lattice*0
# seperating the voxels with neighbours and excluding the occupied agents to retain only facade voxels
facade_lattice[occupied_lat_sten>0]=1
facade_lattice[occupied_lattice==1]=0
# collecting center points of facade voxels
fcd_cens = facade_lattice.centroids_threshold(0)
# initiating the plotter
p = pv.Plotter()
# fast visualization of the facade lattice
facade_lattice.fast_vis(p)
# fast visualization of the facade points
fcd_cens.fast_vis(p)
p.add_mesh(tri_to_pv(agent_mesh), color='#FAE31D')
# plotting
p.show(use_ipyvtk=True)
adding context mesh to the agent mesh
# combining the surrounding building meshes with the agent meshes
context_mesh_complete=context_mesh+agent_mesh
Putting it all together in a single visualization
# initiating the plotter
p = pv.Plotter(notebook=True)
# adding the street grid points
street_g.fast_vis(p)
# facade voxels
facade_lattice.fast_vis(p)
# facade centroids
fcd_cens.fast_vis(p)
# (these aren't the centroids you're looking for)
# not_fcd_cens.fast_vis(p)
# adding the meshes
p.add_mesh(tri_to_pv(context_mesh_complete), color='#aaaaaa')
# plotting
p.show(use_ipyvtk=True)
extracting vector between facade centroid and street points
# Extracting the vector between the facade centroid and the street points.
ray_dir=[]
run_count=0
for i in range(len(fcd_cens)):
for f in range(len(street_g)):
ray_dir1=street_g[f]-fcd_cens[i]
ray_dir.append(ray_dir1)
run_count+=1
ray_dir=np.array(ray_dir)
print(run_count)
sky_view scrpit
# all the ray origins for each sky directions (vectorization format)
ray_srcs = np.tile(fcd_cens, [1, len(street_g)]).reshape(-1, 3)
# all the ray directions for each centroid
ray_dirs = np.tile(street_g, [len(fcd_cens), 1])
print(len(ray_dirs))
print(len(ray_srcs))
# computing the intersections of rays with the context mesh
tri_id, ray_id = context_mesh_complete.ray.intersects_id(ray_origins=ray_srcs, ray_directions=ray_dir, multiple_hits=False)
intersection post process
# initialize the 'hits' array with 0
hits = np.array([0] * len(ray_dirs))
# rays that hit the context mesh are set to 1 and to floats
hits[ray_id] = 1
hits=hits.astype(float)
# we collect only the rays that are unobstructed and we call them the good rays.
# we collect the ID of the rays to use in restoring them back into the right location later
good_ray=[]
good_ray_id=[]
for i in range(len(hits)):
if hits[i]==0:
good_ray.append(ray_dir[i])
good_ray_id.append(i)
good_ray=np.array(good_ray)
good_ray_id=np.array(good_ray_id)
print(good_ray)
print(good_ray_id)
Scoring the streetview rays
# we create a vector in the negatice z direction
z_n=[0,0,-1]
# we computer the sin(theta) of the fot product of the two vectors
# we collect the ray length
# we compute the contribution the length gives
sin_product=[]
ray_len=[]
length_contribution=[]
for i in range(len(good_ray)):
ray_l=np.linalg.norm(good_ray[i])
ray_len.append(ray_l)
dot_prod=np.dot(good_ray[i]/ray_l,z_n)
sin_prod=np.sqrt(1-(dot_prod)**2)
sin_product.append(sin_prod)
len_contr=np.exp((-ray_l**2)/(40**2))
length_contribution.append(len_contr)
sin_product=np.array(sin_product)
ray_len=np.array(ray_len)
length_contribution=np.array(length_contribution)
# we then provide a score to each ray
total_score=100*sin_product*length_contribution
replacing the hit information with the scoring information
# we copy the hits and turn the values to 0,
# we then replace where hits previously had 0 to the total_score of the ray and leave what used to be 1 as 0
score_hits=hits*0
for i in range(len(good_ray_id)):
score_hits[good_ray_id[i]]=total_score[i]
# reshape the 'hits' array to (len(facade centroids), len(directions))
score = score_hits.reshape(len(fcd_cens), -1)
# sum up all the scores per centroid
vox_score=np.sum(score,axis=1)
# normalize the results
vox_score_norm=(vox_score-vox_score.min())/(vox_score.max()-vox_score.min())
# flatten the facade_lattice to retrieve the 1 dimensional lattice id
facade_lattice_flat=facade_lattice.flatten()
# replacing the values of 1 with the normalized score per centroid.
facade_lattice_flat=facade_lattice_flat.astype(float)
fcd_id=np.array(np.argwhere(facade_lattice_flat==1))
for i in range(len(fcd_id)):
facade_lattice_flat[fcd_id[i]]=vox_score_norm[i]
# reshaping and converting to lattice
facade_streetview = facade_lattice_flat.reshape(facade_lattice.shape)
facade_streetview_lat= tg.to_lattice(facade_streetview, facade_lattice)
# initiating the plotter
p = pv.Plotter(notebook=True)
l = facade_streetview_lat
# 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["Sreet visibilw points"] = l.flatten(order="F") # Flatten the Lattice
p.add_mesh(tri_to_pv(agent_mesh), color='#aaaaaa')
street_g.fast_vis(p)
# adding the boundingbox wireframe
p.add_mesh(grid.outline(), color="grey", label="Domain")
# 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)
# plotting
p.show(use_ipyvtk=True)