Skip to content

Facade closeness notebook

In this notebook a csv file is created that maps closeness to the facade voxels to the rest of the voxels

import os
import topogenesis as tg
import pyvista as pv
import trimesh as tm
import numpy as np
import networkx as nx
from scipy.sparse import csr_matrix

1.1. Specify the facade

# creating neighborhood definition
s_sum = tg.create_stencil("von_neumann", 1, 1)
s_sum2 = tg.create_stencil("von_neumann", 1, 1)

# NW
# s_sum.set_index([0,0,0], 0)
# s_sum.set_index([0,0,1], 0)
# s_sum.set_index([0,-1,0], 0)
# s_sum.set_index([0,1,0], 0)
# s_sum.set_index([1,0,0], 0)

# s_sum2.set_index([0,0,0], 0)
# s_sum2.set_index([0,0,1], 0)
# s_sum2.set_index([0,1,0], 0)
# s_sum2.set_index([1,0,0], 0)

# SW
# s_sum.set_index([0,0,0], 0)
# s_sum.set_index([0,0,1], 0)
# s_sum.set_index([0,1,0], 0)
# s_sum.set_index([1,0,0], 0)
# s_sum.set_index([-1,0,0], 0)

# #SE
s_sum.set_index([0,0,0], 0)
s_sum.set_index([0,0,1], 0)
s_sum.set_index([0,-1,0], 0)
s_sum.set_index([0,1,0], 0)
s_sum.set_index([-1,0,0], 0)

# NE
# s_sum.set_index([0,0,0], 0)
# s_sum.set_index([0,0,1], 0)
# s_sum.set_index([0,-1,0], 0)
# s_sum.set_index([1,0,0], 0)
# s_sum.set_index([-1,0,0], 0)

# s_sum2.set_index([0,0,0], 0)
# s_sum2.set_index([0,0,1], 0)
# s_sum2.set_index([0,-1,0], 0)
# s_sum2.set_index([1,0,0], 0)

# Full facade
# s_sum.set_index([0,0,0], 0)
# s_sum.set_index([0,0,1], 0)
# s_sum.set_index([0,0,-1], 0)


s_sum.function = tg.sfunc.sum
s_sum2.function = tg.sfunc.sum


s_plus = tg.create_stencil("von_neumann", 1, 1)
s_plus.set_index([0,0,0], 0)
s_plus.set_index([0,0,1], 0)
s_plus.set_index([0,0,-1], 0)
# s_plus.set_index([0,1,0], 0)
# s_plus.set_index([0,-1,0], 0)
# s_plus.set_index([1,0,0], 0)

print(s_plus)

1.2. Loading environmental data

# loading the lattice from csv
lattice_path = os.path.relpath('../data/solar_envelope_324.csv')
solar_envelope = tg.lattice_from_csv(lattice_path)
avail_lattice = tg.lattice_from_csv(lattice_path)

#Load context mesh from csv
context_path = os.path.relpath("../data/immediate_context.obj")
context_mesh = tm.load(context_path)

# construct full latice
full_lattice = solar_envelope * 0 + 1

2.1. Finding the facade voxels based on the stencils

neighbour_sum = solar_envelope.apply_stencil(s_sum, border_condition='pad_outside', padding_value=0)
exterior_condition = neighbour_sum < 1 #1 for specific facade, 4 for full facade

# for specific facades a second sum function is added to remove unwanted voxels that would otherwise count as 'facades'
neighbour_sum2 = solar_envelope.apply_stencil(s_sum2, border_condition='pad_outside', padding_value=0)
exterior_condition2 = neighbour_sum2 < 1 

# flattening the lattices
exterior_condition_flat = exterior_condition.flatten()
exterior_condition2_flat = exterior_condition2.flatten()

# removing weird corner and edge voxels that would otherwise count as a facade voxel
exterior_condition2_id = np.array(np.where(exterior_condition2_flat == 1)).T
exterior_condition_flat[exterior_condition2_id] = 0

#reshaping the list
exterior_condition = exterior_condition_flat.reshape(exterior_condition.shape)

#finding the indices
avail_index = np.array(np.where(avail_lattice == 1)).T

# creating the border lattice
border_lattice = exterior_condition
border_lattice = solar_envelope * exterior_condition

# setting everything above ground level to 0 (for the all facades function)
border_lattice[:,:,2:] = 0

2.2. Visualizing the border voxels

pv.set_plot_theme("document")

# 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
border_lattice.fast_vis(p)

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

# plotting
cpos = [(785.8704805788776, 708.4540755788776, 741.8613927288776),
 (65.08283250000001, -12.333572500000002, 21.07374465),
 (0.0, 0.0, 1.0)]
p.camera_position = cpos
p.window_size = 2000, 2000
p.show(use_ipyvtk=True)
# p.screenshot("facade_1")
print(p.camera_position)

3.1. Adjacency network creation and calculation

# creating an adjacency list for the creation of a sparce adjacency matrix

avail_lattice_flat = avail_lattice.flatten()

# find the number of all voxels
vox_count = full_lattice.size 
print(vox_count)

# initialize the adjacency list
adj_list = []

# 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(s_plus, 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])

#shaping the list into an array
adj_array = np.array(adj_list).T

#creating a sparce row matrix from the array
adj_matrix_sparse =  csr_matrix((adj_array[0], (adj_array[1], adj_array[2])), shape=(vox_count,vox_count)) 

#turning the sparce adjacency matrix into a nx node network
g = nx.from_scipy_sparse_matrix(adj_matrix_sparse)

border_lattice_flat = border_lattice.flatten()

#finding all the indices that count as 'facade voxels'
id_list = np.array(np.where(border_lattice_flat == 1)).T

#creating a dict that only holds all the lowest key values

counter = 0
for id_border_vox in id_list:
    shortest_path = nx.shortest_path_length(g, id_border_vox[0])
    for key in shortest_path:
        if counter == 0:
            shortest_paths = shortest_path
        else:
            shortest_paths[key] = min(shortest_paths[key], shortest_path[key])
    counter += 1

3.2. Mapping the distance field

# mapping all the distance values to closeness values and reshaping to the lattice shape

all_values = shortest_paths.values()
max_value = max(all_values)

#initializing a 1d distance list
base_lattice_flat = np.zeros(vox_count)

for key in shortest_paths:  
    mapped_value = 1 - shortest_paths[key] / max_value #mapping value between 0 and 1
    base_lattice_flat[key] = mapped_value

# creating the 2D lattice slice
vox_interest_acc_lattice = tg.to_lattice(base_lattice_flat.reshape(avail_lattice.shape), avail_lattice)

3.3. Visualizing the 2D facade closeness slice

# 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

# load the mesh from file
context_path = os.path.relpath('../data/immediate_context.obj')
context_mesh = tm.load(context_path)

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

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

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

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

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

# adding the volume
opacity = np.array([0.0,0.6,0.6,0.6,0.6,0.6,0.6])
p.add_volume(grid, cmap="coolwarm", clim=[0.0, 1.0] ,opacity=opacity)

# plotting
cpos = [(785.8704805788776, 708.4540755788776, 741.8613927288776),
 (65.08283250000001, -12.333572500000002, 21.07374465),
 (0.0, 0.0, 1.0)]
p.camera_position = cpos
p.window_size = 2000, 2000
p.show(use_ipyvtk=True)
p.screenshot("facade_2")
print(p.camera_position)

4.1. Mapping the 2D closeness slice to the 3D lattice

# mapping the 2d facade distance field to the 3D lattice shape

for i in range(len(vox_interest_acc_lattice[0][0])):
    vox_interest_acc_lattice[:,:,i] = vox_interest_acc_lattice[:,:,1]

# shaping the lattice to the solar anvelope
vox_interest_acc_lattice = vox_interest_acc_lattice * avail_lattice

4.2. Visualizing the facade closeness lattice

# 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

# load the mesh from file
context_path = os.path.relpath('../data/immediate_context.obj')
context_mesh = tm.load(context_path)

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

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

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

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

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

# adding the volume
opacity = np.array([0.0,0.6,0.6,0.6,0.6,0.6,0.6])
p.add_volume(grid, cmap="coolwarm", clim=[0.0, 1.0] ,opacity=opacity)

# plotting
cpos = [(785.8704805788776, 708.4540755788776, 741.8613927288776),
 (65.08283250000001, -12.333572500000002, 21.07374465),
 (0.0, 0.0, 1.0)]
p.camera_position = cpos
p.window_size = 2000, 2000
p.show(use_ipyvtk=True)
p.screenshot("facade_3")
print(p.camera_position)

5.1 .Saving the lattice to a csv

#saving the facade lattice to a csv file
csv_path = os.path.relpath('../data/closeness_to_facade_324.csv')
vox_interest_acc_lattice.to_csv(csv_path)