I want to know the code for weighted histogram seperation

weighted histogram seperation matlab 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;