0

I have a 3 X 1000 (and later 3 X 10 000) matrix cord given, which contains the three dimensional coordinates for my pixels.

My intention is to calculate the distance between all the pixels, and I do it with a for loop (see below), but I will have to calculate this for huge matrices soon, and am wondering if I could vectorize the code for making it faster...?

dist = zeros(size(cord,2),size(cord,2)); 
for i = 1:size(cord,2)
    for j = 1:size(cord,2)
        dist(i,j) = norm(cord(:,i)-cord(:,j));
        dist(j,i) = dist(i,j);
    end 
end
Luis Mendo
  • 110,752
  • 13
  • 76
  • 147
Pugl
  • 432
  • 2
  • 5
  • 22

3 Answers3

3

pdist does exactly that. squareform is needed to get the result in the form of a square, symmetric matrix:

dist = squareform(pdist(cord.'));
Luis Mendo
  • 110,752
  • 13
  • 76
  • 147
1

Approach 1 (Vectorized apprach with bsxfun ) -

squeeze(sqrt(sum(bsxfun(@minus,cord,permute(cord,[1 3 2])).^2)))

Not sure if this will be faster though.

Approach 2 -

Inspired by this very smart approach and all credits to the poster. The code posted here is just slightly customized for your case and hopefully slightly better in terms of runtime. Here it is -

A = cord'; %//'
numA = size(cord,2);
helpA = ones(numA,9);
helpB = ones(numA,9);
for idx = 1:3
    sqA_idx = A(:,idx).^2;
    helpA(:,3*idx-1:3*idx) = [-2*A(:,idx), sqA_idx ];
    helpB(:,3*idx-2:3*idx-1) = [sqA_idx , A(:,idx)];
end
dist1 = sqrt(helpA * helpB'); %// desired output
Community
  • 1
  • 1
Divakar
  • 218,885
  • 19
  • 262
  • 358
0

From your code, you have recognized that the dist matrix is symmetric

dist(i,j) = norm(cord(:,i)-cord(:,j));
dist(j,i) = dist(i,j);

You could change the inner loop to account for this and reduce by roughly one half the number of calculations needed

for j = i:size(cord,2)

Further, we can avoid the dist(j,i) = dist(i,j); at each iteration and just do that at the end by extracting the upper triangle part of dist and adding its transpose to the dist matrix to account for the symmetry

dist = zeros(size(cord,2),size(cord,2)); 
for i = 1:size(cord,2)
    for j = i:size(cord,2)
        dist(i,j) = norm(cord(:,i)-cord(:,j));
    end 
end
dist = dist + triu(dist)';

The above addition is fine since the main diagonal is all zeros.

It still performs poorly though and so we should take advantage of vectorization. We can do that as follows against the inner loop

dist = zeros(size(cord,2),size(cord,2)); 
for i = 1:size(cord,2)
    dist(i,i+1:end) = sum((repmat(cord(:,i),1,size(cord,2)-i)-cord(:,i+1:end)).^2);
end
dist = dist + triu(dist)';
dist = sqrt(dist);

For every element in cord we need to calculate its distance with all other elements that follow it. We reproduce the element with repmat so that we can subtract it from every element that follows without the need for the loop. The differences are squared and summed and assigned to the dist matrix. We take care of the symmetry and then take the square root of the matrix to complete the norm operation.

With tic and toc, the original distance calculation with a random cord (cord = rand(3,num);) took ~93 seconds. This version took ~2.8.

Geoff
  • 1,603
  • 11
  • 8