function kd = kvdiags(kernel, p, N, width)
%KVDIAGS Arranges kernel in diagonals for efficient output computation.
%
%   KD = KVDIAGS(KERNEL,P,N,WIDTH) arranges in diagonals the coefficients 
%   of a vectorized Pth-order kernel of memory length N.
%   Only diagonals within a maximum distance of WIDTH to the main diagonal 
%   are kept in the new representation. With KD = KVDIAGS(KERNEL,P,N), 
%   a maximum width N is assumed.
%
%   See also VFDIAGS, KERNDIAGS
%
%   CBVK Toolbox
%   JHMG, PMSB 2016


if p == 1
    kd{1} = kernel;
    return;
end;

% Default width 
if nargin < 4
    width = N;
end

kd = cell(1,nchoosek(width+p-2,p-1));
wp = zeros(1,p-1);

indsv = zeros(N,1);

for n=1:length(kd)
    
    Np = N - wp(end);
    
    if n == 1
        indsv(1) = 1;
        for k = 2:N
           indsv(k) = nchoosek(N+p-k,p-1);
        end
        indsv = cumsum(indsv);
    else
        indsv(1:Np) = indsv(1:Np)+1;
    end
    
    kd{n} = kernel(indsv(1:Np));
    kd{n} = kd{n}(:);
    
    wp = nextinds(wp,width);    
end
