diff --git a/coxeter/shapes/_distance2d.py b/coxeter/shapes/_distance2d.py new file mode 100644 index 00000000..428f70e1 --- /dev/null +++ b/coxeter/shapes/_distance2d.py @@ -0,0 +1,796 @@ +# Copyright (c) 2015-2025 The Regents of the University of Michigan. +# This file is from the coxeter project, released under the BSD 3-Clause License. + +import numpy as np +import numpy.linalg as LA + +# TODO: update docstrings? + + +def point_to_edge_distance( + point: np.ndarray, vert: np.ndarray, edge_vector: np.ndarray +) -> np.ndarray: + """ + Calculates the distances between several points and several varying lines. + + n is the total number of distance calculations that are being made. For example, let's say + we have points: A, B, C & D, and edges: U, V & W, and we want to calculate the distances between: + - A & U, A & W, B & U, C & V, C & W, D & U, D & V, and D & W + n = 8 for this example, and point = [A,A,B,C,C,D,D,D] and edge_vector = [U,W,U,V,W,U,V,W] + + Args: + point (np.ndarray): the positions of the points [shape = (n, 3)] + vert (np.ndarray): positions of the points that lie on each corresponding line [shape = (n, 3)] + edge_vector (np.ndarray): the vectors that describe each line [shape = (n, 3)] + + Returns + ------- + np.ndarray: distances [shape = (n,)] + """ + edge_unit = edge_vector / np.expand_dims( + LA.norm(edge_vector, axis=1), axis=1 + ) # unit vectors of the edges + + dist = LA.norm( + ( + (vert - point) + - ( + np.expand_dims(np.sum((vert - point) * edge_unit, axis=1), axis=1) + * edge_unit + ) + ), + axis=1, + ) # distances + return dist + + +def point_to_edge_displacement( + point: np.ndarray, vert: np.ndarray, edge_vector: np.ndarray +) -> np.ndarray: + """ + Calculates the displacements between several points and several varying lines. + + n is the total number of displacement calculations that are being made. For example, let's say + we have points: A, B, C & D, and edges: U, V & W, and we want to calculate the displacements between: + - A & U, A & W, B & U, C & V, C & W, D & U, D & V, and D & W + n = 8 for this example, and point = [A,A,B,C,C,D,D,D] and edge_vector = [U,W,U,V,W,U,V,W] + + Args: + point (np.ndarray): the positions of the points [shape = (n, 3)] + vert (np.ndarray): positions of the points that lie on each corresponding line [shape = (n, 3)] + edge_vector (np.ndarray): the vectors that describe each line [shape = (n, 3)] + + Returns + ------- + np.ndarray: displacements [shape = (n, 3)] + """ + edge_unit = edge_vector / np.expand_dims( + LA.norm(edge_vector, axis=1), axis=1 + ) # unit vectors of the edges + + disp = (vert - point) - ( + np.expand_dims(np.sum((vert - point) * edge_unit, axis=1), axis=1) * edge_unit + ) # displacements + return disp + + +def point_to_face_distance( + point: np.ndarray, vert: np.ndarray, face_normal: np.ndarray +) -> np.ndarray: + """ + Calculates the distances between a single point and the plane of the polygon. + + Args: + point (np.ndarray): the positions of the points [shape=(n_points, 3)] + vert (np.ndarray): a point that lies on the plane of the polygon [shape=(3,)] + face_normal (np.ndarray): the normal that describes the plane of the polygon [shape=(3,)] + + Returns + ------- + np.ndarray: distances [shape = (n_points,)] + """ + vert_point_vect = vert - point + face_unit = face_normal / LA.norm( + face_normal + ) # unit vector of the normal of the polygon + dist = abs(vert_point_vect @ np.transpose(face_unit)) + + return dist + + +def point_to_face_displacement( + point: np.ndarray, vert: np.ndarray, face_normal: np.ndarray +) -> np.ndarray: + """ + Calculates the displacements between a single point and the plane of the polygon. + + Args: + point (np.ndarray): the positions of the points (shape=(n_points, 3)) + vert (np.ndarray): a point that lies on the plane of the polygon (shape=(3,)) + face_normal (np.ndarray): the normal that describes the plane of the polygon (shape=(3,)) + + Returns + ------- + np.ndarray: displacements (n_points, 3) + """ + vert_point_vect = vert - point + face_unit = face_normal / LA.norm( + face_normal + ) # unit vector of the normal of the polygon + disp = ( + np.expand_dims(np.sum(vert_point_vect * face_unit, axis=1), axis=1) * face_unit + ) # *(-1) + + return disp + + +def get_vert_zones(shape): + """ + Gets the constraints and bounds needed to partition the volume surrounding a polygon into zones + where the shortest distance from any point that is within a vertex zone is the distance between the + point and the corresponding vertex. + + Args: + #will be just `self` once added into coxeter + + Returns + ------- + dict: "constraint": np.ndarray [shape = (n_verts, 2, 3)], "bounds": np.ndarray [shape = (n_verts, 2)] + """ + vert_constraint = np.append( + np.expand_dims(shape.edge_vectors, axis=1), + -1 + * np.expand_dims( + np.append( + np.expand_dims(shape.edge_vectors[-1], axis=0), + shape.edge_vectors[:-1], + axis=0, + ), + axis=1, + ), + axis=1, + ) + + vert_bounds = np.sum( + vert_constraint * np.expand_dims(shape.vertices, axis=1), axis=2 + ) + + return {"constraint": vert_constraint, "bounds": vert_bounds} + + +def get_edge_zones(shape): + """ + Gets the constraints and bounds needed to partition the volume surrounding a polygon into zones + where the shortest distance from any point that is within an edge zone is the distance between the + point and the corresponding edge. + + Args: + #will be just `self` once added into coxeter + + Returns + ------- + dict: "constraint": np.ndarray [shape = (n_edges, 3, 3)], "bounds": np.ndarray [shape = (n_edges, 3)] + """ + # vectors that are 90 degrees from the edges and point inwards + edges_90 = -1 * np.expand_dims(np.cross(shape.edge_vectors, shape.normal), axis=1) + + # Calculating the constraint [shape = (n_edges, 3, 3)] + edge_constraint = np.append( + -1 * np.expand_dims(shape.edge_vectors, axis=1), + np.expand_dims(shape.edge_vectors, axis=1), + axis=1, + ) + edge_constraint = np.append(edge_constraint, edges_90, axis=1) + + # Bounds [shape = (n_edges, 3)] + edge_bounds = np.zeros((shape.num_vertices, 3)) + edge_bounds[:, 0] = np.sum(edge_constraint[:, 0] * (shape.vertices), axis=1) + edge_bounds[:, 1] = np.sum( + edge_constraint[:, 1] + * ( + np.append( + shape.vertices[1:], np.expand_dims(shape.vertices[0], axis=0), axis=0 + ) + ), + axis=1, + ) + edge_bounds[:, 2] = np.sum( + edge_constraint[:, 2] + * ( + np.append( + shape.vertices[1:], np.expand_dims(shape.vertices[0], axis=0), axis=0 + ) + ), + axis=1, + ) + + return {"constraint": edge_constraint, "bounds": edge_bounds} + + +def get_face_zones(shape): + """ + Gets the constraints and bounds needed to partition the volume surrounding a polygon into zones + where the shortest distance from any point that is within a triangulated face zone is the distance between the + point and the corresponding triangulated face. + + Args: + #will be just `self` once added into coxeter + + Returns + ------- + dict: "constraint": np.ndarray , "bounds": np.ndarray + """ + face_constraint = np.cross( + shape.edge_vectors, shape.normal + ) # only one face zone for a polygon | shape = (n_edges, 3) + face_bounds = np.sum(face_constraint * shape.vertices, axis=1) # shape = (n_edges,) + + # Checking to see if all the vertices are in the face zone. If not, the polygon is nonconvex. + if ( + np.all( + face_constraint @ np.transpose(shape.vertices) + <= np.expand_dims(face_bounds, axis=1) + 5e-6 + ) + == False + ): + # --Polygon is nonconvex and needs to be triangulated-- + triangle_verts = [] + + for tri in shape._triangulation(): + triangle_verts.append(list(tri)) + + triangle_verts = np.asarray(triangle_verts) + tri_edges = ( + np.append( + triangle_verts[:, 1:], + np.expand_dims(triangle_verts[:, 0], axis=1), + axis=1, + ) + - triangle_verts + ) # edges point counterclockwise + + face_constraint = np.cross( + tri_edges, shape.normal + ) # shape = (n_triangles, 3, 3) + face_bounds = np.sum( + face_constraint * triangle_verts, axis=2 + ) # shape = (n_triangles, 3) + + else: + # --Polygon is convex-- + face_constraint = np.expand_dims( + face_constraint, axis=0 + ) # shape = (1, n_edges, 3) + face_bounds = np.expand_dims(face_bounds, axis=0) # shape = (1, n_edges) + + return {"constraint": face_constraint, "bounds": face_bounds} + + +def shortest_distance_to_surface( + shape, + points: np.ndarray, + translation_vector: np.ndarray, +) -> np.ndarray: + """ + Solves for the shortest distance between points and the surface of a polygon. + If the point lies inside the polyhedron, the distance is negative. + + This function calculates the shortest distance by partitioning the space around + a polyhedron into zones: vertex, edge, and face. Determining the zone(s) a + point lies in, determines the distance calculation(s) done. For a vertex zone, + the distance is calculated between a point and the vertex. For an edge zone, the + distance is calculated between a point and the edge. For a face zone, the distance + is calculated between a point and the face. Zones are allowed to overlap, and points + can be in more than one zone. By taking the minimum of all the calculated distances, + the shortest distances are found. + + Args: + points (list or np.ndarray): positions of the points + translation_vector (list or np.ndarray): translation vector of the polyhedron [shape = (3,) or (2,)] + + Returns + ------- + np.ndarray: shortest distances [shape = (n_points,)] + """ + points = np.asarray(points) + translation_vector = np.asarray(translation_vector) + + if len(points.shape) == 1: + points = np.expand_dims(points, axis=0) + + n_points = len(points) # number of inputted points + + if points.shape[1] == 2: + points = np.append(points, np.zeros((n_points, 1)), axis=1) + + if ( + translation_vector.shape[0] > 3 + or len(translation_vector.shape) > 1 + or translation_vector.shape[0] < 2 + ): + raise ValueError( + f"Expected the shape of the polygon's position to be either (2,) or (3,), instead it got {translation_vector.shape}" + ) + + if translation_vector.shape[0] == 2: + translation_vector = np.append(translation_vector, [0]) + + # Updating bounds with the position of the polyhedron + vert_bounds = shape.vertex_zones["bounds"] + ( + shape.vertex_zones["constraint"] @ translation_vector + ) + edge_bounds = shape.edge_zones["bounds"] + ( + shape.edge_zones["constraint"] @ translation_vector + ) + face_bounds = shape.face_zones["bounds"] + ( + shape.face_zones["constraint"] @ translation_vector + ) + + points_trans = np.transpose(points) + + max_value = 3 * np.max( + LA.norm(points - (translation_vector + shape.vertices[0]), axis=1) + ) + + min_dist_arr = np.ones((len(points), 1)) * max_value + + # Solving for the distances between the points and any relevant vertices + vert_bool = np.all( + (shape.vertex_zones["constraint"] @ points_trans) + <= np.expand_dims(vert_bounds, axis=2), + axis=1, + ) # <--- shape = (number_of_vertex_zones, number_of_points) + if np.any(vert_bool): + # v--- shape = (number of True in vert_bool,) ---v + vert_used = np.transpose( + np.tile(np.arange(0, shape.num_vertices, 1), (n_points, 1)) + )[ + vert_bool + ] # Contains the indices of the vertices that hold True for vert_bool + v_points_used = np.tile(np.arange(0, n_points, 1), (shape.num_vertices, 1))[ + vert_bool + ] # Contains the indices of the points that hold True for vert_bool + + vert_dist = np.ones((shape.num_vertices, n_points)) * max_value + vert_dist[vert_bool] = LA.norm( + points[v_points_used] - (shape.vertices[vert_used] + translation_vector), + axis=1, + ) # Distances between two points + vert_dist = np.transpose(vert_dist) # <--- shape = (n_points, n_verts) + + vert_dist_arg = np.expand_dims(np.argmin(abs(vert_dist), axis=1), axis=1) + vert_dist = np.take_along_axis(vert_dist, vert_dist_arg, axis=1) + + min_dist_arr = np.concatenate((min_dist_arr, vert_dist), axis=1) + + # Solving for the distances between the points and any relevant edges + edge_bool = np.all( + (shape.edge_zones["constraint"] @ points_trans) + <= np.expand_dims(edge_bounds, axis=2), + axis=1, + ) # <--- shape = (number_of_edge_zones, number_of_points) + if np.any(edge_bool): + # v--- shape = (number of True in edge_bool,) ---v + edge_used = np.transpose( + np.tile(np.arange(0, shape.num_vertices, 1), (n_points, 1)) + )[edge_bool] # Contains the indices of the edges that hold True for edge_bool + e_points_used = np.tile(np.arange(0, n_points, 1), (shape.num_vertices, 1))[ + edge_bool + ] # Contains the indices of the points that hold True for edge_bool + + vert_on_edge = ( + shape.vertices[shape.edges[edge_used][:, 0]] + translation_vector + ) # Vertices that lie on the needed edges + edge_vectors = ( + np.append( + shape.vertices[1:], np.expand_dims(shape.vertices[0], axis=0), axis=0 + ) + - shape.vertices + ) + + edge_dist = np.ones((shape.num_vertices, n_points)) * max_value + edge_dist[edge_bool] = point_to_edge_distance( + points[e_points_used], vert_on_edge, edge_vectors[edge_used] + ) # Distances between a point and a line + edge_dist = np.transpose(edge_dist) # <--- shape = (n_points, n_edges) + + edge_dist_arg = np.expand_dims(np.argmin(abs(edge_dist), axis=1), axis=1) + edge_dist = np.take_along_axis(edge_dist, edge_dist_arg, axis=1) + + min_dist_arr = np.concatenate((min_dist_arr, edge_dist), axis=1) + + face_bool = np.all( + (shape.face_zones["constraint"] @ points_trans) + <= np.expand_dims(face_bounds, axis=2), + axis=1, + ) # <--- shape = (number_of_face_zones, number_of_points) + if np.any(face_bool): + vert_on_face = shape.vertices[0] + translation_vector + face_dist = point_to_face_distance(points, vert_on_face, shape.normal) + face_dist = face_dist + max_value * (np.any(face_bool, axis=0) == False).astype( + int + ) + face_dist = np.expand_dims(face_dist, axis=1) + + min_dist_arr = np.concatenate((min_dist_arr, face_dist), axis=1) + + true_min_dist = np.min(min_dist_arr, axis=1) + + return true_min_dist + + +def shortest_displacement_to_surface( + shape, + points: np.ndarray, + translation_vector: np.ndarray, +) -> np.ndarray: + """ + Solves for the shortest displacement between points and the surface of a polyhedron. + + This function calculates the shortest displacement by partitioning the space around + a polyhedron into zones: vertex, edge, and face. Determining the zone(s) a + point lies in, determines the displacement calculation(s) done. For a vertex zone, + the displacement is calculated between a point and the vertex. For an edge zone, the + displacement is calculated between a point and the edge. For a face zone, the + displacement is calculated between a point and the face. Zones are allowed to overlap, + and points can be in more than one zone. By taking the minimum of all the distances of + the calculated displacements, the shortest displacements are found. + + Args: + points (list or np.ndarray): positions of the points + translation_vector (list or np.ndarray): translation vector of the polyhedron [shape = (3,) or (2,)] + + Returns + ------- + np.ndarray: shortest displacements [shape = (n_points, 3)] + """ + points = np.asarray(points) + translation_vector = np.asarray(translation_vector) + + if points.shape == (3,) or points.shape == (2,): + points = np.expand_dims(points, axis=0) + + n_points = len(points) # number of inputted points + n_verts = shape.num_vertices # number of vertices = number of vertex zones + n_edges = n_verts # number of edges = number of edge zones + + if points.shape[1] == 2: + points = np.append(points, np.zeros((n_points, 1)), axis=1) + + if ( + translation_vector.shape[0] > 3 + or len(translation_vector.shape) > 1 + or translation_vector.shape[0] < 2 + ): + raise ValueError( + f"Expected the shape of the polygon's position to be either (2,) or (3,), instead it got {translation_vector.shape}" + ) + + if translation_vector.shape[0] == 2: + translation_vector = np.append(translation_vector, [0]) + + # Updating bounds with the position of the polyhedron + vert_bounds = shape.vertex_zones["bounds"] + ( + shape.vertex_zones["constraint"] @ translation_vector + ) + edge_bounds = shape.edge_zones["bounds"] + ( + shape.edge_zones["constraint"] @ translation_vector + ) + face_bounds = shape.face_zones["bounds"] + ( + shape.face_zones["constraint"] @ translation_vector + ) + + points_trans = np.transpose(points) + + max_value = 3 * np.max( + LA.norm(points - (translation_vector + shape.vertices[0]), axis=1) + ) + + min_disp_arr = np.ones((n_points, 1, 3)) * max_value + + # Solving for the distances between the points and any relevant vertices + vert_bool = np.all( + (shape.vertex_zones["constraint"] @ points_trans) + <= np.expand_dims(vert_bounds, axis=2), + axis=1, + ) # <--- shape = (number_of_vertex_zones, number_of_points) + if np.any(vert_bool): + # v--- shape = (number of True in vert_bool,) ---v + vert_used = np.transpose(np.tile(np.arange(0, n_verts, 1), (n_points, 1)))[ + vert_bool + ] # Contains the indices of the vertices that hold True for vert_bool + v_points_used = np.tile(np.arange(0, n_points, 1), (n_verts, 1))[ + vert_bool + ] # Contains the indices of the points that hold True for vert_bool + + vert_disp = np.ones((n_verts, n_points, 3)) * max_value + vert_disp[vert_bool] = ( + shape.vertices[vert_used] + translation_vector + ) - points[v_points_used] # Displacements between two points + vert_disp = np.transpose( + vert_disp, (1, 0, 2) + ) # <--- shape = (n_points, n_verts, 3) + + vert_disp_min = np.expand_dims( + np.argmin(LA.norm(vert_disp, axis=2), axis=1), axis=(1, 2) + ) + vert_disp = np.take_along_axis(vert_disp, vert_disp_min, axis=1) + + min_disp_arr = np.concatenate((min_disp_arr, vert_disp), axis=1) + + # Solving for the distances between the points and any relevant edges + edge_bool = np.all( + (shape.edge_zones["constraint"] @ points_trans) + <= np.expand_dims(edge_bounds, axis=2), + axis=1, + ) # <--- shape = (number_of_edge_zones, number_of_points) + if np.any(edge_bool): + # v--- shape = (number of True in edge_bool,) ---v + edge_used = np.transpose(np.tile(np.arange(0, n_edges, 1), (n_points, 1)))[ + edge_bool + ] # Contains the indices of the edges that hold True for edge_bool + e_points_used = np.tile(np.arange(0, n_points, 1), (n_edges, 1))[ + edge_bool + ] # Contains the indices of the points that hold True for edge_bool + + vert_on_edge = ( + shape.vertices[shape.edges[edge_used][:, 0]] + translation_vector + ) # Vertices that lie on the needed edges + edge_vectors = ( + np.append( + shape.vertices[1:], np.expand_dims(shape.vertices[0], axis=0), axis=0 + ) + - shape.vertices + ) + + edge_disp = np.ones((n_edges, n_points, 3)) * max_value + edge_disp[edge_bool] = point_to_edge_displacement( + points[e_points_used], vert_on_edge, edge_vectors[edge_used] + ) # Displacements between a point and a line + edge_disp = np.transpose( + edge_disp, (1, 0, 2) + ) # <--- shape = (n_points, n_edges, 3) + + edge_disp_arg = np.expand_dims( + np.argmin(LA.norm(edge_disp, axis=2), axis=1), axis=(1, 2) + ) + edge_disp = np.take_along_axis(edge_disp, edge_disp_arg, axis=1) + + min_disp_arr = np.concatenate((min_disp_arr, edge_disp), axis=1) + + face_bool = np.all( + (shape.face_zones["constraint"] @ points_trans) + <= np.expand_dims(face_bounds, axis=2), + axis=1, + ) # <--- shape = (number_of_face_zones, number_of_points) + if np.any(face_bool): + face_disp = point_to_face_displacement( + points, shape.vertices[0] + translation_vector, shape.normal + ) + np.repeat( + np.expand_dims( + (max_value * (np.any(face_bool, axis=0) == False).astype(int)), axis=1 + ), + 3, + axis=1, + ) + + min_disp_arr = np.concatenate( + (min_disp_arr, np.expand_dims(face_disp, axis=1)), axis=1 + ) + + disp_list_bool = np.argmin((LA.norm(min_disp_arr, axis=2)), axis=1).reshape( + n_points, 1, 1 + ) + true_min_disp = np.squeeze( + np.take_along_axis(min_disp_arr, disp_list_bool, axis=1), axis=1 + ) + + return true_min_disp + + +def spheropolygon_shortest_displacement_to_surface( + shape, + radius, + points: np.ndarray, + translation_vector: np.ndarray, +) -> np.ndarray: + """ + Solves for the shortest displacement between points and the surface of a polyhedron. + + This function calculates the shortest displacement by partitioning the space around + a polyhedron into zones: vertex, edge, and face. Determining the zone(s) a + point lies in, determines the displacement calculation(s) done. For a vertex zone, + the displacement is calculated between a point and the vertex. For an edge zone, the + displacement is calculated between a point and the edge. For a face zone, the + displacement is calculated between a point and the face. Zones are allowed to overlap, + and points can be in more than one zone. By taking the minimum of all the distances of + the calculated displacements, the shortest displacements are found. + + Args: + points (list or np.ndarray): positions of the points + translation_vector (list or np.ndarray): translation vector of the polyhedron [shape = (3,) or (2,)] + + Returns + ------- + np.ndarray: shortest displacements [shape = (n_points, 3)] + """ + points = np.asarray(points) + translation_vector = np.asarray(translation_vector) + + if points.shape == (3,) or points.shape == (2,): + points = np.expand_dims(points, axis=0) + + n_points = len(points) # number of inputted points + n_verts = shape.num_vertices # number of vertices = number of vertex zones + n_edges = n_verts # number of edges = number of edge zones + + if points.shape[1] == 2: + points = np.append(points, np.zeros((n_points, 1)), axis=1) + + if ( + translation_vector.shape[0] > 3 + or len(translation_vector.shape) > 1 + or translation_vector.shape[0] < 2 + ): + raise ValueError( + f"Expected the shape of the polygon's position to be either (2,) or (3,), instead it got {translation_vector.shape}" + ) + + if translation_vector.shape[0] == 2: + translation_vector = np.append(translation_vector, [0]) + + # Updating bounds with the position of the polyhedron + vert_bounds = shape.vertex_zones["bounds"] + ( + shape.vertex_zones["constraint"] @ translation_vector + ) + edge_bounds = shape.edge_zones["bounds"] + ( + shape.edge_zones["constraint"] @ translation_vector + ) + face_bounds = shape.face_zones["bounds"] + ( + shape.face_zones["constraint"] @ translation_vector + ) + + points_trans = np.transpose(points) + + max_value = 3 * np.max( + LA.norm(points - (translation_vector + shape.vertices[0]), axis=1) + ) + + min_disp_arr = np.ones((n_points, 1, 3)) * max_value + + # Solving for the distances between the points and any relevant vertices + vert_bool = np.all( + (shape.vertex_zones["constraint"] @ points_trans) + <= np.expand_dims(vert_bounds, axis=2), + axis=1, + ) # <--- shape = (number_of_vertex_zones, number_of_points) + if np.any(vert_bool): + # v--- shape = (number of True in vert_bool,) ---v + vert_used = np.transpose(np.tile(np.arange(0, n_verts, 1), (n_points, 1)))[ + vert_bool + ] # Contains the indices of the vertices that hold True for vert_bool + v_points_used = np.tile(np.arange(0, n_points, 1), (n_verts, 1))[ + vert_bool + ] # Contains the indices of the points that hold True for vert_bool + + vert_disp = np.ones((n_verts, n_points, 3)) * max_value + vert_disp[vert_bool] = ( + shape.vertices[vert_used] + translation_vector + ) - points[v_points_used] # Displacements between two points + vert_disp = np.transpose( + vert_disp, (1, 0, 2) + ) # <--- shape = (n_points, n_verts, 3) + + vert_disp_min = np.expand_dims( + np.argmin(LA.norm(vert_disp, axis=2), axis=1), axis=(1, 2) + ) + vert_disp = np.take_along_axis(vert_disp, vert_disp_min, axis=1) + + # for spheropolygon + vert_projection = np.squeeze( + vert_disp + - np.expand_dims( + vert_disp @ (shape.normal / np.linalg.norm(shape.normal)), axis=1 + ) + * (shape.normal / np.linalg.norm(shape.normal)), + axis=1, + ) + v_projection_bool = np.linalg.norm(vert_projection, axis=1) > radius + vert_projection[v_projection_bool] = ( + radius + * vert_projection[v_projection_bool] + / np.expand_dims( + np.linalg.norm(vert_projection[v_projection_bool], axis=1), axis=1 + ) + ) + vert_disp = vert_disp - vert_projection + + min_disp_arr = np.concatenate((min_disp_arr, vert_disp), axis=1) + + # Solving for the distances between the points and any relevant edges + edge_bool = np.all( + (shape.edge_zones["constraint"] @ points_trans) + <= np.expand_dims(edge_bounds, axis=2), + axis=1, + ) # <--- shape = (number_of_edge_zones, number_of_points) + if np.any(edge_bool): + # v--- shape = (number of True in edge_bool,) ---v + edge_used = np.transpose(np.tile(np.arange(0, n_edges, 1), (n_points, 1)))[ + edge_bool + ] # Contains the indices of the edges that hold True for edge_bool + e_points_used = np.tile(np.arange(0, n_points, 1), (n_edges, 1))[ + edge_bool + ] # Contains the indices of the points that hold True for edge_bool + + vert_on_edge = ( + shape.vertices[shape.edges[edge_used][:, 0]] + translation_vector + ) # Vertices that lie on the needed edges + edge_vectors = ( + np.append( + shape.vertices[1:], np.expand_dims(shape.vertices[0], axis=0), axis=0 + ) + - shape.vertices + ) + + edge_disp = np.ones((n_edges, n_points, 3)) * max_value + edge_disp[edge_bool] = point_to_edge_displacement( + points[e_points_used], vert_on_edge, edge_vectors[edge_used] + ) # Displacements between a point and a line + edge_disp = np.transpose( + edge_disp, (1, 0, 2) + ) # <--- shape = (n_points, n_edges, 3) + + edge_disp_arg = np.expand_dims( + np.argmin(LA.norm(edge_disp, axis=2), axis=1), axis=(1, 2) + ) + edge_disp = np.take_along_axis(edge_disp, edge_disp_arg, axis=1) + + # for spheropolygon + edge_projection = np.squeeze( + edge_disp + - np.expand_dims( + edge_disp @ (shape.normal / np.linalg.norm(shape.normal)), axis=1 + ) + * (shape.normal / np.linalg.norm(shape.normal)), + axis=1, + ) + e_projection_bool = np.linalg.norm(edge_projection, axis=1) > radius + edge_projection[e_projection_bool] = ( + radius + * edge_projection[e_projection_bool] + / np.expand_dims( + np.linalg.norm(edge_projection[e_projection_bool], axis=1), axis=1 + ) + ) + edge_disp = edge_disp - edge_projection + + min_disp_arr = np.concatenate((min_disp_arr, edge_disp), axis=1) + + face_bool = np.all( + (shape.face_zones["constraint"] @ points_trans) + <= np.expand_dims(face_bounds, axis=2), + axis=1, + ) # <--- shape = (number_of_face_zones, number_of_points) + if np.any(face_bool): + face_disp = point_to_face_displacement( + points, shape.vertices[0] + translation_vector, shape.normal + ) + np.repeat( + np.expand_dims( + (max_value * (np.any(face_bool, axis=0) == False).astype(int)), axis=1 + ), + 3, + axis=1, + ) + + min_disp_arr = np.concatenate( + (min_disp_arr, np.expand_dims(face_disp, axis=1)), axis=1 + ) + + disp_list_bool = np.argmin((LA.norm(min_disp_arr, axis=2)), axis=1).reshape( + n_points, 1, 1 + ) + true_min_disp = np.squeeze( + np.take_along_axis(min_disp_arr, disp_list_bool, axis=1), axis=1 + ) + + return true_min_disp diff --git a/coxeter/shapes/_distance3d.py b/coxeter/shapes/_distance3d.py new file mode 100644 index 00000000..9df0bb09 --- /dev/null +++ b/coxeter/shapes/_distance3d.py @@ -0,0 +1,1127 @@ +# Copyright (c) 2015-2025 The Regents of the University of Michigan. +# This file is from the coxeter project, released under the BSD 3-Clause License. + +import numpy as np +import numpy.linalg as LA + +# TODO: update docstrings? + + +def get_edge_face_neighbors(shape) -> np.ndarray: + """ + Gets the indices of the faces that are adjacent to each edge. + + Args: + shape (Polyhedron): the polyhedron that is being looked at (can be convex or concave) + + Returns + ------- + np.ndarray: the indices of the nearest faces for each edge [shape = (n_edges, 2)] + """ + faces_len = shape.num_faces + num_edges = shape.num_edges + + # appending the vertex list of each face and a -1 to the end of each list of vertices, then flattening the awkward array (the -1 indicates a change in face) + faces_flat = [] + + for face in shape.faces: + face_extra = np.array([face[0], -1]) + + faces_flat = np.append(faces_flat, face) + faces_flat = np.append(faces_flat, face_extra) + + faces_flat = np.asarray(faces_flat) + + # creating a matrix where each row corresponds to an edge that contains the indices of its two corresponding vertices (a -1 index indicates a change in face) + list_len = len(faces_flat) + face_edge_mat = np.block( + [ + np.expand_dims(faces_flat[:-1], axis=1), + np.expand_dims(faces_flat[1:], axis=1), + ] + ) + + # finding the number of edges associated with each face + fe_mat_inds = np.arange(0, list_len - 1, 1) + find_num_edges = fe_mat_inds[ + (fe_mat_inds == 0) + (np.any(face_edge_mat == -1, axis=1)) + ] + find_num_edges[:][0] = -1 + find_num_edges = find_num_edges.reshape(faces_len, 2) + face_num_edges = find_num_edges[:, 1] - find_num_edges[:, 0] - 1 + + # repeating each face index for the number of edges that are associated with it; length equals num_edges * 2 + face_correspond_inds = np.repeat(np.arange(0, faces_len, 1), face_num_edges) + + # shape.edges lists the indices of the edge vertices lowest to highest. edges1_reshape lists the indices of the edge vertices highest to lowest + edges1_reshape = np.hstack( + ( + np.expand_dims(shape.edges[:, 1], axis=1), + np.expand_dims(shape.edges[:, 0], axis=1), + ) + ) + + # For the new_edge_ind_bool: rows correspond with the face_correspond_inds and columns correpond with the edge index; finding the neighboring faces for each edge + true_face_edge_mat = np.tile( + np.expand_dims(face_edge_mat[np.all(face_edge_mat != -1, axis=1)], axis=1), + (1, num_edges, 1), + ) + new_edge_ind_bool0 = np.all( + true_face_edge_mat == np.expand_dims(shape.edges, axis=0), axis=2 + ).astype(int) # faces to the LEFT of edges if edges are oriented pointing upwards + new_edge_ind_bool1 = np.all( + true_face_edge_mat == np.expand_dims(edges1_reshape, axis=0), axis=2 + ).astype(int) # faces to the RIGHT of edges if edges are oriented pointing upwards + + # tiling face_correspond_inds so it can be multiplied to the new_edge_ind_bool0 and new_edge_ind_bool1 + new_face_corr_inds = np.tile( + np.expand_dims(face_correspond_inds, axis=1), (1, num_edges) + ) + + # getting the face indices and completing the edge-face neighbors + ef_neighbor0 = np.expand_dims( + np.sum(new_face_corr_inds * new_edge_ind_bool0, axis=0), axis=1 + ) # faces to the LEFT + ef_neighbor1 = np.expand_dims( + np.sum(new_face_corr_inds * new_edge_ind_bool1, axis=0), axis=1 + ) # faces to the RIGHT + ef_neighbor = np.hstack((ef_neighbor0, ef_neighbor1)) + + return ef_neighbor + + +def point_to_edge_distance( + point: np.ndarray, vert: np.ndarray, edge_vector: np.ndarray +) -> np.ndarray: + """ + Calculates the distances between several points and several varying lines. + + n is the total number of distance calculations that are being made. For example, let's say + we have points: A, B, C & D, and edges: U, V & W, and we want to calculate the distances between: + - A & U, A & W, B & U, C & V, C & W, D & U, D & V, and D & W + n = 8 for this example, and point = [A,A,B,C,C,D,D,D] and edge_vector = [U,W,U,V,W,U,V,W] + + Args: + point (np.ndarray): the positions of the points [shape = (n, 3)] + vert (np.ndarray): positions of the points that lie on each corresponding line [shape = (n, 3)] + edge_vector (np.ndarray): the vectors that describe each line [shape = (n, 3)] + + Returns + ------- + np.ndarray: distances [shape = (n,)] + """ + edge_unit = edge_vector / np.expand_dims( + LA.norm(edge_vector, axis=1), axis=1 + ) # unit vectors of the edges + + dist = LA.norm( + ( + (vert - point) + - ( + np.expand_dims(np.sum((vert - point) * edge_unit, axis=1), axis=1) + * edge_unit + ) + ), + axis=1, + ) # distances + return dist + + +def point_to_edge_displacement( + point: np.ndarray, vert: np.ndarray, edge_vector: np.ndarray +) -> np.ndarray: + """ + Calculates the displacements between several points and several varying lines. + + n is the total number of displacement calculations that are being made. For example, let's say + we have points: A, B, C & D, and edges: U, V & W, and we want to calculate the displacements between: + - A & U, A & W, B & U, C & V, C & W, D & U, D & V, and D & W + n = 8 for this example, and point = [A,A,B,C,C,D,D,D] and edge_vector = [U,W,U,V,W,U,V,W] + + Args: + point (np.ndarray): the positions of the points [shape = (n, 3)] + vert (np.ndarray): positions of the points that lie on each corresponding line [shape = (n, 3)] + edge_vector (np.ndarray): the vectors that describe each line [shape = (n, 3)] + + Returns + ------- + np.ndarray: displacements [shape = (n, 3)] + """ + edge_unit = edge_vector / np.expand_dims( + LA.norm(edge_vector, axis=1), axis=1 + ) # unit vectors of the edges + + disp = (vert - point) - ( + np.expand_dims(np.sum((vert - point) * edge_unit, axis=1), axis=1) * edge_unit + ) # displacements + return disp + + +def point_to_face_distance( + point: np.ndarray, vert: np.ndarray, face_normal: np.ndarray +) -> np.ndarray: + """ + Calculates the distances between several points and several varying planes. + + n is the total number of distance calculations that are being made. For example, let's say + we have points: A, B, C & D, and faces: P, Q & R, and we want to calculate the distances between: + - A & P, A & R, B & P, C & Q, C & R, D & P, D & Q, and D & R + n = 8 for this example, and point = [A,A,B,C,C,D,D,D] and edge_vector = [P,R,P,Q,R,P,Q.R] + + Args: + point (np.ndarray): the positions of the points [shape = (n, 3)] + vert (np.ndarray): points that lie on each corresponding plane [shape = (n, 3)] + face_normal (np.ndarray): the normals that describe each plane [shape = (n, 3)] + + Returns + ------- + np.ndarray: distances [shape = (n,)] + """ + vert_point_vect = ( + vert - point + ) # displacements between the points and relevent vertices + face_unit = face_normal / np.expand_dims( + LA.norm(face_normal, axis=1), axis=1 + ) # unit vectors of the normals of the faces + + dist = np.sum(vert_point_vect * face_unit, axis=1) * (-1) # distances + + return dist + + +def point_to_face_displacement( + point: np.ndarray, vert: np.ndarray, face_normal: np.ndarray +) -> np.ndarray: + """ + Calculates the displacements between several points and several varying planes. + + n is the total number of displacement calculations that are being made. For example, let's say + we have points: A, B, C & D, and faces: P, Q & R, and we want to calculate the displacements between: + - A & P, A & R, B & P, C & Q, C & R, D & P, D & Q, and D & R + n = 8 for this example, and point = [A,A,B,C,C,D,D,D] and edge_vector = [P,R,P,Q,R,P,Q.R] + + Args: + point (np.ndarray): the positions of the points [shape = (n, 3)] + vert (np.ndarray): points that lie on each corresponding plane [shape = (n, 3)] + face_normal (np.ndarray): the normals that describe each plane [shape = (n, 3)] + + Returns + ------- + np.ndarray: displacements [shape = (n, 3)] + """ + vert_point_vect = ( + vert - point + ) # displacements between the points and relevent vertices + face_units = face_normal / np.expand_dims( + LA.norm(face_normal, axis=1), axis=1 + ) # unit vectors of the normals of the faces + + disp = ( + np.expand_dims(np.sum(vert_point_vect * face_units, axis=1), axis=1) + * face_units + ) # *(-1) #displacements + + return disp + + +def get_vert_zones(shape): + """ + Gets the constraints and bounds needed to partition the volume surrounding a polyhedron into zones + where the shortest distance from any point that is within a vertex zone is the distance between the + point and the corresponding vertex. + + Args: + shape (Polyhedron): the polyhedron that is being looked at (can be convex or concave) + + Returns + ------- + dict: "constraint": np.ndarray [shape = (n_verts, n_edges, 3)], "bounds": np.ndarray [shape = (n_verts, n_edges)] + """ + # For a generalized shape, we cannot assume that every vertex has the same number of edges connected to it + # (EX:vertices in a cube have 3 connected edges each, and for an icosahedron, vertices have 5 conncected edges). + # This would result in a ragged list for the constraint and bounds, which is not ideal. + + # v--- This `for`` loop is used to build and fill that ragged list with zeros, so that it makes an easy to work with array. ---v + for v_i in range(len(shape.vertices)): + pos_adj_edges = shape.edge_vectors[ + shape.edges[:, 0] == v_i + ] # edges that point away from v_i + neg_adj_edges = (-1) * shape.edge_vectors[ + shape.edges[:, 1] == v_i + ] # edges that point towards v_i, so have to multiply by -1 + adjacent_edges = np.append(pos_adj_edges, neg_adj_edges, axis=0) + + if v_i == 0: # initial vertex, start of building the constraint + vert_constraint = np.asarray([adjacent_edges]) + + else: + difference = len(adjacent_edges) - vert_constraint.shape[1] + # ^---difference between the # of edges v_i has and the max # of edges from a previous vertex + + if ( + difference < 0 + ): # adjacent_edges needs to be filled with zeros to match the length of axis=1 for vert_constraint + adjacent_edges = np.append( + adjacent_edges, np.zeros((abs(difference), 3)), axis=0 + ) + + if ( + difference > 0 + ): # vert_constraint needs to be filled with zeros to match the # of edges for v_i + vert_constraint = np.append( + vert_constraint, + np.zeros((len(vert_constraint), difference, 3)), + axis=1, + ) + + vert_constraint = np.append( + vert_constraint, np.asarray([adjacent_edges]), axis=0 + ) + + vert_bounds = np.sum( + vert_constraint * np.expand_dims(shape.vertices, axis=1), axis=2 + ) + + return {"constraint": vert_constraint, "bounds": vert_bounds} + + +def get_edge_zones(shape): + """ + Gets the constraints and bounds needed to partition the volume surrounding a polyhedron into zones + where the shortest distance from any point that is within an edge zone is the distance between the + point and the corresponding edge. + + Args: + shape (Polyhedron): the polyhedron that is being looked at (can be convex or concave) + + Returns + ------- + dict: "constraint": np.ndarray [shape = (n_edges, 4, 3)], "bounds": np.ndarray [shape = (n_edges, 4)] + """ + # Set up + edge_constraint = np.zeros((shape.num_edges, 4, 3)) + edge_bounds = np.zeros((shape.num_edges, 4)) + + # Calculating the normals of the plane boundaries + edge_constraint[:, 0] = shape.edge_vectors + edge_constraint[:, 1] = -1 * shape.edge_vectors + edge_constraint[:, 2] = np.cross( + shape.edge_vectors, shape.normals[shape.edge_face_neighbors[:, 1]] + ) + edge_constraint[:, 3] = -1 * np.cross( + shape.edge_vectors, shape.normals[shape.edge_face_neighbors[:, 0]] + ) + # Constraint shape = (n_edges, 4, 3) + + # Bounds [shape = (n_edges, 4)] + edge_verts = np.zeros((shape.num_edges, 2, 3)) + edge_verts[:, 0] = shape.vertices[shape.edges[:, 0]] + edge_verts[:, 1] = shape.vertices[shape.edges[:, 1]] + + edge_bounds[:, 0] = np.sum(edge_constraint[:, 0] * (edge_verts[:, 1]), axis=1) + edge_bounds[:, 1] = np.sum(edge_constraint[:, 1] * (edge_verts[:, 0]), axis=1) + edge_bounds[:, 2] = np.sum(edge_constraint[:, 2] * (edge_verts[:, 0]), axis=1) + edge_bounds[:, 3] = np.sum(edge_constraint[:, 3] * (edge_verts[:, 0]), axis=1) + + return {"constraint": edge_constraint, "bounds": edge_bounds} + + +def get_face_zones(shape): + """ + Gets the constraints and bounds needed to partition the volume surrounding a polyhedron into zones + where the shortest distance from any point that is within a triangulated face zone is the distance between the + point and the corresponding triangulated face. + + Args: + shape (Polyhedron): the polyhedron that is being looked at (can be convex or concave) + + Returns + ------- + dict: "constraint": np.ndarray [shape = (n_tri_faces, 3, 3)], "bounds": np.ndarray [shape = (n_tri_faces, 3)], + "face_points": np.ndarray [shape= (n_tri_faces, 3)], "normals": np.ndarray [shape=(n_tri_faces, 3)] + """ + # ----- Triangulating the surface of the shape ----- + try: + # checking to see if faces are already triangulated + something = np.asarray(shape.faces).reshape(shape.num_faces, 3) + + except: + # triangulating faces + triangle_verts = [] + for triangle in shape._surface_triangulation(): + triangle_verts.append(list(triangle)) + + triangle_verts = np.asarray( + triangle_verts + ) # vertices of the triangulated faces + tri_edges = ( + np.append( + triangle_verts[:, 1:], + triangle_verts[:, 0].reshape(len(triangle_verts), 1, 3), + axis=1, + ) + - triangle_verts + ) # edges point counterclockwise + tri_face_normals = np.cross( + tri_edges[:, 0], tri_edges[:, 1] + ) # normals of the triangulated faces + + else: + triangle_verts = shape.vertices[ + shape.faces + ] # vertices of the triangulated faces + tri_edges = ( + np.append( + triangle_verts[:, 1:], + triangle_verts[:, 0].reshape(len(triangle_verts), 1, 3), + axis=1, + ) + - triangle_verts + ) # edges point counterclockwise + tri_face_normals = shape.normals # normals of the triangulated faces + + face_constraint = np.cross( + tri_edges, np.expand_dims(tri_face_normals, axis=1) + ) # shape = (n_tri_faces, 3, 3) + face_bounds = np.sum( + face_constraint * triangle_verts, axis=2 + ) # shape = (n_tri_faces, 3) + + face_one_vertex = triangle_verts[ + :, 0 + ] # a point (vertex) that lies on each of the planes of the triangulated faces + + return { + "constraint": face_constraint, + "bounds": face_bounds, + "face_points": face_one_vertex, + "normals": tri_face_normals, + } + + +def get_edge_normals(shape) -> np.ndarray: + """ + Gets the analogous normals of the edges of the polyhedron. The normals point outwards from the polyhedron + and are used to determine whether an edge zone is outside or inside the polyhedron. + + Args: + shape (Polyhedron): the polyhedron that is being looked at (can be convex or concave) + + Returns + ------- + np.ndarray: analogous edge normals [shape = (n_edges, 3)] + """ + face_unit = shape.normals / np.expand_dims( + LA.norm(shape.normals, axis=1), axis=1 + ) # unit vectors of the face normals + face_unit1 = face_unit[shape.edge_face_neighbors[:, 0]] + face_unit2 = face_unit[shape.edge_face_neighbors[:, 1]] + + edge_normals = ( + face_unit1 + face_unit2 + ) # sum of the adjacent face normals for each edge + + # returning the unit vectors of the edge normals + return edge_normals / np.expand_dims(LA.norm(edge_normals, axis=1), axis=1) + + +def get_vert_normals(shape) -> np.ndarray: + """ + Gets the analogous normals of the vertices of the polyhedron. The normals point outwards from the polyhedron + and are used to determine whether a vertex zone is outside or inside the polyhedron. + + Args: + shape (Polyhedron): the polyhedron that is being looked at (can be convex or concave) + + Returns + ------- + np.ndarray: analogous vertex normals [shape = (n_verts, 3)] + """ + n_edges = len(shape.edge_normals) + n_verts = np.max(shape.edges) + 1 + + # Tiling for set up + nverts_edge_vert0 = np.tile(shape.edges[:, 0], (n_verts, 1)) + nverts_edge_vert1 = np.tile(shape.edges[:, 1], (n_verts, 1)) + vert_inds = np.arange(0, n_verts, 1).reshape((n_verts, 1)) + nverts_tile_edges = np.tile(shape.edge_normals, (n_verts, 1)).reshape( + (n_verts, n_edges, 3) + ) + + # Creating the bools needed to get the edges that correspond to each vertex + evbool0 = (np.expand_dims(nverts_edge_vert0 == vert_inds, axis=2)).astype(int) + evbool1 = (np.expand_dims(nverts_edge_vert1 == vert_inds, axis=2)).astype(int) + + # Applying the bools to find the corresponding edges + vert_edges = nverts_tile_edges * evbool0 + nverts_tile_edges * evbool1 + + vert_normals = np.sum( + vert_edges, axis=1 + ) # sum of the adjacent edge normals for each vertex + + # returning the unit vectors of the vertex normals + return vert_normals / np.expand_dims(LA.norm(vert_normals, axis=1), axis=1) + + +def get_weighted_edge_normals(shape) -> np.ndarray: + """ + Gets the weighted normals of the edges of the polyhedron. The normals point outwards from the polyhedron. + + Args: + shape (Polyhedron): the polyhedron that is being looked at (can be convex or concave) + + Returns + ------- + np.ndarray: analogous edge normals [shape = (n_edges, 3)] + """ + face_1 = shape.normals[shape.edge_face_neighbors[:, 0]] + face_2 = shape.normals[shape.edge_face_neighbors[:, 1]] + + edge_normals = face_1 + face_2 # sum of the adjacent face normals for each edge + + return edge_normals + + +def get_weighted_vert_normals(shape) -> np.ndarray: + """ + Gets the weighted normals of the vertices of the polyhedron. The normals point outwards from the polyhedron. + + Args: + shape (Polyhedron): the polyhedron that is being looked at (can be convex or concave) + + Returns + ------- + np.ndarray: analogous vertex normals [shape = (n_verts, 3)] + """ + n_edges = len(shape.edge_normals) + n_verts = np.max(shape.edges) + 1 + + # Tiling for set up + nverts_edge_vert0 = np.tile(shape.edges[:, 0], (n_verts, 1)) + nverts_edge_vert1 = np.tile(shape.edges[:, 1], (n_verts, 1)) + vert_inds = np.arange(0, n_verts, 1).reshape((n_verts, 1)) + nverts_tile_edges = np.tile(shape.weighted_edge_normals, (n_verts, 1)).reshape( + (n_verts, n_edges, 3) + ) + + # Creating the bools needed to get the edges that correspond to each vertex + evbool0 = (np.expand_dims(nverts_edge_vert0 == vert_inds, axis=2)).astype(int) + evbool1 = (np.expand_dims(nverts_edge_vert1 == vert_inds, axis=2)).astype(int) + + # Applying the bools to find the corresponding edges + vert_edges = nverts_tile_edges * evbool0 + nverts_tile_edges * evbool1 + + vert_normals = np.sum( + vert_edges, axis=1 + ) # sum of the adjacent weighted edge normals for each vertex + + return vert_normals + + +def shortest_distance_to_surface( + shp, + points: np.ndarray, + translation_vector: np.ndarray, +) -> np.ndarray: + """ + Solves for the shortest distance between points and the surface of a polyhedron. + If the point lies inside the polyhedron, the distance is negative. + + This function calculates the shortest distance by partitioning the space around + a polyhedron into zones: vertex, edge, and face. Determining the zone(s) a + point lies in, determines the distance calculation(s) done. For a vertex zone, + the distance is calculated between a point and the vertex. For an edge zone, the + distance is calculated between a point and the edge. For a face zone, the distance + is calculated between a point and the face. Zones are allowed to overlap, and points + can be in more than one zone. By taking the minimum of all the calculated distances, + the shortest distances are found. + + Args: + points (list or np.ndarray): positions of the points [shape = (n_points, 3)] + translation_vector (list or np.ndarray): translation vector of the polyhedron [shape = (3,)] + + Returns + ------- + np.ndarray: shortest distances [shape = (n_points,)] + """ + points = np.asarray(points) + translation_vector = np.asarray(translation_vector) + + if translation_vector.shape[0] != 3 or len(translation_vector.shape) > 1: + raise ValueError( + f"Expected the shape of the polygon's position to be (3,), instead it got {translation_vector.shape}" + ) + + if points.shape == (3,): + points = points.reshape(1, 3) + + atol = 1e-8 + n_points = len(points) # number of inputted points + n_verts = len(shp.vertices) # number of vertices = number of vertex zones + n_edges = len(shp.edges) # number of edges = number of edge zones + n_tri_faces = len( + shp.face_zones["bounds"] + ) # number of triangulated faces = number of triangulated face zones + + # arrays consisting of 1 or -1, and used to determine if a point is inside the polyhedron + vert_inside_mult = ( + np.diag( + np.all( + ( + shp.vertex_zones["constraint"] + @ np.transpose(shp.vertex_normals + shp.vertices) + ) + <= np.expand_dims(shp.vertex_zones["bounds"], axis=2), + axis=1, + ) + ).astype(int) + * 2 + - 1 + ) + edge_inside_mult = ( + np.diag( + np.all( + ( + shp.edge_zones["constraint"] + @ np.transpose( + shp.edge_normals + + 0.5 + * ( + shp.vertices[shp.edges[:, 0]] + + shp.vertices[shp.edges[:, 1]] + ) + ) + ) + <= np.expand_dims(shp.edge_zones["bounds"], axis=2), + axis=1, + ) + ).astype(int) + * 2 + - 1 + ) + + # Updating bounds with the position of the polyhedron + vert_bounds = shp.vertex_zones["bounds"] + ( + shp.vertex_zones["constraint"] @ translation_vector + ) + edge_bounds = shp.edge_zones["bounds"] + ( + shp.edge_zones["constraint"] @ translation_vector + ) + face_bounds = shp.face_zones["bounds"] + ( + shp.face_zones["constraint"] @ translation_vector + ) + + points_trans = np.transpose( + points + ) # Have to take the transpose so that 'constraint @ points_trans' returns the right shape and values + max_value = ( + 3 * np.max(LA.norm(points - (translation_vector + shp.vertices[0]), axis=1)) + ) # Placeholder value, it is large so that it is not chosen when taking the min of the distances + + # Calculating the distances + + # Solving for the distances between the points and any relevant vertices + vert_dist = LA.norm( + np.repeat(np.expand_dims(points, axis=1), n_verts, axis=1) + - np.expand_dims(shp.vertices + translation_vector, axis=0), + axis=2, + ) * np.expand_dims(vert_inside_mult, axis=0) # Distances between two points + + # Taking the minimum of the distances for each point + vert_dist_arg = np.expand_dims(np.argmin(abs(vert_dist), axis=1), axis=1) + min_dist_arr = np.take_along_axis(vert_dist, vert_dist_arg, axis=1) + + # Solving for the distances between the points and any relevant edges + edge_bool = np.all( + (shp.edge_zones["constraint"] @ points_trans) + <= (np.expand_dims(edge_bounds, axis=2) + atol), + axis=1, + ) # <--- shape = (n_edges, n_points) + # edge_bool = edge_bool + np.allclose((shp.edge_zones["constraint"] @ points_trans), np.expand_dims(edge_bounds, axis=2), atol=1e-6) + if np.any(edge_bool): + # v--- shape = (number of True in edge_bool,) ---v + edge_used = np.transpose(np.tile(np.arange(0, n_edges, 1), (n_points, 1)))[ + edge_bool + ] # Contains the indices of the edges that hold True for edge_bool + e_points_used = np.tile(np.arange(0, n_points, 1), (n_edges, 1))[ + edge_bool + ] # Contains the indices of the points that hold True for edge_bool + + vert_on_edge = ( + shp.vertices[shp.edges[edge_used][:, 0]] + translation_vector + ) # Vertices that lie on the needed edges + + # Calculating the distances + edge_dist = np.ones((n_edges, n_points)) * max_value + edge_dist[edge_bool] = ( + point_to_edge_distance( + points[e_points_used], vert_on_edge, shp.edge_vectors[edge_used] + ) + * edge_inside_mult[edge_used] + ) # Distances between a point and a line + edge_dist = np.transpose(edge_dist) # <--- shape = (n_points, n_edges) + + # Taking the minimum of the distances for each point + edge_dist_arg = np.expand_dims(np.argmin(abs(edge_dist), axis=1), axis=1) + edge_dist = np.take_along_axis(edge_dist, edge_dist_arg, axis=1) + + min_dist_arr = np.concatenate((min_dist_arr, edge_dist), axis=1) + + # Solving for the distances between the points and any relevant faces + face_bool = np.all( + (shp.face_zones["constraint"] @ points_trans) + <= (np.expand_dims(face_bounds, axis=2) + atol), + axis=1, + ) # <--- shape = (n_tri_faces, n_points) + # face_bool = face_bool + np.allclose((shp.face_zones["constraint"] @ points_trans), np.expand_dims(face_bounds, axis=2), atol=1e-6) + if np.any(face_bool): + # v--- shape = (number of True in face_bool,) ---v + face_used = np.transpose(np.tile(np.arange(0, n_tri_faces, 1), (n_points, 1)))[ + face_bool + ] # Contains the indices of the triangulated faces that hold True for face_bool + f_points_used = np.tile(np.arange(0, n_points, 1), (n_tri_faces, 1))[ + face_bool + ] # Contains the indices of the points that hold True for face_bool + + vert_on_face = ( + (shp.face_zones["face_points"][face_used]) + translation_vector + ) # Vertices that lie on the needed faces + + # Calculating the distances + face_dist = np.ones((n_tri_faces, n_points)) * max_value + face_dist[face_bool] = point_to_face_distance( + points[f_points_used], vert_on_face, shp.face_zones["normals"][face_used] + ) # Distances between a point and a plane + face_dist = np.transpose(face_dist) # <--- shape = (n_points, n_tri_faces) + + # Taking the minimum of the distances for each point + face_dist_arg = np.expand_dims(np.argmin(abs(face_dist), axis=1), axis=1) + face_dist = np.take_along_axis(face_dist, face_dist_arg, axis=1) + + min_dist_arr = np.concatenate((min_dist_arr, face_dist), axis=1) + + min_dist_arg = np.expand_dims( + np.argmin(abs(min_dist_arr), axis=1), axis=1 + ) # determining the distances that are the shortest + true_min_dist = np.take_along_axis(min_dist_arr, min_dist_arg, axis=1).flatten() + + return true_min_dist + + +def shortest_displacement_to_surface( + shp, points: np.ndarray, translation_vector: np.ndarray +) -> np.ndarray: + """ + Solves for the shortest displacement between points and the surface of a polyhedron. + + This function calculates the shortest displacement by partitioning the space around + a polyhedron into zones: vertex, edge, and face. Determining the zone(s) a + point lies in, determines the displacement calculation(s) done. For a vertex zone, + the displacement is calculated between a point and the vertex. For an edge zone, the + displacement is calculated between a point and the edge. For a face zone, the + displacement is calculated between a point and the face. Zones are allowed to overlap, + and points can be in more than one zone. By taking the minimum of all the distances of + the calculated displacements, the shortest displacements are found. + + Args: + points (list or np.ndarray): positions of the points [shape = (n_points, 3)] + translation_vector (list or np.ndarray): translation vector of the polyhedron [shape = (3,)] + + Returns + ------- + np.ndarray: shortest displacements [shape = (n_points, 3)] + """ + points = np.asarray(points) + translation_vector = np.asarray(translation_vector) + + if translation_vector.shape[0] != 3 or len(translation_vector.shape) > 1: + raise ValueError( + f"Expected the shape of the polygon's position to be (3,), instead it got {translation_vector.shape}" + ) + + if points.shape == (3,): + points = points.reshape(1, 3) + + atol = 1e-8 + n_points = len(points) # number of inputted points + n_verts = len(shp.vertices) # number of vertices = number of vertex zones + n_edges = len(shp.edges) # number of edges = number of edge zones + n_tri_faces = len( + shp.face_zones["bounds"] + ) # number of triangulated faces = number of triangulated face zones + + # Updating bounds with the position of the polyhedron + vert_bounds = shp.vertex_zones["bounds"] + ( + shp.vertex_zones["constraint"] @ translation_vector + ) + edge_bounds = shp.edge_zones["bounds"] + ( + shp.edge_zones["constraint"] @ translation_vector + ) + face_bounds = shp.face_zones["bounds"] + ( + shp.face_zones["constraint"] @ translation_vector + ) + + coord_trans = np.transpose( + points + ) # Have to take the transpose so that 'constraint @ coord_trans' returns the right shape and values + max_value = ( + 3 * np.max(LA.norm(points - (translation_vector + shp.vertices[0]), axis=1)) + ) # Placeholder value, it is large so that it is not chosen when taking the min of the distances + + # Calculating the displacements + + # Solving for the displacements between the points and any relevant vertices + vert_disp = ( + -1 * np.repeat(np.expand_dims(points, axis=1), n_verts, axis=1) + ) + np.expand_dims( + shp.vertices + translation_vector, axis=0 + ) # Displacements between two point + + # Taking the minimum of the displacements for each point + vert_disp_min = np.expand_dims( + np.argmin(LA.norm(vert_disp, axis=2), axis=1), axis=(1, 2) + ) + min_disp_arr = np.take_along_axis(vert_disp, vert_disp_min, axis=1) + + # Solving for the displacements between the points and any relevant edges + edge_bool = np.all( + (shp.edge_zones["constraint"] @ coord_trans) + <= (np.expand_dims(edge_bounds, axis=2) + atol), + axis=1, + ) # <--- shape = (n_edges, n_points) + if np.any(edge_bool): + # v--- shape = (number of True in edge_bool,) ---v + edge_used = np.transpose(np.tile(np.arange(0, n_edges, 1), (n_points, 1)))[ + edge_bool + ] # Contains the indices of the edges that hold True for edge_bool + ecoords_used = np.tile(np.arange(0, n_points, 1), (n_edges, 1))[ + edge_bool + ] # Contains the indices of the points that hold True for edge_bool + + vert_on_edge = ( + shp.vertices[shp.edges[edge_used][:, 0]] + translation_vector + ) # Vertices that lie on the needed edges + + # Calculating the displacements + edge_disp = np.ones((n_edges, n_points, 3)) * max_value + edge_disp[edge_bool] = point_to_edge_displacement( + points[ecoords_used], vert_on_edge, shp.edge_vectors[edge_used] + ) # Displacements between a point and a line + edge_disp = np.transpose( + edge_disp, (1, 0, 2) + ) # <--- shape = (n_points, n_edges, 3) + + # Taking the minimum of the displacements for each point + edge_disp_arg = np.expand_dims( + np.argmin(LA.norm(edge_disp, axis=2), axis=1), axis=(1, 2) + ) + edge_disp = np.take_along_axis(edge_disp, edge_disp_arg, axis=1) + + min_disp_arr = np.concatenate((min_disp_arr, edge_disp), axis=1) + + # Solving for the displacements between the points and any relevant faces + face_bool = np.all( + (shp.face_zones["constraint"] @ coord_trans) + <= (np.expand_dims(face_bounds, axis=2) + atol), + axis=1, + ) # <--- shape = (n_tri_faces, n_points) + if np.any(face_bool): + # v--- shape = (number of True in face_bool,) ---v + face_used = np.transpose(np.tile(np.arange(0, n_tri_faces, 1), (n_points, 1)))[ + face_bool + ] # Contains the indices of the triangulated faces that hold True for face_bool + fcoords_used = np.tile(np.arange(0, n_points, 1), (n_tri_faces, 1))[ + face_bool + ] # Contains the indices of the points that hold True for face_bool + + vert_on_face = ( + (shp.face_zones["face_points"][face_used]) + translation_vector + ) # Vertices that lie on the needed faces + + # Calculating the displacements + face_disp = np.ones((n_tri_faces, n_points, 3)) * max_value + face_disp[face_bool] = point_to_face_displacement( + points[fcoords_used], vert_on_face, shp.face_zones["normals"][face_used] + ) # Displacements between a point and a plane + face_disp = np.transpose( + face_disp, (1, 0, 2) + ) # <--- shape = (n_points, n_tri_faces, 3) + + # Taking the minimum of the displacements for each point + face_disp_arg = np.expand_dims( + np.argmin(LA.norm(face_disp, axis=2), axis=1), axis=(1, 2) + ) + face_disp = np.take_along_axis(face_disp, face_disp_arg, axis=1) + + min_disp_arr = np.concatenate((min_disp_arr, face_disp), axis=1) + + disp_arr_bool = np.expand_dims( + np.argmin((LA.norm(min_disp_arr, axis=2)), axis=1), axis=(1, 2) + ) # determining the displacements that are shortest + true_min_disp = np.squeeze( + np.take_along_axis(min_disp_arr, disp_arr_bool, axis=1), axis=1 + ) + + return true_min_disp + + +def spheropolyhedron_shortest_displacement_to_surface( + shp, radius, points: np.ndarray, translation_vector: np.ndarray +) -> np.ndarray: + """ + Solves for the shortest displacement between points and the surface of a polyhedron. + + This function calculates the shortest displacement by partitioning the space around + a polyhedron into zones: vertex, edge, and face. Determining the zone(s) a + point lies in, determines the displacement calculation(s) done. For a vertex zone, + the displacement is calculated between a point and the vertex. For an edge zone, the + displacement is calculated between a point and the edge. For a face zone, the + displacement is calculated between a point and the face. Zones are allowed to overlap, + and points can be in more than one zone. By taking the minimum of all the distances of + the calculated displacements, the shortest displacements are found. + + Args: + points (list or np.ndarray): positions of the points [shape = (n_points, 3)] + translation_vector (list or np.ndarray): translation vector of the polyhedron [shape = (3,)] + + Returns + ------- + np.ndarray: shortest displacements [shape = (n_points, 3)] + """ + points = np.asarray(points) + translation_vector = np.asarray(translation_vector) + + if translation_vector.shape[0] != 3 or len(translation_vector.shape) > 1: + raise ValueError( + f"Expected the shape of the polygon's position to be (3,), instead it got {translation_vector.shape}" + ) + + if points.shape == (3,): + points = points.reshape(1, 3) + + atol = 1e-8 + n_points = len(points) # number of inputted points + n_verts = len(shp.vertices) # number of vertices = number of vertex zones + n_edges = len(shp.edges) # number of edges = number of edge zones + n_tri_faces = len( + shp.face_zones["bounds"] + ) # number of triangulated faces = number of triangulated face zones + + # Updating bounds with the position of the polyhedron + vert_bounds = shp.vertex_zones["bounds"] + ( + shp.vertex_zones["constraint"] @ translation_vector + ) + edge_bounds = shp.edge_zones["bounds"] + ( + shp.edge_zones["constraint"] @ translation_vector + ) + face_bounds = shp.face_zones["bounds"] + ( + shp.face_zones["constraint"] @ translation_vector + ) + + coord_trans = np.transpose( + points + ) # Have to take the transpose so that 'constraint @ coord_trans' returns the right shape and values + max_value = ( + 3 * np.max(LA.norm(points - (translation_vector + shp.vertices[0]), axis=1)) + ) # Placeholder value, it is large so that it is not chosen when taking the min of the distances + min_disp_arr = np.ones((n_points, 1, 3)) * max_value # Initial min_disp_arr + + # Calculating the displacements + + # Solving for the displacements between the points and any relevant vertices + vert_bool = np.all( + (shp.vertex_zones["constraint"] @ coord_trans) + <= np.expand_dims(vert_bounds, axis=2), + axis=1, + ) # <--- shape = (n_verts, n_points) + if np.any(vert_bool): + # v--- shape = (number of True in vert_bool,) ---v + vert_used = np.transpose(np.tile(np.arange(0, n_verts, 1), (n_points, 1)))[ + vert_bool + ] # Contains the indices of the vertices that hold True for vert_bool + vcoords_used = np.tile(np.arange(0, n_points, 1), (n_verts, 1))[ + vert_bool + ] # Contains the indices of the points that hold True for vert_bool + + # Calculating the displacements + vert_disp = np.ones((n_verts, n_points, 3)) * max_value + vert_disp[vert_bool] = (shp.vertices[vert_used] + translation_vector) - points[ + vcoords_used + ] # Displacements between two points + vert_disp = np.transpose( + vert_disp, (1, 0, 2) + ) # <--- shape = (n_points, n_verts, 3) + + # TODO: subtract radius*unit_displacement -- unless displacement is zero, then subtract radius*vert_normal + vert_zero_disp_bool = np.all(vert_disp == 0, axis=2) + vert_disp[vert_zero_disp_bool] = ( + vert_disp[vert_zero_disp_bool] + + radius + * ( + np.repeat(np.expand_dims(shp.vertex_normals, axis=0), n_points, axis=0)[ + vert_zero_disp_bool + ] + ) + ) + vert_disp[np.invert(vert_zero_disp_bool)] = vert_disp[ + np.invert(vert_zero_disp_bool) + ] - radius * ( + vert_disp[np.invert(vert_zero_disp_bool)] + / np.expand_dims( + np.linalg.norm(vert_disp[np.invert(vert_zero_disp_bool)], axis=1), + axis=1, + ) + ) + + # Taking the minimum of the displacements for each point + vert_disp_min = np.expand_dims( + np.argmin(LA.norm(vert_disp, axis=2), axis=1), axis=(1, 2) + ) + vert_disp = np.take_along_axis(vert_disp, vert_disp_min, axis=1) + + min_disp_arr = np.concatenate((min_disp_arr, vert_disp), axis=1) + + # Solving for the displacements between the points and any relevant edges + edge_bool = np.all( + (shp.edge_zones["constraint"] @ coord_trans) + <= (np.expand_dims(edge_bounds, axis=2) + atol), + axis=1, + ) # <--- shape = (n_edges, n_points) + if np.any(edge_bool): + # v--- shape = (number of True in edge_bool,) ---v + edge_used = np.transpose(np.tile(np.arange(0, n_edges, 1), (n_points, 1)))[ + edge_bool + ] # Contains the indices of the edges that hold True for edge_bool + ecoords_used = np.tile(np.arange(0, n_points, 1), (n_edges, 1))[ + edge_bool + ] # Contains the indices of the points that hold True for edge_bool + + vert_on_edge = ( + shp.vertices[shp.edges[edge_used][:, 0]] + translation_vector + ) # Vertices that lie on the needed edges + + # Calculating the displacements + edge_disp = np.ones((n_edges, n_points, 3)) * max_value + edge_disp[edge_bool] = point_to_edge_displacement( + points[ecoords_used], vert_on_edge, shp.edge_vectors[edge_used] + ) # Displacements between a point and a line + edge_disp = np.transpose( + edge_disp, (1, 0, 2) + ) # <--- shape = (n_points, n_edges, 3) + + # TODO: subtract radius*unit_displacement -- unless displacement is zero, then subtract radius*vert_normal + edge_zero_disp_bool = np.all(edge_disp == 0, axis=2) + edge_disp[edge_zero_disp_bool] = ( + edge_disp[edge_zero_disp_bool] + + radius + * ( + np.repeat(np.expand_dims(shp.edge_normals, axis=0), n_points, axis=0)[ + edge_zero_disp_bool + ] + ) + ) + edge_disp[np.invert(edge_zero_disp_bool)] = edge_disp[ + np.invert(edge_zero_disp_bool) + ] - radius * ( + edge_disp[np.invert(edge_zero_disp_bool)] + / np.expand_dims( + np.linalg.norm(edge_disp[np.invert(edge_zero_disp_bool)], axis=1), + axis=1, + ) + ) + + # Taking the minimum of the displacements for each point + edge_disp_arg = np.expand_dims( + np.argmin(LA.norm(edge_disp, axis=2), axis=1), axis=(1, 2) + ) + edge_disp = np.take_along_axis(edge_disp, edge_disp_arg, axis=1) + + min_disp_arr = np.concatenate((min_disp_arr, edge_disp), axis=1) + + # Solving for the displacements between the points and any relevant faces + face_bool = np.all( + (shp.face_zones["constraint"] @ coord_trans) + <= (np.expand_dims(face_bounds, axis=2) + atol), + axis=1, + ) # <--- shape = (n_tri_faces, n_points) + if np.any(face_bool): + # v--- shape = (number of True in face_bool,) ---v + face_used = np.transpose(np.tile(np.arange(0, n_tri_faces, 1), (n_points, 1)))[ + face_bool + ] # Contains the indices of the triangulated faces that hold True for face_bool + fcoords_used = np.tile(np.arange(0, n_points, 1), (n_tri_faces, 1))[ + face_bool + ] # Contains the indices of the points that hold True for face_bool + + vert_on_face = ( + (shp.face_zones["face_points"][face_used]) + translation_vector + ) # Vertices that lie on the needed faces + + # Calculating the displacements + face_disp = np.ones((n_tri_faces, n_points, 3)) * max_value + face_disp[face_bool] = point_to_face_displacement( + points[fcoords_used], vert_on_face, shp.face_zones["normals"][face_used] + ) # Displacements between a point and a plane + + # TODO: subtract radius*unit_displacement -- unless displacement is zero, then subtract radius*vert_normal + # TODO: if point is inside, add radius*unit_displacement instead + point_inside = (-1) * np.ones((n_tri_faces, n_points)) + point_inside[face_bool] = ( + point_to_face_distance( + points[fcoords_used], vert_on_face, shp.face_zones["normals"][face_used] + ) + < 0 + ).astype(int) * 2 - 1 # (+1) outside, (-1) inside + + face_zero_disp_bool = np.all(face_disp == 0, axis=2) + face_disp[face_zero_disp_bool] = ( + face_disp[face_zero_disp_bool] + + radius + * ( + np.repeat( + np.expand_dims( + ( + shp.face_zones["normals"] + / np.expand_dims( + np.linalg.norm(shp.face_zones["normals"], axis=1), + axis=1, + ) + ), + axis=1, + ), + n_points, + axis=1, + )[face_zero_disp_bool] + ) + ) + face_disp[np.invert(face_zero_disp_bool)] = face_disp[ + np.invert(face_zero_disp_bool) + ] + radius * np.expand_dims( + point_inside[np.invert(face_zero_disp_bool)], axis=1 + ) * ( + face_disp[np.invert(face_zero_disp_bool)] + / np.expand_dims( + np.linalg.norm(face_disp[np.invert(face_zero_disp_bool)], axis=1), + axis=1, + ) + ) + + face_disp = np.transpose( + face_disp, (1, 0, 2) + ) # <--- shape = (n_points, n_tri_faces, 3) + # Taking the minimum of the displacements for each point + face_disp_arg = np.expand_dims( + np.argmin(LA.norm(face_disp, axis=2), axis=1), axis=(1, 2) + ) + face_disp = np.take_along_axis(face_disp, face_disp_arg, axis=1) + + min_disp_arr = np.concatenate((min_disp_arr, face_disp), axis=1) + + disp_arr_bool = np.expand_dims( + np.argmin((LA.norm(min_disp_arr, axis=2)), axis=1), axis=(1, 2) + ) # determining the displacements that are shortest + true_min_disp = np.squeeze( + np.take_along_axis(min_disp_arr, disp_arr_bool, axis=1), axis=1 + ) + + return true_min_disp diff --git a/coxeter/shapes/convex_polyhedron.py b/coxeter/shapes/convex_polyhedron.py index 12afc3b7..8b5b57fd 100644 --- a/coxeter/shapes/convex_polyhedron.py +++ b/coxeter/shapes/convex_polyhedron.py @@ -118,6 +118,14 @@ def __init__(self, vertices): self._sort_simplices() self.sort_faces() + # For shortest distance functions + self._edge_face_neighbors = None + self._vertex_zones = None + self._edge_zones = None + self._face_zones = None + self._vertex_normals = None + self._edge_normals = None + def _consume_hull(self, hull): """Extract data from ConvexHull. diff --git a/coxeter/shapes/convex_spheropolygon.py b/coxeter/shapes/convex_spheropolygon.py index 098d8061..fce0112e 100644 --- a/coxeter/shapes/convex_spheropolygon.py +++ b/coxeter/shapes/convex_spheropolygon.py @@ -9,6 +9,9 @@ import numpy as np +from ._distance2d import ( + spheropolygon_shortest_displacement_to_surface, +) from .base_classes import Shape2D from .convex_polygon import ConvexPolygon, _is_convex from .polygon import _align_points_by_normal @@ -57,6 +60,10 @@ def __init__(self, vertices, radius, normal=None): if not _is_convex(self.vertices, self._polygon.normal): raise ValueError("The vertices do not define a convex polygon.") + self._vertex_zones = None + self._edge_zones = None + self._face_zones = None + @property def polygon(self): """:class:`~coxeter.shapes.ConvexPolygon`: The underlying polygon.""" @@ -107,6 +114,9 @@ def _rescale(self, scale): """ self.polygon._vertices *= scale self.radius *= scale + self._vertex_zones = None + self._edge_zones = None + self._face_zones = None @property def signed_area(self): @@ -302,3 +312,76 @@ def to_hoomd(self): self._polygon.centroid = old_centroid return hoomd_dict + + def shortest_distance_to_surface( + self, points, translation_vector=np.array([0, 0, 0]) + ): + """ + Solves for the shortest distance (magnitude) between points and + the surface of a polygon. + + This function calculates the shortest distance by partitioning + the space around a polygon into zones: vertex, edge, and face. + Determining the zone(s) a point lies in, determines the distance + calculation(s) done. For a vertex zone,the distance is calculated + between a point and the vertex. For an edge zone, the distance is + calculated between a point and the edge. For a face zone, the + distance is calculated between a point and the face. Zones are + allowed to overlap, and points can be in more than one zone. By + taking the minimum of all the calculated distances, the shortest + distances are found. + + Args: + points (list or :class:`numpy.ndarray`): + positions of the points [shape = (n_points,3) or (n_points,2)] + translation_vector (list or :class:`numpy.ndarray`): + translation vector of the polygon [shape = (3,) of (2,)] + (Default value: [0,0,0]) + + Returns + ------- + :class:`numpy.ndarray`: + the shortest distance of each point to the surface + [shape = (n_points,)] + """ + return np.linalg.norm( + spheropolygon_shortest_displacement_to_surface( + self._polygon, self.radius, points, translation_vector + ), + axis=1, + ) + + def shortest_displacement_to_surface( + self, points, translation_vector=np.array([0, 0, 0]) + ): + """ + Solves for the shortest displacement (vector) between points and + the surface of a polygon. + + This function calculates the shortest displacement by partitioning + the space around a polygon into zones: vertex, edge, and face. + Determining the zone(s) a point lies in, determines the displacement + calculation(s) done. For a vertex zone, the displacement is + calculated between a point and the vertex. For an edge zone, the + displacement is calculated between a point and the edge. For a face + zone, the displacement is calculated between a point and the face. + Zones are allowed to overlap, and points can be in more than one + zone. By taking the minimum of all the distances of the calculated + displacements, the shortest displacements are found. + + Args: + points (list or :class:`numpy.ndarray`): + positions of the points [shape = (n_points,3) or (n_points,2)] + translation_vector (list or :class:`numpy.ndarray`): + translation vector of the polygon [shape = (3,) or (2,)] + (Default value: [0,0,0]) + + Returns + ------- + :class:`numpy.ndarray`: + the shortest displacement of each point to the surface + [shape = (n_points, 3)] + """ + return spheropolygon_shortest_displacement_to_surface( + self._polygon, self.radius, points, translation_vector + ) diff --git a/coxeter/shapes/convex_spheropolyhedron.py b/coxeter/shapes/convex_spheropolyhedron.py index 6a4f2808..fa6337cd 100644 --- a/coxeter/shapes/convex_spheropolyhedron.py +++ b/coxeter/shapes/convex_spheropolyhedron.py @@ -9,6 +9,10 @@ import numpy as np +from ._distance3d import ( + shortest_distance_to_surface, + spheropolyhedron_shortest_displacement_to_surface, +) from .base_classes import Shape3D from .convex_polyhedron import ConvexPolyhedron from .utils import _hoomd_dict_mapping, _map_dict_keys @@ -69,6 +73,13 @@ def __init__(self, vertices, radius): self._polyhedron = ConvexPolyhedron(vertices) self.radius = radius + self._edge_face_neighbors = None + self._vertex_zones = None + self._edge_zones = None + self._face_zones = None + self._vertex_normals = None + self._edge_normals = None + @property def gsd_shape_spec(self): """dict: Get a :ref:`complete GSD specification `.""" # noqa: D401 @@ -97,6 +108,9 @@ def _rescale(self, scale): """ self.polyhedron._rescale(scale) self.radius *= scale + self._vertex_zones = None + self._edge_zones = None + self._face_zones = None @property def volume(self): @@ -364,3 +378,75 @@ def to_hoomd(self): self._polyhedron.centroid = old_centroid return hoomd_dict + + def shortest_distance_to_surface( + self, points, translation_vector=np.array([0, 0, 0]) + ): + """ + Solves for the shortest distance (magnitude) between points and + the surface of a spheropolyhedron. If the point lies inside the + spheropolyhedron, the distance is negative. + + This function calculates the shortest distance by partitioning + the space around a spheropolyhedron into zones: vertex, edge, and face. + Determining the zone(s) a point lies in, determines the distance + calculation(s) done. For a vertex zone,the distance is calculated + between a point and the vertex. For an edge zone, the distance is + calculated between a point and the edge. For a face zone, the + distance is calculated between a point and the face. Zones are + allowed to overlap, and points can be in more than one zone. By + taking the minimum of all the calculated distances, the shortest + distances are found. + + Args: + points (list or :class:`numpy.ndarray`): + positions of the points [shape = (n_points, 3)] + translation_vector (list or :class:`numpy.ndarray`): + translation vector of the spheropolyhedron [shape = (3,)] + (Default value: [0,0,0]) + + Returns + ------- + :class:`numpy.ndarray`: + the shortest distance of each point to the surface + [shape = (n_points,)] + """ + return ( + shortest_distance_to_surface(self._polyhedron, points, translation_vector) + - self.radius + ) + + def shortest_displacement_to_surface( + self, points, translation_vector=np.array([0, 0, 0]) + ): + """ + Solves for the shortest displacement (vector) between points and + the surface of a spheropolyhedron. + + This function calculates the shortest displacement by partitioning + the space around a spheropolyhedron into zones: vertex, edge, and face. + Determining the zone(s) a point lies in, determines the displacement + calculation(s) done. For a vertex zone, the displacement is + calculated between a point and the vertex. For an edge zone, the + displacement is calculated between a point and the edge. For a face + zone, the displacement is calculated between a point and the face. + Zones are allowed to overlap, and points can be in more than one + zone. By taking the minimum of all the distances of the calculated + displacements, the shortest displacements are found. + + Args: + points (list or :class:`numpy.ndarray`): + positions of the points [shape = (n_points, 3)] + translation_vector (list or :class:`numpy.ndarray`): + translation vector of the spheropolyhedron [shape = (3,)] + (Default value: [0,0,0]) + + Returns + ------- + :class:`numpy.ndarray`: + the shortest displacement of each point to the surface + [shape = (n_points, 3)] + """ + return spheropolyhedron_shortest_displacement_to_surface( + self._polyhedron, self.radius, points, translation_vector + ) diff --git a/coxeter/shapes/polygon.py b/coxeter/shapes/polygon.py index ea750091..9a1b834a 100644 --- a/coxeter/shapes/polygon.py +++ b/coxeter/shapes/polygon.py @@ -10,6 +10,13 @@ from ..extern.bentley_ottmann import poly_point_isect from ..extern.polytri import polytri +from ._distance2d import ( + get_edge_zones, + get_face_zones, + get_vert_zones, + shortest_displacement_to_surface, + shortest_distance_to_surface, +) from .base_classes import Shape2D from .circle import Circle from .utils import ( @@ -132,6 +139,10 @@ class Polygon(Shape2D): def __init__(self, vertices, normal=None, planar_tolerance=1e-5, test_simple=True): vertices = np.array(vertices, dtype=np.float64) + self._vertex_zones = None + self._edge_zones = None + self._face_zones = None + if len(vertices.shape) != 2 or vertices.shape[1] not in (2, 3): raise ValueError("Vertices must be specified as an Nx2 or Nx3 array.") if len(vertices) < 3: @@ -215,6 +226,9 @@ def _rescale(self, scale): Scale factor. """ self._vertices *= scale + self._vertex_zones = None + self._edge_zones = None + self._face_zones = None @property def perimeter(self): @@ -418,6 +432,19 @@ def centroid(self): def centroid(self, value): self._vertices += np.asarray(value) - self.centroid + if self._vertex_zones is not None: + self._vertex_zones["bounds"] += self._vertex_zones["constraint"] @ ( + np.asarray(value) - self.centroid + ) + if self._edge_zones is not None: + self._edge_zones["bounds"] += self._edge_zones["constraint"] @ ( + np.asarray(value) - self.centroid + ) + if self._face_zones is not None: + self._face_zones["bounds"] += self._face_zones["constraint"] @ ( + np.asarray(value) - self.centroid + ) + @property def edges(self): """:class:`numpy.ndarray`: Get the polygon's edges. @@ -838,3 +865,102 @@ def to_hoomd(self): self.centroid = old_centroid return hoomd_dict + + @property + def vertex_zones(self): + """dict: Get the constraints and bounds needed to partition the + volume surrounding a polygon into zones where the shortest + distance from any point that is within a vertex zone is the + distance between the point and the corresponding vertex. + """ + if self._vertex_zones is None: + self._vertex_zones = get_vert_zones(self) + return self._vertex_zones + + @property + def edge_zones(self): + """dict: Get the constraints and bounds needed to partition + the volume surrounding a polygon into zones where the + shortest distance from any point that is within an edge zone + is the distance between the point and the corresponding edge. + """ + if self._edge_zones is None: + self._edge_zones = get_edge_zones(self) + return self._edge_zones + + @property + def face_zones(self): + """dict: Get the constraints and bounds needed to partition + the volume surrounding a polygon into zones where the shortest + distance from any point that is within the face zone + is the distance between the point and the face of the polygon. + """ + if self._face_zones is None: + self._face_zones = get_face_zones(self) + return self._face_zones + + def shortest_distance_to_surface( + self, points, translation_vector=np.array([0, 0, 0]) + ): + """ + Solves for the shortest distance (magnitude) between points and + the surface of a polygon. + + This function calculates the shortest distance by partitioning + the space around a polygon into zones: vertex, edge, and face. + Determining the zone(s) a point lies in, determines the distance + calculation(s) done. For a vertex zone,the distance is calculated + between a point and the vertex. For an edge zone, the distance is + calculated between a point and the edge. For a face zone, the + distance is calculated between a point and the face. Zones are + allowed to overlap, and points can be in more than one zone. By + taking the minimum of all the calculated distances, the shortest + distances are found. + + Args: + points (list or :class:`numpy.ndarray`): + positions of the points [shape = (n_points,3) or (n_points,2)] + translation_vector (list or :class:`numpy.ndarray`): + translation vector of the polygon [shape = (3,) of (2,)] + (Default value: [0,0,0]) + + Returns + ------- + :class:`numpy.ndarray`: + the shortest distance of each point to the surface + [shape = (n_points,)] + """ + return shortest_distance_to_surface(self, points, translation_vector) + + def shortest_displacement_to_surface( + self, points, translation_vector=np.array([0, 0, 0]) + ): + """ + Solves for the shortest displacement (vector) between points and + the surface of a polygon. + + This function calculates the shortest displacement by partitioning + the space around a polygon into zones: vertex, edge, and face. + Determining the zone(s) a point lies in, determines the displacement + calculation(s) done. For a vertex zone, the displacement is + calculated between a point and the vertex. For an edge zone, the + displacement is calculated between a point and the edge. For a face + zone, the displacement is calculated between a point and the face. + Zones are allowed to overlap, and points can be in more than one + zone. By taking the minimum of all the distances of the calculated + displacements, the shortest displacements are found. + + Args: + points (list or :class:`numpy.ndarray`): + positions of the points [shape = (n_points,3) or (n_points,2)] + translation_vector (list or :class:`numpy.ndarray`): + translation vector of the polygon [shape = (3,) or (2,)] + (Default value: [0,0,0]) + + Returns + ------- + :class:`numpy.ndarray`: + the shortest displacement of each point to the surface + [shape = (n_points, 3)] + """ + return shortest_displacement_to_surface(self, points, translation_vector) diff --git a/coxeter/shapes/polyhedron.py b/coxeter/shapes/polyhedron.py index 73bc865c..e09ec4f9 100644 --- a/coxeter/shapes/polyhedron.py +++ b/coxeter/shapes/polyhedron.py @@ -12,6 +12,18 @@ from .. import io from ..extern.polytri import polytri +from ._distance3d import ( + get_edge_face_neighbors, + get_edge_normals, + get_edge_zones, + get_face_zones, + get_vert_normals, + get_vert_zones, + get_weighted_edge_normals, + get_weighted_vert_normals, + shortest_displacement_to_surface, + shortest_distance_to_surface, +) from .base_classes import Shape3D from .convex_polygon import ConvexPolygon, _is_convex from .polygon import Polygon, _is_simple @@ -138,6 +150,12 @@ def __init__(self, vertices, faces, faces_are_convex=None): self._faces_are_convex = faces_are_convex self._find_equations() self._find_neighbors() + self._edge_face_neighbors = None + self._vertex_zones = None + self._edge_zones = None + self._face_zones = None + self._vertex_normals = None + self._edge_normals = None def _find_equations(self): """Find the plane equations of the polyhedron faces.""" @@ -204,6 +222,9 @@ def _rescale(self, scale): """ self._vertices *= scale self._equations[:, 3] *= scale + self._vertex_zones = None + self._edge_zones = None + self._face_zones = None def merge_faces(self, atol=1e-8, rtol=1e-5): """Merge coplanar faces to a given tolerance. @@ -591,6 +612,19 @@ def centroid(self, value): self._vertices += np.asarray(value) - self.centroid self._find_equations() + if self._vertex_zones is not None: + self._vertex_zones["bounds"] += self._vertex_zones["constraint"] @ ( + np.asarray(value) - self.centroid + ) + if self._edge_zones is not None: + self._edge_zones["bounds"] += self._edge_zones["constraint"] @ ( + np.asarray(value) - self.centroid + ) + if self._face_zones is not None: + self._face_zones["bounds"] += self._face_zones["constraint"] @ ( + np.asarray(value) - self.centroid + ) + @property def bounding_sphere(self): """:class:`~.Sphere`: Get the polyhedron's bounding sphere.""" @@ -1057,3 +1091,155 @@ def save(self, filetype, filename): "filetype must be one of the following: OBJ, OFF, " "STL, PLY, VTK, X3D, HTML" ) + + @property + def edge_face_neighbors(self): + """:class:`numpy.ndarray`: Get the indices of the faces that + are adjacent to each edge. + + For a given edge vector oriented pointing upwards and from + an outside perspective of the polyhedron, the index of the + face to the left of the edge is given by the first column, + and the index of the face to the right of the edge is + given by the second column. + """ + if self._edge_face_neighbors is None: + self._edge_face_neighbors = get_edge_face_neighbors(self) + return self._edge_face_neighbors + + @property + def vertex_zones(self): + """dict: Get the constraints and bounds needed to partition the + volume surrounding a polyhedron into zones where the shortest + distance from any point that is within a vertex zone is the + distance between the point and the corresponding vertex. + """ + if self._vertex_zones is None: + self._vertex_zones = get_vert_zones(self) + return self._vertex_zones + + @property + def edge_zones(self): + """dict: Get the constraints and bounds needed to partition + the volume surrounding a polyhedron into zones where the + shortest distance from any point that is within an edge zone + is the distance between the point and the corresponding edge. + """ + if self._edge_zones is None: + self._edge_zones = get_edge_zones(self) + return self._edge_zones + + @property + def face_zones(self): + """dict: Get the constraints and bounds needed to partition + the volume surrounding a polyhedron into zones where the shortest + distance from any point that is within a triangulated face zone + is the distance between the point and the corresponding + triangulated face. + """ + if self._face_zones is None: + self._face_zones = get_face_zones(self) + return self._face_zones + + @property + def vertex_normals(self): + """:class:`numpy.ndarray`: Get the unit vector normals of vertices. + + The normals point outwards from the polyhedron. + """ + if self._vertex_normals is None: + self._vertex_normals = get_vert_normals(self) + return self._vertex_normals + + @property + def edge_normals(self): + """:class:`numpy.ndarray`: Get the unit vector normals of edges. + + The normals point outwards from the polyhedron. + """ + if self._edge_normals is None: + self._edge_normals = get_edge_normals(self) + return self._edge_normals + + @property + def weighted_vertex_normals(self): + """:class:`numpy.ndarray`: Get the weighted normals of vertices. + + The normals point outwards from the polyhedron. + """ + return get_weighted_vert_normals(self) + + @property + def weighted_edge_normals(self): + """:class:`numpy.ndarray`: Get the weighted normals of edges. + + The normals point outwards from the polyhedron. + """ + return get_weighted_edge_normals(self) + + def shortest_distance_to_surface( + self, points, translation_vector=np.array([0, 0, 0]) + ): + """ + Solves for the shortest distance (magnitude) between points and + the surface of a polyhedron. If the point lies inside the + polyhedron, the distance is negative. + + This function calculates the shortest distance by partitioning + the space around a polyhedron into zones: vertex, edge, and face. + Determining the zone(s) a point lies in, determines the distance + calculation(s) done. For a vertex zone,the distance is calculated + between a point and the vertex. For an edge zone, the distance is + calculated between a point and the edge. For a face zone, the + distance is calculated between a point and the face. Zones are + allowed to overlap, and points can be in more than one zone. By + taking the minimum of all the calculated distances, the shortest + distances are found. + + Args: + points (list or :class:`numpy.ndarray`): + positions of the points [shape = (n_points, 3)] + translation_vector (list or :class:`numpy.ndarray`): + translation vector of the polyhedron [shape = (3,)] + (Default value: [0,0,0]) + + Returns + ------- + :class:`numpy.ndarray`: + the shortest distance of each point to the surface + [shape = (n_points,)] + """ + return shortest_distance_to_surface(self, points, translation_vector) + + def shortest_displacement_to_surface( + self, points, translation_vector=np.array([0, 0, 0]) + ): + """ + Solves for the shortest displacement (vector) between points and + the surface of a polyhedron. + + This function calculates the shortest displacement by partitioning + the space around a polyhedron into zones: vertex, edge, and face. + Determining the zone(s) a point lies in, determines the displacement + calculation(s) done. For a vertex zone, the displacement is + calculated between a point and the vertex. For an edge zone, the + displacement is calculated between a point and the edge. For a face + zone, the displacement is calculated between a point and the face. + Zones are allowed to overlap, and points can be in more than one + zone. By taking the minimum of all the distances of the calculated + displacements, the shortest displacements are found. + + Args: + points (list or :class:`numpy.ndarray`): + positions of the points [shape = (n_points, 3)] + translation_vector (list or :class:`numpy.ndarray`): + translation vector of the polyhedron [shape = (3,)] + (Default value: [0,0,0]) + + Returns + ------- + :class:`numpy.ndarray`: + the shortest displacement of each point to the surface + [shape = (n_points, 3)] + """ + return shortest_displacement_to_surface(self, points, translation_vector) diff --git a/tests/test_polygon.py b/tests/test_polygon.py index 40205d8c..b1b61379 100644 --- a/tests/test_polygon.py +++ b/tests/test_polygon.py @@ -641,3 +641,183 @@ def test_edge_lengths(poly): np.diff(poly.vertices[expected_edges], axis=1).squeeze(), axis=1 ) np.testing.assert_allclose(poly.edge_lengths, lengths) + + +def test_shortest_distance_convex(): + tri_verts = np.array( + [[0, 0.5], [-0.25 * np.sqrt(3), -0.25], [0.25 * np.sqrt(3), -0.25]] + ) + triangle = ConvexPolygon(vertices=tri_verts) + + test_points = np.array( + [[3.5, 3.25, 0], [3, 3.75, 0], [3, 3.25, 0], [3, 3, 1], [3.25, 3.5, -1]] + ) + + distances = triangle.shortest_distance_to_surface( + test_points, translation_vector=np.array([3, 3, 0]) + ) + displacements = triangle.shortest_displacement_to_surface( + test_points, translation_vector=np.array([3, 3, 0]) + ) + + true_distances = np.array([0.3080127018, 0.25, 0, 1, 1.0231690965]) + true_displacements = np.array( + [ + [-0.2667468246, -0.1540063509, 0], + [0, -0.25, 0], + [0, 0, 0], + [0, 0, -1], + [-0.1875, -0.1082531755, 1], + ] + ) + + np.testing.assert_allclose(distances, true_distances) + np.testing.assert_allclose(displacements, true_displacements) + + +def test_shortest_distance_concave(): + verts = np.array( + [ + [0, 0.5], + [-0.125, 0.75], + [-0.25 * np.sqrt(3), -0.25], + [0.25 * np.sqrt(3), -0.25], + [0.25 * np.sqrt(3), 0.75], + ] + ) + concave_poly = Polygon(vertices=verts) + + test_points = np.array( + [ + [3.5, 3.25, 0], + [3, 3.75, 0], + [3, 3.25, 0], + [3, 3, -1], + [3.25, 3.5, -1], + [3 + 0.25 * np.sqrt(3), 4, 0], + ] + ) + + distances = concave_poly.shortest_distance_to_surface( + test_points, translation_vector=np.array([3, 3, 0]) + ) + displacements = concave_poly.shortest_displacement_to_surface( + test_points, translation_vector=np.array([3, 3, 0]) + ) + + true_distances = np.array( + [abs(0.25 * np.sqrt(3) - 0.5), np.sqrt(0.0125), 0, 1, 1, 0.25] + ) + true_displacements = np.array( + [ + [0.25 * np.sqrt(3) - 0.5, 0, 0], + [-0.1, -0.05, 0], + [0, 0, 0], + [0, 0, 1], + [0, 0, 1], + [0, -0.25, 0], + ] + ) + + np.testing.assert_allclose(distances, true_distances) + np.testing.assert_allclose(displacements, true_displacements) + + +def test_shortest_distance_general(): + # Creating a random polygon (nonconvex usually) + # np.random.seed(3) + random_angles = np.random.rand(15) * 2 * np.pi # angles + sorted_angles = np.sort(random_angles) + random_dist = np.random.rand(15) * 10 # from origin + + vertices = np.zeros((15, 2)) + vertices[:, 0] = random_dist * np.cos(sorted_angles) # x + vertices[:, 1] = random_dist * np.sin(sorted_angles) # y + + poly = Polygon(vertices=vertices, normal=[0, 0, 1]) + + points2d = np.random.rand(100, 2) * 20 - 10 + points3d = np.random.rand(150, 3) * 20 - 10 + + distances2d = poly.shortest_distance_to_surface(points2d) + distances3d = poly.shortest_distance_to_surface(points3d) + displacements2d = poly.shortest_displacement_to_surface(points2d) + displacements3d = poly.shortest_displacement_to_surface(points3d) + + np.testing.assert_allclose(distances2d, np.linalg.norm(displacements2d, axis=1)) + np.testing.assert_allclose(distances3d, np.linalg.norm(displacements3d, axis=1)) + + triangle_verts = [] + for tri in poly._triangulation(): + triangle_verts.append(list(tri)) + + triangle_verts = np.asarray(triangle_verts) + + tri_edges = ( + np.append( + triangle_verts[:, 1:], np.expand_dims(triangle_verts[:, 0], axis=1), axis=1 + ) + - triangle_verts + ) # edges point counterclockwise + + edges_90 = np.cross(tri_edges, poly.normal) # point outwards (n_triangles, 3, 3) + upper_bounds = np.sum(edges_90 * triangle_verts, axis=2) # (n_triangles, 3) + + def scipy_closest_point(point, edges_90, upper_bounds): + from scipy.optimize import LinearConstraint, minimize + + all_tri_distances = [] + all_tri_displacements = [] + tmps = [] + for triangle in zip(edges_90, upper_bounds, strict=True): + tri_min_point = minimize( + fun=lambda pt: np.linalg.norm(pt - point), # Function to optimize + x0=np.zeros(3), # Initial guess + constraints=[ + LinearConstraint( + np.append( + triangle[0].squeeze(), [[0, 0, 1], [0, 0, -1]], axis=0 + ), + -np.inf, + np.append(triangle[1].squeeze(), [0, 0]), + ) + ], + tol=1e-11, + ) + triangle_distance = np.linalg.norm(tri_min_point.x - point) + all_tri_distances.append(triangle_distance) + all_tri_displacements.append(tri_min_point.x - point) + + return np.min(all_tri_distances), all_tri_displacements[ + np.argmin(all_tri_distances) + ] + + scipy_distances2d = [] + scipy_displacements2d = [] + for point in points2d: + point = np.append(point, [0]) + + scipy_dist2d, scipy_displace2d = scipy_closest_point( + point, edges_90, upper_bounds + ) + scipy_distances2d.append(scipy_dist2d) + scipy_displacements2d.append(scipy_displace2d) + + scipy_distances3d = [] + scipy_displacements3d = [] + for point in points3d: + scipy_dist3d, scipy_displace3d = scipy_closest_point( + point, edges_90, upper_bounds + ) + scipy_distances3d.append(scipy_dist3d) + scipy_displacements3d.append(scipy_displace3d) + + scipy_distances2d = np.asarray(scipy_distances2d) + scipy_displacements2d = np.asarray(scipy_displacements2d) + scipy_distances3d = np.asarray(scipy_distances3d) + scipy_displacements3d = np.asarray(scipy_displacements3d) + + np.testing.assert_allclose(distances2d, scipy_distances2d, atol=2e-8) + np.testing.assert_allclose(displacements2d, scipy_displacements2d, atol=2e-5) + np.testing.assert_allclose(distances3d, scipy_distances3d, atol=2e-8) + np.testing.assert_allclose(displacements3d, scipy_displacements3d, atol=2e-5) diff --git a/tests/test_polyhedron.py b/tests/test_polyhedron.py index d3d1ab84..5188a427 100644 --- a/tests/test_polyhedron.py +++ b/tests/test_polyhedron.py @@ -954,3 +954,350 @@ def test_to_hoomd(poly): for i, face in enumerate(poly.faces): assert np.allclose(face, hoomd_dict["faces"][i]) + + +def test_shortest_distance_convex(): + cube_verts = np.array( + [ + [1, 1, 1], + [-1, 1, 1], + [1, -1, 1], + [1, 1, -1], + [-1, -1, 1], + [-1, 1, -1], + [1, -1, -1], + [-1, -1, -1], + ] + ) + cube = ConvexPolyhedron( + vertices=cube_verts + ) # , faces=[[1,5,7,4],[0,2,6,3],[2,4,7,6],[3,6,7,5],[0,1,4,2],[1,0,3,5]]) + + test_points = np.array( + [[3, 3, 3], [3, 3, 5], [5, 5, 1], [5, 4, 2], [3, 5, 5], [3, 4, 5], [3, 3, 6]] + ) + + distances = cube.shortest_distance_to_surface( + test_points, translation_vector=np.array([3, 3, 3]) + ) + displacements = cube.shortest_displacement_to_surface( + test_points, translation_vector=np.array([3, 3, 3]) + ) + np.testing.assert_allclose(np.abs(distances), np.linalg.norm(displacements, axis=1)) + + cube_surface_distance = cube.shortest_distance_to_surface( + test_points + displacements, translation_vector=np.array([3, 3, 3]) + ) + cube_surface_displacement = cube.shortest_displacement_to_surface( + test_points + displacements, translation_vector=np.array([3, 3, 3]) + ) + np.testing.assert_allclose(cube_surface_distance, np.zeros(len(test_points))) + np.testing.assert_allclose( + cube_surface_displacement, np.zeros((len(test_points), 3)) + ) + + true_distances = np.array([-1, 1, np.sqrt(3), 1, np.sqrt(2), 1, 2]) + true_displacements = np.array( + [[0, 0, -1], [-1, -1, 1], [-1, 0, 0], [0, -1, -1], [0, 0, -1], [0, 0, -2]] + ) + + np.testing.assert_allclose(distances, true_distances) + np.testing.assert_allclose(displacements[1:], true_displacements) + + +def test_shortest_distance_concave(): + test_points = np.array( + [ + [5, 5, 1], + [5, 4, 2], + [3, 5, 5], + [3, 4, 5], + [3, 3, 6], + [3.5, 2.5, 3], + [4.5, 3, 5], + [3, 3, 3], + [3, 3.5, 3.5], + [3, 3, 2.25], + ] + ) + + # --- PYRAMID Point on Cube Case --- + pyramidcube_verts = np.array( + [ + [1, 1, 1], + [-1, 1, 1], + [1, -1, 1], + [1, 1, -1], + [-1, -1, 1], + [-1, 1, -1], + [1, -1, -1], + [-1, -1, -1], + [0, 0, 3], + [0, 3, 0], + ] + ) + pyramid_faces = [ + [0, 1, 8], + [1, 4, 8], + [4, 2, 8], + [2, 0, 8], + [1, 0, 9], + [5, 1, 9], + [3, 5, 9], + [0, 3, 9], + [1, 5, 7, 4], + [0, 2, 6, 3], + [2, 4, 7, 6], + [3, 6, 7, 5], + ] + pyramidcube = Polyhedron(vertices=pyramidcube_verts, faces=pyramid_faces) + + pyramid_distances = pyramidcube.shortest_distance_to_surface( + test_points, translation_vector=np.array([3, 3, 3]) + ) + pyramid_displacements = pyramidcube.shortest_displacement_to_surface( + test_points, translation_vector=np.array([3, 3, 3]) + ) + np.testing.assert_allclose( + np.abs(pyramid_distances), np.linalg.norm(pyramid_displacements, axis=1) + ) + + pyramid_surface_distance = pyramidcube.shortest_distance_to_surface( + test_points + pyramid_displacements, translation_vector=np.array([3, 3, 3]) + ) + pyramid_surface_displacement = pyramidcube.shortest_displacement_to_surface( + test_points + pyramid_displacements, translation_vector=np.array([3, 3, 3]) + ) + np.testing.assert_allclose(pyramid_surface_distance, np.zeros(len(test_points))) + np.testing.assert_allclose( + pyramid_surface_displacement, np.zeros((len(test_points), 3)) + ) + + pyramid_true_distances = np.array( + [ + np.sqrt(3), + 1, + 3 / np.sqrt(5), + 1 / np.sqrt(5), + 0, + -0.5, + 2 / np.sqrt(5), + -1, + -1 * np.sqrt(0.5), + -0.25, + ] + ) + + np.testing.assert_allclose(pyramid_distances, pyramid_true_distances) + + # --- PRISM Point on Cube Case --- + prismcube_verts = np.array( + [ + [1, 1, 1], + [-1, 1, 1], + [1, -1, 1], + [1, 1, -1], + [-1, -1, 1], + [-1, 1, -1], + [1, -1, -1], + [-1, -1, -1], + [1, 0, 3], + [-1, 0, 3], + [1, 3, 0], + [-1, 3, 0], + ] + ) + prism_faces = [ + [0, 1, 9, 8], + [2, 8, 9, 4], + [2, 4, 7, 6], + [3, 6, 7, 5], + [3, 5, 11, 10], + [1, 0, 10, 11], + [0, 8, 2, 6, 3, 10], + [1, 11, 5, 7, 4, 9], + ] + prismcube = Polyhedron(vertices=prismcube_verts, faces=prism_faces) + + prism_distances = prismcube.shortest_distance_to_surface( + test_points, translation_vector=np.array([3, 3, 3]) + ) + prism_displacements = prismcube.shortest_displacement_to_surface( + test_points, translation_vector=np.array([3, 3, 3]) + ) + np.testing.assert_allclose( + np.abs(prism_distances), np.linalg.norm(prism_displacements, axis=1) + ) + + prism_surface_distance = prismcube.shortest_distance_to_surface( + test_points + prism_displacements, translation_vector=np.array([3, 3, 3]) + ) + prism_surface_displacement = prismcube.shortest_displacement_to_surface( + test_points + prism_displacements, translation_vector=np.array([3, 3, 3]) + ) + np.testing.assert_allclose(prism_surface_distance, np.zeros(len(test_points))) + np.testing.assert_allclose( + prism_surface_displacement, np.zeros((len(test_points), 3)) + ) + + prism_true_distances = np.array( + [ + np.sqrt(70) / 5, + 1, + 3 / np.sqrt(5), + 1 / np.sqrt(5), + 0, + -0.5, + 0.5, + -1, + -1 * np.sqrt(0.5), + -0.25, + ] + ) + + np.testing.assert_allclose(prism_distances, prism_true_distances) + + # --- INDENTED Cube Case --- + indented_cube_verts = np.array( + [ + [1, 1, 1], + [-1, 1, 1], + [1, -1, 1], + [1, 1, -1], + [-1, -1, 1], + [-1, 1, -1], + [1, -1, -1], + [-1, -1, -1], + [0, 0, -0.5], + ] + ) + indented_faces = [ + [0, 3, 5, 1], + [0, 2, 6, 3], + [1, 5, 7, 4], + [2, 4, 7, 6], + [3, 6, 7, 5], + [0, 1, 8], + [0, 8, 2], + [2, 8, 4], + [1, 4, 8], + ] + indented_cube = Polyhedron(vertices=indented_cube_verts, faces=indented_faces) + + indented_distances = indented_cube.shortest_distance_to_surface( + test_points, translation_vector=np.array([3, 3, 3]) + ) + indented_displacements = indented_cube.shortest_displacement_to_surface( + test_points, translation_vector=np.array([3, 3, 3]) + ) + np.testing.assert_allclose( + np.abs(indented_distances), np.linalg.norm(indented_displacements, axis=1) + ) + + indented_surface_distance = indented_cube.shortest_distance_to_surface( + test_points + indented_displacements, translation_vector=np.array([3, 3, 3]) + ) + indented_surface_displacement = indented_cube.shortest_displacement_to_surface( + test_points + indented_displacements, translation_vector=np.array([3, 3, 3]) + ) + np.testing.assert_allclose( + indented_surface_distance, np.zeros(len(test_points)), atol=2e-10 + ) + np.testing.assert_allclose( + indented_surface_displacement, np.zeros((len(test_points), 3)), atol=2e-10 + ) + + indented_true_distances = np.array( + [ + np.sqrt(3), + 1, + np.sqrt(2), + 1, + np.sqrt(5), + -0.1714985851, + np.sqrt(1.25), + 1 / np.sqrt(13), + 0.5 / np.sqrt(13), + -0.25, + ] + ) + + np.testing.assert_allclose(indented_distances, indented_true_distances) + + +def test_shortest_distance_convex_general(): + # Creating a random convex polyhedron + # np.random.seed(6) + random_theta = np.random.rand(20) * np.pi + random_phi = np.random.rand(20) * 2 * np.pi + radius = np.random.rand(1) * 5 + + vertices = np.zeros((20, 3)) + vertices[:, 0] = radius * np.sin(random_theta) * np.cos(random_phi) # x + vertices[:, 1] = radius * np.sin(random_theta) * np.sin(random_phi) # y + vertices[:, 2] = radius * np.cos(random_theta) + + poly = ConvexPolyhedron(vertices=vertices) + + points = np.random.rand(1500, 3) * 20 - 10 + + distances = poly.shortest_distance_to_surface(points) + displacements = poly.shortest_displacement_to_surface(points) + + np.testing.assert_allclose(np.abs(distances), np.linalg.norm(displacements, axis=1)) + + # Verifying that the displacements will move the points onto the surface + poly_surface_distance = poly.shortest_distance_to_surface(points + displacements) + poly_surface_displacement = poly.shortest_displacement_to_surface( + points + displacements + ) + + np.testing.assert_allclose(poly_surface_distance, np.zeros(len(points)), atol=2e-8) + np.testing.assert_allclose( + poly_surface_displacement, np.zeros((len(points), 3)), atol=2e-8 + ) + + def scipy_closest_point(point, surface_constraint, surface_bounds): + from scipy.optimize import LinearConstraint, minimize + + tri_min_point = minimize( + fun=lambda pt: np.linalg.norm(pt - point), # Function to optimize + x0=np.zeros(3), # Initial guess + constraints=[LinearConstraint(surface_constraint, -np.inf, surface_bounds)], + tol=1e-12, + ) + + distance = np.linalg.norm(tri_min_point.x - point) + displacement = tri_min_point.x - point + + return distance, displacement + + poly_constraint = poly.normals + poly_bounds = np.sum(poly_constraint * poly.face_centroids, axis=1) + + scipy_distances = [] + scipy_displacements = [] + for point in points: + outside_distance, outside_displacement = scipy_closest_point( + point, poly_constraint, poly_bounds + ) + + scipy_distances.append(outside_distance) + scipy_displacements.append(outside_displacement) + + scipy_distances = np.asarray(scipy_distances) + scipy_displacements = np.asarray(scipy_displacements) + + # Setting points inside the polyhedron to have 0 distance + is_zero_bool = distances <= 0 + + zero_inside_distances = distances + zero_inside_distances[is_zero_bool] = 0 + + zero_inside_displacements = displacements + zero_inside_displacements[is_zero_bool] = np.array([0, 0, 0]) + + np.testing.assert_allclose(zero_inside_distances, scipy_distances, atol=2e-8) + np.testing.assert_allclose( + zero_inside_displacements, scipy_displacements, atol=2e-5 + ) diff --git a/tests/test_spheropolygon.py b/tests/test_spheropolygon.py index 986f404c..7462be44 100644 --- a/tests/test_spheropolygon.py +++ b/tests/test_spheropolygon.py @@ -234,3 +234,147 @@ def test_to_hoomd(unit_rounded_square): hoomd_dict = shape.to_hoomd() for key, val in zip(dict_keys, dict_vals): assert np.allclose(hoomd_dict[key], val), f"{key}" + + +def test_shortest_distance_convex(): + tri_verts = np.array( + [[0, 0.5], [-0.25 * np.sqrt(3), -0.25], [0.25 * np.sqrt(3), -0.25]] + ) + triangle = ConvexSpheropolygon(vertices=tri_verts, radius=0.25) + + test_points = np.array( + [ + [3.5, 3.25, 0], + [3, 3.75, 0], + [3, 3.25, 0], + [3, 3, 1], + [3.25, 3.5, -1], + [3.5, 3.75, 1], + [3 - 0.25 * np.sqrt(3), 2.65, 0], + [3, 4, -1], + ] + ) + + distances = triangle.shortest_distance_to_surface( + test_points, translation_vector=np.array([3, 3, 0]) + ) + displacements = triangle.shortest_displacement_to_surface( + test_points, translation_vector=np.array([3, 3, 0]) + ) + + true_distances = np.array([0.0580127018, 0, 0, 1, 1, 1.0463612304, 0, 1.0307764064]) + true_displacements = np.array( + [ + [-0.0502404735, -0.0290063509, 0], + [0, 0, 0], + [0, 0, 0], + [0, 0, -1], + [0, 0, 1], + [-0.2667468244, -0.1540063509, -1], + [0, 0, 0], + [0, -0.25, 1], + ] + ) + + np.testing.assert_allclose(distances, true_distances) + np.testing.assert_allclose(displacements, true_displacements) + + +def test_shortest_distance_general(): + # Creating a random convex spheropolygon + # np.random.seed(3) + random_angles = np.random.rand(10) * 2 * np.pi # angles + sorted_angles = np.sort(random_angles) + random_dist = np.random.rand(1) * 10 # from origin + radius = np.random.rand(1) * 5 + + vertices = np.zeros((10, 2)) + vertices[:, 0] = random_dist * np.cos(sorted_angles) # x + vertices[:, 1] = random_dist * np.sin(sorted_angles) # y + + poly = ConvexSpheropolygon(vertices=vertices, radius=radius, normal=[0, 0, 1]) + + points2d = np.random.rand(100, 2) * 20 - 10 + points3d = np.random.rand(150, 3) * 20 - 10 + + distances2d = poly.shortest_distance_to_surface(points2d) + distances3d = poly.shortest_distance_to_surface(points3d) + displacements2d = poly.shortest_displacement_to_surface(points2d) + displacements3d = poly.shortest_displacement_to_surface(points3d) + + np.testing.assert_allclose(distances2d, np.linalg.norm(displacements2d, axis=1)) + np.testing.assert_allclose(distances3d, np.linalg.norm(displacements3d, axis=1)) + + edges_90 = np.cross( + poly._polygon.edge_vectors, poly.normal + ) # point outwards (10, 3) + upper_bounds = np.sum(edges_90 * poly.vertices, axis=1) # (10,) + + def scipy_closest_point(point, edges_90, upper_bounds): + from scipy.optimize import LinearConstraint, minimize + + tri_min_point = minimize( + fun=lambda pt: np.linalg.norm(pt - point), # Function to optimize + x0=np.zeros(3), # Initial guess + constraints=[ + LinearConstraint( + np.append(edges_90, [[0, 0, 1], [0, 0, -1]], axis=0), + -np.inf, + np.append(upper_bounds, [0, 0]), + ) + ], + tol=1e-10, + ) + + distance = np.linalg.norm(tri_min_point.x - point) + displacement = tri_min_point.x - point + + return distance, displacement + + # --- 2D --- + scipy_distances2d = [] + scipy_displacements2d = [] + for point in points2d: + point = np.append(point, [0]) + + scipy_dist2d, scipy_displace2d = scipy_closest_point( + point, edges_90, upper_bounds + ) + scipy_distances2d.append(scipy_dist2d) + scipy_displacements2d.append(scipy_displace2d) + + # --- 3D --- + scipy_distances3d = [] + scipy_displacements3d = [] + for point in points3d: + scipy_dist3d, scipy_displace3d = scipy_closest_point( + point, edges_90, upper_bounds + ) + + inplane_disp = scipy_displace3d - (scipy_displace3d @ poly.normal) * poly.normal + inplane_dist = np.linalg.norm(inplane_disp) + if inplane_dist < radius: + subtract_vector = inplane_disp + else: + subtract_vector = radius * inplane_disp / np.linalg.norm(inplane_disp) + + scipy_displace3d = scipy_displace3d - subtract_vector + scipy_displacements3d.append(scipy_displace3d) + scipy_distances3d.append(np.linalg.norm(scipy_displace3d)) + + scipy_distances2d = np.asarray(scipy_distances2d) - radius + is_zero2d = scipy_distances2d < 0 + scipy_distances2d[is_zero2d] = 0 + + scipy_displacements2d = np.asarray(scipy_displacements2d) + scipy_displacements2d = scipy_displacements2d - ( + radius + * scipy_displacements2d + / np.expand_dims(np.linalg.norm(scipy_displacements2d, axis=1), axis=1) + ) + scipy_displacements2d[is_zero2d] = np.array([0, 0, 0]) + + np.testing.assert_allclose(distances2d, scipy_distances2d, atol=2e-8) + np.testing.assert_allclose(displacements2d, scipy_displacements2d, atol=2e-5) + np.testing.assert_allclose(distances3d, scipy_distances3d, atol=2e-8) + np.testing.assert_allclose(displacements3d, scipy_displacements3d, atol=2e-5) diff --git a/tests/test_spheropolyhedron.py b/tests/test_spheropolyhedron.py index b51b1916..85d9fb14 100644 --- a/tests/test_spheropolyhedron.py +++ b/tests/test_spheropolyhedron.py @@ -158,3 +158,142 @@ def test_to_hoomd(poly, r): hoomd_dict = poly.to_hoomd() for key, val in zip(dict_keys, dict_vals): assert np.allclose(hoomd_dict[key], val), f"{key}" + + +def test_shortest_distance_convex(): + radius = 0.5 # np.random.rand(1) + verts = np.array( + [ + [1, 1, 1], + [-1, 1, 1], + [1, -1, 1], + [1, 1, -1], + [-1, -1, 1], + [-1, 1, -1], + [1, -1, -1], + [-1, -1, -1], + ] + ) + poly = ConvexSpheropolyhedron(vertices=verts, radius=radius) + + test_points = np.array( + [ + [3, 3, 3], + [3, 3, 5], + [5, 5, 1], + [5, 4, 2], + [3, 5, 5], + [3, 4, 5], + [3, 3, 6], + [4, 4, 4], + [4, 4, 3], + [4, 3, 3], + ] + ) + + distances = poly.shortest_distance_to_surface( + test_points, translation_vector=np.array([3, 3, 3]) + ) + displacements = poly.shortest_displacement_to_surface( + test_points, translation_vector=np.array([3, 3, 3]) + ) + + print(distances) + print(displacements) + + np.testing.assert_allclose(np.abs(distances), np.linalg.norm(displacements, axis=1)) + + poly_surface_distance = poly.shortest_distance_to_surface( + test_points + displacements, translation_vector=np.array([3, 3, 3]) + ) + poly_surface_displacement = poly.shortest_displacement_to_surface( + test_points + displacements, translation_vector=np.array([3, 3, 3]) + ) + np.testing.assert_allclose( + poly_surface_distance, np.zeros(len(test_points)), atol=1e-10 + ) + np.testing.assert_allclose( + poly_surface_displacement, np.zeros((len(test_points), 3)), atol=1e-10 + ) + + true_distances = ( + np.array([-1, 1, np.sqrt(3), 1, np.sqrt(2), 1, 2, 0, 0, 0]) - radius + ) + true_displacements = np.array( + [[0, 0, -1], [-1, -1, 1], [-1, 0, 0], [0, -1, -1], [0, 0, -1], [0, 0, -2]] + ) + true_displacements = true_displacements - radius * ( + true_displacements + / np.expand_dims(np.linalg.norm(true_displacements, axis=1), axis=1) + ) + + np.testing.assert_allclose(distances, true_distances) + np.testing.assert_allclose(displacements[1:7], true_displacements) + + +def test_shortest_distance_convex_general(): + # Creating a random convex spheropolyhedron + # np.random.seed(6) + random_theta = np.random.rand(20) * np.pi + random_phi = np.random.rand(20) * 2 * np.pi + radius = np.random.rand(1) * 5 # radius for the convex polyhedron + sph_poly_radius = np.random.rand(1) * 2 # radius for the spheropolyhedron + + vertices = np.zeros((20, 3)) + vertices[:, 0] = radius * np.sin(random_theta) * np.cos(random_phi) # x + vertices[:, 1] = radius * np.sin(random_theta) * np.sin(random_phi) # y + vertices[:, 2] = radius * np.cos(random_theta) + + poly = ConvexSpheropolyhedron(vertices=vertices, radius=sph_poly_radius) + + points = np.random.rand(1500, 3) * 20 - 10 + + distances = poly.shortest_distance_to_surface(points) + displacements = poly.shortest_displacement_to_surface(points) + + np.testing.assert_allclose(np.abs(distances), np.linalg.norm(displacements, axis=1)) + + # Verifying that the displacements will move the points onto the surface + poly_surface_distance = poly.shortest_distance_to_surface(points + displacements) + poly_surface_displacement = poly.shortest_displacement_to_surface( + points + displacements + ) + + np.testing.assert_allclose(poly_surface_distance, np.zeros(len(points)), atol=2e-8) + np.testing.assert_allclose( + poly_surface_displacement, np.zeros((len(points), 3)), atol=2e-8 + ) + + def scipy_closest_point(point, surface_constraint, surface_bounds): + from scipy.optimize import LinearConstraint, minimize + + tri_min_point = minimize( + fun=lambda pt: np.linalg.norm(pt - point), # Function to optimize + x0=np.zeros(3), # Initial guess + constraints=[LinearConstraint(surface_constraint, -np.inf, surface_bounds)], + tol=1e-12, + ) + + distance = np.linalg.norm(tri_min_point.x - point) + + return distance + + poly_constraint = poly._polyhedron.normals + poly_bounds = np.sum(poly_constraint * poly._polyhedron.face_centroids, axis=1) + + scipy_distances = [] + for point in points: + outside_distance = scipy_closest_point(point, poly_constraint, poly_bounds) + + scipy_distances.append(outside_distance) + + scipy_distances = np.asarray(scipy_distances) - sph_poly_radius + scipy_zero = scipy_distances < 0 + scipy_distances[scipy_zero] = 0 + + is_zero_bool = distances <= 0 + + zero_inside_distances = distances + zero_inside_distances[is_zero_bool] = 0 + + np.testing.assert_allclose(zero_inside_distances, scipy_distances, atol=2e-8)