Abstract
We propose a framework for indexing of grain and subgrain structures in electron backscatter diffraction (EBSD) images of polycrystalline materials. The framework is based on a previously introduced physicsbased forward model by Callahan and De Graef (2013) relating measured patterns to grain orientations (Euler angle). The forward model is tuned to the microscope and the sample symmetry group. We discretize the domain of the forward model onto a dense grid of Euler angles and for each measured pattern we identify the most similar patterns in the dictionary. These patterns are used to identify boundaries, detect anomalies, and index crystal orientations. The statistical distribution of these closest matches is used in an unsupervised binary decision tree (DT) classifier to identify grain boundaries and anomalous regions. The DT classifies a pattern as an anomaly if it has an abnormally low similarity to any pattern in the dictionary. It classifies a pixel as being near a grain boundary if the highly ranked patterns in the dictionary differ significantly over the pixel’s neighborhood. Indexing is accomplished by computing the mean orientation of the closest dictionary matches to each pattern. The mean orientation is estimated using a maximum likelihood approach that models the orientation distribution as a mixture of Von Mises–Fisher distributions over the quaternionic sphere. The proposed dictionary matching approach permits segmentation, anomaly detection, and indexing to be performed in a unified manner with the additional benefit of uncertainty quantification. We demonstrate the proposed dictionarybased approach on a Nibase IN100 alloy. ^{c1}^{c1}c1 Part of this work was reported in the Proceedings of the IEEE International Conference on Image Processing (ICIP), Melbourne Australia, Sept 2013.
Keywords
electron backscatter diffraction pattern;
EBSD;
dynamical electron scattering;
dictionary matching;
maximum likelihood orientation estimates;
Von Mises–Fisher mixture distribution
1 Introduction
Electron backscatter diffraction, EBSD, is used to perform quantitative microstructure analysis of polycrystalline materials on a millimeter to nanometer scale (Schwartz et al. (2009)). Most current EBSD segmentation (i.e., delineation of individual grains by determination of the grain boundary locations) and indexing (i.e., orientation determination) methods extract orientations and widths of Kikuchi bands in a measured pattern by using a modified Hough transform, implemented using image processing tools such as butterfly convolution, Gaussian filtering, binning, peak detection, and image quality maps to gauge indexing and segmentation accuracy (Tao and Eades (2005)), (Wright and Nowell (2006)). By comparing the measured diffraction line parameters to a precomputed database, indexing yields the crystal orientation, commonly described by three Euler angles with respect to a reference frame, for the volume illuminated by the beam. By repeating the process on a grid of scanning locations on the sample, an orientation map or image is produced. The image is then segmented into grains by thresholding normed differences between the Euler angles (misorientations). The accuracy of the Hough approach to EBSD indexing depends to a large extent on the visibility of the Kikuchi bands, which is often represented in terms of an ”image quality” parameter.
In this paper we propose an alternative indexing approach that uses a physicsbased forward model for the full diffraction patterns and does not require application of the Hough transform or other image processing tools. The proposed approach exploits the known physics of electron scattering phenomena that underlies the image formation process. With this additional information, grain boundaries and anomalous points can be detected as explicit classes at the same time as grains are segmented. Anomaly detection (i.e., the automated detection of abnormal or unexpected diffraction patterns) is an important capability, since anomalies may correspond to defects or contaminants that affect the material properties. In addition, unlike methods based on the Hough transform, the proposed pattern dictionary approach differentiates grain interiors from grain boundaries without requiring additional processing of the measured patterns.
Automated indexing of electron diffraction patterns by means of pattern matching techniques is not new, and was proposed in 1991 by Wright and coworkers (
Wright et al. (1991)) for backscattered Kikuchi diffraction patterns, now commonly known as electron backscatter diffraction patterns or EBSPs. Their approach involved automated comparisons between a series of experimental patterns and idealized patterns, created from a set of orientations that uniformly covers the asymmetric region (or fundamental zone) of Euler space. While the technique showed promising results, computer limitations prevented the approach from gaining widespread acceptance. In the context of precession electron diffraction in the transmission electron microscope, Rauch and coworkers (Rauch and Dupuy (2005); Rauch et al. (2008)) proposed a templatebased pattern matching approach for the automatic orientation determination of quasikinematical diffraction patterns. The templates are computed from kinematical structure factors, and the proper asymmetric part of Euler space is sampled uniformly to generate a template collection, which is then compared quantitatively against the experimental patterns. In this template matching process, only the top match is considered in the determination of the crystal orientation. As we will explain in detail in the present paper, our dictionary appproach uses, in addition to a physicsbased model for the generation of the dictionary, a statistical model for the orientation distribution of all the highly ranked pattern matches to provide both a stable and robust estimate of orientation, as well as a quantitative statistical uncertainty; the method proposed in
Rauch and Dupuy (2005) does not perform such a statistical analysis. The current commercially available EBSD indexing suites also lack a statistical determination of the uncertainties in the orientation determination.The proposed indexing framework relies on two components: offline dictionary generation and online dictionary matching to the sample. Offline dictionary generation is accomplished as follows. First, a dictionary of raw diffraction patterns is generated using the forward model of Callahan and De Graef (2013)
. This dictionary is tuned to the parameters of the microscope and the crystal symmetry group(s) of the sample. Second, the singular value decomposition (SVD) of the dictionary is computed and a second dictionary, called the dictionary of background compensated patterns, is generated by projection of the raw dictionary onto the space orthogonal to the first principal component; this is essentially a background subtraction process.
The first step in the online dictionary matching algorithm is to compute normalized inner products between the uncompensated (raw) sample patterns and the uncompensated dictionary. This is repeated on the compensated sample patterns and the compensated dictionary patterns. The basis for the online dictionary matching algorithm is the construction of a pair of bipartite graphs, uncompensated and compensated, respectively, using these normalized inner products. A bipartite graph is a graph connecting vertices or nodes in two disjoint sets (Diestel (2005)); in our case, the first set contains the experimental patterns at different locations on the sample, the second the dictionary patterns for different orientations. For each spatial location on the sample, the top (normalized) inner products between the sample diffraction pattern and a dictionary determine the
nearestneighbor (kNN) neighborhood of the pattern in the dictionary. These kNN neighborhoods form the backbone of the proposed online dictionary matching algorithm.
From the bipartite graphs, sample patterns can be classified as grain interiors, grain boundaries, or anomalies, using an unsupervised decision tree (DT) classifier defined on the graphs. Specifically, the classifier uses the shapes of these kNN neighborhoods to discriminate between these types of patterns. A tightly clustered and connected kNN neighborhood indicates a grain interior sample pattern. A kNN neighborhood that forms two or three clusters in the dictionary indicates a grain boundary sample pattern. kNN neighborhoods that are very spreadout and scattered indicate anomalous sample patterns. The unsupervised DT classifier uses the uncompensated dictionary to distinguish unusual (anomalous) background patterns. The compensated dictionary is used to distinguish nonanomalous patterns as interior to a grain or on the boundary of a grain. A pixel is classified as ”anomalous”, ”grain interior” or ”grain boundary” using an unsupervised decision tree (DT) classifier to test homogeneity of the pattern matches over a spatial patch centered at the pixel. The effect of surface roughness, and the resulting shifts in the EBSD background intensity, is a topic of ongoing analysis.
Indexing of crystal orientation is performed using a maximum likelihood (ML) estimation strategy for determining the Euler angles. The ML strategy fits a Von MisesFisher mixture (VMFm) density model to the observed distribution of the Euler angles of the top dictionary matches; the von MisesFisher distribution is a probability distribution on a sphere in a
dimensional space. An iterative expectationmaximization (EM) estimation algorithm is used to perform the fit. The use of the VMFm model allows one to account for the symmetryinduced ambiguities of the crystal orientation and produces an estimate of the Euler angles in the desired fundamental zone. A side benefit is that the algorithm yields an estimate of the spread of the VMFm model that can be used as an a posteriori confidence measure on the orientation estimate.
To the best of our knowledge, the framework proposed here is the first EBSD indexing approach that uses a dictionary generated by a physicsbased forward model. Some advantages of our modelbased approach are: 1) it incorporates the physics of dynamical electron scattering; 2) it unifies segmentation, indexing and anomaly detection; 3) it incorporates a statistical model that naturally generates both an estimate of orientation and a measure of confidence in the estimate; 4) it involves parallelizable operations relying on simple inner products and nearest neighbor search. At the same time, the large size of the dictionary, together with the high dimension of the diffraction patterns, create computational challenges as discussed in Sec. 5.
The outline of the chapter is as follows. Section 2 presents the proposed dictionary model for EBSD pattern classification and indexing. Section 3 develops the statistical algorithms for anomaly detection, segmentation, and indexing based on the dictionary model. Section 4 describes the dictionary generation process. Section 6 presents the experimental methods for generating the EBSD samples used in Sec. 7. Section 7 presents the results of applying the proposed dictionarybased classification and indexing to a Nibase IN100 alloy. Finally, Section 8 summarizes the paper and points to future directions.
2 Dictionary Model
This section describes the nonlinear forward model of Callahan and De Graef (2013). It then shows how a physicsbased dictionary of diffraction patterns can be used to approximately linearize the forward model. For descriptive economy, throughout this paper we denote by a ”pixel” a particular scan location on the sample surface; with each pixel, there is an associated EBSD pattern acquired at the sample location.
2.1 Forward model
When the beam is focused on a grain within the sample, the measured backscatter diffraction pattern, , on the detector surface can be expressed as a function of the crystal orientation , parameterized by an Euler angle triplet , the incident electron energy, , and interaction depth, , as:
(1) 
where is a forward model for the backscatter process and is detector noise. Here is a measurement operator that accounts for the instrument geometry and sensitivity, and represents the thickness and energy averaged mean backscatter yield. An accurate forward model for this yield was proposed by Callahan and De Graef (2013). The model employs Monte Carlo simulations to obtain the energy, spatial, and exit depth distributions of the backscattered electrons. This statistical information is then used to compute a series of dynamical ”EBSD master patterns” as a function of the electron exit energies and directions. The master pattern represents all possible EBSD patterns for a given exit energy, and specification of the grain orientation along with the detector geometry then leads to an actual EBSD pattern for that orientation. These three steps constitute the forward model . The measurement operator includes the scintillatortoCCD conversion process in the form of a point spread function for the coupling optics, as well as detector quantum efficiency and CCD binning mode. In the remainder of this paper, we will refer to the depth and energy averaged diffraction pattern generated by the forward model described above as the mean diffraction pattern.
As the electron beam is scanned across the sample surface, the diffraction pattern will change due to changes in the local crystal orientation between homogeneous grains. In grain interiors, the forward model (1) can be used to produce estimates of crystal orientation at each scan location. Elsewhere in the sample, e.g., near boundary regions or near locations of anomalous features, the model (1) will no longer be a good fit to the measured patterns. Thus, the goodness of fit of the forward model can be used to classify grains, grain boundaries, and anomalies in the sample. This forms the basis for our proposed use of the forwardmodel to perform classification and indexing.
2.2 Sparse dictionarybased forward model
In principle one could formulate classification and indexing as a nonlinear inverse problem using the full forward model (1). For example, given a noise model for one could perform maximum likelihood estimation to solve the indexing problem and likelihoodratio testing to solve the classification problem. In the special case of a Gaussian noise model both solutions would require solving the nonlinear leastsquares problem , in which the Euclidean norm squared of the residual fitting errors is to be minimized with respect to the orientation . Here, we take a simpler approach to the inverse problem that leads directly to tractable indexing and classification algorithms.
Let denote a precomputed dictionary of mean diffraction patterns obtained by densely sampling the function over the range of orientations , keeping and fixed. Assume that the size of is . Then, for sufficiently dense samples, the model (1) can be approximated by
(2) 
where are dictionary elements and are coefficients. When corresponds to the measured pattern at a location within a grain, one might expect that only a few ’s will be nonzero, i.e., the representation (2) is sparse. In particular, as the dictionary becomes increasingly dense the sparsity of the representation will also increase and, in the limit, will become a delta function ^{c1}^{c1}c1The Delta function is equal to if and is equal to for other ., where is the index of the true orientation at that location. Note that when , in this limiting case (of a fully dense dictionary) errorless estimation of can be accomplished by finding the index which yields the largest normalized inner product , where, for two diffraction pattern and , , where and index the vertical and horizontal locations on the photodetector. The significance of this fact is that we have simplified the solution of a complicated nonlinear least squares problem to the solution of a linear least squares problem followed by a table lookup (matching an index of the dictionary to the associated Euler angle).
In the practical case of a finite dictionary there will not be an exact match to the true diffraction pattern and (2
) is interpreted as a model that interpolates over the patterns in the dictionary, with interpolation coefficients
. For the purposes of indexing and classification we will restrict ourselves to sparse models, where only a few () of these coefficients are nonzero, i.e., using only a small number of dictionary elements to fit (2). This sparse approximation problem is a well studied mathematical problem with many different iterative algorithms available for identifying the few nonzero coefficients. The brute force algorithm that tries to find the best fit over any set of dictionary elements is intractable except for very small values of . Alternatives include basis pursuit methods such as orthogonal matching pursuit (OMP), stepwise OMP (stOMP), compressive sampling OMP, iterative soft thresholding (IST), and convex optimization relaxation methods such as minimization using active set, interior point, or subgradient methods (Tropp and Wright (2010)). In the EBSD application, the sparse approximation has to be performed for each and every pattern measured on the sample. Even with a relatively modest dictionary size and low sample resolution, e.g., elements and scan locations, these methods are computationally heavy.We have adopted a simpler correlation matching approach with significantly reduced computation requirements. Instead of fitting the model (2) through least squares we simply use inner products to find the top matches between the dictionary and the observations. Specifically, for each measured pattern we compute the normalized inner products (correlation) between and the dictionary and rank them in order of decreasing magnitude. The top correlation matches are the dictionary patterns having the highest inner products^{c1}^{c1}c1Note that these top inner product matches will be identical to the dictionary elements selected by iterations of OMP in the case that the dictionary is orthogonal.. These patterns constitute and pixel’s nearestneighbor (NN) neighborhood in the dictionary. The collections of NN’s define a bipartite graph connecting measured patterns to patterns in the dictionary. In general, these connections will be onetomany, i.e., for each experimental pattern there will be a small number of near matches in the dictionary. The NN neighborhoods will be used to perform classification and indexing as described in Sec. 3.
In addition to the dictionary , referred to as the uncompensated dictionary, we will use a derived dictionary of compensated patterns to cluster the observed patterns into classes that can then be used for anomaly detection, segmentation and indexing. The compensated dictionary consists of patterns in after projecting away the background. This is performed by applying the singular value decomposition to determine the principal component, which closely resembles the population mean of the patterns in , and projecting the dictionary onto the space orthogonal to the principal component (i.e., removing the mean pattern from each individual pattern). This compensation process, which simply removes the background common to all patterns, improves the dictionary’s ability to discriminate between diffraction patterns. However, the uncompensated dictionary will also be used since it can better discriminate anomalies that are primarily manifested in the background.
3 Classification and Indexing
The proposed correlation matching approach to classification procedes in two steps. First inner products between the observed patterns and patterns in the dictionaries and are computed. For each observed pattern we only store its closest dictionary matches, i.e., those dictionary patterns with the highest inner products with respect to the observed pattern. Then a pixel is classified as a grain interior, a grain boundary, or an anomaly using a decision tree classifier applied to the set of top pattern matches in the dictionary. The proposed correlation matching approach to indexing pixel orientation computes an estimate of the mean orientation over the top matching patterns in the dictionary. In order to account for noise and the noneuclidean nature of the Euler sphere, the mean angles are estimated using a specially adapted maximum likelihood estimator, introduced in the indexing subsection below. Pixel classification and indexing are performed independently and are discussed separately.
3.1 Classification
Classification of a pixel is performed by evaluating the inner products between the pixel’s uncompensated and compensated patterns and patterns stored in the dictionaries and , respectively. Anomalies are detected as abnormally low average inner products between an uncompensated pixel pattern and patterns in . Specifically the average inner product similarity measure is defined as
(3) 
where is the average of the normalized inner products between pattern at pixel and the pattern in the dictionary .
Boundaries between grains are detected based on the lack of homogeneity of the matches to the dictionary of a pixel and its adjacent pixels, which we call a spatial patch. Specifically, for pixel we define the neighborhood similarity measure as the average amount of overlap between the NN neighborhood in of the pixel and the NN neighborhoods, with , in of the adjacent pixels on the patch:
(4) 
where are the indices of the neighbors of pixel . The neighborhood similarity will have value close to when the image patch is located in a grain. Its value will be close to zero when the image patch is centered on an anomaly. Its value will be between zero and one when the image patch is at a grain boundary (See Fig. 5).
The inner product similarities (for anomalies) and the neighborhood similarities (for grain boundaries) are used by a pattern classifier to assign each pixel to one of four classes. While many different types of unsupervised classifiers could be used, e.g., means, linear discriminant analysis (Hastie et al. (2005)
) or deep learning networks (
Hinton et al. (2006)), here we propose to cluster patterns using an unsupervised Decision Tree (DT) classifier (Hastie et al. (2005)) whose classification boundaries are determined so that they separate the modes (regions of concentration) of the histograms of innerproduct similarity and neighborhood similarity over the sample. The nonoverlapping modes can easily be separated by thresholding of the similarity value while the others can be estimated using a mode decomposition method such as Gaussian mixture modeling, also called mixture of Gaussian (MoG) modeling (
Figueiredo and Jain (2002)). This is illustrated in Fig. 1. The unsupervised DT classifier is illustrated in Fig. 2 in Sec. 7 for the IN100 sample considered. Four clusters were discovered by the model: anomalous pixels, which divided into two subclusters of shifted background and noisy background, and normal pixels, divided into grain boundary and grain interior.Decision tree classifiers have been previously applied to many imaging applications, e.g., land cover classification in remote sensing (Friedl and Brodley (1997)),(Pal and Mather (2003)). However, there are significant differences between the proposed DT classifier and those previously applied. First, the proposed classifier is a hybrid DT that uses special features (background compensated and uncompensated patterns) and similarity measures (inner products and neighborhood intersections) specific to materials microanalysis. Second, unlike standard nonparametric DT classifiers, the proposed DT is informed by a physics model through the generated dictionary of diffraction patterns. Third, our use of unsupervised Gaussian mixture models to determine the classification thresholds means that the DT classifier threshold parameters are determined by the Gaussian mixture models and do not need to be tuned, thus eliminating the need for labeled training data and time consuming crossvalidation.
3.2 Indexing
The proposed pixel indexing method is formulated as a statistical estimation problem. The pixel’s crystal orientation is estimated via maximum likelihood under a Von MisesFisher mixture density model for the EUler angles of the top dictionary matches. Note that the used in the maximum likelihood model may be different from the used to compute the neighborhood similarity for the DT classifier. We motivate the Von MisesFisher model as follows. Recall that the dictionary is generated for a set of predetermined orientations . Hence using simple table lookup, the indices of dictionary patterns found in the NN neighborhood of a pixel can be mapped to a set of orientations . If the pixel is in a grain then these orientation will be clustered around a true crystal orientation, that we call , at the pixel location. We extract a maximum likelihood estimate of this orientation using a statistical model for the variation of ’s.
We assume that the best matching orientations of a pixel form a random sample from an underlying marginal density , supported on the orientation sphere. Then the maximum likelihood estimator of is
(5) 
As is customary in the theory of directional statistics (Mardia and Jupp (1999)) we parameterize the density of dimensional orientations by their equivalent dimensional unit length quaternions (). Due to crystal symmetry, there are many () orientations that are equivalent to each other, i.e., the Euler angle representation is not unique. For any quaternion , the set of symmetry equivalent quaternions can be represented as where is a quaternion (rotation) matrix and
is the identity matrix. The matrices
are symmetric and constitute an abelian group of actions on the dimensional sphere. For example, Nickel has facecentered cubic symmetry, hence there are symmetryequivalent orientations. Since the quaternion representation of orientations associates two quaternions ( and ) to each orientation, we need matrices to establish a D representation of the point group (De Graef and McHenry (2007)).The proposed model for orientations is based on a generalization of the Von MisesFisher density (Mardia and Jupp (1999)) to group structured domains on the sphere. The standard Von MisesFisher density over the dimensional sphere with location parameter and precision parameter is defined as
(6) 
where and is the modified Bessel function of the first kind. Here and control the location of the mode (maximum) and spread of the density over the sphere, respectively. The natural generalization to the group structured domain is the periodic mixture density
(7) 
This Von MisesFIsher mixture (VMFm) density contains replicates of the Von MisesFisher density over the sphere centered at all the symmetryequivalent values of the location parameters .
Substitution of (7) into (5) and use of the wellknown invariance property of maximum likelihood estimation (Lehmann and Casella (1998)) gives a form for the maximum likelihood estimator of grain orientation in terms of the joint maximum likelihood estimators and :
(8) 
where . Even though this maximization problem appears daunting, it can be iteratively and efficiently computed by applying the well known expectationmaximization (EM) procedure for constrained parameter estimation in mixture models (McLachlan and Peel (2004)). For a full account of this procedure we refer the interested reader to Chen et al. (2015). Note that gives an empirical estimate of the degree of spread of the density about the orientation estimate Thus, is a measure of confidence of this estimate.
4 Generation of the Dictionary
The dictionary approach requires a uniform sampling of orientation space . Several sampling schemes are available in the literature; among the most popular schemes are a deterministic sampling method based on the Hopf fibration (Yershova and LaValle (2004),Yershova et al. (2009)) and the HEALPix framework (Hierarchical Equal Area isoLatitude Pixelization) (Gorski et al. (2005)). Neither of these approaches is easily adaptable for integration with crystallographic symmetry. Instead, we employ a recently developed strategy that starts from a simple D cubic grid which is mapped uniformly onto (Roşca et al. (2014)). This cubochoric mapping is uniform, refinable, and isolatitudinal, and consists of three steps:

a uniform cubic grid of grid points is generated inside a cube of edge length ;

the cube is divided into six pyramids with apex at the cube center and the cube faces as base, and each pyramid is mapped uniformly onto a sextant of a ball, using a generalization of the mapping of a square onto a curved square (Roşca and Plonka (2011));

all points inside the ball are then transformed, using a generalized inverse equalarea Lambert mapping, to the unit quaternion Northern hemisphere, which is isomorphous with .
From the quaternion representation one can readily derive other orientation parameterizations; the Rodrigues parameterization is most suitable for the determination of the orientations that belong to the fundamental zone (FZ) for a given crystal symmetry, because the boundaries of the FZ are planar. The more conventional Euler angle representation typically has curved surfaces as the boundaries of the FZ, so that Euler angles are less useful for uniform sampling approaches. It should be noted that lower crystal symmetry implies a larger dictionary, since the fundamental zone size increases with a reduction in symmetry; acceleration of the dot product calculations by means of a GPU (graphical processing unit) is a topic of ongoing research. The use of the Rodrigues representation to determine the dictionary elements has computational advantages, but care must be taken in the case of crystal symmetries with a single diad axis, for which the Rodrigues fundamental zone becomes infinite in the direction normal to the axis. Our approach still produces a uniform sampling of orientation space, although all rotations by an angle of ° are represented by points at infinity (which correspond to points on the outer cube surface in the cubochoric representation).
The dictionary used for the remainder of this paper was generated by setting
, and keeping only those orientations for which the corresponding Rodrigues vector lies inside the fundamental zone for the octahedral
point group. The patterns in the dictionary were down sampled to by pixel images. This results in a dictionary with elements of dimension . A representative (random) selection of dictionary elements in and is shown in Fig. 3. The left panel of Fig. 4 shows how the sampling points are distributed inside the octahedral Rodrigues FZ. The right panel of the figure shows the rate of dropoff of the top inner products in the compensated dictionary for randomly selected reference elements. This decay rate is used to select the number of nearest neighbors () used for the classifier described in Sec. 3 and implemented in Sec. 7.5 Computational Considerations
The online dictionary matching algorithm requires inner product evaluation between the measured sample diffraction patterns and the dictionary diffraction patterns. Let denote the number of patterns in the dictionary (dictionary size), denote the number of pixels on the photodetector (pattern size), and denote the number of pixels on the sample (sample size). The time complexity of calculating the mean inner product over the entire measured sample is .
For the indexing method, to determine the NN dictionary patterns for a given pixel the largest inner products need to be determined from the set of all inner products between the pixel and patterns in the dictionary. The time complexity of the whole process is , assuming . The computation time and space grow rapidly when the dictionary size and sample size become large.
Computational challenges can be addressed in several ways. The simplest approach is to use parallelization and distributed computation. All of the algorithms introduced here can be parallelized over the spatial domain since they involve local operations. For example, ML orientation estimation is applied independently to each pixel and DT classification is applied to each spatial patch in the sample. To speed up the NN dictionary search one can use methods from information retrieval such as dictionary caching and KD trees to accelerate the inner product evaluation and ranking process. These methods rely on the similarity of diffraction patterns over the NN neighborhood. However, as they rely on approximation, these methods may also introduce indexing errors. A study of these and other computational tradeoffs is important but is outside the scope of this paper.
6 Experimental Methods
To test the dictionary approach against an experimental data set, a polycrystalline IN100 Nickelbased superalloy sample was selected. The sample was polished using a multiplaten Robomet.3D, using a grit of micron diamond slurry on a TexMet cloth and finished with a nm colloidal silica slurry on a ChemoMet cloth. Between polishing steps, a water clean was used on a ChemoMet cloth.
A backscattered electron (BSE) image was recorded using a Tescan Vega XMH scanning electron microscope outfitted with a filament. An EBSD map was obtained with the same SEM and a Bruker eFlash1000 system. A tilt angle of °, voltage of kV, working distance of mm, and emission current of about nA were used to collect the EBSD map. The spatial resolution in both and directions was nm, and a Kikuchi pattern of pixels was acquired and stored at every point in a map.
7 Results
A dictionary was designed as described in Sec. 4 for the octahedral crystal symmetry group to match the known characteristics of the IN100 sample and the SEM system, as described in Sec. 6. Dictionary innerproducts and NN neighborhoods were computed from the detected pattern at each scan location.
We indicate four different scan locations (pixels) with qualitatively different patterns in Fig. 6. These locations are representative of the four different clusters of pattern, described below (see Fig. 2), that were discovered by the unsupervised DT classifier. In Fig. 7 the histograms of the dictionary innerproducts are shown for each of these pixels. The ”shifted background” and ”noisy background” pixels have innerproducts that are well separated from each other in addition to being separated from the innerproducts of the ”grain interior” and ”grain boundary” pixels. Thus one might expect that the group of ”anomalous” pixels, represented by the former two, could be easily separated from the group of ”normal” pixels, represented by the latter two, using any reasonable clustering technique based on the innerproducts. On the other hand, the innerproduct histograms for the grain interior pixel and the grain boundary pixel are overlapping. This overlap makes it difficult to distinguish these two types of pixels and justifies the need for the more sensitive neighborhood similarity measure that is better able to separate them.
Figure 8 shows the fullsample histograms of innerproducts with respect to and neighborhood similarities with respect to , respectively, over all pixels and over all patterns in the dictionaries. The left panel of Fig. 8 can be interpreted as the addition of all other innerproduct histogram to the four histograms shown in Fig. 7. Similarly to Fig. 7, the fullsample innerproduct histogram exhibits three well separated modes, which confirms that anomalous pixels can easily be separated from the normal pixels on the basis of thresholding each pixel’s average inner product measure. The three modes correspond to anomalous pixels with innerproducts clustered around and (not visible in the range of plotted) and normal pixels with innerproducts clustered around .
The right panel of Fig. 8 shows the histogram of all neighborhood similarities computed with neighborhood size . The latter histogram is bimodal and asymmetric about its mean. The higher mode located near corresponds to pixels whose NN neighborhood in has high overlap with the NN neighborhoods of its adjacent pixels. Such pixels are likely to be interior to a grain. The lower mode located near corresponds to patches of pixels near grain boundaries, patches that have less similar NN neighborhoods than ingrain pixels. To separate these tow modes, we fitted a two component mixtureofGaussian model to the histogram in the middle panel using the MoG EM algorithm (McLachlan and Peel (2004)). The result of this fit is shown in Fig. 1. The point of intersection of each of the fitted Gaussian densities (shown in the Figure by vertical dotted line) is used as the DT classification threshold for discriminating between grain boundaries and grain interiors.
Figure 2 shows the unsupervised DT classifier used to cluster observed patterns. The lower four nodes are leaf nodes while the upper three nodes are decision nodes for which thresholds are used to assign labels. These thresholds were determined from the observed histograms shown in Fig. 8 as described above. The DT classifies each pixel on the sample based on the pattern matches in the dictionaries. The top node, labeled ”observation” classifies pixel as a ”anomaly” or as ”normal” by thresholding the average innerproduct between the pattern at the center pixel of the patch and the patterns of the dictionary. Any threshold between and would separate the normal pair (grain interior, grain boundary) from the anomalous pair (noisy background, shifted background) and we selected the midpoint. The anomalous patterns are further subclassified on the left branch of the DT by applying a threshold between and to . The DT classifies normal pixels as either graininterior or grainboundary pixels by applying the threshold to the neighborhood similarities . Representative patterns are shown at the bottom of Fig. 2 that have been identified as belonging to each of the respective four clusters.
Figure 9 shows the pixel neighborhood similarities (left panel) and the pixel classifications (right panel) as images as determined by the unsupervised DT classifier. In the classification image the blue/red regions and black regions respectively correspond to pixels classified as anomalies and boundaries. These class labels can be used for segmentation of the grains and identification of the anomalies. A blowup of these images in Fig. 9 is shown in Fig. 10 for a small region.
Next we illustrate the use of the dictionary for estimation of the Euler angles in the sample. For the same subregion as in Fig. 10, Fig. 11 illustrates the OEM (original equipment manufacturer) orientation estimates (top left image), an estimate equal to the orientation of the element of the dictionary having largest normalized inner product with the pixel pattern (top right image), the proposed ML estimator of the orientation (VMFm location parameter) computed from the top matches in the dictionary (bottom left image), the proposed ML VMFm estimator computed with (bottom right image). The ML estimates of the VMFm scale parameters were computed for each pixel. This parameter is inversely proportional to the width of the estimated VMFm model over the dimensional quaternion sphere. Figure 12 shows images of these ML estimates translated into angular uncertainty (in degrees) by using the transformation , which is the width of the VMFm distribution in any fundamental zone. Note that the OEM image in Fig. 11 has many spurious orientation estimates within grains unlike the proposed dictionary based methods. Note also that the ML orientation estimates produce smoother ingrain orientations. The and
ML orientation estimates have low confidence (high variance) at some locations on grain boundaries and in anomalous region at bottom right. This low confidence is quantified by the ML estimator of the scale parameter
of the VMFm model, shown in Fig. 12 for and .8 Conclusion and Future Work
We have introduced a novel method for indexing polycrystalline materials that uses both mathematicalphysics modeling and mathematicalstatistics modeling. The physicsbased forward model is discretized into a dense dictionary of diffraction patterns that are indexed by Euler angle triplets. The dictionary is fixed for each crystal symmetry group and each SEM instrument configuration. The statistical model is based on the group symmetry of quaternion representation of the Euler angles on the sphere in dimensions. A feature of this method is that it performs classification, segmentation, and indexing in the unified framework of dictionary matching. A feature of the indexing method is that it incorporates a concentration parameter that can be estimated jointly with the Euler angles of a pixel or of a grain. This concentration parameter can be used to report the degree of confidence one can have in the Euler estimates. An iterative maximum likelihood estimator was proposed for estimating the orientation and associated confidence parameters in a statistical Von MisesFisher mixture model. The method was illustrated on a single sectional slice of a Nickel alloy sample.
As the proposed indexing method is pixel driven it is directly applicable to indexing over dimensional volumes. Future work will include algorithm acceleration to make full volumetric indexing fast enough to be practical. Potential acceleration methods include multiresolution and multiscale trees, fast coordinate ascent ML optimization, and parallelization. Other areas for future work include robustification of the dictionary to model mismatch, sensitivity to reductions in detector image resolution, and extensions of the dictionary approach to other electron diffraction modalities, such as electron channeling patterns and precession electron diffraction. Preliminary investigations indicate that, while dictionarybased classification appears to be robust to model mismatch, the proposed dictionarybased indexing algorithm is somewhat sensitive to model mismatch. This suggests that the dictionary design may need to be fine tuned to the SEM instrument in addition to the sample’s crystal symmetry group. The extension to other SEM modalities such as EDS is also possible but would require development of dictionaries that capture other types of data (e.g., spectra).
Acknowledgements
The authors thank Megna Shah of Bluequartz software for having provided the IN100 sample. We also thank Michael Groeber and Michael Uchic at AFRL for their comments on this work. The CMU portion of this work is supported by the Office of Naval Research on contract no. N000141210075. The UM portion of this work was partially supported by USAF/AFMC grant FA86509D5037/04 and by Air Force Office of Scientific Research grant FA95501310043.
References
 (1)
 Callahan and De Graef (2013) Callahan, P. G. and De Graef, M. (2013), ‘Dynamical electron backscatter diffraction patterns. Part I: Pattern simulations’, Microscopy and Microanalysis 19(05), 1255–1265.
 Chen et al. (2015) Chen, Y., Wei, D., Newstadt, G., De Graef, M., Simmons, J. and Hero, A. (2015), ‘Parameter Estimation in Spherical Symmetry Groups’, Signal Processing Letters, IEEE PP(99), 1–1.
 De Graef and McHenry (2007) De Graef, M. and McHenry, M. E. (2007), Structure of materials: an introduction to crystallography, diffraction and symmetry, Cambridge University Press.
 Diestel (2005) Diestel, R. (2005), ‘Graph theory (Graduate texts in mathematics)’.

Figueiredo and Jain (2002)
Figueiredo, M. A. and Jain, A. K. (2002), ‘Unsupervised learning of finite mixture models’,
Pattern Analysis and Machine Intelligence, IEEE Transactions on 24(3), 381–396.  Friedl and Brodley (1997) Friedl, M. A. and Brodley, C. E. (1997), ‘Decision tree classification of land cover from remotely sensed data’, Remote sensing of environment 61(3), 399–409.
 Gorski et al. (2005) Gorski, K. M., Hivon, E., Banday, A., Wandelt, B. D., Hansen, F. K., Reinecke, M. and Bartelmann, M. (2005), ‘HEALPix: a framework for highresolution discretization and fast analysis of data distributed on the sphere’, The Astrophysical Journal 622(2), 759.
 Hastie et al. (2005) Hastie, T., Tibshirani, R., Friedman, J. and Franklin, J. (2005), ‘The elements of statistical learning: data mining, inference and prediction’, The Mathematical Intelligencer 27(2), 83–85.
 Hinton et al. (2006) Hinton, G., Osindero, S. and Teh, Y.W. (2006), ‘A fast learning algorithm for deep belief nets’, Neural computation 18(7), 1527–1554.
 Lehmann and Casella (1998) Lehmann, E. L. and Casella, G. (1998), Theory of point estimation, Vol. 31, Springer Science & Business Media.

Mardia and Jupp (1999)
Mardia, K. V. and Jupp, P. E. (1999), ‘Directional statistics’.
http://cds.cern.ch/record/1254260 
McLachlan and Peel (2004)
McLachlan, G. and Peel, D. (2004),
‘Finite Mixture Models’.
http://cds.cern.ch/record/997567  Pal and Mather (2003) Pal, M. and Mather, P. M. (2003), ‘An assessment of the effectiveness of decision tree methods for land cover classification’, Remote sensing of environment 86(4), 554–565.
 Rauch and Dupuy (2005) Rauch, E. and Dupuy, L. (2005), ‘Rapid spot diffraction patterns idendification through template matching’, Archives of Metallurgy and Materials 50, 87–99.
 Rauch et al. (2008) Rauch, E., Veron, M., Portillo, J., Bultreys, D., Maniette, Y. and Nicolopoulos, S. (2008), ‘Automatic crystal orientation and phase mapping in TEM by precession diffraction’, Microscopy and AnalysisUK (128), S5.
 Roşca et al. (2014) Roşca, D., Morawiec, A. and De Graef, M. (2014), ‘A new method of constructing a grid in the space of 3d rotations and its applications to texture analysis’, Modelling and Simulation in Materials Science and Engineering 22(7), 075013.
 Roşca and Plonka (2011) Roşca, D. and Plonka, G. (2011), ‘Uniform spherical grids via equal area projection from the cube to the sphere’, Journal of Computational and Applied Mathematics 236(6), 1033–1041.
 Schwartz et al. (2009) Schwartz, A. J., Kumar, M., Adams, B. L. and Field, D. P. (2009), Electron backscatter diffraction in materials science, Vol. 2, Springer.
 Tao and Eades (2005) Tao, X. and Eades, A. (2005), ‘Errors, artifacts, and improvements in EBSD processing and mapping’, Microscopy and Microanalysis 11(01), 79–87.
 Tropp and Wright (2010) Tropp, J. A. and Wright, S. J. (2010), ‘Computational methods for sparse solution of linear inverse problems’, Proceedings of the IEEE 98(6), 948–958.
 Wright and Nowell (2006) Wright, S. I. and Nowell, M. M. (2006), ‘EBSD image quality mapping’, Microscopy and Microanalysis 12(01), 72–84.
 Wright et al. (1991) Wright, S. I., Zhao, J.W. and Adams, B. L. (1991), ‘Automated determination of lattice orientation from electron backscattered Kikuchi diffraction patterns’, Texture, Stress, and Microstructure 13(23), 123–131.
 Yershova et al. (2009) Yershova, A., Jain, S., Lavalle, S. M. and Mitchell, J. C. (2009), ‘Generating uniform incremental grids on SO (3) using the Hopf fibration’, The International journal of robotics research .
 Yershova and LaValle (2004) Yershova, A. and LaValle, S. M. (2004), Deterministic sampling methods for spheres and SO (3), in ‘Robotics and Automation, 2004. Proceedings. ICRA’04. 2004 IEEE International Conference on’, Vol. 4, IEEE, pp. 3974–3980.