function allNeighbors = computeNeighbors(V,dist,nPerBatch,distinct)

nTotal = size(V,1);
nBatches = ceil(nTotal/nPerBatch);
start = 1;
neighbors = cell(1,nTotal); maxNeighbors = -inf;
for b=1:nBatches
[b nBatches]
n = min(nPerBatch,nTotal-start+1);
inds = start:(start+n-1);
neighborhood = repmat(sum(V(inds,:).^2,2),1,nTotal) - 2*V(inds,:)*V' + repmat(sum(V.^2,2),1,n)' <= dist;
for i=1:length(inds)
   neighbors{inds(i)} = find(neighborhood(i,:));
   if distinct
      neighbors{inds(i)}(neighbors{inds(i)} <= inds(i)) = [];
   else
      neighbors{inds(i)}(neighbors{inds(i)} == inds(i)) = [];
   end
   if length(neighbors{inds(i)}) > maxNeighbors
      maxNeighbors = length(neighbors{inds(i)});
   end
end
start = start + n;
clear neighborhood inds n inds;
end

allNeighbors = -1*ones(length(neighbors),maxNeighbors);
for i=1:length(neighbors)
   allNeighbors(i,1:length(neighbors{i})) = neighbors{i};
end

