%==========================================================================
%==========================================================================
%
%                                    /$$                                    
%                                   | $$                                    
%                    /$$$$$$/$$$$  /$$$$$$    /$$$$$$  /$$   /$$            
%                   | $$_  $$_  $$|_  $$_/   /$$__  $$|  $$ /$$/            
%                   | $$ \ $$ \ $$  | $$    | $$$$$$$$ \  $$$$/             
%                   | $$ | $$ | $$  | $$ /$$| $$_____/  >$$  $$             
%                   | $$ | $$ | $$  |  $$$$/|  $$$$$$$ /$$/\  $$            
%                   |__/ |__/ |__/   \___/   \_______/|__/  \__/            
%                                                                           
%                                                                           
%                                                                           
%                              /$$                       /$$                
%                             | $$                      | $$                
%                     /$$$$$$ | $$$$$$$   /$$$$$$$  /$$$$$$$                
%                    /$$__  $$| $$__  $$ /$$_____/ /$$__  $$                
%                   | $$$$$$$$| $$  \ $$|  $$$$$$ | $$  | $$                
%                   | $$_____/| $$  | $$ \____  $$| $$  | $$                
%                   |  $$$$$$$| $$$$$$$/ /$$$$$$$/|  $$$$$$$                
%                    \_______/|_______/ |_______/  \_______/                
%                                                                           
%                                                                           
%                                                                           
%                                                    /$$                    
%                                                   |__/                    
%   /$$$$$$  /$$   /$$  /$$$$$$   /$$$$$$   /$$$$$$$ /$$  /$$$$$$$  /$$$$$$ 
%  /$$__  $$|  $$ /$$/ /$$__  $$ /$$__  $$ /$$_____/| $$ /$$_____/ /$$__  $$
% | $$$$$$$$ \  $$$$/ | $$$$$$$$| $$  \__/| $$      | $$|  $$$$$$ | $$$$$$$$
% | $$_____/  >$$  $$ | $$_____/| $$      | $$      | $$ \____  $$| $$_____/
% |  $$$$$$$ /$$/\  $$|  $$$$$$$| $$      |  $$$$$$$| $$ /$$$$$$$/|  $$$$$$$
%  \_______/|__/  \__/ \_______/|__/       \_______/|__/|_______/  \_______/
%                                                                           
%                                                                           
%==========================================================================
%==========================================================================
 
% The aim of this exercise is to get familiar with mtex look at some ebsd
% data

% Notes:
% - usually you can find help by typing `help` command
%   if a command with an identical name might be available for different
%   data types try e.g. `help orientation/calcMDF` or `help odf/calcMDF` or
%   help orientation/mean` (`help mean` will be Matlabs mean for numbers)
% - using TAB-completion makes many tasks faster, prevents typos and helps
%   to explore functions and properties available to an object
%   e.g. ebsd.[hit the tabulator key] will present a list of features
%  available with the ebsd object
%% install/startup mtex (in case it isn't already)
addpath /Path/to/mtex 
startup_mtex

%cd('/Volumes/Lola/digitale bildanalyse/workshops_n_classes/___workshop_chemnitz18/chemnitz18_prep')
%% 0) load some data -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
% we can use the import_wizard and generate a script which imports the ebsd
% data

import_wizard

% selcte EBSD and click Next ...
% when asked for the reference frame, x should be east and y south when
% using convertEuler2SpatialReferenceFrame


%% 1) Explore the ebsd data *-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*

% %########################################################################
% a) Question: What can you tell about the ebsd data?
%    How many indexed pixels ,how many non-indexed?
%    What properties do they have?
%    How can you access the property of a specific pixel?
% %########################################################################
% hint: using ebsd.[TAB] will show you some properties and functions



%       mad = mean angular deviation:misfit between detected bands and
%       simulation (in degree, high mad = 'bad' solution)
%       bc = band contrast: ratio of actual intensity in hough space vs
%       average intensity of kikuchi bands (high value = high contrast)
%       bs =band slope : higest intensity gradient at the morgin of kikuchi
%       bands (high value = sharp band)

%% 2) Look at the map -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*

% You can use the plot command to plot almost anything in mtex
plot(ebsd)
% maybe not so spectacular at first sight, plot(ebsd) without a 
% second argument will only plot the phase

% In general mtex will always honor the convention
% plot(where,what)
% In the case of EBSD data, the "where" are x,y coordiantes of each
% measurement point and the "what" can be anything as long as it has the 
% same size/length as the first argument.

% %########################################################################
% b) Question:
% Colorcode the band contrast in a grey-scale color-coding
% and plot non-indexed pixels in red.
% hint1: when plotting (where,what) values, both need to have the same
% length
% hint2: to overlay two plots, use 'hold on' and 'hold off'
% %########################################################################




%% 3) inspect orientations -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*

% %########################################################################
%  Question:
%  a) access orientations and the first 100 Euler angles
%    note: In mtex, you can access orientations for a single phase only
% %########################################################################




%%
% %########################################################################
% Question:
% b) plot a polefigure (point plot)
%    hint: - define crystal directions (@Miller)
%          - either use plotPDF or remember 
%            that orientiation*Miller = specimen direction
%          - some possibly intersting directions for quartz
%                  [Miller(0,0,0,1,cs,'hkl'),... %c
%                   Miller(1,1,-2,0,cs,'hkl'),... %a
%                   Miller(1,0,-1,0,cs,'hkl')... %m
%                   Miller(0,1,-1,1,cs,'hkl'),... %z
%                   Miller(1,0,-1,1,cs,'hkl')] %r
%
% %########################################################################
% cs = ebsd('q').CS;



%%
% %########################################################################
% Question:
% c) plot an inverse pole figure (point plot)
%    hint: - define a specimen direction (@vector3d)
%          - either use plotIPDF or remember that 
%            orientation*vector3d = crystal direction
%          - possibly interesitng (?) directions for this sample coule be
%            0,0,0 -5,1,0 and 2,-3,0
% %########################################################################



%%
% %########################################################################
% Question:
% d) plot an orientation mapping  
%    Is one mapping sufficient to understand the texture?
%    hint: ipdfHSVOrientationMapping is a nice one!
% %########################################################################



%%
close all
%%
% %########################################################################
% Question:
% e) colorcode points in a polefigure and in an inverse polefigure indentical
%    to the colorcoding in the map
% %########################################################################



%%
close all

%%
% %########################################################################
% Bonus Question
% f1) plot a c-axis mapping
%    hint - the c-axis is just a direction, use .Miller2color and try a
%    and try a mapping for '1' and '-1'
% %########################################################################

%%
% %########################################################################
% Bonus Question
% f2) plot a sigma section, remember that this will give you a c-axis pole
%     figures for different a-axis directions. Does it make sense?
% %########################################################################



%%
% %########################################################################
% Question
% g) Are point plots useful to understand e.g. densities in pdfs or ipdfs?
%    Compare with a controured plot
%    Hint: you plot orientations calculate an odf
%    What kernel size do you want to use?
%    Is it useful to compute if from every pixel?
% %########################################################################





%%
close all
%%
% %########################################################################
% h) Question:
% Do pixels which wereindexed with just few bands (e.g. <=6) show a different 
% polefigure compared to those with more bands (e.g. >9) ?
% %########################################################################




%%
close all
%%
% %########################################################################
% i) Question:
% Compare the <1,1,-2,0> distributions of those orientation which contain 
% c-axes at the periphery of the pole figure to those with c-axes in the center 
% %########################################################################



%% 
% %########################################################################
% j) Question: some more orientation mappings
%  Plot a fibre mapping , color all orientations with <a> up to some angular
%  threshold (e.g. 10*degree) away from the a-maximum
%  hint - find the you can find the maximum using the second output argument
%         of the max function for S2FunHarmonic [~,maxpos]=max()
%       - you can calculate the <a> polefigure of type S2FunHarmonic using 
%         calcDensity(all<a>directions)
%       - fibremapping is called using ipdfCenterOrientationMapping
% %########################################################################




% how many of all orientations are in there? hint fibreVolume


% %########################################################################
% Bonus Question: Why do both methods do not give exactly the same result?
% %########################################################################
%%
% for later use, let's save the ebsd variable
save('ebsd123.mat','ebsd')
%%
 close all

%% -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
% 
%              /$$$$$$$$ /$$$$$$$   /$$$$$$  /$$$$$$$                       
%             | $$_____/| $$__  $$ /$$__  $$| $$__  $$                      
%             | $$      | $$  \ $$| $$  \__/| $$  \ $$                      
%             | $$$$$   | $$$$$$$ |  $$$$$$ | $$  | $$                      
%             | $$__/   | $$__  $$ \____  $$| $$  | $$                      
%             | $$      | $$  \ $$ /$$  \ $$| $$  | $$                      
%             | $$$$$$$$| $$$$$$$/|  $$$$$$/| $$$$$$$/                      
%             |________/|_______/  \______/ |_______/         

%                             /$$                                           
%                            | $$                                           
%                          /$$$$$$$$                                        
%                         |__  $$__/                                        
%                            | $$                                           
%                            |__/          

%         /$$$$$$  /$$$$$$$   /$$$$$$  /$$$$$$ /$$   /$$  /$$$$$$           
%        /$$__  $$| $$__  $$ /$$__  $$|_  $$_/| $$$ | $$ /$$__  $$          
%       | $$  \__/| $$  \ $$| $$  \ $$  | $$  | $$$$| $$| $$  \__/          
%       | $$ /$$$$| $$$$$$$/| $$$$$$$$  | $$  | $$ $$ $$|  $$$$$$           
%       | $$|_  $$| $$__  $$| $$__  $$  | $$  | $$  $$$$ \____  $$          
%       | $$  \ $$| $$  \ $$| $$  | $$  | $$  | $$\  $$$ /$$  \ $$          
%       |  $$$$$$/| $$  | $$| $$  | $$ /$$$$$$| $$ \  $$|  $$$$$$/          
%        \______/ |__/  |__/|__/  |__/|______/|__/  \__/ \______/       

%                                                    /$$                    
%                                                   |__/                    
%   /$$$$$$  /$$   /$$  /$$$$$$   /$$$$$$   /$$$$$$$ /$$  /$$$$$$$  /$$$$$$ 
%  /$$__  $$|  $$ /$$/ /$$__  $$ /$$__  $$ /$$_____/| $$ /$$_____/ /$$__  $$
% | $$$$$$$$ \  $$$$/ | $$$$$$$$| $$  \__/| $$      | $$|  $$$$$$ | $$$$$$$$
% | $$_____/  >$$  $$ | $$_____/| $$      | $$      | $$ \____  $$| $$_____/
% |  $$$$$$$ /$$/\  $$|  $$$$$$$| $$      |  $$$$$$$| $$ /$$$$$$$/|  $$$$$$$
%  \_______/|__/  \__/ \_______/|__/       \_______/|__/|_______/  \_______/
%                                                                           
%                                                                           
%                                                                           
%                                                                      
% The aim of this exercise is to become more familiar with grains in general,
% boundaries, logical indexing and property mapping
%% load the data we previously saved and generate some variable
load('ebsd123.mat')
% in case load did not work, use import_wizard and load the ebsd.ctf
cs=ebsd('q').CS;       

%% 1) make grains -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-
% let's calculate some grains:
[grains,ebsd.grainId,ebsd.mis2mean] = calcGrains(ebsd,'angle',10*degree);
% grains are the grains
% grainId tells every ebsd pixel to which grain it belongs
% mis2mean tells each ebsd pixel it's misorientation to the mean grain orientation
% the 'angle' or 'threshold' is the misorientation angle

% since a short while  it is possible to use any other criterion to segment grains
% e.g. only the c-axis directions - not that it'll be particlualarly
% usefull for some applications but maybe intereting to compare with e.g.
% optical data or consider segmentation procedures which do not depend on
% orientation at all
caxgrains = calcGrains(ebsd('q'),'custom',ebsd('q').orientations*Miller(0,0,0,1,cs),'delta',6*degree)

%% explore the properties of grains: --------------------------------------
% Grain have properties and one in the same way as for ebsd, we can plot
% those e.g.
% %########################################################################
% a) Question:  Plot the frequency distribution of the grain orientation 
%               spread (.GOS) and colorcode each grains acording to the GOS
% %########################################################################



%% -------------------------------------------------------------------------
% %########################################################################
% b) Question: What is the grain size? 
%              plot a histogram of the area equivalent diameter of grains, 
%              what are the frequency weighted mean/mode grain size?
%              hint - to get some area/volume idea use frequency*bincenters*r^2 or
%              r^3 and normalize
% %########################################################################








% So far we did not consider twinning, so keep those numbers in mind
% and go on to see if we were wrong with our gs estimate

%%
%########################################################################
% c) Question: Is there any relation between the GOS and the grain size?
%########################################################################




%% 
% grains have boundaries: and we can access those with grains.boundary
% %########################################################################
% d) Question: How long is the total qtz-qtz boundary in mm
%    hint - inspect grains.boundary('q','q').[TAB]
% %########################################################################




% note these are stepped boundaries and whne considering true length
% one should probably consider the direct path
% compare with smooth boundaries grains_s=smooth(grains('q'),5,'exp');


%% ------------------------------------------------------------------------
% Get rid of twins:
% Before we proceed, we should make sure that we are looking at grains
% and not twins: A very common twin in quartz is the Dauphine twin, 
% which is defined by a 60 degree totation around 0001

% it becomes pretty obvious when plotting e.g. the misorientation angles.
% axis-angle pairs etc. more on that later this week. For now we will just
% look at the angle histogram of grain boundary misorientations
% option 1)
hist(grains('q').boundary('q','q').misorientation.angle./degree,1:110)
% option 2)
plotAngleDistribution(grains('q').boundary('q','q').misorientation, 'resolution',1*degree)

%% Merge grains------------------------------------------------------------
% To get rid of twins e can merge grains which are connected about a twin
% boundary.

%% Option 1:  Merging
% 1) get the grain boundaries
gb_qtz = grains.boundary('q','q');
% 2) define a misorientation assosciated with dauphine twinning
cs=ebsd('q').CS;
twinmis = orientation('axis',Miller(0,0,0,1,cs),'angle',60*degree,cs,cs);
% 3) calculate the angle between the boundary misorientation and the twinning misorientation 
ang_r = angle(gb_qtz('q','q').misorientation,twinmis);
% condition/index all angles below a certain (magic) threshold
ind = ang_r < 6.*degree; % what angle should we choose here?
%select the boundaries based on this condition                                                                        
twinBoundary = gb_qtz(ind);
% merge grains connected with that twin boundary
[grains_merged,parentID_grains]= merge(grains,twinBoundary);
% make a copy of the ebsd dataset
ebsd_merged = ebsd;
% update the grainIds to the parentIds, so each ebsd_merged pixel knows
% again to which grain it belongs to
ebsd_merged.grainId  = parentID_grains(ebsd.grainId);
% plot it
figure
plot(ebsd('q'),ebsd('q').orientations,'faceAlpha',0.15)
hold on
plot(grains_merged.boundary('q','q'),'linewidth',2)
hold on
plot(twinBoundary,'linecolor','r','linewidth',2)

% note that now, grains do not carry any meanOrientation - since it is
% probably nonsense!

%% Option 2:  Temporarily increase symmetry
% just for this particular twinning relation on can also temporarily
% increase the crystal symmetry to hexagonal - the twin dioes not exist - 
% segment grains and set back the symmetry. Now the meanOrientation esists
% and represents something like the modal orientation of the grains

% option 2) change the symmetry
csq= ebsd('q').CS;
cshq = crystalSymmetry('622','mineral', 'Quartz');
ebsd('q').CS=cshq;
[grainsD,ebsd.grainId,ebsd.mis2mean] = calcGrains(ebsd,'angle',10*degree);
ebsd('q').CS=csq;
grainsD('q').CS=csq;
hold on
plot(grainsD.boundary('q','q'),'linecolor','b','linewidth',2);
hold off

% from here, go on with grains_merged / ebsd_merged or grainsD/ebsd
figure
plotAngleDistribution(grainsD('q').boundary('q','q').misorientation, 'resolution',1*degree);


% compare grainsD and grains_merged
grainsD
grains_merged

% %########################################################################
% Question: Do they have the same number of grains and if not, why and what
% could we do about it?
% %########################################################################

%%
close all
%% ------------------------------------------------------------------------
% Now we can plot a more "real" equivalent diameter map and a histogram
plot(grainsD('q'),2.*grainsD('q').equivalentRadius)

bins=0:3:30;
[gs,bc] = hist(grains('q').equivalentRadius.*2,bins)
[gsD,bcD] = hist(grainsD('q').equivalentRadius.*2,bins)

subplot(2,1,1)
bar(bc', [gs' gsD'])
legend(['no Dauphine twins considered: dmean =' num2str(mean(grains('q').equivalentRadius.*2))], ...
       ['Dauphine twins considered: dmean =' num2str(mean(grainsD('q').equivalentRadius.*2))])

subplot(2,1,2)
vol = gs.*bc.^2;
volD = gsD.*bc.^2;
 
bar(bc', [vol/sum(vol); volD/sum(volD)]')
legend(['no Dauphine twins considered'], ...
       ['Dauphine twins considered '])
title('area weighted')
   
   
%% 
close all
%%
% %########################################################################
% e) Question: plot a map of grains which have an aspect ratio larger 2.5

%   for simplicity create a variable just containing quartz grains
%   and use interpolated shapes and discard grain which are too small
%   note: with serious, work-on data, you probably wouldn't want to trust
%   shapes of grains smaller 10,20,50,100 ... pixels, depending on the
%   property one is interested in

gmqtz=grainsD(grainsD.grainSize> 5, 'q')
gmqtz = smooth(gmqtz,5,'exp')

%   hint: use principalComponents to get the parameters for best fit ellipses
[rho,a,b]= gmqtz.principalComponents('hull'); %
rho =  mod(rho./degree,180);
% %########################################################################





%% 
% ########################################################################
% f) Question:  plot only grains with an aspect ratio larger 2.5 AND
%               with a long axis direction +/- 10 degree from ESE-WNW
% %########################################################################

% plot all laxdir with a circular colorcoding




%% 
close all
%% 
% %########################################################################
% g) Question: Compare the texture of grains with high aspect ratios and 
% long axes +/- parallel to ESE-WNW to grains with lower aspect ratios. 
% Also compare  number weighted (1 point per grains) and area weighted
% (1 point per pixel) measures of texture strength.
% hint: to access the ebsd of grains, you can index ebsd by grains, for
% example ebsd(grains([1 7 8])) returns all the ebsd points contrained in
% grains 1,7 and 8.





%%
 close all
% 