160222

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):
         """

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]
+
+

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.

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))

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]])