function [spot, grain, conf, tol, spin] = gtINDEXIndividualSpots(grain, spin, tol, parameters)
% GTINDEXINDIVIDUALSPOTS Finds grain ID-s for all diffraction spots and 
% detects conflicts (overlaps and segmentation faults).
%
% gid = gtIndexDifSpots(grain, spin, tol, parameters)
%
% ---------------------------------------------------------
%
% It assigns a grain ID to each diffraction spot. Those for which 
% no acceptable grain is found or compatible with more than one grain 
% (potential overlap) or too close to another spot (spot potentially 
% segmented in multiple parts) will have zero assigned.
%
% A summary of spot and grain conflicts are printed at the end. If a grain
% needs to be merged with another one, just set that grain's "active"
% property to false, and run this function again. The grain will then be
% ignored in the reindexing.
%
%
% Run these functions beforehand:
%   index = load('4_grains/phase_##/index.mat')
%   grain = gtIndexAddInfo(index.grain);
%
% 
% INPUT
%   grain      - all grain data (1xn or nx1 cell array) as in index.mat
%
%   OPTIONAL:
%   spin       - input diffraction spot data; if undefined or empty, 
%                will be loaded from database
%   tol        - tolerances
%                  .u, .v, .w  - maximum absolut deviation between the 
%                                predicted and simulated spot positions 
%                                [pixel, pixel, image]; if undefined, 
%                                they will be computed from the data 
%                                using tol.std;
%                  .hklmaxmult - maximum multiplicity for a given (hkl);
%                                for monochromatic scans and a full 360deg
%                                rotation it's normally 4; {4}
%                  .equigrain  - threshold for grains to be signalled as 
%                                potential equivalents - min. no. of 
%                                overlapping spots between two grains; {3}                              
%                  .std        - max. standard deviations of spot position
%                                and spot properties; {4}
%   parameters - as in parameters file; if undefined or empty, will be
%                loaded from parameters file
%
% OUTPUT
%   spot       - list of diffraction spots with the assigned grain ID 
%                and other info
%   grain      - 'spot' field added with the assigned spot info
%   spin       - the original input diff. spot info (or loaded from database)
%   conf       - lists of potentially conflicting spots and grain combinations, 
%                and the frequency of conflicts (no of spots in conflict 
%                between the grains)
%                  spotoverlap    - true for overlapping spots
%                  allgrains      - all unique grain combinations with overlapping
%                                   spots
%                  allgrainsfreq  - no. of overlapping spots for each 'allgrains'
%                  equigrains     - equivalent grains (no. of overlaps higher than limit)
%                  equigrainsfreq - no. of overlapping spots for each 'equigrains'
%
%
% Version 001 18-04-2014 by P.Reischig
%


if ~exist('tol','var')
    tol.hklmaxmult = [];
end

% Limit to the multiplicity of detection of a given (hkl); 
% for monochromatic scans and a full 360deg rotation, it's normally 4
if ~isfield(tol,'hklmaxmult') || isempty(tol.hklmaxmult)
    tol.hklmaxmult = 4;
end

% No. of overlapping spots as the threshold for equivalent grains  
if ~isfield(tol,'equigrain') || isempty(tol.equigrain)
    tol.equigrain = 3;
end

% Tolerance limits for (u,v,w) spot position deviation from expected; 
% determined in std of the deviations measured for the indexed Friedel pairs
if ~isfield(tol,'std') || isempty(tol.std)
    tol.std = 4;
end

phaseID = grain{1}.phaseid;
fprintf('Processing phase #%d ...\n', phaseID)

%%%%%%%%%%%%%%%%%%
%% Prepare input
%%%%%%%%%%%%%%%%%%

% Load parameters
if ~exist('parameters','var') || isempty(parameters)
    parameters = load('parameters.mat');
    parameters = parameters.parameters;
end

% Get spot and pair data
if ~exist('spin','var') || isempty(spin)
    if isfield(parameters.acq, 'difA_name')
        difspottable = [parameters.acq.difA_name 'difspot'];
    else
        difspottable = [parameters.acq.name 'difspot'];
    end
    
    fprintf('Loading diff. spot data from difspottable %s ...\n', difspottable)
    gtDBConnect;
    mysqlcmd = sprintf([ ...
        'SELECT difspotID, CentroidX, CentroidY, CentroidImage, ' ...
        '   BoundingBoxXsize, BoundingBoxYsize, BoundingBoxXOrigin, ' ...
        '   BoundingBoxYOrigin, Integral ' ...
        'FROM %s ' ...
        'ORDER BY difspotID'], difspottable);
    [spin.id, spin.cu, spin.cv, spin.cw, spin.bbx, spin.bby, spin.bbxo, spin.bbyo, spin.int] = mym(mysqlcmd);
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Vectorise fwd simulated grain spots 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Active grains
if isfield('grain','active')
    actgr = gtFitAllGrainValues(grain,'active',[],1);
	actgr = logical(actgr);
else
    actgr = true(length(grain), 1);
end


gr.cu      = [];
gr.cv      = [];
gr.cw      = [];
gr.ind     = [];
gr.fsimind = [];
%gr.hklspind = [];

for ii = find(actgr)'%1:length(grain)
   gr.cu      = [gr.cu; grain{ii}.fsim.cu'];
   gr.cv      = [gr.cv; grain{ii}.fsim.cv'];
   gr.cw      = [gr.cw; grain{ii}.fsim.cw'];
   %gr.id      = [gr.id; grain{ii}.id(ones(1,length(grain{ii}.fsim.cu)))'];
   gr.ind     = [gr.ind; ii(ones(length(grain{ii}.fsim.cu),1))];
   gr.fsimind = [gr.fsimind; (1:length(grain{ii}.fsim.cu))'];
   %gr.hklspind = [gr.hklspind; grain{ii}.fsim.hklspind'];
end
    
%%%%%%%%%%%%%%%%%%%%%%%
%% Set tolerance limits
%%%%%%%%%%%%%%%%%%%%%%%

% Gather all grain data ; only for statistics
bbusrat = gtFitAllGrainValues(grain(actgr),'stat','bbxsrat',[],[]);
bbvsrat = gtFitAllGrainValues(grain(actgr),'stat','bbysrat',[],[]);
intrat  = gtFitAllGrainValues(grain(actgr),'stat','intrat',[],[]);
uvwA    = gtFitAllGrainValues(grain(actgr),'fsimA','uvw',[],[]);
uvwB    = gtFitAllGrainValues(grain(actgr),'fsimB','uvw',[],[]);
centA   = gtFitAllGrainValues(grain(actgr),'centA',[],[],[]);
centB   = gtFitAllGrainValues(grain(actgr),'centB',[],[],[]);
centimA = gtFitAllGrainValues(grain(actgr),'centimA',[],[],[]);
centimB = gtFitAllGrainValues(grain(actgr),'centimB',[],[],[]);
npairs  = gtFitAllGrainValues(grain(actgr),'nof_pairs',[],[],[]);
nspots  = 2*sum(npairs);

% Deviations between fwd simulation and measured spot positions
du = [uvwA(1,:) - centA(1,:), uvwB(1,:) - centB(1,:)];
dv = [uvwA(2,:) - centA(2,:), uvwB(2,:) - centB(2,:)];
dw = [uvwA(3,:) - centimA,    uvwB(3,:) - centimB];

% Standard deviations
std_du  = std(du(~isnan(du)));
std_dv  = std(dv(~isnan(dv)));
std_dw  = std(dw(~isnan(dw)));
std_bbu = std(bbusrat);
std_bbv = std(bbvsrat);
std_int = std(intrat);
mean_bbu = mean(bbusrat);
mean_bbv = mean(bbvsrat);
mean_int = mean(intrat);

% Set limits within 3std; the abs(mean) value is added to allow for a
% bias (e.g. due to a geometry problem) that might affect just 
% part of the spots (it may be unnecessary but allows for larger errors)
if ~isfield(tol,'u') || isempty(tol.u)
    tol.u = abs(mean(du(~isnan(du)))) + tol.std*std_du;
end
if ~isfield(tol,'v') || isempty(tol.v)
    tol.v = abs(mean(dv(~isnan(dv)))) + tol.std*std_dv;
end
if ~isfield(tol,'w') || isempty(tol.w)
    tol.w = abs(mean(dw(~isnan(dw)))) + tol.std*std_dw;
end

% The bounding box size ratios are always >1. 
lim_bbuhi = mean_bbu + tol.std*std_bbu;
lim_bbvhi = mean_bbv + tol.std*std_bbv;
lim_inthi = mean_int + tol.std*std_int;
lim_bbulo = 1/lim_bbuhi; 
lim_bbvlo = 1/lim_bbvhi; 
lim_intlo = 1/lim_inthi; 


disp(' ')
disp('Tolerance limits used:')
fprintf('  multiplicity of an (hkl) per grain:             %g\n', tol.hklmaxmult)
fprintf('  deviation of spot properties in std:            %g\n', tol.std)
fprintf('  min. no. of shared spots for equivalent grains: %g\n', tol.equigrain)

disp(' ')
disp('Deviation of measured from expected spot positions')
disp('based on Friedel pairs (in std):')
fprintf('  U (pixels):      %g\n', std_du)
fprintf('  V (pixels):      %g\n', std_dv)
fprintf('  W (images):      %g\n', std_dw)

disp(' ')
disp('Bounding box size and intensity ratios in the grains')
disp('based on Friedel pairs:')
disp('  Mean:')
fprintf('    bbox U size: %g\n', mean_bbu)
fprintf('    bbox V size: %g\n', mean_bbv)
fprintf('    intensity:   %g\n', mean_int)
disp('  Std:')
fprintf('    bbox U size: %g\n', std_bbu)
fprintf('    bbox V size: %g\n', std_bbv)
fprintf('    intensity:   %g\n', std_int)

disp(' ')
disp('Tolerance limits used for reindexing individual spots:')
fprintf('  U pos. (pixels):   %g\n', tol.u)
fprintf('  V pos. (pixels):   %g\n', tol.v)
fprintf('  W pos. (images):   %g\n', tol.w)
fprintf('  bbox U size ratio: %g\n', lim_bbuhi)
fprintf('  bbox V size ratio: %g\n', lim_bbvhi)
fprintf('  intensity ratio:   %g\n', lim_inthi)

disp(' ')
disp('All tolerance parameters:')
disp(tol)


%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Loop to find grain ID-s
%%%%%%%%%%%%%%%%%%%%%%%%%%

spin.grainind = zeros(length(spin.cu), 1);  % grain ID-s found
spin.fsimind  = zeros(length(spin.cu), 1);  % linear index in grain.fsim
ok = false(length(gr.cu),3);

graincomb           = cell(0);
conf.spotoverlap    = false(1,length(spin.cu));
conf.equigrains     = [];
conf.equigrainsfreq = [];
cnt                 = 0;


disp(' ')
fprintf('Processing %d spots: 0000000',length(spin.cu))

tic

% Loop through spots - could be a parallel parfor loop
for ii = 1:length(spin.cu)
    
    % Display progress
    cnt = cnt + 1;
    if cnt == 100
       fprintf('\b\b\b\b\b\b\b%07d',ii)
       cnt = 0;
    end
    
    % Check if position is within limits
    ok(:,1) = abs(gr.cu - spin.cu(ii)) < tol.u;
    ok(:,2) = abs(gr.cv - spin.cv(ii)) < tol.v;
    ok(:,3) = abs(gr.cw - spin.cw(ii)) < tol.w;
        
    oks = all(ok, 2);
    sumoks = sum(oks);
    
    % Check metadata, if exactly one candidate grain found
    if sum(oks) == 1
       
        indok = find(oks,1,'first');
        grind = gr.ind(indok);
        
        bbxrat = spin.bbx(ii)/grain{grind}.stat.bbxsmean;
        bbyrat = spin.bby(ii)/grain{grind}.stat.bbysmean;
        intrat = spin.int(ii)/grain{grind}.stat.intmean;
        
        ok2(1) = bbxrat <= lim_bbuhi;
        ok2(2) = bbyrat <= lim_bbvhi;
        ok2(3) = intrat <= lim_inthi;
        ok2(4) = bbxrat >= lim_bbulo;
        ok2(5) = bbyrat >= lim_bbvlo;
        ok2(6) = intrat >= lim_intlo;
        
        if all(ok2)
            % !!! adding results to the spin output is not good; consider
            % using another variable name 
            spin.grainind(ii) = grind;
            spin.fsimind(ii)  = gr.fsimind(indok);
            %hklspind(ii) = gr.hklspind(indok);
        end
    elseif sumoks > 1
        % spot in conflict (potential overlap)
        conf.spotoverlap(ii) = true; 
        confgid       = gr.ind(oks)';
        %allgrainid  = [allgrainid, confgid];
        graincomb{end+1} = confgid;
    end         
end

fprintf('\b\b\b\b\b\b\bdone.\n')
                            
toc


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Conflicts between grains (overlaps)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Find grains in conflict

% Unique conflicting grain combinations
for ii = 1:length(graincomb)
    gcombarray(ii,1:length(graincomb{ii})) = graincomb{ii};
end
unigcarray = unique(gcombarray, 'rows');

% Frequency of conflicts
for ii = 1:size(unigcarray, 1)
    allgrains = unigcarray(ii,:);
    conf.allgrains{ii} = allgrains(allgrains > 0);
    
    equi = bsxfun(@eq, unigcarray(ii,:), gcombarray);    
    conf.allgrainsfreq(ii) = sum(all(equi,2));
end
conf.equigrains     = conf.allgrains(conf.allgrainsfreq > tol.equigrain);
conf.equigrainsfreq = conf.allgrainsfreq(conf.allgrainsfreq > tol.equigrain);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Create output, detect broken and overlapping spots
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

disp('Creating output...')

% Cartesian coordinates of plane normals in the unstrained reference crystal
% (output plane normals are normalized)
Bmat  = gtCrystHKL2CartesianMatrix(parameters.cryst(phaseID).latticepar);
plref = gtCrystHKL2Cartesian(parameters.cryst(phaseID).hklsp, Bmat);

% D-spacings in the undeformed reference crystal
dspref = parameters.cryst(phaseID).dspacing;

% Bins for hklsp indices
bins = 1:size(parameters.cryst(phaseID).hklsp,2);

spot.grainind = [];
spot.hklind   = [];
spot.hklspind = [];
spot.omind    = [];
spot.simu     = [];
spot.simv     = [];
spot.simw     = [];
spot.mesu     = [];
spot.mesv     = [];
spot.mesw     = [];
spot.bbxo     = [];
spot.bbyo     = [];
indspot       = 0;
ovspot        = 0;
okspot        = 0;

for ii = 1:length(grain)
    
    % OVERLAPPING OR BROKEN SPOTS

    % Detect overlapping or broken spots that are too close to each other.
    % After this check the multiplicity of a given (hkl) is also 
    % expected to be ok, so it's not dealt with later on.
    
    % Spot indices assigned to grain ii
    spind = find(spin.grainind==ii);
      
    mesu = spin.cu(spind);
    mesv = spin.cv(spind);
    mesw = spin.cw(spind);

    % Check if spots are too close in (U,V,W) -> overlapping or broken; 
    % use matrix form for speed
    ch(:,:,3) = bsxfun(@(a,b) (abs(a-b) <= 2*tol.w), mesw', mesw);
    ch(:,:,2) = bsxfun(@(a,b) (abs(a-b) <= 2*tol.v), mesv', mesv);
    ch(:,:,1) = bsxfun(@(a,b) (abs(a-b) <= 2*tol.u), mesu', mesu);
    ovm       = all(ch,3);
    ovm(1:(size(ch,1)+1):end) = false;
    ov        = any(ovm, 1);
    clear ch
        
    % Remove overlapping/broken spots from the grain's list
    spind(ov) = [];
    
    ovspot = ovspot + sum(ov);       % total number of overlapping spots
    okspot = okspot + length(spind); % total number of non-overlapping spots
    
    
    % GRAIN OUTPUT
    
    grain{ii}.spot.numspots  = length(spind);
    grain{ii}.spot.difspotid = spin.id(spind)';
    grain{ii}.spot.fsimind   = spin.fsimind(spind);

    fsimind = grain{ii}.spot.fsimind;
    
    grain{ii}.spot.hklind     = grain{ii}.fsim.hklind(fsimind);
    grain{ii}.spot.hklspind   = grain{ii}.fsim.hklspind(fsimind);
    grain{ii}.spot.plref      = grain{ii}.fsim.plref(:,fsimind);
    grain{ii}.spot.dspref     = dspref(grain{ii}.spot.hklind);
    %grain{ii}.spot.sinthundef = grain{ii}.fsim.sinth(fsimind);
    grain{ii}.spot.plsam0     = gtFedRod2RotTensor(grain{ii}.R_vector) *...
                                grain{ii}.spot.plref;
    
    % The final spot index containing only non-overlapping spots
    grain{ii}.spot.spotind = (indspot + 1) : (indspot + length(spind));
    indspot  = indspot + length(spind);
    
    grain{ii}.spot.simu = grain{ii}.fsim.cu(fsimind);
    grain{ii}.spot.simv = grain{ii}.fsim.cv(fsimind);
    grain{ii}.spot.simw = grain{ii}.fsim.cw(fsimind);    
    grain{ii}.spot.mesu = spin.cu(spind)';
    grain{ii}.spot.mesv = spin.cv(spind)';
    grain{ii}.spot.mesw = spin.cw(spind)';


    % SPOT OUTPUT

    spot.grainind = [spot.grainind; ii(ones(length(spind),1),1)];
    spot.hklind   = [spot.hklind;   grain{ii}.fsim.hklind(fsimind)'];
    spot.hklspind = [spot.hklspind; grain{ii}.fsim.hklspind(fsimind)'];
    spot.omind    = [spot.omind;    grain{ii}.fsim.omind(fsimind)'];
    spot.simu     = [spot.simu;     grain{ii}.fsim.cu(fsimind)'];
    spot.simv     = [spot.simv;     grain{ii}.fsim.cv(fsimind)'];
    spot.simw     = [spot.simw;     grain{ii}.fsim.cw(fsimind)'];
    spot.mesu     = [spot.mesu;     spin.cu(spind)];
    spot.mesv     = [spot.mesv;     spin.cv(spind)];
    spot.mesw     = [spot.mesw;     spin.cw(spind)];
    spot.bbxo     = [spot.bbxo;     spin.bbxo(spind)];
    spot.bbyo     = [spot.bbyo;     spin.bbyo(spind)];
    
    % Check multiplicity
    mult     = histc(grain{ii}.spot.hklspind, bins);
    multhigh = find(mult > tol.hklmaxmult); 
    
    if ~isempty(multhigh)
        for jj = 1:length(multhigh)
            fprintf('Multiplicity too high for hklsp %d in grain %d\n',...
                multhigh(jj), ii)
        end
    end
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Print summary of results
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

fprintf('\nTotal number of active grains: %d\n', sum(actgr))

fprintf('\nTotal number of diff. spots:\n')
fprintf('  analysed in reindexing:         %d\n', length(spin.grainind))
fprintf('  indexed in input Friedel pairs: %d\n\n', nspots)
fprintf('  reindexed initially:            + %d\n', sum(spin.grainind > 0) + sum(conf.spotoverlap))
fprintf('  fitted to multiple grains:      - %d\n', sum(conf.spotoverlap))
disp('    (potential spot overlap)')
fprintf('  overlapped within the grains:   - %d\n', ovspot)
disp('    (potential segmentation fault)')
fprintf('  indexed in final output:        = %d\n', okspot)

disp(' ')
fprintf('Potentially equivalent grain combinations: %d\n',length(conf.equigrains))
disp('  Grain indices | No. of indexed spots | No. of overlapping spots')
for ii = 1:length(conf.equigrains)
%     fprintf('  %5d -%5d       %3d :%3d          %3d\n', ...
%         conf.equigrains(ii,1), ...
%         conf.equigrains(ii,2), ...
%         grain{conf.equigrains(ii,1)}.spot.numspots, ...
%         grain{conf.equigrains(ii,2)}.spot.numspots, ...
%         conf.equigrainsfreq(ii))
    fprintf('  %s   |   %s   |   %d\n', num2str(conf.equigrains{ii}), ...
        num2str(gtFitAllGrainValues(grain(conf.equigrains{ii}),'spot','numspots',[],[])), ...
        conf.equigrainsfreq(ii))
end
disp(' ')


% FIGURES
norig = 2*gtFitAllGrainValues(grain,'nof_pairs');
nnew  = gtFitAllGrainValues(grain,'spot','numspots');

figure('name','Number of indexed spots per grain')
bar([norig(:), nnew(:)])
xlabel('Grain ID')
title('Total no. of indexed spots per grain')
legend('Indexed in Friedel pairs', 'Re-indexed individually','Location','NorthEast')

figure('name','Correlation of number of indexed spots per grain')
ha = axes;
hold on
im = accumarray([nnew(:), norig(:)]+1, ones(1,length(norig)), [max(nnew), max(norig)]+1);
imagesc(im, 'Parent',ha);
xlabel('No. of spots indexed originally as Friedel pairs')
ylabel('No. of reindexed spots')
title('Correlation of indexed spots per grain')
axis equal
set(ha, 'ydir','normal')
plot([0 max(norig)], [0 max(norig)], 'w-.')


end % of function