Fast Marching method

Fast Marching method 跟 dijkstra 方法类似,只不过dijkstra方法的路径只能沿网格,而Fast Marching method的方法可以沿斜线.
[Level Set Methods and Fast Marching Methods p94-95

]

这里u理解为到达点的时间, Fijk理解为在点ijk的流速.
然后就可以跟Boundary Value Formulation对应起来了.

[Level Set Methods and Fast Marching Methods p7

]

本例,首先加载一张灰度图片, 将里面像素的值W作为该点的流速F. 然后算到达该点的时间,也就是程序里面的距离D.

matlab部分源代码如下:
example1.m

n = 400;
[M,W] = load_potential_map('road2', n);
%start_point = [16;219];
%end_point = [394;192];
% You can use instead the function
[start_point,end_point] = pick_start_end_point(W);

clear options;
options.nb_iter_max = Inf;
options.end_points = end_point; % stop propagation when end point reached
[D,S] = perform_fast_marching(W, start_point, options);
% nicde color for display
A = convert_distance_color(D);
imageplot({W A}, {'Potential map' 'Distance to starting point'}); 
colormap gray(256);

perform_front_propagation_2d_slow.m

function [D,S,father] = perform_front_propagation_2d_slow(W,start_points,end_points,nb_iter_max,H)

%   [D,S] = perform_front_propagation_2d_slow(W,start_points,end_points,nb_iter_max,H);
%
%   [The mex function is perform_front_propagation_2d]
%
%   'D' is a 2D array containing the value of the distance function to seed.
%   'S' is a 2D array containing the state of each point : 
%       -1 : dead, distance have been computed.
%        0 : open, distance is being computed but not set.
%        1 : far, distance not already computed.
%   'W' is the weight matrix (inverse of the speed).
%   'start_points' is a 2 x num_start_points matrix where k is the number of starting points.
%   'H' is an heuristic (distance that remains to goal). This is a 2D matrix.
%   
%   Copyright (c) 2004 Gabriel Peyr?

data.D = W.*0 + Inf; % action 先把所有点的距离标为Inf
start_ind = sub2ind(size(W), start_points(1,:), start_points(2,:));
data.D( start_ind ) = 0; %将起点的距离设置为0
data.O = start_points; % 将起点加入Open list 
% S=1 : far,  S=0 : open,  S=-1 : close
data.S = ones(size(W));% 将所点的状态设为Far
data.S( start_ind ) = 0; % 将起点的状态设为open(trial)
data.W = W;
data.father = zeros( [size(W),2] );% father维度400*400*2,父节点有两个,因为走斜线

verbose = 1;

if nargin<3
    end_points = [];
end
if nargin<4
    nb_iter_max = round( 1.2*size(W,1)^2 );
end
if nargin<5
    data.H = W.*0;
else
    if isempty(H)
        data.H = W.*0;
    else
        data.H = H;
    end
end

if ~isempty(end_points)
    end_ind = sub2ind(size(W), end_points(1,:), end_points(2,:));
else
    end_ind = [];
end

str = 'Performing Fast Marching algorithm.';
if verbose
    h = waitbar(0,str);
end

i = 0; 
while i0 && jj>0 && ii<=n && jj<=p

        f = [0 0];  % current father

        %%% update the action using Upwind resolution
        P = h/W(ii,jj);
        % neighbors values
        a1 = Inf;
        if ii1 && D( ii-1,jj )1 && D( ii,jj-1 )a2    % swap to reorder
            tmp = a1; a1 = a2; a2 = tmp;
            f = f([2 1]);
        end
        % now the equation is   (a-a1)^2+(a-a2)^2 = P^2, with a >= a2 >= a1.
        % 书上95页公式为:(ux^2 + uy^2)^(1/2)=1/Fijk
        % u理解为到达点的时间,Fijk理解为在点ijk处的流速
        if P^2 > (a2-a1)^2%delta 大于0
            delta = 2*P^2-(a2-a1)^2;
            A1 = (a1+a2+sqrt(delta))/2;
        else%否则用dijkstra方法,沿着格子走,公式为:max|ux,uy|=1/Fijk
            A1 = a1 + P;
            f(2) = 0;%将第2个父节点设为0
        end

        switch S(ii,jj)
            case kClose%闭集不用更新
                % check if action has change. Should not appen for FM
                if A1<d(ii,jj) %="" warning('fastmarching:nonmonotone',="" 'the="" update="" is="" not="" monotone');="" pop="" from="" close="" and="" add="" to="" open="" if="" 0="" don't="" reopen="" points="" o(:,end+1)="[ii;jj];" s(ii,jj)="kOpen;" d(ii,jj)="A1;" end="" case="" kopen%开集才更新="" check="" action="" has="" change.="" a1运行结果: