2D Fast Marching Computations

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

]

这里uu理解为到达点的时间, FijkF_{ijk}理解为在点ijkijk的流速.
然后就可以跟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<3end_points = [];
end
if nargin<4nb_iter_max = round( 1.2*size(W,1)^2 );
end
if nargin<5data.H = W.*0;
elseif isempty(H)data.H = W.*0;elsedata.H = H;end
endif ~isempty(end_points)end_ind = sub2ind(size(W), end_points(1,:), end_points(2,:));
elseend_ind = [];
endstr = 'Performing Fast Marching algorithm.';
if verboseh = waitbar(0,str);
endi = 0;
while i<nb_iter_max && ~isempty(data.O) && isempty( find( data.S(end_ind)==-1 ) )i = i+1;data = perform_fast_marching_step(data);if verbosewaitbar(i/nb_iter_max, h, sprintf('Performing Fast Marching algorithm, iteration %d.', i) );end
endif verboseclose(h);
endD = data.D;
S = data.S;
father = data.father;function data1 = perform_fast_marching_step(data)%有多个起点也是一样的,只需将他们的距离都设为0即可% perform_fast_marching_step - perform one step in the Fast Marching algorithm
%
%   data1 = perform_fast_marching_step(data);
%
%   Data is a structure that records the state before/after a step
%   of the FM algorithm.
%
%   Copyright (c) 2004 Gabriel Peyr?% some constant
kClose = -1;
kOpen = 0;
kFar = 1;D = data.D; % action, a 2D matrix
O = data.O; % open list
S = data.S; % state, either 'O' or 'C', a 2D matrix
H = data.H; % Heuristic
W = data.W; % weight matrix, a 2D array (speed function)
father = data.father;[n,p] = size(D);  % size of the grid,   n,p都为400% step size
h = 1/n;if isempty(O)%看开集是否为空data1 = data;return;
endind_O = sub2ind(size(D), O(1,:), O(2,:));%获取开集里面的顶点的索引[m,I] = min( D(ind_O)+H(ind_O) ); %m里面是最小值,I里面是该最小值在被检测矩阵里面的索引
I = I(1);%取第一个索引
% selected vertex
% 取开集中的第I个点的索引
i = O(1,I);
j = O(2,I);
O(:,I) = [];  % pop from open ,将此点从开集中移除
S(i,j) = kClose; % now its close, 将此点加入闭集(known set)中% its neighbor
% 准备遍历他的右,上,左,下的邻近点
nei = [1,0; 0,1; -1,0; 0,-1 ];for k = 1:4ii = i+nei(k,1);jj = j+nei(k,2);if ii>0 && jj>0 && ii<=n && jj<=pf = [0 0];  % current father%%% update the action using Upwind resolutionP = h/W(ii,jj);% neighbors valuesa1 = Inf;if ii<na1 = D( ii+1,jj );f(1) = sub2ind(size(W), ii+1,jj);endif ii>1 && D( ii-1,jj )<a1a1 = D( ii-1,jj );f(1) = sub2ind(size(W), ii-1,jj);enda2 = Inf;if jj<na2 = D( ii,jj+1 );f(2) = sub2ind(size(W), ii,jj+1);endif jj>1 && D( ii,jj-1 )<a2a2 = D( ii,jj-1 );f(2) = sub2ind(size(W), ii,jj-1);endif a1>a2    % swap to reordertmp = 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处的流速% 那么 ux = (a-a1)/(1/400), uy = (a-a2)/(1/400)% 所以方程变为:((a-a1)^2/(1/400))^2+((a-a2)^2/(1/400))^2 = (1/Wij)^2% 把1/400移到右边,则得Pif P^2 > (a2-a1)^2%delta 大于0delta = 2*P^2-(a2-a1)^2;A1 = (a1+a2+sqrt(delta))/2;else%否则用dijkstra方法,沿着格子走,公式为:max|ux,uy|=1/Fijk% (a-a1) / (1/400) = 1 / WijA1 = a1 + P;f(2) = 0;%将第2个父节点设为0endswitch S(ii,jj)case kClose%闭集不用更新% check if action has change. Should not appen for FMif A1<D(ii,jj)% warning('FastMarching:NonMonotone', 'The update is not monotone');% pop from Close and add to Openif 0        % don't reopen close pointsO(:,end+1) = [ii;jj];S(ii,jj) = kOpen;D(ii,jj) = A1;endendcase kOpen%开集才更新% check if action has change.if A1<D(ii,jj)D(ii,jj) = A1;father(ii,jj,:) = f;endcase kFar%远集不仅更新,而且加入开集if D(ii,jj)~=Infwarning('FastMarching:BadInit', 'Action must be initialized to Inf');  end    % add to openO(:,end+1) = [ii;jj];S(ii,jj) = kOpen;% action must have change.D(ii,jj) = A1;father(ii,jj,:) = f;otherwiseerror('Unknown state');endend
enddata1.D = D;
data1.O = O;
data1.S = S;
data1.W = W;
data1.H = H;
data1.father = father;

运行结果:

matlab完整源代码

2D Fast Marching Computations相关推荐

Fast Marching on 3D Meshes
3D mesh的fast marching 跟2D图片的fast marching类似. 2D图片是规则的平面网格,点的ux,uyu_x, u_y是通过上或下,左或右(具体哪个,是通过距离小的点去确定 ...
fast marching matlab,Fast Marching method
function [D,S,father] = perform_front_propagation_2d_slow(W,start_points,end_points,nb_iter_max,H) % ...
Fast marching on 3D meshes with diffusion distance
与fast marching on 3D meshes with euclidian distance 不同(http://blog.csdn.net/seamanj/article/details/ ...
fast marching method 计算内波相速度
在计算内孤立波传播轨迹的时候,可以应用(Jackson,2009)(J09)提出的经验公式 C(x,y)=Cmaxtanh[B1+H(x,y)B2]−−−−−−−−−−−−−−−−√ C(x,y)=C ...
Fast Marching算法及其在点云测地线计算中的应用
1. 前言: 在解离散测地线问题中,Fast Marching算法被广泛使用.其最大的优点是可以直接应用于点云数据.要知道,大部分内蕴几何算法需要原始数据提供连续的网格信息,才能够被使用.Fast M ...
3D 快速行进算法(Fast Marching Method in 3D case)
Intro 快速行进算法用来高效求解程函方程(Eikonal Equation) F ∥ ∇ U ∥ = 1 F\Vert\nabla U\Vert = 1 F∥∇U∥=1, F为速度场函数( Vel ...
MITK医学Python开发入门详细版
1.关于MITK: MITK的全称是"The Medical Imaging Interaction Toolkit".它是一款开源的交互式医学图像处理软件开发和应用平台.MITK ...
快速行进算法（fast_marching_kroon）的matlab代码
快速行进算法(fast_marching_kroon)的matlab代码快速行进算法用于求解程函方程得到走时场,源代码转载于link 文章目录快速行进算法(fast_marching_kroon) ...
Farthest Point Sampling on 2d image
Farthest Point Sampling的原理是,先随机选一个点,然后呢选择离这个点距离最远的点(D中值最大的点)加入起点,然后继续迭代,直到选出需要的个数为止其主要代码如下: %main.m ...

2D Fast Marching Computations

2D Fast Marching Computations相关推荐

最新文章

热门文章