How to vectorize for loops?

1 次查看(过去 30 天)
Hi Everybody, I have three for loops and their processing is very slow, I need to speed up the process. For that purpose we need to convert it to vectors . Any help will be strongly encouraged. Below is the code:
for k = 1:size_glcm_3
glcm_sum(k) = sum(sum(glcm(:,:,k)));
glcm(:,:,k) = glcm(:,:,k)./glcm_sum(k); % Normalize each glcm
glcm_mean(k) = mean2(glcm(:,:,k)); % compute mean after norm
glcm_var(k) = (std2(glcm(:,:,k)))^2;
for i = 1:size_glcm_1
for j = 1:size_glcm_2
p_x(i,k) = p_x(i,k) + glcm(i,j,k);
p_y(i,k) = p_y(i,k) + glcm(j,i,k); % taking i for j and j for i
p_xplusy((i+j)-1,k) = p_xplusy((i+j)-1,k) + glcm(i,j,k);
p_xminusy((abs(i-j))+1,k) = p_xminusy((abs(i-j))+1,k) +...
glcm(i,j,k);
end
end
end
All arrays are pre-allocated, size of size_glcm_1 and size_glcm_2 is 512 and size of size_glcm_3 is 1 .
  3 个评论
azizullah khan
azizullah khan 2015-12-12
编辑:azizullah khan 2015-12-12
@Jan Thank you for your response. Yes, i pre-allocated all arrays, both the size of size_glcm_1 and size_glcm_2 is 512 and size_glcm_3 is 1. the time it is taking is due the two internal for loops ( i & j). Waiting for your response. thx
Image Analyst
Image Analyst 2015-12-13
Why not use var() instead of std2() and squaring?

请先登录,再进行评论。

采纳的回答

dpb
dpb 2015-12-13
编辑:dpb 2015-12-15
To simplify, I presume copy the plane slice to a 2D array to eliminate the k index from the expressions. The following duplicated your results for a trial array of random values...
px=sum(glcm,2);
py=sum(glcm,1);
N=size_glcm_1-1;
j=0;
for i=-N:N
j=j+1;
pxp(j)=sum(diag(fliplr(glcm),i));
end
pxm(1)=sum(diag(glcm)); % is only one main diagonal
for i=1:N
pxm(i+1)=sum(diag(glcm),-i)+sum(diag(glcm),i); % +/- off-diagonals
end
NB: You may want a temporary for fliplr(glcm); I'm not sure if the JIT optimizer will avoid doing the operation every pass or not; you can test and adjust as seems necessary.
You can also 'spearmint w/ accumarray and friends to see about eliminating the remaining loops; it wasn't patently apparent to me it would help altho for a fixed size you perhaps could build the needed indexing arrays a priori.
ADDENDUM
OK, the accumarray solution isn't as bad as I thought it might have been...see comments below on "how it works".
Alternate solution--
% the preliminaries...
N=size_glcm_1-1;
[i j]=ind2sub(size(glcm),1:numel(glcm));
idx=i+j-1; idx(end)=1;
% calculations
px=sum(glcm,2);
py=sum(glcm,1);
pxp=accumarray(idx.',glcm);
pxm(1)=sum(diag(glcm)); % is only one main diagonal
for i=1:N
pxm(i+1)=sum(diag(glcm),-i)+sum(diag(glcm),i); % +/- off-diagonals
end
[i j]=ind2sub(size(glcm),1:numel(glcm));
idx=i+j-1; idx(end)=1;
pxp=accumarray(idx.',glcm); % the result
  12 个评论
dpb
dpb 2015-12-14
编辑:dpb 2015-12-15
"... fliplr is time consumable"
Did you follow the suggestion of making the result of the operation a temporary outside the loop?
Else, I'd have to think about it some; it's the way you've defined the subscripts unless you could live with the alternate diagonals in the "normal" direction.
ADDENDUM
Actually, on reflection, there is a rather cute relationship that can be exploited that should be reasonably fast.
[i j]=ind2sub(size(g),1:numel(g))-1; % the indices of the array less one
idx=i+j; idx(end)=1; % the combination of the two, correct last position
pxp=accumarray(idx.',g); % the result
It turns out the "positive" diagonals have corresponding indices whose sum is a constant excepting for the first and last which are combined as the first element; hence the fixup for the last position in the summation.
The index vector can be computed outside any loop, of course, and if they were to always be the same could be precomputed and simply read or passed in saving even the generation time.
One would have to test whether this beats the fliplr or not; I'm guessing as it also eliminates the double call to diag it will...
azizullah khan
azizullah khan 2015-12-16
Thank you @dpb , I didn't test flip outside the loop, No, i tested it outside and result is tremendous, Computataional time is decreased now !!!!!

请先登录,再进行评论。

更多回答(0 个)

类别

Help CenterFile Exchange 中查找有关 Loops and Conditional Statements 的更多信息

Community Treasure Hunt

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

Start Hunting!

Translated by