deconvolveIRF.m

% Results = deconvolveIRF(Detected, IRF, Method, varargin)
function Results = deconvolveIRF(Detected, IRF, Method, varargin)
    %%
    if strcmpi(Method, 'BackSubstitution')
        IRFFirstNonZeroIdx = find(IRF, 1);
        ReducedDimension = length(IRF) - IRFFirstNonZeroIdx + 1;
        U = zeros(ReducedDimension);
        for i = 1 : ReducedDimension
            U(i : -1 : 1, i) = IRF(IRFFirstNonZeroIdx : ...
                IRFFirstNonZeroIdx + i - 1);
        end

        b = reshape(...
            Detected(end : -1 : end - length(IRF) + IRFFirstNonZeroIdx), ...
            ReducedDimension, 1);
        x = solveUxEqualsb(U, b);
        Results.Decay = flipud([zeros(IRFFirstNonZeroIdx - 1, 1); x]);
    end
    %%
    if strcmpi(Method, 'Fourier')
        if (nargin == 3)
            Results.Decay = ifft(fft(Detected) ./ ...
                fft(IRF, length(Detected)));
        else
            Results.Decay = ifft(fft(Detected, varargin{1}) ./ ...
                fft(IRF, varargin{1}), length(Detected));
        end
    end
    %%
    if strcmpi(Method, 'deconv')
        FirstNonZeroIdx = find(IRF, 1);
        ZeroPadding = zeros(1, length(IRF(FirstNonZeroIdx : end)) - 1);
        [Results.Decay, ~] = deconv(...
            [Detected(FirstNonZeroIdx : end), ZeroPadding], ...
            IRF(FirstNonZeroIdx : end));
        Results.Decay = [zeros(1, FirstNonZeroIdx - 1), Results.Decay]; 
    end
    %%
    if strcmpi(Method, 'Fitting1') || strcmpi(Method, 'Fitting2')
        X0 = [1, 2, 1, 100, 0];
        if strcmpi(Method, 'Fitting1')
            XLowerBound = [1, 0, 0, 0, 0];
        else
            XLowerBound = [0, 0, 0, 0, 0];
        end
        XUpperBound = [1, 10, 1000, Inf, 100];
        Options = optimoptions('fmincon', 'Algorithm', 'sqp', ...
            'ConstraintTolerance', 1e-10, 'StepTolerance', 1e-10, ...
            'MaxFunctionEvaluations', 1e6);
        [Results.OptimX, ~, ~, Results.Output] = fmincon(...
            @(X) chiSquareInitialAmplitudeWeight(Detected, IRF, varargin{1}, X), ...
            X0, [], [], [], [], XLowerBound, XUpperBound, [], Options);
        Results.Decay = Results.OptimX(4) * (...
            Results.OptimX(1)       * diff(- exp(- varargin{1} / Results.OptimX(2))) * Results.OptimX(2) + ...
            (1 - Results.OptimX(1)) * diff(- exp(- varargin{1} / Results.OptimX(3))) * Results.OptimX(3) ...
            ); % Background noise is not included in decay curves
    end
    %%
    if strcmpi(Method, 'Fitting1S') || strcmpi(Method, 'Fitting2S')
        X0 = [1, 2, 1, sum(Detected), 0];
        if strcmpi(Method, 'Fitting1S')
            XLowerBound = [1, 0, 0, 0.95 * sum(Detected), 0];
        else
            XLowerBound = [0, 0, 0, 0.95 * sum(Detected), 0];
        end
        XUpperBound = [1, 10, 1000, 1.05 * sum(Detected), 100];
        Options = optimoptions('fmincon', 'Algorithm', 'sqp', ...
            'ConstraintTolerance', 1e-10, 'StepTolerance', 1e-10, ...
            'MaxFunctionEvaluations', 1e6);
        [Results.OptimX, ~, ~, Results.Output] = fmincon(...
            @(X) residualNormTotalIntensityWeight(Detected, IRF, varargin{1}, X), ...
            X0, [], [], [], [], XLowerBound, XUpperBound, [], Options);
        Results.Decay = Results.OptimX(4) * (...
            Results.OptimX(1)       * diff(- exp(- varargin{1} / Results.OptimX(2))) + ...
            (1 - Results.OptimX(1)) * diff(- exp(- varargin{1} / Results.OptimX(3))) ...
            ); % Background noise is not included in decay curves
    end
end

function x = solveUxEqualsb(U, b)
    x = U \ b;
end

function ChiSquare = chiSquareInitialAmplitudeWeight(Detected, IRF, ...
    EvalTEdges, X)
    % For Chi-Square stablity
    Threshold = 25;
    
    FullConvolution = conv(IRF, ...
        X(1)       * diff(- exp(- EvalTEdges / X(2))) * X(2) * X(4) + ...
        (1 - X(1)) * diff(- exp(- EvalTEdges / X(3))) * X(3) * X(4), 'full') ...
        + X(5); % X(4) = (A + A'), and X(5) is camera noise
    
    ChiSquare = sum( ...
        (reshape(Detected, 1, []) - ...
        reshape(FullConvolution(1 : length(IRF)), 1, [])) .^ 2 ./ ...
        max([reshape(FullConvolution(1 : length(IRF)), 1, []); ...
        Threshold * ones(1, length(IRF))]) ...
    );
end

function Residual1Norm = residualNormTotalIntensityWeight(Detected, IRF, ...
    EvalTEdges, X)
    FullConvolution = conv(IRF, ...
        X(1)       * diff(- exp(- EvalTEdges / X(2))) * X(4) + ...
        (1 - X(1)) * diff(- exp(- EvalTEdges / X(3))) * X(4), 'full') ...
        + X(5); % X(5) is camera noise
    
    Residual1Norm = norm(reshape(Detected, 1, []) - ...
        reshape(FullConvolution(1 : length(IRF)), 1, []), ...
        1); % 1-norm
end