www.pudn.com > level_set_methods_1.1.zip > evolve2D.m, change:2005-05-12,size:11010b

```function [phi] = evolve2D(phi, dx, dy, alpha, iterations, accuracy, is_signed_distance, normal_evolve, Vn, vector_evolve, u, v, kappa_evolve, b)
%
% function [phi] = evolve2D(phi, dx, dy, iterations, accuracy, ...
%   is_signed_distance, normal_evolve, Vn, vector_evolve, u, v, kappa_evolve, b)
%
% Calculates evolution for a 2D curve (3D level set function) phi.
% phi is the input level set function.
%
% dx and dy are the resolution of the grid at x and y dimensions.
% alpha is a constant for calculating the euler step (dt). Should
% be between 0 and 1. 0.5 is quite safe whereas 0.9 can be risky.
% iterations specifies the number of iterations before the function returns.
% normal_evolve, vector_evolve, kappa_evolve should be either set to 0 and 1. This
% indicates if these forces are present or not. If any of these are set to 1,
% corresponding  variables (e.g. Vn, u, v, b) should not be an empty array []. if
% either normal_evolve or vector_evolve are set to 1, accuracy needs to be specified,
% otherwise accuracy can be set to an empty array. Allowed values for accuracy are
% 'ENO1', 'ENO2', 'ENO3', 'WENO'. These correspond to 1st, 2nd, 3rd and 5th order
% accurate schemes for calculating the derivative of phi. if both normal_evolve and
% vector_evolve are set to 1, then is_signed_distance needs to be specified. If phi is
% approximately a signed distance function, set this variable to 1. If
% is_signed_distance == 1, Godunov scheme will be used, otherwise Stencil Local
% Lax-Friedrichs (SLLF) scheme will be used (See Osher&Fedkiw section 6.4).
%
% Other variables (these are either a scalar or of the same size as phi):
% Vn: Force in the normal direction.
% u: x component of the vector field
% v: y component of the vector field
% b: this specifies the weighting for the curvature (always positive everywhere).
%
% Author: Baris Sumengen  sumengen@ece.ucsb.edu
% http://vision.ece.ucsb.edu/~sumengen/
%

tic

if alpha <= 0 | alpha >= 1
error('alpha needs to be between 0 and 1!');
end

% expand parameters if they are scalar
if length(Vn) == 1
Vn = Vn.*ones(size(phi));
end
if length(u) == 1
u = u.*ones(size(phi));
end
if length(v) == 1
v = v.*ones(size(phi));
end
if length(b) == 1
b = b.*ones(size(phi));
end

% Understand what kind of evolution the user is interested in.
if (normal_evolve == 1) & (vector_evolve == 1) & (kappa_evolve == 1)
evolution_type = 'normal_vector_kappa';
if isempty(Vn) | isempty(u) | isempty(v) | isempty(b)
error('Some of the variables should not be empty!')
end
switch(accuracy)
case 'ENO1'
init_normal = @init_normal_ENO1;
init_vector = @init_vector_ENO1;
if is_signed_distance == 1
evolve_normal_vector = @evolve_normal_vector_ENO1_SD;
else
evolve_normal_vector = @evolve_normal_vector_ENO1;
end
case 'ENO2'
init_normal = @init_normal_ENO2;
init_vector = @init_vector_ENO2;
if is_signed_distance == 1
evolve_normal_vector = @evolve_normal_vector_ENO2_SD;
else
evolve_normal_vector = @evolve_normal_vector_ENO2;
end
case 'ENO3'
init_normal = @init_normal_ENO3;
init_vector = @init_vector_ENO3;
if is_signed_distance == 1
evolve_normal_vector = @evolve_normal_vector_ENO3_SD;
else
evolve_normal_vector = @evolve_normal_vector_ENO3;
end
case 'WENO'
init_normal = @init_normal_WENO;
init_vector = @init_vector_WENO;
if is_signed_distance == 1
evolve_normal_vector = @evolve_normal_vector_WENO_SD;
else
evolve_normal_vector = @evolve_normal_vector_WENO;
end
otherwise
error('Desired type of the accuracy is not correctly specified!');
end
elseif (normal_evolve == 1) & (vector_evolve == 1) & (kappa_evolve == 0)
evolution_type = 'normal_vector';
if isempty(Vn) | isempty(u) | isempty(v)
error('Some of the variables should not be empty!')
end
switch(accuracy)
case 'ENO1'
init_normal = @init_normal_ENO1;
init_vector = @init_vector_ENO1;
if is_signed_distance == 1
evolve_normal_vector = @evolve_normal_vector_ENO1_SD;
else
evolve_normal_vector = @evolve_normal_vector_ENO1;
end
case 'ENO2'
init_normal = @init_normal_ENO2;
init_vector = @init_vector_ENO2;
if is_signed_distance == 1
evolve_normal_vector = @evolve_normal_vector_ENO2_SD;
else
evolve_normal_vector = @evolve_normal_vector_ENO2;
end
case 'ENO3'
init_normal = @init_normal_ENO3;
init_vector = @init_vector_ENO3;
if is_signed_distance == 1
evolve_normal_vector = @evolve_normal_vector_ENO3_SD;
else
evolve_normal_vector = @evolve_normal_vector_ENO3;
end
case 'WENO'
init_normal = @init_normal_WENO;
init_vector = @init_vector_WENO;
if is_signed_distance == 1
evolve_normal_vector = @evolve_normal_vector_WENO_SD;
else
evolve_normal_vector = @evolve_normal_vector_WENO;
end
otherwise
error('Desired type of the accuracy is not correctly specified!');
end
elseif (normal_evolve == 1) & (vector_evolve == 0) & (kappa_evolve == 1)
evolution_type = 'normal_kappa';
if isempty(Vn) | isempty(b)
error('Some of the variables should not be empty!')
end
switch(accuracy)
case 'ENO1'
init_normal = @init_normal_ENO1;
evolve_normal = @evolve_normal_ENO1;
case 'ENO2'
init_normal = @init_normal_ENO2;
evolve_normal = @evolve_normal_ENO2;
case 'ENO3'
init_normal = @init_normal_ENO3;
evolve_normal = @evolve_normal_ENO3;
case 'WENO'
init_normal = @init_normal_WENO;
evolve_normal = @evolve_normal_WENO;
otherwise
error('Desired type of the accuracy is not correctly specified!');
end
elseif (normal_evolve == 1) & (vector_evolve == 0) & (kappa_evolve == 0)
evolution_type = 'normal';
if isempty(Vn)
error('Some of the variables should not be empty!')
end
switch(accuracy)
case 'ENO1'
init_normal = @init_normal_ENO1;
evolve_normal = @evolve_normal_ENO1;
case 'ENO2'
init_normal = @init_normal_ENO2;
evolve_normal = @evolve_normal_ENO2;
case 'ENO3'
init_normal = @init_normal_ENO3;
evolve_normal = @evolve_normal_ENO3;
case 'WENO'
init_normal = @init_normal_WENO;
evolve_normal = @evolve_normal_WENO;
otherwise
error('Desired type of the accuracy is not correctly specified!');
end
elseif (normal_evolve == 0) & (vector_evolve == 1) & (kappa_evolve == 1)
evolution_type = 'vector_kappa';
if isempty(u) | isempty(v) | isempty(b)
error('Some of the variables should not be empty!')
end
switch(accuracy)
case 'ENO1'
init_vector = @init_vector_ENO1;
evolve_vector = @evolve_vector_ENO1;
case 'ENO2'
init_vector = @init_vector_ENO2;
evolve_vector = @evolve_vector_ENO2;
case 'ENO3'
init_vector = @init_vector_ENO3;
evolve_vector = @evolve_vector_ENO3;
case 'WENO'
init_vector = @init_vector_WENO;
evolve_vector = @evolve_vector_WENO;
otherwise
error('Desired type of the accuracy is not correctly specified!');
end
elseif (normal_evolve == 0) & (vector_evolve == 1) & (kappa_evolve == 0)
evolution_type = 'vector';
if isempty(u) | isempty(v)
error('Some of the variables should not be empty!')
end
switch(accuracy)
case 'ENO1'
init_vector = @init_vector_ENO1;
evolve_vector = @evolve_vector_ENO1;
case 'ENO2'
init_vector = @init_vector_ENO2;
evolve_vector = @evolve_vector_ENO2;
case 'ENO3'
init_vector = @init_vector_ENO3;
evolve_vector = @evolve_vector_ENO3;
case 'WENO'
init_vector = @init_vector_WENO;
evolve_vector = @evolve_vector_WENO;
otherwise
error('Desired type of the accuracy is not correctly specified!');
end
elseif (normal_evolve == 0) & (vector_evolve == 0) & (kappa_evolve == 1)
evolution_type = 'kappa';
if isempty(b)
error('Some of the variables should not be empty!')
end
else
error('Incorrect input parameters! normal_evolve, vector_evolve, kappa_evolve can only be 0 or 1 and all of them cannot be 0 at the same time');
end

% Now the main part of the program: Evolution
disp('Evolution type: ');
disp(evolution_type);
switch(evolution_type)
case 'normal_vector_kappa'
Vn_ext = feval(init_normal, Vn);
[u_ext, v_ext] = feval(init_vector, u, v);
[dx2, dy2] = init_kappa(dx, dy);
it=0;
t = 0;
while(it < iterations)
[delta_normal_vector, H1_abs, H2_abs] = feval(evolve_normal_vector, ...
phi, dx, dy, Vn_ext, u_ext, v_ext);
delta_kappa = evolve_kappa(phi, dx, dy, b, dx2, dy2);
dt = get_dt_normal_vector_kappa(alpha, dx, dy, H1_abs, H2_abs, b, dx2, dy2);
phi = phi + dt*(delta_kappa - delta_normal_vector);
it = it+1
t = t+dt;
end
case 'normal_vector'
Vn_ext = feval(init_normal, Vn);
[u_ext, v_ext] = feval(init_vector, u, v);
it=0;
t = 0;
while(it < iterations)
[delta_normal_vector, H1_abs, H2_abs] = feval(evolve_normal_vector, ...
phi, dx, dy, Vn_ext, u_ext, v_ext);
dt = get_dt_normal_vector(alpha, dx, dy, H1_abs, H2_abs);
phi = phi - dt*delta_normal_vector;
it = it+1
t = t+dt;
end
case 'normal_kappa'
Vn_ext = feval(init_normal, Vn);
[dx2, dy2] = init_kappa(dx, dy);
it=0;
t = 0;
while(it < iterations)
[delta_normal, H1_abs, H2_abs] = feval(evolve_normal, phi, dx, dy, Vn_ext);
delta_kappa = evolve_kappa(phi, dx, dy, b, dx2, dy2);
dt = get_dt_normal_kappa(alpha, dx, dy, H1_abs, H2_abs, b, dx2, dy2);
phi = phi + dt*(delta_kappa - delta_normal);
it = it+1
t = t+dt;
end
case 'normal'
Vn_ext = feval(init_normal, Vn);
it=0;
t = 0;
while(it < iterations)
[delta_normal, H1_abs, H2_abs] = feval(evolve_normal, phi, dx, dy, Vn_ext);
dt = get_dt_normal(alpha, dx, dy, H1_abs, H2_abs);
phi = phi - dt*delta_normal;
it = it+1
t = t+dt;
end
case 'vector_kappa'
[u_ext, v_ext] = feval(init_vector, u, v);
[dx2, dy2] = init_kappa(dx, dy);
dt = get_dt_vector_kappa(alpha, dx, dy, u, v, b, dx2, dy2);
it=0;
t = 0;
while(it < iterations)
delta_vector = feval(evolve_vector, phi, dx, dy, u_ext, v_ext);
delta_kappa = evolve_kappa(phi, dx, dy, b, dx2, dy2);
phi = phi + dt*(delta_kappa - delta_vector);
it = it+1
t = t+dt;
end
case 'vector'
[u_ext, v_ext] = feval(init_vector, u, v);
dt = get_dt_vector(alpha, dx, dy, u, v);
it=0;
t = 0;
while(it < iterations)
delta_vector = feval(evolve_vector, phi, dx, dy, u_ext, v_ext);
phi = phi - dt*delta_vector;
it = it+1
t = t+dt;
end
case 'kappa'
[dx2, dy2] = init_kappa(dx, dy);
dt = get_dt_kappa(alpha, dx, dy, b, dx2, dy2);
it=0;
t = 0;
while(it < iterations)
delta_kappa = evolve_kappa(phi, dx, dy, b, dx2, dy2);
phi = phi + dt*delta_kappa;
it = it+1;
t = t+dt;
end
otherwise
error('Something went wrong when checking the evolution type!');
end

t

toc

```