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

Question

我将基于简单的线性反馈控制系统 (LFCS) 实现显着目标检测方法。控制系统模型表示为如下方程​​:

我想出了以下程序代码，但结果不是应该的。具体来说，输出应该类似于下图:

但是代码产生了这样的输出:

代码如下。

%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

以上代码是后续部分的先决条件(实际上，它们为下一节生成必要的数据和矩阵。提供上述代码是为了使整个过程可重现)。

具体来说，我认为这部分没有给出预期的结果。恐怕我没有正确模拟使用 for 循环的并行性。此外，永远不会满足终止条件(与 for 和 if 语句一起使用以模拟 do-while 循环)，并且循环会一直持续到最后一次迭代(而不是在指定条件出现时终止)。这里的一个主要问题是终止条件是否得到正确实现。以下代码的伪算法如下图:

 %%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

Answer 1

您的部分终止标准是:

max(abs(Y_a(:,t)-Y_a(:,t-1)))<=eps

说 Y_a 趋向于 2。你真的很接近......事实上，如果后续值不相同，你可以获得的最接近的值是 Y_a(t)-Y_a(t-1 ) == 4.4409e-16。如果这两个值更接近，则它们的差将为 0，因为这是表示浮点值的精度。所以你已经达到了离你的目标如此接近的程度。随后的迭代以尽可能小的量更改目标值，4.4409e-16。但是你的测试返回错误!为什么？因为eps == 2.2204e-16!

eps 是 eps(1) 的简写，是 1 和下一个可表示的较大值之间的差值。由于浮点值的表示方式，此差值是 2 与下一个可表示的较大值(由 eps(2) 给出)之间差值的一半。

但是，如果 Y_a 趋向于 1e-16，后续迭代可能会使 Y_a 的值加倍或减半，您仍然会遇到停止标准!

因此，您需要提出一个合理的停止标准，该标准是目标值的一小部分，如下所示:

max(abs(Y_a(:,t)-Y_a(:,t-1))) <= 1e6 * eps(max(abs(Y_a(:,t))))

---

不请自来的建议

您真的应该研究一下 MATLAB 中的矢量化运算。例如，

for y=1:N
   Temp_sig_a=0;
   for z=1:N
      Temp_sig_a = Temp_sig_a + abs(Y_a(y,t)- Y_a(z,t));
   end
   sigma_a(t-1,y)= Temp_sig_a;
end

可以写成

for y=1:N
   sigma_a(t-1,y) = sum(abs(Y_a(y,t) - Y_a(:,t)));
end

这又可以写成

sigma_a(t-1,:) = sum(abs(Y_a(:,t).' - Y_a(:,t)));

避免循环通常不仅效率更高，而且代码更短更易于阅读。

另外，这个:

P_FB_a = diag(sigma_a(t-1,:));
u_a(:,t) = u_a(:,t-1) + P_FB_a * Y_a(:,t);

与

相同

u_a(:,t) = u_a(:,t-1) + sigma_a(t-1,:).' .* Y_a(:,t);

当然，创建一个对角矩阵并用如此多的零进行矩阵乘法比直接计算逐元素乘法要昂贵得多。