function outputParams = Image2Mesh(inputParams)

%% NIRFAST PATH

% Define new path
if isfield(inputParams,'nirfastDir') && exist(inputParams.nirfastDir, 'dir')
    NIRFASTPath = inputParams.nirfastDir;
else
    errordlg('Select a NIRFAST Matlab directory', 'I/O Error')
    error('I/O Error: missing input (nirfast matlab path)')
end

% Add path
addpath(genpath(NIRFASTPath));

%% CHECK I/O

if isfield(inputParams,'labelmap')
    inputpath = inputParams.labelmap;
else
    errordlg('Select an input label map', 'I/O Error')
    error('I/O Error: missing input (label map)')
end
if isfield(inputParams,'fiducials')
    fiducials = inputParams.fiducials;
else
    errordlg('Select a fiducials list', 'I/O Error')
    error('I/O Error: missing input (fiducials)')
end
if strcmp(inputParams.meshdir,'/')
    errordlg('Select an output mesh directory', 'I/O Error')
    error('I/O Error: missing output mesh directory)')
end
if strcmp(inputParams.meshdir,'~')
    if ispc
        meshdir = [getenv('HOMEDRIVE') getenv('HOMEPATH')];
    else
        meshdir = getenv('HOME');
    end
else
    meshdir = inputParams.meshdir;
end
if ~isfield(inputParams,'meshname')
    errordlg('Select an output mesh name', 'I/O Error')
    error('I/O Error: missing output mesh name)')
end

%% READ SD

if iscell(fiducials)
% %     executed if more than one fiducials
    sdcoords = cell2mat(fiducials)';
else
% %     executed if only one fiducial
    sdcoords = fiducials';
end

%% READ IMAGE
img = cli_imageread(inputpath);
mask = img.pixelData;

%% READ PARAMETERS FROM FILE
% % File info
% param.Filename = inputpath;
% param.Format = 'NRRD';
% param.DataType = img.metaData.type; % 'short'

% % Compression info
% % TODO : understand those
% param.CompressedData = 'false';
% param.ObjectType = 'image';
% param.BinaryData = 'true';
% param.ByteOrder = 'false';
% param.DataFile = 'LOCAL';
% param.BitDepth = 32;
% param.HeaderSize = 359;

% % Volume Information
% % !!!!!!!!!!!!!!!!!!!!!!!!!!!! WE DO NOT KNOW A PRIORI!!!!!!!!!!!!!!!!!!!!!!!!
% % if the image is oriented differently (as my image is) this will not run!
% % the real orientation can be read from img.metaData.space
% param.AnatomicalOrientation = 'LPS'; % 'LPI'; %
% % noot needed
% param.CenterOfRotation = [0,0,0];
param.NumberOfDimensions = str2double(img.metaData.dimension); % 3
% % Next calculation steps can handle only 3D models. We better check it
% % now.
if param.NumberOfDimensions ~= 3
    error('We do not support 2D meshing yet.')
end
% % not needed and warning: it has 4 numbers not 3! if we have more than one
% % segmented structure. It holds sizes in [structures x y z], e.g. [5 219 239 164].
% dim = sscanf(img.metaData.sizes, '%f')';
% param.Dimensions = dim; % [103,115,130]
% % pffset not needed since we use affine transform now
% % BTW: offset = sscanf(img.metaData.space_origin,'(%f,%f,%f)');
% offset = strrep(img.metaData.space_origin,'(',' ');
% offset = strrep(offset,',',' ');
% offset = strrep(offset,')',' ');
% offset = sscanf(offset, '%f')';
% param.Offset = offset; % [126.9220,2.6726,-78.2673]

% % Transform setup
% % Why not to take the whole affine matrix?
% param.TransformMatrix = img.ijkToLpsTransform(1:3, 1:3);
% % This is the affine transform matrix
param.TransformMatrix = img.ijkToLpsTransform;
% % should be always the same as img.ijkToLpsTransform(1:3,1:3). It's true
% % on all images I have used. So, its better to use transformarion matrix
% % instead of reading a new string value. Do we do not trust the
% % img.ijkToLpsTransformv returned by NIRFASTClicer? We can try to
% % compare both and raise a warning/error just in case?
% % BTW: spacedir = sscanf(img.metaData.space_directions,'(%f,%f,%f) (%f,%f,%f) (%f,%f,%f)')';
% spacedir = strrep(img.metaData.space_directions,'(',' ');
% spacedir = strrep(spacedir,',',' ');
% spacedir = strrep(spacedir,')',' ');
% spacedir = sscanf(spacedir, '%f')';
% spacedir = reshape(spacedir, param.NumberOfDimensions, []);
spacedir = param.TransformMatrix(1:3,1:3);
% % !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
% %  If the anatomical orientation is not 'LPS' as in my case ('right-anterior-superior'),
% %  abs(spacedir(1,1)) can be 0 and this division by zero ruins further
% %  calculations! My 'spacedir' is:
% %      0     0     1
% %     -1     0     0
% %      0     1     0
% % So the following produces NaN
% param.TransformMatrix(1,1) = param.TransformMatrix(1,1) / abs(spacedir(1,1));
% param.TransformMatrix(2,2) = param.TransformMatrix(2,2) / abs(spacedir(2,2));
% param.TransformMatrix(3,3) = param.TransformMatrix(3,3) / abs(spacedir(3,3));
% % This is not needed since wee use affine transform
% param.Offset(1) = param.TransformMatrix(1,1)*param.Offset(1);
% param.Offset(2) = param.TransformMatrix(2,2)*param.Offset(2);
% param.Offset(3) = param.TransformMatrix(3,3)*param.Offset(3);
% % This one i do not understand. Is it some sort of hand-coded fix or to match the fiducials space?
% % Not needed if affine transform is used
% param.TransformMatrix(1,1) = -1 * param.TransformMatrix(1,1);
% param.TransformMatrix(2,2) = -1 * param.TransformMatrix(2,2);

% % In my orientation of ('right-anterior-superior')this will produce 0s.
% % And this does not have to be aligned with xyz, so not only diaonal
% % values.
% % axes
% x=abs(spacedir(1,1));
% y=abs(spacedir(2,2));
% z=abs(spacedir(3,3));
% % This is the way to calculate resolution/voxel_size, i believe.
x = sqrt(sum(spacedir(:,1).^2));
y = sqrt(sum(spacedir(:,2).^2));
z = sqrt(sum(spacedir(:,3).^2));

param.PixelDimensions(1) = x; %(0.73884);
param.PixelDimensions(2) = y; %(0.73884);
param.PixelDimensions(3) = z; %(1)
param.PixelSpacing(1) = x; %(0.73884);
param.PixelSpacing(2) = y; %(0.73884);
param.SliceThickness  = z; %(1);

% Mesh Parameters
param.cell_size = inputParams.cell_size;
param.cell_radius_edge = inputParams.cellradiusedge; % (3.0)
param.facet_size = inputParams.facet_size;
param.facet_angle = inputParams.facetangle; % (25.0)
param.facet_distance = inputParams.facetdistance; % (3.0)
param.special_subdomain_label = (0); % TODO ? label of region to be refined. If this is set to zero, no subdomain resizing will be performed
param.special_subdomain_size  = (0); % TODO ? size of tetrahedron for the special subdomain label
% param.medfilter = 0; % TODO ? apply a median filter to iron out speckle noises of images
% param.pad = 0; % TODO ? Add a zero padding to all XY dimensions

%% WRITE MESHES
% Write CGAL mesh
param.tmppath = meshdir;
% % clean up leftovers
% param.delmedit = 0;
param.delmedit = 1;
% % this formulation is faster
if (exist('RunCGALMeshGenerator','file') > 1)
%     if (exist('RunCGALMeshGenerator'))
% %     see for changes inside RunCGALMeshGenerator
    [e, p] = RunCGALMeshGenerator(mask, param);
    
% %     This is rotation only if param.TransformMatrix = img.ijkToLpsTransform(1:3, 1:3) as in the original script.
%     p= p * param.TransformMatrix;
% %     Originally the translation is carried out inside RunCGALMeshGenerator function
% %     Now this script carries out the the affine transform. All in one
% %     place, close to the affine matrix definition.

% %     The affine transform:
    p = [p ones(size(p,1),1)]';
    p = param.TransformMatrix * p;
    p = p(1:3,:)';
    
% %     It would be nice to group together mesh and fiducials transforms. Here or within
% %     RunCGALMeshGenerator (this needs some redesign of parameters passing).
% %     Further the fiducial coordinates always have the
% %     following rotation matrix from patient to xyz space:
% %     -1     0     0
% %      0    -1     0
% %      0     0     1
% %     no mather what is the orientation of the model
% %     (param.TransformMatrix(1:3,1:3)). So this will fail:
% %     sdcoords = (param.TransformMatrix(1:3,1:3) * sdcoords')';
% %     This works:
    sdcoords = ([-1 0 0;0 -1 0; 0 0 1]*sdcoords')';
% %     However, I see all sort of things that can go wrong with this
% %     hard-coded rotation matrix.
else
    rmpath(genpath(NIRFASTPath))
    errordlg('RunCGALMeshGenerator function does not exist. Check that your path to NIRFAST Matlab is correct.', 'Matlab Error')
    error('Matlab Error: RunCGALMeshGenerator function does not exist. Check that your path to NIRFAST Matlab is correct.')
end

% Optimize mesh
% % not used
% q1 = simpqual(p, e, 'Min_Sin_Dihedral');
if isfield(inputParams,'optimizemesh')
    [genmesh.ele, genmesh.node, mat] = call_improve_mesh_use_stellar(e, p);
% %     not used
%     q2 = simpqual(genmesh.node, genmesh.ele, 'Min_Sin_Dihedral');
else
    genmesh.ele = e;
    genmesh.node = p;
    if size(e,2) > 4
        mat = e(:,5);
    else
        mat = ones(size(e,1),1);
    end
end
genmesh.ele(:,5) = mat;
genmesh.nnpe = 4;
genmesh.dim = 3;

% Write Nirfast mesh
switch inputParams.meshtype
    case 'Standard'
        meshtype = 'stnd';
    case 'Fluorescence'
        meshtype = 'fluor';
    case 'Spectral'
        meshtype = 'spec';
    case 'BEM'
        meshtype = 'stnd_bem';
    case 'BEM Fluorescence'
        meshtype = 'fluor_bem';
    case 'BEM Spectral'
        meshtype = 'spec_bem';
    case 'SPN'
        meshtype = 'stnd_spn';
    otherwise
        error('This type of Nirfast mesh is not supported!');
end

nirfastMeshPath = fullfile(meshdir,inputParams.meshname);
solidmesh2nirfast(genmesh, nirfastMeshPath, meshtype);

%% S/D WINDOW
% % not needed
% mesh = load_mesh(nirfastMeshPath);

h = gui_place_sources_detectors('mesh', nirfastMeshPath);
data = guidata(h);

set(data.sources,  'String', cellstr(num2str(sdcoords,'%.8f %.8f %.8f')));
set(data.detectors, 'String', cellstr(num2str(sdcoords,'%.8f %.8f %.8f')));
plot3(data.mesh, sdcoords(:,1),sdcoords(:,2),sdcoords(:,3),'ro');
plot3(data.mesh, sdcoords(:,1),sdcoords(:,2),sdcoords(:,3),'bx');

%% CLEAN
% Delete files
% TODO ?

%% Help for reading/writing parameters
%
% Reading input parameters
%
%  integer, integer-vector, float, float-vector, double, double-vector, string, string-enumeration, file:
%    value=inputParams.name;
%  string-vector:
%    value=cli_stringvectordecode(inputParams.name);
%  point-vector:
%    value=cli_pointvectordecode(inputParams.name);
%  boolean:
%    value=isfield(inputParams,'name');
%  image:
%    value=cli_imageread(inputParams.name);
%  transform:
%    value=cli_lineartransformread(inputParams.name);
%   or (for generic transforms):
%    value=cli_transformread(inputParams.name);
%  point:
%    values_LPS=cli_pointvectordecode(inputParams.name);
%    [pointDim pointCount]=size(values_LPS);
%  region:
%    values=cli_pointvectordecode(inputParams.name);
%    [regionDim regionCount]=size(values_LPS);
%    center_LPS=[values(1:3,regionIndex); 1];
%    radius_LPS=abs([values(4:6,regionIndex); 1]);
%  measurement:
%    value=cli_measurementread(inputParams.name);
%  geometry:
%    value=cli_geometryread(inputParams.name);
%    Important: in the CLI definition file the following attribute shall be added to the geometry element: fileExtensions=".ply"
%    See https://subversion.assembla.com/svn/slicerrt/trunk/MatlabBridge/src/Examples/MeshScale/ for a complete example.
%
%  Notes:
%    - Input and file (image, transform, measurement, geometry) parameter names are defined by the <longflag> element in the XML file
%    - Output parameter names are defined by the <name> element in the XML file
%    - For retrieving index-th unnamed parameter use inputParams.unnamed{index+1} instead of inputParams.name
%
%
% Writing output parameters
%
%  integer, integer-vector, float, float-vector, double, double-vector, string, string-enumeration, file:
%    outputParams.name=value;
%  image:
%    cli_imagewrite(inputParams.name, value);
%  transform:
%    cli_lineartransformwrite(inputParams.name, value);
%   or (for generic transforms):
%    cli_transformwrite(inputParams.name, value);
%  measurement:
%    cli_measurementwrite(inputParams.name, value);
%  geometry:
%    cli_geometrywrite(inputParams.name, value);
%    Important: in the CLI definition file the following attribute shall be added to the geometry element: fileExtensions=".stl"
%    See https://subversion.assembla.com/svn/slicerrt/trunk/MatlabBridge/src/Examples/MeshScale/ for a complete example.
%
