八领域搜索，水平集分割，分水岭分割，
首先是水平集方法，效果并不好
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)
转载于:https://www.cnblogs.com/natalie/p/5631100.html
相关资源：基于python opencv的目标图像自动识别提取感兴趣区域（本项目提取矩形目标区域）.rar