MATLAB: How to vectorize for loops

for loopvectorization

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 .

Best Answer

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
Related Question