From 7d8ab800ec95a05124fea33f7a2ea2ae06224cec Mon Sep 17 00:00:00 2001
From: Nicola Vigano <nicola.vigano@esrf.fr>
Date: Wed, 6 May 2015 17:18:07 +0200
Subject: [PATCH] EBSD: added a few functions for interacting with EBSD data
Signed-off-by: Nicola Vigano <nicola.vigano@esrf.fr>
---
zUtil_Deformation/gtDefCalibrateEBSDMap.m | 68 ++++++++++++++
zUtil_Deformation/gtReadEBSDMapEulerCSVFile.m | 54 +++++++++++
.../plotting/gt6DRenderEBSDMapIPF.m | 91 +++++++++++++++++++
3 files changed, 213 insertions(+)
create mode 100755 zUtil_Deformation/gtDefCalibrateEBSDMap.m
create mode 100755 zUtil_Deformation/gtReadEBSDMapEulerCSVFile.m
create mode 100644 zUtil_Deformation/plotting/gt6DRenderEBSDMapIPF.m
diff --git a/zUtil_Deformation/gtDefCalibrateEBSDMap.m b/zUtil_Deformation/gtDefCalibrateEBSDMap.m
new file mode 100755
index 00000000..191e61fa
--- /dev/null
+++ b/zUtil_Deformation/gtDefCalibrateEBSDMap.m
@@ -0,0 +1,68 @@
+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
diff --git a/zUtil_Deformation/gtReadEBSDMapEulerCSVFile.m b/zUtil_Deformation/gtReadEBSDMapEulerCSVFile.m
new file mode 100755
index 00000000..1ca09ecb
--- /dev/null
+++ b/zUtil_Deformation/gtReadEBSDMapEulerCSVFile.m
@@ -0,0 +1,54 @@
+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
diff --git a/zUtil_Deformation/plotting/gt6DRenderEBSDMapIPF.m b/zUtil_Deformation/plotting/gt6DRenderEBSDMapIPF.m
new file mode 100644
index 00000000..66ace7d3
--- /dev/null
+++ b/zUtil_Deformation/plotting/gt6DRenderEBSDMapIPF.m
@@ -0,0 +1,91 @@
+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
--
GitLab