function cbpxs = cbpmat( cbpcenters,grid,delta )
% cbpxs = cbpmat( cbpcenters,grid,delta )
%  Generates design matrix used for first-order CBP
%
% Inputs:
%  - cbpcenters, 3xK matrix: initial grid of parameters
%  - grid, 3xN matrix: a grid=(v_1,...,v_n) of directions on the sphere
%     where diffusion measurements are taken
%  - delta, scalar: minimum space of cbpcenters
%
% Outputs:
%  - cbpxs, Nx(8*K) matrix: CBP design matrix
%     Use by calling lsqnonneg(cbpxs, y)
[xsc,d1sc,~] = ste_tan0(cbpcenters,grid);
cbpxs = [xsc+delta.*d1sc{1}+delta.*d1sc{2}+delta.*d1sc{3},...
    xsc+delta.*d1sc{1}+delta.*d1sc{2}-delta.*d1sc{3},...
    xsc+delta.*d1sc{1}-delta.*d1sc{2}+delta.*d1sc{3},...
    xsc+delta.*d1sc{1}-delta.*d1sc{2}-delta.*d1sc{3},...
    xsc-delta.*d1sc{1}+delta.*d1sc{2}+delta.*d1sc{3},...
    xsc-delta.*d1sc{1}+delta.*d1sc{2}-delta.*d1sc{3},...
    xsc-delta.*d1sc{1}-delta.*d1sc{2}+delta.*d1sc{3},...
    xsc-delta.*d1sc{1}-delta.*d1sc{2}-delta.*d1sc{3}];

end

