一、简介
1 基本思想： 令（n+1)维函数 t = phi(x(t)) 等于0，即水平集。
2 基本问题：
2.1 初始化
定义 level set 函数 phi(x) 为符号距离函数，界面线上距离为0。一般为隐函数。
2.1.1 符号判定， 界面闭合线内/外
2.1.1.1 2D留数算法： 将边界C离散成n段，c_1, c_2, … c_n 为结点， oc_i 表示某点o 指向c_i的向量
sum ( atan(( oc_i x oc_(i+1)) / (oc_i . oc_(i+1))), i = 1: n)
if sum == 0 —> 点o 在边界C外
if sum == 2 * pi —> 点o在边界C内
2.1.1.2 交点判定算法: 给定已知符号的点P0, 待定符号的点Px 与 P0 之间连线穿过边界C，
如果穿过次数是偶数表示Px与P0符号相同；
如果是奇数，则符号相反。
实现：
将边界C离散成n段，初始定义所有网格点（2D面）符号为正（C外），当离散的边界线段，与某条垂向网格线相交，则将该网格线上，在边界线段上侧的点全部变负号。
2.1.2 距离计算
2.1.2.1 边界C可由参数方程表示，先搜索距离边界线最近的网格点，计算距离；
2.1.2.2 边界C分段连续，直接计算各个端点到网格点的距离
2.2 演化方程
2.2.1 边界线随时间的变化，采用PDE描述 phi(x(t), t) = 0
d phi(x)/ dt + grad_phi (a. n + b.t) = 0
其中n是法向单位向量，t是切向单位向量。a是法向速度，b是切向速度。n影响边界线拓扑。
<- - n = grad_ phi / norm( grad_phi)
–> d phi(x)/ dt = a norm(grad_phi) = 0
2. 2.2 演化方程的数值解
演化方程具有 Hamilton-Jacobi形式， d phi /dt + H(grad_phi) = 0 。 这类方程不存在经典解，只有弱解。解C0连续但不保证C1连续，即使初值C1连续。守恒型双曲PDE也有类似的解性质。比如超声速压缩流动，给定连续的初状态值，仍然可能演化出不连续的间断激波。所以演化方程的高阶精度数值解可借鉴守恒型双曲方程的ENO， WENO等解方案。LSM演化方程实际上是就是一个包含源项的守恒型双曲方程。
常用一阶精度的演化方程常采用迎风格式，Godunov格式， Lax_Friedrich格式等。
2.2.2.1 2D Godunov格式
norm(grad_phi) ~= ( max(-D_{+x}，D_{-x}, 0) ^2 + max(-D_{+y}, D_{-y}, 0)^2 ) ^.5;
D_{+x} forward差分, D_{-x} backward差分。同理 D_{+y}, D_{-y}。 物理意义上要求保证信息必须从已知传向未知，所以对于每个当前要求的phi， 必须选择传播正方向上的差分格式。
一个简单的判断：当地(local) phi_x (i,j) 大于还是小于0. 如果 大于0，表示phi(i,j)沿x轴正方向传播，显然先要经过phi(i-1, j)，故应该选用 D_{-x}差分。而且此时还满足
abs( D_{-x}) > abs(D_{+x})， 所以才有统一的上式样子。同理，如果local phi_x(i,j) 小于0， 则表示 phi(i,j )沿x轴负方向传播，即先经过 phi(i+1, j)， 故采用 D_{+x}。同理y方向。
MATALB CODE here
2. 2.2.2 ENO格式
参见 力学所李新亮老师CFD讲义， 激波平滑。
2.2 法向速度延扩
从边界线到全水平集面，任意垂直边界线的直线上速度恒定 --》 所以水平集面上任意一点的速度，可由两条垂直边界线的直线l1, l2唯一确定
2.3 边界设定
布置ghost point， 速度边界Ok, 位移边界可能出问题 – > 重初始化(重新选择一个level set funtion)
假定周期边界 --》 保证单元边界符号连续
2.3 窄带level set method
只跟新距离边界线一定范围的网格点，远场点保持不变。节省计算/存储。
二、源代码
unction [phi L_z] = ls_sparse(phi, L_z, h, iter)
% Whitaker "A Level-Set Approach to 3D Reconstruction..." IJCV 1998
CFL = 0.15;
% shift and peek operations
peekD = @(M,rr,cc) M(safe_sub2ind(size(phi), rr-1, cc));
peekU = @(M,rr,cc) M(safe_sub2ind(size(phi), rr+1, cc));
peekR = @(M,rr,cc) M(safe_sub2ind(size(phi), rr, cc-1));
peekL = @(M,rr,cc) M(safe_sub2ind(size(phi), rr, cc+1));
shiftD = @(M) M([1 1:end-1],:);
shiftL = @(M) M(:,[2:end end]);
shiftR = @(M) M(:,[1 1:end-1]);
shiftU = @(M) M([2:end end],:);
%- gather layers
Nd = shiftD(phi); Nu = shiftU(phi); Nr = shiftR(phi); Nl = shiftL(phi);
L_i = find(phi < -.5 & (Nu >= -.5 | Nd >= -.5 | Nl >= -.5 | Nr >= -.5));
L_o = find(phi > .5 & (Nu <= .5 | Nd <= .5 | Nl <= .5 | Nr <= .5));
%- layer reverse lookup
L_rev = 2*ones(size(phi), 'int8'); % default: outside
L_rev(phi < 0) = -2; % inside
L_rev(L_i) = -1;
L_rev(L_o) = 1;
L_rev(L_z) = 0;
for i = 1:iter
% record sign
L_z_ = L_z;
L_z_phi_ = phi(L_z);
%- 1. update value of zero layer (with uncontrollable force)
F = h.init_iteration(phi, L_z);
phi(L_z) = phi(L_z) - CFL*F/max(abs(F)); % control that force
%- 3. update status of zero layer
S_i = L_z(phi(L_z) < -.5 );
S_o = L_z( .5 < phi(L_z));
L_z(phi(L_z) < -.5 | .5 < phi(L_z)) = []; % drop
%- 2. update value of nonzero layers
[rr cc] = ind2sub(size(phi), L_i);
N = [peekR(phi,rr,cc) peekL(phi,rr,cc) peekU(phi,rr,cc) peekD(phi,rr,cc)]';
phi(L_i) = max(N) - 1;
[rr cc] = ind2sub(size(phi), L_o);
N = [peekR(phi,rr,cc) peekL(phi,rr,cc) peekU(phi,rr,cc) peekD(phi,rr,cc)]';
phi(L_o) = min(N) + 1;
%- 4. update status of nonzero layers
S_zi = L_i( -.5 <= phi(L_i));
S_zo = L_o(phi(L_o) <= .5 );
L_rev(L_i(phi(L_i) < -1.5)) = -2;
L_rev(L_o(phi(L_o) > 1.5)) = 2;
L_i(phi(L_i) < -1.5 | -.5 <= phi(L_i)) = [];
L_o(phi(L_o) <= .5 | 1.5 < phi(L_o)) = [];
%- 5. update reverse lookup
L_rev(S_i) = -1;
L_rev(S_o) = 1;
L_rev([S_zi; S_zo]) = 0;
%- 7. add new neighbors
Nl = shiftL(L_rev); Nr = shiftR(L_rev);
Nu = shiftU(L_rev); Nd = shiftD(L_rev);
is_near = Nl==0 | Nr==0 | Nu==0 | Nd==0;
S_oo = find(L_rev == 2 & is_near);
S_ii = find(L_rev == -2 & is_near);
L_o = [L_o; S_o; S_oo];
L_i = [L_i; S_i; S_ii];
L_z = [L_z; S_zi; S_zo];
L_rev(S_ii) = -1;
L_rev(S_oo) = 1;
%- 8. update values of new neighbors
[rr cc] = ind2sub(size(phi), S_ii);
N = [peekR(phi,rr,cc) peekL(phi,rr,cc) peekU(phi,rr,cc) peekD(phi,rr,cc)]';
phi(S_ii) = max(N) - 1;
[rr cc] = ind2sub(size(phi), S_oo);
N = [peekR(phi,rr,cc) peekL(phi,rr,cc) peekU(phi,rr,cc) peekD(phi,rr,cc)]';
phi(S_oo) = min(N) + 1;
%- 9. update bookkeeping
L_z_phi = phi(L_z_);
idx_out = L_z_(L_z_phi_ >= 0 & L_z_phi < 0);
idx_in = L_z_(L_z_phi_ < 0 & L_z_phi >= 0);
h.move_in(idx_in);
h.move_out(idx_out);
end
end
function [phi L_z] = ls_sparse3(phi, L_z, h, iter,img)
% Whitaker "A Level-Set Approach to 3D Reconstruction..." IJCV 1998
CFL = .4;
% shift and peek operations
peekD = @(M,rr,cc,dd) M(safe_sub2ind(size(phi), rr-1, cc, dd));
peekU = @(M,rr,cc,dd) M(safe_sub2ind(size(phi), rr+1, cc, dd));
peekR = @(M,rr,cc,dd) M(safe_sub2ind(size(phi), rr, cc-1, dd));
peekL = @(M,rr,cc,dd) M(safe_sub2ind(size(phi), rr, cc+1, dd));
peekF = @(M,rr,cc,dd) M(safe_sub2ind(size(phi), rr, cc, dd+1));
peekB = @(M,rr,cc,dd) M(safe_sub2ind(size(phi), rr, cc, dd-1));
shiftD = @(M) M([1 1:end-1],:,:);
shiftL = @(M) M(:,[2:end end],:);
shiftR = @(M) M(:,[1 1:end-1],:);
shiftU = @(M) M([2:end end],:,:);
shiftF = @(M) M(:,:,[2:end end]);
shiftB = @(M) M(:,:,[1 1:end-1]);
%- gather layers
Nd = shiftD(phi); Nu = shiftU(phi); Nr = shiftR(phi); Nl = shiftL(phi); Nf = shiftF(phi); Nb = shiftB(phi);
L_i = find(phi < -.5 & (Nu >= -.5 | Nd >= -.5 | Nl >= -.5 | Nr >= -.5 | Nf >= -.5 | Nb >= -.5));
L_o = find(phi > .5 & (Nu <= .5 | Nd <= .5 | Nl <= .5 | Nr <= .5 | Nf <= .5 | Nb <= .5));
%- layer reverse lookup
L_rev = 2*ones(size(phi), 'int8'); % default: outside
L_rev(phi < 0) = -2; % inside
L_rev(L_i) = -1;
L_rev(L_o) = 1;
L_rev(L_z) = 0;
for i = 1:iter
% record sign
L_z_ = L_z;
L_z_phi_ = phi(L_z);
%- 1. update value of zero layer (with uncontrollable force)
F = h.init_iteration(phi, L_z,img);%speed1%speed2
% F = h.init_iteration(phi, L_z);%mean throshhold..
phi(L_z) = phi(L_z) - CFL*F/max(abs(F)); % control that force
%- 3. update status of zero layer
S_i = L_z(phi(L_z) < -.5 );
S_o = L_z( .5 < phi(L_z));
L_z(phi(L_z) < -.5 | .5 < phi(L_z)) = []; % drop
%- 2. update value of nonzero layers
[rr cc dd] = ind2sub(size(phi), L_i);
N = [peekR(phi,rr,cc,dd) peekL(phi,rr,cc,dd) ...
peekU(phi,rr,cc,dd) peekD(phi,rr,cc,dd) ...
peekF(phi,rr,cc,dd) peekB(phi,rr,cc,dd)];
phi(L_i) = max(N,[],2) - 1;
[rr cc dd] = ind2sub(size(phi), L_o);
N = [peekR(phi,rr,cc,dd) peekL(phi,rr,cc,dd) ...
peekU(phi,rr,cc,dd) peekD(phi,rr,cc,dd) ...
peekF(phi,rr,cc,dd) peekB(phi,rr,cc,dd)];
phi(L_o) = min(N,[],2) + 1;
%- 4. update status of nonzero layers
S_zi = L_i( -.5 <= phi(L_i));
S_zo = L_o(phi(L_o) <= .5 );
L_rev(L_i(phi(L_i) < -1.5)) = -2;
L_rev(L_o(phi(L_o) > 1.5)) = 2;
L_i(phi(L_i) < -1.5 | -.5 <= phi(L_i)) = [];
L_o(phi(L_o) <= .5 | 1.5 < phi(L_o)) = [];
%- 5. update reverse lookup
L_rev(S_i) = -1;
L_rev(S_o) = 1;
L_rev([S_zi; S_zo]) = 0;
%- 7. add new neighbors
Nl = shiftL(L_rev) == 0; Nr = shiftR(L_rev) == 0; Nu = shiftU(L_rev) == 0;
Nd = shiftD(L_rev) == 0; Nf = shiftF(L_rev) == 0; Nb = shiftB(L_rev) == 0;
is_near = Nl | Nr | Nu | Nd | Nf | Nb;
S_oo = find(L_rev == 2 & is_near);
S_ii = find(L_rev == -2 & is_near);
L_o = [L_o; S_o; S_oo];
L_i = [L_i; S_i; S_ii];
L_z = [L_z; S_zi; S_zo];
L_rev(S_ii) = -1;
L_rev(S_oo) = 1;
%- 8. update values of new neighbors
[rr cc dd] = ind2sub(size(phi), S_ii);
N = [peekR(phi,rr,cc,dd) peekL(phi,rr,cc,dd) ...
peekU(phi,rr,cc,dd) peekD(phi,rr,cc,dd) ...
peekF(phi,rr,cc,dd) peekB(phi,rr,cc,dd)];
phi(S_ii) = max(N,[],2) - 1;
[rr cc dd] = ind2sub(size(phi), S_oo);
N = [peekR(phi,rr,cc,dd) peekL(phi,rr,cc,dd) ...
peekU(phi,rr,cc,dd) peekD(phi,rr,cc,dd) ...
peekF(phi,rr,cc,dd) peekB(phi,rr,cc,dd)];
phi(S_oo) = min(N,[],2) + 1;
%- 9. update bookkeeping
L_z_phi = phi(L_z_);
idx_out = L_z_(L_z_phi_ > 0 & L_z_phi <= 0);
idx_in = L_z_(L_z_phi_ <= 0 & L_z_phi > 0);
h.move_in(idx_in);
h.move_out(idx_out);
end
end
三、运行结果
四、备注
版本：2014a
