Clear Filters
Clear Filters

for loop vectorization of De Boor algorithm

11 views (last 30 days)
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
% Copyright 2010 Levente Hunyadi
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

Answers (1)

Bruno Luong
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')

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!