Efficient computation of the sum of pairwise absolute differences

7 次查看（过去 30 天）
Santiago Benito
Santiago Benito 2021-10-14
评论： Santiago Benito 2021-10-15
Ran in:
Hi,
I am working on a project that requires the computation of the sum of all pairwise absolute differences between elements at either end of a randomly placed vector with coordinates . For instance, starting with the following 3x3 matrix:
A = magic(3)
A = 3×3
8 1 6 3 5 7 4 9 2
A vector with coordinates will produce the absolute differences:
b = sum(abs([A(1,1)-A(1,2),A(1,2)-A(1,3),A(2,1)-A(2,2),A(2,2)-A(2,3),...
A(3,1)-A(3,2),A(3,2)-A(3,3)]))
b = 28
Note that I took the differences between all "horizontal" matrix elements separated by a distance of .
The idea is to save this value in a matrix, let's say, C, in the position corresponding to the vector coordinate differences, i.e., with respect to the center of the matrix A. Thus, if we were to do the same with the vector defined by , we get:
C = zeros(size(A));
C(1,3) = sum(abs([A(2,1)-A(1,2),A(2,2)-A(1,3),A(3,1)-A(2,2),A(3,2)-A(2,3)]));
And adding in the previously calculated value:
C(2,3) = b
C = 3×3
0 0 6 0 0 28 0 0 0
As a proof of concept, I borrowed and adapted a function from another answer that does the computation. I have called it diffMap and is at the bottom of this question/script.
dMap = diffMap(A)
dMap = 5×5
6 2 16 6 2 2 18 20 6 6 8 28 0 28 8 6 6 20 18 2 2 6 16 2 6
As one can probably easily notice, this is just a modification of a rudimentary for-loopy cross-correlation.
I have been giving some thought to this problem in the hope of making this work with xcorr2 and discrete Fourier transforms. This is what I have got until now:
h = ones(size(A));
dMap_ = abs(xcorr2(A,h)-xcorr2(h,A))
dMap_ = 5×5
6 2 0 6 2 2 6 0 2 6 0 0 0 0 0 6 2 0 6 2 2 6 0 2 6
It would be great if this worked, because it could be easily adapted to work with fft2, ifft2, and fftshift thanks to the convolution theorem. However, this does not return what i was hoping for: it computes the absolute value of the sum of the differences and NOT the sum of the absolute differences.
The question: is there any way to make this work efficiently, ideally with conv2 or xcorr2? With big matrices four for-loops takes forever, obviously.
Thanks!!
Here's the diffMap function:
%% Function
function dMap = diffMap(A)
%DIFFMAP Compute the absolute differences map
% Get the size of the matrix A
[r,c] = size(A);
% Initialize and index matrix
h = ones([r,c]);
% Padding
Rep = zeros(r + r*2-2, c + c*2-2);
Rep_ = zeros(r + r*2-2, c + c*2-2);
for x = r : r+r-1
for y = c : c+r-1
Rep(x,y) = A(x-r+1, y-c+1);
Rep_(x,y) = h(x-r+1, y-c+1);
end
end
dMap = zeros(r+r-1,c+c-1);
for x = 1 : r+r-1
for y = 1 : c+c-1
for i = 1 : r
for j = 1 : c
dMap(x, y) = dMap(x, y) + abs((Rep(x+i-1, y+j-1) * h(i, j)) -...
(Rep_(x+i-1, y+j-1) * A(i, j)));
end
end
end
end
end

请先登录，再回答此问题。

采纳的回答

Matt J
Matt J 2021-10-15
编辑：Matt J 2021-10-15
I think parallelization might be the only way to accelerate things. Parfor seems to work well with the rewritten version of diffMap() below. With 12 cores, I can do a 250x250 matrix in about 2.7 sec., whereas the original diffMap code takes about 60 sec.
A = magic(250);
tic;
dMap = diffMap(A);
toc
%% Function
function dMap = diffMap(A)
%DIFFMAP Compute the absolute differences map
[r,c] = size(A);
[D1,D2]=deal(zeros(r,c));
parfor k=1:r*c
[i,j]=ind2sub([r,c],k);
if i==1 && j==1, continue; end
kernel=zeros(i,j);
kernel([1,end])=[-1,1];
D1(k)=sum(abs(conv2(A,kernel,'valid')),'all');
D2(k)=sum(abs(conv2(A,flipud(kernel),'valid')),'all');
end
dMap=[flipud(D2);D1];
dMap(end/2,:)=[];
dMap=[rot90(dMap,2),dMap];
dMap(:,end/2)=[];
end
  1 个评论
Santiago Benito
Santiago Benito 2021-10-15
Thank you very much for this version of the code! I tried replacing conv2 with conv2fft by Bruno Luong, but that actually takes longer. Probably because the matrices for the convolution computation in the initial iterations of the parfor loop are really small, or because that function is not really parallelizable.
In my current setup, a 480x680 matrix takes around 340 seconds to compute with conv2 and with conv2fft it took so long, that I had to stop it after 4000 seconds. For more info, I have four cores (AMD Ryzen 3 3200G).

请先登录，再进行评论。

更多回答（2 个）

Matt J
Matt J 2021-10-14
编辑：Matt J 2021-10-15
The question: is there any way to make this work efficiently, ideally with conv2 or xcorr2?
If you were taking the sum of squared differences, it could be done with xcorr2 in just a few lines:
W=ones(size(A));
tmp=xcorr2(A.^2,W);
dMap=tmp+rot90(tmp,2)-2*xcorr2(A);
However, because you are taking the sum of absolute differences, there is no simple connection with xcorr2.
With big matrices four for-loops takes forever, obviously.
Is there an upper bound on the size of the matrix that you will need to process?
  1 个评论
Santiago Benito
Santiago Benito 2021-10-15
Hello Matt, thank you very much for all of your answers, I really appreciate it!
I just read your comments and functions and tried some things on my own. I will post comments to your answers separately.
Firstly, I have at the moment matrices ranging from 480x640 to around 1000x800; they are relatively big for this kind of thing. 1200x1200 might be as big as it gets in this case.
Yeah, I wish I could stick to squared differences. I have that working relatively quickly already. I was trying to find a way to either:
  • Compute the absolute differences independently but also using convolutions.
  • Compute the absolute differences as a derivation of the squared differences map.
Arguably, these two options can be seen as the same thing, though.

请先登录，再进行评论。

Matt J
Matt J 2021-10-15
This GPU implementation may also be useful. I was able to process a 500x500 matrix in 20 seconds on the GTX 1080 Ti.
A = magic(500);
gd=gpuDevice;
tic;
dMap = diffMap(A);
wait(gd);
toc
%% Function
function dMap = diffMap(A)
%DIFFMAP Compute the absolute differences map
% Get the size of the matrix A
[r,c] = size(A);
A=gpuArray(A);
T=gpuArray(toeplitz(1:r,[1,zeros(1,r-1)]));
t=nonzeros(T);
W=reshape(logical(T),r,1,[]);
R=double(W);
% Initialize and index matrix
% Padding
dMap = gpuArray.zeros(r,2*c-1);
for j=1:c
col=A(:,j);
R(W)=col(t);
tmp=permute(sum( abs(R-A).*W,1) ,[3,2,1]);
%ctr=c
dMap(:,c-(j-1):2*c-j) = dMap(:,c-(j-1):2*c-j) +tmp;
end
dMap=[rot90(dMap,2);dMap];
dMap(end/2,:)=[];
end
  1 个评论
Santiago Benito
Santiago Benito 2021-10-15
Again, thanks for the help! This function is, because of my setup, slower for me than the parallel one; I only have a GT 710. Just for reference, it takes 610 seconds to run, a little less than double the time of the parfor version.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Startup and Shutdown 的更多信息

标签

产品


版本

R2019b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by