From e8d0c90441041710c483d6e85c4e0d4839353f6e Mon Sep 17 00:00:00 2001 From: Simon Date: Wed, 16 Feb 2022 18:06:25 +0100 Subject: [PATCH] 160222 --- PySubdiv/PySubdiv.py | 58 +++++++++++++------ PySubdiv/backend/automatization.py | 16 +++-- PySubdiv/data/data_structure.py | 3 + .../optimization/variational_minimazation.py | 19 +++--- PySubdiv/subdivision_algorithms/loop.py | 1 - 5 files changed, 64 insertions(+), 33 deletions(-) diff --git a/PySubdiv/PySubdiv.py b/PySubdiv/PySubdiv.py index 7eac664..684ef2e 100644 --- a/PySubdiv/PySubdiv.py +++ b/PySubdiv/PySubdiv.py @@ -122,7 +122,7 @@ def faces(self, values): self.data['mesh_type'] = utils.mesh_type(values) self.data['faces'] = values - def add_vertex(self, vertex): + def add_vertex_old(self, vertex): """ Add a vertex to the end of the vertex array @@ -140,6 +140,44 @@ def add_vertex(self, vertex): else: raise ValueError("Vertices are not in the shape (n,3)") + def add_vertex(self, face_index, position=None): + """ + Add a vertex at the passed position or centroid of the passed faces (index). Deletes the passed face and create three new + faces to fill the hole in the mesh. + + Parameters + -------------- + face_index : int + index of face in self.faces + position: [x, y, z] float or None + when position is new vertex will be at that position otherwise at the centroid of the face + + """ + # check if passed face index is in the range of the meshes' faces + if face_index > len(self.faces) - 1: + print(f"face index {face_index} out of Range for mesh with {len(self.faces)} faces") + else: + if position is None: + position = self.face_centroids[face_index] # get centroid position of face + self.vertices = np.append(self.vertices, [position], axis=0) # add centroid to the vertices + new_faces = [] # list to store the three new faces + # create three new faces with the reference from the new vertex and the old vertices building the face + for counter, idx_vertex in enumerate(self.faces[face_index]): + new_face = [idx_vertex, len(self.vertices) - 1] + if counter < 2: + new_face.append(self.faces[face_index][counter + 1]) + else: + new_face.append(self.faces[face_index][0]) + new_faces.append(new_face) + # appending the new faces to self.faces + self.faces = np.append(self.faces, new_faces, axis=0) + # deleting the old face + self.faces = np.delete(self.faces, face_index, axis=0) + + self.define_face_centroids() + self.derive_edges() + self.edges_unique() + @property def ctr_points(self): """ @@ -837,7 +875,6 @@ def visualize_mesh(self): self.define_face_normals() self.define_face_centroids() self.define_vertex_normals() - model = visualize.print_model(self) return model @@ -862,23 +899,6 @@ def visualize_mesh_interactive(self, iteration=1, additional_meshes=None): subdivided_mesh = visualize.visualize_subdivision(self, iteration, additional_meshes) return subdivided_mesh - def pick_geodesic_path(self): - """ - Opens an interactive pyvista model of the mesh, where vertices can be picked and the shortest path are - calculated with Dijkstras's algorithm. - Returns a list of arrays with the indices of the vertices on the path - Parameters - ---------- - self: mesh - - Returns - --------- - nested list: - arrays of the vertices indices of the path - - """ - self.data['geodesic_path'] = visualize.geodesic_path(self) - return self.data['geodesic_path'] def save_mesh(self, filename): """ diff --git a/PySubdiv/backend/automatization.py b/PySubdiv/backend/automatization.py index 6ec758e..2092616 100644 --- a/PySubdiv/backend/automatization.py +++ b/PySubdiv/backend/automatization.py @@ -5,6 +5,7 @@ import math + def find_boundary_vertices(mesh, return_index=False, return_edge=False): """ Find the coordinates of the vertices that are on the boundary in a 2D plane @@ -29,11 +30,14 @@ def find_boundary_vertices(mesh, return_index=False, return_edge=False): List of vertex indices making up the boundary_edges """ - if "faces_edges_matrix" in mesh.data: + if "edge_faces_dictionary" in mesh.data: pass else: - mesh.faces_edges_incidence() - index_boundary_edges = np.nonzero(np.sum(mesh.data['faces_edges_matrix'], axis=0) == 1)[0] + mesh.edges_faces_connected() + index_boundary_edges = [] + for edge_idx in mesh.data["edge_faces_dictionary"]: + if len(mesh.data["edge_faces_dictionary"][edge_idx]) < 2: + index_boundary_edges.append(edge_idx) boundary_edges = mesh.data['unique_edges'][index_boundary_edges] index_boundary_vertices = np.unique(boundary_edges) boundary_vertices = mesh.data['unique_edges'][index_boundary_vertices] @@ -83,8 +87,10 @@ def find_corner_vertices(mesh, index_boundary_vertices=None, return_index=False) if index_boundary_vertices is None: index_boundary_vertices = np.unique(find_boundary_vertices(mesh, return_index=True)[0]) + if "vertices_connected_dictionary" not in mesh.data: + mesh.vertices_connected_dictionary() - connected_vertices_matrix = mesh.vertices_connected() # build vertex adjacency matrix + connected_vertices_matrix = mesh.data["vertices_connected_dictionary"] # build vertex adjacency matrix num_boundary_vertices = len(index_boundary_vertices) # number of vertices on the boundary of the mesh # create array to store connected vertices connected_vert_boundary_arr = np.ones((num_boundary_vertices, 4), dtype=int) @@ -316,3 +322,5 @@ def change_tolerance(tolerance): pl.add_mesh(model_2, style='wireframe', color='red') pl.show() return intersection_points[0] + + diff --git a/PySubdiv/data/data_structure.py b/PySubdiv/data/data_structure.py index 55c5739..4810df6 100644 --- a/PySubdiv/data/data_structure.py +++ b/PySubdiv/data/data_structure.py @@ -2,6 +2,7 @@ import numpy as np from PySubdiv.backend import calculation from scipy.sparse import coo_matrix +from scipy.sparse import csr_matrix import math @@ -164,9 +165,11 @@ def face_edges_incident(unique_edges, sorted_edges, faces, mesh_type): row = calculation.faces_to_edges(faces, mesh_type, return_index=True)[1] data = np.ones(len(col)) coo = coo_matrix((data, (row, col)), shape=(len(faces), len(unique_edges))).toarray() + return coo + def edge_faces_dict(unique_edges, sorted_edges): """ Returns incidence edges and faces. The dictionary key is the edge index and the values are the face indices. diff --git a/PySubdiv/optimization/variational_minimazation.py b/PySubdiv/optimization/variational_minimazation.py index 2eddde5..75f51ef 100644 --- a/PySubdiv/optimization/variational_minimazation.py +++ b/PySubdiv/optimization/variational_minimazation.py @@ -4,6 +4,7 @@ from PySubdiv.backend import optimization from PySubdiv.backend import automatization from scipy import optimize +from scipy import linalg import pyswarms as ps # Variational approach for fitting subdivision surfaces after Wu et. al. [http://dx.doi.org/10.1007/s41095-017-0088-2] @@ -212,20 +213,20 @@ def optimize(self, number_iteration=5, epsilon_0=1e-5, iterations_swarm=500, nr_ epsilon = np.inf iteration = 0 + A = self.control_cage.subdivide().data['subdivision_weight_matrix'].toarray() if len(self.variable_edges) == 0: print('It seems that all crease sharpness values of the edges are constrained, only the control cage can be' ' optimized') else: init_crease = self._control_cage.creases[self.variable_edges] while epsilon > epsilon_0 and iteration < number_iteration: - p_0 = self.p - result_p = optimize.minimize(objective_functions.objective_p, self.p, method='SLSQP', - args=(self._control_cage, self.v, self.z, self.lambda_z, self.a_z, - self.variable_vertices_idx, self.iterations_subdivision), - bounds=self.bounds_p) - - self.p = result_p.x.reshape(-1, 3) - self._control_cage.vertices[self.variable_vertices_idx] = self.p + p_0 = self.control_cage.vertices + # result_p = optimize.minimize(objective_functions.objective_p, self.p, method='SLSQP', + # args=(self._control_cage, self.v, self.z, self.lambda_z, self.a_z, + # self.variable_vertices_idx, self.iterations_subdivision), + # bounds=self.bounds_p) + self.p = linalg.inv(A.T @ A) @ (A.T @ self.v - A.T @ self.z - (A.T @ self.lambda_z)/self.a_z) + #self.p = result_p.x.reshape(-1, 3) if len(self.variable_edges) == 0: pass @@ -261,7 +262,7 @@ def optimize(self, number_iteration=5, epsilon_0=1e-5, iterations_swarm=500, nr_ self.z = result_z subdivided_mesh = self._control_cage.subdivide(self.iterations_subdivision) - + A = subdivided_mesh.data['subdivision_weight_matrix'].toarray() mp = subdivided_mesh.data['vertices'] self.lambda_z += self.a_z * (self.z - (self.v - mp)) diff --git a/PySubdiv/subdivision_algorithms/loop.py b/PySubdiv/subdivision_algorithms/loop.py index cdf55d7..18365f2 100644 --- a/PySubdiv/subdivision_algorithms/loop.py +++ b/PySubdiv/subdivision_algorithms/loop.py @@ -204,7 +204,6 @@ def generate_edge_points(self, subdivison_weight_matrix): -1, 2) # finding vertices opposite to the new edge point opposite_verts_zero_crease = np.zeros((len(connected_tris_zero_crease_masked), 2), dtype=int) - for i in range(len(connected_tris_zero_crease_masked)): connected_face_1 = set(self.data['faces'][connected_tris_zero_crease_masked[i][0]]) connected_face_2 = set(self.data['faces'][connected_tris_zero_crease_masked[i][1]])