一、简介

1 原理分析\ \ \ \ \ 2 编程过程\ 2.1 准备工作\ \ 2.2 迭代运算\

二、源代码

``` %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clear; clc;

%读图像并将rgb图像转化为灰度图象，为了节省cpu时间，将图像变小为原来大小的1/2

im=imread('3.bmp'); im=rgb2gray(im); im=double(im); im = imresize( im, 0.5 );
figure(1);imshow(uint8(im)); [nrow,ncol]=size(im);

%为了计算函数g，先对图像作预平滑。 J= gauss( im,3,2 ); %guassian平滑

%求图像梯度模值 Jx = (J(:,[2:ncol ncol])-J(:,[1 1:ncol-1]))/2; Jy = (J([2:nrow nrow],:)-J([1 1:nrow-1],:))/2;
gradim = (Jx.^2 + Jy.^2).^0.5;
%计算边缘函数g kk=5; %contrast parameter g=1./(1+(gradim/kk).^2);

%初始化曲线为一个图，初始化u为带符号距离函数 curvIndex=zeros(5nrowncol,2); curvImag=zeros([nrow,ncol]); num=0; for i = 2 : nrow - 1 for j = 2 : ncol - 1 if ((i==8|i==nrow-8)&(j>=8&j<=(ncol-8)))|((i>=8&i<=(nrow-8))&(j==8|j==ncol-8)) num=num+1; InitCurvImag(i,j)=255; curvIndex(num,1)=i; curvIndex(num,2)=j; for k=1:num-1 %check is it the last point of the curve if curvIndex(k,1)==i & curvIndex(k,2)==j num=num-1;break end end end end end

u = zeros(nrow,ncol); dist=zeros(1,num); for j=1:ncol for i=1:nrow for k=1:num dist(k)=sqrt((i-curvIndex(k,1)).^2+(j-curvIndex(k,2)).^2); end u(i,j)=min(dist); if i>8&i<(nrow-8)&j>8&j<(ncol-8) u(i,j)=-u(i,j); end end end

%将当前曲线加入到原图像中,然后写入文件 newim=createimage(im,u,0);
figure(2);imshow(uint8(newim));

%选定迭代步长 delta_t=20; %选定常数速度 c=0.2;

%迭代开始 for iterations=1:4000 %%% compute Roe upwind gradient
uxe = u(:,[2:ncol,ncol])-u; uye = (u([2:nrow,nrow],:)+u([2:nrow,nrow],[2:ncol,ncol])-u([1 1:nrow-1],:)-u([1 1:nrow-1],[2:ncol ncol]))/4; uGe = sqrt(uxe .^2+uye .^2); ge = 0.5*(g(:,[2:ncol,ncol])+g); Terme = ge.*uxe./(uG_e+eps);

uxw = u-u(:,[1 1:ncol-1]); uyw = (u([2:nrow,nrow],:)+u([2:nrow,nrow],[1 1:ncol-1])-u([1,1:nrow-1],:)-u([1,1:nrow-1],[1 1:ncol-1]))/4; uGw = sqrt(uxw.^2+uyw.^2); gw = 0.5*(g(:,[1 1:ncol-1])+g); Termw = gw.*uxw./(uG_w+eps);

uys = u([2:nrow,nrow],:)-u; uxs = (u(:,[2:ncol,ncol])+u([2:nrow,nrow],[2:ncol,ncol])-u(:,[1 1:ncol-1])-u([2:nrow,nrow],[1 1:ncol-1]))/4; uGs = sqrt(uys.^2+ uxs.^2); gs = 0.5*(g([2:nrow,nrow],:)+g); Terms = gs.*uys./(uG_s+eps);

uyn = u-u([1 1:nrow-1],:); uxn = (u(:,[2:ncol,ncol])+u([1 1:nrow-1],[2:ncol,ncol])-u(:,[1 1:ncol-1])-u([1 1:nrow-1],[1 1:ncol-1]))/4; uGn = sqrt(uyn.^2+uxn.^2); gn = 0.5*(g([1,1:nrow-1],:)+g); Termn = gn.*uyn./(uG_n+eps); function Ig=gauss(I,ks,sigma2) %private function: gauss (by Guy Gilboa): % Ig=gauss(I,ks,sigma2) % ks - kernel size (odd number) % sigma2 - variance of Gaussian

[Ny,Nx]=size(I); hks=(ks-1)/2; % half kernel size if (Ny else %% 2-d convolution x=ones(ks,1)*(-hks:hks); y=x'; flt=exp(-(x.^2+y.^2)/(2*sigma2)); % 2D gaussian flt=flt/sum(sum(flt)); % normalize % expand if (hks>1) xL=mean(I(:,1:hks)')'; xR=mean(I(:,Nx-hks+1:Nx)')'; else xL=I(:,1); xR=I(:,Nx); end eI=[xLones(1,hks) I xRones(1,hks)]; if (hks>1) xU=mean(eI(1:hks,:)); xD=mean(eI(Ny-hks+1:Ny,:));
else xU=eI(1,:); xD=eI(Ny,:);
end ```

三、运行结果

\ \

四、备注