% Tuomo Nyyss?nen 2018
% Scripts to run the algorithm functions without GUI.

%% Import Script for EBSD Data

clear all
close all

%% Specify Crystal and Specimen Symmetries

% crystal symmetry
CS = {... 
  'notIndexed',...
  crystalSymmetry('m-3m', [2.866 2.866 2.866], 'mineral', 'Iron bcc (old)', 'color', 'light blue'),...
  crystalSymmetry('m-3m', [3.6599 3.6599 3.6599], 'mineral', 'Iron fcc', 'color', 'light green')};

% plotting convention
setMTEXpref('xAxisDirection','west');
setMTEXpref('zAxisDirection','outofPlane');

%% Specify File Names

% which files to be imported
fname = ['demo_martensite.cpr'];

%% Import the Data

% create an EBSD variable containing the data
ebsd = loadEBSD(fname,CS,'interface','crc',...
  'convertEuler2SpatialReferenceFrame');

region = [70 45 70 70];
ebsd = ebsd(inpolygon(ebsd,region));
ebsd.id = 1:length(ebsd);

clearvars -except ebsd

%% The previous lines can be generated with the import_wizard feature of MTEX. Now begins the code for running the functions.

%% Create the initial martensite grain map:

%Identify strings for your phases:
string_ferrite = 'Iron bcc';
string_austenite = 'Iron fcc';

%Identify symmetries:
%Iron bcc:
symmetries{1} = ebsd.CSList{2};
%Iron fcc:
symmetries{2} = ebsd.CSList{3};

%%
%Identify grains with 3 degree boundary
[grains,ebsd(string_ferrite).grainId] = calcGrains(ebsd(string_ferrite),'angle',3*degree);

%Remove small grains
ind = grains.grainSize < 4;
ebsd(grains(ind)).phase = 0;
ebsd(grains(ind)).grainId = 0;

%Reidentify grains with small grains removed:
[grains,ebsd.grainId] = calcGrains(ebsd(string_ferrite),'angle',3*degree);

%% Initial map with K-S

% Determine the initial orientation relationship as Kurdjumov-Sachs:
fcc2bcc_KS = (orientation('map',Miller(1,1,1,symmetries{2}),Miller(0,1,1,symmetries{1}),Miller(-1,0,1,symmetries{2}),Miller(-1,-1,1,symmetries{1})));

% Symmetry order correction according to Morito et al.:
cs_symmetries = symmetries{1}.properGroup;

symrots_true(1,1) = cs_symmetries(1);
symrots_true(2,1) = cs_symmetries(21);
symrots_true(3,1) = cs_symmetries(11);
symrots_true(4,1) = cs_symmetries(16);
symrots_true(5,1) = cs_symmetries(24);
symrots_true(6,1) = cs_symmetries(8);
symrots_true(7,1) = cs_symmetries(18);
symrots_true(8,1) = cs_symmetries(10);
symrots_true(9,1) = cs_symmetries(20);
symrots_true(10,1) = cs_symmetries(3);
symrots_true(11,1) = cs_symmetries(7);
symrots_true(12,1) = cs_symmetries(23);
symrots_true(13,1) = cs_symmetries(19);
symrots_true(14,1) = cs_symmetries(2);
symrots_true(15,1) = cs_symmetries(9);
symrots_true(16,1) = cs_symmetries(22);
symrots_true(17,1) = cs_symmetries(17);
symrots_true(18,1) = cs_symmetries(12);
symrots_true(19,1) = cs_symmetries(5);
symrots_true(20,1) = cs_symmetries(15);
symrots_true(21,1) = cs_symmetries(4);
symrots_true(22,1) = cs_symmetries(14);
symrots_true(23,1) = cs_symmetries(6);
symrots_true(24,1) = cs_symmetries(13);

symrots_true = rotation(symrots_true);

% Determine the full orientation relationship according to cubic symmetry
% as by Morito et al.:
fcc2bcc_true = orientation(symrots_true*rotation(fcc2bcc_KS),symmetries{1},1);
l = 1;
    % Determine intermartensitic misorientations:
    bcc_trans =  inv(orientation(rotation(fcc2bcc_true),symmetries{1},1)) * (orientation(rotation(fcc2bcc_true(1)),symmetries{1},1));
    
    % Determine the interferritic grain boundaries from ebsd data:
    bound_bcc = grains.boundary('Iron bcc','Iron bcc');
    
    % Determine all neighboring grain pairs:
    bound_bcc_ids_sorted = sort(bound_bcc.grainId,2);
    [bound_bcc_ids,~,bbcc_ic] = unique(bound_bcc_ids_sorted,'rows');

        % Determine the misorientations between all grain pairs using mean orientations:
    misos_bcc_grains = inv(grains(bound_bcc_ids(:,1)).meanOrientation).*grains(bound_bcc_ids(:,2)).meanOrientation;
    
    % Determine the angles between the possible misorientation set and the
    % grain pair misorientations:
    angles = angle_outer(misos_bcc_grains,bcc_trans) / degree;
    
    % Determine the minimum angle and corresponding index between the
    % misorientation set and each grain pair misorientation:
    [omega, omega_I] = min(angles,[],2);

    [hist_omega_I, bins_omega_I] = hist(omega_I,1:24);
        
    % Determine cutoff as value encompassing 90 % of misorientations
    [hist_omega, bins_omega] = hist(omega,100);
    hist_omega = cumsum(hist_omega)/sum(hist_omega);
    I_hist_omega = find(hist_omega>0.9,1,'first');
    cutoff = bins_omega(I_hist_omega);
    cutoff(cutoff>5) = 5;
    all_cutoffs(l) = cutoff;

    % Limit the list of indexes of the misorientations to those whose minimum angle between the
    % misorientation set is below the cutoff value:
    omega_I_cutoff = omega_I(omega < cutoff);
    
    % Get the misorientations between intermartensitic grain pairs.
    misos_bcc_grains_omega_cutoff = misos_bcc_grains(omega < cutoff);
    
    % Get the grain pairs whose minimum angle between the misorientation set is
    % below the cutoff value (intramartensitic grain pairs).
    OR_rows = bound_bcc_ids(omega < cutoff,:);

 
    
figure
plot(ebsd,ebsd.bc)
colormap gray
hold on;
plot(bound_bcc(ismember(bound_bcc_ids_sorted,OR_rows,'rows')),'linecolor','r','linewidth',2,'faceAlpha',0.5)
plot(bound_bcc(ismember(bound_bcc_ids_sorted,OR_rows_omega_I_cutoff,'rows')),'linecolor','g','linewidth',2,'faceAlpha',0.5)
plot(bound_bcc(not(ismember(bound_bcc_ids_sorted,OR_rows,'rows'))),'linecolor','k','linewidth',3)


%% Iteratively determine the correct orientation relationship and generate figures for checking:

% cond = 0 if you do not want an updating graph, 1 if yes
cond = 1;

[fcc2bcc,omega,bound_bcc_ids,bound_bcc,OR_rows,OR_rows_omega_I_cutoff,bound_bcc_ids_sorted,hist_all_omega_I,symrots_true,bcc2fcc_odf] = determineor(symmetries,grains,cond)

%% Create figure showing band contrast in the background and block (green), packet (red) and PAG (black) boundaries:

figure
plot(ebsd,ebsd.bc)
colormap gray
hold on;
plot(bound_bcc(ismember(bound_bcc_ids_sorted,OR_rows,'rows')),'linecolor','r','linewidth',2,'faceAlpha',0.5)
plot(bound_bcc(ismember(bound_bcc_ids_sorted,OR_rows_omega_I_cutoff,'rows')),'linecolor','g','linewidth',2,'faceAlpha',0.5)
plot(bound_bcc(not(ismember(bound_bcc_ids_sorted,OR_rows,'rows'))),'linecolor','k','linewidth',3)

%% Create a graph describing the grain map and segment it using MCL:

%Determine inflation power:
inflationpower = 1.6;

%Use the Burr survival algorithm to determine 0-1 similarity between each grain pair (based on misorientation angle):
surv_val = 1-cdf('Burr',omega,2,5,1);
%Prune surv_val below certain value
surv_val(surv_val<0.001) = 0;
%Construct the sparse undirected graph required by MCL:
mcl = [bound_bcc_ids, surv_val];

%Create the initial, unsymmetric matrix:
z = sparse(mcl(:,1),mcl(:,2),mcl(:,3));
zz = z;

%Symmetrise the matrix
zzz = [zz;zeros((length(zz) - min(size(zz))),length(zz))];
zzz = zzz' + zzz;

%Put 1s on the diagonal:
zzzz = zzz + speye(length(zzz));
% 
% %Make sure the matrix is indeed sparse:
% zzzzzz = sparse(zzzz);

%Feed the sparse matrix to MCL:
clusters = full(mcl_func(inflationpower,zzzz));

%Change MCL output to row-by-row grain lists (each row is one cluster):
[row_mcl,col_mcl] = find(clusters);
[a,b] = ismember(row_mcl,1:length(grains));
idx = col_mcl;
ib = accumarray(nonzeros(b), idx(a), [], @(x){x});
%Find empty cells
emptyCells = cellfun(@isempty,ib);
%Remove empty cells
ib(emptyCells) = [];

%% Reconstruct the austenite orientation map cluster by cluster:

[ebsd_aus,devis] = recaus(ebsd,grains,symmetries,fcc2bcc,ib);

%% Plot the results

%Plot the martensitic orientations, along with incomplete PAG boundaries:
figure
plot(ebsd(string_ferrite),ebsd(string_ferrite).orientations)
hold on;
plot(bound_bcc(not(ismember(bound_bcc_ids_sorted,OR_rows,'rows'))),'linecolor','k','linewidth',2)

%Plot the austenitic orientations, along with incomplete PAG boundaries:
figure
plot(ebsd_aus(string_austenite),ebsd_aus(string_austenite).orientations)
hold on;
plot(bound_bcc(not(ismember(bound_bcc_ids_sorted,OR_rows,'rows'))),'linecolor','k','linewidth',2)

%% Variant map, entire ebsd data

fccs2bccs = inv(project2FundamentalRegion(ebsd_aus(grains).orientations)).*project2FundamentalRegion(ebsd(grains).orientations);

%Define the fcc to bcc orientation relationship as by Morito et al.:
fcc2bcc_true = orientation(fcc2bcc,symmetries{1},symmetries{2})

%Calculate the angles between the misorientations and the components of the
%full orientation relationship:
variant_angles = angle_outer(orientation(rotation(fccs2bccs),symmetries{1},1),orientation(rotation(fcc2bcc_true),symmetries{1},1)) / degree;

%Determine variant numbers corresponding to best match:
[M_variants, I_variants] = min(variant_angles,[],2);

figure
plot(ebsd(grains),I_variants)
colormap jet(24)
hold on
plot(bound_bcc(not(ismember(bound_bcc_ids_sorted,OR_rows,'rows'))),'linecolor','k','linewidth',3)