accessing all elements along the diagonals and rows (simple backprojection algorithm for image reconstruction)

Hi, I'm wondering how to effectively sum all elements along every diagonal and row that touches an element in a given matrix. For example, in the matrix below

1 2 3

4 5 6

7 8 9

I would like to compute 1 + (1+2+3)/3 + (1+5+9)/3 + (1+4+7)/3. This value should replace the element 1 at (1,1). To replace the element 2 at (1,2), it would be (2+4)/2 + (2+6)/2 + (1+2+3)/3 + (2+5+8)/3. Similarly, element 6 at (2,3) should be (2+6)/2 + (6+8)/2 + (4+5+6)/3 + (3+6+9)/3.

It's a bit like a Sudoku calculation. The algorithm needs to be flexible for any square matrix, of any size (my data set is 360x360). Any help is much appreciated!

[KSSV]

You can access diagonal elements using diag. Row elements by A(i,:).

[Image Analyst]

For a360x360 matrix, does the sum just go over the 3 nearest elements inside the matrix, or does it go all the way to the outer edge of the matrix?

Explain why you want to do this thing. What's the use case?

Answer 1

[Guillaume]

You can indeed use convolutions:

m = reshape(1:100, 10, 10);  %demo matrix

convlength = 2 * size(m, 1);  %also 2 * size(m, 2) since m is square
%size of convolution matrices is twice the size of the input matrix to make sure every element is  always taken into account
rowmean = conv2(m, ones(1, convlength) / convlength * 2, 'same');
colmean = conv2(m, ones(convlength, 1) / convlength * 2, 'same');
diagsum = conv2(m, eye(convlength), 'same');
diagcount = toeplitz(size(m, 1):-1:1);
diagmean = diagsum ./ diagcount;
antidiagsum = conv2(fliplr(m), eye(convlength), 'same');
antidiagmean = fliplr(antidiagsum ./ diagcount);
result = rowmean + colmean + diagmean + antidiagmean

[Bruno Luong]

I think diagcount must be different for diag and antidiag.

Shouldn't flip m but eye for antidiag.

Not sure either why you multiply by 2 for row and col.

[Guillaume]

I think diagcount must be different for diag and antidiag.

Shouldn't flip m but eye for antidiag.

Since I flip m, the antidiagonals are now diagonals, and hence diagcount still apply. Note that the whole mean matrix is then flipped back. You could flip eye instead indeed, and use a flipped diagcount.

Not sure either why you multiply by 2 for row and col.

The convolution matrix is twice as big as the input matrix (to make sure the convolution covers all elements at the edges). However, for the mean you want to divide by the number of elements, which is convlengh/2 (or multiply by 2/convlength). Maybe it'd be clearer if I'd written:

msize = size(m, 1);
rowmean = conv2(m, ones(1, 2*msize) / msize, 'same');

However, I was so focussed on convolutions that I missed that both rowmean and colmean could be calculated more easily:

rowmean = repmat(mean(m, 2), 1, size(m, 2));
colmean = repmat(mean(m, 1), size(m, 1), 1));

[Bruno Luong]

Since I flip m, the antidiagonals are now diagonals, and hence diagcount still apply. Note that the whole mean matrix is then flipped back. You could flip eye instead indeed, and use a flipped diagcount.

But when you combine later the results (of antidia with the rest) you combine the flip one with the right one.

[Guillaume]

But when you combine later the results (of antidia with the rest) you combine the flip one with the right one

I'm not sure I understand. I may very well have made a mistake, I haven't tested this very thoroughly. However, for the example matrix in the question, my algorithm produces the desired result at the 3 example points.

[Bruno Luong]

conv2(flip(A),K,'same')

is equal to

flip(conv2(A,flip(K),'same')

For the correct anti-diagonal result, OP wants

K = eye(...)
conv2(A,flip(K),'same')

and not the flip of it.

[Guillaume]

Hum, no, the sum of the antidiagonals, replicated along the antidiagonals is

flip(conv2(flip(A), K, 'same'))

or am I missing something?

[Andrei Bobrov]

+1. Nice idea!

[Guillaume]

Andrei's answer reminded me of spdiag which is another way of calculating the diagonal and antidiagonal means. So the whole thing can be calculated without convolutions:

m = reshape(1:100, 10, 10);  %demo matrix

rowmean = repmat(mean(m, 2), 1, size(m, 2));
colmean = repmat(mean(m, 1), size(m, 1), 1);
diagmean = sum(spdiags(m)) ./ sum(spdiags(ones(size(m))));
diagmean = diagmean(toeplitz(size(m, 1):-1:1, size(m, 1):size(m, 1)+size(m, 2)-1));
antidiagmean = sum(spdiags(flip(m))) ./ sum(spdiags(ones(size(m))));
antidiagmean = antidiagmean(hankel(1:size(m, 1), size(m, 1):size(m, 1)+size(m, 2)-1));
result = rowmean + colmean + diagmean + antidiagmean

[Bruno Luong]

Oh OK I didn't see you flip it back, my bad

Answer 2

[Image Analyst]

You can do this with conv2(). If you can't figure it out, let me know.

[kvk]

I'm not familiar using conv2() here. Could you please elaborate more? What is convolving with what? Thank you.

I'm trying to use another approach other than diag, as it seems to be inefficient accessing elements along every diagonal for every element in the 360x360 matrix.

Answer 3

[Andrei Bobrov]

A = [1 4 7 10 13;
     2 5 8 11 14;
     3 6 9 12 15];

s = size(A);
n = ones(s(1),1);
a = n*sum(A)/s(1);
b = sum(A,2)*ones(1,s(2))/s(2);
k = sum(triu(rot90(triu(ones(s + [0, s(1)-1])),2)));
c = sum(spdiags(A))./k;
d = sum(spdiags(rot90(A)))./k;
dlr = full(spdiags(n*c,1-s(1):s(2)-1,s(1),s(2)));
drl = flip(full(spdiags(n*d,1-s(1):s(2)-1,s(1),s(2))),1);
out = a + b + dlr + drl;

variant by Guillaume's idea

s = size(A);
n = min(s)*2-1;
o = ones(s);
Iud = eye(n);
Idu = flip(Iud,2);
nn = conv2(o,Iud,'same');
out1 = conv2(A,Iud,'same')./nn + ...
       conv2(A,Idu,'same')./flip(nn,2) + ...
       o(:,1)*sum(A)/s(1) + sum(A,2)*o(1,:)/s(2);

Answer 4

[Bruno Luong]

OK, here is the conv method as proposed by IA and started by Guillaume, it works for rectangular input array as well. Some convolution kernels are computed larger than necessary for clarity.

A = [1 2 3
     4 5 6
     7 8 9];

[m,n] = size(A);
p = 2*m-1;
q = 2*n-1;
r = (1:p)';
c = (1:q)';
d = min(m,n)-1;
k = (-d:d)';
Kh = accumarray([m+0*c, c],1,[p,q]);
Kv = accumarray([r, n+0*r],1,[p,q]);
Kd = accumarray([m+k, n+k],1,[p,q]);
Ka = accumarray([m+k, n-k],1,[p,q]);

I = ones(size(A));
cfun = @(K) conv2(A,K,'same')./conv2(I,K,'same');
% alternative way of calculation of h and v:
% [h,v] = ndgrid(mean(A,2),mean(A,1))
h = cfun(Kh);
v = cfun(Kv);
d = cfun(Kd);
a = cfun(Ka);

Result = h+v+d+a