happy 发表于 2006-9-6 15:42

水平集算法matlab实现

function phiy = activecontourCV( u0, center, radius, isinside, d_it, m_it, m_name )
% 用主动轮廓线CV算法对输入图像u0实现图像边缘提取
% 输入图像为double型,灰度为1—256的图像。选用圆形起始轮廓线
% center为起始轮廓线原点,radius为起始轮廓线半径。isinside 表示边缘目标在起始轮廓线外还是内,=1表示目标在起始轮廓线内,=0表示在外

% 初始化参数
ITERATIONS = 500;%迭代次数
delta_t = 0.1;%时间步长
%轮廓内外能量参数
lambda1 = 1;
lambda2 = 1;
nu = 0;
%曲率项参数
h = 1; h_sq = h^2;
epsilon = 1;
mu = 0.01 * 255^2;

% 初始化符号距离函数
phi = initsdf( size( u0 ), center, radius, isinside );

for ii = 1 : ITERATIONS;

% 显示当前迭代次数
fprintf( 1, '%d\n', ii );

% 每d_it显示一次图像
if( mod( ii - 1, d_it ) == 0 )
    disp( 'Displaying Segmented Image' );
    segim = createim( u0, phi );
    clf; imshow( segim );
    drawnow;
end;

% 每m_it次保存一次图像
   if( mod( ii - 1, m_it ) == 0 )
    segim = createim( u0, phi );
    filename = strcat( m_name, sprintf( '%06d', ( ( ii - 1 )/ m_it ) + 1 ), '.png' );
    imwrite( segim, filename );
end;

%delta_t*狄力克函数
dirac_delta_t = delta_t * diracfunction( phi, epsilon );

% 计算轮廓线内外能量
[ inside, outside ] = calcenergyf( u0, phi, epsilon );
energy_term = -nu - lambda1 .* inside + lambda2 .* outside;%能量项

% 中心差分
dx_central = ( circshift( phi, [ 0, -1 ] ) - circshift( phi, [ 0, 1 ] ) ) / 2;
dy_central = ( circshift( phi, [ -1, 0 ] ) - circshift( phi, [ 1, 0 ] ) ) / 2;

% div(delt_phi/|delta_phi|)
abs_grad_phi = ( sqrt( dx_central.^2 + dy_central.^2 ) + 0.00001 );
x = dx_central ./ abs_grad_phi;
y = dy_central ./ abs_grad_phi;
grad_term = ( mu / h_sq ) .* divergence( x, y );%梯度能量

% phi(n+1)
phi = phi + dirac_delta_t .* ( grad_term + energy_term );
phiy(:,:,ii)=phi;%保存每次迭代结果
end;


%子函数
%初始化水平集函数,符号距离函数
function phi = initsdf( imsize, center, radius, isinside )
m = imsize( 1 ); n = imsize( 2 );
phi = zeros( imsize );
for i = 1 : m;
for j = 1 : n;
   distance = sqrt( sum( ( center - [ i, j ] ).^2 ) );
   phi( i, j ) = distance - radius;
    if( isinside == 0 )
      phi( i, j ) = -phi( i, j );
    end
end
end

%狄力克函数,利于检测出目标的空心部分
function y = diracfunction( x, epsilon )
y = ( 1 ./ pi ) .* ( epsilon ./ ( epsilon.^2 + x.^2 ) );

%用于求轮廓内外能量
function value = heavisidef( z, epsilon )
value = 0.5 .* ( 1 + ( 2 ./ pi ) .* atan( z ./ epsilon ) );

%求轮廓内外能量
function [ inside, outside ] = calcenergyf( u0, phi, epsilon )
H_phi = heavisidef( phi, epsilon );
H_phi_minus = 1 - heavisidef( phi, epsilon );

c1 = sum( sum( u0 .* H_phi ) ) /sum( sum( H_phi ) );
c2 = sum( sum( u0 .* H_phi_minus ) ) / sum( sum( H_phi_minus ) );

inside = ( u0 - c1 ).^2;
outside = ( u0 - c2 ).^2;

%显示代轮廓的图像,轮廓为零水平集,用ifro检测出。将原图和轮廓叠加为rgb图像,轮廓用绿色表示
function newim = createim( im, phi )
newim( :, :, 1 ) = im;
newim( :, :, 3 ) = im;

tempim = im;
tempim( find( isfro( phi ) ) ) = 255;

newim( :, :, 2 ) = tempim;

newim = uint8( newim );

%检测零水平集的位置
function front = isfro( phi )
[ n, m ] = size( phi );
front = zeros( size( phi ) );
for i = 2 : n - 1;
for j = 2 : m - 1;

    maxVal = max( max( phi( i:i+1, j:j+1 ) ) );
    minVal = min( min( phi( i:i+1, j:j+1 ) ) );
    front( i, j ) = ( ( maxVal > 0 ) & ( minVal < 0 ) ) | phi( i, j ) == 0;

end
end

mathfish 发表于 2006-9-8 14:44

谢谢!不过水平集理解起来还是比较困难的

wangyongq246 发表于 2006-9-14 10:08

我同学一直想要这方面的内容,这回可以帮他了!!

bolandi 发表于 2011-2-9 19:36

程序怎么使用啊?

bolandi 发表于 2011-2-9 19:37

m_name是什么啊?
页: [1]
查看完整版本: 水平集算法matlab实现