fm.m

function [T eFlag] = fm(F, SourcePoints, Dxyz, varargin)
% usage:
% [T eFlag] = fm(F, SourcePoints, Dxyz, opts)
%
%
% Implements the Fast-marching method to solve the Eikonal eqn
% in 2 or 3 dimensions:
% || grad(T) ||^2 = 1/F^2
%
% Note that this is an interface only. It treats the input and
% executes the correct functions accordingly.
%
%
% T				= distance (in time) map
% F				= velocity map
% SourcePoints	= matrix with dimensions (nDims * NumOfSPs)
% Lxyz			= Domain length in x y and z directions
%
% opts			= struct with the following case-sensitive fields:
%	'implementation' = 'Matlab', 'noHeap' or 'C++' 
%						'noHeap' is the old Matlab implementation, that 
%						doesn't use a heap to store narrow band pixels.
%	'useMC'		 = 1 or 0. Should the method use multicore?
%   'order'      = 1 or 2. Order accuracy of finite difference approx.
% opts can also be entered as name/value pairs
%
% The meaning of the error flag 'eFlag' is given by the displayed
% warnings (at the bottom of this file).


% Add path of functions
addpath('./functions');



if (nargin<3)
	error('Must supply at least the first 3 input arguments.');
end

% Set default options
opts.implementation = 'C++';
opts.useMC = 0;
opts.order = 2;
	
% Set input options
opts = parsepropval(opts,varargin{:});

% Check options for errors
opts.implementation = validatestring(opts.implementation,...
	{'Matlab','C++','noHeap'});
if ~any(opts.useMC == [0 1])
	error('''useMC'' must be either 0 or 1');
end
if ~any(opts.order == [1 2])
	error('''order'' must be either 1 or 2');
end

% Convert order to type integer, coz thats what the C++ program expects.
opts.order = int32(opts.order);



% Check dimensionality coherence
nDims = ndims(F);
if((nDims ~= 2) && (nDims ~= 3))
	error('Problem must be either 2- or 3-dimensional');
end
if((nDims ~= size(SourcePoints,1)) || ...
		(length(Dxyz) ~= nDims))
	error('Incompatible dimensionality of inputs');
end


% Get n, Lx, etc...
n = size(F,2);
m = size(F,1);
if (nDims==3)
	o = size(F,3);
end
dx = Dxyz(1);
dy = Dxyz(2);
if (nDims==3)
	dz = Dxyz(3);
end


% Check that F>=0
if any(F < 0)
	error('F must be >= 0');
end


% Check Nxyz and Lxyz conditions
if ((n < 1) || (m < 1) || (exist('o','var') && (o < 1)))
	error('n,m,o must all be > 0');
end
if ((dx <= 0) || (dy <= 0) || (exist('dz','var') && (dz <= 0)))
	error('Dxyz must all be > 0');
end


% Grid spacing
Lx = n*dx;
Ly = m*dy;
if (nDims==3)
	Lz = o*dz;
end


% Check that all input is real
if(~isreal(F) || ~isreal(SourcePoints) || ~isreal([dx dy]) || ...
		(exist('dz','var') && (~isreal(dz))))
	error('All inputs must be real');
end


% Check that SourcePoints are in the domain
if (	any(SourcePoints(:) < 0) || ...
		any(SourcePoints(1,:) > Lx) || ...
		any(SourcePoints(2,:) > Ly) || ...
		((nDims == 3) && any(SourcePoints(3,:) > Lz)))
	error('SourcePoints supplied are out of bounds');
end

% Make SourcePoints coherent with matlab dimension ordering.
% I.e. swap x and y rows (1 and 2). Convert to integer data type
SourcePoints(1:2,:) = flipud(SourcePoints(1:2,:));


% Run Multicore processes
if opts.useMC
	% Add multicore system path
	addpath('../Multicore');
	
	% Adjust nCores according to the number of cores in computer.
	nCores = 2;
	nSPs = size(SourcePoints,2);
	
	% make sure EvalTimeSingle is big enough...otherwise everything
	% will be done by one core. It has to be greater than the time
	% taken to calculate one time-distance map
	MCsettings.maxEvalTimeSingle = 200;
	
	% Should this be adjusted in case mod(nSPs,nCores) != 0 ???
	MCsettings.nrOfEvalsAtOnce = max(1,floor(nSPs/nCores));
	
	% Other settings
	MCsettings.multicoreDir = 'mcDir';
	MCsettings.useWaitbar = false;
	MCsettings.masterIsWorker = true;
	
	% Put correct parameters in the parameter cell
	parameterCell = cell(1,nSPs);
	if (nDims==2)
		switch opts.implementation
			case 'Matlab'
				funhandle = @fm2d;
			case 'C++'
				funhandle = @fm2dc;
			case 'noHeap'
				funhandle = @fm2d_noHeap;
		end
		for iSP=1:nSPs
			parameterCell{iSP} = ...
				{funhandle,F,SourcePoints(:,iSP),dy,dx,opts.order};
		end
	elseif (nDims==3)
		switch opts.implementation
			case 'Matlab'
				funhandle = @fm3d;
			case 'C++'
				funhandle = @fm3dc;
			case 'noHeap'
				error('noHeap has not been implemented in 3d');
		end
		for iSP=1:nSPs
			parameterCell{iSP} = ...
				{funhandle,F,SourcePoints(:,iSP),dy,dx,dz,opts.order};
		end
	end
	
	% Run multicore processes
	resultCells = startmulticoremaster(@encapsulateOutput, parameterCell, MCsettings);
	
	% Split result cells into T and eFlag cells. Create maxEFlag for
	% use with the warning messages given here.
	maxEFlag = 0;
	T = cell(nSPs,1);
	eFlag = cell(nSPs,1);
	for iSP=1:nSPs
		T{iSP} = resultCells{iSP}{1};
		eFlag{iSP} = resultCells{iSP}{2};
		if (eFlag{iSP} > maxEFlag)
			maxEFlag = eFlag{iSP};
		end
	end
	
% Run regular (non-multicore) processes
else
	if (nDims==2)
		switch opts.implementation
			case 'Matlab'
			[T eFlag] = fm2d(F,SourcePoints,dy,dx,opts.order);
			case 'C++'
			[T eFlag] = fm2dc(F,SourcePoints,dy,dx,opts.order);
			case 'noHeap'
				error('noHeap in 2D is outdated.');
% 			T = fm2d_noHeap(F,SourcePoints,dy,dx,opts.order);
		end
	elseif(nDims==3)
		switch opts.implementation
			case 'Matlab'
				[T eFlag] = fm3d(F,SourcePoints,dy,dx,dz,opts.order);
			case 'C++'
				[T eFlag] = fm3dc(F,SourcePoints,dy,dx,dz,opts.order);
			case 'noHeap'
				error('noHeap has not been implemented in 3D.');
		end
	end
	maxEFlag = eFlag;
end


% Display warning according to error flag.
if (maxEFlag==1)
	display([10 'Warning 1: Had to revert to order 1 scheme at least' 10 ...
		'once because the 2nd order scheme resultet in complex' 10 ...
		'values. Note that this does not mean the upwind' 10 ...
		'causality requirement was violated.']);
elseif (maxEFlag==2)
	display([10 'Warning 2: The upwind causality requirement has' 10 ...
		'been violated somehow. The algorithm had to use a' 10 ...
		'brute fix which might lead to large errors.' 10 ...
		'Too extreme speed map values input?']);
end

end