%
%
%                  /$$$$$$$$ /$$$$$$$   /$$$$$$  /$$$$$$$                       
%                 | $$_____/| $$__  $$ /$$__  $$| $$__  $$                      
%                 | $$      | $$  \ $$| $$  \__/| $$  \ $$                      
%                 | $$$$$   | $$$$$$$ |  $$$$$$ | $$  | $$                      
%                 | $$__/   | $$__  $$ \____  $$| $$  | $$                      
%                 | $$      | $$  \ $$ /$$  \ $$| $$  | $$                      
%                 | $$$$$$$$| $$$$$$$/|  $$$$$$/| $$$$$$$/                      
%                 |________/|_______/  \______/ |_______/                       
% 
% 
% 
%                                                        /$$                    
%                                                       |__/                    
%       /$$$$$$  /$$   /$$  /$$$$$$   /$$$$$$   /$$$$$$$ /$$  /$$$$$$$  /$$$$$$ 
%      /$$__  $$|  $$ /$$/ /$$__  $$ /$$__  $$ /$$_____/| $$ /$$_____/ /$$__  $$
%     | $$$$$$$$ \  $$$$/ | $$$$$$$$| $$  \__/| $$      | $$|  $$$$$$ | $$$$$$$$
%     | $$_____/  >$$  $$ | $$_____/| $$      | $$      | $$ \____  $$| $$_____/
%     |  $$$$$$$ /$$/\  $$|  $$$$$$$| $$      |  $$$$$$$| $$ /$$$$$$$/|  $$$$$$$
%      \_______/|__/  \__/ \_______/|__/       \_______/|__/|_______/  \_______/
% 
% 
% The aim of this exercise is to become a bit more familiar mtex and
% handling of ebsd data

% Quite often, 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)
%
%% 0) load some data -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
% change directory to the place where the sample data is located
% cd workshopfolder

% we will use a cropped part of a larger ebsd map, hence at the moment we
% will not use the import_wizard for now
load ebsd17L.mat

% the crystal symmetry as it would be produced by the import_wizard
% template
% CS = {... 
%   'notIndexed',...
%   crystalSymmetry('-3m1', [4.913 4.913 5.501], 'X||a*', 'Y||b', 'Z||c', 'mineral', 'Quartz', 'color', 'purple')};
% this dataset originates from an hkl/oxford system, hence X||a*,
% originally they set the crystalSymmetry to be -3m and the c-axis as 5.5,
% however it should be 321 and 5.401, we will change this later
% % we can retrieve the crystal symmetry as stored in the dataset from the ebsd 
% CS = ebsd.CSList

% plotting convention
setMTEXpref('xAxisDirection','east');
setMTEXpref('zAxisDirection','intoPlane');

%% 1) explore the ebsd data *-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
% hint: using ebsd.[TAB] will show you some properties and functions
% available for the ebsd object
% a) Which properties does it have, how big is it (in pixel, in micron), which phases are in there ...
        ebsd %will display the ebsd oject
        ext = ebsd.extend
        ext(2)-ext(1)
        ext(4)-ext(3)
        % which is identical to...
        max(ebsd.x)-min(ebsd.x)
        max(ebsd.y)-min(ebsd.y)
        % number of x,y pixels
        (max(ebsd.y)-min(ebsd.y))/(2*ebsd.unitCell(1,1))
        % we seem to a little rounding issue with this dataset
        
%% -------------------------------------------------------------------------
% b) What is the cystal symmetry, which parameters are used?
        ebsd('indexed').CS
% Quartz is trigonal with a symmetry of '321', and cell parameters of a,b= 4.913 c= 5.405
% Oxford/HKL decides decided to use -3m, is that a proper group?   
          
% Can we change the symmetry to a proper group?   

 if  ~ebsd('q').CS.isProper	% if it's not a proper group, make it proper
        ebsd('q').CS=ebsd('q').CS.properGroup
 end      
 
% The c-axis should be 5.401, define the correct cs with a correct correct
% c-length
        cs=crystalSymmetry('321', [4.913 4.913 5.401],...
            'X||a*', 'Y||b', 'Z||c', 'mineral', 'quartz', 'color', 'light blue')
        ebsd('q').CS=cs
        

%% -------------------------------------------------------------------------
% c) inspect properties

% some properties that are in this ebsd object
ebsd.prop

%       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)
%      disregard grainID,mis2mean (we will come to that later)

% so to plot any property
plot(ebsd,ebsd.prop.mad)
% %########################################################################
% Question:
% What is the mean/min/max quartz mad?
% %########################################################################
       mean(ebsd('q').prop.mad)
       % will return nan since there were some non-indexed pixels which were
       % filled but no mad value was generated
       nanmean(ebsd('q').prop.mad)
% %########################################################################       
% Question:       
%   -  What is the percentage of pixels with a less reliable indexing?
        % query points with a specific condition e.g. % of pixels with mad > 1 
        sum(ebsd('indexed').prop.mad >1)/length(ebsd('indexed'))*100
%   -  Is there any relation between mad and band contrast?
        plot(ebsd('indexed').prop.mad,ebsd('indexed').prop.bc,'*','MarkerSize',2)
        ylabel('band contrast'); xlabel('mad')
% %########################################################################        
%% % 2) plot properties in a map -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*

% use the plot command to plot almost anything in mtex
plot(ebsd)
% maybe not so spectacular at first sight
% mtex will (should) always honor the convention
% plot(where,what) where inthsi case are x,y coordiantes of the ebsd map
% %########################################################################
% Question:
% a) Colorcode the number of bands used for indexing of the quartz plase
% in a grey-scale color-coding and non-indexed pixels in red
% hint: when plotting (where,what) values, both need to have the same
% length
% %########################################################################
plot(ebsd('q'),ebsd('q').prop.bands)
mtexColorMap white2black
c = mtexColorbar
c.Label.String='bands'
% xlabel(c,'bands') % on older matlab versions
hold on  % this keeps the present figure axis
plot(ebsd('notIndexed'),'FaceColor','r')
hold off 

%% -------------------------------------------------------------------------
% b) plot the band contrast and plot band contrast lower 50 as black
% hint: you could either plot two maps on top of each other using hold
% on/off, you can set the colorrange using CLim and you can get the
% colormap cm=colormap of the current plot and maniplulate the values
% another hint: use logical indexing: e.g. a=1:9; a(a > 4) returns 5 6 7 8 9
% ebsd data can be querried by logical indexing
plot(ebsd,ebsd.bc)
CLim(gcm,[50  max(ebsd.bc)]) % set the colorraneg for the current axis
mtexColorbar
hold on
ind = ebsd.bc < 50;
plot(ebsd(ind),'FaceColor','k','DisplayName','bc < 50')
hold off


% another option in case min(ebsd.bc)==0 and you liek to overwrite the
% colormap
plot(ebsd,ebsd.bc)
cm = colormap;
id50=round(length(cm)/max(ebsd.bc)*50)
cm((1:id50),:)= repmat([0 0 0],length((1:id50)),1);
colormap(cm)
c = mtexColorbar;
c.Label.String='band contrast'
%% -------------------------------------------------------------------------
% c) select ebsd data based on their position in the map
% for example: define a property and colorcode
%  all quartz points in the dataset as a function of distance from
% an arbitrary line

% hint: distance from a line defined by two points, 
% e.g. http://mathworld.wolfram.com/Point-LineDistance2-Dimensional.html
% using [pq]= ginput(2) will let you select two points on a map
% this data can also be added to the ebsd dataset by 
% ebsd.prop.distancefromline =distancefromline

plot(ebsd)
[pq]= ginput(2)
xy=[ebsd('q').prop.x ebsd('q').prop.y];
p=pq(1,:); q=pq(2,:);
d =(abs((q(2)-p(2)).*xy(:,1) - (q(1)-p(1)).*xy(:,2) + q(1).*p(2) -q(2).*p(1)))./sqrt((q(2)-p(2)).^2 + (q(1)-p(1)).^2);
plot(ebsd('q'),d)
hold on
line(pq(:,1),pq(:,2),'linewidth',2)
hold off
%

% %########################################################################
% Question:
% colorcode all ebsd data points as a function of the distance from of an
% a point of your choice
% hint: choose a point using xy= ginput(1)
% %########################################################################
plot(ebsd)
xy= ginput(1)

% distance from the point
pointdist=sqrt(sum([ebsd.x-xy(1)  ebsd.y-xy(2)].^2,2));
plot(ebsd,pointdist)

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

% a) access orientations and the first 100 Euler angles
ebsd('q').orientations
ebsd('q').orientations(1:100).Euler

% ------------------------------------------------------------------------- 
% b) display orientations as polefigure plot for c (0,0,0,1) and a
%  (1,1,-2,0), learn that orientation*Miller gives a vector in sample
%  reference which can be plotted using plot() (see also help vector3d/plot)
%  and also use plotPDF() (see help orientation/plotPDF)
% 

%define a crystal direction
h = [Miller(0,0,0,1,cs,'uvw') Miller(1,1,-2,0,cs,'uvw')]
% option 1: calculate directions and sue plot - not conventient but you'll
% see how it works

o=ebsd('q').orientations;
% note this will likely plot a lot of points
c_dir=o *h(1)
a_dirs=o*symmetrise(h(2))
mf=mtexFigure;
plot(c_dir,'antipodal','MarkerSize',1,'parent',mf.nextAxis)
mtexTitle('0001')
plot(a_dirs,'antipodal','MarkerSize',1,'MarkerFaceColor','r','parent',mf.nextAxis)
mtexTitle('11-20')
mf.drawNow

% just plot 1000 points
o=discreteSample(o,1000);
c_dir=o *h(1)
a_dirs=o*symmetrise(h(2))
plot(c_dir,'antipodal','MarkerSize',2)
hold on
plot(a_dirs,'antipodal','MarkerSize',2)
hold off


% option 2: use plotPDF - much more convenient
plotPDF(ebsd('q').orientations,h,'points',1000,'markerSize',2)

%% -------------------------------------------------------------------------
% %########################################################################
% c) % Question: plot orientations as sigma sections
% hint: have a look at help plotSection - this helps to study the help
% %########################################################################
plotSection(ebsd('q').orientations,'points',4000,'sections',6,'sigma','MarkerSize',3)
%% -------------------------------------------------------------------------
% d) plot orientations in an inverse pole figure wrt to x,z and a direction 
% oriented NW-SE 
% option 1: learn that orientation\vector or inv(o)*vector gives a Miller
% and option 2, use plotIDPF (see help orientation/plotIPDF)

dir=[xvector zvector vector3d(1,1,0)]
for i=1:length(dir)
m(:,i) =inv(ebsd('q').orientations).*dir(i);
end
m=project2FundamentalRegion(m);

% to plot just 1000 points e.g.
sel=discretesample(1:length(ebsd('q')),1000);
mf=mtexFigure
for i=1:length(dir)
plot(m(sel,i),'FundamentalRegion','MarkerSize',2,'parent',mf.nextAxis)
end

% option 2: use plotIPDF
plotIPDF(ebsd('q').orientations,dir,'points',1000,'MarkerSize',2)
%% ------------------------------------------------------------------------
% %########################################################################
% Question:
% e) do pixels indexed with just few bands (e.g. <=6) show a different 
% pf than those with many (e.g. >9) ?
% 
% Color-code each point according to the mad.
% hint: plotPDF(where, what, Miller)
% %########################################################################

ebsdQ=ebsd('q')
ind6=ebsdQ.prop.bands <= 6;
ind9=ebsdQ.prop.bands > 9;

figure
plotPDF(ebsdQ(ind6).orientations,ebsdQ(ind6).prop.mad,h,'MarkerSize',3)
c = mtexColorbar; ylabel(c,'mad')
CLim(gcm,[min(ebsd('q').prop.mad) max(ebsd('q').prop.mad)])
figure
plotPDF(ebsdQ(ind9).orientations,ebsdQ(ind9).prop.mad,h,'MarkerSize',3)
c = mtexColorbar; ylabel(c,'mad')
CLim(gcm,[min(ebsd('q').prop.mad) max(ebsd('q').prop.mad)])

figure
plotPDF(ebsdQ(ind6).orientations,h,'MarkerSize',3)
hold on
plotPDF(ebsdQ(ind9).orientations,h,'MarkerSize',3)
hold off

%% -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
% 4) orientation/direction/angle mapping 
% -------------------------------------------------------------------------
% a) plot ebsd data in a map with ipdf colorcoding wrt to Z
% and a direction paralle to E-W (X)  (z in our example is the kinematic/strain y-direction,
% and the flow plane is oriented E-W
% plot also a pole figure of c,a,-a,m,r,z where individual points are colorcoded
% with the ipdf color-coding

oM = ipdfHSVOrientationMapping(ebsd('q').CS);
oM.inversePoleFigureDirection = zvector;
plot(oM)
figure;
plot(ebsd('q'),oM.orientation2color(ebsd('q').orientations))
hs= [... % planes
 Miller(0,0,0,1,cs,'hkl'),... %c
 Miller(1,1,-2,0,cs,'hkl'),... %a
 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
 ];
figure
plotPDF(ebsd('q').orientations,oM.orientation2color(ebsd('q').orientations),hs)

oM.inversePoleFigureDirection = vector3d(1,0,0);
figure;
plot(ebsd('q'),oM.orientation2color(ebsd('q').orientations))
figure
plotPDF(ebsd('q').orientations,oM.orientation2color(ebsd('q').orientations),hs)

%% -------------------------------------------------------------------------
% b) plot ebsd data color-coded with respect to a fibre e.g. with 0,0,0,1 
% parallel to the mean c-axis direction in a range of 20 degree
% hint: get the mean direction of the c-axis and use it as inversePoleFigureDirection
% in the orientation mapping
cs=ebsd('q').CS;
h=Miller(0,0,0,1,cs);

c_dir = ebsd('q').orientations.*h;
oM = ipdfCenterOrientationMapping(cs);
oM.inversePoleFigureDirection = mean(c_dir,'antipodal'); %specimendirection
oM.center = h; % fibre
oM.color = [1 0 1];
conewidth = 10*degree;
oM.psi = deLaValeePoussinKernel('halfwidth',conewidth);

figure;
plot(oM)
hold on
circle(oM.center,2*conewidth,'linewidth',2)
hold off
figure;
plot(ebsd('q'),oM.orientation2color(ebsd('q').orientations))

% how many of all orientations are in there?
fibreVolume(ebsd('q').orientations,h,mean(c_dir,'antipodal'),10*degree)
%% -------------------------------------------------------------------------
% c) plot the direction of the c-axis in a map in a "nice" color coding
% hint: the c-axis is just a direction, use .oM.Miller2color
oMs = TSLOrientationMapping(crystalSymmetry('-1'))
plot(oMs,'zAxisDirection','intoPlane','xAxisDirection','east')
% you will be informed that this color-coding is not topologically continuous
c_dir = ebsd('q').orientations.*Miller(0,0,0,1,ebsd('q').CS)


figure
ind= discretesample(1:length(ebsd('q')),500);
plot(c_dir(ind),oMs.Miller2color(c_dir(ind)),'antipodal','MarkerSize',5)

figure
plot(ebsd('q'),oMs.Miller2color(c_dir))

%% -------------------------------------------------------------------------
% %########################################################################
% Question:
% d) plot the angle between flow pane (NW-SW) and an a-axis and plot only those
% within a range < 30 degree
% hint: first we need to get all eq a-axes and we can loop through all
% angles between the a-axes and the reference direction and find the minaml
% value
% %########################################################################
% rotate the vector into crystal coordinates and calc the angle with the
% ref. crystal direction
aang = angle(Miller(1,1,-2,0,ebsd('q').CS),(inv(ebsd('q').orientations).*vector3d(1,1,0)),'antipodal')./degree;

figure
plot(ebsd('q'),aang)
c=mtexColorbar; ylabel(c,'min. angle between flow plane and <a>')
%plot only those within a specific range
figure
ebsdq = ebsd('q');
plot(ebsdq(aang < 30),aang(aang < 30))
c=mtexColorbar; ylabel(c,'angles between <a> and flow plane < 30 degree')
%%
 close all
 
 %% 5) manually select data - ebsd area
%select an ebsd area by clicking around it, then hit enter
plot(ebsd('q'),ebsd('q').orientations)
ebsd_xy = ginput; % click around a polygon
ebsd_sel = ebsd(inpolygon(ebsd,ebsd_xy))
plot(ebsd_sel('q'),ebsd_sel('q').orientations)
%%
 close all
%% -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
% 
%              /$$$$$$$$ /$$$$$$$   /$$$$$$  /$$$$$$$                       
%             | $$_____/| $$__  $$ /$$__  $$| $$__  $$                      
%             | $$      | $$  \ $$| $$  \__/| $$  \ $$                      
%             | $$$$$   | $$$$$$$ |  $$$$$$ | $$  | $$                      
%             | $$__/   | $$__  $$ \____  $$| $$  | $$                      
%             | $$      | $$  \ $$ /$$  \ $$| $$  | $$                      
%             | $$$$$$$$| $$$$$$$/|  $$$$$$/| $$$$$$$/                      
%             |________/|_______/  \______/ |_______/         

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

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

%                                                    /$$                    
%                                                   |__/                    
%   /$$$$$$  /$$   /$$  /$$$$$$   /$$$$$$   /$$$$$$$ /$$  /$$$$$$$  /$$$$$$ 
%  /$$__  $$|  $$ /$$/ /$$__  $$ /$$__  $$ /$$_____/| $$ /$$_____/ /$$__  $$
% | $$$$$$$$ \  $$$$/ | $$$$$$$$| $$  \__/| $$      | $$|  $$$$$$ | $$$$$$$$
% | $$_____/  >$$  $$ | $$_____/| $$      | $$      | $$ \____  $$| $$_____/
% |  $$$$$$$ /$$/\  $$|  $$$$$$$| $$      |  $$$$$$$| $$ /$$$$$$$/|  $$$$$$$
%  \_______/|__/  \__/ \_______/|__/       \_______/|__/|_______/  \_______/
%                                                                           
%                                                                           
%                                                                           
%                                                                      
% The aim of this exercise is to become more familiar with ebsd data and
% grains, access and plot properties of grains

                                                       
%% 5) grains -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*

% let's use cleaned data
% (cleanup procedure was: grains were calculated, non-indexed areas <= 2px have been filled, grains 
% calculated again and dauphine twins were merged, qtz-grains of 1 px have been deleted and filled
% been deleted and filled)
load ebsd17L.mat

% plotting conventions
setMTEXpref('xAxisDirection','east');
setMTEXpref('zAxisDirection','intoPlane');


% make grains inspect grains and access grains properties
[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

%% -------------------------------------------------------------------------
% a) explore the properties of grains:
% plot the frequency distribution of the grain orientation spread (.GOS) and
% colorcode each grains acording to the GOS
hist(grains('q').GOS./degree,20)
figure
plot(grains('q'),grains('q').GOS./degree)
c = mtexColorbar
C.Label.String='GOS'

%% -------------------------------------------------------------------------
% b) plot a histogram of the area equivalent diameter of grains, what is
% the frequency weighted mean/mode grain size? (for simplicty we skip the
% area/volume conversion)
figure
hist(grains('q').equivalentRadius.*2,20)
mean(grains('q').equivalentRadius.*2)

% 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

%% ------------------------------------------------------------------------
% %########################################################################
% Question:
% c) grains have boundaries: How long is the total qtz-qtz boundary in mm
% hint: inspect grains.boundary('q','q').[TAB]
% %########################################################################
sum(grains('q').boundary('q','q').segLength)/1000

%note these are stepped boundaries so maybe not so close to true
% comarpe with smooth boundaries
% hint:
% grains_smoothed=smooth(grains('q'),5,'exp');
sum(grains_smoothed('q').boundary('q','q').segLength)/1000

%% -------------------------------------------------------------------------
% d) plot the misorientation angle distribution and the misorientation axis
% hint: inspect grains.boundary('q','q').[TAB]
% and also consider using plotAngleDistribution
% option 1)
hist(grains('q').boundary('q','q').misorientation.angle./degree,1:110)
% option 2)
plotAngleDistribution(grains('q').boundary('q','q').misorientation)

% in each case, we not that there is a huge peak at 60degree,
% this is due to dauphine twinning, a 60 degree rotation around the c-axis

%% -------------------------------------------------------------------------
% e) get rid of dauphine twins
% hint: you can use merge, however it comes with a new meanorientation of
% grains which might not be any of the twinned portions, hence setting a
% higer symmetry might be an option

% option 1) merge the grains
% get the grain boundaries
gb_qtz = grains.boundary('q','q');
 %define a rotation assosciated with dauphine twinning
rot = rotation('axis',Miller(0,0,0,1,ebsd('q').CS),'angle',60*degree);
% calculates the angle between two trotations - the misorientation and the twinning rotation 
ang_r = angle(gb_qtz('q','q').misorientation,rot);
% condition/index all angles below a certain threshold
ind = ang_r < 4.*degree;
%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.3)
hold on
plot(grains_merged('q').boundary,'linewidth',2)
hold on
plot(twinBoundary,'linecolor','r','linewidth',2)

% option 2) temporarily switch to higher 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


% form here, go on with grains_merged / ebsd_merged or grainsD/ebsd

% !!!!!!!!!
% note that merged grains will NOT have a mean orienation unless you use
% 'calcMeanOrientation' since in many cases this meanorientation will be
% nonsense
% for this particular twinnign realtion, the hexagonal trick produces the
% modal orientatation for each grain with twins

%%
close all
%% ------------------------------------------------------------------------
% e) plot "real" equivalent diameter in a map and in a histogram
plot(grainsD('q'),2.*grainsD('q').equivalentRadius)

bins=0:3:30;
[gs,~] = hist(grains('q').equivalentRadius.*2,bins)
[gsD,~] = hist(grainsD('q').equivalentRadius.*2,bins)
figure
bar(bins', [gs' gsD'])
legend(['no Dauphine twins considered: dmean =' num2str(mean(grains('q').equivalentRadius.*2))], ...
       ['Dauphine twins considered: dmean =' num2str(mean(grainsD('q').equivalentRadius.*2))])

%% -------------------------------------------------------------------------
% f) other built-in parameters:
% plot a map of grains which have an aspect ratio larger 2
% create a new variable of qtz-only ebsd (makes some things easier)
gmqtz = grainsD('q');

%define an index for a condition
condA = gmqtz.aspectRatio > 2;
figure
plot(gmqtz(condA),gmqtz(condA).aspectRatio)
c= mtexColorbar; ylabel(c,'aspect ratio')


%% -------------------------------------------------------------------------
% g) define your own parameters and learn about logical &
% ########################################################################
% Question:
% plot only grains with an aspect ratio larger 2 AND with a long axis
% direction +/- 10 degree from ESE-WNW, plot ipdf standar colorcoding
%get long axis direction
[ev,ew]= gmqtz.principalComponents; %ev being the angle as complex, ew the length of axs
laxdir = atan2(ev(1,2,:),ev(1,1,:)).*-1;
laxdir=laxdir(:);
%dir could go ccw with 0 in the east
laxdir(laxdir<0) = pi + laxdir(laxdir<0);
% %########################################################################

% plot all laxdir with a reasonable colorcoding
plot(gmqtz,laxdir./degree)
mtexColorMap('hsv')
c = mtexColorbar; ylabel(c,'trend of long axis')

% get the angle between vector3d(2,1,0) and laxdir within range of 10
% degree, in our case 2,1,0 is eqal to 153.4 (otherwise one can use polar, take care for the correct direction
[~, azi]= polar(vector3d(2,1,0));
condB = abs(azi+pi-laxdir) <= 10;
figure
plot(gmqtz(condA & condB),gmqtz(condA & condB).meanOrientation)

%% -------------------------------------------------------------------------
% h) ebsd of a subset of grains
% %########################################################################
% Question:
% compare the area weighted pole figures of large and small grains and
% compare the pole figures obtained by grains.meanOrientation
% hint: ebsd(grains,'phase') gives all ebsd data points which belong to
% grains and 'phase'

gmqtz = grainsD('q');
condB = gmqtz.equivalentRadius.*2 < 5;
% %########################################################################


hs= [... % planes
 Miller(0,0,0,1,cs,'hkl'),... %c
 Miller(1,1,-2,0,cs,'hkl'),... %a
 Miller(1,0,-1,1,cs,'hkl')... %r
 ];
figure
plotPDF(ebsd(gmqtz(condB),'q').orientations,hs,'contourf','halfwidth',10*degree,'minmax','colorrange','equal')
mtexColorbar
CLim(gcm,[0 10])
figure
plotPDF(ebsd(gmqtz(~condB),'q').orientations,hs,'contourf','halfwidth',10*degree,'minmax','colorrange','equal')
mtexColorbar
CLim(gcm,[0 10])
figure
plotPDF(gmqtz(condB).meanOrientation,hs,'contourf','halfwidth',10*degree,'minmax','colorrange','equal')
mtexColorbar
CLim(gcm,[0 10])
figure
plotPDF(gmqtz(~condB).meanOrientation ,hs,'contourf','halfwidth',10*degree,'minmax','colorrange','equal')
mtexColorbar
CLim(gcm,[0 10])

%%
 close all
 
 
%% -------------------------------------------------------------------------
% i)  manually select data - grains
%select a grain by clicking at it, when done, hit enter
plot(grainsD('q'),grainsD('q').meanOrientation)
grain_xy= ginput;% close on those grains you want to select, hit ENTER when done
selected_grains = grainsD(findByLocation(grainsD,grain_xy))
plot(ebsd(selected_grains),ebsd(selected_grains).orientations)

%%
 close all

%% 6) misorientations-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*

% have a look at grans.boundary.misorientations.[TAB]

% a) colorcode grain boundaries according to the misorientation angle on
% top of a semi-tranparent map showing grain meanorientation
% 'FaceAlpha' 0-1 should take case of transparency
plot(grainsD('q'),grainsD('q').meanOrientation,'FaceAlpha',0.2)
hold on
plot(grainsD('q').boundary('q','q'),grainsD('q').boundary('q','q').misorientation.angle./degree,'linewidth',2)
hold off
c = mtexColorbar;ylabel(c,'misorientation angle')
%% -------------------------------------------------------------------------
% b) plot and colorcode the misorientation axes along boundaries with respect 
% to crystal ref. frame. Those will be of class Miller, so you can use an
% orinetation mapping with Miller2Color
oMC = TSLOrientationMapping(grainsD('q').CS.Laue);
plot(grainsD('q').boundary('q','q'),oMC.Miller2color(grainsD('q').boundary('q','q').misorientation.axis), 'linewidth',2)
figure
plot(oMC)
%% -------------------------------------------------------------------------
% %########################################################################
% Question:
% c) plot a map with only those boundaries which have misorientation axes 
% close to 10-10
% %########################################################################

gb= grainsD('q').boundary('q','q');
misax = gb.misorientation.axis;
cond = angle(misax,Miller(1,0,-1,0,misax.CS),'antipodal') < 15*degree;
plot(gb(cond),oMC.Miller2color(misax(cond)), 'linewidth',2)
%% -------------------------------------------------------------------------
% d) plot and colorcode the misorientation axes along boundaries with
% respect to sample ref. frame. Those will be of class vector3d
% get orientations along boundary
% hint: you can get the orientation along the boundary from
% ebsd(grains.boundary.ebsdID).orientations
% see also help orientation/axis
% note: misorientation is based in crystal coordinates, and
% misorientation.axis will be in crystalcords.
% orientations are based in specimen coordinates, hence axis(o1,o2) will be
% based in specimen coordinates

% get all orientations from ebsd along boundary
gbo = ebsd(grainsD.boundary('q','q').ebsdId,'q').orientations;    
% get the misorientation axis
misaxS = axis(gbo(:,1),gbo(:,2));
% what is the symmetry of misaxS
oMS = TSLOrientationMapping(crystalSymmetry('-1')); 
    plot(oMS,'zAxisDirection','intoPlane','xAxisDirection','east')
    hold on
    plot(misaxS,'antipodal','contour','contours',0:1:3)
    hold off
figure    
plot(grainsD('q').boundary('q','q'),oMS.Miller2color(misaxS), 'linewidth',2)
%% -------------------------------------------------------------------------
% %########################################################################
% Question:
% e) plot those misorientation axes in a pf/ipf with misorientation angles between 10 and
% 20degree
% %########################################################################

misAng = grainsD('q').boundary('q','q').misorientation.angle./degree;
cond = misAng >=10 & misAng <=20;
figure
plot(misax(cond),'contourf','FundamentalRegion','antipodal','minmax')
figure
plot(misaxS(cond),'antipodal','contourf','halfwidth',10*degree,'minmax')

%%
close all
%% -------------------------------------------------------------------------
% f) comapre c-axis "misorientation" and real misorientation angles

% plot the angles
h=Miller(0,0,0,1,ebsd('q').CS)
gb= grainsD('q').boundary('q','q');
gbo = ebsd(gb.ebsdId,'q').orientations;  

c_axis_angle =   angle(gbo(:,1).*h,gbo(:,2).*h,'antipodal');

plot (c_axis_angle./degree,...
     gb.misorientation.angle./degree, 'o','MarkerSize',1)
xlabel('c-axis angle'); ylabel('misor angle');

% plot the c-axis angle
figure
plot(ebsd('q'),ebsd('q').orientations,'FaceAlpha',0.1)
hold on
plot(gb,c_axis_angle./degree,'linewidth',2)
hold off
c = mtexColorbar;
ylabel(c,'c axis "misor" angle')


% plot the difference betwen c-axis misorientation angles and real
% misorientation angles 
figure
plot(ebsd('q'),ebsd('q').orientations,'FaceAlpha',0.1)
hold on
plot(gb,(gb.misorientation.angle - c_axis_angle)./degree ,'linewidth',2)
hold off
c = mtexColorbar;
ylabel(c,'diff angle')

% but maybe more interesting;
% plot the misorientation angles for those boundaries 
% with c-axis "misorientations" smaller e.g. 3degree
% this can be used when comaring ebsd based segmentations with optically
% derived segmentations

lim=3
ind=c_axis_angle<lim*degree;

figure
plot(ebsd('q'),ebsd('q').orientations,'FaceAlpha',0.1)
hold on
plot(gb(ind),gb(ind).misorientation.angle./degree,'linewidth',2)
hold off
c = mtexColorbar;
ylabel(c,['true misor angle (for c-ax misor angle < ' num2str(lim) ')'])

%%
close all

%% %% -------------------------------------------------------------------------
% g) manually select data along lines -> misor profiles
%select ebsd along a line
plot(ebsd('q'),ebsd('q').orientations)
lineX = ginput(2) % select a line change (2) if your profile should have more than two points
[ebsd_lineX, dist ]= spatialProfile(ebsd('q'),lineX);
%plot pdf with points colorcoded after distance to origin
plotPDF(ebsd_lineX('q').orientations,dist, hs,'points','all','MarkerSize',5)

% ebsd along the line can be used for (mis)orientation profiles, eg
% misorinetation angles
o=ebsd_lineX('q').orientations
m = o(1:end-1).\o(2:end) %misor to neigbour
mo = o(1).\o(2:end) %misor to origin

%
figure
plot(dist(1:end-1),angle(mo)/degree,'linewidth',2)
hold on
plot(dist(1:end-1), angle(m)/degree,'color','r','linewidth',2)
xlabel('distance','FontSize',16); ylabel('misorientation angle in degree','FontSize',16)
legend('to reference orientation','to neighbour')
title(['misor angles (all angles)'],'FontSize',16)
%% -------------------------------------------------------------------------
% %########################################################################
% h) Big-Exercise: Compare misorientations of low angle boundaries (2-10 degree)
% of two types of grains, those with a flat (peripheral) and those with a steep
% (central) c-axis orientation, plot misorientation axes in sample and in 
% crystal reference frame, plot maps of grain boundaries, with a colorcoding
% for misorientation axes in crystal and sample reference frame

% Note: if you want to access boundaries with smaller misorientation angles, 
% you have to segment grains with a smaller threshold and select relevant angles above noise
% ebsd_s= ebsd
% [subgrains,ebsd_S.grainId,ebsd_s.mis2mean] = calcGrains(ebsd,'angle',0.01*degree);
%
% Suggestion: 
% 1) get caxis directio of all grains
% 2) discrimniate two groups of grains based on angle between c-axis direction and a vector
% 3) create 2 ebsd subsets from those grains
% 4) calculate two sets of new_grains with a very small segmentation angle
% 5) caluclate misorientation angle, axes etc from both new subsets of new_grains
% 6) use the misorientation angle to index reasonable values e.g. 2-10
% degree
% 7) plot grains of both c-axes subsets, plot boundary misorientations of new_grains
% ...

% %########################################################################

gmqtz = grainsD('q');
c_dir =gmqtz.meanOrientation.*Miller(0,0,0,1,ebsd('q').CS);
ebsdq =ebsd('q');

% 1) get an index for grains with the c-axis close to the center
cond1 = angle(c_dir, zvector,'antipodal') < 20*degree;
% 2) with the c-axis close to the periphery
cond2 = angle(c_dir, zvector,'antipodal') > 70*degree;

% check selection
figure
plot(gmqtz(cond1),gmqtz(cond1).meanOrientation)
figure
plot(gmqtz(cond2),gmqtz(cond2).meanOrientation)

%get the ebsd of selected  grains
ebsd1 = ebsdq(gmqtz(cond1));
ebsd2 = ebsdq(gmqtz(cond2));

% make new grains with a very small threshold
[subgrains_1, ebsd1.grainId,ebsd1.mis2mean] = calcGrains(ebsd1,'angle',0.01*degree,'unitcell');
[subgrains_2, ebsd2.grainId,ebsd2.mis2mean] = calcGrains(ebsd2,'angle',0.01*degree,'unitcell');


gb1 = subgrains_1.boundary('q','q');
gb2 = subgrains_2.boundary('q','q');


% misorientation angle
misor_angle1 = gb1.misorientation.angle./degree;
misor_angle2 = gb2.misorientation.angle./degree;

% misorientation axis (crytal ref. frame)
misor_ax_crystal1 = gb1.misorientation.axis;
misor_ax_crystal2 = gb2.misorientation.axis;


% misor axes (sample ref. frame)
%get the boundary orientations
gbo1 = ebsd1(gb1.ebsdId).orientations;    
gbo2 = ebsd2(gb2.ebsdId).orientations;
% calc the axes
misor_ax_sample1 = axis(gbo1(:,1),gbo1(:,2));
misor_ax_sample2 = axis(gbo2(:,1),gbo2(:,2));

% we can all the above to a specific range of angles
min_angle = 2;
max_angle = 10;

conda1 = misor_angle1 > min_angle & misor_angle1 < max_angle;
conda2 = misor_angle2 > min_angle & misor_angle2 < max_angle;


% set 1
% make a pleasant colorcoding for misor axes
oMC = TSLOrientationMapping(misor_ax_crystal1.CS,misor_ax_crystal1.CS);
f=newMtexFigure
plot(ebsd1,ebsd1.orientations,'faceAlpha',0.3)
hold on
plot(gb1(conda1),oMC.Miller2color(misor_ax_crystal1(conda1)),'linewidth',2)
hold off

figure
plot(misor_ax_crystal1(conda1),misor_angle1(conda1),'MarkerSize',3,'antipodal','FundamentalRegion')
c = mtexColorbar(gcm);
ylabel(c,'misor angle')

figure
plot(misor_ax_crystal1(conda1),'contourf','antipodal','FundamentalRegion','minmax')
mtexColorbar;


figure
plot(misor_ax_sample1(conda1),oMC.Miller2color(misor_ax_crystal1(conda1)) ,'antipodal','minmax','MarkerSize',2)
hold on
plot(misor_ax_sample1(conda1),'contour','contours',1:1:10,'antipodal','minmax','linewidth',2,'halfwidth',10*degree,'LineColor','k')
hold off



% set 2
figure
plot(ebsd2('q'),ebsd2('q').orientations,'faceAlpha',0.3)
hold on
plot(gb2(conda2),oMC.Miller2color(misor_ax_crystal2(conda2)),'linewidth',2)
hold off


figure
plot(misor_ax_crystal2(conda2),misor_angle2(conda2),'MarkerSize',3,'antipodal','FundamentalRegion')
c = mtexColorbar(gcm);
ylabel(c,'misor angle')

figure
plot(misor_ax_crystal2(conda2),'contourf','antipodal','FundamentalRegion','minmax')
mtexColorbar;

figure
plot(misor_ax_sample2(conda2),oMC.Miller2color(misor_ax_crystal2(conda2)),'antipodal','minmax','MarkerSize',2)
hold on
plot(misor_ax_sample2(conda2),'contour','contours',1:1:10,'antipodal','minmax','linewidth',2,'halfwidth',10*degree,'LineColor','k')
hold off


% 