89815618levelsetImage.rar_level set_图像分割_水平集_水平集方法

function phi = activecontourpaper( u0, center, radius, isinside, d_it, m_it, m_name ) % ACTIVECONTOURPAPER Segment an image using active contours % ACTIVECONTOURPAPER( u0, center, radius, isinside, d_it, m_it, m_name ) % segments the image u0 with initial segmented region as a % circle with argument center and radius. The current segmented % image is displayed every d_it iterations and written to disk % every m_iterations with name m_name*.png. Uses Jacobi Method. % setup constants ITERATIONS = 25000; INNER_ITERATIONS = 100; delta_t = 0.1; lambda1 = 1; lambda2 = 1; nu = 0; h = 1; h_sq = h^2; epsilon = 1; mu = 0.0001 * 255^2; % initialize phi to signed distance function from circle phi = initphi( size( u0 ), center, radius, isinside ); for ii = 1 : ITERATIONS; % display current iteration fprintf( 1, '%d\n', ii ); % display the segmented image every 'd_it' iterations if( mod( ii - 1, d_it ) == 0 ) disp( 'Displaying Segmented Image' ); segim = createimage( u0, phi ); clf; imshow( segim ); drawnow; end; % write current segmented image to file every 'm_it' % iterations if( mod( ii - 1, m_it ) == 0 ) segim = createimage( u0, phi ); filename = strcat( m_name, sprintf( '%06d', ( ( ii - 1 )/ m_it ) + 1 ), '.png' ); imwrite( segim, filename ); end; % compute dirac(phi) * delta_t factor dirac_delta_t = delta_t * dirac( phi, epsilon ); % calculate inside and outside curve terms [ inside, outside ] = calcenergy( u0, phi, epsilon ); energy_term = -nu - lambda1 .* inside + lambda2 .* outside; % calculate forward, backward and central differences dx_plus = circshift( phi, [ 0, -1 ] ) - phi; dy_plus = circshift( phi, [ -1, 0 ] ) - phi; dx_minus = phi - circshift( phi, [ 0, 1 ] ); dy_minus = phi - circshift( phi, [ 1, 0 ] ); dx_central = ( dx_plus + dx_minus ) ./ 2; dy_central = ( dy_plus + dy_minus ) ./ 2; % calculate sqrt terms sqrt_term1 = sqrt( ( dx_plus.^2 ./ h_sq ) + ( dy_central.^2 ./ h_sq ) ); sqrt_term2 = sqrt( ( dx_central.^2 ./ h_sq ) + ( dy_plus.^2 ./ h_sq ) ); % calculate main terms mult_term = mu / h_sq; first_term = mult_term ./ sqrt_term1; second_term = mult_term ./ sqrt_term2; % intialize phi_k to phi phi_k = phi; % Use Jacobi method to estimate forward differences of phi^n+1 for k = 1 : INNER_ITERATIONS; % determine forward differences dx_plus_k = circshift( phi_k, [ 0, -1 ] ) - phi_k; dy_plus_k = circshift( phi_k, [ -1, 0 ] ) - phi_k; grad_term = dx_plus_k .* first_term + dy_plus_k .* second_term; % calculate phi_k+1 old_phi_k = phi_k; phi_k = phi + dirac_delta_t .* ( ( mult_term .* grad_term ) + energy_term ); % if we have converged then break if( max( max( abs( old_phi_k - phi_k ) ) ) <= 0.000001 ) break; end; end; % update phi^n to phi^n+1 phi = phi_k; end;