gtINDEXMatchGrains.m

function [match, dev, dcm, drotm, conflicts, toBeChecked] = gtINDEXMatchGrains(...
                      grain1, grain2, ph1, dc, drot, tol)
% GTINDEXMATCHGRAINS  Matches multiple indexed grains between two datasets; 
%                     can also be used to find unmerged grains.
%
% [match, dev, dcm, drotm, conflicts, toBeChecked] = gtINDEXMatchGrains(grain1, grain2, ph1, dc, drot, tol)
% --------------------------------------------------------------------------------------------------------
%
% For each grain in the reference set 'grain1' it tries to find the 
% corresponding grain in the set 'grain2'. It uses tolerances in 'tol' 
% or default ones to find matching candidates. The best candidate is 
% selected based on the smallest average deviations of its parameters 
% from the reference grain. Uniqueness of the match is not fulfilled, 
% more than one matches may be found for a grain in 'grain2'.
%
% When used with grain1 = grain2, conflicting grains which should be merged 
% can be found in a dataset. Essentially the same criteria are used as in 
% gtIndexter but the results may be slightly different for matches close 
% to the tolerance limits.  
%
% It considers translation and rotation between the two datasets. 
% Optimized values for the translation and rotation can be achieved by
% using the function iteratively. The first guess by the user has to be 
% close enough to find correct matching grains. Then the returned 
% parameters may be used to update the guess.
% 
% 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
% according to the input translation and rotation. The new Rodrigues 
% vectors for grain2 are sought in the fundamental zone that is extended
% to take into account inaccuracies of the measurements. Candidates for 
% each grain1 is found and the best fitting one is selected, based on a 
% figure of merit that is weighted by the actual tolerance values. Median 
% deviations of the grain centers and Rodrigues vectors are calculated, 
% from which updated values of the misfit between the two datasets are 
% given assuming a small misorientation. The grain centers are not
% but only the grain orientations are used to determine the 
% average misorientation of the two datasets.
% Suggestion: to compare orientations, instead of Rodrigues coordinates, 
% it might be easier to use the misorientation matrices and symmetry 
% operators.
% 
%
% OPTIONAL INPUT: 
%   grain1  = indexed grain data of reference set as output from Indexter
%             {if undefined, it is loaded for the current dataset}
%   grain2  = indexed grain data to be matched as output from Indexter
%             {if undefined, it is same as grain1}
%   ph1     = phase ID in grain1 {default is 1}
%   dc      = global translation (of the origin) 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 = drot*p2 + dc;
%   tol     = tolerances for finding matching grains and their default values
%               tol.distf = 0.5   - distance between center of mass;
%                                   mult. factor of smaller bounding box size
%                                     toldist = tol.distf * min(bbxs,bbys)
%               tol.Rdist = 0.02  - absolute 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
%   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
%   toBeChecked = contains the two grain structures with R_vector lists updated
%                 (moved to the fundamental zone)
%
%
% 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);
%
%   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
%


% Attepts were to determine the average misorientation matrix: 
%   1. Calculating one from the median deviation of the Rodrigues vectors.
%      It works well (small rotations), and is very simple, so this is used. 
%   2. Composed a true rotation matrix from small rotations around the 
%      base axes with the symmetric part of the average misorientation 
%      matrix. It also worked well, but more complicated.
%   3. n-th root mean of misorientation matrices multiplied together. 
%      Didn't work.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Initialize
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


if ~exist('ph1','var') || isempty(ph1)
  ph1 = 1;
end

if ~exist('dc','var') || isempty(dc)
  dc = [0 0 0];
end

if ~exist('drot','var') || isempty(drot)
  drot = [];
end

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') || isempty(tol)
  tol.distf = 0.5;          % 0.5
  tol.Rdist = 0.02;         % 0.05
  tol.bbxs  = 2;            % 3 
  tol.bbys  = 2;            % 3   
  tol.int   = 1e10;         % 1e10
end

nof_grains1 = length(grain1);
nof_grains2 = length(grain2);


%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Prepare grain data
%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Sort grains by Z position and only keep one third at the center 
% if cent
% 	cz1  = gtIndexAllGrainValues(grain1,'center',[],3);
% 	ids1 = (1:length(grain1))';
% 
% 	S = sortrows([ids1,cz1],2);
% 
% 	cgrain1 = [];
% 	for ii = floor(nof_grains1/3):nof_grains1*2/3
% 		cgrain1 = [cgrain1,grain1(S(ii,1))];
% 	end
% 
% 	grain1 = cgrain1;
% 	nof_grains1 = length(grain1);
% end

disp('Loading pixelsize and crystallographic info from parameters.mat ...')
parameters = [];
load parameters
cryst   = parameters.cryst(ph1);
pixsize = 0.5*(parameters.labgeo.pixelsizeu + parameters.labgeo.pixelsizev);


match      = NaN(nof_grains1,2);
match(:,1) = (1:nof_grains1)';


% Load relevant info into the gr1 structure (vectors and matrices):
%gr1.id            = gtIndexAllGrainValues(grain1,'id',[],1);
gr1.ind            = (1:nof_grains1)';

gr1.center(:,1)   = gtIndexAllGrainValues(grain1,'center',[],1);
gr1.center(:,2)   = gtIndexAllGrainValues(grain1,'center',[],2);
gr1.center(:,3)   = gtIndexAllGrainValues(grain1,'center',[],3);

gr1.R_vector(:,1) = gtIndexAllGrainValues(grain1,'R_vector',[],1);
gr1.R_vector(:,2) = gtIndexAllGrainValues(grain1,'R_vector',[],2);
gr1.R_vector(:,3) = gtIndexAllGrainValues(grain1,'R_vector',[],3);

gr1.R_onedge(:,1) = gtIndexAllGrainValues(grain1,'R_onedge',[],1);
gr1.bbxs(:,1)     = gtIndexAllGrainValues(grain1,'stat','bbxsmean',1);
gr1.bbys(:,1)     = gtIndexAllGrainValues(grain1,'stat','bbysmean',1);
gr1.int(:,1)      = gtIndexAllGrainValues(grain1,'stat','intmean',1);


% Load relevant info into the gr2 structure (vectors and matrices):
%gr2.id            = gtIndexAllGrainValues(grain2,'id',[],1);
gr2.ind            = (1:nof_grains2)';

gr2.center(:,1)   = gtIndexAllGrainValues(grain2,'center',[],1);
gr2.center(:,2)   = gtIndexAllGrainValues(grain2,'center',[],2);
gr2.center(:,3)   = gtIndexAllGrainValues(grain2,'center',[],3);

gr2.R_vector(:,1) = gtIndexAllGrainValues(grain2,'R_vector',[],1);
gr2.R_vector(:,2) = gtIndexAllGrainValues(grain2,'R_vector',[],2);
gr2.R_vector(:,3) = gtIndexAllGrainValues(grain2,'R_vector',[],3);

gr2.R_onedge(:,1) = gtIndexAllGrainValues(grain2,'R_onedge',[],1);
gr2.bbxs(:,1)     = gtIndexAllGrainValues(grain2,'stat','bbxsmean',1);
gr2.bbys(:,1)     = gtIndexAllGrainValues(grain2,'stat','bbysmean',1);
gr2.int(:,1)      = gtIndexAllGrainValues(grain2,'stat','intmean',1);


gr2_rot = NaN(3,3,nof_grains2);
dev.rot = NaN(3,3,nof_grains1);


%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Reference planes
%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Extension to Rodrigues fundamental zone
% Rext = 0;
Rext = tol.Rdist;

if (any(gtIndexAllGrainValues(grain1,'R_onedge',[],1,[])) || ~isempty(drot))
        
    Bmat = gtCrystHKL2CartesianMatrix(cryst.latticepar);
    spacegroup = cryst.spacegroup;
    
    fzone_acc = gtCrystRodriguesFundZone(spacegroup,0);
    fzone_ext = gtCrystRodriguesFundZone(spacegroup,Rext);
    
    
    if (isfield(cryst,'hklsp') && ~isempty(cryst.hklsp))
        
        % Take the first three plane normals
        hklsp = cryst.hklsp(:,1:3);

        % All signed hkl-s for the given family
        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];
        shkl2 = [shkl2, -shkl2];
        shkl3 = [shkl3, -shkl3];
  
    else
        % Take the {1,0,0} planes (doesn't matter if reflections really exist)
        if size(cryst.hkl,1) == 3
            hklsp = [1 0 0; 0 1 0; 0 0 1]';
        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.')
        end
        
        % All signed hkl-s for the given family (both hkl and -h-k-l)
        shkl1 = gtCrystSignedHKLs(hklsp(:,1)',spacegroup)';
        shkl2 = gtCrystSignedHKLs(hklsp(:,2)',spacegroup)';
        shkl3 = gtCrystSignedHKLs(hklsp(:,3)',spacegroup)';
    end
    
    
    % Plane normals with Cartesian Crystal coordinates
    pl_cry = gtCrystHKL2Cartesian(hklsp,Bmat);  % takes and returns column vectors
    
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Grain2 grain centers
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if isempty(drot)
    rotmat = eye(3);
else
    rotmat = drot;   % for column vectors ! for row vectors, use transpose!
end

% Rotate and translate grain centers
gr2.center = gr2.center*rotmat' + dc(ones(nof_grains2,1),:);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Grain2 Rodrigues coordinates
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculate new Rodrigues lines and vectors inside the fundamental zone 
% for the new, rotated positions

gr2.Rlines = cell(nof_grains2,1);
gr2.Rids = cell(nof_grains2,1);

% All orientation matrices in dataset2
all_g2 = gtMathsRod2RotMat(gr2.R_vector);

if (any(gtIndexAllGrainValues(grain1,'R_onedge',[],1,[])) || ~isempty(drot))
   
    disp('Calculating Rodrigues coordinates ...');

    % Generate a prime number larger than 3 for each grain to be used to 
    % generate distinct reflection id-s later 
    pr = primes(1e6);
    pr(1:2) = [];
    pr(nof_grains2+1:end) = [];

    % Plane normals in the Sample ref. of dataset2
    all_pl_sam2 = gtVectorCryst2Lab(pl_cry.', all_g2);

    for ii = 1:nof_grains2
        % Plane normals rotated to the position in Sample ref. of dataset1
        pl_sam1 = rotmat*all_pl_sam2(:, :, ii).';

        % Lines in Rodrigues space (might be empty)
        Rline1 = gtCrystRodriguesVector(pl_sam1(:,1)',shkl1',Bmat,fzone_ext);
        Rline2 = gtCrystRodriguesVector(pl_sam1(:,2)',shkl2',Bmat,fzone_ext);
        Rline3 = gtCrystRodriguesVector(pl_sam1(:,3)',shkl3',Bmat,fzone_ext);

        Rline = [Rline1; Rline2; Rline3];

        % Vector of line indices to identify which plane they belong to
        Rid = [1*ones(size(Rline1,1),1); 2*ones(size(Rline2,1),1); 3*ones(size(Rline3,1),1)]; 
        
        % Find new Rodrigues vector
        Rvec = gtCrystRodriguesTestCore(Rline, Rid, fzone_acc, fzone_ext, Rext);
        
        if isempty(Rvec)
            disp('WARNING! Rodrigues vector could not be fitted. Possible bug.' )
            gr2.R_vector(ii,:) = NaN;
            gr2_rot(:, :, ii)  = NaN;
        else
            gr2.R_vector(ii,:) = Rvec;
            gr2_rot(:, :, ii)  = all_rot(:, :, ii);
        end           

        % Rodrigues line and distinct line ID (Rid contains (1,2,3))
        gr2.Rlines{ii} = Rline;
        gr2.Rids{ii}   = Rid*pr(ii);

        % For debugging: 
        %   Rvec = gtCrystRodriguesTest(pl_sam1', [shkl1(:,1), shkl2(:,1), shkl3(:,1)]', spacegroup, Bmat, 0, 0, 0);
       
    end
else
    gr2_rot = all_rot;
end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Grain1 Rodrigues coordinates
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculate Rodrigues lines inside the fundamental zone

gr1.Rlines = cell(nof_grains1,1);
gr1.Rids   = cell(nof_grains1,1);

% All orientation matrices in dataset1
all_g1 = gtMathsRod2RotMat(gr1.R_vector);

% Plane normals in the Sample ref. of dataset2
all_pl_sam1 = gtVectorCryst2Lab(pl_cry.', all_g1);

for ii = 1:nof_grains1
  
    if gr1.R_onedge(ii)
        
        % Plane normals in the Sample ref. of dataset1
        pl_sam1 = all_pl_sam1(:, :, ii);
        
        % Lines in Rodrigues space (some might be empty)
        Rline1 = gtCrystRodriguesVector(pl_sam1(1,:), shkl1', Bmat, fzone_ext);
        Rline2 = gtCrystRodriguesVector(pl_sam1(2,:), shkl2', Bmat, fzone_ext);
        Rline3 = gtCrystRodriguesVector(pl_sam1(3,:), shkl3', Bmat, fzone_ext);
        
        gr1.Rlines{ii} = [Rline1; Rline2; Rline3];
        
        % Vector of line indices to identify which plane they belong to
        gr1.Rids{ii} = [1*ones(size(Rline1,1),1); 2*ones(size(Rline2,1),1); 3*ones(size(Rline3,1),1)];
        
        % For debugging: find Rodrigues vector
        Rvec = gtCrystRodriguesTestCore(gr1.Rlines{ii}, gr1.Rids{ii}, fzone_acc, fzone_ext, Rext);
        %   Rvec = gtCrystRodriguesTest(pl_sam1', [shkl1(:,1), shkl2(:,1), shkl3(:,1)]', spacegroup, Bmat, 0, 0, 0);
        if ~isempty(Rvec)
            gr1.Rvec(ii,:) = Rvec;
        end
    end
    
end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Match grains
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

disp('Matching grains ...')

conflicts = cell(0,2);


for ii = 1:nof_grains1

    % Actual grain data
	act = sfSelectGrains(ii,gr1);      
    
    % Get candidate grain indices to be matched with the actual grain
	tomatch = sfCheckMatching(act, gr2, tol, fzone_acc, fzone_ext, pixsize); 

	if ~isempty(tomatch)

		ltm = length(tomatch);
		
		if ltm > 1

            toldist    = min(act.bbxs, act.bbys) * pixsize * tol.distf;

			dcenter    = sum((act.center(ones(ltm,1),:)   - gr2.center(tomatch,:)).^2,2)   / (toldist^2);
			dR_vector  = sum((act.R_vector(ones(ltm,1),:) - gr2.R_vector(tomatch,:)).^2,2) / (tol.Rdist^2);
			dbbxs      = sum((act.bbxs - gr2.bbxs(tomatch,:)).^2,2)/(act.bbxs^2) / (tol.bbxs^2);
			dbbys      = sum((act.bbys - gr2.bbys(tomatch,:)).^2,2)/(act.bbys^2) / (tol.bbys^2);
			dint       = sum((act.int  - gr2.int(tomatch,:)).^2,2) / (tol.int^2);

			dsum = dcenter + dR_vector + dbbxs + dbbys + dint;
			
			sortM      = [];
			sortM(:,1) = dsum;
			sortM(:,2) = tomatch;
			sortM(:,3) = gr2.ind(tomatch);
			
            % Find best candidate by sorting
			sortM = sortrows(sortM,1);

			tomatch = sortM(:,2);

			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);

	end

end

toBeChecked.gr1 = gr1;
toBeChecked.gr2 = gr2;
    

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Evalute grain match
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

dev.xyz   = NaN(nof_grains1,3);
dev.Rvec  = NaN(nof_grains1,3);
dev.dist  = NaN(nof_grains1,1);
dev.Rdist = NaN(nof_grains1,1);
dev.int   = NaN(nof_grains1,1);
dev.bbys  = NaN(nof_grains1,1);
dev.bbxs  = NaN(nof_grains1,1);
dev.rot   = NaN(3,3,nof_grains1);

for ii = 1:size(match,1)
    jj = match(ii, 2);
    if ~isnan(jj)
        dev.xyz(ii,:)   = gr2.center(jj,:) - gr1.center(ii,:);
        dev.dist(ii,:)  = sqrt(dev.xyz(ii,:) * dev.xyz(ii,:)');
        dev.Rvec(ii,:)  = gr2.R_vector(jj,:) - gr1.R_vector(ii,:);
        dev.Rdist(ii,:) = sqrt(dev.Rvec(ii,:) * dev.Rvec(ii,:)');
        dev.bbxs(ii,1)  = gr2.bbxs(jj) / gr1.bbxs(ii);
        dev.bbys(ii,1)  = gr2.bbys(jj) / gr1.bbys(ii);
        dev.int(ii,1)   = gr2.int(jj)  / gr1.int(ii);
        dev.rot(:,:,ii) = all_g1(:, :, ii).' - gr2_rot(:, :, jj);
    end

end


ok = ~isnan(dev.dist(:,1));

Rvec_med = median(dev.Rvec(ok,:),1);
dxyz_med = median(dev.xyz(ok,:));

ddist  = sqrt(sum(dev.dist(ok,:).^2,1)/sum(ok));
dRdist = sqrt(sum(dev.Rdist(ok,:).^2,1)/sum(ok));


dcm = dc - dxyz_med;


% rotasymm = 0.5*(rot_med0 - rot_med0');
% 
% rmc  = gtMathsRotationMatrixComp([1; 0; 0],'col');
% rotx = gtMathsRotationTensor(asind(rotasymm(3,2)),rmc);
% 
% rmc  = gtMathsRotationMatrixComp([0; 1; 0],'col');
% roty = gtMathsRotationTensor(asind(rotasymm(1,3)),rmc);
% 
% rmc  = gtMathsRotationMatrixComp([0; 0; 1],'col');
% rotz = gtMathsRotationTensor(asind(rotasymm(2,1)),rmc);
% 
% rot_med = rotx*roty*rotz;


rot_med = gtMathsRod2RotMat(Rvec_med);

if isempty(drot)
    drotm = rot_med;
else
    drotm = rot_med*drot;
end


fprintf('\nRemaining median deviation of Rodrigues vectors:\n  [%s]',num2str(Rvec_med));

fprintf('\nRemaining median deviation of grain centers:\n  [%s]',num2str(dxyz_med));

fprintf('\nRemaining root mean square deviation of Rodrigues vectors:\n  %g',dRdist);

fprintf('\nRemaining root mean square deviation of grain center distances:\n  %g',ddist);

%disp('Remaining median misorientation of grains:')
%disp(rot_med)

fprintf('\n\nRecommended total displacement correction (dc):\n  [%s]\n',num2str(dcm));

disp('Recommended total rotation correction (drot):');
disp(drotm);

fprintf('Number of matching grains found:\n  %d \n\n',sum(~isnan(match(:,2))));


end % of main function


%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUB-FUNCTIONS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function tomatch = sfCheckMatching(act, cand, tol, fzone_acc, fzone_ext, pixsize)
        
    tomatch = (1:length(cand.ind))';  % initial indices of grain candidates to be merged
    
    % Delete those indices that don't meet the following constraints:
    
    % Bounding box Y size close enough?
    cons = (cand.bbys > act.bbys/tol.bbys) & ...
        (cand.bbys < tol.bbys*act.bbys);
    cand = sfSelectGrains(cons, cand);
    tomatch(~cons) = [];
    
    % Distance of centers close enough?
    toldist = min(act.bbxs, act.bbys) * pixsize * tol.distf;
    dvec    = repmat(act.center,length(cand.ind),1) - cand.center;
    dist    = sqrt(sum(dvec.*dvec,2));
    cons    = dist < toldist;
    cand    = sfSelectGrains(cons, cand);
    tomatch(~cons) = [];
    
    % Bounding box X size close enough?
    cons = (cand.bbxs > act.bbxs/tol.bbxs) & ...
        (cand.bbxs < tol.bbxs*act.bbxs);
    cand = sfSelectGrains(cons, cand);
    tomatch(~cons) = [];
    
    % Intensity close enough?
    cons = (cand.int > act.int/tol.int) & ...
        (cand.int < tol.int*act.int);
    cand = sfSelectGrains(cons,cand);
    tomatch(~cons) = [];
    
    if isempty(tomatch)
        return;
    end
    
    % Rodriguez vector close enough?
    
    % If any of them is on edge of fundamental zone
    if any([act.R_onedge; cand.R_onedge])  
        
        % Do Rodrigues test from scratch
        
        Rlines = vertcat(act.Rlines{:}, cand.Rlines{:});
        Rids   = vertcat(act.Rids{:}, cand.Rids{:});
        
        candinds = zeros(length(act.Rids{:}),1);
        for jj = 1:length(tomatch)
            candinds = [candinds; jj(ones(length(cand.Rids{jj}),1))];
        end
        
        % Find new Rodrigues vector
        [~, ~, goodlines] = gtCrystRodriguesTestCore(Rlines, Rids, ...
                               fzone_acc, fzone_ext, tol.Rdist);
        
        
        % Check which candidate has all of its 3 reflections accepted
        
        candinds = candinds(goodlines);
        candinds(candinds==0) = [];
        
        unicandinds = unique(candinds);
        
        okinds = false(size(tomatch));
        
        for jj = 1:length(unicandinds)
            if (sum(candinds==unicandinds(jj)) == 3)
                okinds(unicandinds(jj)) = true;
            end
        end
        
        tomatch = tomatch(okinds);
        
    else
        dRvec = repmat(act.R_vector,length(cand.ind),1) - cand.R_vector;
        Rdist = sqrt(sum(dRvec.*dRvec,2));
        cons  = Rdist < tol.Rdist;
        tomatch(~cons) = [];
    end
    
end


function grk = sfSelectGrains(keeplist,gr)

    grk.ind      = gr.ind(keeplist);
    grk.int      = gr.int(keeplist);
    grk.center   = gr.center(keeplist,:);
    grk.R_vector = gr.R_vector(keeplist,:);
    grk.R_onedge = gr.R_onedge(keeplist,:);
    grk.bbys     = gr.bbys(keeplist);
    grk.bbxs     = gr.bbxs(keeplist);
    grk.Rlines   = gr.Rlines(keeplist);
    grk.Rids     = gr.Rids(keeplist);
 
end