function [perpADC, sigFrac, parADC] = SMT3(PA, diffC, G, exclusionIndices, shellVec)


% -- SMT-3 -- %

% This function fits a perpendicular apparent diffusion coefficient (ADC),
% signal fraction and parallel ADC to a powder-averaged diffusion MRI signal,
% assuming a single restricting compartment. The restricting compartment is
% assumed to consist of cylindrical micro-geometries e.g. cylinders, axons.


% -- Inputs & Outputs -- %

% PA - Normalised, powder averaged signal. Dimensions: substrates (rows) x shells
% (columns) x repeats. The third dimension "repeats" contains of repeated acquisitions
% for the same substrate and shell.

% G - Scheme file

% diffC - Intrinsic diffusivity in units of m^2/s. Used to provide bounds for the fit only.

% exclusionIndices - In case the scheme file includes non-isotropically
% distributed directions that are not included in the powder average,
% provide a list of which rows in the scheme file these appear on.
% Otherwise, provide an empty vector [].

% shellVec - Specify which acquisition shells to use in the calculation of
% perpendicular ADC. Dimensions: 1 x N, where N is the number of
% shells.

% sigFrac - Signal (or volume) fraction of the restricted compartment.

% perpADC - Fitted perpendicular ADC in units of m^2/s.

% parADC - Parallel diffusivity in units of m^2/s.


% -- Author -- %

% Mariam Andersson mariama@drcmr.dk, 2021


% -- If used, please cite -- %

% Mariam Andersson, Marco Pizzolato, Hans Martin Kjer, Katrine Forum Skodborg, Henrik Lundell, Tim B. Dyrby,
% Does powder averaging remove dispersion bias in diffusion MRI diameter estimates within real 3D axonal architectures?,
% NeuroImage,
% 2021,
% 118718,
% ISSN 1053-8119,
% https://doi.org/10.1016/j.neuroimage.2021.118718.

format long

% Calculate b-values of the different shells
gyro = 42.58e6 *2 * pi;
bCalc = gyro^2 .* G(:,4).^2 .* G(:,6).^2 .* [G(:,5) - (G(:,6)./3)];
bCalc = reshape(bCalc, length(bCalc)/size(PA,2), size(PA,2));
bCalc(exclusionIndices, :) = [];
bCalc = squeeze(bCalc);
bCalc = bCalc(1,:);
bCalc = repmat(bCalc, size(PA,1),1);
bCalc = bCalc(:,shellVec);



for r = 1 : size(PA,3)
    
    for j = 1 : size(PA, 1)
        
        
        
        tmpSig = PA(j,shellVec, r);
        
        % - PARAMETER RESCALING - %
        diffC_s = diffC * 1e10; % Rescale diffusion coefficients by 1e10
        tmp_b = bCalc(j,:) .* 1e-10; % Rescale b-values by 1e-10
        
        
       %options = optimoptions(@lsqnonlin, 'DiffMinChange', 0.001, 'FunctionTolerance',1e-16,'OptimalityTolerance',1e-10, 'StepTolerance',1e-7, 'MaxIterations',10000, 'Display','off', 'FiniteDifferenceType', 'central');
       options = optimoptions(@lsqnonlin, 'MaxIterations',10000, 'Display','off');

        fun = @(params)  (params(2) * (exp(-tmp_b * params(1)) .* sqrt(pi./(4*tmp_b * (params(3) - params(1))))  .* erf(sqrt(abs(tmp_b*(params(3) - params(1))))))) - tmpSig;
        params_0 = [0,0,diffC_s]; %[perpendicular ADC, volume fraction, parallel ADC]
        
        [params] = lsqnonlin(fun,params_0,[0, 0, diffC_s*0.5], [diffC_s*1.5, 1, diffC_s*1.5], options);
        
        perpADC(j,r) = params(1) * 1e-10;
        
        sigFrac(j,r) = params(2);
        
        parADC(j,r) = params(3);
        
    end
end

end

