% lets's have some data
mtexdata olivine
% Let's take a crop of the data
% This can also be accomplished interactively by:
% plot(ebsd)
% [ebsd,rec] = selectInteractive(ebsd);
ebsd = ebsd(inpolygon(ebsd,[152 296 427 313]))
plot(ebsd)
%% calculating grains
% In mtex, there are currently three different approaches to compute grains
% based on: 
% 1) a misorientation angle threshold (that's probably what most people use)
%    Bachmann et al, 2011, https://doi.org/10.1016/j.ultramic.2011.08.002               
% 2) a threshold of an arbitrary property which is assigned to each pixel
% 3) fast multiscale clustering, FMC 
%    McMahon et al 2013, http://dx.doi.org/10.1016/j.ultramic.2013.04.009
%
% In each case, calcGrains produces 3 outputs:
% - grains : grain object with a lot of different properties
% - ebsd.grainsId :  grain ID to tell each pixel to which grain it belongs
%                    this allows for indexing of the type ebsd(grains)
% - ebsd.mis2mean :  misorientation of each pixel to the mean grain
%                    orientation

%% "default" misorientation angle based segmentation (default is 10 degree)
% The following computes grains from the entire ebsd set
[grains,ebsd.grainId]=ebsd.calcGrains('angle',5*degree)
% plot the grains
plot(grains)

%% calculating grains of a single phase
%  (equivalent to an incomplete list of pixels)
[grains,ebsd('o').grainId]=calcGrains(ebsd('o'),'angle',5*degree)
% plot the grains
plot(grains)
hold on
% and plot the pixels associated with those grains
plot(ebsd('o'),'faceAlpha',0.3,'faceColor','r')
hold off

% Anything notably different?

%% in case filling non-'o' phase-areas is not desired 
% a) use 'unitcell' option
[grains,ebsd('o').grainId]=calcGrains(ebsd('o'),'angle',5*degree,'unitCell');
plot(grains)
hold on
plot(ebsd('o'),'faceAlpha',0.3,'faceColor','r')
hold off
% or 
%% b) assign non-'o' pixels the notIndexed phase
% check available phase options stored in ebsd
ebsd.CSList
ebsdO = ebsd('o')
% a simple way to assign non-indexed pixels
ebsdO = ebsdO.gridify; % this creates an n-by-m matrix, non-existent 
                       % points will be filled with nan
ebsdO(isnan(ebsdO.rotations)).phaseId=1;

% and now calculate grains
[grains,ebsdO.grainId]=calcGrains(ebsdO,'angle',5*degree);
plot(grains)

%% Taking advantage of the grains "ignorance" of non-existent 
% (in contrast to nonIndexed) pixels 
% EXAMPLE-> dealing with single pixel grains

[grains,ebsd.grainId]=calcGrains(ebsd,'angle',5*degree);
plot(grains)
hold on
% plot grains with 1px size
plot(grains(grains.grainSize==1),'faceAlpha',0.7,'faceColor',[0 0 0])
hold off

%% delete ebsd information from 1px grains
ebsd(grains(grains.grainSize ==1)) = [] % <- delete those pixels
% re-construct grains again
[grains,ebsd.grainId]=calcGrains(ebsd,'angle',5*degree);


plot(grains.boundary)
hold on
plot(ebsd('indexed'),'faceAlpha',0.7)
hold off

%% and lets decide that large areas which were previously occupied by 1px
% grains should actually be assigned notIndexed
ebsdg = ebsd.gridify;
ebsdg(isnan(ebsdg.rotations)).phaseId=1;

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

ebsdg(grains(grains.grainSize <= 3)) = [] 

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

% only plot the boundaries
plot(ebsdg,'faceAlpha',0.7)
hold on
% plot grains with 1px size
plot(grains.boundary)
hold off

%% 2) using any arbitrary property
% e.g. iq, c-axis direction, ...

plot(ebsd,ebsd.prop.iq)
mtexColorbar

[grains,ebsd.grainId]=calcGrains(ebsd,'custom',ebsd.prop.iq,'delta',40000);
% some vector colorkey

hold on
plot(grains.boundary)
hold off


%%  effects of 'unitcell', and 'tight'
%% prep: delete grains at the edge
% reconstruct grains
[grains_nob,ebsd.grainId]=ebsd.calcGrains('angle',5*degree);
% find those boundary ids which are 0 ==> outer boundary
outerBoundary_id = any(grains_nob.boundary.grainId==0,2);
% grainId of grains which are at edge
grain_id = grains_nob.boundary(outerBoundary_id).grainId;
% check if selected grains are correct
plot(grains_nob(grain_id(:,2)))
grain_id(grain_id==0) = [] ; % delete column with zeros
grains_nob(grain_id) = [] % delete boundary grains
ebsd_nob = ebsd(grains_nob) % set ebsd according to the remaining grains
% and delete all non-indexed pixels
ebsd_nob('n')=[];
%% effect of 'unitCell, 'tight'
grainsN=ebsd_nob.calcGrains;
grainsU=ebsd_nob.calcGrains('unitCell');
grainsT=ebsd_nob.calcGrains('boundary','tight');
plot(ebsd_nob)
hold on
plot(grainsN.boundary,'linewidth',2,'displayName','default')
plot(grainsU.boundary,'lineColor',[1 0 0],'displayName','unitCell')
plot(grainsT.boundary,'linewidth',2,'lineColor',[0 1 0],'displayName','tight')
hold off
%% above result might only be semi-nice

%% let's replace all non-existent pixels with non-indexed
% and remove nonIndexed "inclusion"
ebsd_nob=ebsd_nob.gridify;
ebsd_nob(isnan(ebsd_nob.prop.iq)).phase=0;
ebsd_nob(isnan(ebsd_nob.prop.iq)).phaseId=1;
ebsd_nob = EBSD(ebsd_nob);

% and find inclusions
[grainsU,ebsd_nob.grainId]=ebsd_nob.calcGrains('unitCell');
isInside = checkInside(grainsU('indexed'),grainsU);
plot(grainsU(any(isInside,2)),'FaceColor',[0 0 0])
hold on
plot(grainsU.boundary,'lineColor',[1 0 0])
plot(ebsd_nob(grainsU(any(isInside,2) & grainsU.phase==0)),'FaceColor',[1 0 1]) % <- those we might want to get rid of
hold off

% delete the isolated pixels of phase non-indexed inside indexed grains
ebsd_nob(grainsU(any(isInside,2) & grainsU.phase==0))=[];
% and delete single pixel, non-indexed grains
ebsd_nob(grainsU(grainsU.grainSize==1),'n')=[];

grainsT=ebsd_nob.calcGrains('boundary','tight');
plot(ebsd_nob)
hold on
plot(ebsd_nob('n'),'FaceColor',[1 1 0])
plot(grainsT.boundary,'linewidth',2)
hold off

%% FMC
% alternatively FMC (fast multiscale clustering) segmentation
% sometimes good for intragranular segmentation
[grains,ebsd.grainId,ebsd.mis2mean]=calcGrains(ebsd,'FMC',2);
% "the larger the parameter, the more grains", where the parameter relates
% to the Mahalanobis distance for clustering, other parameters described in
% McMahon et al., 2013 // http://dx.doi.org/10.1016/j.ultramic.2013.04.009
% are hard coded in @EBSD/private/gbc_FMC.m

plot(ebsd,ebsd.mis2mean.angle/degree)
hold on
plot(grains.boundary,'linewidth',2)
hold off
CLim(gcm,[0, 5])