support vertex attributes in remeshing

src/cpp/binding_remesher_botsch.cpp

@@ -13,13 +13,14 @@ using EigenDRef = Ref<MatrixType, 0, EigenDStride>; //allows passing column/row
 
 void binding_remesh_botsch(py::module& m) {
     m.def("_remesh_botsch_cpp_impl",[](EigenDRef<MatrixXd> v,
-                         EigenDRef<MatrixXi> f, int i, double h, EigenDRef<VectorXi> feature, bool project)
+                         EigenDRef<MatrixXi> f, int i, double h, EigenDRef<VectorXi> feature, bool project, EigenDRef<MatrixXd> va)
         {
             Eigen::MatrixXd V(v);
             Eigen::MatrixXi F(f);
             Eigen::VectorXi FEATURE(feature);
-            remesh_botsch(V, F, h, i, FEATURE, project);
-            return std::make_tuple(V, F);
+            Eigen::MatrixXd VA(va);
+            remesh_botsch(V, F, h, i, FEATURE, project, VA);
+            return std::make_tuple(V, F, VA);
         });
 
 }

src/cpp/remesher/collapse_edges.cpp

@@ -28,7 +28,7 @@
 #include <igl/decimate_callback_types.h>
 using namespace std;
 
-void collapse_edges(Eigen::MatrixXd & V,Eigen::MatrixXi & F, Eigen::VectorXi & feature, Eigen::VectorXd & high, Eigen::VectorXd & low){
+void collapse_edges(Eigen::MatrixXd & V,Eigen::MatrixXi & F, Eigen::VectorXi & feature, Eigen::VectorXd & high, Eigen::VectorXd & low, Eigen::MatrixXd & VA){
         using namespace Eigen;
     MatrixXi E,uE,EI,EF;
     VectorXi EMAP,I,J;
@@ -189,13 +189,17 @@ std::vector<int> N;
 
     Eigen::VectorXd high_new,low_new;
     Eigen::VectorXi feature_new;
+    Eigen::MatrixXd VA_new;
+
     feature_new.resize(num_feature);
     high_new.resize(U.rows());
     low_new.resize(U.rows());
+    VA_new.resize(U.rows(),VA.cols());
     int j = 0;
     for (int s = 0; s<U.rows(); s++) {
         high_new(s) = high(I(s));
         low_new(s) = low(I(s));
+        VA_new.row(s) = VA.row(I(s));
         if (is_feature_vertex[I(s)]) {
             feature_new(j) = s;
             j = j+1;
@@ -209,7 +213,7 @@ std::vector<int> N;
     high = high_new;
     low = low_new;
     feature = feature_new;
-
+    VA = VA_new;
 
 
 

src/cpp/remesher/collapse_edges.h

@@ -5,7 +5,7 @@
 
 #include <Eigen/Core>
 
-void collapse_edges(Eigen::MatrixXd & V,Eigen::MatrixXi & F, Eigen::VectorXi & feature, Eigen::VectorXd & high, Eigen::VectorXd & low);
+void collapse_edges(Eigen::MatrixXd & V,Eigen::MatrixXi & F, Eigen::VectorXi & feature, Eigen::VectorXd & high, Eigen::VectorXd & low, Eigen::MatrixXd & VA);
 
 
 #endif

src/cpp/remesher/remesh_botsch.cpp

@@ -11,33 +11,42 @@
 #include <igl/avg_edge_length.h>
 #include <iostream>
 
-void remesh_botsch(Eigen::MatrixXd & V,Eigen::MatrixXi & F, Eigen::VectorXd & target,int iters, Eigen::VectorXi & feature, bool project){
+void remesh_botsch(Eigen::MatrixXd & V,Eigen::MatrixXi & F, Eigen::VectorXd & target,int iters, Eigen::VectorXi & feature, bool project, Eigen::MatrixXd & VA){
     Eigen::MatrixXd V0;
     Eigen::MatrixXi F0;
+    Eigen::MatrixXd VA0;
 
     Eigen::VectorXd high,low,lambda;
     high = 1.4*target;
     low = 0.7*target;
 
 	F0 = F;
 	V0 = V;
+	VA0 = VA;
     // Iterate the four steps
     for (int i = 0; i<iters; i++) {
-    	split_edges_until_bound(V,F,feature,high,low); // Split
-    	collapse_edges(V,F,feature,high,low); // Collapse
+    	split_edges_until_bound(V,F,feature,high,low,VA); // Split
+    	collapse_edges(V,F,feature,high,low,VA); // Collapse
     	equalize_valences(V,F,feature); // Flip
     	int n = V.rows();
     	lambda = Eigen::VectorXd::Constant(n,1.0);
 	if(!project){
 		V0 = V;
 		F0 = F;
+		VA0 = VA;
 	}
-	tangential_relaxation(V,F,feature,V0,F0,lambda); // Relax
+	tangential_relaxation(V,F,feature,V0,F0,lambda,VA,VA0); // Relax
     }
 }
 
+void remesh_botsch(Eigen::MatrixXd & V,Eigen::MatrixXi & F, Eigen::VectorXd & target,int iters, Eigen::VectorXi & feature, bool project){
+    Eigen::MatrixXd VA;
+	VA.resize(V.rows(), 1);
+	remesh_botsch(V,F,target,iters,feature,false,VA);
+}
+
 void remesh_botsch(Eigen::MatrixXd & V,Eigen::MatrixXi & F, Eigen::VectorXd & target,int iters, Eigen::VectorXi & feature){
-remesh_botsch(V,F,target,iters,feature,false);
+	remesh_botsch(V,F,target,iters,feature,false);
 }
 
 void remesh_botsch(Eigen::MatrixXd & V,Eigen::MatrixXi & F, Eigen::VectorXd & target,int iters){
@@ -76,10 +85,17 @@ void remesh_botsch(Eigen::MatrixXd & V,Eigen::MatrixXi & F, double target_double
 }
 
 void remesh_botsch(Eigen::MatrixXd & V,Eigen::MatrixXi & F, double target_double,int iters, Eigen::VectorXi feature, bool project){
+    Eigen::VectorXd target;
+    int n = V.rows();
+    target = Eigen::VectorXd::Constant(n,target_double);
+    remesh_botsch(V,F,target,iters,feature,project);
+}
+
+void remesh_botsch(Eigen::MatrixXd & V,Eigen::MatrixXi & F, double target_double,int iters, Eigen::VectorXi feature, bool project, Eigen::MatrixXd & VA){
 	Eigen::VectorXd target;
 	int n = V.rows();
 	target = Eigen::VectorXd::Constant(n,target_double);
-	remesh_botsch(V,F,target,iters,feature,project);
+	remesh_botsch(V,F,target,iters,feature,project,VA);
 }
 
 void remesh_botsch(Eigen::MatrixXd & V,Eigen::MatrixXi & F, double target_double){

src/cpp/remesher/remesh_botsch.h

@@ -5,8 +5,12 @@
 
 #include <Eigen/Core>
 
+void remesh_botsch(Eigen::MatrixXd & V,Eigen::MatrixXi & F,Eigen::VectorXd & target,int iters, Eigen::VectorXi & feature, bool project, Eigen::MatrixXd & VA);
+
 void remesh_botsch(Eigen::MatrixXd & V,Eigen::MatrixXi & F,Eigen::VectorXd & target,int iters, Eigen::VectorXi & feature, bool project);
 
+void remesh_botsch(Eigen::MatrixXd & V,Eigen::MatrixXi & F, double target_double,int iters, Eigen::VectorXi feature, bool project, Eigen::MatrixXd & VA);
+
 void remesh_botsch(Eigen::MatrixXd & V,Eigen::MatrixXi & F, double target_double,int iters, Eigen::VectorXi feature, bool project);
 
 void remesh_botsch(Eigen::MatrixXd & V,Eigen::MatrixXi & F,Eigen::VectorXd & target,int iters, Eigen::VectorXi & feature);

src/cpp/remesher/split_edges.cpp

@@ -25,7 +25,7 @@
 #include <igl/infinite_cost_stopping_condition.h>
 using namespace std;
 
-void split_edges(Eigen::MatrixXd & V, Eigen::MatrixXi & F, Eigen::MatrixXi & E0, Eigen::MatrixXi & uE, Eigen::VectorXi & EMAP0, std::vector<std::vector<int>> & uE2E,Eigen::VectorXd & high, Eigen::VectorXd & low,const std::vector<int> & edges_to_split){
+void split_edges(Eigen::MatrixXd & V, Eigen::MatrixXi & F, Eigen::MatrixXi & E0, Eigen::MatrixXi & uE, Eigen::VectorXi & EMAP0, std::vector<std::vector<int>> & uE2E,Eigen::VectorXd & high, Eigen::VectorXd & low,const std::vector<int> & edges_to_split, Eigen::MatrixXd & VA){
     using namespace Eigen;
 
     Eigen::VectorXi EMAP;
@@ -72,6 +72,7 @@ void split_edges(Eigen::MatrixXd & V, Eigen::MatrixXi & F, Eigen::MatrixXi & E0,
 
     F.conservativeResize(num_faces,3);
     V.conservativeResize(num_vertices,3);
+    VA.conservativeResize(num_vertices,VA.cols());
     high.conservativeResize(num_vertices);
     low.conservativeResize(num_vertices);
     uE.conservativeResize(num_uE,2);
@@ -176,6 +177,7 @@ void split_edges(Eigen::MatrixXd & V, Eigen::MatrixXi & F, Eigen::MatrixXi & E0,
 
         // naïve: use mid-point
         V.row(n+i) = (V.row(v1)+V.row(v2))/2;
+        VA.row(n+i) = (VA.row(v1)+VA.row(v2))/2;
         high(n+i) = (high(v1)+high(v2))/2;
         low(n+i) = (low(v1)+low(v2))/2;
         // *** UPDATE F ***

src/cpp/remesher/split_edges.h

@@ -5,7 +5,7 @@
 
 #include <Eigen/Core>
 
-void split_edges(Eigen::MatrixXd & V,Eigen::MatrixXi & F, Eigen::MatrixXi & E0, Eigen::MatrixXi & uE, Eigen::VectorXi & EMAP0, std::vector<std::vector<int>> & uE2E,Eigen::VectorXd & high, Eigen::VectorXd & low,const std::vector<int> & edges_to_split);
+void split_edges(Eigen::MatrixXd & V,Eigen::MatrixXi & F, Eigen::MatrixXi & E0, Eigen::MatrixXi & uE, Eigen::VectorXi & EMAP0, std::vector<std::vector<int>> & uE2E,Eigen::VectorXd & high, Eigen::VectorXd & low,const std::vector<int> & edges_to_split, Eigen::MatrixXd & VA);
 
 
 #endif

src/cpp/remesher/split_edges_until_bound.cpp

@@ -28,7 +28,7 @@
 #include "split_edges.h"
 using namespace std;
 
-void split_edges_until_bound(Eigen::MatrixXd & V,Eigen::MatrixXi & F, Eigen::VectorXi & feature, Eigen::VectorXd & high, Eigen::VectorXd & low){
+void split_edges_until_bound(Eigen::MatrixXd & V,Eigen::MatrixXi & F, Eigen::VectorXi & feature, Eigen::VectorXd & high, Eigen::VectorXd & low, Eigen::MatrixXd & VA){
 
     using namespace Eigen;
     int m = F.rows();
@@ -82,7 +82,7 @@ void split_edges_until_bound(Eigen::MatrixXd & V,Eigen::MatrixXi & F, Eigen::Vec
 
             //std::cout << "Before call to split_edges" << std::endl;
             //std::cout << edges_to_split.size() << std::endl;
-            split_edges(V,F,E,uE,EMAP,uE2E,high,low,edges_to_split);
+            split_edges(V,F,E,uE,EMAP,uE2E,high,low,edges_to_split,VA);
             //igl::writeOBJ("test.obj",V,F);
             //igl::unique_edge_map(F,E,uE,EMAP,uE2E);
             //std::cout << igl::is_edge_manifold(F) << std::endl;

src/cpp/remesher/split_edges_until_bound.h

@@ -5,7 +5,7 @@
 
 #include <Eigen/Core>
 
-void split_edges_until_bound(Eigen::MatrixXd & V,Eigen::MatrixXi & F, Eigen::VectorXi & feature, Eigen::VectorXd & high, Eigen::VectorXd & low);
+void split_edges_until_bound(Eigen::MatrixXd & V,Eigen::MatrixXi & F, Eigen::VectorXi & feature, Eigen::VectorXd & high, Eigen::VectorXd & low, Eigen::MatrixXd & VA);
 
 
 #endif

src/cpp/remesher/tangential_relaxation.cpp

@@ -5,6 +5,8 @@
 #include <igl/adjacency_list.h>
 #include <igl/per_face_normals.h>
 #include <igl/barycenter.h>
+#include <igl/barycentric_coordinates.h>
+#include <igl/barycentric_interpolation.h>
 #include <igl/pinv.h>
 #include <igl/writeOBJ.h>
 #include <igl/edges.h>
@@ -26,7 +28,7 @@
 using namespace std;
 
 void tangential_relaxation(Eigen::MatrixXd & V,Eigen::MatrixXi & F, Eigen::VectorXi & feature,
-        Eigen::MatrixXd & V0 ,Eigen::MatrixXi & F0, Eigen::VectorXd & lambda){
+        Eigen::MatrixXd & V0 ,Eigen::MatrixXi & F0, Eigen::VectorXd & lambda, Eigen::MatrixXd & VA, Eigen::MatrixXd & VA0){
     using namespace Eigen;
         MatrixXd Q,P,N,V_projected,V_fixed;
         VectorXd dblA,sqrD;
@@ -37,7 +39,9 @@ void tangential_relaxation(Eigen::MatrixXd & V,Eigen::MatrixXi & F, Eigen::Vecto
         Eigen::MatrixXd SV;
         Eigen::MatrixXi SVI,SVJ;
 
-
+        MatrixXd VA_projected;
+        VA_projected.resize(VA.rows(), VA.cols());
+
 
 
 
@@ -54,9 +58,9 @@ void tangential_relaxation(Eigen::MatrixXd & V,Eigen::MatrixXi & F, Eigen::Vecto
 
 
         //for (int j = 0; j < m; j++) {
-        //    vertex_areas[F(j,0)] = vertex_areas[F(j,0)] + (abs(dblA(j))/6);
-        //    vertex_areas[F(j,1)] = vertex_areas[F(j,1)] + (abs(dblA(j))/6);
-        //    vertex_areas[F(j,2)] = vertex_areas[F(j,2)] + (abs(dblA(j))/6);
+            //    vertex_areas[F(j,0)] = vertex_areas[F(j,0)] + (abs(dblA(j))/6);
+            //    vertex_areas[F(j,1)] = vertex_areas[F(j,1)] + (abs(dblA(j))/6);
+            //    vertex_areas[F(j,2)] = vertex_areas[F(j,2)] + (abs(dblA(j))/6);
         //}
 
 
@@ -116,8 +120,16 @@ void tangential_relaxation(Eigen::MatrixXd & V,Eigen::MatrixXi & F, Eigen::Vecto
 //
 	igl::point_mesh_squared_distance(V,V0,F0,sqrD,sqrI,V_projected);
 //
+    // Use the barycentric coordinates associated with the reprojections to interpolate the vertex attributes
+    MatrixXd B;
+    MatrixXd vA = V0(F0(sqrI, 0), all);
+    MatrixXd vB = V0(F0(sqrI, 1), all);
+    MatrixXd vC = V0(F0(sqrI, 2), all);
+    igl::barycentric_coordinates(V_projected, vA, vB, vC, B);
+    igl::barycentric_interpolation(VA0, F0, B, sqrI, VA_projected);
 //
     V = V_projected;
+    VA = VA_projected;
 //	igl::writeOBJ("post-project.obj",V,F);
 //    igl::remove_duplicate_vertices(V,0,SV,SVI,SVJ);
 //    std::cout << V.rows()-SV.rows() << std::endl;

src/cpp/remesher/tangential_relaxation.h

@@ -6,7 +6,7 @@
 #include <Eigen/Core>
 
 void tangential_relaxation(Eigen::MatrixXd & V,Eigen::MatrixXi & F, Eigen::VectorXi & feature,
-Eigen::MatrixXd & V0 ,Eigen::MatrixXi & F0, Eigen::VectorXd & lambda);
+Eigen::MatrixXd & V0 ,Eigen::MatrixXi & F0, Eigen::VectorXd & lambda, Eigen::MatrixXd & VA, Eigen::MatrixXd & VA0);
 
 
 #endif

src/gpytoolbox/remesh_botsch.py

@@ -5,7 +5,7 @@
 from gpytoolbox.non_manifold_edges import non_manifold_edges
 
 
-def remesh_botsch(V, F, i=10, h=None, project=True, feature=np.array([], dtype=int)):
+def remesh_botsch(V, F, i=10, h=None, project=True, feature=np.array([], dtype=int), vertex_attrs=None):
     """Remesh a triangular mesh to have a desired edge length
     
     Use the algorithm described by Botsch and Kobbelt's "A Remeshing Approach to Multiresolution Modeling" to remesh a triangular mesh by alternating iterations of subdivision, collapse, edge flips and collapses.
@@ -24,14 +24,17 @@ def remesh_botsch(V, F, i=10, h=None, project=True, feature=np.array([], dtype=i
         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).
+    vertex_attrs : numpy double array
+        Matrix of per-vertex attributes of arbitrary (flat) size
 
     Returns
     -------
     U : numpy double array
         Matrix of output triangle mesh vertex coordinates
     G : numpy int array
         Matrix of output triangle mesh indices into U
-
+    A : numpy double array
+        Matrix of upsampled per-vertex attributes
 
     Notes
     -----
@@ -63,6 +66,12 @@ def remesh_botsch(V, F, i=10, h=None, project=True, feature=np.array([], dtype=i
 
     feature = np.concatenate((feature, boundary_vertices(F)), dtype=np.int32)
 
+    if vertex_attrs is None:
+        return_va = False
+        vertex_attrs = np.zeros((V.size(0), 1), dtype=np.float64)
+    else:
+        return_va = True
+
     # reorder feature nodes to the beginning of the array (contributed by Michael Jäger)
     if feature.shape[0] > 0:
         # feature indices need to be unique (including the boundary_vertices)
@@ -82,6 +91,7 @@ def remesh_botsch(V, F, i=10, h=None, project=True, feature=np.array([], dtype=i
 
         # reorder vertex coordinates
         V = V[order]
+        vertex_attrs = vertex_attrs[order]
         # reorder faces
         F = tmp[F]
         # features are now 0 to n_features
@@ -97,6 +107,5 @@ def remesh_botsch(V, F, i=10, h=None, project=True, feature=np.array([], dtype=i
         # return error
         raise ValueError("Input mesh is non-manifold.")
 
-    v, f = _remesh_botsch_cpp_impl(V, F.astype(np.int32), i, h, feature, project)
-
-    return v, f
+    v, f, va = _remesh_botsch_cpp_impl(V, F.astype(np.int32), i, h, feature, project, vertex_attrs)
+    return (v,f,va) if return_va else (v, f)