gtEBSDLoadMapEulerCTFFile.m

function EBSD_struct = gtEBSDLoadMapEulerCTFFile(filename, varargin)
% function [EBSD_e_map, EBSD_r_map] = gtEBSDLoadMapEulerCTFFile(filename, varargin)
%   INPUT:
%       axes_perm       : [1:3]   - vector<1x3>
%       stretch_axes    : [1 1]   - vector<1x2>
%       axes_convention : {'a'} | 'b'
%       phase_id        : 1       - scalar
%       rotate_vecs     : [0 0 0] - vector<1x3>
%       permute_maps    : {false} - boolean
%
%   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!
%

    conf = struct( ...
        'axes_perm', 1:3, ...
        'stretch_axes', [1 1], ...
        'phase_id', 1, ...
        'axes_convention', 'a', ...
        'rotate_vecs', [0 0 0], ...
        'permute_maps', false, ...
        'euler_in_degrees', true, ...
        'euler_convention', 'ZXZ' );
    conf = parse_pv_pairs(conf, varargin);

    c = tic();
    fprintf('Reading from file: "%s"..', filename);

    fid = fopen(filename);
    fgetl(fid);
    fgetl(fid);
    fgetl(fid);
    fgetl(fid);
    dims = fscanf(fid, 'XCells\t%d\nYCells\t%d\nXStep\t%f\nYStep\t%f\n');
    fscanf(fid, 'AcqE1\t%d\nAcqE2\t%d\nAcqE3\t%d\n'); % Unknown parameters
    fgetl(fid);
    fgetl(fid);
    fgetl(fid);
    fgetl(fid);
    input = fscanf(fid, '%d\t%f\t%f\t%d\t%d\t%f\t%f\t%f\t%f\t%d\t%d');
    fclose(fid);

    fprintf('\b\b, Done in %f seconds.\n', toc(c));
    c = tic();
    fprintf('Converting to Euler and Rodriguez maps..');

    input = reshape(input, 11, [])';

    scales = dims(3:4);

%     xx = (input(:, 2) ./ scales(1)) + 1;
%     yy = (input(:, 3) ./ scales(2)) + 1;

    ea1 = input(:, 6);
    ea2 = input(:, 7);
    ea3 = input(:, 8);

    sizes = dims(1:2)';
    EBSD_e_map = cat(3, reshape(ea1, sizes), reshape(ea2, sizes), reshape(ea3, sizes));
%     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_mad = reshape(input(:, 9), sizes);
    EBSD_bc = reshape(input(:, 10), sizes);

    if (~conf.euler_in_degrees)
        EBSD_e_map = EBSD_e_map * 180 / pi;
    end

    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);
    switch (upper(conf.euler_convention))
        case 'ZXZ'
            gvdm_EBSD_r_map = gtMathsEuler2Rod(gvdm_EBSD_e_map);
        case 'XYZ'
            gvdm_EBSD_r_map = gtMathsOriMat2Rod(gtMathsEuler2OriMat(gvdm_EBSD_e_map, 'XYZ'));
    end

    fprintf('\b\b, Done in %f seconds.\n', toc(c));
    c = tic();
    fprintf('Applying reference transformations..');

    gvdm_EBSD_r_map = gvdm_EBSD_r_map(abs(conf.axes_perm), :);
    negative = conf.axes_perm < 0;
    gvdm_EBSD_r_map(negative, :) = - gvdm_EBSD_r_map(negative, :);

    if (conf.axes_convention == 'b')
        conv_rotate_vecs = [0 0 tand(30/2)]';
        gvdm_EBSD_r_map = gtMathsRodSum(gvdm_EBSD_r_map, conv_rotate_vecs);
    end
    if (any(conf.rotate_vecs ~= [0 0 0]))
        gvdm_EBSD_r_map = gtMathsRodSum(conf.rotate_vecs', gvdm_EBSD_r_map);
    end

    fprintf('\b\b, Done in %f seconds.\n', toc(c));
    c = tic();
    fprintf('Mapping Rodrigues vectors to fundamental zone..');

    p = gtLoadParameters();
    cryst_system = p.cryst(conf.phase_id).crystal_system;
    cryst_spacegroup = p.cryst(conf.phase_id).spacegroup;
    symm = gtCrystGetSymmetryOperators(cryst_system, cryst_spacegroup);

    gvdm_EBSD_r_map = gtMathsRod2RodInFundZone(gvdm_EBSD_r_map, symm);
    dmvol_EBSD_r_map = gtDefGvdm2Dmvol(gvdm_EBSD_r_map, size_dmvol_EBSD_map);
    EBSD_r_map = reshape(dmvol_EBSD_r_map, [sizes 3]);

    fprintf('\b\b, Done in %f seconds.\n', toc(c));
    c = tic();
    fprintf('Applying requested transformations..');

    EBSD_mask = input(:, 1);
    EBSD_mask = reshape(EBSD_mask, sizes);

    EBSD_mask_inds = find(~EBSD_mask);
    EBSD_mask_inds = cat(1, EBSD_mask_inds, ...
        EBSD_mask_inds + numel(EBSD_mask), ...
        EBSD_mask_inds + 2 * numel(EBSD_mask));
    EBSD_r_map(EBSD_mask_inds) = 0;

    if (any(conf.stretch_axes ~= 1))
        EBSD_e_map = imresize(EBSD_e_map, 'nearest', 'scale', conf.stretch_axes);
        EBSD_r_map = imresize(EBSD_r_map, 'nearest', 'scale', conf.stretch_axes);
        EBSD_mask = imresize(EBSD_mask, 'nearest', 'scale', conf.stretch_axes);
        EBSD_mad = imresize(EBSD_mad, 'nearest', 'scale', conf.stretch_axes);
        EBSD_bc = imresize(EBSD_bc, 'nearest', 'scale', conf.stretch_axes);
    end

    EBSD_mask = logical(EBSD_mask);

    if (conf.permute_maps)
        EBSD_e_map = permute(EBSD_e_map, [2 1 3]);
        EBSD_r_map = permute(EBSD_r_map, [2 1 3]);
        EBSD_mask = EBSD_mask';
        EBSD_mad = EBSD_mad';
        EBSD_bc = EBSD_bc';
    end

    EBSD_struct = conf;
    EBSD_struct.map_e = EBSD_e_map;
    EBSD_struct.map_r = EBSD_r_map;
    EBSD_struct.mask = EBSD_mask;
    EBSD_struct.mad = EBSD_mad;
    EBSD_struct.bc = EBSD_bc;
    EBSD_struct.pixel_size = reshape(scales, 1, []);

    fprintf('\b\b, Done in %f seconds.\n', toc(c));
end