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