用matlab编写维尔斯特拉斯图像的程序

elseif nargin<2

N=128

end

image_input=imread(image_name)

subplot(2,2,1)

imshow(image_input)

title('原图像')

[size_m,size_n]=size(image_input)

matrix_temp=zeros(size_m,size_n)

% 求水平集

for row=1:size_m

for col=1:size_n

if image_input(row,col) >N

matrix_temp(row,col)=1

end

end

end

subplot(2,2,2)

imshow(matrix_temp,[])

title('图像的水平集')

imwrite(matrix_temp,'level_setzhan.bmp')

% 图像矩阵扩展 赋值 便于处理边界

matrix_ex=zeros(size_m+2,size_n+2)

for row=1:size_m

for col=1:size_n

matrix_ex(row+1,col+1)=matrix_temp(row,col)

end

end

% 四邻域反填充 得水平线

matrix_new=matrix_temp

for row=2:size_m+1

for col=2:size_n+1

if matrix_ex(row+1,col)==0 &matrix_ex(row-1,col)==0 &matrix_ex(row,col+1)==0 &matrix_ex(row,col-1)==0

matrix_new(row-1,col-1)=1

end

end

end

subplot(2,2,3)

imshow(matrix_new,[])

title('图像的水平线')

%imwrite(matrix_new,'level_line.bmp')

% 求图像的等高线

contour=zeros(size_m,size_n)

for row=1:size_m

for col=1:size_n

if image_input(row,col)==N

contour(row,col)=1

end

contour(row,col)=1-contour(row,col)

end

end

subplot(2,2,4)

imshow(contour,[])

title('图像的等高线')

%imwrite(contour,'contour.bmp')

function levelset2image(out_filename)

% 水平集检验 体现图象与水平集的关系

% out_filename: 输出文件名

cd('level_set')

% 生成各层水平集－－－－－－－－－－－－－－－－－－－－－－－

temp=imread('level_set_0.bmp')

[size_r,size_c]=size(temp)

level_image=zeros(size_r,size_c,255)

for N=0:254

level_image(:,:,N+1)=imread(strcat('level_set_',num2str(N),'.bmp'))

end

% 求各层矩阵最大值

conduct_image=zeros(size_r,size_c)

% 由各层水平集重构图像－－－－－－－－－－－－－－－－－－－－

for row=1:size_r

for col=1:size_c

for i=255:-1:1

if level_image(row,col,i)==255

conduct_image(row,col)=i %%%%%%

break

end

end

end

end

imshow(conduct_image,[])

title('重构出的图像')

imwrite(uint8(conduct_image),out_filename)

cd('..')

function u = EVOLUTION(u0, g, lambda, mu, alf, epsilon, delt, numIter)

% EVOLUTION(u0, g, lambda, mu, alf, epsilon, delt, numIter) updates the level set function

% according to the level set evolution equation in Chunming Li et al's paper:

% "Level Set Evolution Without Reinitialization: A New Variational Formulation"

% in Proceedings CVPR'2005,

% Usage:

% u0: level set function to be updated

% g: edge indicator function

% lambda: coefficient of the weighted length term L(\phi)

% mu: coefficient of the internal (penalizing) energy term P(\phi)

% alf: coefficient of the weighted area term A(\phi), choose smaller alf

% epsilon: the papramater in the definition of smooth Dirac function, default value 1.5

% delt: time step of iteration, see the paper for the selection of time step and mu

% numIter: number of iterations.

%

u=u0

[vx,vy]=gradient(g)

for k=1:numIter

u=NeumannBoundCond(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 EVOLUTION

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 = NeumannBoundCond(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])