function [D,V] = sdistm(X,Y,r,v) % SDISTM create sparse distance matrix. % D = SDISTM(X,Y,r) is the (sparse) matrix of pairwise euclidean distances % between the points in the rows of X and Y not exceeding distance r. % % [D,V] = SDISTM(X,Y,r,v) additionally returns the sparse matrix V with entries % V_ij = v(Y(j,:),X(i,:)). The sparsity structure of V and D are the same. % % Uses the Matlab kdtree package by Guy Shechter. N = size(X,1); M = size(Y,1); SD = zeros(3,25*max(M,N)); % row index, column index and entry of distance matrix if nargin == 4 SV = zeros(3,25*max(M,N)); end if N <= M, % distinction for efficiency reasons j = 1; [~,~,t] = kdtree(Y,[]); for i = 1:N [~,dist,jdx] = kdrangequery(t,X(i,:),r); dj = length(jdx); % + realmin in order to preserve sparsity structure SD(:,j:j+dj-1) = [i*ones(1,dj); jdx'; (dist+realmin)']; if nargin == 4 SV(:,j:j+dj-1) = [i*ones(1,dj); jdx'; v(Y(jdx,:), X(i*ones(1,dj),:))' + realmin]; end j = j + dj; end else i = 1; [~,~,t] = kdtree(X,[]); for j = 1:M [~,dist,idx] = kdrangequery(t,Y(j,:),r); di = length(idx); % + realmin in order to preserve sparsity structure SD(:,i:i+di-1) = [idx'; j*ones(1,di); (dist+realmin)']; if nargin == 4 SV(:,i:i+di-1) = [idx'; j*ones(1,di); v(Y(j*ones(1,di),:), X(idx,:))' + realmin]; end i = i + di; end end I = find(SD(1,:)); D = sparse(SD(1,I),SD(2,I),SD(3,I),N,M); if nargin == 4 I = find(SV(1,:)); V = sparse(SV(1,I),SV(2,I),SV(3,I),N,M); end kdtree([],[],t);