%% Import Script for EBSD Data (MTEX 4)
%
% David Mainprice 27/10/2014
%
% version simplified for Paper on Spherulites
%
% Sample D85 map 1 (polished surface)
%****************************
%  D85-map1-cl.ctf
%****************************
%
%% clear memory and close plots - for new analysis of EBSD data
clear all
close all
%% Specify Crystal and Specimen Symmetries

% crystal symmetry
CS = {... 
  'notIndexed',...
  crystalSymmetry('12/m1', [9.6809 8.847 5.2169], [90,105.57,90]*degree, 'X||a*', 'Y||b*', 'Z||c', 'mineral', 'Diopside', 'color', 'light blue'),...
  crystalSymmetry('6/m', [9.4202 9.4202 6.9157], 'X||a*', 'Y||b', 'Z||c*', 'mineral', 'Apatite', 'color', 'light green'),...
  crystalSymmetry('m-3m', [8.3958 8.3958 8.3958], 'mineral', 'Magnetite', 'color', 'light red'),...
  crystalSymmetry('12/m1', [8.3482 12.98 7.1582], [90,116.109,90]*degree, 'X||a*', 'Y||b*', 'Z||c', 'mineral', 'Anorthoclase', 'color', 'cyan'),...
  crystalSymmetry('12/m1', [5.357 9.245 20.234], [90,94.98,90]*degree, 'X||a*', 'Y||b*', 'Z||c', 'mineral', 'Biotite - C 2/c', 'color', 'magenta')};

% specimen symmetry : the MTEX4 default is triclinic so do not need this line
SS = specimenSymmetry('triclinic'); 
% plotting convention
setMTEXpref('xAxisDirection','east');
setMTEXpref('zAxisDirection','outOfPlane');

%% Specify File Names

% path to files
pname = '/MatLab_Programs/D85_Rita_Fabio';

% which files to be imported
fname = [pname '/D85-map1-cl.ctf'];

%% Import the Data

% create an EBSD variable containing the data
ebsd = loadEBSD(fname,CS,'interface','ctf',...
  'convertSpatial2EulerReferenceFrame', 'wizard')
%% Edit RGB Red Green Blue values (tested with MTEX4)
% Find the RGB values for mineral colour on EBSD phase map
RGB_Anorthoclase = ebsd('Anorthoclase').color;
% Replace with new RGB values,  the R,G and B Have values between zero and one,
% they can be fraction like 0.25
% Set [R G B] to [0 1 1] R=0 G=1 B=1 this gives 'cyan'
ebsd('Anorthoclase').color = [0 1 1];
%%
%***********************************************************************
% Area fractions of all phases in the ORIGINAL EBSD map (tested with MTEX4)
%***********************************************************************
% phase names
Phase_names = ebsd.mineralList;
Area_Fractions = zeros(numel(Phase_names));
% number of grid points in the map
Map_points = length(ebsd);
% print results to screen
fprintf(' \n')
fprintf(' Area fractions for all phases in the EBSD object \n')
fprintf(' N.B. Includes NON-indexed points \n')
fprintf(' \n')
fprintf('   #      Phase       Points  Area fraction \n')
for i=1:numel(Phase_names);
% indexed point for mineral
N_Points = length(ebsd(Phase_names(i)));
% area fraction for mineral
Area_Fractions(i) = N_Points/Map_points;
% retain first part of mineral name
P_Name = strtok(char(Phase_names{i}),' ');
      fprintf(' %3i %12s %10i %11.4f \n',...
      i,P_Name,N_Points,Area_Fractions(i))
end
AF_total = sum(sum(Area_Fractions));
fprintf(' \n')
fprintf(' Total area fraction of whole map = %8.4f \n',AF_total)
fprintf(' \n')
%
% total number of indexed points in map
Indexed_map_points = sum(length(ebsd(Phase_names(2:numel(Phase_names)))));
% percentage of indexed point in map
Percent_indexed = 100*Indexed_map_points/Map_points;
fprintf(' \n')
fprintf(' Percentage indexed points in map = %5.1f \n',Percent_indexed)
fprintf(' \n')
%
% print results to screen
fprintf(' \n')
fprintf(' Area fractions for INDEXED phases in the EBSD object \n')
fprintf(' \n')
fprintf('   #      Phase       Points  Area fraction \n')
for i=2:numel(Phase_names);
% indexed point for mineral
N_Points = length(ebsd(Phase_names(i)));
% area fraction for mineral
Area_Fractions(i) = N_Points/Indexed_map_points;
% retain first part of mineral name
P_Name = strtok(char(Phase_names{i}),' ');
      fprintf(' %3i %12s %10i %11.4f \n',...
      i,P_Name,N_Points,Area_Fractions(i))
end
% total fraction area of index points
AF_total = sum((Area_Fractions(2:numel(Phase_names))));
%
fprintf(' \n')
fprintf(' Total area fraction of INDEXED map points = %8.4f \n',AF_total)
fprintf(' \n')

%% Generate simple Crystal Symmetry variables for all phases (tested with MTEX4)
fprintf(' \n');
fprintf(' Crystal Symmetry (CS) : names of variables \n');
fprintf(' \n');
% phase names
Phase_names = ebsd.mineralList;
for i=1:length(CS)
    if(~strcmpi(Phase_names{i},'notIndexed'))
% retain first part of mineral name
      P_Name = strtok(char(Phase_names{i}),' ');
% make variables for CS with mineral names
      myvariable = strcat(P_Name,'_CS');
      datavalues = CS{i};
% print variable names for CS
      eval([sprintf(myvariable) ' = CS{i};'])
      fprintf(' %s %s \n','Variable name =',myvariable);
    end
end
disp(' ')
disp(' Use mineral_CS to defined your CS in m-file ')
disp(' ')
%% Calculate an ODF (tested with MTEX4)
odf_Anorthoclase = calcODF(ebsd('Anorthoclase').orientations, 'halfwidth',10*degree)
% Texture index
Jindex_Anorthoclase = textureindex(odf_Anorthoclase)
%%
%**************************************************************************
% Create list of pole figures with hkls 100,010,001     (tested with MTEX4)
%**************************************************************************
% MTEX cannot mix hkls and uvws
% pole to planes can be specified by the option 'pole' or 'hkl' 
%
PFs_Anorthoclase_hkl = [ ...
    Miller(1,0,0,Anorthoclase_CS,'hkl','phase','Anorthoclase'), ...
    Miller(0,1,0,Anorthoclase_CS,'hkl','phase','Anorthoclase'), ...
    Miller(0,0,1,Anorthoclase_CS,'hkl','phase','Anorthoclase'), ...
    ]
PFs_Anorthoclase_uvw = [ ...
    Miller(1,0,0,Anorthoclase_CS,'uvw','phase','Anorthoclase'), ...
    Miller(0,1,0,Anorthoclase_CS,'uvw','phase','Anorthoclase'), ...
    Miller(0,0,1,Anorthoclase_CS,'uvw','phase','Anorthoclase'), ...
    ]
%%
%**************************************************************************
% Plot Anorthoclase Pole Figures (PFs)  in Sample co-ordinates (X,Y,Z) (tested with MTEX4)
%**************************************************************************
figure
setMTEXpref('FontSize',25)
plotPDF(odf_Anorthoclase,PFs_Anorthoclase_hkl,'antipodal','resolution',5*degree,'contourf');
colorbar
% add labels for X,Y and Z sample directions on PFs
annotate([xvector,yvector,-zvector],'label',{'X','Y','Z'},'backgroundcolor','w','MarkerSize',15,'FontSize',20)
% gray for publications
%colormap(white2blackColorMap)
% file name
saveFigure('/MatLab_Programs/PFs_contour_hkls_Anorthoclase_area_wt_D85_map1_cl.pdf')
%close;
%%
%**************************************************************************
% Plot Anorthoclase Pole Figures (PFs)  in Sample co-ordinates (X,Y,Z) (tested with MTEX4)
%**************************************************************************
figure
setMTEXpref('FontSize',25)
plotPDF(odf_Anorthoclase,PFs_Anorthoclase_uvw,'antipodal','resolution',5*degree,'contourf');
colorbar
% add labels for X,Y and Z sample directions on PFs
annotate([xvector,yvector,-zvector],'label',{'X','Y','Z'},'backgroundcolor','w','MarkerSize',20,'FontSize',20)
% gray for publications
%colormap(white2blackColorMap)
% file name
saveFigure('/MatLab_Programs/PFs_contour_uvws_Anorthoclase_area_wt_D85_map1_cl.pdf')
%%
%**************************************************************************
% EBSD Maps
%**************************************************************************
%% Default map is the phase map you can chose the size 'small', 'normal' or 'large' (tested with MTEX4)
% By default MTEX plots a phase map for EBSD data
close all
figure
plot(ebsd,'position','small')
%% plot an EBSD property  - Properties: bands, bc, bs, error, mad, x, y (tested with MTEX4)
% like mean angular deviation (mad)
close all
figure
plot(ebsd,ebsd.mad,'position','small')
colorbar
%% plot an EBSD property - band contrast (bc) (tested with MTEX4)
close all
figure
plot(ebsd,ebsd.bc,'position','large')
% you can select the colormap to suit you requirements
% 'gray' scale from black (0) to white (255)
%colormap gray
% or 'white2black' black (255) to white (0)
mtexColorMap white2black
colorbar
%% orientation map using the rotation angle in specimen coordinates to define the colormap (tested with MTEX4)
% The default reference is reference orientation in Bunge Euler angles (0,0,0) degrees
plot(ebsd('Anorthoclase'),ebsd('Anorthoclase').orientations.angle./degree)
colorbar
%% orientation map using the rotation angle in specimen coordinates to define the colormap (tested with MTEX4)
% User defines reference orientation Bunge Euler angles (17,86,5) degrees
% what else
reference_o = orientation('Euler',17*degree,86*degree,5*degree,ebsd('Anorthoclase').CS)
plot(ebsd('Anorthoclase'),angle(ebsd('Anorthoclase').orientations,reference_o)./degree)
colorbar
%% orientation map one mineral with Euler RGB colormap (tested with MTEX4)
close all
figure
% define the RGB to Bunge Euler angles to RGB orientation colormap
oM = BungeRGBOrientationMapping(ebsd('Anorthoclase'))
plot(ebsd('Anorthoclase'),oM.orientation2color(ebsd('Anorthoclase').orientations))
%% Orientation colour bar plot - euler angle ODF sections (tested with MTEX4)
close all
figure
plot(oM)
%% orientation map one mineral with IPF colorbar for r = xvector (tested with MTEX4)
close all
figure
% select IPF colorbar for specific mineral (a,b,c,alpha,beta,gamma) and symmetry
oM = ipdfHSVOrientationMapping(ebsd('Anorthoclase'))
% define the specimen direction for IPF orientation colormap
oM.inversePoleFigureDirection = xvector;
plot(ebsd('Anorthoclase'),oM.orientation2color(ebsd('Anorthoclase').orientations))
% The 'notindexed' phase are by default plotted as white, you can change it to black 
hold on
plot(ebsd('notIndexed'),'FaceColor','black')
hold off
%% Orientation IPF colourbar plot - monoclinic with Anorthoclase symmetry and r = x = 1,0,0 (tested with MTEX4)
close all
figure
plot(oM)
%% orientation map with colourbar with restricted symmetry (tested with MTEX4)
close all
figure
color = oM.orientation2color(ebsd('Anorthoclase').orientations);
plot(ebsd('Anorthoclase'),color)
%% orientation map with colourbar with all ebsd data (*** does not work bcause multi-phase **)
oM = ipdfHSVOrientationMapping(ebsd);
plot(ebsd,oM.orientation2color(ebsd.orientations))
%% Orientation colourbar with restricted symmetry RGB mix vs number of orientations
close all
figure
plot(color)
%**************************************************************************
%% Correcting EBSD data use mad (tested with MTEX4)
%**************************************************************************
close all
figure
plot(ebsd,ebsd.mad)
colorbar
% plot a histogram
close all
figure
hist(ebsd.mad)
%%
% take only those measurements with MAD smaller then 1.3
ebsd_corrected = ebsd(ebsd.mad<1.3)
close all
figure
hist(ebsd_corrected.mad)
%%
close all
figure
plot(ebsd_corrected,ebsd_corrected.mad)
colorbar
%**************************************************************************
%% Grain analysis  (tested with MTEX4)
%**************************************************************************
%
% You can loop over this cell changing segmentation angle and the
% the 'keepNotIndexed' options
% You should compared the resulting grain map with EBSD (pixel) phase
% or orientation (pixel) map or Band Contrast map
%**************************************************************************
%
%%
disp(' ')
disp(' Grain segmentation angle option ')
disp(' Choose a high angle typically between 15 to 10 degrees for geological samples')
disp(' OR choose low angle of 2 degrees if you want to detect sub-grains')
segmentation_angle = input('The segmentation angle (e.g. 2-15):');
segAngle = segmentation_angle*degree;
%
disp(' ')
disp(' Keep non-indexed points option ')
disp('*1= Scientifically correct, not extrapolating raw indexed data')
disp('    model grains BUT keep non-index points')
disp(' 2= May be more geologically correct in some cases, use with care')
disp('    model grains AND include non-index points within grains boundaries')
disp('    N.B. this option does NOT ADD map pixels with neighbouring ORIENTATIONS')
non_indexed_option = input('Option an integer  (1-2):');
%
% keep non-indexed
if(non_indexed_option == 1)
    [grains,ebsd.grainId,ebsd.mis2mean] = calcGrains(ebsd,'angle',segAngle,'keepNotIndexed')
end
% remove non-indexed points
if(non_indexed_option == 2)
    [grains,ebsd.grainId,ebsd.mis2mean] = calcGrains(ebsd,'angle',segAngle)
end 
% number of orientations (*** does not give number of orientations in
% grains)
n_Anorthoclase_orientations = numel(grains('Anorthoclase'))
% number of grains (*** works***)
n_Anorthoclase_grains = grains('Anorthoclase').length
% plot 'grain' phase map
figure
% Correct MatLab problem with colour buffer
set(gcf,'renderer','zbuffer')
% Default phase map
plot(grains)
%
%**************************************************************************
%% Removing small grains - not representative small grains, may be errors
%**************************************************************************
%
% You can loop over this cell changing number of indexed points per grain
% You should compared the resulting grain map with EBSD (pixel) phase
% or orientation (pixel) map
%**************************************************************************
%
disp(' ')
disp(' Small grains option ')
disp(' Remove small grains containing less than a critical')
disp(' number of indexed points as they error prone or ')
disp(' If you require an accurate grain size and shape analysis')
disp(' the recommended minimum number indexed points per grain size is 10')
disp(' You can decide to keep all grain by accepting all grains with 0')
small_grains_option = input('Indexed points per grain an integer (e.g. 0-10):');
%
% remove grains containing less than critical number of indexed points, 
selected_grains = grains(grains.grainSize > small_grains_option);
% number of orientations  (*** does not work **)
n_Anorthoclase_orientations = length(selected_grains('Anorthoclase'))
% number of grains using numel
n_Anorthoclase_grains = selected_grains('Anorthoclase').length
figure
% Correct MatLab problem with colour buffer
set(gcf,'renderer','zbuffer')
%plot(selected_grains,'property','phase')
% MTEX 4 default map is phase map
plot(selected_grains)
%% Orientation from reference
plot(ebsd('Anorthoclase'),ebsd('Anorthoclase').orientations.angle./degree)
colorbar
hold on
% plot grain boundaries
plot(grains.boundary,'position','large','linewidth',1)
%% Anorthoclase grains coloured by phase (default map) (*** works  **)
figure
plot(grains('Anorthoclase'))
%**************************************************************************
% Size and shape of grains : area,grainSize,perimeter,
% equivalent radius & perimeter,aspectRatio,shapeFactor
%**************************************************************************
%% Anorthoclase grains coloured by area (tested with MTEX4)
plot(grains('Anorthoclase'),grains('Anorthoclase').area)
colorbar
%% Anorthoclase grains coloured by grainSize (tested with MTEX4)
plot(grains('Anorthoclase'),grains('Anorthoclase').grainSize)
colorbar
%% Anorthoclase grains coloured by perimeter (tested with MTEX4)
plot(grains('Anorthoclase'),grains('Anorthoclase').perimeter)
colorbar
%% Anorthoclase grains coloured by equivalentPerimeter (tested with MTEX4)
plot(grains('Anorthoclase'),grains('Anorthoclase').equivalentPerimeter)
colorbar
%% Anorthoclase grains coloured by equivalentRadius (tested with MTEX4)
plot(grains('Anorthoclase'),grains('Anorthoclase').equivalentRadius)
colorbar
%% Anorthoclase grains coloured by aspect ratio (tested with MTEX4)
figure
plot(grains('Anorthoclase'),grains('Anorthoclase').aspectRatio)
colorbar
%% Anorthoclase grains coloured by shape factor (tested with MTEX4)
figure
plot(grains('Anorthoclase'),grains('Anorthoclase').shapeFactor)
colorbar
%**************************************************************************
% Orientation, misorientation of grains : meanOrientation,KAM,GOS,Angle,Axis
%**************************************************************************
%% Anorthoclase grains coloured by mean grain orientation (tested with MTEX4)
% default specimen direction (r = z = 0,0,1) for IPF orientation colormap
color = oM.orientation2color(ebsd('Anorthoclase').orientations);
plot(grains('Anorthoclase'),grains('Anorthoclase').meanOrientation,color)
colorbar
%% Anorthoclase grains coloured by mean grain orientation (tested with MTEX4)
% select IPF colorbar for specific mineral (a,b,c,alpha,beta,gamma) and symmetry
oM = ipdfHSVOrientationMapping(ebsd('Anorthoclase'))
% define the specimen direction (r = x = 1,0,0) for IPF orientation colormap
oM.inversePoleFigureDirection = xvector;
plot(grains('Anorthoclase'),oM.orientation2color(grains('Anorthoclase').meanOrientation))
%% Anorthoclase grains coloured by Grain Orienation Spread (GOS) (tested with MTEX4)
plot(grains('Anorthoclase'),grains('Anorthoclase').GOS)
colorbar
%% Anorthoclase grains coloured by Kernel Average Misorientation (KAM)(tested with MTEX4)
% Note KAM is defined using EBSD data, it is the kernel average
% misorientation for each EBSD measurement
% options : 'threshold',10*degree and 'secondOrder'
% Kernels for square EBSD grids
% $ = center point *= points used for Kernel
% 1st order KAM (default)    2nd order KAM (option 'secondOrder')
%                                      *
%           *                         ***
%          *$*                       **$**
%           *                         ***
%                                      *
plot(ebsd('Anorthoclase'),KAM(ebsd('Anorthoclase'),'secondOrder','threshold',15*degree)./degree)
colorbar
%% Rotation ANGLE of mean orientation of each grain from reference orientation (tested with MTEX4)
% Uniform colour for each grain mean orientation
% get all rotations of mean rotation from the reference orientation
% using all ebsd indexed points inside grains
Mean_rotations = grains('Anorthoclase').meanRotation;
figure
plot(grains('Anorthoclase'),angle(Mean_rotations)/degree)
% colorbar
colorbar
%% Misorientation ANGLE from mean within a grain (tested with MTEX4)
% get all misorientations from the mean orientation
% using all ebsd indexed points inside grains
% error message non-existent field 'mis2meanRotation' if you have not used
% new syntax [grains ebsd] = calcGrains(ebsd)
figure
plot(ebsd('Anorthoclase'),ebsd('Anorthoclase').mis2mean.angle./degree)
% caxis - allows you focus the scale of the angle range of interest
% example 0 to 5 degrees
caxis([0 5])
%% Misorientation AXIS from mean within a grain (tested with MTEX4 scale 0-1 is this normal ?)
% get all misorientations from the mean orientation
% using all ebsd indexed points inside grains
% colourmap
oM = ipdfHSVOrientationMapping(ebsd('Anorthoclase'))
% plot
figure
plot(ebsd('Anorthoclase'),oM.Miller2color(ebsd('Anorthoclase').mis2mean.axis))
colorbar
% caxis - allows you focus the scale of the range of interest
% does not work with this plot, must be normalied colour scale ...
%caxis([0 30])

