www.pudn.com > MATLAB.zip > demo3.m, change:2006-08-22,size:2512b


clear all; 
close all; 
 
[X,Y] = meshgrid(0:.01:2, 0:.01:2); 
Z = zeros(size(X,1),size(X,2)); 
% I = find((X-1).^2+(Y-1).^2>0.25); 
 
%D2 
% Z(I) = 1 - sqrt(1 - (sqrt((X(I)-1).^2+(Y(I)-1).^2)-0.5).^2); 
 
% Z(I) = 1 - sqrt(1 - (sqrt((X(I)-1).^2+(Y(I)-1).^2)-0.6).^2); 
 
%D1 
I = find((X-1).^2+(Y-1).^2<0.25); 
% Z(I) = sqrt(0.25-(X(I)-1).^2-(Y(I)-1).^2); 
% Z(I) = 2.25 - sqrt(2.5-2*((X(I)-1).^2+(Y(I)-1).^2)); 
 
%D0 
Z(I) = 0.5 + sqrt(0.25-(X(I)-1).^2-(Y(I)-1).^2); 
 
zErrB = 0.0; 
 
for i=1:size(Z,1), 
    for j=1:size(Z,2), 
        Z(i,j) = Z(i,j) + 2*zErrB*(rand(1)-0.5); 
    end 
end 
 
surf(X,Y,Z) 
shading interp; 
axis equal; 
figure; 
 
plot3(X(:),Y(:),Z(:),'.'); 
axis equal; 
figure; 
 
Img = Z; 
 
[Ix,Iy]=gradient(Img); 
imshow(uint8(Img)); 
 
f=Ix.^2+Iy.^2; 
g=1./(1+f);  % edge indicator function. 
 
figure; 
g = mat2gray(g); 
imshow(g); 
 
Img = mat2gray(Z); 
 
lambda=30; % coefficient of the weighted length term L(\phi) 
mu=0.04;  % coefficient of the internal (penalizing) energy term P(\phi) 
alf=3;  % coefficient of the weighted area term A(\phi) 
epsilon=1.5; % the papramater in the definition of smooth Dirac function 
timestep=5;  % time step 
 
% define initial level set function (LSF) as -c, 0, c at points outside, on 
% the boundary, and inside of a region R, respectively. 
[nrow, ncol]=size(Img); 
c=10; 
initialLSF=c*ones(nrow,ncol); 
d=30; 
initialLSF(d+1:end-d, d+1:end-d)=0;  % zero level set on the boundary of R 
initialLSF(d+2:end-d-1, d+2: end-d-1)=-c; % negative constant inside of R 
u=initialLSF; 
figure;imagesc(Img);colormap(gray);hold on; 
[c,h] = contour(u,[0 0],'r','LineWidth',2); 
 
% pause(3); 
 
title('Initial contour'); 
% start level set evolution 
axis equal; axis off; 
for n=1:280 
    u=EVOLUTION_PLSE(u, g ,lambda, mu, alf, epsilon, timestep, 1); 
    pause(0.01); 
    if mod(n,20)==0 
        imagesc(Img);colormap(gray);hold on; 
        axis equal; axis off; 
        [c,h] = contour(u,[0 0],'g','LineWidth',3); 
        iterNum=[num2str(n), ' iterations']; 
        title(iterNum); 
        hold off; 
    end 
end 
 
imagesc(Img);colormap(gray);hold on; 
[c,h] = contour(u,[0 0],'g','LineWidth',4); 
totalIterNum=[num2str(n), ' iterations']; 
title(['Final contour, ', totalIterNum]); 
axis equal; axis off; 
 
%Mean square error 
 
r = size(Img,1)/4; 
center = size(Img)/2; 
error = 0; 
for i=1:size(c,2), 
    error = error + abs((c(1,i)-center(1))^2+(c(2,i)-center(2))^2-r^2)/size(Img,1)^2; 
end 
error = sqrt(error/size(c,2));