gtMatchInitialize.m

function handles = gtMatchInitialize(handles)

% Sets initial values for the main GUI's constants, handles, tables
% and figures. Loads difspot data from database.
%
% It should not modify labgeo or any other parameter! All required 
% parameters should be set upstream or from within the GUI.

gtMatchSetBusy(handles, 'on')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Parameters and handles
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Maximum no. of candidates considered (only to preallocate a blank array)
handles.maxcands    = 100;

% Plotting in 3D
handles.plot3D      = true;

% Working directory
handles.workingdir  = pwd;

% Load parameters
parameters = gtLoadParameters();
if (isfield(parameters, 'detgeo'))
    parameters.labgeo = gtGeoConvertDetgeo2LegacyLabgeo(parameters.detgeo, parameters.labgeo);
end

check = inputwdefault('Do you want to reset all the parameters for matching? [y/n]', 'n');
if (strcmpi(check, 'y'))
     parameters.match = gtMatchDefaultParametersGUI();
     parameters.match.thetalimits = 0.5 * ...
         [parameters.labgeo.detanglemin, parameters.labgeo.detanglemax];
     % Do not change parameters file sneakily! (PR)
     % gtSaveParameters(parameters);
     disp('Match parameters have been reset to default values.')
end

if 0
    % checking parameters first: adds missing parameters as new fields
    parameters = gtCheckParameters(parameters, 'acq');
    parameters = gtCheckParameters(parameters, 'match');
    parameters = gtCheckParameters(parameters, 'labgeo');
    % Do not change parameters file sneakily! (PR)
    % gtSaveParameters(parameters);
else
    disp('Skipping gtCheckParameters - very slow (!?)')
end

handles.parfilename = [];

% Store main parameters in handles
handles.parameters = parameters;
handles.matchpars  = rmfield(parameters.match, {'thr_meanerror', 'thetalimits'});
handles.labgeo = parameters.labgeo;
handles.nof_phases = parameters.acq.nof_phases;
handles.nproj      = parameters.acq.nproj;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Matching parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

handles.ErrorFilter = parameters.match.thr_meanerror;
set(handles.MatchPairs_Edit_Filter, 'String', num2str(handles.ErrorFilter))

% Bring theta thresholds to top of list
parnames = fieldnames(handles.matchpars);
perm     = 1:length(parnames);
perm0    = perm;
indA     = find(strcmp('thr_theta', parnames));
indB     = find(strcmp('thr_theta_scale', parnames));

perm(1)    = indA;
perm(2)    = indB;
perm(indA) = perm0(1);
perm(indB) = perm0(2);

handles.matchpars     = orderfields(handles.matchpars, perm);
parnames              = fieldnames(handles.matchpars);
handles.matchparnames = parnames;

% Get parameters description
parlist = build_list_v2();

% Set table fields
partable  = cell(length(parnames), 3);

jj = 1;

for ii = 1:length(parnames)
    partable{jj, 1} = parnames{ii};

    partable{jj, 2} = handles.matchpars.(parnames{ii});

    rowinds        = strcmp(parnames{ii}, parlist.match(:, 1));
    partable{jj, 3} = ['  ', parlist.match{rowinds, 2}];

    jj = jj + 1;
end

partable{strcmp('addconstr', partable(:, 1)), 3} = ' Not used currently';

set(handles.MatchParams_Table, 'RowName', partable(:, 1))
set(handles.MatchParams_Table, 'Data', partable(:, [2, 3]))


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Phases
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for ip = 1:handles.nof_phases
  handles.PhaseNames{ip} = handles.parameters.cryst(ip).name;
  phasetable{ip, 1} = handles.PhaseNames{ip};
  phasetable{ip, 2} = [];
end

handles.UsedPhases = 1:handles.nof_phases;

set(handles.Phases_Listbox_Select, 'String', handles.PhaseNames)
set(handles.Phases_Listbox_Active, 'String', handles.PhaseNames)
set(handles.Phases_Listbox_Active, 'Min', 1);
set(handles.Phases_Listbox_Active, 'Max', length(handles.PhaseNames));


% Theta center line colors
linecolors = ...
    [  0, 200,   0;...
       0,   0, 255;...
     255,   0,   0;...
     200,   0, 200;...
       0, 200, 200;...
     200, 200,   0]/255;

handles.ThetaLineColor = linecolors(mod((0:handles.nof_phases-1), 6)+1, :);


% Theta rectangle colors
rectcolors = ...
    [200, 255, 200;...
     220, 220, 255;...
     255, 220, 220;...
     255, 220, 255;...
     160, 255, 225;...
     225, 235, 100]/255;

handles.ThetaRectColor = rectcolors(mod((0:handles.nof_phases-1), 6)+1, :);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% labgeo.detnorm, labgeo.detorig
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Original detector normal; (re)calculate handles.labgeo.detnorm
origdetnorm = cross(handles.parameters.labgeo.detdiru, ...
                    handles.parameters.labgeo.detdirv);
handles.labgeo.detnorm = origdetnorm/sqrt(origdetnorm*origdetnorm');

% Detector origin: pixel (0, 0) coordinates in LAB ref.
handles.labgeo.detorig = handles.labgeo.detrefpos ...
    - handles.labgeo.detdiru * handles.labgeo.detrefu * handles.labgeo.pixelsizeu ...
    - handles.labgeo.detdirv * handles.labgeo.detrefv * handles.labgeo.pixelsizev;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Labgeobase for fitting
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% U and V direction as a base for fitting
handles.fitting.labgeobase.detdiru = handles.parameters.labgeo.detdiru;
handles.fitting.labgeobase.detdirv = handles.parameters.labgeo.detdirv;

handles.fitting.labgeobase.detnorm = handles.labgeo.detnorm;

% Reference position U, V coordinates
handles.fitting.labgeobase.detrefu = handles.parameters.labgeo.detrefu;
handles.fitting.labgeobase.detrefv = handles.parameters.labgeo.detrefv;

% Beam and rot. axis
handles.fitting.labgeobase.beamdir = handles.parameters.labgeo.beamdir;
handles.fitting.labgeobase.rotpos  = handles.parameters.labgeo.rotpos;
handles.fitting.labgeobase.rotdir  = handles.parameters.labgeo.rotdir;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Detector parameters and table
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Row names in table
rownames{1} = 'Position X (mm)';
rownames{2} = 'Position Y (mm)';
rownames{3} = 'Position Z (mm)';
rownames{4} = 'Tilt in plane (deg)';
rownames{5} = 'Angle U-V (deg)';
rownames{6} = 'Tilt around U (deg)';
rownames{7} = 'Tilt around V (deg)';
rownames{8} = 'Pixel size mean (um)';
rownames{9} = 'Pixel size ratio';

set(handles.DetPars_Table, 'RowName', rownames);


% Initial, saved values
handles.DetParsCorrSaved = zeros(1, 9);
handles.DetParsCorrSaved([1, 2, 3]) = handles.labgeo.detrefpos(1:3);
handles.DetParsCorrSaved([4, 6, 7]) = 0;
handles.DetParsCorrSaved(5) = acosd(handles.labgeo.detdiru * handles.labgeo.detdirv');
handles.DetParsCorrSaved(8) = 1000*0.5*(handles.labgeo.pixelsizeu + handles.labgeo.pixelsizev);
handles.DetParsCorrSaved(9) = handles.labgeo.pixelsizeu / handles.labgeo.pixelsizev;

% Actual and last fitted values
handles.DetParsCorrLast = handles.DetParsCorrSaved;
handles.DetParsCorrAct  = handles.DetParsCorrSaved;

% Default fitted parameters: detector positions and tilts
handles.DetParsTicked = false(1, length(handles.DetParsCorrAct));
handles.DetParsTicked([1 2 4 6 7]) = true;

% Set table for display
dettable = cell(length(handles.DetParsCorrAct), 6);
dettable(:, 1) = num2cell(handles.DetParsCorrSaved);
dettable(:, 2) = num2cell(handles.DetParsTicked);
dettable(:, 4) = num2cell(handles.DetParsCorrAct);
dettable(:, [5, 6]) = {false};
set(handles.DetPars_Table, 'Data', dettable);


% Parameters increment values in table :
% Detector position
handles.DetParsIncrement([1, 2, 3]) = handles.labgeo.detsizeu*handles.labgeo.pixelsizeu*0.01;

% Tilts increment
handles.DetParsIncrement([4, 6, 7]) = 0.5;

% UV angle in degrees
handles.DetParsIncrement(5) = 0.5;

% Pixel sizes
handles.DetParsIncrement(8) = handles.DetParsCorrSaved(8)*0.1;
handles.DetParsIncrement(9) = 0.1;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Lattice parameters and table
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

handles.LattParsTicked  = false(handles.nof_phases, 6);
handles.LattParsToFit   = false(handles.nof_phases, 6);
handles.LatticeParTable = cell(6, 6, handles.nof_phases);
handles.LatticeParTable(:, [2, 5, 6], :) = {false};

for ip = 1:handles.nof_phases
    handles.LattParsCorrSaved(ip, :) = parameters.cryst(ip).latticepar;
    handles.LatticeParTable(:, 1, ip) = num2cell(parameters.cryst(ip).latticepar');
    handles.LatticeParTable(:, 4, ip) = num2cell(parameters.cryst(ip).latticepar');
end

handles.LattParsCorrLast = handles.LattParsCorrSaved;
handles.LattParsCorrAct  = handles.LattParsCorrLast;

% Active phase
ph = handles.UsedPhases(get(handles.Phases_Listbox_Active, 'Value'));

set(handles.LattPars_Table, 'Data', handles.LatticeParTable(:, :, ph));

% Parameters increment values in table
handles.LattParsIncrement(1:3) = 0.01;
handles.LattParsIncrement(4:6) = 0.5;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Tolerances for fitting
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Indices of fitted parameters:
%
% 1  = 'Position X';
% 2  = 'Position Y';
% 3  = 'Position Z';
% 4  = 'Tilt in plane';
% 5  = 'Angle U-V';
% 6  = 'Tilt around U';
% 7  = 'Tilt around V';
% 8  = 'Pixel size mean';
% 9  = 'Pixel size ratio';
% 10 = 'lengtha';
% 12 = 'lengthb';
% 13 = 'lengthc';
% 14 = 'alpha';
% 15 = 'beta';
% 16 = 'gamma';

% Typical values of fitted parameters to scale finite differences
% handles.fitting.typicalx([1, 2, 3])    = 1;     % detector position in mm
% handles.fitting.typicalx([4, 6, 7])    = 0.1;   % detector vector rotations in degrees
% handles.fitting.typicalx(5)          = 0.1;   % detector U-V angle in degrees
% handles.fitting.typicalx(8)          = 1;     % mean pixel size in microns
% handles.fitting.typicalx(9)          = 1;     % pixel size ratio
% handles.fitting.typicalx([10, 11, 12]) = 1;     % unit cell size in Angstrom
% handles.fitting.typicalx([13, 14, 15]) = 10;    % unit cell angles in degrees

% Termination tolerance on fitted parameters
handles.fitting.tolx = 1e-6;

% Termination tolerance on the fitting function value
handles.fitting.tolfun = 1e-10;

% Step size of finite differences to calculate derivatives
handles.fitting.derstep = 1e-7;

% Activate use of equal number of pairs for fitting
if ~isfield(parameters.match, 'equalize_hkls')        
    handles.fitting. equalize_hkls = false;
else
    handles.fitting.equalize_hkls = parameters.match.equalize_hkls;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Define blank plot handles
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

handles.ThetaRectPlot{1} = [];
handles.ThetaRectHist{1} = [];
handles.ThetaLinePlot{1} = [];
handles.ThetaLineHist{1} = [];

handles.pairsfilt.eta    = NaN;
handles.pairsfilt.theta  = NaN;

handles.GUIPairFig.Fullim_RotAxis   = [];
handles.GUIPairFig.Fullim_RotAxis3D = [];


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Database; load data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if isfield(parameters.acq, 'difA_name')
  handles.difspottable = [parameters.acq.difA_name 'difspot'];
else
  handles.difspottable = [parameters.acq.name 'difspot'];
end

handles.pairtable = parameters.acq.pair_tablename ;

handles = gtMatchLoadDifspotData(handles);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Theoretical Bragg angles - load crystallographic data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% X-ray beam wavelength in Angstroms
lambda_true = gtConvEnergyToWavelength(handles.parameters.acq.energy);

% {hkl}-s, and theoretical Bragg angles from energy and crystal structure
for ip = 1:handles.nof_phases

    %handles.CompPhaseThetas{ip} = double(parameters.cryst(ip).theta);
    handles.CompPhaseHKLs{ip}   = double(parameters.cryst(ip).hkl);
    handles.CompIntensities{ip} = double(parameters.cryst(ip).int);

    % 'B matrix': HLK -> Cartesian transformation
    Bmat = gtCrystHKL2CartesianMatrix(handles.LattParsCorrAct(ip, :));
    % d-spacings
    dsp_comp = gtCrystDSpacing(handles.CompPhaseHKLs{ip}, Bmat);
    % Bragg angle theta; avoid using the saved theta values in
    % parameters.cryst, so that the angles surely correspond to the energy
    handles.CompPhaseThetas{ip} = asind(0.5 * lambda_true ./ dsp_comp);

    handles.CompThetaIndices{ip} = 1:length(handles.CompPhaseThetas{ip});
    
    % added by LNervo : keeping trace of used families
    % !!! Notes by PReischig:
    % This part is to be improved. The used families should be stored and read from  
    % parameters.match only. Plus, double check the use of CompThetasTicked and others
    % here when multiple phases have overlapping reflections!
    filename = gtLastFileName(fullfile('3_pairmatching', ['match_' date() '_']), 'existing');
    if (isempty(filename))
        list = dir(fullfile('3_pairmatching', 'match_*.mat'));
        [~, ind] = max([list.datenum]);
        if ~isempty(ind)
            filename = fullfile('3_pairmatching', list(ind).name);
        else
            filename = [];
        end
    end
    if (~isempty(filename))
        match = [];
        load(filename);
        handles.CompThetasTicked{ip} = ismember(handles.parameters.cryst(ip).hkl', match.usedhkls', 'rows')';
    elseif isfield(handles.parameters.cryst(ip), 'usedfam')
        handles.CompThetasTicked{ip} = handles.parameters.cryst(ip).usedfam;
    else
        handles.CompThetasTicked{ip} = true(1, length(handles.CompPhaseThetas{ip}));
    end

    for it = 1:size(handles.CompPhaseHKLs{ip}, 2)
        hkltext = num2str(handles.CompPhaseHKLs{ip}(:, it)');
        hkltext(hkltext == ' ') = [];
        hkltext = ['\{', hkltext, '\}'];
        handles.CompPhaseHKLLabels{ip, it} = hkltext;
    end

end

% Tolerance for distinguishing unique {hkl} families in degrees
if (isfield(parameters.match, 'uniquetheta'))
    handles.ThetaUniqueTol = parameters.match.uniquetheta;
else
    handles.ThetaUniqueTol = 0.1;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Theta filter limits
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if isempty(parameters.match.thetalimits)
    if (isempty(handles.labgeo.detanglemin) ...
            || isempty(handles.labgeo.detanglemax))
        handles.ThetaLimits = [0 45];
    else
        handles.ThetaLimits = 0.5*[handles.labgeo.detanglemin, ...
                                   handles.labgeo.detanglemax];
    end
else
    handles.ThetaLimits = parameters.match.thetalimits;
end

handles.ThetaLimLastLow  = handles.ThetaLimits(1);
handles.ThetaLimLastHigh = handles.ThetaLimits(2);

set(handles.ThetaLims_Listbox, 'String', num2str(handles.ThetaLimits));
set(handles.ThetaLims_Listbox, 'Value', 1);

handles.ThetaFilterPlot = [];


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Interruptibles
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Set every callack non-interruptible
allch = get(handles.figure1, 'Children');
set(allch, 'Interruptible', 'off');
set(allch, 'BusyAction', 'cancel');

% Set every callack non-interruptible, except Pre-match and Match
set(handles.MatchPairs_Push_PreMatch, 'Interruptible', 'on')
set(handles.MatchPairs_Push_Match, 'Interruptible', 'on')
set(handles.MatchPairs_Push_StopCont, 'Interruptible', 'on')


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Progress and status display
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Strings for Pause/Cont button
handles.stopcont = {'Pause', 'Cont.'};

% Plot list
handles.PlotList{1} = 'Eta vs. Theta';
handles.PlotList{2} = 'Statistics';
handles.PlotList{3} = 'Pair figure';
set(handles.Plot_Listbox, 'String', handles.PlotList);

% Define status text and color
handles.statusstring{1} = '!!!';  % warning
handles.statusstring{2} = 'OK';   % OK

handles.statuscolor{1} = [1 0 0]; % red
handles.statuscolor{2} = [0 1 0]; % green


% Set status text and color for status indicators
set(handles.Text_FitDone, 'String', handles.statusstring{1});
set(handles.Text_FitDone, 'BackgroundColor', handles.statuscolor{1});

set(handles.Text_FitSaved, 'String', handles.statusstring{1});
set(handles.Text_FitSaved, 'BackgroundColor', handles.statuscolor{1});

set(handles.Text_MatchDone, 'String', handles.statusstring{1});
set(handles.Text_MatchDone, 'BackgroundColor', handles.statuscolor{1});

set(handles.Text_MatchSaved, 'String', handles.statusstring{1});
set(handles.Text_MatchSaved, 'BackgroundColor', handles.statuscolor{1});


% Display spot results as found in difspot table
spottext{1} = handles.nof_allspots;
spottext{2} = handles.nof_spotsA;
spottext{3} = handles.nof_spotsB;
set(handles.Text_SpotsResults, 'String', spottext);

% Pair results
pairtext{1} = '0000';
pairtext{2} = '0000';
pairtext{3} = '0000';
set(handles.Text_PairsResults, 'String', pairtext);

pairtext{1} = '00.0%';
pairtext{2} = '00.0%';
pairtext{3} = '00.0%';
set(handles.Text_PairsResultsPercent, 'String', pairtext);

set(handles.MatchPairs_Push_StopCont, 'Enable', 'off');

if handles.nof_pairs_ini > 0
    set(handles.MatchPairs_Push_Reset, 'Enable', 'on');
else
    set(handles.MatchPairs_Push_Reset, 'Enable', 'off');
end

handles.matchdone  = false;
handles.matchsaved = false;
handles.fitdone    = false;
handles.fitsaved   = false;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
%%% Save parameters backup file
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5

set(handles.Text_Progress, 'String', '-');
gtMatchSetBusy(handles, 'off')

end