快速行进算法（Fast Marching）

快速行进算法（Fast Marching）

Fast Marching方法简介

快速行进算法 （Fast Marching Method） 是求解程函方程 （Eikonal Equation）$F|\nabla T|=1$ 的一种高效数值算法，程函方程属于非线性偏微分方程，可以认为是一种近似波动方程。在界面演化问题中的求解的结果 $T$ 的物理含义是曲线以速度 F ( x ⃗ ) F(\vec x) F(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
这里的 Δ + x {\Delta _{ + x}} Δ+x​ 表示对于自变量 $x$ 的一阶向前差分算子， Δ − x {\Delta _{ - x}} Δ−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 T_1 T1​是$x$方向值较小的一个邻居， T 2 T_2 T2​是$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 ≥ T 2 T_1 \ge T_2 T1​≥T2​)：

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

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

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

在第二个条件下，我们还得满足： T 2 + 1 / F ≤ T 1 T_2+1/F \leq T_1 T2​+1/F≤T1​，即 ∣ T 1 − T 2 ∣ ≥ 1 / F |T_1-T_2| \ge 1/F ∣T1​−T2​∣≥1/F，刚好和前面一个条件相补。

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

• ∣ T 1 − T 2 ∣ < 1 / F |T_1-T_2|<1/F ∣T1​−T2​∣<1/F
  t = T 1 + T 2 + 2 / F 2 − ( T 1 − T 2 ) 2 2 t = \frac{T_1+T_2 +\sqrt{2/F^2-(T_1-T_2)^2}}{2} t=2T1​+T2​+2/F2−(T1​−T2​)2 ​​

• ∣ T 1 − T 2 ∣ ≥ 1 / F |T_1-T_2|\ge1/F ∣T1​−T2​∣≥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 i j ≡ 1 F_{ij}\equiv 1 Fij​≡1 。 $S$ 为源点集，在求符号距离函数中，即为零水平离散点集。 T s T_s Ts​ 为源点集对应的函数值，在求符号距离函数中，这个值为零。 $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:dimentiongrid_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;%%让这些源点上的状态值都为0for k=1:source_set_sizecurrent_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) ~= 0points_attribution_status(temp) = 1;%判断是否alive，不是alive置为1T(temp) = h./F(temp) + source_values(k);%当前点周围点的T值初始化.endendenditer = 0;
while 1iter = 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_numcurrent_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:dimentionneighbor_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);endneighbor_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:dimensionTn(u) = min(neighbor_neighbors_T(2*u-1:2*u));endTn = 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))%%%%%%%%%%%%%%%%%%%%%%%addT(updateT_neighbors_ind(t)) = min...(T(updateT_neighbors_ind(t)),T_update);%%%%%%%%%%%%%%%%%%%%%%%%%%endendfunction 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 :dimentionif max(point_set(ix,j) > bounds(j)) || min(point_set(ix,j) < 1)flag (ix) = 1;break;endend
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:dimentionneighbors_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);
endfunction ndx = mysub2ind(siz,point_set)
dimension = size(point_set,2);
switch dimensioncase 1ndx = point_set;%源点位置转全局索引case 2ndx = sub2ind(siz,point_set(:,1),point_set(:,2));%源点位置转全局索引case 3ndx = sub2ind(siz,point_set(:,1),point_set(:,2),point_set(:,3));otherwiseerror('wrong dimension for sub2ind input!');
end
endfunction band_min_T = myind2sub(siz,ndx)
dimension = length(siz);
switch dimensioncase 1band_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);otherwiseerror('wrong dimension for ind2sub output!')
end
endend

本文来自互联网用户投稿，文章观点仅代表作者本人，不代表本站立场，不承担相关法律责任。如若转载，请注明出处。 如若内容造成侵权/违法违规/事实不符，请点击【内容举报】进行投诉反馈！