speeding up code (mean shift colour image segmentation)

I know there are several links out there - but i can't seem to speed up my code (below) -

The code implements *mean-shift algorithm* for colour image segmentation I tried reducing the number of steps, and even calculated *first_exp* outside the while loop

to run 10 iterations with a large *window size* (for e.g 61x61) it's taking about 10 hours.

clc;

clear all;

        tic;

        seg1=read_raw('seg1.raw',321,481,3);

        figure(1)
        imshow(uint8(seg1));

        seg1extended(:,:,1)=extend(seg1(:,:,1),24);
        seg1extended(:,:,2)=extend(seg1(:,:,2),24);
        seg1extended(:,:,3)=extend(seg1(:,:,3),24);

        figure(2)
        imshow(uint8(seg1extended));

        in= seg1extended;    
        out=zeros(size(in));

        sigma_s=256;
        sigma_r=1024;

        iter=1;
        MAXITER=4;
        filter_window_dimension=24;

        ES=filter_window_dimension/2;

        R=size(in,1);
        C=size(in,2);

        while (iter<=MAXITER)
        iter

          for l=1+ES:R-ES
          for m=1+ES:C-ES
        l
        m
                             sum_w=0;
                            sum_nr=0;

                    for win_w=-ES:ES
                    for win_h=-ES:ES

                        yi=in(l,m,:);

                        dist= sqrt((l-(l+win_w))^2 + (m-(m+win_h))^2);
                        first_exp= -((dist^2)/sigma_s); %spatial exponent

                        feat= yi-in(l+win_w,m+win_h,:);
                        feat= sqrt((feat(1)^2) + (feat(2)^2) + (feat(3)^2));

                        second_exp= -((feat^2)/sigma_r); %range exponent

                        w_term=exp(first_exp)*exp(second_exp)*feat;

                        sum_w= sum_w + w_term; %denominator of yi_k
                        sum_nr = sum_nr + (in(l+win_w,m+win_h,:).*w_term); %numerator of yi_k

                        yi=ceil(sum_nr./sum_w);

                    end
                    end
          out(l,m,:)=yi;          
          end
          end

         iter=iter+1;  
        end
        outp=out(1+ES:R-ES,1+ES:C-ES,:);

        figure(4)
        imshow(uint8(outp));

        fid=fopen('seg1result_1632.raw','wb');
        fwrite(fid,outp);

        time_elapsed=toc

*modifications - made *

clc;
    clear all;
    tic;

      seg1=read_raw('seg1.raw',321,481,3);

      figure(1)
      imshow(uint8(seg1));

      seg1extended(:,:,1)=extend(seg1(:,:,1),66);
      seg1extended(:,:,2)=extend(seg1(:,:,2),66);
      seg1extended(:,:,3)=extend(seg1(:,:,3),66);

      figure(2)
      imshow(uint8(seg1extended));

      image= seg1extended;

      in=image;

      out=zeros(size(in));

      sigma_s=256;
      sigma_r=1024;

      iter=1;
      MAXITER=3;

      filter_window_dimension=66;
      ES=filter_window_dimension/2;

      R=size(in,1);
      C=size(in,2);

      for win_w = -ES : ES
      for win_h = - ES : ES
          dist(win_w+ES+1,win_h+ES+1)=sqrt((win_h^2) + (win_w^2));

          first_exp(win_w+ES+1,win_h+ES+1)=-(power(dist(win_w+ES+1,win_h+ES+1),2))/(sigma_s);
      end
      end

      while (iter<=MAXITER)
          iter

          for l=1+ES:R-ES
          for m=1+ES:C-ES
              l
              m

              sum_w=0;
              sum_nr=0;

              for win_w=-ES:ES
                  for win_h=-ES:ES

                      yi=in(l,m,:);

                      feat=sqrt((yi(1)-in(l+win_w,m+win_h,1)).^2  + (yi(2)-in(l+win_w,m+win_h,2)).^2 + (yi(3)-in(l+win_w,m+win_h,3)).^2);
                      second_exp=-((feat^2)/sigma_r);

                      w_term= exp(first_exp(1+win_w+ES,1+win_h+ES))*exp(second_exp)*feat;

                      sum_w= sum_w+w_term;
                      sum_nr=sum_nr + (in(l+win_w,m+win_h,:)*w_term);

                      yi=ceil(sum_nr/sum_w);

                  end
              end
              out(l,m,:)=yi;
          end
          end
          iter=iter+1;
      end

              outp=out(1+ES:R-ES,1+ES:C-ES,:);

              figure(4)
              imshow(uint8(outp));

              toc;

the code is written keeping in mind that we cannot use any inbuilt functions. this is a homework project.

Answer 1

[Image Analyst]

I'm not about to try something for you that takes 10 hours. I suggest you use the timing functions to see where most of the time is taken up.

[pooja]

Dear @image analyst

The problem isn't finding out where it's taking long. I know it's taking long within the nested for loop.

Where I'm computing "second_exp" and until I output to a new image matrix . I'm seeking alternative approaches to frame the algorithm (first and second exp are to get gaussian kernel weights)

The modified code runs in 15 mins for one iteration whereas the original one takes approximately an hour for one iteration.

If you have any hints or you suspect some other modification might work, please let me know.

Basically , the 15 mins code should be reduced to around minutes .

[Image Analyst]

Try switching the order of the for loops so that you're going down rows first before moving over columns. This is the way the array is stored in memory. Your way of going over columns before incrementing row is the slower way.