for loop vectorization of De Boor algorithm

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

Answer 1

[Bruno Luong]

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