Cleaning + improvements: switch to a single function for IVP color map -> gtIVPCmap

2 other new functions are introduced to factorize and speed up the procedure:
- gtVectorLab2Cryst: rotate a lab vector to crystal CS
- gtCrystVector2SST: move a vector from crystal CS to SST zone and give IVP color code

git-svn-id: https://svn.code.sf.net/p/dct/code/trunk@906 4c865b51-4357-4376-afb4-474e03ccb993

6_rendering/gtIVPCmap.m

+function cmap = gtIVPCmap(phaseid, LD, varargin)
+% GTIVPCMAP  Returns a color map for the specified phase and lab direction.
+%
+%     cmap = gtIVPCmap([phaseid, LD, varargin])
+%     ------------------------------------------------------------------------
+%     Makes an IVP colormap of the lab_ direction (LD) in the crystal 
+%     coordinates. This can be specified as a second arguement, otherwise we
+%     assume the z axis [0 0 1]
+%
+%     OPTIONAL INPUT:
+%       phaseid          = <int>     Id of the specific phase {1}
+%       LD               = <double>  Lab direction to study   {[0 0 1]}
+%
+%     OPTIONAL INPUT (varargin as a list of pairs, see parse_pv_pairs.m)
+%       'save'           = <bool>    Save color map to cmap_phase<phaseid>.mat
+%       'crystal_system' = <string>  Legacy to set (force) the crystal system
+%                                    for old datasets
+%   
+%     OUTPUT:
+%       cmap             = <double>  RGB color map, color code is the inverse
+%                                    pole figure in the SST
+%
+%     Version 001 15-11-2012 by YGuilhem
+%       Built by gathering old gtIVPcubCmap and gtIVPhexCmap functions
+
+% Default parameters
+par.save = false;
+par.crystal_system = '';
+par = parse_pv_pairs(par, varargin);
+
+% Set default phaseid
+if ~exist('phaseid', 'var') || isempty(phaseid)
+    phaseid = 1;
+end
+
+% Set default lab direction if LD empty
+if ~exist('LD', 'var') || isempty(LD)
+    LD = [0 0 1];
+end
+
+% Get the r-vectors
+r_vectors = [];
+if exist(fullfile('4_grains', sprintf('phase_%02d', phaseid), 'r_vectors.mat'), 'file')
+    load(fullfile('4_grains', sprintf('phase_%02d', phaseid), 'r_vectors.mat'));
+    r_vectors = r_vectors(:, 2:4);
+elseif exist('r_vectors.mat', 'file')
+    load('r_vectors.mat');
+    r_vectors = r_vectors(:, 2:4);
+elseif exist(fullfile('4_grains', sprintf('phase_%02d', phaseid), 'index.mat'), 'file')
+    grain = [];
+    load(fullfile('4_grains', sprintf('phase_%02d', phaseid), 'index.mat'), 'grain');
+    r_vectors = gtIndexAllGrainValues(grain, 'R_vector', [], 1:3, []);
+else
+    gtError('gtIVPCmap:missing_file', 'Could not find r_vectors or index file!');
+end
+
+% Get the symmetry operators
+if isempty(par.crystal_system)
+    parameters = load('parameters.mat');
+    par.crystal_system = parameters.parameters.cryst.crystal_system;
+end
+tmp = gtCrystGetSymmetryOperators(par.crystal_system);
+if isfield(tmp, 'g3')
+    symm = {tmp.g3};
+else
+    symm = {tmp.g};
+end
+
+% Create the map
+M = length(r_vectors);
+map = zeros(M, 3);
+
+% Get sample to crystal rotation matrices
+sam2crys = gtMathsRod2RotMat(r_vectors);
+
+% Compute LDc = LD expressed in the crystal CS
+LDc = gtVectorLab2Cryst(LD, sam2crys);
+
+% Get the color in the SST triangle
+map = gtCrystVector2SST(LDc, par.crystal_system, symm);
+
+cmap = [ 0 0 0 ; map ];
+
+if par.save
+    fileName = fullfile('6_rendering', sprintf('cmap_phase%02d.mat', phaseid));
+    save(fileName, 'cmap');
+end
+
+end % end of function

6_rendering/gtIVPcubCmap.m

-function cmap = gtIVPcubCmap(grain)
-
-% as for hexagonal materials
-% aim for RED = 001 / GREEN = 101 / BLUE = 111
-
-% get the r-vectors
-r_vectors=gtIndexAllGrainValues(grain, 'R_vector', [], 1:3, []);
-symm = gtCrystGetSymmetryOperators('cubic');
-
-% base the IVP on the lab z axis
-Zlab=[0 0 1];
-
-map = zeros(length(r_vectors+1),3);
-
-for i=1:length(r_vectors)
-    
-    crys2sam=Rod2g(r_vectors(i, :));
-    
-    % z axis in cryst axes \ to go sample->crystal
-    Zc = crys2sam\Zlab';
-    
-    % all equivilents by symm
-    Zc_symm=[];
-    for jj=1:length(symmm)
-        Zc_symm(jj, :)=(symmm(jj).g*Zc)';
-    end
-    Zc_symm=[Zc_symm; -Zc_symm];
-    
-    % define angles phi - rotation around 001 axis, projected on the xy plane
-    % chi - rotation around 010, projected on xz plane
-    Zc_chi=rad2deg(atan2(Zc_symm(:,1), Zc_symm(:,3)));
-    Zc_phi=rad2deg(atan2(Zc_symm(:,2), Zc_symm(:,1)));
-    
-    % pick equiv that is in our triangle
-    ndx=find(Zc_chi>=0 & Zc_chi<45 & Zc_phi>=0 & Zc_phi<45);
-   
-    % Z symm equivilent in the part of the crystal we are interested in
-    Zp(i, :)=Zc_symm(ndx, :);
-    
-    % colourmap
-    chi=Zc_chi(ndx);
-    phi=Zc_phi(ndx);
-    
-    % these shouldn't come out negative!   
-    red=(45-chi)/45;
-    green=(chi/45)*(45-phi)/45;
-    blue=(chi/45)*phi/45;
- 
-    % saturate colours    
-    rgb=[red green blue];
-    rgb=rgb./(max(rgb));
-
-    map(i+1, :)=rgb;
-    
-end
-
-cmap = map;
-
-end % end of function

6_rendering/gtIVPhexCmap.m

-function cmap = gtIVPhexCmap(phaseid, varargin)
-%     cmap = gtIVPhexCmap(phaseid, varargin)
-%     --------------------------------------
-%     Makes an IVP colourmap of the lab_direction (LD) in the crystal coordinates.
-%     this can be specified as a second arguement, otherwise we assume the z
-%     axis [0 0 1]
-
-r_vectors = [];
-if exist(fullfile('4_grains',sprintf('phase_%02d',phaseid),'r_vectors.mat'),'file')
-    load(fullfile('4_grains',sprintf('phase_%02d',phaseid),'r_vectors.mat'));
-    r_vectors = r_vectors(:,2:4);
-elseif exist('r_vectors.mat','file')
-    load('r_vectors.mat');
-    r_vectors = r_vectors(:,2:4);
-elseif exist(fullfile('4_grains',sprintf('phase_%02d',phaseid),'index.mat'),'file')
-    grain = [];
-    load(fullfile('4_grains',sprintf('phase_%02d',phaseid),'index.mat'),'grain');
-    r_vectors = gtIndexAllGrainValues(grain, 'R_vector', [], 1:3, []);
-end
-
-if isempty(varargin)
-    LD=[0 0 1];
-else
-    LD=varargin{1};
-end
-
-
-% make some symmetry operators
-% six fold rotation around z
-for i=1:6
-    a=(i-1)*60;
-    symm(i).g3 = [cosd(a) sind(a) 0;...
-                  -sind(a) cosd(a) 0;...
-                  0 0 1];
-    % include the mirror
-    symm(i+6).g3=[1 0 0; 0 -1 0; 0 0 1]*symm(i).g3;
-end
-
-map = zeros(length(r_vectors+1),3);
-
-for i=1:length(r_vectors)
-    
-    crys2sam=Rod2g(r_vectors(i, :));
-    
-    % RD in cryst axes - \ to do inverse, sample->crystal
-    LDc = crys2sam\LD';
-    % all equivilents by sym
-    LDc_symm=[];
-    for jj=1:12
-        LDc_symm(jj, :)=(symm(jj).g3*LDc)';
-    end
-    LDc_symm=[LDc_symm; -LDc_symm];
-    
-    LDc_psi=acosd(LDc_symm(:,3));
-    LDc_phi=rad2deg(atan2(LDc_symm(:,2), LDc_symm(:,1)));
-    
-    
-    % pick equiv that is in out triangle
-    ndx=find(LDc_psi<90 & LDc_phi>=60 & LDc_phi<90);
-    
-    % RD symm equivilent in the part of the crystal we are interested in
-    LDp(i, :)=LDc_symm(ndx, :);
-    
-    % colourmap
-    psi=LDc_psi(ndx);
-    phi=LDc_phi(ndx);
-    
-    % for hexagonal colour map need only psi and phi
-    red=(90-psi)/90;
-    green=(psi/90)*(90-phi)/30;
-    blue=(psi/90)*(phi-60)/30;
-    
-    rgb=[red green blue];
-    % do we saturate colours?
-    if 1
-        rgb=rgb./(max(rgb));
-    end
-    
-    map(i+1, :)=rgb;
-    
-end
-
-map(1,:) = [0 0 0];
-
-cmap = map;
-
-save('6_rendering/cmap.mat','cmap');
-
-end % end of function
-
-
-
-

zUtil_Cryst/gtCrystGetSymmetryOperators.m

@@ -46,9 +46,10 @@ if ~exist('crystal_system', 'var') || isempty(crystal_system)
     end
 elseif ~exist('spacegroup', 'var') || isempty(spacegroup)
     if exist('parameters.mat', 'file')
-        disp('Guessing spacegroup from parameters file...');
         parameters = load('parameters.mat');
-        spacegroup = parameters.parameters.cryst.spacegroup;
+        if isfield(parameters.parameters, 'cryst')
+            spacegroup = parameters.parameters.cryst.spacegroup;
+        end
     end
 end
 

zUtil_Cryst/gtCrystVector2SST.m

+function [rgb, Vsst] = gtCrystVector2SST(Vc, crystal_system, symm, saturation)
+% Vc             : row vector expressed in crystal frame
+% crystal_system : type of crystal_system ('cubic'|'hexagonal')
+% symm           : cell containing every symmetry operators
+% saturation     : saturate or not the RGB output {true}
+%
+% rgb    : RGB color in the sst
+% Vsst   : Vc moved vector in the SST triangle
+% Vsst_p : Vc moved vector in the SST triangle, projected on equatorial plane
+
+    if ~exist('saturation', 'var') || isempty(saturation)
+        saturation = true;
+    end
+
+    %Vc_norm = sqrt(sum(Vc.^2, 2));
+    %Vc = Vc ./ Vc_norm(:, [1 1 1]);
+
+    % All equivalents by symm
+    M = size(Vc, 1);
+    % New method
+    N = length(symm);
+    for isymm=1:N
+        symm{N+isymm} = -symm{isymm};
+    end
+    Vc_symm = gtVectorLab2Cryst(Vc, symm);
+
+    % Old method
+    %N = 2*length(symm);
+    %Vc_symm = zeros(N, 3, M);
+    %for iV=1:M
+    %    V = Vc(iV, :);
+    %    tmp = gtVectorLab2Cryst(V, symm); 
+    %    Vc_symm_tmp = [ tmp ; -tmp ];
+    %    Vc_symm(:, :, iV) = reshape(Vc_symm_tmp, [N 3 1]);
+    %end
+
+    x = Vc_symm(:, 1, :);
+    y = Vc_symm(:, 2, :);
+    z = Vc_symm(:, 3, :);
+
+    aaa = find(x == y);
+    bbb = find(x == z);
+    ccc = find(y == z);
+
+    eps = 1.e-1;
+    if size(aaa)
+        x(aaa) = x(aaa)+eps;
+        y(aaa) = y(aaa)-eps;
+        %disp('x / y')
+        %disp([x(aaa)'; y(aaa)'; z(aaa)'])
+        %disp([x(aaa)'; y(aaa)'; z(aaa)'])
+    end
+    if size(bbb)
+        x(bbb) = x(bbb)+eps;
+        z(bbb) = z(bbb)-eps;
+        %disp('x / z')
+        %disp([x(bbb)'; y(bbb)'; z(bbb)'])
+        %disp([x(bbb)'; y(bbb)'; z(bbb)'])
+    end
+    if  size(ccc)
+        y(ccc) = y(ccc)+eps;
+        z(ccc) = z(ccc)-eps;
+        %disp('y / z')
+        %disp([x(ccc)'; y(ccc)'; z(ccc)'])
+        %disp([x(ccc)'; y(ccc)'; z(ccc)'])
+    end
+
+    aaa = find(x == y);
+    bbb = find(x == z);
+    ccc = find(y == z);
+
+    % Define angles
+
+    % phi : rotation around 001 axis,
+    %       from 100 axis to Vc vector,
+    %       projected on (100,010) plane
+    %Vc_phi = rad2deg(atan2(Vc_symm(:, 2, :), Vc_symm(:, 1, :)));
+    Vc_phi = rad2deg(atan2(y, x));
+
+    % chi : rotation around 010 axis,
+    %       from 001 axis to Vc vector,
+    %       projected on (100,001) plane
+    %Vc_chi = rad2deg(atan2(Vc_symm(:, 1, :), Vc_symm(:, 3, :)));
+    Vc_chi = rad2deg(atan2(x, z));
+
+    % psi : angle from 001 axis to Vc vector
+    Vc_psi = acosd(z);
+    %Rxy = sqrt(1 - Vc_symm(:, 3, :).*Vc_symm(:, 3, :));
+    %Vc_psi = rad2deg(atan2(Rxy, Vc_symm(:, 3, :)));
+
+    switch crystal_system
+    case {'cubic', 'cub'}
+        anglesR   = 45 - Vc_chi; % red color proportional to (45-chi)
+        minAngleR =  0;
+        maxAngleR = 45;
+        anglesB   = Vc_phi; % blue color proportional to phi
+        minAngleB =  0;
+        maxAngleB = 45;
+    case {'hexagonal', 'hex'}
+        anglesR   = 90 - Vc_psi; % red color proportional to (90 - psi)
+        minAngleR =  0;
+        maxAngleR = 90;
+        anglesB   = Vc_phi; % blue color proportional to phi
+        minAngleB =  0;
+        maxAngleB = 30;
+    otherwise
+        gtError('gtCrystVector2SST:wrong_crystal_system', ...
+                'Unsupported crystal system!');
+    end
+
+    % Pick equiv that is in our triangle
+    iSST = anglesR>=minAngleR & anglesR<maxAngleR & anglesB>=minAngleB & anglesB<maxAngleB;
+
+    ok = squeeze(sum(iSST, 1));
+    ok_no = (ok == 0);
+    ok_too = (ok > 1);
+
+    if any(ok_no)
+        warning('Problem when moving to SST triangle, no indices found for vectors:');
+        %gtError('gtCrystVector2SST:no_solution', ...
+        %    'Problem when moving to SST triangle, no indices found for vectors:');
+        disp('  indices:');
+        disp(find(ok_no)');
+        disp('  vectors:');
+        disp(Vc(find(ok_no), :));
+        iSST(1, :, ok_no) = 1;
+        anglesR(1, :, ok_no) = 0;
+        anglesB(1, :, ok_no) = 0;
+    end
+    if any(ok_too)
+        warning('Might be a problem when moving to SST triangle, too many indices found!');
+        disp('  indices:');
+        disp(find(ok_too)');
+        disp('  vectors:');
+        disp(Vc(find(ok_too), :));
+        i_too = find(iSST(:, :, ok_too));
+        iSST(i_too(2:end), :, ok_too) = 0;
+    end
+
+    angleR = anglesR(iSST);
+    angleB = anglesB(iSST);
+    if nargout > 1
+        Vsst = reshape(Vc_symm(iSST(:, [1 1 1], :)), [3 M]).';
+        if nargout > 2
+            denom = 1 + Vsst(:, 3);
+            Vsst_p = Vsst(:, 1:2) ./ denom(:, [1 1]);
+        end
+    end
+
+    % Color map
+    red   = angleR/maxAngleR;
+    green = (maxAngleR - angleR)/maxAngleR .* (maxAngleB - angleB)/maxAngleB;
+    blue  = (maxAngleR - angleR)/maxAngleR .* angleB/maxAngleB;
+
+    rgb = [red green blue];
+    if saturation
+        sat = max(rgb, [], 2);
+        rgb = rgb./sat(:, [1 1 1]);
+    end
+
+end % end of function

zUtil_Cryst/gtVectorLab2Cryst.m

+function Vc = gtVectorLab2Cryst(Vs, g, diag)
+% GTVECTORLAB2CRYST  Express lab row vector(s) to crystal CS
+%
+%     Vc = gtVectorLab2Cryst(Vs, g)
+%     -------------------------------------------------------------------------
+%
+%     INPUT:
+%       Vs = <double>   Vector(s) in the lab CS.
+%       g  = <cell>     Cell containing crystal orientation matrix(ces).
+%
+%     OUTPUT:
+%       Vc  = <double>  Vector Vs expressed in the crystal CS.
+%
+%     Note:
+%     This function can be used in 4 different modes:
+%       - Vs = 1x3 matrix (1 row vector)  , g = {N}x3x3 matrices
+%         output N vectors expressed in the crystal CS defined by g matrices
+%
+%       - Vs = Mx3 matrix (M row vectors) , g = 3x3 matrix
+%         output M vectors expressed in the crystal CS defined by g matrix
+%
+%       - Vs = Mx3 matrix (M row vectors) , g = {N}x3x3 matrices
+%         output Nx3xM matrix containing NxM row vectors, which are all the Vs
+%         vectors expressed in the crystals CS defined by g matrices
+%
+%       - diag mode (need diag argument, otherwise previous case)
+%         Vs = Mx3 matrix (M row vectors) , g = {N}3x3 matrices (M=N)
+%         output M (equal to N) vectors in the crystal CS defined by g matrices
+%
+%     Version 001 17-10-2012 by YGuilhem
+
+if ~exist('diag', 'var') || isempty(diag)
+    diag = false;
+else
+    diag = logical(diag);
+end
+
+% Get number of vectors and orientation matrices
+M = size(Vs, 1);
+if iscell(g)
+    N = length(g);
+else
+    N = 1;
+end
+
+if M==1 && N >= 1
+    % Convert cell to array of matrices
+    sam2crys = cat(3, g{:});
+
+    % Expand Vs in 3rd dimensiodifferent n
+    Vext = Vs([1 1 1], :, ones(N, 1));
+
+    % Compute Vc and permute to get in in row vectors
+    Vc = permute(sum(sam2crys.*Vext, 2), [3 1 2]);
+
+elseif M > 1 && N == 1
+    % Convert cell to array of matrices
+    sam2crys = g(:, :, ones(M, 1));
+
+    % Expand Vs in 3rd dimension
+    Vext = reshape(Vs.', [1 3 M]);
+
+    % Compute Vc and permute to get in in row vectors
+    Vc = permute(sum(sam2crys.*Vext([1 1 1], :, :), 2), [3 1 2]);
+
+elseif M == N && sync
+    % Convert cell to array of matrices
+    sam2crys = cat(3, g{:});
+
+    % Expand Vs in 3rd dimension
+    Vext = reshape(Vs.', [1 3 M]);
+
+    % Compute Vc and permute to get in in row vectors
+    Vc = permute(sum(sam2crys.*Vext([1 1 1], :, :), 2), [3 1 2]);
+
+elseif M>0 && N>0
+    % Convert cell to array of matrices
+    sam2crys = cat(3, g{:});
+    sam2crys2 = sam2crys(:, :, :, ones(M, 1));
+
+    % Expand Vs in 3rd dimension
+    Vext = reshape(Vs.', [1 3 1 M]);
+    VV = Vext([1 1 1], :, ones(N, 1), :);
+
+    tmp = sum(sam2crys2.*VV, 2);
+    Vc = permute(squeeze(tmp), [2 1 3]);
+    %Vc = permute(squeeze(sum(sam2crys2.*VV, 2)), [2 1 3]);
+
+else
+    gtError('gtVectorLab2Cryst:wrong_input_size', ...
+            ['Cannot handle the case where size(Vs)=' num2str(size(Vs)) ...
+            ' and length(g)=' num2str(N)]);
+end
+
+end % end of function