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 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),...
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(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:
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.
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:
dMap_ = abs(xcorr2(A,h)-xcorr2(h,A))
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.
Here's the diffMap function:
function dMap = diffMap(A)
Rep = zeros(r + r*2-2, c + c*2-2);
Rep_ = zeros(r + r*2-2, c + c*2-2);
Rep(x,y) = A(x-r+1, y-c+1);
Rep_(x,y) = h(x-r+1, y-c+1);
dMap = zeros(r+r-1,c+c-1);
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)));