EBSD: added a few functions for interacting with EBSD data

Signed-off-by: Nicola Vigano <nicola.vigano@esrf.fr>

zUtil_Deformation/gtDefCalibrateEBSDMap.m

+function [rot_angle, rot_axis, rot_tensor] = gtDefCalibrateEBSDMap(EBSD_r_map, points, gr, symm)
+
+    if (~exist('symm', 'var'))
+        symm = gtCrystGetSymmetryOperators();
+    end
+
+    num_points = numel(gr);
+    for ii = 1:num_points
+%         ref_e_vec = extract_from_map(EBSD_e_map, points{ii});
+%         ref_r_vec = gtMathsEuler2RodInFundZone(ref_e_vec', symm);
+        ref_r_vec = extract_from_map(EBSD_r_map, points{ii});
+        ref_r_vec = gtMathsRod2RodInFundZone(ref_r_vec', symm);
+        % Sometimes the axes definition is not the same
+%         ref_r_vec = [-ref_r_vec(2), ref_r_vec(3), -ref_r_vec(1)];
+
+        gr_r_vec = gr(ii).R_vector;
+        gr_r_vec = gtMathsRod2RodInFundZone(gr_r_vec', symm);
+
+        [rot_vec(ii, :)] = get_angle_axis(ref_r_vec, gr_r_vec, symm); %#ok<AGROW>
+%         [rot_angle(ii), rot_axis(ii, :)] = get_angle_axis(ref_r_vec, gr_r_vec, symm); %#ok<AGROW>
+    end
+
+%     rot_angle = sum(rot_angle) / num_points;
+%     rot_axis = sum(rot_axis, 1) / num_points;
+%     rot_axis = rot_axis / norm(rot_axis);
+
+    rot_vec = sum(rot_vec, 1) / num_points;
+    rot_angle = 2 * atand(norm(rot_vec));
+    rot_axis = rot_vec / norm(rot_vec);
+
+    rotcomp = gtMathsRotationMatrixComp(rot_axis', 'col');
+    rot_tensor = gtMathsRotationTensor(rot_angle, rotcomp);
+end
+
+function r_vec_3 = get_angle_axis(r_vec_1, r_vec_2, symm)
+% function [angle, axis] = get_angle_axis(r_vec_1, r_vec_2, symm)
+%     [angle, axis] = gtDisorientation(r_vec_1', r_vec_2', symm);
+
+    g1 = gtMathsRod2OriMat(r_vec_1);
+    g2 = gtMathsRod2OriMat(r_vec_2);
+
+    M = g1' * g2;
+
+    r_vec_3 = gtMathsOriMat2Rod(M);
+    r_vec_3 = gtMathsRod2RodInFundZone(r_vec_3, symm);
+%     angle = 2 * atand(norm(r_vec_3));
+%     axis = r_vec_3 ./ norm(r_vec_3);
+end
+
+function r_vec = extract_from_map(EBSD_r_map, points)
+    kernel = 7;
+    k_dist = (kernel - 1) / 2;
+
+    num_points = size(points, 1);
+    r_vecs = zeros(num_points, 3);
+
+    for ii = 1:num_points
+        point = points(ii, :);
+
+        ranges = [point(2)-k_dist, point(2)+k_dist; point(1)-k_dist, point(1)+k_dist];
+
+        temp_r_vecs = EBSD_r_map(ranges(1, 1):ranges(1, 2), ranges(2, 1):ranges(2, 2), :);
+        temp_r_vecs = sum(sum(temp_r_vecs, 1), 2) / (size(temp_r_vecs, 1) * size(temp_r_vecs, 2));
+        r_vecs(ii, :) = reshape(temp_r_vecs, 1, []);
+    end
+
+    r_vec = sum(r_vecs, 1) / num_points;
+end

zUtil_Deformation/gtReadEBSDMapEulerCSVFile.m

+function [EBSD_e_map, EBSD_r_map] = gtReadEBSDMapEulerCSVFile(filename, axes_perm)
+% function [EBSD_e_map, EBSD_r_map] = gtReadEBSDMapEulerCSVFile(filename, axes_perm)
+%   Axes permutation allows to easily deal with different reference
+%   systems, like for instance:
+%       DCT Axes <-> EBSD Axes ([x y z] <-> [x'y'z'])
+%           x     =     -y'
+%           y     =      z'
+%           z     =     -x'
+%   Translates to: axes_perm = [-2 3 -1]
+%
+%   NOTE: Only the Rodrigues vector map has the axes changed!
+%
+
+    c = tic();
+    fprintf('Reading from file: "%s"..', filename);
+
+    fid = fopen(filename);
+    sizes = fscanf(fid, 'Width,%d\nHeight,%d\nX Position,Y Position,Eul 1,Eul 2,Eul 3\n')';
+    input = fscanf(fid, '%d,%d,%f,%f,%f');
+    fclose(fid);
+
+    fprintf('\b\b, Done in %f seconds.\n', toc(c));
+    c = tic();
+    fprintf('Converting to Euler and Rodriguez maps..');
+
+    xx = input(1:5:end)+1;
+    yy = input(2:5:end)+1;
+    ea1 = input(3:5:end);
+    ea2 = input(4:5:end);
+    ea3 = input(5:5:end);
+
+    sizes = sizes([2 1]);
+    EBSD_e_map = zeros([sizes 3], 'single');
+    ii = sub2ind([sizes, 3], yy, xx, 1 * ones(size(xx)));
+    EBSD_e_map(ii) = ea1;
+    ii = sub2ind([sizes, 3], yy, xx, 2 * ones(size(xx)));
+    EBSD_e_map(ii) = ea2;
+    ii = sub2ind([sizes, 3], yy, xx, 3 * ones(size(xx)));
+    EBSD_e_map(ii) = ea3;
+
+    EBSD_e_map = EBSD_e_map * 180 / pi;
+
+    dmvol_EBSD_e_map = permute(EBSD_e_map, [1 2 4 3]);
+    [gvdm_EBSD_e_map, size_dmvol_EBSD_map] = gtDefDmvol2Gvdm(dmvol_EBSD_e_map);
+    gvdm_EBSD_r_map = gtMathsEuler2Rod(gvdm_EBSD_e_map);
+    dmvol_EBSD_r_map = gtDefGvdm2Dmvol(gvdm_EBSD_r_map, size_dmvol_EBSD_map);
+    EBSD_r_map = reshape(dmvol_EBSD_r_map, [sizes 3]);
+
+    EBSD_r_map = EBSD_r_map(:, :, abs(axes_perm));
+    negative = axes_perm < 0;
+    EBSD_r_map(:, :, negative) = - EBSD_r_map(:, :, negative);
+
+    fprintf('\b\b, Done in %f seconds.\n', toc(c));
+end

zUtil_Deformation/plotting/gt6DRenderEBSDMapIPF.m

+function gt6DRenderEBSDMapIPF(EBSD_r_map, phase_id, varargin)
+
+    conf = struct( ...
+        'ipf_plane_normal', [], ...
+        'flip_ud', false, ...
+        'flip_lr', false, ....
+        'clims', [], ...
+        'lims_x', [], ...
+        'pixel_to_cm', 0.05, ...
+        'borders', 'no', ...
+        'colorbar', 'no' );
+    conf = parse_pv_pairs(conf, varargin);
+
+    if (~isempty(conf.lims_x))
+        % Note the inversion between x/y because of matlab's transpose problem
+        EBSD_r_map = EBSD_r_map(:, conf.lims_x(1):conf.lims_x(2), :);
+    end
+
+    slice_ids = all(EBSD_r_map ~= 0, 3);
+    indx = find(slice_ids);
+    r_vecs = reshape(EBSD_r_map, [], 3);
+    r_vecs = r_vecs(indx, :);
+
+    cmap = get_cmap(r_vecs, phase_id, conf.ipf_plane_normal);
+
+    slice_R = zeros(size(slice_ids));
+    slice_G = zeros(size(slice_ids));
+    slice_B = zeros(size(slice_ids));
+
+    slice_R(indx) = cmap(:, 1);
+    slice_G(indx) = cmap(:, 2);
+    slice_B(indx) = cmap(:, 3);
+
+    if (conf.flip_lr)
+        slice_R = fliplr(slice_R);
+        slice_G = fliplr(slice_G);
+        slice_B = fliplr(slice_B);
+    end
+    if (conf.flip_ud)
+        slice_R = flipud(slice_R);
+        slice_G = flipud(slice_G);
+        slice_B = flipud(slice_B);
+    end
+
+    f = figure();
+    ax = axes('Parent', f);
+    im = imagesc(cat(3, slice_R, slice_G, slice_B), 'Parent', ax);
+
+    set(f, 'Units', 'centimeters')
+    img_size = size(slice_ids) * conf.pixel_to_cm;
+    img_size = img_size([2 1]);
+    if (strcmpi(conf.borders, 'no'))
+        position = [0, 0, img_size];
+        set(f, 'Position', position)
+        set(f, 'Paperposition', position)
+
+        set(ax, 'Units', 'normalized')
+        set(ax, 'Position', [0, 0, 1, 1])
+    else
+        set(ax, 'Units', 'centimeters')
+        if (strcmpi(conf.colorbar, 'no'))
+            position = [0, 0, img_size] + [0, 0, 2, 1.5];
+            set(f, 'Position', position)
+            set(f, 'Paperposition', position)
+        else
+            position = [0, 0, img_size] + [0, 0, 3, 2];
+            set(f, 'Position', position)
+            set(f, 'Paperposition', position)
+
+            cb = colorbar('peer', ax, 'location', 'EastOutside');
+            set(cb, 'Units', 'centimeters')
+            position = [img_size(1)+1.35, 1, 0.65, img_size(2)];
+            set(cb, 'Position', position)
+        end
+        position = [1, 1, img_size];
+        set(ax, 'Position', position)
+    end
+end
+
+function cmap = get_cmap(r_vecs, phase_id, ipf_plane_normal)
+    p = gtLoadParameters();
+    cryst_system = p.cryst(phase_id).crystal_system;
+    cryst_spacegroup = p.cryst(phase_id).spacegroup;
+    symm = gtCrystGetSymmetryOperators(cryst_system, cryst_spacegroup);
+
+    [cmap, ~, ~] = gtIPFCmap(phase_id, ipf_plane_normal, ...
+        'r_vectors', r_vecs, ...
+        'crystal_system', cryst_system, ...
+        'background', false, ...
+        'symm', symm);
+end