gtIndexStatisticalFit.m

function [grain, remaining] = gtIndexStatisticalFit(grain, restids, tot,...
                                                   strategy_s, cryst)
% GTINDEXSTATISTICALFIT Adds unindexed reflections to grains based on 
% statistics.
% 
% [grain, remaining] = gtIndexStatisticalFit(grain, restids, ...
%                                            tot, strategy_s, cryst)
%
% ------------------------------------------------------------------
%
% Adds unindexed reflections to existing grains based on statistical 
% fitting. Those reflections of which intensity, bounding box sizes, 
% distance and angles do not deviate more than the specified tolerance 
% value (in std) from the average grain value will be added to the grain. 
% If a reflection fits to more than one grain, it will be added to the
% grain it fits best. 
%
% INPUT
%   grain      - all grain data
%   restids    - unindexed reflection id-s
%   tot        - all original reflection data
%   strategy_s - parameters for statistical fit
%                (as in parameters.index.strategy.s)
%   cryst      - crystallography and other data
%
% OUTPUT 
%   grain      - all grain data updated with added reflections
%   remaining  - remaining unindedxed reflection ID-s after the fit
%


fprintf('\n\n--------------------------------------------');
fprintf('\n ADDING REFLECTIONS BASED ON STATISTICS...\n\n');
disp(strategy_s);
disp(' ');
tic;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Check fit and goodness of fit for each refleciton
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Remaining reflection data
srefs      = gtIndexSelectRefs(restids, tot);
nof_srefs  = length(srefs.id);
nof_grains = length(grain);
conslist   = cell(nof_srefs,1);
gnesslist  = cell(nof_srefs,1);

% All possible normalised signed hkl vectors in cartesian CRYSTAL coordinates
shkldirs_cry = gtCrystHKL2Cartesian(vertcat(cryst.shklfam{:})', cryst.Bmat)';

% Extract all the Rodrigues vectors into column vectors
R_vectors = gtIndexAllGrainValues(grain, 'R_vector', [], 1, 1:3);

% Compute all orientation matrices g
all_g = gtMathsRod2OriMat(R_vectors.');

% Express normalised signed hkl vectors in cartesian SAMPLE CS
all_shkldirs = gtVectorCryst2Lab(shkldirs_cry, all_g); 

for ii = 1:nof_grains
    % Store normalised signed hkl vectors in cartesian SAMPLE CS
    grain{ii}.shkldirs = all_shkldirs(:, :, ii); 

	% check consistency of all single reflections with the given grain based
	% on standard deviations
	[cons, angs] = gtIndexCheckAddRefs(srefs, grain{ii}, strategy_s, 'std');

	if ~isempty(cons)
	
		srefscons = gtIndexSelectRefs(cons, srefs);
		
		dcoms = gtMathsPointLinesDists(grain{ii}.center, ...
                                       [srefscons.ca srefscons.dir]);
		
		% compute goodness values for all consistent reflections
		gness = sfFitGoodness(srefscons, grain{ii}, angs, dcoms) ;
		
		for jj = 1:length(cons)
			% add grain id and goodness value to list
			conslist{cons(jj)}  = [conslist{cons(jj)} ii];
			gnesslist{cons(jj)} = [gnesslist{cons(jj)} gness(jj)];
		end
		
	end
		
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Find best fitting grain for each reflection
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

todel = false(nof_srefs,1);

% 'addtograins' will tell for each grain which reflections to add to it
addtograins = cell(nof_grains,1);

for ii = 1:nof_srefs
	
    % if reflection fits to some grain(s)
	if ~isempty(gnesslist{ii})
	
        % best goodness value for reflection
		[~, minind] = min(gnesslist{ii});
        
		addtograins{conslist{ii}(minind)} = ...
            [addtograins{conslist{ii}(minind)}, srefs.id(ii)];
		
        todel(ii) = true;
        
	end
end

remaining  = restids;
remaining(todel) = [];


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Add reflections to grains
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Fastest is to do a grain with all its new reflections simultaneously.

% Grains to delete
gtodel = [];

% !!! ? can be a parfor loop ?
for ii = 1:nof_grains
	if ~isempty(addtograins{ii})
        
		fprintf('Pairs fitted to grain #%d    :   %s\n', ii, ...
                num2str(tot.pairid(addtograins{ii})'));
            
            
		[grain{ii}, excrefids] = gtIndexModifyGrain(grain{ii}, [], ...
                                 addtograins{ii}', tot, strategy_s, cryst);
		
        remaining = [remaining; excrefids];

        if isempty(grain{ii})
            gtodel = [gtodel ii];
            fprintf('Grain #%d was found invalid and deleted.', ii);
        end

	end
end


% Delete bad grains
grain(gtodel) = [];


disp(' ');
toc;


end % of main function


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% SUB-FUNCTIONS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function goodness = sfFitGoodness(refs,grain,angs,dcoms)

  goodness(:,1) = abs(angs/grain.stat.dangstd);
  goodness(:,2) = abs(dcoms/grain.stat.dcomstd);
  goodness(:,3) = abs((refs.int-grain.stat.intmean)/grain.stat.intstd);
  goodness(:,4) = abs((refs.bbxs-grain.stat.bbxsmean)/grain.stat.bbxsstd);
  goodness(:,5) = abs((refs.bbys-grain.stat.bbysmean)/grain.stat.bbysstd);

  % !!! Normalised intensity should be used in future
  %goodness(:,3) = abs((refs.int-grain.stat.intmean)/grain.stat.intstd);

  goodness = sqrt(sum(goodness.*goodness,2));
	
end