MCDA: Seed Allocation + Voxel Occupation + Shaft and corridor growth
In this notebook the agents will grow in their desired shape. Afterwards the shafts and corridors are generated inside of the generated lattice.
0. Initialization
0.1. Load required libraries
import os
import topogenesis as tg
import pyvista as pv
import trimesh as tm
import numpy as np
import networkx as nx
import pandas as pd
import matplotlib.pyplot as plt
from scipy.sparse import csr_matrix
from sklearn.cluster import KMeans
np.random.seed(0)
# extra import function
def lattice_from_csv(file_path):
# read metadata
meta_df = pd.read_csv(file_path, nrows=3)
shape = np.array(meta_df['shape'])
unit = np.array(meta_df['unit'])
minbound = np.array(meta_df['minbound'])
# read lattice
lattice_df = pd.read_csv(file_path, skiprows=5)
# create the buffer
buffer = np.array(lattice_df['value']).reshape(shape)
# create the lattice
l = tg.to_lattice(buffer, minbound=minbound, unit=unit)
return l
0.2.1 Define the Neighborhood (Stencil)
# creating neighborhood definition
stencil = tg.create_stencil("von_neumann", 1, 1)
# setting the center to zero
stencil.set_index([0,0,0], 0)
0.2.2 Custom neighbours function to circumvent bugs
def neighbours(loc, avail_lattice):
output = []
output_id = []
my_stencil = np.array([[1, 0, 0],
[-1, 0, 0],
[0, 1, 0],
[0, -1, 0],
[0, 0, 1],
[0, 0, -1]])
for n, i in enumerate(my_stencil):
temp = np.array(loc) + i
try:
if avail_lattice[tuple(temp)] == True and occ_lattice[tuple(temp)] == -1:
output.append(temp)
except:
pass
return output
0.3. Load the envelope lattice as the avialbility lattice
# loading the lattice from csv
lattice_path = os.path.relpath('../data/324/solar_envelope_324.csv')
avail_lattice = lattice_from_csv(lattice_path)
init_avail_lattice = tg.to_lattice(np.copy(avail_lattice), avail_lattice)
0.4. Load Agents Information
# loading program (agents information) from CSV
prgm_path = os.path.relpath('../data/324/Relations - Blad11.csv')
agn_info = np.genfromtxt(prgm_path, delimiter=',')[1:, 1:]
# extract agent ids
length = len(agn_info)
agn_ids = agn_info[:, 0]
agn_prefs = agn_info[:, 1:]
# extract agent preferences
agn_prefs_env = agn_info[:, length + 1:]
0.5. Initialize environment information layers from Sun Access Lattice and Entrance Access Lattice
# loading the lattices from csv
sun_acc_path = os.path.relpath('../data/324/solar_accessibility_324.csv')
sun_acc_lattice = lattice_from_csv(sun_acc_path)
svf_acc_path = os.path.relpath('../data/324/SVF_324.csv')
svf_acc_lattice = lattice_from_csv(svf_acc_path)
gnd_floor_pref_path = os.path.relpath('../data/324/closeness_ground_floor_324.csv')
gnd_floor_pref_lattice = lattice_from_csv(gnd_floor_pref_path)
sub_floor_pref_path = os.path.relpath('../data/324/closeness_subterrain_324.csv')
sub_floor_pref_lattice = lattice_from_csv(sub_floor_pref_path)
coffee_pref_path = os.path.relpath('../data/324/closeness_coffee_corner_324.csv')
coffee_pref_lattice = lattice_from_csv(coffee_pref_path)
vis_ent_pref_path = os.path.relpath('../data/324/closeness_visitors_entrance_324.csv')
vis_ent_pref_lattice = lattice_from_csv(vis_ent_pref_path)
facade_pref_path = os.path.relpath('../data/324/closeness_to_facade_324.csv')
facade_pref_lattice = lattice_from_csv(facade_pref_path)
park_ent_pref_path = os.path.relpath('../data/324/closeness_car_park_entrance_324.csv')
park_ent_pref_lattice = lattice_from_csv(park_ent_pref_path)
NE_pref_path = os.path.relpath('../data/324/closeness_to_NE_facade_324.csv')
NE_pref_lattice = lattice_from_csv(NE_pref_path)
NW_pref_path = os.path.relpath('../data/324/closeness_to_NW_facade_324.csv')
NW_pref_lattice = lattice_from_csv(NW_pref_path)
SE_pref_path = os.path.relpath('../data/324/closeness_to_SE_facade_324.csv')
SE_pref_lattice = lattice_from_csv(SE_pref_path)
SW_pref_path = os.path.relpath('../data/324/closeness_to_SW_facade_324.csv')
SW_pref_lattice = lattice_from_csv(SW_pref_path)
vis_ent_pref_path = os.path.relpath('../data/324/closeness_visitors_entrance_324.csv')
vis_ent_pref_lattice = lattice_from_csv(vis_ent_pref_path)
# list the environment information layers (lattices)
env_info = [sun_acc_lattice, svf_acc_lattice, gnd_floor_pref_lattice, sub_floor_pref_lattice,
coffee_pref_lattice, vis_ent_pref_lattice, park_ent_pref_lattice, NE_pref_lattice,
NW_pref_lattice, SE_pref_lattice, SW_pref_lattice, facade_pref_lattice]
# Fixing a fault in one of the lattices
facade_pref_lattice2 = np.zeros((len(facade_pref_lattice), len(facade_pref_lattice[0]), len(facade_pref_lattice[0][0]) - 2))
for i, x in enumerate(facade_pref_lattice):
for j, y in enumerate(x):
y = np.delete(y, 4)
y = np.delete(y, 0)
facade_pref_lattice2[i][j] = y
facade_pref_lattice = facade_pref_lattice2
env_info = [sun_acc_lattice, svf_acc_lattice, gnd_floor_pref_lattice, sub_floor_pref_lattice,
coffee_pref_lattice, vis_ent_pref_lattice, park_ent_pref_lattice, NE_pref_lattice,
NW_pref_lattice, SE_pref_lattice, SW_pref_lattice, facade_pref_lattice]
1. ABM Simulation
1.1. Initialize the Agents
# initialize the occupation lattice
occ_lattice = avail_lattice * 0 - 1
test_lattice = avail_lattice * 0 - 1
keep_distance_lattice = tg.to_lattice(np.copy(occ_lattice), occ_lattice)
# 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
# count the number of spaces (rows) and intiialize an agent for each space
agn_num = len(agn_info)
# adding the origins to the agents locations
agn_locs = []
# retrieving the entrance access value of the free neighbours
for a_id in agn_ids:
voxel_vals = []
pot_voxels = []
# retrieve agent preferences
a_pref = agn_prefs_env[int(a_id)]
print(agn_prefs_env, a_id)
# Voxel Evaluation Loop
for pot_vox in avail_index:
if avail_lattice[tuple(pot_vox)]: #so if the voxel is available
global_vox_value = 1.0
# for every lattice in the environment informations
for i, info_lattice in enumerate(env_info):
# Here we utilise Fuzzy Logics to be able to compare different layers
# of environmental information and evaluate the voxel for the agent.
# This method is introduced, and generalised in Pirouz Nourian dissertation:
# section 5.7.3, pp. 201-208, eq. 57. You can refer to this section for
# comprehensive mathematical details.
vox_val = info_lattice[tuple(pot_vox)]
agn_vox_val = np.power(vox_val, a_pref[i])
global_vox_value *= agn_vox_val
# add the neighbour value to the list of values
voxel_vals.append(global_vox_value)
pot_voxels.append(pot_vox)
# convert to numpy array
voxel_vals = np.array(voxel_vals)
# convert to numpy array
pot_voxels = np.array(pot_voxels)
# select the neighbour with highest value
selected_int = np.argmax(voxel_vals) #use argsort?
# find 3D intiger index of selected neighbour
selected_neigh_3d_id = tuple(pot_voxels[selected_int].T)
# 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.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
1.2.1 Initialize the distance matrices
# func for distance matrix for every agent seed, can be used in a kind of dynamic relaxation
# find the number of all voxels
vox_count = avail_lattice.size
# initialize the adjacency matrix
adj_mtrx = np.zeros((vox_count,vox_count))
# Finding the index of the available voxels in avail_lattice
avail_index = np.array(np.where(avail_lattice == 1)).T
# fill the adjacency matrix using the list of all neighbours
for vox_loc in avail_index:
# find the 1D id
vox_id = np.ravel_multi_index(vox_loc, avail_lattice.shape)
# retrieve the list of neighbours of the voxel based on the stencil
vox_neighs = avail_lattice.find_neighbours_masked(stencil, loc = vox_loc)
# iterating over the neighbours
for neigh in vox_neighs:
# setting the entry to one
adj_mtrx[vox_id, neigh] = 1.0
# construct the graph
g = nx.from_numpy_array(adj_mtrx)
# insert the distance lattices into env_info
for i in range(agn_num):
env_info.insert(0, avail_lattice * float(0))
1.2.2 Function to update the distance matrices
def update_closeness_mtrx():
global agn_locs, avail_lattice, g
for n, i in enumerate(agn_locs):
dist_mtrx = avail_lattice * float(0)
for j in g.nodes():
if np.unravel_index(j, avail_lattice.shape)[2] == i[0][2]: #Calculate distance on the horizontal axis
temp_vox = np.unravel_index(j, avail_lattice.shape)
dist_vec = np.array(i) - np.array(temp_vox)
#now we have absolute distance, let's make it normalized closeness (0 - 1)
dist_mtrx[temp_vox[0], temp_vox[1], :] = float(np.linalg.norm(dist_vec))
argmax = np.argmax(dist_mtrx)
loc3d = np.unravel_index(argmax, avail_lattice.shape)
maxdist = dist_mtrx[tuple(loc3d)]
closeness_mtrx = 1 - dist_mtrx / maxdist #normalising
env_info[n]= closeness_mtrx
update_closeness_mtrx() #first calculation after the initial seed agent placement
1.3. Defining the depth function
###This function makes sure the building doesn't exceed a maximum depth
def depth(occ_lattice, a_id, avail_lattice, agn_locs, max_depth, max_amount):
global init_avail_lattice
#Loop through the agent locations
for loc in agn_locs[a_id]:
output = 0
my_stencil = np.array([[1, 0, 0],
[-1, 0, 0],
[0, 1, 0],
[0, -1, 0]])
neigh_check = []
forbidden_voxels = []
for i in my_stencil:
temp = loc
# append the neighbours in a given distance to neigh_check
for j in range(max_depth):
temp = tuple(np.add(i, temp))
neigh_check.append(temp)
counter = 0
for n in neigh_check:
try:
#is the neighbours occupied by a voxel?
if avail_lattice[n] == False and init_avail_lattice[n] == True:
counter += 1
except:
pass
counter_forbidden = 0
if counter == len(neigh_check):
#check how many max_dist + 1 voxels aren't occupied
for t in my_stencil:
vox_loc = tuple(np.add((np.array(t) * (max_depth + 1)), loc))
try:
dist_plus_one = occ_lattice[vox_loc]
if dist_plus_one == -1:
forbidden_voxels.append(vox_loc)
except:
pass
#make those not occupied voxels unavailable if the number of these voxels is equal to the max_amount
if len(forbidden_voxels) == max_amount:
for f in forbidden_voxels:
avail_lattice[f[0], f[1], 1:] = False
1.4. Defining the keep_distance function
###This function makes sure other functions don't get closer than 'distance' to this function
def keep_distance(a_id, loc, keep_distance_lattice, distance, avail_lattice, switch):
global occ_lattice
my_stencil = np.array([[1, 0, 0],
[-1, 0, 0],
[0, 1, 0],
[0, -1, 0]])
def recurse(my_stencil, loc, n, keep_distance_lattice): #recursive function to find all voxels within a radius
if n > 0:
for i in my_stencil:
recurse(my_stencil, loc+i, n-1, keep_distance_lattice)
else:
try:
keep_distance_lattice[loc[0], loc[1], 3:] = a_id #set the keep_distance_lattice to the agent id
except:
pass
for i in range(distance):
recurse(my_stencil, loc, i, keep_distance_lattice)
2.1. 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 = 1700
threshold = 100 #for how many iterations should the seed agents walk around
forbidden = []
have_max_depth = [11, 21] # selecting functions that have a max depth
squareness_factor = 1.5
are_square = [0, 3, 4, 5, 6, 9, 10, 11, 12, 13, 14, 15, 16, 18, 21, 22, 24] # selecting specific functions to be more square
most_recent = np.array(agn_num)
# Simulation Loop
# main feedback loop of the simulation (for each time step ...)
while t<n_frames:
print(t)
# Agent Loop
# for each agent ...
for a_id in range(agn_num):
if len(agn_locs[a_id]) <= agn_prefs[a_id][-1]: #
# retrieve the list of the locations of the current agent
a_locs = agn_locs[a_id]
# initialize the list of free neighbours
free_neighs = []
free_neighs_1d = []
switch = True #use the depth function? Mainly for troubleshooting reasons
if switch == True and (a_id in have_max_depth) and t > threshold:
depth(occ_lattice, a_id, avail_lattice, agn_locs, 2, 1)
# Location loop
# for each location of the agent
for loc in a_locs:
# retrieve the list of neighbours of the agent based on the stencil
neighs = neighbours(loc, avail_lattice)
# for each neighbour ...
for n in neighs
# if the neighbour is available
if ((avail_lattice[tuple(n)] or (np.array(loc) == np.array(n)).all()) and
(keep_distance_lattice[tuple(n)] == -1 or (keep_distance_lattice[tuple(n)] == a_id)))
#if the agents are walking around
if t < threshold:
free_neighs.append(n)
free_neighs_1d.append(np.ravel_multi_index(n, avail_lattice.shape))
#else: the growing phase has started
else:
temp_down = [n[0], n[1], n[2] - 1]
#if the function cannot grow underground, the z value should be between 1 and 5 or the voxel below
#should also be occupied
if t >= threshold and not (a_id == 0 or a_id == 1 or a_id == 22) and n[2] > 1 and ((1 < n[2] < 5) or ((occ_lattice[tuple(temp_up)] == a_id) or
(occ_lattice[tuple(temp_down)] != -1))):
free_neighs.append(n)
free_neighs_1d.append(np.ravel_multi_index(n, avail_lattice.shape))
#if the function must grow underground, the z value should be 1
elif t >= threshold and (a_id == 0 or a_id == 1 or a_id == 22) and n[2] == 1:
free_neighs.append(n)
free_neighs_1d.append(np.ravel_multi_index(n, avail_lattice.shape))
# check if found any free neighbour
if len(free_neighs)>0:
# convert free neighbours to a numpy array
free_neighs = np.array(free_neighs)
free_neighs_1d = np.array(free_neighs_1d)
neigh_vals = []
# retrieve agent preferences
a_pref = agn_prefs[a_id]
# Neighbour Evaluation Loop
for neigh in free_neighs:
neigh_value = 1.0
# for every lattice in the environment informations
for i, info_lattice in enumerate(env_info):
if (a_id != i and t < threshold) or t > threshold:
# Here we utilise Fuzzy Logics to be able to compare different layers
# of environmental information and evaluate the voxel for the agent.
# This method is introduced, and generalised in Pirouz Nourian dissertation:
# section 5.7.3, pp. 201-208, eq. 57. You can refer to this section for
# comprehensive mathematical details.
vox_val = info_lattice[tuple(neigh)]
agn_vox_val = np.power(vox_val, a_pref[i])
neigh_value *= agn_vox_val
# add the neighbour value to the list of values
neigh_vals.append(neigh_value)
#squareness weights added to specific functions
if a_id in are_square:
uniq_neighs, unique_to_norm_order, unique_neighs_count = np.unique(free_neighs_1d, return_counts=True, return_inverse=True)
extra_weight = (unique_neighs_count - 1) * squareness_factor
for i in range(len(neigh_vals)):
neigh_vals[i] *= (1 + extra_weight[unique_to_norm_order][i])
# convert to numpy array
neigh_vals = np.array(neigh_vals)
# select the neighbour with highest value
selected_int = np.argmax(neigh_vals)
# find 3D intiger index of selected neighbour
selected_neigh_3d_id = tuple(free_neighs[selected_int].T)
#if the agent needs to have a maximum depth, it also should have a distance to other functions
if t >= threshold and (a_id in have_max_depth):
keep_distance(a_id, selected_neigh_3d_id, keep_distance_lattice, 3, avail_lattice, False)
# 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)
#if the agents are walking, the previous agent must be deleted and that voxel must be made available
if t < threshold and len(agn_locs[a_id]) > 1:
temp = agn_locs[a_id][0]
avail_lattice[tuple(temp)] = True
occ_lattice[tuple(temp)] = -1
del agn_locs[a_id][0]
# set the newly selected neighbour as UNavailable (0) in the availability lattice
avail_lattice[selected_neigh_3d_id] = False
#No voxels can grow above the hub
if (a_id == 2 or a_id == 18) and t >= threshold:
avail_lattice[selected_neigh_3d_id[0], selected_neigh_3d_id[1], 2:] = False
# 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
if t < threshold:
occ_lattice = avail_lattice * 0 - 1
update_closeness_mtrx()
2.2. Visualizing the simulation
pv.set_plot_theme("document")
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")
cmap = plt.cm.get_cmap("viridis", 5)
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, show_scalar_bar=False)
return
p.add_slider_widget(create_mesh, [0, n_frames], title='Time', value=0, event_type="always", style="classic", pointa=(0.1, 0.1), pointb=(0.9, 0.1))
p.window_size = 1000, 1000
# for i in range(171):
# create_mesh(i * 10)
# p.screenshot("growth_" + str(i * 10))
p.show(use_ipyvtk=True)
2.3. Saving lattice frames in CSV
lattice = frames[-1]
csv_path = os.path.relpath('../data/324/growth_output.csv')
lattice.to_csv(csv_path)
3. Shafts and Corridors
3.1. stencil creation
# creating neighborhood definition
stencil_flat = tg.create_stencil("von_neumann", 1, 1)
# setting the center to zero
stencil_flat.set_index([0,0,0], 0)
stencil_flat.set_index([0,0,1], 0)
stencil_flat.set_index([0,0,-1], 0)
print(stencil_flat)
avail_index_bool = np.array((frames[-1] >= 0).tolist()) # making a new availability boolean lattice
# ignoring some functions that do not need a shaft to be grown inside (park for instance)
avail_index_bool[np.where(frames[-1] == 2)] = False
avail_index_bool[np.where(frames[-1] == 13)] = False
avail_index_bool[np.where(frames[-1] == 14)] = False
#lattice creation
avail_lattice_bool = tg.to_lattice(avail_index_bool.reshape(avail_lattice.shape),avail_lattice)
3.2. Visualize the lattice for shaft and corridor growth
# 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
# Visualize the mesh using pyvista plotter
# initiating the plotter
p = pv.Plotter(notebook=True)
# fast visualization of the lattice
avail_lattice_bool.fast_vis(p)
# plotting
p.window_size = 1000, 1000
# p.screenshot("corridors_1")
p.show(use_ipyvtk=True)
4.1. calculating cluster centers for each agent, based on agent size
cluster_centers = []
for func_id in range(0, int(np.max(frames[-1])) + 1):
if not func_id in [2, 13, 14, 19, 20]: #again not caring about specific functions, such as the garden
function_locs = np.array(np.where(frames[-1] == func_id)).T
# for every multiple of 100 voxels a cluster center is created and at least 1 per agent
nclusters = int(len(function_locs) / 100 + 1)
kmeans_model = KMeans(n_clusters= nclusters, random_state=0).fit(function_locs)
cluster_centers.append(np.round(kmeans_model.cluster_centers_).astype(np.int8))
# making shure it's no longer a nested cluster center list
cluster_center_list = []
for i in range(len(cluster_centers)):
cluster_center_list.append(cluster_centers[i][0])
#creating a lattice that shows all cluster centers
cluster_center_lattice = avail_lattice_bool * 0
for cluster_center in cluster_center_list:
cluster_center_lattice[cluster_center[0], cluster_center[1], cluster_center[2]] = 1
4.2. Visualizing the cluster centers
p = pv.Plotter(notebook=True)
base_lattice = cluster_center_lattice
# 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 avilability lattice
avail_lattice_bool.fast_vis(p)
# Add the data values to the cell data
grid.cell_arrays["Agents"] = base_lattice.flatten(order="F").astype(int) # Flatten the array!
# filtering the voxels
threshed = grid.threshold([0.9, 1.1])
# adding the voxels
p.add_mesh(threshed, name='sphere', show_edges=True, opacity=1.0, show_scalar_bar=False, color = [1,1,0])
p.window_size = 1000, 1000
# p.screenshot("corridors_2")
p.show(use_ipyvtk=True)
4.3. Calculating 2D shaft placement
kmeans_model = KMeans(n_clusters=6, random_state=0).fit(cluster_center_list) # select how many shafts
cluster_centers_total = np.round(kmeans_model.cluster_centers_).astype(np.int8)
# init shaft lattice
shft_lattice = avail_lattice_bool * 0
# setting the shafts to 1
for cl_cen in cluster_centers_total:
shft_lattice[cl_cen[0],cl_cen[1],:] = 1
4.4. Visualizing the shaft lattice
p = pv.Plotter(notebook=True)
base_lattice = shft_lattice * avail_lattice_bool
# 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 avilability lattice
avail_lattice_bool.fast_vis(p)
# Add the data values to the cell data
grid.cell_arrays["Agents"] = base_lattice.flatten(order="F").astype(int) # Flatten the array!
# filtering the voxels
threshed = grid.threshold([0.9, 1.1])
# adding the voxels
p.add_mesh(threshed, name='sphere', show_edges=True, opacity=1.0, show_scalar_bar=False)
p.window_size = 1000, 1000
# p.screenshot("corridors_3")
p.show(use_ipyvtk=True)
5.1. Adjacency matrix creation
avail_lattice_flat = avail_lattice_bool.flatten()
# find the number of all voxels
vox_count = avail_lattice_bool.size
# initialize the adjacency matrix
adj_list = []
# Finding the index of the available voxels in avail_lattice
avail_index = np.array(np.where(avail_lattice_bool == 1)).T
# fill the adjacency matrix using the list of all neighbours
for vox_loc in avail_index:
# find the 1D id
vox_id = np.ravel_multi_index(vox_loc, avail_lattice_bool.shape)
# retrieve the list of neighbours of the voxel based on the stencil
vox_neighs = init_avail_lattice.find_neighbours_masked(stencil, loc = vox_loc)
# iterating over the neighbours
for neigh in vox_neighs:
if avail_lattice_flat[neigh] == 1:
adj_list.append([1.0, vox_id, neigh])
#list to array
adj_array = np.array(adj_list).T
#array to sparce matrix
adj_matrix_sparse = csr_matrix((adj_array[0], (adj_array[1], adj_array[2])), shape=(vox_count,vox_count))
#sparce matrix to nx connectivity graph
g = nx.from_scipy_sparse_matrix(adj_matrix_sparse)
5.2. Corridor growth
occ_ind = cluster_center_list
# initialize corridor lattice
cor_lattice = shft_lattice * 0
cor_flat = cor_lattice.flatten()
# for each voxel that needs to have access to shafts
for a_vox in occ_ind:
# slice the corridor lattice horizontally
cor_floor = shft_lattice[:,:,a_vox[2]]
# find the vertical shaft voxel indices
shaft_vox_inds = np.array(np.where(cor_floor > 0)).T
paths = []
path_lens = []
for shft_ind in shaft_vox_inds:
# construct the destination address
dst_vox = np.array([shft_ind[0],shft_ind[1],a_vox[2]])
# construct 1-dimensional indices
src_ind = np.ravel_multi_index(a_vox, shft_lattice.shape)
dst_ind = np.ravel_multi_index(dst_vox, shft_lattice.shape)
# find the shortest path
if nx.algorithms.shortest_paths.generic.has_path(g, src_ind, dst_ind):
path = nx.algorithms.shortest_paths.astar.astar_path(g, src_ind, dst_ind)
paths.append(path)
path_lens.append(len(path))
# find the number of shortest path
for shortest_path_index in np.array(path_lens).argsort()[:1]: #select how many paths to connect to closest shafts
cor_flat[paths[shortest_path_index]] = 1
#set the floor level the shafts need to be connected
connected_floor_levels = [2]
for level in connected_floor_levels:
cur_floor_level = shft_lattice[:,:,level] #finding floor level
shaft_vox_inds = np.array(np.where(cur_floor_level > 0)).T #finding shafts
dst_ind_list = []
source_ind_list = []
for shft_ind in shaft_vox_inds:
# construct the destination address
dst_vox = np.array([shft_ind[0],shft_ind[1],level])
index = np.ravel_multi_index(dst_vox, shft_lattice.shape) #finding the 1D index
#setting the indices as both src and dst
dst_ind_list.append(index)
source_ind_list.append(index)
for dst_ind in dst_ind_list:
paths = []
path_lens = []
for source_ind in source_ind_list: #checking if path to itself
if dst_ind == source_ind:
continue
if nx.algorithms.shortest_paths.generic.has_path(g, source_ind, dst_ind):
path = nx.algorithms.shortest_paths.astar.astar_path(g, source_ind, dst_ind) #finding the closest paths
paths.append(path)
path_lens.append(len(path))
# find the 2 shortest paths
for shortest_path_index in np.array(path_lens).argsort()[:2]: # selecting how many paths to connect to closests shafts
cor_flat[paths[shortest_path_index]] = 1
# reshape the flat lattice
cor_lattice = cor_flat.reshape(cor_lattice.shape)
5.3. Vizualizing the shafts and corridors lattice
p = pv.Plotter(notebook=True)
base_lattice = (shft_lattice + cor_lattice) * avail_lattice_bool
# 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 avilability lattice
avail_lattice_bool.fast_vis(p)
# Add the data values to the cell data
grid.cell_arrays["Agents"] = base_lattice.flatten(order="F").astype(int) # Flatten the array!
# filtering the voxels
threshed = grid.threshold([0.9, 2.1])
# adding the voxels
p.add_mesh(threshed, name='sphere', show_edges=True, opacity=1.0, show_scalar_bar=False)
p.window_size = 1000, 1000
p.screenshot("corridors_4")
p.show(use_ipyvtk=True)
6.1. Setting the shafts as 1 and corridors as 2 for export
shafts_and_corridors_lattice = avail_lattice_bool * 0
shafts_and_corridors_lattice[np.where(cor_lattice == 1)] = 2
shafts_and_corridors_lattice[np.where(shft_lattice == 1)] = 1
shafts_and_corridors_lattice *= avail_lattice_bool
6.2. Exporting the lattices
csv_path = os.path.relpath('../data/324/shafts_and_corridors_lattice.csv')
shafts_and_corridors_lattice.to_csv(csv_path)
csv_path = os.path.relpath('../data/324/avail_lattice_bool.csv')
avail_lattice_bool.to_csv(csv_path)
Credits
__author__ = "Shervin Azadi and Pirouz Nourian"
__license__ = "MIT"
__version__ = "1.0"
__url__ = "https://github.com/shervinazadi/spatial_computing_workshops"
__summary__ = "Spatial Computing Design Studio Workshop on MCDA and Path Finding for Generative Spatial Relations"