gtCrystAxisPoles : changed a bit the function using pars (varargin input),...

gtCrystAxisPoles : changed a bit the function using pars (varargin input), cleaned the code. Added 'axesRot' option to rotate sample CS to match for example pole figures. Added comments and variable names more reasonable.

Signed-off-by: Laura Nervo <laura.nervo@esrf.fr>

6_rendering/gtCrystAxisPoles.m

-function [poles, cmap, valids] = gtCrystAxisPoles(phaseid, pl_norm, r_vectors, axesRot, save_flag)
+function [poles, cmap, valids] = gtCrystAxisPoles(phaseid, pl_norm, varargin)
 % GTCRYSTAXISPOLES  Makes a colourmap according to a crystallographic axis
-%                   orientation and gives the colormap.
-%     [poles, cmap, valids] = gtCrystAxisPoles([phaseid], [pl_norm], r_vectors, [axesRot], [save_flag])
-%     -------------------------------------------------------------------------------------------------
+%                   orientation and gives the colormap
+%     [poles, cmap, valids] = gtCrystAxisPoles([phaseid], [pl_norm], varargin)
+%     ------------------------------------------------------------------------
 %     Calculates poles wrt a crystallographic direction using r_vectors.
-%     Loads parameters.mat:cryst(phaseid):spacegroup, latticepar, crystal_system
+%     Loads parameters.mat:cryst(phaseid):latticepar, crystal_system
 %
 %     INPUT:
-%       phaseid   = <int>        phase number {1}
-%       pl_norm   = <double>     hkil of the plane normal wrt calculate the poles (1,3(4)) 
-%                                {[0 0 0 2]},[2 -1 -1 0] for hexagonal, [0 0 1],[0 1 0] for cubic...
-%       r_vectors = <double>     Rodrigues vectors (N,3)
-%       axesRot   = <double>     vector to trasform sample reference. It rotates
-%                                poles basing on component order in this vector {[1 2 3]}
-%       save_flag = <logical>    flag to save file {false}
+%       phaseid          = <int>        phase number {1}
+%       pl_norm          = <double>     hkil of the plane normal wrt calculate the poles (1,3(4)) 
+%                                       {[0 0 0 2]},[2 -1 -1 0] for hexagonal, [0 0 1],[0 1 0] for cubic...
+%
+%     OPTIONAL INPUT (varargin as a list of pairs, see parse_pv_pairs.m)
+%       'crystal_system' = <string>     Legacy to set (force) the crystal system
+%                                       for old datasets
+%       'latticepar'     = <double>     Lattice parameters to convert from hex to
+%                                       cartesian base
+%       'r_vectors'      = <double>     Rodrigues vectors as N x 3 array
+%       'axesRot'        = <double>     vector to trasform sample reference. It rotates
+%                                       poles basing on component order in this vector {[1 2 3]}
+%       'background'     = <logical>    Keep or remove the background color in first row {true}
+%       'save'           = <logical>    flag to save file {false}
 %
 %     OUTPUT:
-%       poles     = <double>     list of poles direction of all the grains (M,4)
-%                                We keep the index of the reference r_vector (first column)
-%                                NB. M is different from N because of the crystal symmetry.
-%       cmap      = <double>     poles colourmap (N,4)
-%       valids    = <double>     poles list (only one per orientation) (N,3)
+%       poles            = <double>     list of poles direction of all the grains (M,4)
+%                                       We keep the index of the reference r_vector (first column)
+%                                       NB. M is different from N because of the crystal symmetry.
+%       cmap             = <double>     poles colourmap (N,4)
+%       valids           = <double>     poles list (only one per orientation) (N,3)
 %
 %
 %     Version 005 18-09-2013 by LNervo
@@ -36,78 +43,103 @@ function [poles, cmap, valids] = gtCrystAxisPoles(phaseid, pl_norm, r_vectors, a
 %     Version 002 28-06-2012 by LNervo
 %       Update to version 2 of parameters; cleaning formatting
 
+% Default parameters
+par.crystal_system = '';
+par.latticepar     = [];
+par.r_vectors      = [];
+par.axesRot        = [1 2 3];
+par.background     = true;
+par.save           = false;
 
-if ~exist('r_vectors','var') || isempty(r_vectors)
-    gtError('gtCrystAxisPoles:missingArgument','The argument r_vectors is empty...')
-end
+par = parse_pv_pairs(par, varargin);
+
+% Set default phaseid
 if ~exist('phaseid','var') || isempty(phaseid)
     phaseid = 1;
 end
 
-parameters = [];
-load('parameters.mat');
-crystal_system = parameters.cryst(phaseid).crystal_system;
-latticepar = parameters.cryst(phaseid).latticepar;
-
+% Set default crystal direction if pl_norm empty
 if ~exist('pl_norm','var') || isempty(pl_norm)
-  if strcmp(crystal_system, 'hexagonal')
-      pl_norm = [0 0 0 2];
-  else
-      pl_norm = [0 0 1];
-  end
-end
-if ~exist('axesRot','var') || isempty(axesRot)
-    axesRot = [1 2 3];
-end
-if ~exist('save_flag','var') || isempty(save_flag)
-    save_flag = false;
+    pl_norm = [0 0 1];
 end
 
-pl_norm_str = num2str(pl_norm);
-pl_norm_str = strrep(pl_norm_str,' ','');
+r_vectors = [];
+phaseDir = fullfile('4_grains', sprintf('phase_%02d', phaseid));
+if ~isempty(par.r_vectors)
+    r_vectors = par.r_vectors;
+elseif exist(fullfile(phaseDir, 'r_vectors.mat'), 'file')
+    load(fullfile(phaseDir, 'r_vectors.mat'));
+elseif exist(fullfile(phaseDir, 'index.mat'), 'file')
+    grain = [];
+    load(fullfile(phaseDir, 'index.mat'), 'grain');
+    r_vectors = gtIndexAllGrainValues(grain, 'R_vector', [], 1, 1:3);
+else
+    gtError('gtCrystAxisPoles:missing_file', 'Could not find r_vectors.mat or index.mat file!');
+end
 
+% Resize r_vectors if needed
 if size(r_vectors,2) == 4
     r_vectors = r_vectors(:,2:4);
 end
 
-% take all the symmetric equivalents
-all_hkils = gtCrystSignedHKLs(pl_norm, gtCrystGetSymmetryOperators(crystal_system));
-% finds the original pl_norm among the list of its symmetries
-[~,ind,~] = findDifferentRowIntoMatrix(all_hkils, pl_norm);
-% going to cartesian reciprocal space + normalisation of the cartesian space
-normalised_hkls = gtHex2Cart(all_hkils, latticepar);
+% Get the crystal system
+if isempty(par.crystal_system)
+    parameters = load('parameters.mat');
+    par.crystal_system = parameters.parameters.cryst(phaseid).crystal_system;
+end
+% Get the lattice parameters
+if isempty(par.latticepar)
+    parameters = load('parameters.mat');
+    par.latticepar = parameters.parameters.cryst(phaseid).latticepar;
+end
 
-poles = zeros(0,4);
-valid_poles = zeros(0,4);
+% convert to Miller-Bravais notation if the case
+if strcmp(par.crystal_system, 'hexagonal') && size(pl_norm,2) == 3
+    pl_norm(4) = pl_norm(3);
+    pl_norm(3) = -pl_norm(1)-pl_norm(2);
+end
 
 % Compute all orientation matrices g
+%   Gcry = g Gsam   --> Gsam = Gcry inv(g)
+%   g = U = gtMathsRod2OriMat(r_vectors.')
 all_g = gtMathsRod2OriMat(r_vectors.');
 
-% Express normalised hkls in cartesian SAMPLE CS
-all_normals = gtVectorCryst2Lab(normalised_hkls, all_g);
 
+% take all the symmetric equivalents
+pl_norm_symm = gtCrystSignedHKLs(pl_norm, gtCrystGetSymmetryOperators(par.crystal_system));
+% finds the original pl_norm among the list of its symmetries
+[~,ind,~] = findDifferentRowIntoMatrix(pl_norm_symm, pl_norm);
+% going to cartesian reciprocal space + normalisation of the cartesian space
+% PLnorm = pl_norm normalised and symmetrised
+PLnorm = gtHex2Cart(pl_norm_symm, par.latticepar);
+% Express normalised hkls (PLs) in cartesian SAMPLE CS
+PLs = gtVectorCryst2Lab(PLnorm, all_g);
+
+
+poles = zeros(0,4);
+valid_poles = zeros(0,4);
 for ii = 1:length(all_g)
-    normals = all_normals(:, :, ii);
+    normals = PLs(:, :, ii);
 
     n_bkp = normals;
     % rotate
-    if axesRot(1) < 0
-        normals(:,1) = -n_bkp(:,axesRot(1));
+    if par.axesRot(1) < 0
+        normals(:,1) = -n_bkp(:,par.axesRot(1));
     else
-        normals(:,1) = n_bkp(:,axesRot(1));
+        normals(:,1) = n_bkp(:,par.axesRot(1));
     end
-    if axesRot(2) < 0
-        normals(:,2) = -n_bkp(:,axesRot(2));
+    if par.axesRot(2) < 0
+        normals(:,2) = -n_bkp(:,par.axesRot(2));
     else
-        normals(:,2) = n_bkp(:,axesRot(2));
+        normals(:,2) = n_bkp(:,par.axesRot(2));
     end
-    if axesRot(3) < 0
-        normals(:,3) = -n_bkp(:,axesRot(3));
+    if par.axesRot(3) < 0
+        normals(:,3) = -n_bkp(:,par.axesRot(3));
     else
-        normals(:,3) = n_bkp(:,axesRot(3));
+        normals(:,3) = n_bkp(:,par.axesRot(3));
     end
 
-    for jj = 1:size(all_normals, 1)
+    for jj = 1:size(PLs, 1)
         % takes the positive poles and saves the grainid which it refers to
         if (normals(jj, 3) >= 0)
             poles(end+1, 1) = ii;
@@ -130,13 +162,22 @@ clear tmp
 
 if nargout > 1
     cmap = gtCrystAxisCmap(valid_poles);
+    
+    if ~par.background
+        cmap(1,:) = [];
+    end
+
     if nargout > 2
         valids = valid_poles;
     end
 end
 
 
-if save_flag
+if par.save
+  
+    pl_norm_str = num2str(pl_norm);
+    pl_norm_str = strrep(pl_norm_str,' ','');
+
     % saves the list of poles to be used in the pole figure
     fname = fullfile('6_rendering',['poles_' pl_norm_str '.mat']);
     save(fname,'poles'); % list of the poles