水平集算法matlab实现

happy

1. function phiy = activecontourCV( u0, center, radius, isinside, d_it, m_it, m_name )
   
2. % 用主动轮廓线CV算法对输入图像u0实现图像边缘提取
   
3. % 输入图像为double型，灰度为1—256的图像。选用圆形起始轮廓线
   
4. % center为起始轮廓线原点，radius为起始轮廓线半径。isinside 表示边缘目标在起始轮廓线外还是内，=1表示目标在起始轮廓线内，=0表示在外
   
5. 
   
6. % 初始化参数
   
7. ITERATIONS = 500;%迭代次数
   
8. delta_t = 0.1;%时间步长
   
9. %轮廓内外能量参数
   
10. lambda1 = 1;
    
11. lambda2 = 1;
    
12. nu = 0;
    
13. %曲率项参数
    
14. h = 1; h_sq = h^2;
    
15. epsilon = 1;
    
16. mu = 0.01 * 255^2;
    
17. 
    
18. % 初始化符号距离函数
    
19. phi = initsdf( size( u0 ), center, radius, isinside );
    
20. 
    
21. for ii = 1 : ITERATIONS;
    
22. 
    
23.   % 显示当前迭代次数
    
24.   fprintf( 1, '%d\n', ii );
    
25. 
    
26.   % 每d_it显示一次图像
    
27.   if( mod( ii - 1, d_it ) == 0 )
    
28.     disp( 'Displaying Segmented Image' );
    
29.     segim = createim( u0, phi );
    
30.     clf; imshow( segim );
    
31.     drawnow;
    
32.   end;
    
33.   
    
34.   % 每m_it次保存一次图像
    
35.    if( mod( ii - 1, m_it ) == 0 )
    
36.     segim = createim( u0, phi );
    
37.     filename = strcat( m_name, sprintf( '%06d', ( ( ii - 1 )/ m_it ) + 1 ), '.png' );
    
38.     imwrite( segim, filename );
    
39.   end;
    
40. 
    
41.   %delta_t*狄力克函数
    
42.   dirac_delta_t = delta_t * diracfunction( phi, epsilon );
    
43. 
    
44.   % 计算轮廓线内外能量
    
45.   [ inside, outside ] = calcenergyf( u0, phi, epsilon );
    
46.   energy_term = -nu - lambda1 .* inside + lambda2 .* outside;%能量项
    
47. 
    
48.   % 中心差分
    
49.   dx_central = ( circshift( phi, [ 0, -1 ] ) - circshift( phi, [ 0, 1 ] ) ) / 2;
    
50.   dy_central = ( circshift( phi, [ -1, 0 ] ) - circshift( phi, [ 1, 0 ] ) ) / 2;
    
51. 
    
52.   % div(delt_phi/|delta_phi|)
    
53.   abs_grad_phi = ( sqrt( dx_central.^2 + dy_central.^2 ) + 0.00001 );
    
54.   x = dx_central ./ abs_grad_phi;
    
55.   y = dy_central ./ abs_grad_phi;
    
56.   grad_term = ( mu / h_sq ) .* divergence( x, y );%梯度能量
    
57. 
    
58.   % phi(n+1)
    
59.   phi = phi + dirac_delta_t .* ( grad_term + energy_term );
    
60.   phiy(:,:,ii)=phi;%保存每次迭代结果
    
61. end;
    
62. 
    
63. 
    
64. %子函数
    
65. %初始化水平集函数，符号距离函数
    
66. function phi = initsdf( imsize, center, radius, isinside )
    
67. m = imsize( 1 ); n = imsize( 2 );
    
68. phi = zeros( imsize );
    
69. for i = 1 : m;
    
70.   for j = 1 : n;
    
71.      distance = sqrt( sum( ( center - [ i, j ] ).^2 ) );
    
72.      phi( i, j ) = distance - radius;
    
73.     if( isinside == 0 )
    
74.       phi( i, j ) = -phi( i, j );
    
75.     end
    
76.   end
    
77. end
    
78. 
    
79. %狄力克函数，利于检测出目标的空心部分
    
80. function y = diracfunction( x, epsilon )
    
81. y = ( 1 ./ pi ) .* ( epsilon ./ ( epsilon.^2 + x.^2 ) );
    
82. 
    
83. %用于求轮廓内外能量
    
84. function value = heavisidef( z, epsilon )
    
85. value = 0.5 .* ( 1 + ( 2 ./ pi ) .* atan( z ./ epsilon ) );
    
86. 
    
87. %求轮廓内外能量
    
88. function [ inside, outside ] = calcenergyf( u0, phi, epsilon )
    
89. H_phi = heavisidef( phi, epsilon );
    
90. H_phi_minus = 1 - heavisidef( phi, epsilon );
    
91. 
    
92. c1 = sum( sum( u0 .* H_phi ) ) /  sum( sum( H_phi ) );
    
93. c2 = sum( sum( u0 .* H_phi_minus ) ) / sum( sum( H_phi_minus ) );
    
94. 
    
95. inside = ( u0 - c1 ).^2;
    
96. outside = ( u0 - c2 ).^2;
    
97. 
    
98. %显示代轮廓的图像，轮廓为零水平集，用ifro检测出。将原图和轮廓叠加为rgb图像，轮廓用绿色表示
    
99. function newim = createim( im, phi )
    
100. newim( :, :, 1 ) = im;
     
101. newim( :, :, 3 ) = im;
     
102. 
     
103. tempim = im;
     
104. tempim( find( isfro( phi ) ) ) = 255;
     
105. 
     
106. newim( :, :, 2 ) = tempim;
     
107. 
     
108. newim = uint8( newim );
     
109. 
     
110. %检测零水平集的位置
     
111. function front = isfro( phi )
     
112. [ n, m ] = size( phi );
     
113. front = zeros( size( phi ) );
     
114. for i = 2 : n - 1;
     
115.   for j = 2 : m - 1;
     
116. 
     
117.     maxVal = max( max( phi( i:i+1, j:j+1 ) ) );
     
118.     minVal = min( min( phi( i:i+1, j:j+1 ) ) );
     
119.     front( i, j ) = ( ( maxVal > 0 ) & ( minVal < 0 ) ) | phi( i, j ) == 0;
     
120. 
     
121.   end
     
122. end

复制代码

mathfish

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

wangyongq246

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

bolandi

程序怎么使用啊？

bolandi

m_name是什么啊？