%%                 Chemnitz MTEX Workshop 2018
%                     
%%                  EBSD Analysis with MTEX
%                Ralf Hielscher, TU Chemnitz, Germany
%
%%
% 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.properGroup;
cs_aus = ebsd('Aus').CS.properGroup;

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

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

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

grains = smooth(grains,4);

%%
plotx2east
plot(ebsd('fe'),ebsd('fe').orientations)

hold on
plot(ebsd('Aus'),'facecolor','black','DisplayName','Austenite')
plot(grains.boundary,'linewidth',2)
hold off

%% pole figures

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

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

figure(1)
plotPDF(ori,h)


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

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

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

%%

km = cmeans(3);

km.n = 3;
km.m = 2;
[cid,center] = km.doClustering(ori);


%% select orientations by cluster

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

plotPDF(ori,h)

% plot orientations of cluster 1
plot(ori(cid == 1),'add2all','MarkerFaceColor','r','MarkerSize',3)

hold on
annotate(center(1))
hold off

%% plot the clusters in pole figures

% which clusters to plot
outId = 1:max(cid);

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

figure(1)
opt = {'MarkerSize',3,'MarkerFaceColor','none'};
plotPDF(ori(cid==1),h,...
    'MarkerEdgeColor',C(1,:),opt{:})
for i = outId(2:end)
  plot(ori(cid==i),'MarkerEdgeColor',C(i,:),'add2all',opt{:});
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


%% plot the clusters in axis angle space

figure(2)
clf
% 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



%% colorize grains according to cluster

grains_fe = grains('fe');

% 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 
% parent to child misorientation
KS = orientation.KurdjumovSachs(cs_aus,cs_bcc);
%KS = orientation.NishiyamaWassermann(cs_aus,cs_bcc);

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

% all child orientation for oriP
unique(oriP.symmetrise * inv(KS))

% and its variants
% KSV = unique(symmetrise(oriP) * inv(KS))
KSV = KS.variants;

figure(1)
plot(oriP * inv(KSV),'MarkerEdgeColor','k','add2all',opt{:},'MarkerSize',9)

%% 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 * inv(KSV);

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

% find the maximum value
[~,id] = max(prod(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
plot(oriP * inv(KSV),'MarkerColor','none','MarkerEdgeColor','k','add2all')


%% Visualize the result
% skip this 

% 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)*inv(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')

%% 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.Pitsch(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')

%%

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

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

%% 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

% 
oriP2 = mean(ebsd('Au').orientations,'robust')

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')

%%

mdf2 = calcKernelODF(gB.misorientation,'halfwidth',2*degree)

mori = gB.mis

values = mdf2.eval(gB.misorientation);
plot(gB.misorientation,values)
caxis([min(values),max(values)])

%%
% maybe there is a correlation between grain size and chemistry?

figure
scatter(sqrt(grains.area),grains.fe)


%% extract orientations other values 

ori = ebsd('Fe').orientations;
fe = ebsd('Fe').fe;

hist(fe)

%%

ori = grains('Fe').meanOrientation;
fe = grains('Fe').fe;
area = grains('Fe').area;

hist(grains('fe'))

%% Pole figures in dependency of area or chemistry

% set a condition
cond = area<60;

plotPDF(ori(cond),Miller(0,0,1,ebsd('Fe').CS),'points','all')

% superpose both pole figures with different colors
hold on
plotPDF(ori(~cond),Miller(0,0,1,ebsd('Fe').CS),'points','all')
hold off


