%% Import Script for EBSD Data (Tested using MTEX 4)
%**************************************************************************
% Selecting subregions of the EBSD using a polygon of any shape
% or rectangular region (Rubber band method) of any size using your mouse
%**************************************************************************
%
% David Mainprice 14/01/2015
%
% Importing example EBSD dataset used in many examples of the MTEX Help
% file accessible in MATLAB menu Help -> Production Help -> MTEX Toolbox
%
% Reference for EBSD dataset used in this example
%
% Andrea Tommasi, Alain Vauchez, Dmitri A. Ionov (2008)
% Deformation, static recrystallization, and reactive melt transport in shallow
% subcontinental mantle xenoliths (Tok Cenozoic volcanic field, SE Siberia)
% Earth and Planetary Science Letters 272 (2008) 65-77.
%
% Sample
%**************************************************************************
% 9510-16-map-stitch-pscleaned.ctf
%**************************************************************************
% 
% A) Import of Oxford/HKL *.ctf (ASCII) file
%
% Uses *.ctf data file at pname = '/MatLab_Programs/MTEX_Workshop_Distibution';
%
% B) Area fractions of all indexed phases
%
% C) Create simplified "Crystal Symmetry variables" for all phases
%
% D) Create list of pole figures with hkls 100,010,001 for Forsterite
%    Exercise :  create a list for Enstatite and Diopside
%
% E) Calculate the ODF of Forsterite and the Texture Index (= J-index)
%    Exercise :  Vary the hafwidth of ODF from 5 to 20, what effect does
%                it have on texture index and pole figures ?
%                Create an ODF for Enstatite and Diopside and calculate
%
% F) Calculate pole figures of Forsterite
%    Exercise :  Try to plot the polefigures using ODFs of differnt half-width.
%                Try different pole figure options 'contourf', 'antipodal'
%                colorbar % remove colorbars
%                CLim(gcm,'equal');
%                colorbar % add a single colorbar
%
% G) Plot EBSD phase map on orientations (EBSD pixels), this the RAW data.
%
%
% H) Modelling grains using segmentation angle and
%    the 'keepNotIndexed' options
%
% I) Removing small grains - not representative small grains, may be errors
%
%***** MAIN OBJECTIVE OF THIS EXAMPLE Script designed selection ****
%***** of subregions for detailed study *****
% 
% J) Select a region of any shape use the mouse to define polygon
%    Exercise :  Select different areas of your map and plot the pole
%    figures
%    
% K) Select a rectangular region of any size using youe mouse
%    Exercise :  Select different areas of your map and plot the pole
%    figures
%
%% 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('mmm', [4.756 10.207 5.98], 'mineral', 'Forsterite', 'color', 'dark green'),...
  crystalSymmetry('mmm', [18.2406 8.8302 5.1852], 'mineral', 'Enstatite', 'color', 'light red'),...
  crystalSymmetry('12/m1', [9.746 8.99 5.251], [90,105.63,90]*degree, 'X||a*', 'Y||b*', 'Z||c', 'mineral', 'Diopside', 'color', 'dark blue'),...
  crystalSymmetry('m-3m', [5.431 5.431 5.431], 'mineral', 'Silicon', 'color', 'cyan')};

% 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/MTEX_Workshop_Distibution';

% which files to be imported
fname = [pname '/9510-16-map-stitch-pscleaned.ctf'];
%% Import the Data

% create an EBSD variable containing the data
ebsd = loadEBSD(fname,CS,'interface','ctf',...
  'convertSpatial2EulerReferenceFrame', 'wizard')

%%
%***********************************************************************
% Set up area fractions                            (tested with MTEX4)
%***********************************************************************
% list of phase names
Phase_names = ebsd.mineralList;
% zero area fractions
Area_Fractions = zeros(numel(Phase_names));
% number of measurement points in the map
Map_points = length(ebsd);
% default first phase is indexed the perfect sample !
N_first_indexed_phase = 1;
% remove 'nonIndexed' points from total
% if 'notIndexed' is present
if(strcmp(Phase_names(1),'notIndexed'))
    N_first_indexed_phase = 2;
end
%**************************************************************************
% Area fractions for all phases in the EBSD object including non-indexed
%**************************************************************************
% print results to screen
fprintf(' \n')
fprintf(' Area fractions for all phases in the EBSD object \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)));
   if(N_Points > 0)
% 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
end
AF_total = sum(sum(Area_Fractions));
fprintf(' \n')
fprintf(' Total area fraction = %8.4f \n',AF_total)
fprintf(' \n')

% total number of indexed points in map without 'notIndexed' phase
Indexed_map_points = sum(length(ebsd(Phase_names(N_first_indexed_phase:numel(Phase_names)))));
% percentage of indexed point in map
Percent_indexed = 100*Indexed_map_points/Map_points;
% print to screen
fprintf(' \n')
fprintf(' Percentage indexed points in map = %5.1f \n',Percent_indexed)
fprintf(' \n')
%**************************************************************************
% Area fractions based on minerals with EBSD indexed points only
%**************************************************************************
% print results to screen
fprintf(' \n')
fprintf(' Area fractions based on minerals with EBSD indexed points only \n')
fprintf(' \n')
fprintf('   #      Phase       Points  Area fraction \n')
for i=N_first_indexed_phase:numel(Phase_names);
% indexed point for mineral
N_Points = length(ebsd(Phase_names(i)));
   if(N_Points > 0)
% 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
end
AF_total = sum((Area_Fractions(N_first_indexed_phase:numel(Phase_names))));
fprintf(' \n')
fprintf(' Total area fraction = %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)
% indexed point for mineral
     N_Points = length(ebsd(Phase_names(i)));
% print only indexed phases
    if((~strcmpi(Phase_names{i},'notIndexed')) && (N_Points > 0))
% 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_Forsterite = calcODF(ebsd('Forsterite').orientations,'halfwidth',10*degree)
% Texture index
Jindex_Forsterite = textureindex(odf_Forsterite)
%%
%**************************************************************************
% 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' 
% directions are specified by the option 'uvw' 
%
PFs_Forsterite_hkl = [ ...
    Miller(1,0,0,Forsterite_CS,'hkl','phase','Forsterite'), ...
    Miller(0,1,0,Forsterite_CS,'hkl','phase','Forsterite'), ...
    Miller(0,0,1,Forsterite_CS,'hkl','phase','Forsterite'), ...
    ]
PFs_Forsterite_uvw = [ ...
    Miller(1,0,0,Forsterite_CS,'uvw','phase','Forsterite'), ...
    Miller(0,1,0,Forsterite_CS,'uvw','phase','Forsterite'), ...
    Miller(0,0,1,Forsterite_CS,'uvw','phase','Forsterite'), ...
    ]
%%
%**************************************************************************
% Plot Forsterite Pole Figures (PFs)  in Sample co-ordinates (X,Y,Z) (tested with MTEX4)
% Colorbar for each pole figure (default)
%**************************************************************************
figure
setMTEXpref('FontSize',25)
plotPDF(odf_Forsterite,PFs_Forsterite_hkl,'antipodal','resolution',5*degree);
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_Forsterite_colorbar_for_each_default.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,'figSize','large')
%**************************************************************************
%% Grain modelling  (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 : Default in MTEX 4')
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 in now default in MTEX4
if(non_indexed_option == 1)
    [grains,ebsd.grainId,ebsd.mis2mean] = calcGrains(ebsd,'angle',segAngle)
end
% remove non-indexed points, restrict to indexed points only
if(non_indexed_option == 2)
    [grains,ebsd.grainId,ebsd.mis2mean] = calcGrains(ebsd('indexed'),'angle',segAngle)
end
% number of grains all phases 
number_of_calcGrains = grains.length
% plot 'grain' phase map
figure
% Correct MatLab problem with colour buffer
set(gcf,'renderer','zbuffer')
% Default phase map
plot(grains,'figSize','large')
%
%**************************************************************************
%% 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 small grains removed from all phase
number_of_small_grains_removed = number_of_calcGrains - selected_grains.length
%n_Forsterite_grains = selected_grains('Forsterite').length
figure
% Correct MatLab problem with colour buffer
set(gcf,'renderer','zbuffer')
% MTEX 4 default map is phase map
plot(selected_grains,'figSize','large')
%**************************************************************************
%% Select EBSD region of interest with your mouse
%**************************************************************************
close all
figure
plot(selected_grains,'figSize','large') % MTEX4 compatable
% selecting a single grain by x,y coordinates
% Instructions for using the mouse
disp(' ')
disp('Define selected region of map using the curse and mouse')
disp('The last point does not need to exactly the first point')
disp('as the program ensures the polygon is closed')
disp('Left mouse button picks points')
disp('Right mouse button picks last point somewhere near the first point')
disp(' ')
region = selectPolygon;
% select region of polygon with condition
condition = inpolygon(ebsd,region);
% select EBSD data within region and printout to command window
ebsd_subregion = ebsd(condition)
saveFigure('/MatLab_Programs/ebsd_region_selection_polygon_outline_Map.png')
%% plot whole EBSD Orientation map with region of interest marked in red
figure
plot(selected_grains,'figSize','large') % MTEX4 compatable
hold on
plot(region(:,1),region(:,2),'r--','linewidth',2)
hold off
saveFigure('/MatLab_Programs/ebsd_region_polygon_outline_Map.png')
%% plot EBSD Orientation map of region of interest
close all
figure
plot(ebsd_subregion,'figSize','large') % MTEX4 compatable
saveFigure('/MatLab_Programs/ebsd_polygon_subregion_Map.png')
%% reconstruct the grains in the region of interest
grains_subregion = calcGrains(ebsd_subregion,'angle',10*degree)
plot(ebsd,'figSize','large') % MTEX4 compatable
hold on
plot(grains_subregion.boundary,'r--','linewidth',2)
hold off
%%
plot(selected_grains,'figSize','large') % MTEX4 compatable
hold on
plot(grains_subregion.boundary,'color','black')
hold off
%**************************************************************************
%% Create rubberband selection rectangle for EBSD region
%**************************************************************************
close all
figure
plot(selected_grains,'figSize','large') % MTEX4 compatable
% selecting a single grain by x,y coordinates
% Instructions for using the mouse
disp(' ')
disp('Rubberband selection rectangle for EBSD region')
disp('Drag mouse from top right corner with mouse')
disp('button pressed (down), release (up) button to define botton right corner')
disp(' ')
k = waitforbuttonpress; % Wait for key press or mouse-button click
point1 = get(gca,'CurrentPoint');    % button PRESS (down) detected
finalRect = rbbox;                   % return figure units
point2 = get(gca,'CurrentPoint');    % button RELEASE (up) detected
point1 = point1(1,1:2);              % extract x
point2 = point2(1,1:2);              %     and y
p1 = min(point1,point2);             % calculate locations
offset = abs(point1-point2);         % and dimensions
% EBSD map subregion
xy_subregion(1,:) = [p1(1) p1(1)+offset(1) p1(1)+offset(1) p1(1) p1(1)];
xy_subregion(2,:) = [p1(2) p1(2) p1(2)+offset(2) p1(2)+offset(2) p1(2)];
hold on
axis manual
% region = rectangle
% N.B. transpose of xy_subregion
region = xy_subregion';
% draw selection rectangle
plot(region(:,1),region(:,2),'r--','linewidth',2)
hold off
% select region of rectangle with condition
condition = inpolygon(ebsd,region);
% select EBSD data within region and printout to command window
ebsd_subregion = ebsd(condition)
saveFigure('/MatLab_Programs/ebsd_region_selection_rectangular_outline_Map.png')
%% plot whole EBSD Orientation map with region of interest marked in red
close all
figure
plot(selected_grains,'figSize','large') % MTEX4 compatable
hold on
plot(region(:,1),region(:,2),'r--','linewidth',2)
hold off
saveFigure('/MatLab_Programs/ebsd_region_rectangular_outline_Map.png')
%% plot EBSD Orientation map of region of interest
close all
figure
plot(ebsd_subregion,'figSize','large') % MTEX4 compatable
saveFigure('/MatLab_Programs/ebsd_rectangular_subregion_Map.png')