www.pudn.com > level_set_MCM_0_1_image.rar > level_set_MCM_0_1_image.m
clear;
clc;
%%%%%%%%%%%%%% 图像大小变为二分之一 %%%%%%%%%%%%%%%%%
p=1;
K=20;
%%%%%%%%%%%%%% input image %%%%%%%%%%%%%
I0=imread('5jx.bmp');
% Plot the Curve and convert the plot into Image
% h0=subplot(1,1,1);
% imshow(I0);
% axis off;
% f0 = getframe(h0);
% [im0, map] = frame2im(f0);
if( isrgb( I0 ) )
I0 = rgb2gray( I0 );
end;
I0=double(I0);
[nx0,ny0]=size(I0);
figure(1);imshow(I0);
%%%%%%%% 计算函数g及其一阶偏导 %%%%%%%%%%%%%%%
I0_x = (I0(:,[2:ny0 ny0])-I0(:,[1 1:ny0-1]))/2;
I0_y = (I0([2:nx0 nx0],:)-I0([1 1:nx0-1],:))/2;
% for i=1:nx0
% for j=1:ny0
% g0(i,j)=sqrt(I0_x(i,j)^2+I0_y(i,j)^2);
% end
% end
g0=zeros(nx0,ny0);
g0=(I0_x.^2+I0_y.^2).^(1/2);
% for i=1:nx0
% for j=1:ny0
% g(i,j)=1/(1+(g0(i,j)/K)^p);
% end
% end
g=zeros(nx0,ny0);
g=1./(1+(g0./K).^p);
g_x = (g(:,[2:ny0 ny0])-g(:,[1 1:ny0-1]))/2;
g_y = (g([2:nx0 nx0],:)-g([1 1:nx0-1],:))/2;
%%%%%%%%%%%%%%%%% input yuan_hua_image %%%%%%%%%
I=imread('yuan.bmp');
% Plot the Curve and convert the plot into Image
figure(2);
h=subplot(1,1,1);
imshow(I);
axis off;
f = getframe(h);
[im, map] = frame2im(f);
if( isrgb( im ) )
im = rgb2gray( im );
end;
im=double(im);
[nx,ny]=size(im);
figure(2);imshow(im);
%Detect the curve in the half_sized image
% and store in curvIndex
nnx=floor(nx/2);nny=floor(ny/2);
MaxlLengh=2*nnx+2*nny;
curvImag=255*ones(nnx,nny);
curvIndex=zeros(MaxlLengh,2);
num=0;
for i=1:nx
for j=1:ny
if im(i,j)<5 & (im(i-1,j)>120 | im(i+1,j)>120 | im(i,j-1)>120 |im(i,j+1)>120)
num=num+1;ii=floor(i/2);jj=floor(j/2);
curvImag(ii,jj)=0;
curvIndex(num,1)=ii;
curvIndex(num,2)=jj;
for k=1:num-1
if curvIndex(k,1)==ii & curvIndex(k,2)==jj
num=num-1;break
end
end
end
end
end
figure(3);imshow(uint8(curvImag));
%% Initialize U as distence function to the initial curve
U = zeros(nnx,nny);
dist=zeros(1,num);
for j=1:nny
for i=1:nnx
for k=1:num
dist(k)=sqrt((i-curvIndex(k,1)).^2+(j-curvIndex(k,2)).^2);
end
U(i,j)=min(dist);
if im(2*i,2*j)<5 % if (i,j) is inside of the curve,then negtive
U(i,j)=-U(i,j);
end
end
end
figure(4);surf(U); % Display U(i,j)
colormap(jet); %着色hot or cool or pink etc.
shading flat;view([-18,28]); %flat or interp or faceted
%%%%%%%%%%%%%%%%% Iteration began here
% figure(5);imshow(im0);
dt=0.2;
for n=1:2000
%%%%%% 计算第一项有关数据 %%%%%%%%%
% U_x = (U(:,[2:nny nny])-U(:,[1 1:nny-1]))/2;
% U_y = (U([2:nnx nnx],:)-U([1 1:nnx-1],:))/2;
% U_xx = U(:,[2:nny nny])+U(:,[1 1:nny-1])-2*U;
% U_yy = U([2:nnx nnx],:)+U([1 1:nnx-1],:)-2*U;
U_x = (U([2:nny nny],:)-U([1 1:nny-1],:))/2;
U_y = (U(:,[2:nnx nnx])-U(:,[1 1:nnx-1]))/2;
U_xx = U([2:nny nny],:)+U([1 1:nny-1],:)-2*U;
U_yy = U(:,[2:nnx nnx])+U(:,[1 1:nnx-1])-2*U;
Dp = U([2:nnx nnx],[2:nny nny])+U([1 1:nnx-1],[1 1:nny-1]);
Dm = U([1 1:nnx-1],[2:nny nny])+U([2:nnx nnx],[1 1:nny-1]);
U_xy = (Dp-Dm)/4;
Num = U_xx.*U_y.^2-2*U_x.*U_y.*U_xy+U_yy.*U_x.^2;
Den = U_x.^2+U_y.^2;
I_t = Num./(Den+eps);
%%%%%% 计算第二项有关数据 %%%%%%%%%
D_x_L=U(:,:)-U([1 1:nny-1],:);
D_x_R=U([2:nny nny],:)-U(:,:);
D_y_L=U(:,:)-U(:,[1 1:nny-1]);
D_y_R=U(:,[2:nny nny])-U(:,:);
%%%%%% 计算迭代数据 %%%%%%%%%
% for i=1:nnx
% for j=1:nny
% U(i,j)=U(i,j)+dt*(g(i,j)*I_t(i,j)...
% +max(g_x(i,j),0)*D_x_L(i,j)+min(g_x(i,j),0)*D_x_R(i,j)...
% +max(g_y(i,j),0)*D_y_L(i,j)+min(g_y(i,j),0)*D_y_R(i,j));
% end
% end
U=U+dt*(g.*I_t+max(g_x,0).*D_x_L+min(g_x,0).*D_x_R...
+max(g_y,0).*D_y_L+min(g_y,0).*D_y_R);
%%%%%%%%%% 每400次重新初始化一次 %%%%%%%%%%%%%%%%%%%
if mod(n,400)==0
%%%%%%%%%%% Detect current zero_level_set and display it
curvIndex=zeros(MaxlLengh,2);
num=0;
for i = 2 : nnx - 1
for j = 2 : nny - 1
if U(i,j)<0 & (U(i+1,j)>0 | U(i-1,j)>0 | U(i,j+1)>0 | U(i,j-1)>0)
num=num+1; curvIndex(num,1)=i;curvIndex(num,2)=j;
curvImag(i,j)=0;
end
end
end
figure(5);imshow(uint8(curvImag));
%%%%%%%%%%%%%% Then reinitialize the embedding function U
new_u = zeros(nnx,nny);
dist=zeros(1,num);
for j=1:nny
for i=1:nnx
for k=1:num
dist(k)=sqrt((i-curvIndex(k,1)).^2+(j-curvIndex(k,2)).^2);
end
new_u(i,j)=min(dist);
if U(i,j)<0
new_u(i,j)=-new_u(i,j);
end
end
end
U=new_u;
end % end of if (n mod 100)==0
end % end of iteration