From c14fc1e7a12f6c86d14bac4e4ae80010848ac91f Mon Sep 17 00:00:00 2001
From: Peter Reischig <peter.reischig@esrf.fr>
Date: Wed, 12 Sep 2012 10:10:36 +0000
Subject: [PATCH] Bug fix, small changes and addition of conflicts to the
 output.

git-svn-id: https://svn.code.sf.net/p/dct/code/trunk@731 4c865b51-4357-4376-afb4-474e03ccb993
---
 zUtil_Indexter/gtINDEXMatchGrains.m | 132 +++++++++++++++-------------
 1 file changed, 73 insertions(+), 59 deletions(-)

diff --git a/zUtil_Indexter/gtINDEXMatchGrains.m b/zUtil_Indexter/gtINDEXMatchGrains.m
index aef432d2..7188d6ba 100644
--- a/zUtil_Indexter/gtINDEXMatchGrains.m
+++ b/zUtil_Indexter/gtINDEXMatchGrains.m
@@ -1,9 +1,9 @@
-function [match, dev, dcm, drotm, pippo] = gtINDEXMatchGrains(grain1, grain2, ...
-                                    ph1, dc, drot, tol)
+function [match, dev, dcm, drotm, conflicts] = gtINDEXMatchGrains(grain1, grain2, ...
+                                               ph1, dc, drot, tol)
 
 % GTINDEXMATCHGRAINS Matches multiple indexed grains between two datasets.
 %
-% [match, dev, dcm, drotm, pippo] = gtINDEXMatchGrains(grain1, grain2, ...
+% [match, dev, dcm, drotm, conflicts] = gtINDEXMatchGrains(grain1, grain2, ...
 %                                                ph1, dc, drot, tol)
 %
 % For each grain in the reference set 'grain1' it tries to find the 
@@ -19,8 +19,10 @@ function [match, dev, dcm, drotm, pippo] = gtINDEXMatchGrains(grain1, grain2, ..
 % close enough to find correct matching grains. Then the returned 
 % parameters may be used to update the guess.
 % 
-% Loads crstallography info from parameters.mat in the current folder
-% for phase 'ph1'.
+% NOTE !
+% It needs to be launched from the folder belonging to grain1. 
+% Crstallography info is loaded from parameters.mat in the current folder
+% for phase 'ph1'. 
 %
 % ALGORITHM 
 % First the new grain centers and orientations are calculated for grain2
@@ -39,39 +41,48 @@ function [match, dev, dcm, drotm, pippo] = gtINDEXMatchGrains(grain1, grain2, ..
 % operators.
 % 
 %
-% INPUT 
+% OPTIONAL INPUT 
 %   grain1  - indexed grain data of reference set (as output from Indexter)
 %   grain2  - indexed grain data to be matched (as output from Indexter)
-% Optional:
 %   ph1     - phase ID in grain1 (default is 1)
 %   dc      - global translation between datasets from a position given in 
 %             grain2 to the same in grain1 (default is [0 0 0])
 %   drot    - global rotation around the origin of grain2; the rotation 
 %             matrix that transforms a position vector given in grain2 
 %             into grain1 (size 3x3 for column vectors; default is the
-%             identity matrix):   p1 = rotmat*p2 + dc;
-%   tol     - tolerances for finding matching grains (and default values)
-%     tol.dist  (0.01) - distance between center of mass
-%     tol.Rdist (0.01) - distance in Rodrigues space     
-%     tol.bbxs  (2)    - max. bounding box size x ratio
-%     tol.bbys  (2)    - max. bounding box size y ratio      
-%     tol.int   (1e10) - max. intensity ratio      
+%             identity matrix):   p1 = drot*p2 + dc;
+%   tol     - tolerances for finding matching grains {default values}
+%               tol.dist  {0.01} - distance between center of mass
+%               tol.Rdist {0.01} - distance in Rodrigues space     
+%               tol.bbxs  {2}    - max. bounding box size x ratio
+%               tol.bbys  {2}    - max. bounding box size y ratio      
+%               tol.int   {1e10} - max. intensity ratio      
 %
 %
 % OUTPUT
-%   match  - array of indices of matching grains [grain_ind1 grain_ind2]
-%   dev    - vectors of deviations of grain properties
-%   dcm    - updated guess for dc
-%   drotm  - updated guess for drot
+%   match     - array of indices of matching grains [grain_ind1 grain_ind2]
+%   dev       - vectors of deviations of grain properties
+%   dcm       - updated guess for dc
+%   drotm     - updated guess for drot
+%   conflicts - list of conflicts found (same as printed on command line);
+%                 conflicts{:,1}: grains from set1 for which more than one
+%                                 possible matches were found
+%                 conflicts{:,2}: corresponding possible matches from set2
 %
+% USAGE EXAMPLE
+%   Revise tolerances; then run:
+%     [match,dev,dc,drot,conflicts] = gtINDEXMatchGrains(grain1,grain2,1,[],[],tol);
+%   Then run again a few times to iterate:
+%     [match,dev,dc,drot,conflicts] = gtINDEXMatchGrains(grain1,grain2,1,dc,drot,tol);
 %
-% USAGE
-%   First run:
-%     [match,dev,dc,drot] = gtIndexMatchGrains2(grain1,grain2,1);
-%   Then run a few times to iterate:
-%     [match,dev,dc,drot] = gtIndexMatchGrains2(grain1,grain2,1,dc,drot);
+%   For checking unmerged grains in a single dataset (revise the tolerances!):
+%     gtINDEXMatchGrains;
+%     [~,~,~,~,conflicts] = gtINDEXMatchGrains([],[],1,[],[],tol);
 %
 %
+% Version 002 12-09-2012 by PReischig
+%   Added conflicts to output.
+% 
 % Version 001 25-06-2012 by PReischig
 %
 
@@ -107,6 +118,17 @@ if ~exist('cent','var') || isempty(cent)
   cent = false;	
 end
 
+if ~exist('grain1','var') || isempty(grain1)
+    fname = sprintf('4_grains/phase_%02d/index.mat',ph1);
+    fprintf('Loading grain info from %s ...\n',fname)
+    grain1 = load(fname);
+    grain1 = grain1.grain;
+end
+
+if ~exist('grain2','var') || isempty(grain2)
+    grain2 = grain1;
+end
+
 
 if ~exist('tol','var')
   tol.dist  = 0.01;         % 0.01
@@ -198,26 +220,24 @@ if (any(gtIndexAllGrainValues(grain1,'R_onedge',[],1,[])) || ~isempty(drot))
     disp('Loading crystallographic info from parameters.mat...')
     parameters = [];
     load parameters
+    cryst = parameters.cryst(ph1);
     
-    Bmat = gtCrystHKL2CartesianMatrix(parameters.cryst(ph1).latticepar);
-    spacegroup = parameters.cryst(ph1).spacegroup;
+    Bmat = gtCrystHKL2CartesianMatrix(cryst.latticepar);
+    spacegroup = cryst.spacegroup;
     
     fzone_acc = gtCrystRodriguesFundZone(spacegroup,0);
     fzone_ext = gtCrystRodriguesFundZone(spacegroup,Rext);
     
     
-    if (isfield(parameters.cryst,'hklsp') && ~isempty(parameters.cryst.hklsp))
+    if (isfield(cryst,'hklsp') && ~isempty(cryst.hklsp))
         
         % Take the first three plane normals
-        hklsp = parameters.cryst.hklsp(:,1:3);
+        hklsp = cryst.hklsp(:,1:3);
 
         % All signed hkl-s for the given family
-        shkl1 = parameters.cryst.hklsp(:,parameters.cryst(ph1).thetatypesp ==...
-                  parameters.cryst(ph1).thetatypesp(1));
-        shkl2 = parameters.cryst.hklsp(:,parameters.cryst(ph1).thetatypesp ==...
-                  parameters.cryst(ph1).thetatypesp(2));
-        shkl3 = parameters.cryst.hklsp(:,parameters.cryst(ph1).thetatypesp ==...
-                  parameters.cryst(ph1).thetatypesp(3));
+        shkl1 = cryst.hklsp(:,cryst.thetatypesp == cryst.thetatypesp(1));
+        shkl2 = cryst.hklsp(:,cryst.thetatypesp == cryst.thetatypesp(2));
+        shkl3 = cryst.hklsp(:,cryst.thetatypesp == cryst.thetatypesp(3));
     
         % Both hkl and -h-k-l have to be considered:
         shkl1 = [shkl1, -shkl1];
@@ -226,9 +246,9 @@ if (any(gtIndexAllGrainValues(grain1,'R_onedge',[],1,[])) || ~isempty(drot))
   
     else
         % Take the {1,0,0} planes (doesn't matter if reflections really exist)
-        if size(parameters.cryst(ph1).hkl,1) == 3
+        if size(cryst.hkl,1) == 3
             hklsp = [1 0 0; 0 1 0; 0 0 1]';
-        elseif size(parameters.cryst(ph1).hkl,1) == 4
+        elseif size(cryst.hkl,1) == 4
             hklsp = [1 0 -1 0; 0 1 -1 0; 0 0 0 1]';
         else
             error('Field parameters.cryst.hkl is not set correctly.')
@@ -309,10 +329,10 @@ if (any(gtIndexAllGrainValues(grain1,'R_onedge',[],1,[])) || ~isempty(drot))
         if isempty(Rvec)
             disp('WARNING! Rodrigues vector could not be fitted. Possible bug.' )
             gr2.R_vector(ii,:) = NaN;
-            gr2_rot(:,:,ii) = NaN;
+            gr2_rot(:,:,ii)    = NaN;
         else
             gr2.R_vector(ii,:) = Rvec;
-            gr2_rot(:,:,ii) = Rod2g(gr2.R_vector(ii,:));
+            gr2_rot(:,:,ii)    = Rod2g(gr2.R_vector(ii,:));
         end           
 
         % Rodrigues line and distinct line ID (Rid contains (1,2,3))
@@ -375,7 +395,8 @@ end
 
 disp('Matching grains ...')
 
-pippo = zeros(size(nof_grains1),1);
+conflicts = cell(0,2);
+
 
 for ii = 1:nof_grains1
 
@@ -411,38 +432,31 @@ for ii = 1:nof_grains1
 
 			disp(['Multiple matches found for grain #' num2str(ii) ':'])
 			disp(tomatch)
+            
+            conflicts{end+1,1} = ii;
+            conflicts{end,2}   = sort(tomatch)';
+            
 		end
 
 		match(ii,2) = tomatch(1);
 
-        if length(tomatch) > 1
-            pippo(ii) = tomatch(2);
-        else
-            pippo(ii) = NaN;
-        end
-        
-    else
-        pippo(ii) = [];
 	end
 
 end
 
-% removing NaN
-pippo = [find(~isnan(pippo))' pippo(~isnan(pippo))'];
-
-% sorting and unique
-pippo2 =[];
-for iter = 1 : size(pippo,1)
-    pippo2(iter,:) = sort(pippo(iter,:));
-end
 
-pippo = unique(pippo2,'rows');
-pippo_cc =[];
-for i = 1 : size(pippo,1)
-    pippo_cc{i} = pippo(i,:);
+if isempty(conflicts)
+    fprintf('\nNo conflicts have been found with the actual tolerances.\n')
+else
+    if (length([conflicts{:,2}]) == length(unique([conflicts{:,2}])))
+        fprintf('Match conflict list is unique.\n')
+    else
+        fprintf('Match conflict list is NOT unique.\n')
+        fprintf('  Some grain(s) in set2 may be assigned to more than one in set1.\n')
+    end
 end
 
-pippo = pippo_cc;
+    
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% Evalute grain match
-- 
GitLab