lichunming.m

Level Set Evolution Without Re-initialization一文的matlab实现

clear;
I=imread('ChanVese1.jpg');
I=I(:,:,1);
I=double(I);


u0=I;
isinside = 0;
center=[size(I)/2];
radius=min(size(I))/2; 
d_it=10;


%初始化参数
ITERATIONS=500;
delta_t=10;

lambda=5;


ve=-1.5;
epsilon=1.5;
mu=0.04;


%初始化符号函数
imsize=size(I);
m=imsize(1); n=imsize(2);
phi=zeros(imsize);

for i = 1 : m;
  for j = 1 : n;
     distance = sqrt( sum( ( center - [ i, j ] ).^2 ) );
     phi( i, j ) = distance - radius;   
     phi( i, j ) = -phi( i, j );
     if phi(i,j)>0.00001
	phi(i,j)=4;
    	 else if   phi(i,j)<-0.00001
		phi(i,j)=-4;
         end
     end
  end
end

%计算边缘检测函数
h=fspecial('gaussian',3,0.5);
I=imfilter(I,h);
[Ix Iy]=gradient(I);
g=1./(1+Ix.^2+Iy.^2);

%进行迭代运算
for ii=1:ITERATIONS
	%显示当前的迭代次数
	fprintf( 1, '%d\n', ii );
		 
	 % 每d_it显示一次图像
 	 if( mod( ii - 1, d_it ) == 0 )
		pause(0.2);clf;
		imagesc(uint8(I));colormap(gray)
        	hold on;
		contour(phi,[0,0],'r');	
         end;
	
	%狄拉克函数
%	dirac=0.5 .* ( 1 + ( 2 ./ pi ) .* atan( phi ./ epsilon ) );		
	dirac=zeros(m,n);
	for i=1:m,
	   for j=1:n,
		if abs(phi(i,j))<=epsilon
	 	    dirac(i,j)=1./(2*epsilon).*(1+cos(pi*phi(i,j)/epsilon));
		end
	   end
	end	

	%计算第一部分
	  
		% 中心差分
  		dx_central = ( circshift( phi, [ 0, -1 ] ) - circshift( phi, [ 0, 1 ] ) ) / 2;
  		dy_central = ( circshift( phi, [ -1, 0 ] ) - circshift( phi, [ 1, 0 ] ) ) / 2;
		
		abs_grad_phi = ( sqrt( dx_central.^2 + dy_central.^2 ) + 0.00001 );
		x=(1-1./abs_grad_phi).* dx_central;
		y=(1-1./abs_grad_phi).* dy_central;
		
		first=mu.*divergence(x,y);
	
	%计算第二部分
		
		x = dx_central ./ abs_grad_phi;
 		y = dy_central ./ abs_grad_phi;
		x=x.*g;
		y=y.*g;
		second=lambda.*dirac.*divergence(x,y);

	% 计算第三部分
		third=ve.* g.*dirac;

	%进行更新
	phi=phi+delta_t.*(first+second+third);
end;