accessing all elements along the diagonals and rows (simple backprojection algorithm for image reconstruction)
1 view (last 30 days)
Show older comments
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!
2 Comments
Image Analyst
on 19 Oct 2018
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?
Accepted Answer
Guillaume
on 19 Oct 2018
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
9 Comments
Guillaume
on 19 Oct 2018
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
More Answers (3)
Andrei Bobrov
on 19 Oct 2018
Edited: Andrei Bobrov
on 19 Oct 2018
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);
0 Comments
Bruno Luong
on 19 Oct 2018
Edited: Bruno Luong
on 19 Oct 2018
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
0 Comments
See Also
Categories
Find more on Operating on Diagonal Matrices in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!