function T = msfm2dim(F, SourcePoints,source_set_values,usesecond, usecross)
% This function MSFM2D calculates the shortest distance from a list of
% points to all other pixels in an image, using the
% Multistencil Fast Marching Method (MSFM). This method gives more accurate
% distances by using second order derivatives and cross neighbours.
%
% T=msfm2d(F, SourcePoints, UseSecond, UseCross)
%
% inputs,
% F: The speed image. The speed function must always be larger
% than zero (min value 1e-8), otherwise some regions will
% never be reached because the time will go to infinity.
% SourcePoints : A list of starting points [2 x N] (distance zero)
% UseSecond : Boolean Set to true if not only first but also second
% order derivatives are used (default)
% UseCross : Boolean Set to true if also cross neighbours
% are used (default)
% outputs,
% T : Image with distance from SourcePoints to all pixels
%
% Note:
% Compile the c file "mex msfm2d.c" for cpu-effective registration
%
% Literature : M. Sabry Hassouna et Al. Multistencils Fast Marching
% Methods: A Highly Accurate Solution to the Eikonal Equation on
% Cartesian Domains
%
% Example,
% SourcePoint = [51; 51];
% SpeedImage = ones([101 101]);
% [X Y] = ndgrid(1:101, 1:101);
% T1 = sqrt((X-SourcePoint(1)).^2 + (Y-SourcePoint(2)).^2);
%
% % Run fast marching 1th order, 1th order multi stencil
% % and 2th orde and 2th orde multi stencil
%
% tic; T1_FMM1 = msfm2d(SpeedImage, SourcePoint, false, false); toc;
% tic; T1_MSFM1 = msfm2d(SpeedImage, SourcePoint, false, true); toc;
% tic; T1_FMM2 = msfm2d(SpeedImage, SourcePoint, true, false); toc;
% tic; T1_MSFM2 = msfm2d(SpeedImage, SourcePoint, true, true); toc;
%
% % Show results
% fprintf('\nResults with T1 (Matlab)\n');
% fprintf('Method L1 L2 Linf\n');
% Results = cellfun(@(x)([mean(abs(T1(:)-x(:))) mean((T1(:)-x(:)).^2) max(abs(T1(:)-x(:)))]), {T1_FMM1(:) T1_MSFM1(:) T1_FMM2(:) T1_MSFM2(:)}, 'UniformOutput',false);
% fprintf('FMM1: %9.5f %9.5f %9.5f\n', Results{1}(1), Results{1}(2), Results{1}(3));
% fprintf('MSFM1: %9.5f %9.5f %9.5f\n', Results{2}(1), Results{2}(2), Results{2}(3));
% fprintf('FMM2: %9.5f %9.5f %9.5f\n', Results{3}(1), Results{3}(2), Results{3}(3));
% fprintf('MSFM2: %9.5f %9.5f %9.5f\n', Results{4}(1), Results{4}(2), Results{4}(3));
%
% Example multiple starting points,
% SourcePoint=rand(2,100)*255+1;
% SpeedImage = ones([256 256]);
% tic; T1_MSFM2 = msfm2d(SpeedImage, SourcePoint, true, true); toc;
% figure, imshow(T1_MSFM2,[]); colormap(hot(256));
%
% Function is written by D.Kroon University of Twente (June 2009)
% Distance image, also used to store the index of narrowband pixels
% during marching process
T = zeros(size(F))-1; %生成F规模的-1矩阵
% Augmented Fast Marching (For skeletonize)
%Ed=nargout>1;
% Euclidian distance image
%if(Ed), Y = zeros(size(F)); end
% Pixels which are processed and have a final distance are frozen
Frozen = zeros(size(F));%生成0矩阵
% Free memory to store neighbours of the (segmented) region
neg_free = 100000;%分配内存
neg_pos=0;%位置为零
%if(Ed),
% neg_list = zeros(4,neg_free);
%else
neg_list = zeros(3,neg_free);%neg_list是3行,neg_free列的零矩阵。
%end
% (There are 3 pixel classes:
% - frozen (processed)
% - narrow band (boundary) (in list to check for the next pixel with smallest distance)
% - far (not yet used)
% Neighbours
ne =[-1 0;
1 0;
0 -1;
0 1];%用于查找邻居点的矩阵,顺序为左右下上。
SourcePoints=int32(floor(SourcePoints));%假如源点不是格点,则四舍五入取整。
% set all starting points to distance zero and frozen
for z=1:size(SourcePoints,2)%对每个源点做遍历
% starting point
x= SourcePoints(1,z); y=SourcePoints(2,z);
%提取源点的两个坐标
% Set starting point to frozen and distance to zero
Frozen(x,y)=1; %源点的Fro值置位为1,表示alive。
T(x,y)=0;%T(x,y)=source_set_values(z);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%改了
%源点赋值
end
% Add all neighbours of the starting points to narrow list
for z=1:size(SourcePoints,2)
% starting point
x=SourcePoints(1,z);
y=SourcePoints(2,z);
for k=1:4,
% Location of neighbour
i=x+ne(k,1); j=y+ne(k,2);%取源点的邻居坐标
% Check if current neighbour is not yet frozen and inside the
% picture
%判断邻居是否尚未alive且在区域内
if((i>0)&&(j>0)&&(i<=size(F,1))&&(j<=size(F,2))&&(Frozen(i,j)==0))
Tt=1/max(F(i,j),eps);%Tt=1/max(F(i,j),eps) + T(x,y);%;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%改了
%Ty=1;
% Update distance in neigbour list or add to neigbour list
if(T(i,j)>0)%如果邻居大于零,说明邻居是做过的非源点
if (neg_list(1,T(i,j))>Tt)%如果neg_list的第一行第T值位大于这个位上的1/F
neg_list(1,T(i,j))=Tt;%取小
end
% if(Ed)
% neg_list(4,T(i,j))=min(Ty,neg_list(4,T(i,j)));
% end
else%如果邻居不大于零
neg_pos=neg_pos+1;%neg_pos自加1
% If running out of memory at a new block
%如何内存超了,再分配10万内存。neg_list的前第一行10万位置零
if(neg_pos>neg_free), neg_free = neg_free +100000; neg_list(1,neg_free)=0; end
% if(Ed)
% neg_list(:,neg_pos)=[Tt;i;j;Ty];
% else
neg_list(:,neg_pos)=[Tt;i;j];%记下为1/F,并记录位置
% end
T(i,j)=neg_pos;%T值记为标为1/F的序号
end
end
end
end
%此时T中哪些邻居点按序号递增
%Far的信息在T中,alive的信息在frozen中,trial的信息在neg_list中
Ttest = T;
for k = 1:neg_pos
i = neg_list(2,k);
j = neg_list(3,k);
Ttest(i,j) = 1;
end
T = Ttest;
% Loop through all pixels of the image
for itt=1:numel(F)%返回元素的总个数
% Get the pixel from narrow list (boundary list) with smallest
% distance value and set it to current pixel location
[t,index]=min(neg_list(1,1:neg_pos));%返回neg_list中最小值的位置
if(neg_pos==0), break; end %如果位置等于0,就直接break了。
x=neg_list(2,index); y=neg_list(3,index);%找到最小值对应的位置
Frozen(x,y)=1;%置位为alive
fprintf('min_position:%d,%d\n',x,y)
Ttest(x,y) = neg_list(1,index);
T(x,y)=neg_list(1,index);%T值刷新为此最小值,%%%%%%%%%%%%%%%%%%%%%%%%这里是真正的给值
%以上的操作就是找窄带中最小值点,并将T值刷新,改点置入alive
%T中要么是真值,要么是-1,要么就是在neg_list的值对应的指标。
% if(Ed), Y(x,y)=neg_list(4,index); end
% Remove min value by replacing it with the last value in the array
if(index<neg_pos),%如果没做到最后一个
neg_list(:,index)=neg_list(:,neg_pos);%将neg_list的这一列赋值为最后一位
x2=neg_list(2,index); y2=neg_list(3,index);
% T(x2,y2)=index;%最后一位对应的坐标位置的T置位为最小值的位置
end
neg_pos =neg_pos-1;%个数减1
%哪一个最小值点去掉了,把窄带中最后一位补进它的位置
% Loop through all 4 neighbours of current pixel
for k=1:4,%将最小值附近的值更新加入窄带
% Location of neighbour
i=x+ne(k,1); j=y+ne(k,2);%找零点坐标
fprintf('%d %d',ne(k,:));
fprintf('%d,%d neighbor:%d,%d ',x,y,i,j);
% Check if current neighbour is not yet frozen and inside the
% picture
%如果零点坐标没有超过边界,且不是alive
if((i>0)&&(j>0)&&(i<=size(F,1))&&(j<=size(F,2))&&(Frozen(i,j)==0))
Tt=CalculateDistance(T,F(i,j),size(F),i,j,usesecond,usecross,Frozen);
% if(Ed)
% Ty=CalculateDistance(Y,1,size(F),i,j,usesecond,usecross,Frozen);
% end