局部直方图分离能量实现

Question

我正在研究水平集方法，特别是兰克顿方法paper。我尝试实现直方图分离 (HS) 能量问题（第 III.C 部分）。它基于 Bhattacharyya 来控制轮廓的演变。为了理解它，首先我们考虑给定输入图像和轮廓的全局方法。轮廓将图像分为内部和外部区域。 Bhattacharyya 距离的计算公式为

B=sqrt (P_in.*P_out)

其中 P_in 和 P_pout 是内部和外部区域的 pdf。 要将 Bhattacharyya 应用于全局级别集，您可以在 here 处查看源代码。现在我们返回兰克顿论文。它是本地级别的设置。其中，他通过 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;

Answer 1

我试过这个 - 不知道它是否正确 - 它没有直方图平滑（我认为没有必要）

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