gtINDEXIndividualSpots added to repo.

Complete reindexing of all individual diff. spots with previously found grains.

Signed-off-by: preischig <preischig@gmail.com>

zUtil_Indexter/gtINDEXIndividualSpots.m

+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, sp, parameters)
+%
+% ---------------------------------------------------------
+% WARNING! NOT properly tested yet! Watch for bugs!
+%
+% 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 (opt.)   - 
+%                  .hklmaxmult - maximum multiplicity for a given (hkl);
+%                                for monochromatic scans and a full 360deg
+%                                rotation, it's normally 4
+%                  .equigrain  - threshold for grains to be signalled as 
+%                                potential equivalents - min. no. of 
+%                                overlapping spots between two grains                              
+%                  .std        - max. standard deviations of spot properties
+%   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
+
+
+%%%%%%%%%%%%%%%%%%
+%% 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,Integral FROM %s ORDER BY difspotID',difspottable);
+    [spin.id, spin.cu, spin.cv, spin.cw, spin.bbx, spin.bby, spin.int] = mym(mysqlcmd);
+end
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% Vectorise fwd simulated grain spots 
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+% Active grains
+actgr = gtFitAllGrainValues(grain,'active',[],1);
+
+
+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
+if ~isfield(tol,'u') || isempty(tol.u)
+    tol.u = mean(du(~isnan(du))) + tol.std*std_du;
+end
+if ~isfield(tol,'v') || isempty(tol.v)
+    tol.v = mean(dv(~isnan(dv))) + tol.std*std_dv;
+end
+if ~isfield(tol,'w') || isempty(tol.w)
+    tol.w = mean(dw(~isnan(dw))) + tol.std*std_dw;
+end
+
+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.latticepar);
+plref = gtCrystHKL2Cartesian(parameters.cryst.hklsp, Bmat);
+
+% D-spacings in the undeformed reference crystal
+dspref = parameters.cryst.dspacing;
+
+% Bins for hklsp indices
+bins = 1:size(parameters.cryst.hklsp,2);
+
+spot.grainind = [];
+spot.hklind   = [];
+spot.hklspind = [];
+spot.omind    = [];
+spot.simu     = [];
+spot.simv     = [];
+spot.simw     = [];
+spot.mesu     = [];
+spot.mesv     = [];
+spot.mesw     = [];
+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)];
+    
+    
+    % 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('  indexed in input Friedel pairs: %d\n', nspots)
+fprintf('  analysed in reindexing:         %d\n', length(spin.grainind))
+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 index | Indexed spots | 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(' ')
+
+figure('name','Number of indexed spots per grain')
+hold on
+plot(2*gtFitAllGrainValues(grain,'nof_pairs'),'b.')
+plot(gtFitAllGrainValues(grain,'spot','numspots'),'r.')
+xlabel('grain ID')
+ylabel('Total no. of indexed spots per grain')
+title('Total no. of indexed spots per grain')
+legend('Indexed in Friedel pairs', 'Re-indexed individually','Location','NorthEast')
+
+
+end % of function
\ No newline at end of file