%%                MTEX Workshop ICOTOM 2017
%                     
%%                  EBSD Analysis with MTEX
%                Ralf Hielscher, TU Chemnitz, Germany
%                David Mainprice, Geosciences Montpellier, France
%
%%
% The following script demostrates some of the facilities of MTEX to
% analyze orientations and misorientations.
%
% Run this script section by section and follow the output. You are
% encauraged to alter the script to get a better feeling about MTEX.
%
% important shortcuts:
% Strg + Return         - run the current section
% Shift + Strg + Return - run the current section and go to the next one
% F9                    - evaluate the selected command
%

%% data import

% import the ebsd data
ebsd = loadEBSD('Emsland_plessite_500x_15_DS.ctf');

% extract crystal symmetries
cs_bcc = ebsd('Fe').CS;
cs_aus = ebsd('Aus').CS;

% recover grains
ebsd = ebsd('indexed');

[grains,ebsd.grainId,ebsd.mis2mean] = calcGrains(ebsd,'angle',5*degree);
ebsd(grains(grains.grainSize<=2)) = [];
[grains,ebsd.grainId,ebsd.mis2mean] = calcGrains(ebsd,'angle',5*degree)

grains = smooth(grains,4);


%% Cluster analysis
% the orientations nicely seperates into different clusters

ori  = grains('fe').meanOrientation;
weights = grains('fe').grainSize;

% 3 clusters
[cid,center] = calcCluster(ori,weights,'halfwidth',10*degree);

% 12 clusters
% [cid,center] = calcCluster(ori,weights,'halfwidth',2.5*degree);

%% select orientations by cluster

% some crystal planes
h = Miller({1,0,0},{1,1,0},{1,1,1},cs_bcc)

% plot orientations of cluster 1
plotPDF(ori(cid == 1),h)

hold on
annotate(center(1))
hold off

%% visualize the clusters

% which clusters to plot
outId = 1:12;

% select nicely seperating colors
C = linspecer(length(outId));

figure(1)
% plot the cluster centers in axis - angle space
plot(center(outId),'axisAngle','MarkerSize',10,'MarkerColor','k')
hold on

% plot clusters in axis angle space
for i = outId
  plot(ori(cid==i),'MarkerEdgeColor',C(i,:))
end
hold off

% plot the clusters in pole figures
figure(2)
for i = outId
  plotPDF(ori(cid==i),h,...
    'MarkerEdgeColor',C(i,:),'MarkerSize',1)
  hold on
end

% the following line can only be evaluated of oriP has been determined
%plotPDF(oriP*KS.variants,'Marker','x','MarkerColor','k','MarkerSize',5)
hold off

%% colorize grains according to cluster

grains_fe = grains('fe');

% which clusters to plot
outId = 1:12;

% select nicely seperating colors
C = linspecer(length(outId));

% plot clusters in axis angle space
for i = outId
  plot(grains_fe(cid==i),'FaceColor',C(i,:),'micronbar','off')  
  hold on
end
plot(grains.boundary)
hold off


%% The Kurdjumov Sachs misorientation

% define Kurdjumov Sachs 
KS = orientation.KurdjumovSachs(cs_aus,cs_bcc);

% and its variants
KSV = KS.variants;

% assume the parent orientation to be in cube orientation
oriP = orientation.id(cs_aus);

hold on
plot(oriP * KSV,'MarkerColor','none','MarkerEdgeColor','k')
hold off


%% Fit Kurdjumov Sachs misorientation

%odf = calcKernelODF(ori,'weights',weights,'halfwidth',10*degree);
odf = calcKernelODF(ori,'weights',weights,'halfwidth',5*degree);


%% First Step: search on a 2.5 degree grid

% generate a 2.5 degree search grid for the parent orientation
S3G = equispacedSO3Grid(cs_aus,'resolution',2.5*degree)

% compute the corresponding child orientations
S3GV = S3G * KSV;

% evaluate the child ODF at KS-child possitions
value = odf.eval(S3GV);

% find the maximum value
[~,id] = max(sum(value,2))

% check 
%plot(ori,odf.eval(ori),'axisAngle')
%hold on
%plot(S3GV(id,:),'MarkerFaceColor','k','MarkerSize',10)
%hold off

% the guessed parent orientation
oriP = S3G(id);

%% check the fitting

% plot the true orientations
plotPDF(ori,h)

% plot the theoretical child orientations
hold on
plot(oriP * KSV,'MarkerColor','none','MarkerEdgeColor','k')
hold off



%% Visualize the result

% rotate the odf such that parent ori would have been the idential
% orientation
odfId = rotate(odf,inv(rotation(oriP)));

plotSection(odfId,'sigma','sections',28)

% visualize alle KS variants
% we should see 72 points
annotate(orientation.id(cs_aus)*KSV)


%% look in misorientation plot

% compute the misorientations
mori = inv(ori) * oriP;

% visulaize in a axis angle plot
plotSection(mori,'axisangle',(43:0.5:47)*degree,...
  'MarkerFaceColor','none','MarkerEdgeColor','b')

% mark orientation relation ships
annotate(orientation.KurdjumovSachs(cs_aus,cs_bcc),'label','KS','backgroundColor','w')
annotate(orientation.GreningerTrojano(cs_aus,cs_bcc),'label','GT','backgroundColor','w')
annotate(orientation.NishiyamaWassermann(ori.CS,ori.CS),'label','NW','backgroundColor','w')
annotate(orientation.Pitch(ori.CS,ori.CS),'label','P','backgroundColor','w')

%% with misorientation distribution function

% compute a misorientation density function from the individual
% orientations
mdf = calcKernelODF(mori,'weights',weights,'halfwidth',2*degree);

% plot it in axis angle section
% this takes some time
plotSection(mdf,'axisangle',(43:0.5:47)*degree,'resolution',1.5*degree)

% mark known orientation relation ships
annotate(orientation.KurdjumovSachs(cs_aus,cs_bcc),'label','KS','backgroundColor','w')
annotate(orientation.GreningerTrojano(cs_aus,cs_bcc),'label','GT','backgroundColor','w')
annotate(orientation.NishiyamaWassermann(ori.CS,ori.CS),'label','NW','backgroundColor','w')
annotate(orientation.Pitch(ori.CS,ori.CS),'label','P','backgroundColor','w')

%% find the "true" orientation relation ship
%

% find the maximum of the misorientation density function
% this takes some time
[~,KS2] = max(mdf)

% and mark it
annotate(KS2,'MarkerFaceColor','red')

%% compare the fitted and the KS orientation relation ship

% plot the rotated odf
plotSection(odfId,'sigma','sections',28)

% visualize alle KS2 variants
% we should see 72 points
annotate(orientation.id(ori.CS)*KS2.variants,'MarkerFaceColor','r')
annotate(orientation.id(ori.CS)*KS.variants,'MarkerFaceColor','k')

%% what about the true parent orientation??

% visualize the remaining austenite orientations
plot(ebsd('Au').orientations,'axisAngle')

% compare with the computed parent orientation
hold on
plot(oriP,'MarkerFaceColor','red','MarkerSize',20)
hold off

% take the ROBUST mean of all austenite orientations
oriP2 = mean(ebsd('Au').orientations,'robust')

% the angle between the "true" and the fitted misorientation relation ship
angle(oriP,oriP2)./degree

%%

% consider Ferrite to Austenite boundaries
gB = grains.boundary('fe','aus')

% and plot them in axis angle section 
plotSection(gB.misorientation,'axisangle',(43:0.5:47)*degree,...
  'MarkerFaceColor','none','MarkerEdgeColor','b')

%% 
% and in 3d space
plot(gB.misorientation,'axisangle')

%%

% compute a misorientation function from the boundary misorientations
% this takes a lot of time :)
mdf2 = calcKernelODF(gB.misorientation,'halfwidth',2*degree)

% visualize the misorientation density function
values = mdf2.eval(gB.misorientation);
plot(gB.misorientation,values)
caxis([min(values),max(values)])

%% investiage orientation relation ships

% this computes planes and directions which are approximately transfered
% into each other by a misorientation
[n1,n2,d1,d2] = round2Miller(KS)

% define the corresponding misorientation
m  = orientation('map',n1,n2,d1,d2)

% lets look how close it is to the initial misorientation
angle(KS,m)

%%
% lets do the same for the fitted orientation relation ship

[n1,n2,d1,d2] = round2Miller(KS2,'penalty',0.001)

% define the corresponding misorientation
m  = orientation('map',n1,n2,d1,d2)

% lets look how close it is to the initial misorientation
angle(KS2,m) ./ degree
