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);

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