1
2
3
4
5
6
7
8
9
10
11
| clc clear all; close all; im = imread('coins.png'); figure(1),imshow(im);
dimension=ndims(im); if dimension==3 im = rgb2gray(im); end Img1 = im; Img=double(Img1); [nrow, ncol]=size(Img); BW = roipoly(Img); figure(8),imshow(BW); title('Polygon'); u=8*(0.5-BW); figure(9);imshow(He, [0, 255]);hold on; [c,h] = contour(u,[0 0],'r'); [Ix,Iy]=gradient(Img); f=Ix.^2+Iy.^2; g=1./(1+f); lambda=5; delt=3; iteration=200; for n=1:iteration u=levelset(u, g, lambda, delt); if mod(n,25)==0 pause(1); figure(10),imshow(He, [0, 255]); hold on; [c,h] = contour(u,[0 0],'r'); iterNum=[num2str(n), ' iterations']; title(iterNum); hold off; end end
function u = levelset(u0, g, lambda, delt) u=u0; [vx,vy]=gradient(g); epsilon=1.5; mu=0.2/5; alf=5;
for k=1:1 u=weightvalue(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 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 = weightvalue(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]); |