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)
axis equal;
figure;

plot3(X(:),Y(:),Z(:),'.');
axis equal;
figure;

Img = Z;

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));