keep indices of feature nodes, closes #30 (#45)

src/gpytoolbox/remesh_botsch.py

@@ -1,4 +1,5 @@
 import numpy as np
+import warnings
 from gpytoolbox.boundary_vertices import boundary_vertices
 from gpytoolbox.halfedge_lengths import halfedge_lengths
 
@@ -18,7 +19,7 @@ def remesh_botsch(V,F,i=10,h=None,project=True,feature = np.array([],dtype=int))
     h : double, optional (default 0.1)
         Desired edge length (if None, will pick average edge length)
     feature : numpy int array, optional (default np.array([],dtype=int))
-        List of indeces of feature vertices that should not change (i.e., they will also be in the output) 
+        List of indices of feature vertices that should not change (i.e., they will also be in the output). They will be placed at the beginning of the output array in the same order (as long as they were unique).
     project : bool, optional (default True)
         Whether to reproject the mesh to the input (otherwise, it will smooth over iterations).
 
@@ -29,6 +30,11 @@ def remesh_botsch(V,F,i=10,h=None,project=True,feature = np.array([],dtype=int))
     G : numpy int array
         Matrix of output triangle mesh indices into U
 
+
+    Notes
+    -----
+    The ordering of the output can be somewhat confusing. The output vertices are ordered as follows: [feature vertices, boundary vertices, other vertices]. If a vertex is both a feature and boundary one, it is treated as a feature vertex for the purposes of the ordering. For a more in-depth explanation see [PR #45](https://github.com/sgsellan/gpytoolbox/pull/45).
+
     Examples
     --------
     ```python
@@ -46,11 +52,36 @@ def remesh_botsch(V,F,i=10,h=None,project=True,feature = np.array([],dtype=int))
     if (h is None):
         h = np.mean(halfedge_lengths(V,F))
 
+    # check that feature is unique
+    if feature.shape[0] > 0:
+        if np.unique(feature).shape[0] != feature.shape[0]:
+            warnings.warn("Feature array is not unique. We will compute its unique entries and use those as an input. We recommend you do this yourself to avoid this warning.")
+
     feature = np.concatenate((feature,boundary_vertices(F)),dtype=np.int32)
-    # print(feature)
-    # print(boundary_vertices(F))
-    # bV = boundary_vertices(F)
 
+    # reorder feature nodes to the beginning of the array
+    if feature.shape[0] > 0:
+        # feature indices need to be unique (including the boundary_vertices)
+        tmp, ind = np.unique(feature, return_index=True)
+        # unique feature array while preserving the order [feature, boundary_vertices]
+        feature = tmp[np.argsort(ind)]
+
+        # number of vertices
+        n_vertices = V.shape[0]
+        # 0 ... n_vertices array
+        old_order = np.arange(n_vertices, dtype=np.int32)
+        # new order
+        order = np.concatenate((feature, np.delete(old_order, feature)), dtype=np.int32)
+        # generate tmp array for reordering mesh indices
+        tmp = np.empty(n_vertices, dtype=np.int32)
+        tmp[order] = old_order  # this line will fail if features are not unique
+
+        # reorder vertex coordinates
+        V = V[order]
+        # reorder faces
+        F = tmp[F]
+        # features are now 0 to n_features
+        feature = old_order[:feature.shape[0]]
 
     v,f = _remesh_botsch_cpp_impl(V,F.astype(np.int32),i,h,feature,project)
 

test/test_remesh_botsch.py

@@ -47,6 +47,48 @@ def test_with_boundary(self):
             dist = np.min(np.linalg.norm(np.tile(boundary_verts[i,:][None,:],(boundary_verts_output.shape[0],1)) - boundary_verts_output,axis=1))
             self.assertTrue(dist==0.0)
 
+    def test_with_unique_features_closed(self):
+        np.random.seed(0)
+        v,f = gpytoolbox.read_mesh("test/unit_tests_data/bunny_oded.obj")
+        # pick random faces of the model that are fixed
+        feature = f[np.random.choice(range(f.shape[0]), v.shape[0]//1000, replace=False)].flatten()
+        u,g = gpytoolbox.remesh_botsch(v,f.astype(np.int32),20,0.01,True,feature=feature)
+        self.assertTrue(np.allclose(v[feature], u[:feature.shape[0]]))
+
+    def test_with_unique_features(self):
+        np.random.seed(0)
+        v,f = gpytoolbox.read_mesh("test/unit_tests_data/bunny.obj")
+        # pick random faces of the model that are fixed
+        feature = f[np.random.choice(range(f.shape[0]), v.shape[0]//1000, replace=False)].flatten()
+        u,g = gpytoolbox.remesh_botsch(v,f.astype(np.int32),20,0.01,True,feature=feature)
+        self.assertTrue(np.allclose(v[feature], u[:feature.shape[0]]))
+
+    def test_with_not_unique_features(self):
+        np.random.seed(8)
+        v,f = gpytoolbox.read_mesh("test/unit_tests_data/bunny.obj")
+        # pick random faces of the model that are fixed
+        feature = f[np.random.choice(range(f.shape[0]), v.shape[0]//1000, replace=False)].flatten()
+        # check that they are not unique
+        self.assertFalse(feature.shape[0] == np.unique(feature).shape[0])
+        u,g = gpytoolbox.remesh_botsch(v,f.astype(np.int32),20,0.01,True,feature=feature)
+        # unique feature nodes
+        tmp, ind = np.unique(feature, return_index=True)
+        feature_unique = tmp[np.argsort(ind)]
+        self.assertTrue(np.allclose(v[feature_unique], u[:feature_unique.shape[0]]))
+
+    def test_with_not_unique_features_and_boundary(self):
+        np.random.seed(8)
+        v,f = gpytoolbox.read_mesh("test/unit_tests_data/bunny.obj")
+        # pick random faces of the model that are fixed and add some boundary nodes
+        feature = f[np.random.choice(range(f.shape[0]), v.shape[0]//1000, replace=False)].flatten()
+        feature = np.concatenate((feature, np.random.choice(gpytoolbox.boundary_vertices(f), 20, replace=False)))
+        # check that they are not unique
+        self.assertFalse(feature.shape[0] == np.unique(feature).shape[0])
+        u,g = gpytoolbox.remesh_botsch(v,f.astype(np.int32),20,0.01,True,feature=feature)
+        # unique feature nodes
+        tmp, ind = np.unique(feature, return_index=True)
+        feature_unique = tmp[np.argsort(ind)]
+        self.assertTrue(np.allclose(v[feature_unique], u[:feature_unique.shape[0]]))
 
     # def test_github_issue_30(self):
     #     np.random.seed(0)