Local histogram

I am working in a level method, especially the Lankton document method . I am trying to implement the (HS) Energy Histogram problem (Part III.C). It is based on Bhattacharya to control the evolution of the circuit. To understand this, let's first consider a global method that defines an input image and an outline. The outline divides the image into inner and outer regions. Bhattacharya distance calculated for

B=sqrt (P_in.*P_out)

      

        
        
        
      

    

where P_in and P_pout are pdfs inside and outside regions. To apply Bhattacharyya globally, you can see the source code here . Now we are returning Lentenko paper. This is the local level. In which he divides the image into a small area using the Ball function. The outline will then divide these regions into an inner and outer region. Each small region has a P_in and a P_out. And we can calculate the distance of Bhattacarya. I took this step. But I cannot implement the final step as formal. Could you help me??? and Av and Au are the area inside and outside these regions. This is my main code. You can download the source code

  for its = 1:max_its   % Note: no automatic convergence test

%-- get the curve 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;