# Histogram separation energy equation based on Bhattacharyya

I am working in level set method, specially Lankton method paper. I try to implement Histogram Separation (HS) Energy problem (Part III.C). It based on Bhattacharyya to control the evolution of contour.

To understand it, the first consider the global method in which we are given an input image and a contour. The contour divides the image into an inside and outside region. The Bhattacharyya distance is calculated by

$$B_{global}=\sqrt{P_{u}*P_{v}}$$

where ${P_{u}}$ and ${P_{v}}$ are pdf of inside and outside regions.

Now we return the Lankton paper which is a local level set. He divides the image into small region by Ball function B(x,y). Then, the contour will separate these regions into inside and outside region. Each small regions have ${P_{u}}$ and ${P_{v}}$. And we can calculate Bhattacharyya distance.

$$F_{HS}=B_{local}=\sqrt{P_{u,x}*P_{v,x}}$$

where $x$ is center point of $B(x,y)$ in image.

I understand that that term. But I cannot understand second term in the formula (integration term, in $\color{red}{\textrm{red}}$ below). Can I represent it by convolution operator? Here is my source code. It works but the result is not correct.

\begin{aligned} \frac{d\phi}{dt}(x) =\ &\delta\phi(x) \int_{\Omega_y}\frac{\mathcal{B}(x,y)\delta\phi(y)}{2} \times \left[F_{HS}\left(\frac{1}{A_v} - \frac{1}{A_u}\right)\right. \\ &+\left.\color{red}{\int_z K(z - I(y)) \times \left(\frac{1}{A_u}\sqrt\frac{P_{v,x}(z)}{P_{u,x}(z)} - \frac{1}{A_v}\sqrt\frac{P_{u,x}(z)}{P_{v,x}(z)}\right) dz} \right] dy \\ &+ \lambda\delta\phi(x)\ \mathrm{div}\left(\frac{\nabla\phi(x)}{\left|\nabla\phi(x)\right|}\right) \end{aligned}

where $K$ is a Gaussian kernel and $A_v$ and $A_u$ are areas inside and outside of these regions.

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(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(img(upts))/ Ain(i) + eps; % leave histograms unsmoothed
q = imhist(img(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(:));%integeration can do by this code?
end

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


## migrated from codereview.stackexchange.comAug 20 '14 at 7:52

This question came from our site for peer programmer code reviews.

• What does this have to do with code or reviewing? – syb0rg Aug 20 '14 at 1:27
• This question appears to be off-topic because it has nothing to do with reviewing code. – Jamal Aug 20 '14 at 1:35
• @syb0rg and Jamal: Sorry, My brower is fail, so I does not upload code successfully – johnmark Aug 20 '14 at 2:14
• Are you just looking for a general review? Also, does this code work as intended? – Jamal Aug 20 '14 at 2:21
• Hi johnmark. Here on Code Review we require that your code be working as intended before we take a look. Your comment I don't know how to code the integral term leads me to believe this is not the case, and so I've voted to leave your question closed. – Schism Aug 20 '14 at 2:51