diff -r 78d7e4e28e3e -r 00203fcc962f +parametrization/Curve.m
--- a/+parametrization/Curve.m	Mon Apr 08 20:15:37 2019 +0000
+++ b/+parametrization/Curve.m	Tue Apr 09 20:24:21 2019 +0000
@@ -103,7 +103,10 @@
             % Construct arcLength function using splines
             tvec = linspace(0,1,N);
             arcVec = obj.arcLength(0,tvec);
-            tFunc = spline(arcVec,tvec); % t as a function of arcLength
+
+            % t as a function of arcLength. Monotonicity-preserving cubic splines.
+            tFunc = @(arcLen) pchip(arcVec,tvec,arcLen);
+
             L = obj.arcLength(0,1);
             arcPar = @(s) tFunc(s*L);
 
@@ -349,8 +352,6 @@
     end
 end
 
-
-
 function g_norm = normalize(g0)
     g1 = g0(1,:);
     g2 = g0(2,:);

diff -r 78d7e4e28e3e -r 00203fcc962f +parametrization/dataSpline.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+parametrization/dataSpline.m	Tue Apr 09 20:24:21 2019 +0000
@@ -0,0 +1,19 @@
+% dataSpline calculates a Curve through the points f_i using cubic spline interpolation.
+% The spline curve is parametrized with the arc length parametrization
+% to facilitate better grids.
+%
+% f 	- m x D matrix of m points in D dimensions
+function C = dataSpline(f)
+	m = size(f, 1);
+
+	t = linspace(0,1,m);
+
+	pp_g = spapi(4, t, f');
+	pp_gp = fnder(pp_g);
+
+	g  = @(t) fnval(pp_g, t);
+	gp = @(t) fnval(pp_gp, t);
+
+	C = parametrization.Curve(g, gp);
+	C = C.arcLengthParametrization();
+end