diff --git a/src/cpp/binding_remesher_botsch.cpp b/src/cpp/binding_remesher_botsch.cpp
index ba540e30..32025930 100644
--- a/src/cpp/binding_remesher_botsch.cpp
+++ b/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);
         });
     
 }
\ No newline at end of file
diff --git a/src/cpp/remesher/collapse_edges.cpp b/src/cpp/remesher/collapse_edges.cpp
index 583b5651..800c9f39 100644
--- a/src/cpp/remesher/collapse_edges.cpp
+++ b/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;
 
 
 
diff --git a/src/cpp/remesher/collapse_edges.h b/src/cpp/remesher/collapse_edges.h
index 0873557e..27ed2659 100644
--- a/src/cpp/remesher/collapse_edges.h
+++ b/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
diff --git a/src/cpp/remesher/remesh_botsch.cpp b/src/cpp/remesher/remesh_botsch.cpp
index 8f6f8b44..2263fb7e 100644
--- a/src/cpp/remesher/remesh_botsch.cpp
+++ b/src/cpp/remesher/remesh_botsch.cpp
@@ -11,9 +11,10 @@
 #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;
@@ -21,23 +22,31 @@ void remesh_botsch(Eigen::MatrixXd & V,Eigen::MatrixXi & F, Eigen::VectorXd & ta
 
 	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){
diff --git a/src/cpp/remesher/remesh_botsch.h b/src/cpp/remesher/remesh_botsch.h
index 7ef4b261..4217989f 100644
--- a/src/cpp/remesher/remesh_botsch.h
+++ b/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);
diff --git a/src/cpp/remesher/split_edges.cpp b/src/cpp/remesher/split_edges.cpp
index e952b123..0f9b8bd5 100644
--- a/src/cpp/remesher/split_edges.cpp
+++ b/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 ***
diff --git a/src/cpp/remesher/split_edges.h b/src/cpp/remesher/split_edges.h
index 1dd99216..2833ff86 100644
--- a/src/cpp/remesher/split_edges.h
+++ b/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
diff --git a/src/cpp/remesher/split_edges_until_bound.cpp b/src/cpp/remesher/split_edges_until_bound.cpp
index 999b1d78..31910b72 100644
--- a/src/cpp/remesher/split_edges_until_bound.cpp
+++ b/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;
diff --git a/src/cpp/remesher/split_edges_until_bound.h b/src/cpp/remesher/split_edges_until_bound.h
index 13493b16..c3839f9f 100644
--- a/src/cpp/remesher/split_edges_until_bound.h
+++ b/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
diff --git a/src/cpp/remesher/tangential_relaxation.cpp b/src/cpp/remesher/tangential_relaxation.cpp
index ba111410..34330e35 100644
--- a/src/cpp/remesher/tangential_relaxation.cpp
+++ b/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;
diff --git a/src/cpp/remesher/tangential_relaxation.h b/src/cpp/remesher/tangential_relaxation.h
index 40bb8fb4..dc869514 100644
--- a/src/cpp/remesher/tangential_relaxation.h
+++ b/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
diff --git a/src/gpytoolbox/remesh_botsch.py b/src/gpytoolbox/remesh_botsch.py
index 7bf576cc..26ae8ab5 100644
--- a/src/gpytoolbox/remesh_botsch.py
+++ b/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,6 +24,8 @@ 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
     -------
@@ -31,7 +33,8 @@ def remesh_botsch(V, F, i=10, h=None, project=True, feature=np.array([], dtype=i
         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)
\ No newline at end of file