关于感兴趣区域提取

八领域搜索，水平集分割，分水岭分割，

首先是水平集方法，效果并不好

clear all;
close all;
Img = imread(‘index_1.bmp‘);  % The same cell image in the paper is used here
Img=double(Img(:,:,1));
sigma=1.5;    % scale parameter in Gaussian kernel for smoothing.
G=fspecial(‘gaussian‘,15,sigma);
Img_smooth=conv2(Img,G,‘same‘);  % smooth image by Gaussiin convolution
[Ix,Iy]=gradient(Img_smooth);
f=Ix.^2+Iy.^2;
g=1./(1+f);  % edge indicator function.

epsilon=1.5; % the papramater in the definition of smoothed Dirac function

timestep=5;  % time step
mu=0.2/timestep;  % coefficient of the internal (penalizing) energy term P(\phi)
          % Note: the product timestep*mu must be less than 0.25 for stability!

lambda=5; % coefficient of the weighted length term Lg(\phi)
alf=1.5;  % coefficient of the weighted area term Ag(\phi);
          % Note: Choose a positive(negative) alf if the initial contour is outside(inside) the object.

% define initial level set function (LSF) as -c0, 0, c0 at points outside, on
% the boundary, and inside of a region R, respectively.
[nrow, ncol]=size(Img);
c0=4;
initialLSF=c0*ones(nrow,ncol);
w=8;
initialLSF(w+1:end-w, w+1:end-w)=0;  % zero level set is on the boundary of R.
                                     % Note: this can be commented out. The intial LSF does NOT necessarily need a zero level set.

initialLSF(w+2:end-w-1, w+2: end-w-1)=-c0; % negative constant -c0 inside of R, postive constant c0 outside of R.
u=initialLSF;
figure;imagesc(Img);colormap(gray);hold on;
[c,h] = contour(u,[0 0],‘r‘);
title(‘Initial contour‘);

% start level set evolution
for n=1:300
    u=EVOLUTION(u, g ,lambda, mu, alf, epsilon, timestep, 1);
    if mod(n,20)==0
        pause(0.001);
        imagesc(Img);colormap(gray);hold on;
        [c,h] = contour(u,[0 0],‘r‘);
        iterNum=[num2str(n), ‘ iterations‘];
        title(iterNum);
        hold off;
    end
end
imagesc(Img);colormap(gray);hold on;
[c,h] = contour(u,[0 0],‘r‘);
totalIterNum=[num2str(n), ‘ iterations‘];
title([‘Final contour, ‘, totalIterNum]);

function u = EVOLUTION(u0, g, lambda, mu, alf, epsilon, delt, numIter)
%  EVOLUTION(u0, g, lambda, mu, alf, epsilon, delt, numIter) updates the level set function
%  according to the level set evolution equation in Chunming Li et al‘s paper:
%      "Level Set Evolution Without Reinitialization: A New Variational Formulation"
%       in Proceedings CVPR‘2005,
%  Usage:
%   u0: level set function to be updated
%   g: edge indicator function
%   lambda: coefficient of the weighted length term L(\phi)
%   mu: coefficient of the internal (penalizing) energy term P(\phi)
%   alf: coefficient of the weighted area term A(\phi), choose smaller alf
%   epsilon: the papramater in the definition of smooth Dirac function, default value 1.5
%   delt: time step of iteration, see the paper for the selection of time step and mu
%   numIter: number of iterations.
%

u=u0;
[vx,vy]=gradient(g);

for k=1:numIter
    u=NeumannBoundCond(u);
    [ux,uy]=gradient(u);
    normDu=sqrt(ux.^2 + uy.^2 + 1e-10);
    Nx=ux./normDu;
    Ny=uy./normDu;
    diracU=Dirac(u,epsilon);
    K=curvature_central(Nx,Ny);
    weightedLengthTerm=lambda*diracU.*(vx.*Nx + vy.*Ny + g.*K);
    penalizingTerm=mu*(4*del2(u)-K);
    weightedAreaTerm=alf.*diracU.*g;
    u=u+delt*(weightedLengthTerm + weightedAreaTerm + penalizingTerm);  % update the level set function
end

% the following functions are called by the main function EVOLUTION
function f = Dirac(x, sigma)   %水平集狄拉克计算
f=(1/2/sigma)*(1+cos(pi*x/sigma));
b = (x<=sigma) & (x>=-sigma);
f = f.*b;

function K = curvature_central(nx,ny);  %曲率中心
[nxx,junk]=gradient(nx);
[junk,nyy]=gradient(ny);
K=nxx+nyy;

function g = NeumannBoundCond(f)
% Make a function satisfy Neumann boundary condition
[nrow,ncol] = size(f);
g = f;
g([1 nrow],[1 ncol]) = g([3 nrow-2],[3 ncol-2]);
g([1 nrow],2:end-1) = g([3 nrow-2],2:end-1);
g(2:end-1,[1 ncol]) = g(2:end-1,[3 ncol-2]);

　　%%下面是八领域搜索算法代码，没搞明白怎么回事，先贴上来吧。。、、

clear all;
close all;
clc;
%外边界
img=imread(‘rice.png‘);
img=img>128;
imshow(img);
[m n]=size(img);

imgn=zeros(m,n);        %边界标记图像
ed=[-1 -1;0 -1;1 -1;1 0;1 1;0 1;-1 1;-1 0]; %从左上角像素判断
for i=2:m-1
    for j=2:n-1
        if img(i,j)==1      %如果当前像素是前景像素

            for k=1:8
                ii=i+ed(k,1);
                jj=j+ed(k,2);
                if img(ii,jj)==0    %当前像素周围如果是背景，边界标记图像相应像素标记
                    imgn(ii,jj)=1;
                end
            end

        end
    end
end

figure;
imshow(imgn,[]);

%不过要是真取二值图像外边界，通常是原图膨胀图减去原图就行了
se = strel(‘square‘,3);
imgn=imdilate(img,se)-img;
figure;
imshow(imgn)