DigitalBreastTomosynthesis/Fitting.m

%% Author: Aleks Prša
% Date: July, 2022
% alex.prsa@gmail.com
% =========================================================================
%{
% 
%     DESCRIPTION:
% 
%     Program to obtain search parameters from projection images using
%     optimisation function and least-square method.
%
%     Faculty of Mathematic and Physics
%     University in Ljubljana
%     Ljubljana, Slovenia
%
%}
% =========================================================================
%% Optimisation function

It = [5.1; 5.5];                            % Anode current of each projection images

Sb = [440.8, 479.2];                        % Signal of dark side of detector
Sigmab = [12.95, 18.31];                    % Standard deviation od a signal of dark side of detector

Sw = [27918, 30200];                        % Signal of bright side of detector
Sigmaw = [335.18, 350.5];                   % Standard deviation od a signal of bright side of detector

% Optimisation function - chi-square method
R = @(x)( ...
         (Sb(1)-x(1)*x(2)*x(5)*It(1)-x(3))^2/1.692^2  + (Sb(2)-x(1)*x(2)*x(5)*It(2)-x(3))^2/3.741^2 + ...
         (Sw(1)-x(1)*x(2)*It(1)-x(3))^2/25^2 + (Sw(2)-x(1)*x(2)*It(2)-x(3))^2/25.823^2 + ...
         (Sigmab(1)^2-x(1)^2*x(2)*x(5)*It(1)-x(4)^2)^2/0.247^2 + (Sigmab(2)^2-x(1)^2*x(2)*x(5)*It(2)-x(4)^2)^2/0.614^2 + ...
         (Sigmaw(1)^2-x(1)^2*x(2)*It(1)-x(4)^2)^2/8.527^2 + (Sigmaw(2)^2-x(1)^2*x(2)*It(2)-x(4)^2)^2/3.1^2 ...
         );

x0 = [4, 1300, 400, 10, 0.0003];            % Initial parameters

[x, r] = fminsearch(R, x0);

Sb1 = x(1)*x(2)*x(5)*It(1) + x(3)
Sb2 = x(1)*x(2)*x(5)*It(2) + x(3)
Sw1 = x(1)*x(2)*It(1) + x(3)
Sw2 = x(1)*x(2)*It(2) + x(3)

Sigmab1 = sqrt(x(1)^2*x(2)*x(5)*It(1) + x(4)^2)
Sigmab2 = sqrt(x(1)^2*x(2)*x(5)*It(2) + x(4)^2)
Sigmaw1 = sqrt(x(1)^2*x(2)*It(1) + x(4)^2)
Sigmaw2 = sqrt(x(1)^2*x(2)*It(2) + x(4)^2)

fprintf('Optimal parameters:\n')
fprintf('\ne0:') 
disp(x(1))
fprintf('\nalpha:') 
disp(x(2))
fprintf('\np:') 
disp(x(3))
fprintf('\nsigmap:') 
disp(x(4))
fprintf('\nepsilon:') 
disp(x(5))
fprintf('\nOstanek:') 
disp(r)