Optimizing nested loop for calculating geodesic distance I segmented an image into N superpixels and constructed a graph based on that in which each superpixel is considered as a node. The information about neighbor superpixels is stored in glcms array. The weight between each pair of neighboring superpixels is stored in matrix W. 

Finally, I want to compute geodesic distance between non-adjacent superpixels using graphshortestpath function. The following code does the mentioned process, however it is very time-consuming. Specifically, the last section of the code that computes geodesic distance takes more time than expected (more than 15 sec). 

    Img=imread('input.jpg');
    [rows, columns, numberOfColorChannels] = size(Img);
    [L,N] = superpixels(Img,250);

    %Identifying neighborhood relationships
    glcms = graycomatrix(L,'NumLevels',N,'GrayLimits',[1,N],'Offset',  [0,1;1,0]); %Create gray-level co-occurrence matrix from image
    glcms = sum(glcms,3);    % add together the two matrices
    glcms = glcms + glcms.'; % add upper and lower triangles together, make it symmetric
    glcms(1:N+1:end) = 0;    % set the diagonal to zero, we don't want to see "1 is neighbor of 1"

    data = zeros(N,3);
    for labelVal = 1:N
    redIdx = idx{labelVal};
    greenIdx = idx{labelVal}+numRows*numCols;
    blueIdx = idx{labelVal}+2*numRows*numCols;
    data(labelVal,1) = mean(Img(redIdx));
    data(labelVal,2) = mean(Img(greenIdx));
    data(labelVal,3) = mean(Img(blueIdx));

    end    

   Euc=zeros(N);
% Euclidean Distance
for i=1:N
    for j=1:N
        if glcms(i,j)~=0
            Euc(i,j)=sqrt(((data(i,1)-data(j,1))^2)+((data(i,2)-data(j,2))^2)+((data(i,3)-data(j,3))^2));
        end
    end
end


W=zeros(N);
W_num=zeros(N);

W_den=zeros(N);
OMG1=0.1;
for i=1:N
    for j=1:N
        if(Euc(i,j)~=0)
         W_num(i,j)=exp(-OMG1*(Euc(i,j)));
         W_den(i,i)=W_num(i,j)+W_den(i,i);
        end
    end
end

for i=1:N
    for j=1:N
         if(Euc(i,j)~=0)
             W(i,j)=(W_num(i,j))/(W_den(i,i));   % Connectivity Matrix W

         end
    end
end


s_star_temp=zeros(N);   %temporary variable for geodesic distance measurement
W_sparse=zeros(N);
W_sparse=sparse(W);
for i=1:N
    for j=1:N
        if W(i,j)==0 &amp; i~=j;
            s_star_temp(i,j)=graphshortestpath(W_sparse,i,j); % Geodesic Distance
        end
    end
end




The question is how to optimize the code to be more efficient , i.e. less time-consuming.
I'll first go through the different section and outline possible improvements, then give some general comments at the end. I am not checking the algorithm for correctness, as I don't know it.

Img=imread('input.jpg');
[rows, columns, numberOfColorChannels] = size(Img);
[L,N] = superpixels(Img,250);

%Identifying neighborhood relationships
glcms = graycomatrix(L,'NumLevels',N,'GrayLimits',[1,N],'Offset',  [0,1;1,0]); %Create gray-level co-occurrence matrix from image
glcms = sum(glcms,3);    % add together the two matrices
glcms = glcms + glcms.'; % add upper and lower triangles together, make it symmetric
glcms(1:N+1:end) = 0;    % set the diagonal to zero, we don't want to see "1 is neighbor of 1"


You need to have defined idx by now, but you haven't, so the code doesn't run. Based on the superpixels documentation, you need to add idx = label2idx(L);

data = zeros(N,3);
for labelVal = 1:N
    redIdx = idx{labelVal};
    greenIdx = idx{labelVal}+numRows*numCols;
    blueIdx = idx{labelVal}+2*numRows*numCols;
    data(labelVal,1) = mean(Img(redIdx));
    data(labelVal,2) = mean(Img(greenIdx));
    data(labelVal,3) = mean(Img(blueIdx));
end    


This is fine, it probably is fine to just write data(labelVal,1) = mean(Img(idx{1:N})), etc. but the performance should be similar.

Euc = zeros(N);
% Euclidean Distance
for i=1:N
    for j=1:N
        if glcms(i,j)~=0
            Euc(i,j)=sqrt(((data(i,1)-data(j,1))^2)+((data(i,2)-data(j,2))^2)+((data(i,3)-data(j,3))^2));
        end
    end
end


This can be replaced by Euc = pdist2(data,data).*(glcms~=0);

W=zeros(N);
W_num=zeros(N);

W_den=zeros(N);
OMG1=0.1;
for i=1:N
    for j=1:N
        if(Euc(i,j)~=0)
         W_num(i,j)=exp(-OMG1*(Euc(i,j)));
         W_den(i,i)=W_num(i,j)+W_den(i,i);
        end
    end
end


This can be replaced by W_num = exp(-OMG1*Euc).*(Euc~=0); and W_den = sum(W_num); (here W_den is a vector, not a matrix, but that's what you want because you were only using diagonal elements of your W_den matrix).

for i=1:N
    for j=1:N
         if(Euc(i,j)~=0)
             W(i,j)=(W_num(i,j))/(W_den(i,i));   % Connectivity Matrix W
         end
    end
end


This can be replace by W = W_num./W_den.'; and W(isnan(W)) = 0;.

s_star_temp=zeros(N);   %temporary variable for geodesic distance measurement
W_sparse=zeros(N);
W_sparse=sparse(W);
for i=1:N
    for j=1:N
        if W(i,j)==0 &amp; i~=j;
            s_star_temp(i,j)=graphshortestpath(W_sparse,i,j); % Geodesic Distance
        end
    end
end


Here you define W_sparse twice, you don't need to initialise it with zeros(N). This part is finding all shortest paths, and there is a builtin for that (which is probably a lot more efficient that generating the paths independently). Replace this whole section with s_star_temp = graphallshortestpaths(W).*(W==0); This seems to speed it up from ~15 seconds to 0.005 seconds.



So this is how I would write your algorithm

Img = imread('input.jpg');
[rows, columns, numberOfColorChannels] = size(Img);
[L,N] = superpixels(Img,250);

idx = label2idx(L);

%Identifying neighborhood relationships
glcms = graycomatrix(L,'NumLevels',N,'GrayLimits',[1,N],'Offset',  [0,1;1,0]); %Create gray-level co-occurrence matrix from image
glcms = sum(glcms,3);    % add together the two matrices
glcms = glcms + glcms.'; % add upper and lower triangles together, make it symmetric
glcms(1:N+1:end) = 0;    % set the diagonal to zero, we don't want to see "1 is neighbor of 1"

data = zeros(N,3);
for labelVal = 1:N
    data(labelVal,1) = mean(Img(idx{labelVal}));
    data(labelVal,2) = mean(Img(idx{labelVal}+rows*columns));
    data(labelVal,3) = mean(Img(idx{labelVal}+2*rows*columns));
end

Euc = sparse(pdist2(data,data).*(glcms~=0));

OMG1 = 0.1;

W_num = exp(-OMG1*Euc).*(Euc~=0);
W_den = sum(W_num);

W = W_num./W_den.';
W(isnan(W)) = 0;

s_star_temp = graphallshortestpaths(W).*(W==0);


Notice that I made Euc a sparse matrix. This makes the resulting matrices (W_num, W_den, W) also sparse. Your original code left them as full matrices, which was unnecessary. This code runs in about 1.3 seconds on my machine, compared to about 16 for your original code (with the picture I used).



Minor general comments on your code:


Use spaces around your = signs, it makes things easier to read
Try to avoid very long lines, for example the graycomatrix line with its comment will go over the end of the window and be lost
Write more comments
Matlab best practice is to avoid using i and j for loop variables, prefer ii and jj. It's not very important, but i and j can also mean the imaginary unit.