gtINDEXIndividualSpots: modified tolerance computation; other small changes; new figure;

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

zUtil_Indexter/gtINDEXIndividualSpots.m

@@ -2,10 +2,9 @@ function [spot, grain, conf, tol, spin] = gtINDEXIndividualSpots(grain, spin, to
 % GTINDEXINDIVIDUALSPOTS Finds grain ID-s for all diffraction spots and 
 % detects conflicts (overlaps and segmentation faults).
 %
-% gid = gtIndexDifSpots(grain, sp, parameters)
+% gid = gtIndexDifSpots(grain, spin, tol, 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 
@@ -22,28 +21,34 @@ function [spot, grain, conf, tol, spin] = gtINDEXIndividualSpots(grain, spin, to
 %   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:
+%   OPTIONAL:
 %   spin       - input diffraction spot data; if undefined or empty, 
 %                will be loaded from database
 %   tol        - tolerances
-%                  .u (opt.)   - 
+%                  .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
+%                                rotation it's normally 4; {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
+%                                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
+%   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 
@@ -81,6 +86,8 @@ if ~isfield(tol,'std') || isempty(tol.std)
     tol.std = 4;
 end
 
+phaseID = grain{1}.phaseid;
+fprintf('Processing phase #%d ...\n', phaseID)
 
 %%%%%%%%%%%%%%%%%%
 %% Prepare input
@@ -119,8 +126,9 @@ end
 % Active grains
 if isfield('grain','active')
     actgr = gtFitAllGrainValues(grain,'active',[],1);
+	actgr = logical(actgr);
 else
-    actgr = ones(length(grain), 1);
+    actgr = true(length(grain), 1);
 end
 
 
@@ -174,17 +182,20 @@ mean_bbu = mean(bbusrat);
 mean_bbv = mean(bbvsrat);
 mean_int = mean(intrat);
 
-% Set limits within 3std
+% 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 = mean(du(~isnan(du))) + tol.std*std_du;
+    tol.u = abs(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;
+    tol.v = abs(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;
+    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;
@@ -339,14 +350,14 @@ 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);
+Bmat  = gtCrystHKL2CartesianMatrix(parameters.cryst(phaseID).latticepar);
+plref = gtCrystHKL2Cartesian(parameters.cryst(phaseID).hklsp, Bmat);
 
 % D-spacings in the undeformed reference crystal
-dspref = parameters.cryst.dspacing;
+dspref = parameters.cryst(phaseID).dspacing;
 
 % Bins for hklsp indices
-bins = 1:size(parameters.cryst.hklsp,2);
+bins = 1:size(parameters.cryst(phaseID).hklsp,2);
 
 spot.grainind = [];
 spot.hklind   = [];
@@ -459,18 +470,18 @@ end
 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))
+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)
+fprintf('  overlapped within the grains:   - %d\n', ovspot)
 disp('    (potential segmentation fault)')
-fprintf('  indexed in final output:        %d\n', okspot)
+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')
+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), ...
@@ -484,14 +495,28 @@ for ii = 1:length(conf.equigrains)
 end
 disp(' ')
 
+
+% FIGURES
+norig = 2*gtFitAllGrainValues(grain,'nof_pairs');
+nnew  = gtFitAllGrainValues(grain,'spot','numspots');
+
 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')
+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
\ No newline at end of file