我正在使用级别设置方法,特别是Lankton方法paper。我尝试解决直方图分离(HS)能耗问题(第III.C部分)。它基于Bhattacharyya来控制轮廓的演变。为了理解它,首先我们考虑全局方法,其中给定输入图像和轮廓。轮廓线将图像分为内部和外部区域。 Bhattacharyya距离的计算公式为

B=sqrt (P_in.*P_out)

其中P_in和P_pout是内部和外部区域的pdf。
要将Bhattacharyya应用于全局级别集,可以在here上查看源代码。现在我们返回Lankton论文。它是本地级别设置的。其中,他通过Ball函数将图像划分为较小的区域。然后,轮廓会将这些区域分为内部区域和外部区域。每个小区域都有P_in和P_out。我们可以计算出Bhattacharyya距离。我完成了这一步。但是我无法正式执行最后一步。你能帮助我吗???

Av和Au是这些区域内部和外部的区域。这是我的主要代码。您可以在source code下载
  for its = 1:max_its   % Note: no automatic convergence test

%-- get the curve's narrow band
idx = find(phi <= 1.2 & phi >= -1.2)';
[y x] = ind2sub(size(phi),idx);

%-- get windows for localized statistics
xneg = x-rad; xpos = x+rad;      %get subscripts for local regions
yneg = y-rad; ypos = y+rad;
xneg(xneg<1)=1; yneg(yneg<1)=1;  %check bounds
xpos(xpos>dimx)=dimx; ypos(ypos>dimy)=dimy;

%-- re-initialize u,v,Ain,Aout
Ain=zeros(size(idx)); Aout=zeros(size(idx));
B=zeros(size(idx));integral=zeros(size(idx));
%-- compute local stats
for i = 1:numel(idx)  % for every point in the narrow band
  img = I(yneg(i):ypos(i),xneg(i):xpos(i)); %sub image
  P = phi(yneg(i):ypos(i),xneg(i):xpos(i)); %sub phi

  upts = find(P<=0);            %local interior
  Ain(i) = length(upts)+eps;

  vpts = find(P>0);             %local exterior
  Aout(i) = length(vpts)+eps;

  %% Bha distance
  p = imhist(I(upts))/ Ain(i) + eps; % leave histograms unsmoothed
  q = imhist(I(vpts)) / Aout(i) + eps;
  B(i) = sum(sqrt(p.* q));
  term2= sqrt(p./q)/Aout(i) - sqrt(q./p)/Ain(i); %Problem in here===I don't know how to code the integral term
  integral(i) =sum(term2(:));
end

F =-B./2.*(1./Ain - 1./Aout) - integral./2;

最佳答案

我试过了-不知道它是否正确-它没有直方图平滑化(我认为没有必要)

if type==3  % Set up for bhatt

F=zeros(size(idx,1),2);

for i = 1:numel(idx)

img2 = img(yneg(i):ypos(i),xneg(i):xpos(i));
P = phi(yneg(i):ypos(i),xneg(i):xpos(i));

upts = find(P<=0);            %local interior
Ain = length(upts)+eps;
[u,~] = hist(img2(upts),1:256);

vpts = find(P>0);             %local exterior
Aout = length(vpts)+eps;
 [v,~] = hist(img2(vpts),1:256);

   Ap = Ain;
   Aq = Aout;
   In=Ap;
   Out=Aq;
    try
     p = ((u))  ./ Ap + eps;
     q = ((v)) ./ Aq + eps;
     catch
     g
    end


 B = sum(sqrt(p .* q));

 F(i)=B.*((1/numel(In))-(1/numel(Out)))+(0.5.*(1/numel(In)*(q(img(idx(i))+1)...      /p(img(idx(i))+1))))-(0.5.*(1/numel(Out)*(p(img(idx(i))+1)/q(img(idx(i))+1))));
end

09-07 06:55