%-------------------------------------------------------------------------%
%         This software is licensed by Creative Commons BY-NC-SA:         %
%      http://creativecommons.org/licenses/by-nc-sa/3.0/it/legalcode      %
%-------------------------------------------------------------------------%
%
% file BlockBased3D_DistanceMatrix(dsites,rbfctrs,ep)
%
% Goal: forms the distance matrix of two sets of points in R^3 for
%       compactly supported RBFs, i.e.,
%       DM(i,j) = || datasite_i - centre_j ||_2.
%
% Inputs: dsites:   NX3 matrix representing a set of N data sites
%         rbfctrs:  MX3 matrix representing a set of M centres for RBFs
%         ep:       shape parameter
%
% Outputs: DM: NXM SPARSE matrix that contains the Euclidean u-distance
%              (u=1-e*r) between the i-th data site and the j-th center
%              in the i,j position
%
% Calls on: BlockBased3D_Structure, BlockBased3D_ContainingQuery, 
%           BlockBased3D_RangeSearch
%
% Remark:  The original routine comes from the book: 
%          [G. E. Fasshauer, Meshfree approximation methods with Matlab, 
%          World Scientific, Singapore, 2007].
%
%-------------------------------------------------------------------------%
function DM = BlockBased3D_DistanceMatrix(dsites,rbfctrs,ep)
N = size(dsites,1); M = size(rbfctrs,1); rbfradius = 1/ep;                                     
nzmax = 25*N; rowidx = zeros(1,nzmax); colidx = zeros(1,nzmax);
validx = zeros(1,nzmax); istart = 1; iend = 0;
max_clmns_ds = max(dsites,[],1); min_clmns_ds = min(dsites,[],1);
min_dsites = min(min_clmns_ds); max_dsites = max(max_clmns_ds);
max_clmns_rb = max(rbfctrs,[],1); min_clmns_rb = min(rbfctrs,[],1);
min_rbfctrs = min(min_clmns_rb); max_rbfctrs = max(max_clmns_rb);
min_box = min(min_rbfctrs,min_dsites);
max_box = max(max_rbfctrs,max_dsites);
% Define the number of square cells in one direction
q = ceil((max_box-min_box)./(rbfradius));
if M > N % Faster if more centres than data sites
    % Build cell-based partitioning structure
    if q>2
        idx_ct = BlockBased3D_Structure(rbfctrs,q,rbfradius,min_box);
    end
    for i = 1:N
        if q>2
            % Find the square cell containing the i-th centre
            index = BlockBased3D_ContainingQuery(dsites(i,:),...
                q,rbfradius,min_box);
            % For each centre (basis function), find the data sites in its
            % support
            [idx dist] = BlockBased3D_RangeSearch(dsites(i,:),rbfradius,...
                rbfctrs(idx_ct{index},:),idx_ct{index});
        else
            [idx dist] = BlockBased3D_RangeSearch(dsites(i,:),rbfradius,...
                rbfctrs);
        end
        newentries = length(idx); iend = iend + newentries;
        rowidx(istart:iend) = repmat(i,1,newentries);
        colidx(istart:iend) = idx'; validx(istart:iend) = 1-ep*dist';
        istart = istart + newentries;
    end
else
    if q>2
        % Build cell-based partitioning structure
        idx_ds = BlockBased3D_Structure(dsites,q,rbfradius,min_box);
    end
    for j = 1:M
        if q>2
            % Find the square cell containing the j-th centre
            index = BlockBased3D_ContainingQuery(rbfctrs(j,:),q,...
                rbfradius,min_box);
            % For each centre, find the data sites in its support
            [idx dist] = BlockBased3D_RangeSearch(rbfctrs(j,:),...
                rbfradius,dsites(idx_ds{index},:),idx_ds{index});
        else
            [idx dist] = BlockBased3D_RangeSearch(rbfctrs(j,:),...
                rbfradius,dsites);
        end
        newentries = length(idx); iend = iend + newentries;
        rowidx(istart:iend) = idx';
        colidx(istart:iend) = repmat(j,1,newentries);
        validx(istart:iend) = 1-ep*dist'; istart = istart + newentries;
    end
end
idx = find(rowidx); DM = sparse(rowidx(idx),colidx(idx),validx(idx),N,M);