问matlab中的并行循环和图像处理

EN

%Calculation of euclidian distance between adjacent superpixels stores in variable of Euc

  A = imread('aa.jpg'); 
  [rows, columns, cnumberOfColorChannels] = size(A);
  [L,N] = superpixels(A,400);

  %% Determination of adjacent superpixels
  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"

  idx = label2idx(L);    % Convert label matrix to cell array of linear indices
  numRows = size(A,1);
  numCols = size(A,2);

 %%Mean color in Lab color space for each channel

 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(A(redIdx));
data(labelVal,2) = mean(A(greenIdx));
data(labelVal,3) = mean(A(blueIdx));

end    

Euc=zeros(N);

  %%Calculation of euclidian distance between adjacent superpixels stores in Euc

for a=1:N
for b=1:N
    if glcms(a,b)~=0
        Euc(a,b)=sqrt(((data(a,1)-data(b,1))^2)+((data(a,2)-data(b,2))^2)+((data(a,3)-data(b,3))^2));
    end
end
end


 %%Creation of Connectivity matrix "W" between adjacent superpixels

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

 W_den=zeros(N);
 OMG1=0.1;
 for c=1:N
 for d=1:N
    if(Euc(c,d)~=0)
     W_num(c,d)=exp(-OMG1*(Euc(c,d)));

      W_den(c,c)=W_num(c,d)+W_den(c,c);  % 

    end
end
end

%Connectivity matrix W between adjacent superpixels 

for e=1:N
for f=1:N
     if(Euc(e,f)~=0)
         W(e,f)=(W_num(e,f))/(W_den(e,e));

     end
end
end


   %%calculation of geodesic distance between nonadjacent superpixels  stores in variable "s_star_temp"

  s_star_temp=zeros(N);   %temporary variable for geodesic distance measurement
  W_sparse=zeros(N);
  W_sparse=sparse(W);
  for g=1:N
  for h=1:N
    if W(g,h)==0 & g~=h;
        s_star_temp(g,h)=graphshortestpath(W_sparse,g,h,'directed',false); 
    end
end
end


  %%Calculation of connectivity matrix for nonadjacent superpixels stores in "S_star" variable" 

  S_star=zeros(N);
  OMG2=8;   
  for i=1:N
  for j=1:N
    if s_star_temp(i,j)~=0
        S_star(i,j)=exp(-OMG2*s_star_temp(i,j));
    end
end
end


  %%Calculation of connectivity matrix "S" for measuring connectivity between all superpixels

 S=zeros(N);

 S=S_star+W;


 %% Defining non-isolation level for connectivity matrix "W" 
 g_star=zeros(N);

 for k=1:N
 g_star(k,k)=max(W(k,:));   
 end


   %%Limiting the range of g_star and calculation of isolation cue matrix "G"

  alpha1=0.15;
  alpha2=0.85;
  G=zeros(N);
  for l=1:N
  G(l,l)=alpha1*(g_star(l,l)- min(g_star(:)))/(max(g_star(:))- min(g_star(:)))+(alpha2 - alpha1);
  end



  %%Determining the supperpixels that surrounding the image boundary
  lr = L([1,end],:);

  tb = L(:,[1,end]);

  labels = unique([lr(:);tb(:)]);



  %% Calculation of background likelihood for each superpixels stores in"BgLike"
 sum_temp=0;
 temp=zeros(1,N);
 BgLike=zeros(N,1);
 BgLike_num=zeros(N);
 BgLike_den=zeros(N);

for m=1:N
for n=1:N
    if ismember(n,labels)==1

        BgLike_num(m,m)=S(m,n)+ BgLike_num(m,m);

    end
   end
  end

 for o=1:N
 for p=1:N
    for q=1:N
        if W(p,q)~=0
            temp(q)=S(o,p)-S(o,q);
        end
    end
          sum_temp=max(temp)+sum_temp;
          temp=0;
end
BgLike_den(o,o)=sum_temp;
sum_temp=0;
end


for r=1:N

    BgLike(r,1)= BgLike_num(r,r)/BgLike_den(r,r); 

end


  %%%%Calculation of Foreground likelihood for each superpixels stores in "FgLike"

 FgLike=zeros(N,1);

 for s=1:N
 for t=1:N
    FgLike(s,1)=(exp(-BgLike(t,1))) * Euc(s,t)+ FgLike(s,1); 
 end
 end

 %%parallel operations for background and foreground  implemented  here
 T0 = 0 ;
 Tf = 20 ;
 Ts = 0.1 ;
 Ti = T0:Ts:Tf ;
 Nt=numel(Ti);
 Y_Bg=zeros(N,Nt);
 Y_Fg=zeros(N,Nt);

 P_Back_Bg=zeros(N,N);
 P_Back_Fg=zeros(N,N);
 u_Bg=zeros(N,Nt);
 u_Fg=zeros(N,Nt);
 u_Bg_Star=zeros(N,Nt);
 u_Fg_Star=zeros(N,Nt);
 u_Bg_Normalized=zeros(N,Nt);
 u_Fg_Normalized=zeros(N,Nt);
 tau=0.1;
 sigma_Bg=zeros(Nt,N);

Temp_Bg=0;
Temp_Fg=0;

C_Bg=zeros(Nt,N);
C_Fg=zeros(Nt,N);

 %%System Initialization

for u=1:N
u_Bg(u,1)=(BgLike(u,1)- min(BgLike(:)))/(max(BgLike(:))- min(BgLike(:)));
u_Fg(u,1)=(FgLike(u,1)- min(FgLike(:)))/(max(FgLike(:))- min(FgLike(:)));
end

%% P_state and P_input
P_state=G*W;           
P_input=eye(N)-G;

% State Initialization


X_Bg=zeros(N,Nt);
X_Fg=zeros(N,Nt);


   for v=1:20   % v starts from 1 because we have no matrices with 0th column number
           %The first column of X_Bg and X_Fg is 0 for system initialization     
       X_Bg(:,v+1)=P_state*X_Bg(:,v) + P_input*u_Bg(:,v);
       X_Fg(:,v+1)=P_state*X_Fg(:,v) + P_input*u_Fg(:,v);
  v=v+1;  
  if v==2
  C_Bg(1,:)=1;       
 C_Fg(1,:)=1;   
 else
       for w=1:N

           for x=1:N

      Temp_Fg=S(w,x)*X_Fg(x,v-1)+Temp_Fg;
      Temp_Bg=S(w,x)*X_Bg(x,v-1)+Temp_Bg;
           end
       C_Fg(v-1,w)=inv(X_Fg(w,v-1)+((Temp_Bg)/(Temp_Fg)*(1-X_Fg(w,v-1))));    
       C_Bg(v-1,w)=inv(X_Bg(w,v-1)+((Temp_Fg)/(Temp_Bg))*(1-X_Bg(w,v-1)));    
       Temp_Bg=0;
       Temp_Fg=0;
       end
 end
 P_Bg=diag(C_Bg(v-1,:));  
 P_Fg=diag(C_Fg(v-1,:));  
 Y_Bg(:,v)= P_Bg*X_Bg(:,v);
 Y_Fg(:,v)= P_Fg*X_Fg(:,v);

 for y=1:N
 Temp_sig_Bg=0;
 Temp_sig_Fg=0;
 for z=1:N
  Temp_sig_Bg = Temp_sig_Bg +S(y,z)*abs(Y_Bg(y,v)- Y_Bg(z,v));
  Temp_sig_Fg = Temp_sig_Fg +S(y,z)*abs(Y_Fg(y,v)- Y_Fg(z,v));
 end
 if Y_Bg(y,v)>= Y_Bg(y,v-1)
    sign_Bg=1;
 else
   sign_Bg=-1;
 end

 if Y_Fg(y,v)>= Y_Fg(y,v-1)
   sign_Fg=1;
 else
   sign_Fg=-1;
 end
 sigma_Bg(v-1,y)=sign_Bg*Temp_sig_Bg;
 sigma_Fg(v-1,y)=sign_Fg*Temp_sig_Fg;
 end

  %Calculation of P_Back for background and foreground
  P_Back_Bg=tau*diag(sigma_Bg(v-1,:));  
  P_Back_Fg=tau*diag(sigma_Fg(v-1,:));

 u_Bg_Star(:,v)=u_Bg(:,v-1)+P_Back_Bg*Y_Bg(:,v);
 u_Fg_Star(:,v)=u_Fg(:,v-1)+P_Back_Fg*Y_Fg(:,v);
 for aa=1:N   %Normalization of u_Bg and u_Fg

 u_Bg(aa,v)=(u_Bg_Star(aa,v)- min(u_Bg_Star(:,v)))/(max(u_Bg_Star(:,v))-min(u_Bg_Star(:,v)));
  u_Fg(aa,v)=(u_Fg_Star(aa,v)- min(u_Fg_Star(:,v)))/(max(u_Fg_Star(:,v))-min(u_Fg_Star(:,v)));

end

if (max(abs(Y_Fg(:,v)-Y_Fg(:,v-1)))<=0.0118) &&(max(abs(Y_Bg(:,v)-Y_Bg(:,v-1)))<=0.0118)  %% epsilon= 0.0118
 break;
 end 
 end

K=4;
T=0.4;
phi_1=(2-(1-T)^(K-1))/((1-T)^(K-2));
phi_2=(1-T)^(K-1);
phi_3=1-phi_1;

for bb=1:N
Y_Output_Preliminary(bb,1)=Y_Fg(bb,v)/((Y_Fg(bb,v)+Y_Bg(bb,v)));
end

for hh=1:N
 Y_Output(hh,1)=(phi_1*(T^K))/(phi_2*(1-Y_Output_Preliminary(hh,1))^K+(T^K))+phi_3;
 end


   V_rs=zeros(N);
   V_Final=zeros(rows,columns);
   for cc=1:rows
   for dd=1:columns
    V_rs(cc,dd)=Y_Output(L(cc,dd),1); 
   end
  end

  maxDist = 10;      % Maximum chessboard distance from image

  wSF=zeros(rows,columns);
  wSB=zeros(rows,columns);

  % Get the range of x and y indices who's chessboard distance from pixel (0,0) are less than 'maxDist'
  xRange = (-(maxDist-1)):(maxDist-1);
  yRange = (-(maxDist-1)):(maxDist-1);

  % Create a mesgrid to get the pairs of (x,y) of the pixels
  [pointsX, pointsY] = meshgrid(xRange, yRange);
  pointsX = pointsX(:);
  pointsY = pointsY(:);

  % Remove pixel (0,0)
  pixIndToRemove = (pointsX == 0 & pointsY == 0);
  pointsX(pixIndToRemove) = [];
  pointsY(pixIndToRemove) = [];

  for ee=1:rows
  for ff=1:columns
    % Get a shifted copy of 'pointsX' and 'pointsY' that is centered
    % around (x, y)
    pointsX1 = pointsX + ee;
    pointsY1 = pointsY + ff;

    % Remove the the pixels that are out of the image bounds        
    inBounds =...
        pointsX1 >= 1 & pointsX1 <= rows &...
        pointsY1 >= 1 & pointsY1 <= columns;

    pointsX1 = pointsX1(inBounds);
    pointsY1 = pointsY1(inBounds);

    % Do stuff with 'pointsX1' and 'pointsY1'

    wSF_temp=0;
    wSB_temp=0;

    for gg=1:size(pointsX1)


        Temp=exp(-OMG1*(sqrt(double(A(pointsX1(gg),pointsY1(gg),1))-double(A(ee,ff,1)))^2+(double(A(pointsX1(gg),pointsY1(gg),2))-double(A(ee,ff,2)))^2 + (double(A(pointsX1(gg),pointsY1(gg),3))-double(A(ee,ff,3)))^2));
        wSF_temp=wSF_temp+(Temp*V_rs(pointsX1(gg),pointsY1(gg)));
        wSB_temp=wSB_temp+(Temp*(1-V_rs(pointsX1(gg),pointsY1(gg))));


    end
    wSF(ee,ff)= wSF_temp;   
    wSB(ee,ff)= wSB_temp;   
    V_Final(ee,ff)=V_rs(ee,ff)/(V_rs(ee,ff)+(wSB(ee,ff)/wSF(ee,ff))*(1-V_rs(ee,ff))); 

end
end

imshow(V_Final,[]);    %%Saliency map of the image