# for loop vectorization of De Boor algorithm

11 views (last 30 days)
Maximilian Hoffmann on 4 Jul 2017
Edited: Bruno Luong on 22 Mar 2024
Hi there, I'm currently using Levente Hunyadi's DeBoor algorithm to determine 1000 Points on a 1D-B-Spline-Curve. This is my script:
if true
k=5;
u=[0 0 0 0 0 0.125 0.25 0.375 0.5 0.625 0.75 0.875 1 1 1 1 1];
d=[0 1 2 3 4 5 6 7 8 9 10 11;
0 0.06 0.2 0.3 0.4 0.5 0.55 0.6 0.7 0.8 0.9 1];
tic;
[C,U]=DeBoor(k,u,d,1000);
dura1=toc;
disp(dura1);
end
And this is his code:
if true
function [C,U] = DeBoor(n,t,P,U)
% Evaluate explicit B-spline at specified locations.
%
% Input arguments:
% n:
% B-spline order (2 for linear, 3 for quadratic, etc.)
% t:
% knot vector
% P:
% control points, typically 2-by-m, 3-by-m or 4-by-m (for weights)
% u (optional):
% values where the B-spline is to be evaluated, or a positive
% integer to set the number of points to automatically allocate
%
% Output arguments:
% C:
% points of the B-spline curve
validateattributes(n, {'numeric'}, {'positive','integer','scalar'});
d = n-1; % B-spline polynomial degree (1 for linear, 2 for quadratic, etc.)
validateattributes(t, {'numeric'}, {'real','vector'});
assert(all( t(2:end)-t(1:end-1) >= 0 ), 'bspline:deboor:InvalidArgumentValue', ...
'Knot vector values should be nondecreasing.');
validateattributes(P, {'numeric'}, {'real','2d'});
nctrl = numel(t)-(d+1);
assert(size(P,2) == nctrl, 'bspline:deboor:DimensionMismatch', ...
'Invalid number of control points, %d given, %d required.', size(P,2), nctrl);
if nargin < 4
U = linspace(t(d+1), t(end-d), 10*size(P,2)); % allocate points uniformly
elseif isscalar(U) && U > 1
validateattributes(U, {'numeric'}, {'positive','integer','scalar'});
U = linspace(t(d+1), t(end-d), U); % allocate points uniformly
else
validateattributes(U, {'numeric'}, {'real','vector'});
assert(all( U >= t(d+1) & U <= t(end-d) ), 'bspline:deboor:InvalidArgumentValue', ...
'Value outside permitted knot vector value range.');
end
m = size(P,1); % dimension of control points
t = t(:).'; % knot sequence
U = U(:);
S = sum(bsxfun(@eq, U, t), 2); % multiplicity of u in t (0 <= s <= d+1)
I = bspline_deboor_interval(U,t);
Pk = zeros(m,d+1,d+1);
a = zeros(d+1,d+1);
C = zeros(size(P,1), numel(U));
for j = 1 : numel(U)
u = U(j);
s = S(j);
ix = I(j);
Pk(:) = 0;
a(:) = 0;
% identify d+1 relevant control points
Pk(:, (ix-d):(ix-s), 1) = P(:, (ix-d):(ix-s));
h = d - s;
if h > 0
% de Boor recursion formula
for r = 1 : h
q = ix-1;
for i = (q-d+r) : (q-s);
a(i+1,r+1) = (u-t(i+1)) / (t(i+d-r+1+1)-t(i+1));
Pk(:,i+1,r+1) = (1-a(i+1,r+1)) * Pk(:,i,r) + a(i+1,r+1) * Pk(:,i+1,r);
end
end
C(:,j) = Pk(:,ix-s,d-s+1); % extract value from triangular computation scheme
elseif ix == numel(t) % last control point is a special case
C(:,j) = P(:,end);
else
C(:,j) = P(:,ix-d);
end
end
function ix = bspline_deboor_interval(u,t)
% Index of knot in knot sequence not less than the value of u.
% If knot has multiplicity greater than 1, the highest index is returned.
i = bsxfun(@ge, u, t) & bsxfun(@lt, u, [t(2:end) 2*t(end)]); % indicator of knot interval in which u is
[row,col] = find(i);
[row,ind] = sort(row); %#ok<ASGLU> % restore original order of data points
ix = col(ind);
end
I would like to know, if there's a way to change the for loop with "for j = 1 : numel(U)" by vectorizing it and replace the (slow) if-else-conditions using logical indexing, making the code faster for a large number of points. I completely failed to do this after 12 hours of trying... I'm really looking forward to hearing from you. Regards, Max

Bruno Luong on 21 Mar 2024
Edited: Bruno Luong on 22 Mar 2024
Old question, but I'll reply it for future readers that wonder how to tackle the same problem.
I did a quick and dirdty vectorization of DeBoor.m extracted from this FEX the speed up is not significant (therefore the vectorized code is not provided)
However I have another code that I took from my FEX BSFK named Bernstein.m that can do the same evaluation of interpolation spline. For middle size of query points, 1000 here, the speed up is about 4-7 times.
The algorithm in Bernstein.m is not the so-Called "De Boor algoritm" but based on Cox-Deboor formula on "non-normalized" B-spline, from Carl De Boor 1972 paper. This primitive version seems to be more suiitable for MATLAB vectorization.
EDIT 22/Mar/2024: new version of Bernstein uploaded
k=4;
uu=cumsum(rand(1,10));
[a,b] = bounds(uu);
u=[a+zeros(1,k-1), uu, b+zeros(1,k-1)];
d=rand(1,length(uu)+k-2);
[C,U]=DeBoor(k,u,d,1000);
t0=timeit(@()DeBoor(k,u,d,1000),1);
% Vectorized DeBoor
% B = DBoorBS(U, u, [], k, d);
% t1=timeit(@()DBoorBS(U, u, [], k, d),1);
% t1r = t1/t0;
D = Bernstein(U, u, [], k, d);
t2=timeit(@()Bernstein(U, u, [], k, d),1);
t2r = t2/t0;
if t2r > 1
fprintf("Bernstein is %g time slower than DeBoor\n", t2r)
else
fprintf("Bernstein is %g time faster than DeBoor\n", 1/t2r)
end
Bernstein is 6.92918 time faster than DeBoor
close all
plot(uu, d((k-2)+(1:length(uu))),'ro')
hold on
plot(U, C','k.')
plot(U, D, '-g')
legend('Conrol points', 'DeBoor', 'Bernstein', ...
'Location', 'best')