Mercurial > repos > public > sbplib

comparison +grid/bspline.m @ 197:4f7930d2d2c4

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Added function to calculate bspline.
author Jonatan Werpers <jonatan@werpers.com>
date Thu, 26 May 2016 16:45:24 +0200
parents
children
comparison
equal deleted inserted replaced
196:289fd8b1635c 197:4f7930d2d2c4
1 % Calculates a D dimensional p-order bspline at t with knots T and control points P.
2 % T = [t0 t1 t2 ... tm] is a 1 x (m+1) vector with non-decresing elements and t0 = 0 tm = 1.
3 % P = [P0 P1 P2 ... Pn] is a D x (n+1) matrix.
4
5 % knots p+1 to m-p-1 are the internal knots
6
7 % Implemented from: http://mathworld.wolfram.com/B-Spline.html
8 function C = bspline(t,p,P,T)
9 m = length(T) - 1;
10 n = size(P,2) - 1;
11 D = size(P,1);
12
13 assert(p == m - n - 1);
14
15 C = zeros(D,length(t));
16
17 for i = 0:n
18 for k = 1:D
19 C(k,:) = C(k,:) + P(k,1+i)*B(i,p,t,T);
20 end
21 end
22
23 % Curve not defined for t = 1 ? Ugly fix:
24 I = find(t == 1);
25 C(:,I) = repmat(P(:,end),[1,length(I)]);
26 end
27
28 function o = B(i, j, t, T)
29 if j == 0
30 o = T(1+i) <= t & t < T(1+i+1);
31 return
32 end
33
34 if T(1+i+j)-T(1+i) ~= 0
35 a = (t-T(1+i))/(T(1+i+j)-T(1+i));
36 else
37 a = t*0;
38 end
39
40 if T(1+i+j+1)-T(1+i+1) ~= 0
41 b = (T(1+i+j+1)-t)/(T(1+i+j+1)-T(1+i+1));
42 else
43 b = t*0;
44 end
45
46 o = a.*B(i, j-1, t, T) + b.*B(i+1, j-1, t, T);
47 end