快速行进算法（Fast Marching）

快速行进算法（Fast Marching）

Fast Marching方法简介

快速行进算法 （Fast Marching Method） 是求解程函方程 （Eikonal Equation）$F|\nabla T|=1$ 的一种高效数值算法，程函方程属于非线性偏微分方程，可以认为是一种近似波动方程。在界面演化问题中的求解的结果 $T$ 的物理含义是曲线以速度 $F(\vec{x})$ 到达计算域每一点所需要消耗的最短时间。寻找 $T$ 在不同高度上的水平集，它总能给出时间演化曲线的一个位置。当 $F=1$ 时，方程解就代表计算域中的距离场，也就是所谓的符号距离函数。

快速行进算法使用了一个特殊的差商代替微商，比如在二维情况，使用的是如下近似：
∣ ∇ T ∣ 2 ≈ ( max ⁡ { − Δ + x T , Δ − x T , 0 } ) 2 + ( max ⁡ { − Δ + y T , Δ − y T , 0 } ) 2 |\nabla T{|^2} \approx {(\max \{ - {\Delta _{ + x}}T,{\Delta _{ - x}}T,0\} )^2} + {(\max \{ - {\Delta _{ + y}}T,{\Delta _{ - y}}T,0\} )^2} ∣∇T∣2≈(max{−Δ+x​T,Δ−x​T,0})2+(max{−Δ+y​T,Δ−y​T,0})2
这里的 ${\Delta }_{+x}$ 表示对于自变量 $x$ 的一阶向前差分算子， ${\Delta }_{-x}$ 表示对于自变量 $x$ 的一阶向后差分算子。

算法推导

同差分思想， $F|\nabla T|=1$的一阶近似为：

∣ ∇ T ∣ 2 ≈ ( max ⁡ { − Δ + x T , Δ − x T , 0 } ) 2 + ( max ⁡ { − Δ + y T , Δ − y T , 0 } ) 2 |\nabla T{|^2} \approx {(\max \{ - {\Delta _{ + x}}T,{\Delta _{ - x}}T,0\} )^2} + {(\max \{ - {\Delta _{ + y}}T,{\Delta _{ - y}}T,0\} )^2} ∣∇T∣2≈(max{−Δ+x​T,Δ−x​T,0})2+(max{−Δ+y​T,Δ−y​T,0})2

假设$T_1$是$x$方向值较小的一个邻居，$T_2$是$y$方向值较小的一个邻居，$t$是当前待求点的值，为了方便，假设网格步长为1（因为算法是基于网格的，若步长不是为1，最后的结果乘以网格步长h即可），那么有：

( max ⁡ { t − T 1 , 0 } ) 2 + ( max ⁡ { t − T 2 , 0 } ) 2 = 1 / F 2 {(\max \{ t-T_1,0\} )^2} + {(\max \{ t-T_2,0\} )^2}= 1/F^2 (max{t−T1​,0})2+(max{t−T2​,0})2=1/F2

我们现在所要做的，就是解这个方程，分几种情况来考虑(不妨假设$T_1\ge T_2$)：

• t > max ⁡ { T 1 , T 2 } = T 1 t>\max\{T_1,T_2\}=T_1 t>max{T1​,T2​}=T1​
  求解 $(t-T_1{)}^2+(t-T_2{)}^2=1/F^2$，可得(若有实根)：
  $$t=\frac{T_1+T_2+\sqrt{2/F^2-(T_1-T_2{)}^2}}{2}$$
  这里舍去了另外一个根是因为另外一个根不满足满足 t > max ⁡ { T 1 , T 2 } t>\max\{T_1,T_2\} t>max{T1​,T2​}。
• $T_1\ge t>T_2$
  求解$(t-T_2{)}^2=1/F^2$，因F为正，可得：
  $$t=T_2+1/F$$

在这样的两种条件下，我们忽略的一个东西，就是在第一种情况下，有可能判别式小于零，导致无解。另外忽略的一个东西就是有可能我们按上面的方式求出来的t不满足前面的关于$t$的范围的条件假设。

也就说，第一个条件下，我们还得满足：
$$\frac{T_1+T_2+\sqrt{2/F^2-(T_1-T_2{)}^2}}{2}>T_1$$
即$|T_1-T_2|<1/F$，显然，这个条件满足了，判别式自然已经大于零了。

在第二个条件下，我们还得满足：$T_2+1/F\le T_1$，即$|T_1-T_2|\ge 1/F$，刚好和前面一个条件相补。

综上，我们去掉开始时的那个“不妨假设”，我们可以得到如下描述：

• $|T_1-T_2|<1/F$
  $$t=\frac{T_1+T_2+\sqrt{2/F^2-(T_1-T_2{)}^2}}{2}$$

• $|T_1-T_2|\ge 1/F$
  t = min ⁡ { T 1 , T 2 } + 1 / F t=\min\{T_1,T_2\}+1/F t=min{T1​,T2​}+1/F

算法描述

这里给出二维情况下的算法描述：

输入的 $F$ 给出了速度场，在求符号距离函数中， $F_{ij}\equiv 1$ 。 $S$ 为源点集，在求符号距离函数中，即为零水平离散点集。 $T_s$ 为源点集对应的函数值，在求符号距离函数中，这个值为零。 $h$ 为离散空间步长。输出的 $T$ 给出了时间标量场，当求符号距离函数时，它所求解得到的就是符号距离函数的离散表示。对于一个 $N\times N$ 的网格， FMM 的总的算术复杂度为 $\mathrm{O}(N\text{log}N)$ ，误差阶为 $\mathrm{O}(\Delta x)$ 。

快速行进算法和 Dijkstra 算法思想相似，不同之处在于 Dijkstra 算法利用节点之间的欧式距离进行更新，而 FMM 算法利用由程函方程化简得到的近似偏微分方程进行更新。快速行进算法有很多特质：首先，它是“一次通过”的算法，减少了不必要的计算；其次，它计算量小，一般只要计算不超过 $N\text{log}N$ 次。

在水平集方法中，无符号距离函数的精度是至关重要的。因此，人们针对符号距离函数的重构和初始化提出了很多方法、算法和离散方式。

程序代码

提供一个自用的程序代码如下：

function T = fast_marching(source_set,source_values,bounds,h)
%bounds = [Ns Ns]. Ns表示一共有几个点. 点的下标从1开始.
%如若[0,1]N等分，那么Ns = N+1
dimention = size(source_set,2);
grid_num = zeros();
for i = 1:dimention
    grid_num(i) = bounds(i);
end
%算全局坐标，算矩阵的时候，网格要统一往上挪一格。
source_set_size = size(source_set,1);%源点的个数
if bound_exceed_judging(source_set,bounds)
    error('source_set out of bounds！');
end   %源点不能超出那个区域,超出报错


T = ones(grid_num)*inf; %网格点上值都标为无穷
points_attribution_status = 2*ones(grid_num); %网格上点的状态都先标为2
%2表示far，1表示trial，0表示alive
F = ones(grid_num);%特殊考虑
source_set_ind = mysub2ind(grid_num,source_set);%源点位置转全局索引
%对于矩阵取二维下标操作是不可行的，所以要转为全局坐标
T(source_set_ind) = source_values;%让这些格点上的T值都为给定的原始值
points_attribution_status(source_set_ind) = 0;%%让这些源点上的状态值都为0



for k=1:source_set_size
    current_source = source_set(k,:);
    current_source_neighbors = neighbors_seeking(current_source);
    current_source_neighbors_ind = mysub2ind (grid_num,current_source_neighbors);
    %把当前点周围点的全局坐标都找出来
    for temp = current_source_neighbors_ind'
        if points_attribution_status(temp) ~= 0
            points_attribution_status(temp) = 1;
            %判断是否alive，不是alive置为1
            T(temp) = h./F(temp) + source_values(k);
            %当前点周围点的T值初始化.
        end
    end
    
end



iter = 0;
while 1
    iter = iter + 1;
    narrow_band_ind = find(points_attribution_status == 1);
    %标出窄带元素
    if isempty(narrow_band_ind)
        break;
    end
    %如果窄带空了，不再继续
    [~,band_min_T_position] = min(T(narrow_band_ind));
    %找到窄带中第一个T最小值的在窄带中的位置
    band_min_T_ind = narrow_band_ind(band_min_T_position);
    %找到窄带中最小T值全局坐标
    points_attribution_status(band_min_T_ind) = 0;
    %送到alive中去
    
    
    
    %对激活元素周围点的T值进行更新，并加入到窄带中
    band_min_T = myind2sub(grid_num,band_min_T_ind);
    %最小T值(坐标形式)
    
    %找出最小T值的邻居,alive邻居就不需要找了
    band_min_T_neighbors = neighbors_seeking(band_min_T);
    band_min_T_neighbors_ind = mysub2ind(grid_num,band_min_T_neighbors);
    updateT_neighbors_position = points_attribution_status...
        (band_min_T_neighbors_ind) ~= 0;
    updateT_neighbors_ind = band_min_T_neighbors_ind(updateT_neighbors_position);
    points_attribution_status(updateT_neighbors_ind) = 1;
    updateT_neighbors = myind2sub(grid_num , updateT_neighbors_ind);
    updateT_neighbors_num = size (updateT_neighbors,1);
    %要更新T值的元素个数的个数
    
    for t = 1:updateT_neighbors_num
        current_updateT_neighbor = updateT_neighbors(t,:);
        F_current=F(mysub2ind(grid_num,current_updateT_neighbor));
        neighbor_matrix =...
            [-1 0 0;
            1 0 0;
            0 -1 0;
            0 1 0;
            0 0 -1;
            0 0 1];
        for jx = 1:dimention
            neighbor_neighbors_aug(:,jx) = current_updateT_neighbor(jx) ...
                + neighbor_matrix(:,jx);
            bound_exceed_position = find(neighbor_neighbors_aug(:,jx) >...
                bounds(jx) | neighbor_neighbors_aug(:,jx) < 1) ;
            neighbor_neighbors_aug(bound_exceed_position ,jx) =...
                2*current_updateT_neighbor(jx) - ...
                neighbor_neighbors_aug(bound_exceed_position ,jx);
        end
        neighbor_neighbors = neighbor_neighbors_aug(1:2*dimension,:);
        neighbor_neighbors_ind = mysub2ind(grid_num,neighbor_neighbors);
        neighbor_neighbors_T = T(neighbor_neighbors_ind);
        Tn = zeros();
        for u = 1:dimension
            Tn(u) = min(neighbor_neighbors_T(2*u-1:2*u));
        end
        Tn = sort(Tn);
        T_update = Tn(1) + h/F_current;
        ups = 2;
        while ups <= dimension && T_update > Tn(ups)
            T_update = (sum(Tn)+sqrt(sum(Tn)^2 +...
                ups*h^2/F_current^2-ups*(Tn*Tn')))/ups;
            ups = ups + 1;
        end
        %                 band_min_T
        %                 myind2sub([9,9],updateT_neighbors_ind(t))%%%%%%%%%%%%%%%%%%%%%%%add
        T(updateT_neighbors_ind(t)) = min...
            (T(updateT_neighbors_ind(t)),T_update);%%%%%%%%%%%%%%%%%%%%%%%%%%
        
        
    end
    
    
    
    
    
    
end








function flag = bound_exceed_judging(point_set,bounds)
dimention = size(point_set,2);
flag = zeros(size(point_set,1),1);
for ix = 1:size(point_set,1)
    for j = 1 :dimention
        if max(point_set(ix,j) > bounds(j)) || min(point_set(ix,j) < 1)
            flag (ix) = 1;
            break;
        end
    end
end
end
%越界判断函数,只要到边界上都判定为越界



function [neighbors,neighbors_num]= neighbors_seeking(point)
dimention = size(point,2);
neighbor_matrix =...
    [-1 0 0;
    1 0 0;
    0 -1 0;
    0 1 0;
    0 0 -1;
    0 0 1];
%neighbors_aug = zeros();
for a = 1:dimention
    neighbors_aug(:,a)= neighbor_matrix(:,a) + point(a);
end
neighbors(1:dimention*2,1:dimention) = neighbors_aug(1:dimention*2,1:dimention);
neighbors = neighbors(bound_exceed_judging(neighbors,bounds) == 0,:);
%超出边界的邻居都不算是邻居。
neighbors_num = size(neighbors,1);
end



function ndx = mysub2ind(siz,point_set)
dimension = size(point_set,2);
switch dimension
    case 1
        ndx = point_set;%源点位置转全局索引
    case 2
        ndx = sub2ind(siz,point_set(:,1),point_set(:,2));%源点位置转全局索引
    case 3
        ndx = sub2ind(siz,point_set(:,1),point_set(:,2),point_set(:,3));
    otherwise
        error('wrong dimension for sub2ind input!');
end
end

function band_min_T = myind2sub(siz,ndx)
dimension = length(siz);
switch dimension
    case 1
        band_min_T = ind2sub(siz,ndx);
    case 2
        [band_min_T(:,1),band_min_T(:,2)] = ind2sub(siz,ndx);
    case 3
        [band_min_T(:,1),band_min_T(:,2),band_min_T(:,3)] = ind2sub(siz,ndx);
    otherwise
        error('wrong dimension for ind2sub output!')
end
end



end