Code covered by the BSD License

# Accurate Fast Marching

### Dirk-Jan Kroon (view profile)

23 Jun 2009 (Updated )

Multistencils second order Fast Marching 2D and 3D including rk4 shortest path and skeletonize

```function  EndPoint = e1(StartPoint, GradientVolume, StepSize)
%  E1 is a function which performs one step of the Euler ray tracing
%
%   EndPoint = E1(StartPoint, GradientVolume, StepSize);
%
%    inputs :
%        StartPoint: 2D or 3D location in vectorfield
%        Stepsize : The stepsize
%
%   outputs :
%        EndPoint : The new location (zero if outside image)
%
% Function is written by D.Kroon University of Twente (Oct 2010)

if(numel(StartPoint)==2)
%Linear interpolation of current location
xBas=[0 0 1 1]+floor(StartPoint(1));
yBas=[0 1 0 1]+floor(StartPoint(2));
xCom=StartPoint(1)-floor(StartPoint(1));
yCom=StartPoint(2)-floor(StartPoint(2));

% Linear interpolation percentages.
perc=[(1-xCom) * (1-yCom);
(1-xCom) * yCom    ;
xCom     * (1-yCom);
xCom     * yCom    ];

% Split in Gradient Volumes

xBas=min(max(xBas,1),size(Gx,1));
yBas=min(max(yBas,1),size(Gx,2));

ind=sub2ind(size(Gx),xBas(:),yBas(:));

% Set a step in the direction of the gradient.
check=(EndPoint(1)<1)||(EndPoint(2)<1)||(EndPoint(1)>size(Gx,1))||(EndPoint(2)>size(Gx,2));
if(check), EndPoint=[0;0]; end
else
%Linear interpolation of current location
xBas=[0 0 0 0 1 1 1 1]+floor(StartPoint(1));
yBas=[0 0 1 1 0 0 1 1]+floor(StartPoint(2));
zBas=[0 1 0 1 0 1 0 1]+floor(StartPoint(3));
xCom=StartPoint(1)-floor(StartPoint(1));
yCom=StartPoint(2)-floor(StartPoint(2));
zCom=StartPoint(3)-floor(StartPoint(3));

% Linear interpolation percentages.
perc=[(1-xCom) * (1-yCom) * (1-zCom);
(1-xCom) * (1-yCom) * zCom;
(1-xCom) * yCom     * (1-zCom);
(1-xCom) * yCom     * zCom;
xCom     * (1-yCom) * (1-zCom);
xCom     * (1-yCom) * zCom;
xCom     * yCom     * (1-zCom);
xCom     * yCom     * zCom;];

% Split in Gradient Volumes

xBas=min(max(xBas,1),size(Gx,1));
yBas=min(max(yBas,1),size(Gx,2));
zBas=min(max(zBas,1),size(Gx,3));

ind=sub2ind(size(Gx),xBas(:),yBas(:),zBas(:));

% Set a step in the direction of the gradient.
check=(EndPoint(1)<1)||(EndPoint(2)<1)||(EndPoint(3)<1)||(EndPoint(1)>size(Gx,1))||(EndPoint(2)>size(Gx,2))||(EndPoint(3)>size(Gx,3));
if(check), EndPoint=[0;0;0]; end
end
```