activecontourpaper.m

这是一个利用水平集方法进行分割的方法

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;