Module geodome.tessellation
Expand source code
import numpy as np
from numba import njit
from numba.typed import Dict
from typing import Union
import math
try:
from geodome.neighbourhood_search import *
except:
from neighbourhood_search import *
SCALE = 1
@njit
def find_adjacent_triangles(triangles: np.ndarray, vertices: np.ndarray) -> np.ndarray:
"""Finds adjacent triangles for a given point
Args:
triangles (np.ndarray): The list of triangles in the dome
vertices (np.ndarray): The vertices (index) to find the adjacent triangles for
Returns:
np.ndarray: List of indices for triangles in the array
"""
found_triangles = np.zeros((len(triangles)), dtype=np.int64)
found_index = 0
for i in range(len(triangles)):
t = triangles[i]
if t[0] in vertices or t[1] in vertices or t[2] in vertices:
found_triangles[found_index] = i
found_index += 1
resized = np.zeros(found_index, dtype=np.int64)
for i in range(found_index):
resized[i] = found_triangles[i]
return resized
@njit
def is_zero(coord, message="") -> None:
"""Debug method, checks if a coord is at origin
Args:
coord (np array): coordinate to check
message (str, optional): descriptive message when triggered. Defaults to "".
"""
if (coord == np.array([0, 0, 0])).all():
print("Something is zero that shouldn't be: " + message)
@njit
def normalise_length(coords: np.ndarray) -> np.ndarray:
"""Normalises the distance from origin of a coord. Multiplies by the
frequency the icosphere to avoid floating point precision errors
Args:
coords (np.ndarray): coordinate to normalise
Returns:
np.ndarray: normalised coordinate
"""
length = math.sqrt(
math.pow(coords[0], 2) + math.pow(coords[1], 2) + math.pow(coords[2], 2)
)
is_zero(coords, "normalise")
return np.array(
[
(coords[0] / length) * SCALE,
(coords[1] / length) * SCALE,
(coords[2] / length) * SCALE,
]
)
@njit
def get_middle_coords(v1: np.ndarray, v2: np.ndarray) -> np.ndarray:
"""Gets the midpoint between two coords
Args:
v1 (np.ndarray): coord 1
v2 (np.ndarray): coord 2
Returns:
np.ndarray: the midpoint (not normalised)
"""
ret = np.array(
[(v2[0] + v1[0]) / 2, (v2[1] + v1[1]) / 2, (v2[2] + v1[2]) / 2],
dtype=np.float64,
)
return ret
@njit
def old_add_middle_get_index(
matrix: np.ndarray,
new_vertices: np.ndarray,
old_vertices: np.ndarray,
v_index: np.int64,
v1: np.int64,
v2: np.int64,
) -> Union[np.int64, np.int64]:
"""Matrix based method for inserting new vertices with the guarantee of uniqueness
Args:
matrix (np.ndarray): (Nvertices x Nvertices) matrix
new_vertices (np.ndarray): Array to store new vertices
old_vertices (np.ndarray): Array of old vertices
v_index (np.int64): Current index of new vertices
v1 (np.int64): Index of outer vertex
v2 (np.int64): Index of other outer vertex
Returns:
Union[np.int64, np.int64]: index of inserted vertex, new v_index
"""
if matrix[v1][v2] == -1:
matrix[v1][v2] = v_index
matrix[v2][v1] = v_index
new_vertices[v_index] = get_middle_coords(old_vertices[v1], old_vertices[v2])
v_index += 1
return matrix[v1][v2], v_index
@njit
def add_vertex_get_index(mid: np.ndarray, vertices: np.ndarray, v_dict: Dict):
"""Adds a given midpoint to a list of new vertices
Args:
mid (np.ndarray): the midpoint to add
vertices (np.ndarray): an array of new vertices, to be concatenated with
existing vertices
v_index (Dict): dictionary containing the indexes of existing midpoints,
to prevent duplicates
Returns:
int: the index of the midpoint that was added to vertices
"""
# Creating (hopefully) unique key for each coordinate
mid_sum = mid[0] * 100 + mid[1] * 10 + mid[2]
# Adding key to dictionary of coords with index as value
inserted = False
while not inserted:
if mid_sum not in v_dict:
v_dict[mid_sum] = len(v_dict)
inserted = True
else:
if (mid != vertices[v_dict[mid_sum]]).any():
mid_sum += 100
else:
break
index = v_dict[mid_sum]
# Add new midpoint to new vertices array
if (vertices[index] == np.array([0, 0, 0])).all():
vertices[index] = mid
return index
@njit
def calc_dist(points) -> list:
"""Calculates the distance of each point in the Dome from the origin
Args:
points (list): list of points
Returns:
list: List of distances of each point in the Dome
"""
distances = []
for p in points:
dist = math.sqrt(p[0] * p[0] + p[1] * p[1] + p[2] * p[2])
distances.append(round(dist, 2))
return distances
@njit
def tessellate_geodesic_dome(
vertices: np.ndarray,
triangles: np.ndarray,
target_triangles: np.ndarray = np.zeros(0, dtype=np.int64),
) -> Union[np.ndarray, np.ndarray]:
"""Tessellates the target triangles in the dome
Args:
vertices (np.ndarray): Array of vertices
triangles (np.ndarray): Array of triangles
target_triangles (np.ndarray): Target triangles to tessellate
(optional: tessellates all triangles if not given)
Returns:
Union[np.ndarray, np.ndarray, np.ndarray]: Array of new vertices,
triangles, and new adjacency list
"""
# Fill target with all triangles if not given
if len(target_triangles) == 0:
target_triangles = np.arange(len(triangles))
# First calculate the number of edges covered by the current triangles
# create new array for new triangles
n_new_triangles = len(triangles) + 3 * len(target_triangles)
new_triangles = np.zeros((n_new_triangles, 3), dtype=np.int64)
n_old_vertices = len(vertices)
generous_edge_count = len(target_triangles) * 3 + n_old_vertices
# create new array for new vertices
new_vertices = np.zeros((generous_edge_count, 3), dtype=np.float64)
v_dict = Dict.empty(key_type=np.float64, value_type=np.int64)
i = 0
hit_triangles = np.zeros((len(triangles)), dtype=np.int8)
# Add all existing vertices to dictionary
for v in vertices:
add_vertex_get_index(v, new_vertices, v_dict)
for t in target_triangles:
tri = triangles[t]
hit_triangles[t] = 1
v0 = vertices[tri[0]]
v1 = vertices[tri[1]]
v2 = vertices[tri[2]]
# Get midpoints for each edge of the triangle and normalise
mid01 = normalise_length(get_middle_coords(v0, v1))
mid12 = normalise_length(get_middle_coords(v1, v2))
mid02 = normalise_length(get_middle_coords(v0, v2))
index01 = add_vertex_get_index(mid01, new_vertices, v_dict)
index12 = add_vertex_get_index(mid12, new_vertices, v_dict)
index02 = add_vertex_get_index(mid02, new_vertices, v_dict)
# Create new triangles
new_triangles[i] = [tri[0], index01, index02]
new_triangles[i + 1] = [tri[1], index12, index01]
new_triangles[i + 2] = [tri[2], index02, index12]
new_triangles[i + 3] = [index01, index12, index02]
i += 4
# Add all the untargeted triangles
for j in range(len(triangles)):
if hit_triangles[j] == 0:
t = triangles[j]
new_triangles[i] = t
i += 1
# Create array to concatenate old vertices with new midpoints
vertices = np.zeros((len(v_dict), 3), dtype=np.float64)
for i in range(len(v_dict)):
vertices[i] = new_vertices[i]
# i = 0
# # Add old vertices
# for v in old_vertices:
# vertices[i] = v
# i += 1
# # Add new midpoints
# for j in range(v_index):
# vertices[i] = new_vertices[j]
# i += 1
new_adj_list = create_adj_list(vertices, new_triangles)
return vertices, new_triangles, new_adj_list
@njit
def create_geodesic_dome(freq=0) -> Union[np.ndarray, np.ndarray, np.ndarray]:
"""Creates an geodesic dome of a given frequency
Args:
freq (int, optional): the frequency of the dome. Defaults to 0.
Returns:
Union[np.ndarray, np.ndarray]: the array of vertices, the array of
triangles
Vertices = [[x,y,z], ... , [x,y,z]]
Triangles = [[v1, v2, v3], ...] where vx is the index of a vertex in the vertices array
Adjacency list = [[v1, ..., v5, v6?], ...] where vx is the index of a vertex. v6 may not exist for some vertices
"""
# Set normalised scaling
if freq != 0:
SCALE = freq
g_ratio = (1 + math.sqrt(5)) / 2
# creating initial icosahedron vertices
icosa_vertices = np.array(
[
(-1, g_ratio, 0),
(1, g_ratio, 0),
(-1, -(g_ratio), 0),
(1, -(g_ratio), 0),
(0, -1, g_ratio),
(0, 1, g_ratio),
(0, -1, -(g_ratio)),
(0, 1, -(g_ratio)),
(g_ratio, 0, -1),
(g_ratio, 0, 1),
(-(g_ratio), 0, -1),
(-(g_ratio), 0, 1),
],
dtype=np.float64,
)
# creating initial icosahedron edges
triangles = np.array(
[
(0, 11, 5),
(0, 5, 1),
(0, 1, 7),
(0, 7, 10),
(0, 10, 11),
(1, 5, 9),
(5, 11, 4),
(11, 10, 2),
(10, 7, 6),
(7, 1, 8),
(3, 9, 4),
(3, 4, 2),
(3, 2, 6),
(3, 6, 8),
(3, 8, 9),
(4, 9, 5),
(2, 4, 11),
(6, 2, 10),
(8, 6, 7),
(9, 8, 1),
],
dtype=np.int64,
)
# Array for normalised vertices
vertices = np.zeros((len(icosa_vertices), 3), dtype=np.float64)
# Normalise all icosahedron vertices
for i in range(len(icosa_vertices)):
vertices[i] = normalise_length(icosa_vertices[i])
# Tessellate icosahedron
adj_list = create_adj_list(vertices, triangles)
for i in range(freq):
vertices, triangles, adj_list = tessellate_geodesic_dome(vertices, triangles)
return vertices, triangles, adj_list
Functions
def add_vertex_get_index(mid: numpy.ndarray, vertices: numpy.ndarray, v_dict: numba.typed.typeddict.Dict)-
Adds a given midpoint to a list of new vertices
Args
mid:np.ndarray- the midpoint to add
vertices:np.ndarray- an array of new vertices, to be concatenated with
- existing vertices
v_index:Dict- dictionary containing the indexes of existing midpoints,
to prevent duplicates
Returns
int- the index of the midpoint that was added to vertices
Expand source code
@njit def add_vertex_get_index(mid: np.ndarray, vertices: np.ndarray, v_dict: Dict): """Adds a given midpoint to a list of new vertices Args: mid (np.ndarray): the midpoint to add vertices (np.ndarray): an array of new vertices, to be concatenated with existing vertices v_index (Dict): dictionary containing the indexes of existing midpoints, to prevent duplicates Returns: int: the index of the midpoint that was added to vertices """ # Creating (hopefully) unique key for each coordinate mid_sum = mid[0] * 100 + mid[1] * 10 + mid[2] # Adding key to dictionary of coords with index as value inserted = False while not inserted: if mid_sum not in v_dict: v_dict[mid_sum] = len(v_dict) inserted = True else: if (mid != vertices[v_dict[mid_sum]]).any(): mid_sum += 100 else: break index = v_dict[mid_sum] # Add new midpoint to new vertices array if (vertices[index] == np.array([0, 0, 0])).all(): vertices[index] = mid return index def calc_dist(points) ‑> list-
Calculates the distance of each point in the Dome from the origin
Args
points:list- list of points
Returns
list- List of distances of each point in the Dome
Expand source code
@njit def calc_dist(points) -> list: """Calculates the distance of each point in the Dome from the origin Args: points (list): list of points Returns: list: List of distances of each point in the Dome """ distances = [] for p in points: dist = math.sqrt(p[0] * p[0] + p[1] * p[1] + p[2] * p[2]) distances.append(round(dist, 2)) return distances def create_geodesic_dome(freq=0) ‑> numpy.ndarray-
Creates an geodesic dome of a given frequency
Args
freq:int, optional- the frequency of the dome. Defaults to 0.
Returns
Union[np.ndarray, np.ndarray]- the array of vertices, the array of
triangles
Vertices = [[x,y,z], … , [x,y,z]] Triangles = [[v1, v2, v3], …] where vx is the index of a vertex in the vertices array
Adjacency list = [[v1, …, v5, v6?], …] where vx is the index of a vertex. v6 may not exist for some vertices
Expand source code
@njit def create_geodesic_dome(freq=0) -> Union[np.ndarray, np.ndarray, np.ndarray]: """Creates an geodesic dome of a given frequency Args: freq (int, optional): the frequency of the dome. Defaults to 0. Returns: Union[np.ndarray, np.ndarray]: the array of vertices, the array of triangles Vertices = [[x,y,z], ... , [x,y,z]] Triangles = [[v1, v2, v3], ...] where vx is the index of a vertex in the vertices array Adjacency list = [[v1, ..., v5, v6?], ...] where vx is the index of a vertex. v6 may not exist for some vertices """ # Set normalised scaling if freq != 0: SCALE = freq g_ratio = (1 + math.sqrt(5)) / 2 # creating initial icosahedron vertices icosa_vertices = np.array( [ (-1, g_ratio, 0), (1, g_ratio, 0), (-1, -(g_ratio), 0), (1, -(g_ratio), 0), (0, -1, g_ratio), (0, 1, g_ratio), (0, -1, -(g_ratio)), (0, 1, -(g_ratio)), (g_ratio, 0, -1), (g_ratio, 0, 1), (-(g_ratio), 0, -1), (-(g_ratio), 0, 1), ], dtype=np.float64, ) # creating initial icosahedron edges triangles = np.array( [ (0, 11, 5), (0, 5, 1), (0, 1, 7), (0, 7, 10), (0, 10, 11), (1, 5, 9), (5, 11, 4), (11, 10, 2), (10, 7, 6), (7, 1, 8), (3, 9, 4), (3, 4, 2), (3, 2, 6), (3, 6, 8), (3, 8, 9), (4, 9, 5), (2, 4, 11), (6, 2, 10), (8, 6, 7), (9, 8, 1), ], dtype=np.int64, ) # Array for normalised vertices vertices = np.zeros((len(icosa_vertices), 3), dtype=np.float64) # Normalise all icosahedron vertices for i in range(len(icosa_vertices)): vertices[i] = normalise_length(icosa_vertices[i]) # Tessellate icosahedron adj_list = create_adj_list(vertices, triangles) for i in range(freq): vertices, triangles, adj_list = tessellate_geodesic_dome(vertices, triangles) return vertices, triangles, adj_list def find_adjacent_triangles(triangles: numpy.ndarray, vertices: numpy.ndarray) ‑> numpy.ndarray-
Finds adjacent triangles for a given point
Args
triangles:np.ndarray- The list of triangles in the dome
vertices:np.ndarray- The vertices (index) to find the adjacent triangles for
Returns
np.ndarray- List of indices for triangles in the array
Expand source code
@njit def find_adjacent_triangles(triangles: np.ndarray, vertices: np.ndarray) -> np.ndarray: """Finds adjacent triangles for a given point Args: triangles (np.ndarray): The list of triangles in the dome vertices (np.ndarray): The vertices (index) to find the adjacent triangles for Returns: np.ndarray: List of indices for triangles in the array """ found_triangles = np.zeros((len(triangles)), dtype=np.int64) found_index = 0 for i in range(len(triangles)): t = triangles[i] if t[0] in vertices or t[1] in vertices or t[2] in vertices: found_triangles[found_index] = i found_index += 1 resized = np.zeros(found_index, dtype=np.int64) for i in range(found_index): resized[i] = found_triangles[i] return resized def get_middle_coords(v1: numpy.ndarray, v2: numpy.ndarray) ‑> numpy.ndarray-
Gets the midpoint between two coords
Args
v1:np.ndarray- coord 1
v2:np.ndarray- coord 2
Returns
np.ndarray- the midpoint (not normalised)
Expand source code
@njit def get_middle_coords(v1: np.ndarray, v2: np.ndarray) -> np.ndarray: """Gets the midpoint between two coords Args: v1 (np.ndarray): coord 1 v2 (np.ndarray): coord 2 Returns: np.ndarray: the midpoint (not normalised) """ ret = np.array( [(v2[0] + v1[0]) / 2, (v2[1] + v1[1]) / 2, (v2[2] + v1[2]) / 2], dtype=np.float64, ) return ret def is_zero(coord, message='') ‑> None-
Debug method, checks if a coord is at origin
Args
coord:np array- coordinate to check
message:str, optional- descriptive message when triggered. Defaults to "".
Expand source code
@njit def is_zero(coord, message="") -> None: """Debug method, checks if a coord is at origin Args: coord (np array): coordinate to check message (str, optional): descriptive message when triggered. Defaults to "". """ if (coord == np.array([0, 0, 0])).all(): print("Something is zero that shouldn't be: " + message) def normalise_length(coords: numpy.ndarray) ‑> numpy.ndarray-
Normalises the distance from origin of a coord. Multiplies by the frequency the icosphere to avoid floating point precision errors
Args
coords:np.ndarray- coordinate to normalise
Returns
np.ndarray- normalised coordinate
Expand source code
@njit def normalise_length(coords: np.ndarray) -> np.ndarray: """Normalises the distance from origin of a coord. Multiplies by the frequency the icosphere to avoid floating point precision errors Args: coords (np.ndarray): coordinate to normalise Returns: np.ndarray: normalised coordinate """ length = math.sqrt( math.pow(coords[0], 2) + math.pow(coords[1], 2) + math.pow(coords[2], 2) ) is_zero(coords, "normalise") return np.array( [ (coords[0] / length) * SCALE, (coords[1] / length) * SCALE, (coords[2] / length) * SCALE, ] ) def old_add_middle_get_index(matrix: numpy.ndarray, new_vertices: numpy.ndarray, old_vertices: numpy.ndarray, v_index: numpy.int64, v1: numpy.int64, v2: numpy.int64) ‑> numpy.int64-
Matrix based method for inserting new vertices with the guarantee of uniqueness
Args
matrix:np.ndarray- (Nvertices x Nvertices) matrix
new_vertices:np.ndarray- Array to store new vertices
old_vertices:np.ndarray- Array of old vertices
v_index:np.int64- Current index of new vertices
v1:np.int64- Index of outer vertex
v2:np.int64- Index of other outer vertex
Returns
Union[np.int64, np.int64]- index of inserted vertex, new v_index
Expand source code
@njit def old_add_middle_get_index( matrix: np.ndarray, new_vertices: np.ndarray, old_vertices: np.ndarray, v_index: np.int64, v1: np.int64, v2: np.int64, ) -> Union[np.int64, np.int64]: """Matrix based method for inserting new vertices with the guarantee of uniqueness Args: matrix (np.ndarray): (Nvertices x Nvertices) matrix new_vertices (np.ndarray): Array to store new vertices old_vertices (np.ndarray): Array of old vertices v_index (np.int64): Current index of new vertices v1 (np.int64): Index of outer vertex v2 (np.int64): Index of other outer vertex Returns: Union[np.int64, np.int64]: index of inserted vertex, new v_index """ if matrix[v1][v2] == -1: matrix[v1][v2] = v_index matrix[v2][v1] = v_index new_vertices[v_index] = get_middle_coords(old_vertices[v1], old_vertices[v2]) v_index += 1 return matrix[v1][v2], v_index def tessellate_geodesic_dome(vertices: numpy.ndarray, triangles: numpy.ndarray, target_triangles: numpy.ndarray = array([], dtype=int64)) ‑> numpy.ndarray-
Tessellates the target triangles in the dome
Args
vertices:np.ndarray- Array of vertices
triangles:np.ndarray- Array of triangles
target_triangles:np.ndarray- Target triangles to tessellate (optional: tessellates all triangles if not given)
Returns
Union[np.ndarray, np.ndarray, np.ndarray]- Array of new vertices,
triangles, and new adjacency list
Expand source code
@njit def tessellate_geodesic_dome( vertices: np.ndarray, triangles: np.ndarray, target_triangles: np.ndarray = np.zeros(0, dtype=np.int64), ) -> Union[np.ndarray, np.ndarray]: """Tessellates the target triangles in the dome Args: vertices (np.ndarray): Array of vertices triangles (np.ndarray): Array of triangles target_triangles (np.ndarray): Target triangles to tessellate (optional: tessellates all triangles if not given) Returns: Union[np.ndarray, np.ndarray, np.ndarray]: Array of new vertices, triangles, and new adjacency list """ # Fill target with all triangles if not given if len(target_triangles) == 0: target_triangles = np.arange(len(triangles)) # First calculate the number of edges covered by the current triangles # create new array for new triangles n_new_triangles = len(triangles) + 3 * len(target_triangles) new_triangles = np.zeros((n_new_triangles, 3), dtype=np.int64) n_old_vertices = len(vertices) generous_edge_count = len(target_triangles) * 3 + n_old_vertices # create new array for new vertices new_vertices = np.zeros((generous_edge_count, 3), dtype=np.float64) v_dict = Dict.empty(key_type=np.float64, value_type=np.int64) i = 0 hit_triangles = np.zeros((len(triangles)), dtype=np.int8) # Add all existing vertices to dictionary for v in vertices: add_vertex_get_index(v, new_vertices, v_dict) for t in target_triangles: tri = triangles[t] hit_triangles[t] = 1 v0 = vertices[tri[0]] v1 = vertices[tri[1]] v2 = vertices[tri[2]] # Get midpoints for each edge of the triangle and normalise mid01 = normalise_length(get_middle_coords(v0, v1)) mid12 = normalise_length(get_middle_coords(v1, v2)) mid02 = normalise_length(get_middle_coords(v0, v2)) index01 = add_vertex_get_index(mid01, new_vertices, v_dict) index12 = add_vertex_get_index(mid12, new_vertices, v_dict) index02 = add_vertex_get_index(mid02, new_vertices, v_dict) # Create new triangles new_triangles[i] = [tri[0], index01, index02] new_triangles[i + 1] = [tri[1], index12, index01] new_triangles[i + 2] = [tri[2], index02, index12] new_triangles[i + 3] = [index01, index12, index02] i += 4 # Add all the untargeted triangles for j in range(len(triangles)): if hit_triangles[j] == 0: t = triangles[j] new_triangles[i] = t i += 1 # Create array to concatenate old vertices with new midpoints vertices = np.zeros((len(v_dict), 3), dtype=np.float64) for i in range(len(v_dict)): vertices[i] = new_vertices[i] # i = 0 # # Add old vertices # for v in old_vertices: # vertices[i] = v # i += 1 # # Add new midpoints # for j in range(v_index): # vertices[i] = new_vertices[j] # i += 1 new_adj_list = create_adj_list(vertices, new_triangles) return vertices, new_triangles, new_adj_list