%-------------------------------------------------------------------------%
%         This software is licensed by Creative Commons BY-NC-SA:         %
%      http://creativecommons.org/licenses/by-nc-sa/3.0/it/legalcode      %
%-------------------------------------------------------------------------%
%
% File: DistanceMatrixCSRBF(A,ctrs,ep)
%
% Goal: forms the distance matrix of two sets of points in R^3
%       for compactly supported radial basis functions, i.e.,
%       DM(i,j) = || datasite_i - center_j ||_2.
%
% Inputs: A:      Mx3 matrix representing a set of M points
%         ctrs:   Nx3 matrix representing a set of N centers for RBFs 
%         ep:     determines size of support of basis function.
%  
% Outputs: DM:    MxN 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: kD-tree package by Pramod Vemulapalli from
%           MATLAB Central File Exchange
%
% Remark:  the original routine comes from the book: 
%          [G. E. Fasshauer, Meshfree approximation methods with Matlab, 
%          World Scientific, Singapore, 2007]. The program has been 
%          modified for the use of the Kd-tree package written by 
%          P. Vemulapalli.
%
%-------------------------------------------------------------------------%
function DM = DistanceMatrixCSRBF(A,ctrs,ep)
N = size(A,1);
M = size(ctrs,1);
% For each center (basis function), find the data sites
% in its support along with u-distance
support = 1/ep;                                     
nzmax = 25*N; rowidx = zeros(1,nzmax); colidx = zeros(1,nzmax);
validx = zeros(1,nzmax); istart = 1; iend = 0;
if M > N % Faster if more centers than data sites
    [Tree] = kd_buildtree(ctrs,[]);
    for i = 1:N
         [idx,dist] = kd_rangequery_ball(Tree,A(i,:),support);
        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
     [Tree] = kd_buildtree(A,[]);
     for j = 1:M
        [idx,dist] = kd_rangequery_ball(Tree,ctrs(j,:),support);
        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);