Mercurial > repos > public > sbplib

changeset 1098:36d092a00040 feature/dataspline

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Bugfix Curve.arcLengthParametrization() by using linear interpolation instead of higher-order spline for parameter as a funciton of arclength. Generalize dataSpline() to curves in higher dimensions.
author Martin Almquist <malmquist@stanford.edu>
date Tue, 09 Apr 2019 09:50:44 -0700
parents eec03e78c6f2
children 7f4aae76e06a
files +parametrization/Curve.m +parametrization/dataSpline.m
diffstat 2 files changed, 13 insertions(+), 19 deletions(-) [+]
line wrap: on
line diff
--- a/+parametrization/Curve.m	Mon Apr 08 22:30:47 2019 +0200
+++ b/+parametrization/Curve.m	Tue Apr 09 09:50:44 2019 -0700
@@ -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. Use linear interpolation to preserve monotonicity.
+            tFunc = @(arcLen)interp1(arcVec,tvec,arcLen, 'linear', 'extrap');
+
             L = obj.arcLength(0,1);
             arcPar = @(s) tFunc(s*L);
 
--- a/+parametrization/dataSpline.m	Mon Apr 08 22:30:47 2019 +0200
+++ b/+parametrization/dataSpline.m	Tue Apr 09 09:50:44 2019 -0700
@@ -1,28 +1,19 @@
-% Accepts data points (t_i, f_i) and returns a Curve object,
-% using spline interpolation.
+% 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.
 %
-% t_data 	- vector
-% f_data 	- vector
-% C 	- curve object
-function C = dataSpline(t_data, f_data)
+% f 	- m x D matrix of m points in D dimensions
+function C = dataSpline(f)
+	m = size(f, 1);
 
-	assert(length(t_data)==length(f_data),'Vectors must be same length');
-	m_data = length(t_data);
+	t = linspace(0,1,m);
 
-	pp_g = spapi(4, t_data, f_data);
+	pp_g = spapi(4, t, f');
 	pp_gp = fnder(pp_g);
 
-	% Reparametrize with parameter s from 0 to 1 to use Curve class
-    tmin = min(t_data);
-    tmax = max(t_data);
-    t = @(s) tmin + s*(tmax-tmin);
-    dt_ds = @(s) 0*s + (tmax-tmin);
-	g = @(s) [t(s); fnval(pp_g, t(s))];
-	gp = @(s) [dt_ds(s); fnval(pp_gp, t(s)).*dt_ds(s)];
+	g  = @(t) fnval(pp_g, t);
+	gp = @(t) fnval(pp_gp, t);
 
-	% Create Curve object and reparametrize with arclength parametrization
 	C = parametrization.Curve(g, gp);
-	C = C.arcLengthParametrization(m_data);
+	C = C.arcLengthParametrization();
 end