Dictionary-Based Indexing Approach

Dictionary-based indexing refers to a wide class of algorithms that employ a library of precomputed patterns to index experimental diffraction patterns. Instead of extracting features, such as bands, spots etc., from the experimental patterns, the dictionary approach uses all the available pixels in the pattern to find a best match with respect to the dictionary patterns. There are three necessary ingredients to generate this library: (1) a generalized forward projector (GFP), which is a physics-based forward model of the scattering process; (2) an accurate detector and noise model; and (3) a method to uniformly sample orientation space, SO(3). Although the GFP and the detector and noise model will change for different diffraction modalities, the method to uniformly sample orientation space is common to all modalities; this has been discussed in detail in (Singh & De Graef, Reference Singh and De Graef2016). The three ingredients combined with a pattern matching engine (Goshtasby, Reference Goshtasby2012) (cross-correlation, mutual information, etc.) form the basis of the dictionary-based approach to indexing diffraction patterns. In the remainder of this section, we will discuss in detail each aspect of the dictionary-based approach as it is adapted to the ECP modality. A schematic for dictionary-based indexing of ECPs is shown in Figure 1.

Figure 1 Schematic representation of the electron channeling patterns (ECP) dictionary indexing process. Arrows represent the flow of data/information. GFP, generalized forward projector.

Forward Model for ECPs

The GFP described in this section has three main components: dynamical scattering of electrons, Monte Carlo (MC) simulations to generate spatial and depth statistics, and a realistic detector, noise and distortion model. In this section, we will briefly consider each of these components to generate a computational framework for the simulation of ECPs. The dynamical scattering accounts for the elastic channeling of the electrons by the crystal potential, whereas the MC simulation accounts for the stochastic nature of inelastic scattering. The MC simulations also provide the backscattering depth statistics, which are needed to set the integration depth for the dynamical simulations. The two approaches are merged together to compute a “master channeling pattern” that represents the diffracted intensity distribution on a spherical surface surrounding the crystal. The pattern is stored as a square intensity grid, using a modified Lambert equal-area projection (Roşca, Reference Roşca and Plonka2010). The master pattern can be sampled for a given sample and detector geometry, as well as a given crystal orientation, which can be specified as a triplet of Euler angles, a unit quaternion, etc. A similar computational approach has already been used for the simulation of EBSD patterns (Callahan & De Graef, Reference Callahan and De Graef2013). Such a GFP is also necessary to make the dictionary indexing method computationally tractable as a large number of dynamical diffraction patterns need to be calculated for reliable indexing. The master pattern for Silicon using the modified Lambert projection along with the stereographic projection are shown in Figures 2a and 2b, respectively. In the next subsection, we will briefly outline the dynamical scattering process and the Bloch wave method used to compute the BSE amplitude for different channeling directions, as well as the detector model. Finally, we conclude the section by discussing the relevant noise model and electron optical aberrations.

Figure 2 Master electron channeling pattern for 20 keV beam acceleration voltage for Si using (a) modified Lambert projection and (b) stereographic projection.

Dynamical Diffraction

A number of theoretical models have been proposed to explain channeling contrast as a function of incoming beam direction with respect to the crystal reference frame (Spencer & Humphreys, Reference Spencer and Humphreys1980; Marthinsen & Høier Reference Marthinsen and Høier1986; Rossouw et al., Reference Rossouw, Miller, Josefsson and Allen1994; Dudarev et al., Reference Dudarev, Rez and Whelan1995; Winkelmann et al., Reference Winkelmann, Schröter and Richter2003). For the present study, we focus on the theoretical model outlined in Picard et al. (Reference Picard, Liu, Lammatao, Kamaladasa and De Graef2014), which distinguishes between BSE1 and BSE2 electrons; BSE1 electrons undergo a backscatter event as their very first scattering event and are responsible for the channeling contrast. These electrons have nearly the same energy as the incident electrons as they immediately backscatter. The BSE2s, on the other hand, backscatter after undergoing multiple elastic and inelastic scattering events. These electrons carry no channeling information and form the background intensity. It has been shown in (Picard et al., Reference Picard, Liu, Lammatao, Kamaladasa and De Graef2014) that for a reasonable detector geometry, the typical ratio of BSE1s to BSE2s is of the order of 10⁻³–10⁻⁴. Therefore, to detect the weak BSE1 signal superimposed on top of the large BSE2 background signal, the brightness and contrast settings of the microscope have to be adjusted accordingly. Most of the BSE2 signal is usually removed by adjusting pattern brightness and contrast such that the dominant signal results from BSE1 electrons. Therefore, in the remainder of this paper, we will only focus on the BSE1 electrons.

The BSE1s can be generated from a range of depths inside the sample. Consequently, the overall signal will be a sum of the signals coming from different depths. Note that the incident electron channels on its way into the crystal, which is different from the EBSD technique, where the electron channels on its way out of the crystal. The range of depths from which the BSE1 signal arises can be estimated by means of MC simulations. There will also be normal and anomalous absorption of the BSE1 electrons as they channel through the crystal. This absorption is described phenomenologically by the imaginary part of the lattice potential and a depth dependent weight factor, whereas integrating the signal from various depths. For a crystal structure with N _a atoms in the unit cell distributed over n sites in the asymmetric unit, with a set S _n of equivalent atom positions, we can write the overall BSE1 signal for a channeling direction described by the wave vector k ₀ as

(1) $${\cal P}({\bf{k}} _{0} ){\equals}\mathop{\sum}\limits_n \mathop{\sum}\limits_{i\in S_{n} } {{Z_{n}^{2} D_{n} } \over {z_{0} }}{\int}_0^{z_{0} } \lambda (z)\mid\Psi _{{{\bf{k}} _{0} }} ({\bf{r}} _{i} )\!\mid^{2} dz,$$

where Z _n and D _n are the atomic number and Debye–Waller factor, respectively, of the atom at site n; z ₀ is the depth of integration; λ(z) is a depth dependent weight factor, and Ψ _k0(r _i ) is the wave function of the electron which is channeling in the direction k ₀. The electron wave function is evaluated using either the Bloch wave or the scattering matrix approach (De Graef, Reference De Graef2003); in this paper, we restrict our discussion to the Bloch wave approach, in which the electron wave function is written as a superposition of Bloch waves traveling in different directions k ^(j)

(2) $$\Psi ({\bf r}){\equals}\mathop{\sum}\limits_j \alpha ^{{(j)}} \mathop{\sum}\limits_g C_{{\bf g}}^{{(j)}} e^{{2\pi i({\bf k}^{{(j)}} {\plus}{\bf g}) \cdot {\bf r}}} .$$

Here, α ^(j) is the excitation amplitude of the jth Bloch wave and $C_{{\bf g}}^{{(j)}} $ is the Bloch wave coefficient of the g beam for that wave. The above two equations can be combined to obtain the overall BSE1 signal for an electron incident in the k ₀ direction. It can be shown that the result can be expressed as the sum of elements in a Hadamard product (denoted by the symbol °) of two matrices, S and L. The matrix S can be thought of as a structure factor, whereas the matrix L encodes the dynamical diffraction information. The individual elements of S and L are given by

(3) $$\eqalignno{ & S_{{{\bf gh}}} {\equals}\mathop{\sum}\limits_n {\mathop{\sum}\limits_{{\rm i}\inS_{n} } {Z_{n}^{2} D_{n} {\rm e}^{{2\pi {\rm i}(h{\minus}g) \cdot {\rm r}_{i} }} ^; } } \cr & L_{{{\bf gh}}} {\equals}\mathop{\sum}\limits_j {\mathop{\sum}\limits_k {C_{{\rm g}}^{{(j){\asterisk}}} \alpha ^{{(j){\asterisk}}} {\scr I}_{{jk}} } } \alpha ^{{(k)}} C_{h}^{{(k)}} . $$

Here, ${\scr I}_{{jk}} $ is defined as

(4) $${\scr I}_{{jk}} {\equals}{1 \over {z_{0} }}{\int}_0^{z_{0} } {\lambda (z){\rm e}^{{{\minus}2\pi (\alpha _{{jk}} {\plus}{\rm i}\beta _{{jk}} )z}} } dz,$$

where α _jk =q ^j +q ^k and β _jk =γ ^j −γ ^k ; the complex numbers λ ^j ≡γ ^j +iq ^j are the eigenvalues of the dynamical Bloch matrix with absorption taken into account. Finally, all of this can be put together to obtain the depth-integrated probability of backscattering of the electron traveling in the direction k ₀, which is given by

(5) $${\cal P}({\bf k}_{0} ){\equals}\langle {\bf u}^{T} \!\mid\!S \circ L \!\mid\!{\bf u}\rangle .$$

Here, u is a column vector with all its entries equal to unity and ° represents the Hadamard (element-wise) product of two matrices. The above equation is nothing but the sum of elements of the Hadamard product of S and L.

Detector Model

A channeling pattern is recorded by measuring the backscattered signal as a function of the incident beam direction. MC simulations for BSE1 electrons reveal that the majority of BSE1s are scattered at a large angle with respect to the incident beam direction, typically in the range 50–70°. Figure 3a shows the stereographic projection of the BSE1 exit direction distribution for 20 keV electrons scattered from a Si sample. A polar plot of the BSE1 intensity as a function of angle with respect to the sample normal is shown in Figure 3b. The maximum intensity is at an exit angle of 58.1°. As it is desirable to capture a large portion of the backscattered electrons to enhance the already weak BSE1 cumulative signal, an annular integrating backscatter detector is used, as shown in the schematic of Figure 4a. The inner and outer radii of the detector and the sample to detector distance, ξ are adjusted such that the detector captures a significant fraction of the total backscattered electron flux. For defect imaging, on the other hand, the sample is usually tilted to be close to a Bragg orientation. As a result, the background intensity will no longer have rotational symmetry around the incident beam direction. This asymmetry can be taken into account when the dynamical backscatter intensities are merged with the BSE1 results from stochastic MC simulations.

Figure 3 MC simulation for 20 keV beam acceleration voltage for Si (a) stereographic projection of BSE1 electrons and (b) polar plot of backscatter yield as a function of angle from sample normal with the maximum at 58.1°.

Figure 4 a: Schematic of the electron channeling patterns setup. Only a small fraction (black dashed) of the BSE1 electrons get captured by the annular detector and (b) stereographic projection of BSE1 exit direction for 20 keV acceleration voltage on Si sample superimposed with the stereographic projection of all exit directions for which the electron is captured by the detector (detector geometry in text).

The MC simulations are performed for a range of incident beam angles, typically from 0–20° in steps of 1°, resulting in the BSE1 backscatter yield as a function of exit direction and incidence angle, ${\bf Y}(\theta ,{\bf \hat{n}})$ . The inner and outer radius of the annular detector, together with the working distance, are used to calculate the solid angle range over which the MC simulations need to be integrated. The detector is discretized with a polar and azimuthal angular step size of 1°, producing a set of exit directions {e _i }. For every incidence angle, θ _j from 0–20°, the effective weight factor is then calculated by summing the backscatter yield over the calculated solid angle range; mathematically, we have

(6) $$w(\theta; \,R_{{{\rm in}}} ,\,R_{{{\rm out}}} ,\,\xi ){\equals}\mathop{\sum}\limits_{\{ {\bf e}_{i} \} } {{\bf Y}(\theta ,\,{\vskip -3pt \bf \hat}{\hskip -3pt \bf n}})} .$$

For a realistic detector geometry with inner radius, R _in=5.0 mm, outer radius, R _out=13.0 mm and a working distance of WD=5.0 mm, the fraction of electrons captured by the detector is shown in Figure 4b. The weight factor for an arbitrary incidence angle is calculated by linear interpolation using weight factors of the known angles. For a grain with orientation described by a unit quaternion, $$\tilde{q}$$ , each pixel (i, j) in the channeling pattern corresponds to an incident wave vector k _ij in the microscope reference frame; the intensity in the corresponding direction is interpolated from the MC and master patterns and is given by the expression

(7) $$I_{{{\rm BSE}}} {\equals}w(\theta;\,R_{{{\rm in}}} ,R_{{{\rm out}}} ,\xi ){\cal M}({\bf k}_{{ij}};\tilde{q}),$$

where θ is the angle between the wave vector k _ij and the sample normal and ${\cal M}({\bf k}_{{ij}} )$ is extracted from the master pattern using bilinear interpolation.

Lens Aberrations and Noise

Aberrations in the channeling patterns arise from the fact that the incident electron beam is scanned at a relatively large angle (10–15°) with respect to the optical axis. The primary (Seidel) aberrations of the objective lens (spherical aberration, coma, astigmatism, field curvature, and distortion) combine with the aberrations of the deflection system to introduce a significant amount of distortion in the high-angle portion of the ECPs. We model the distortion in the following way: if P _i represents a point in the image plane corresponding to a point P ₀ in the object plane, then we assume that the location of P _i is given by u _i +Δu _i , where ∆u _i =(D+id)|u _o|² u _o; u ₀ is a complex number representing the location of point P ₀, and D and d are the real and imaginary parts of the distortion coefficient. Depending on the sign of D, the distortion is either referred to as a barrel distortion (D>0) or a pin-cushion distortion (D<0). Note that the displacement field only depends on the position in the object plane, which makes the computation of the distortion field relatively straightforward. A full determination of the aberration field requires consideration of higher order terms. More details about this and various other SEM lens aberrations can be found in (Szilagyi, Reference Szilagyi1988; Hawkes & Kasper, 1989 Reference Hawkes and Kasperb , 1989 Reference Hawkes and Kaspera ).

Figure 5a shows an experimental channeling pattern for a semiconductor grade poly-Si sample with a large grain size (500 μm), acquired on the TESCAN MIRA 3 FE-SEM equipped with a rocking beam setup; the acquisition parameters were: acceleration voltage V=30 kV, opening angle of the incident beam cone, θ _c =14.35°, working distance, 6.2 mm. A pattern fit resulted in the following Euler angles: (φ ₁, Φ, φ ₂)≡(121.1°, 46.2°, 203.7°). Figure 5b shows the corresponding simulated ECP without distortion. Figures 5c and 5d show the addition of barrel distortion and Poisson noise, respectively, to the simulated pattern. The distortion coefficient, D=(0.4, 1.0)×10⁻⁶ was found to give the best match between experimental and simulated patterns; it is clear from a visual comparison that barrel distortion alone does not fully reproduce the distortions present in the experimental pattern. An adaptive histogram equalization (Pizer et al., Reference Pizer, Amburn, Austin, Cromartie, Geselowitz, Greer, Romney, Zimmerman and Zuiderveld1987) has been performed on both the experimental and simulated patterns to enhance and match the contrast. Note that there is significant bending of the Kikuchi bands close to the border of the pattern. We shall see later that even though the noise and distortion are important for improving the overall agreement between simulated and experimental patterns, the dictionary approach is relatively insensitive to these second-order corrections.

Figure 5 a: Experimental channeling pattern. b: Simulated channeling pattern for correct detector parameters and Euler angle triplet. c: Distortion added to the simulated pattern (d) both distortion and Poisson noise added to simulated pattern. The contrast for both experimental and simulated patterns have been enhanced using the adaptive histogram equalization algorithm.

Estimating Detector Parameters

For the dictionary approach to work correctly, the right set of geometric detector parameters for a particular experiment must be estimated. For ECPs, the relevant parameters are the working distance, the sample tilt and the incident beam cone semi-angle. The problem of finding the correct set of detector parameters can be reformulated in terms of an optimization problem, where the best set of parameters will maximize either the mutual information or the dot product between the simulated and experimental images (Goshtasby, Reference Goshtasby2012). As the dependence of these metrics on the detector parameters is not known analytically, DFO algorithms are well suited to this optimization problem. The performance of a number of DFO algorithms, in terms of finding the global optimum and refining near optimal solutions, has been documented in (Rios & Sahinidis, Reference Rios and Sahinidis2013). For the present paper, two of these algorithms, namely the Nelder–Mead simplex (Nelder & Mead, Reference Nelder and Mead1965) and the bound optimization by quadratic approximation (BOBYQA) (Powell, Reference Powell2009) were evaluated. It is also worthwhile to note that further refinement of the Euler angles obtained from the dictionary approach can also be performed using this method, with the Euler angles obtained from the dictionary method serving as a starting point for the refinement algorithm. The details of obtaining the detector parameters using dynamically simulated diffraction patterns will be published elsewhere.

Dictionary Generation and Indexing

The dictionary approach requires a uniform sampling of orientation space. This is essential so that different regions of SO(3) are represented with equal weight in the pattern dictionary. It is also important to note that for a given crystal structure, only the relevant FZ needs to be sampled as all possible unique orientations for the given crystal symmetry are represented by the FZ. For the present study, we will use the Rodrigues FZ because, unlike in other orientation representations, the FZ in Rodrigues space has planar boundaries and is convenient for sampling. Although there are other sampling methods in the literature, such as the Hopf method described in Yershova et al. (Reference Yershova, Jain, LaValle and Mitchell2010) and the HEALPix framework described in Górski et al. (Reference Górski, Hivon, Banday, Wandelt, Hansen, Reinecke and Bartelmann2005), it is not straightforward to adapt them for integration with crystallographic symmetry. In this study, we have employed the cubochoric sampling method (Roşca et al., Reference Roşca, Morawiec and De Graef2014), which is an equal-volume mapping of a uniform three-dimensional (3D) cubic grid onto a grid on the unit quaternion sphere. Thus, the sampling is uniform in the equal-volume sense. This cubochoric mapping is uniform, hierarchical, and isolatitudinal, and is carried out in three steps:

• ∙ Generate a uniform grid inside a cube of edge length π^2/3. With N being the number of sampling points along a semi-edge of the cube, there are two choices for the grid. The first choice, S ₀₀₀ is a grid of (2N+1)³ points containing the identity orientation. The other choice is the dual grid, $S_{{{1 \over 2}{1 \over 2}{1 \over 2}}} $ which has 8N ³ points and does not contain the identity orientation. The fineness of the resulting orientation sampling depends on the value of N (Roşca et al., Reference Roşca, Morawiec and De Graef2014).

• ∙ Mapping the uniform grid on the cube to a uniform grid on an homochoric ball of radius (3π/4)^1/3. This is done by dividing the cube into six square pyramids with apex at the center of the cube and mapping each pyramid to a sextant of the ball (Roşca & Plonka, Reference Roşca2011).

• ∙ Finally, the uniform grid on the ball is mapped onto the Northern hemisphere of the unit quaternion sphere, isomorphic with the space of 3D orientations, using the generalization of the azimuthal equal-area Lambert projection, also referred to as the homochoric projection (Roşca et al., Reference Roşca, Morawiec and De Graef2014).

The dictionary can be generated by combining all the steps in the previous section and serves as a precomputed look-up table against which the experimental patterns are matched. As the dictionary is generated on a discrete grid of orientations, it is not necessary and highly unlikely that the exact orientation corresponding to the experimental pattern is present in the dictionary. For the dictionary approach to be effective, the pattern matching step needs to pick out the pattern that is closest in misorientation to the experimental pattern. Therefore, the criterion used for pattern matching must serve as a proxy for the misorientation between the dictionary and experimental patterns. Using uniform sampling of misorientation iso-surfaces, it has been shown in Singh & De Graef (Reference Singh and De Graef2016) that the dot product method is a good proxy for misorientations up to about 3–5° for electron backscatter diffraction patterns (EBSPs). As ECPs and EBSPs are very similar, the dot product was used as the similarity metric in this work as well. Each experimental channeling pattern is arranged in a column vector, e _i (i∈{1, 2, … , N _e }) and normalized to give unit vectors $${\bf {\vskip -3pt \hat} {\hskip -2pt e}} {\equals}{\bf e}_{i} \,/\mid\mid\,{\bf e}_{i} \mid\mid\,$$ , where ||·|| represents the Euclidean norm of the vector. These patterns are then matched with each normalized dictionary member, ${{\bf \vskip -4pt \hat} {\hskip -4pt \bf d}}_{j} $ (j∈{1, 2, … , N _d }) using the dot product metric given by

(8) $$\alpha _{{ij}} {\equals}\,{\bf {\vskip -3pt \hat} {\hskip -2pt e}} \cdot {\bf {{\vskip -3pt \hat} {\hskip -2pt d}}_{j} .$$

For each experimental pattern ${\bf \vskip -2.7pt \hat} {\hskip -2.3pt \bf e}}_{i} $ , the dictionary element with the highest dot product, σ _ij is considered to be the best match and the orientation corresponding to that dictionary element is assigned as the orientation for the experimental pattern. To achieve an average grid spacing of about 1.4°, the value of N=100 must be used in the S ₀₀₀ grid. Even for this modest angular step size, the number of entries in the dictionary is 333,227 for cubic (octahedral) rotational symmetry. The number of patterns can very easily become greater than a few million for crystals with lower symmetry, with more than eight million entries in the absence of any symmetry (for N=100); this makes the dictionary approach computationally intensive. The approach is made computationally tractable by using a combination of OpenMP parallellization for the dictionary generation with a massively parallel graphics processing unit platform using the OpenCL approach for the computation of the dot products between the dictionary and experimental patterns.

Summary

In this paper, we have introduced a new dictionary-based approach to indexing ECPs. The basic ingredients of this method include a GFP, which is a physics-based forward model for the scattering process, a detector and noise model and a method to uniformly sample orientation space. The deterministic forward model is merged with stochastic MC simulations to calculate the scattering master pattern. Individual channeling patterns are extracted from the master pattern using bilinear interpolation. The cubochoric representation is introduced as a means to uniformly sample orientation space. The forward model coupled with a DFO algorithm are used to evaluate the correct detector parameters and to refine the orientations obtained from the dictionary approach. In contrast to the Hough transform-based method which extracts features from diffraction patterns, our approach uses all available pixels. Pattern matching is performed using the normalized dot product as the similarity metric. The method is applied to poly-Si as the model system. The results are in excellent agreement with the calibrated values. Finally, an error analysis is performed using a set of 1,000 simulated patterns with random orientations. The results show that the method produces a mean orientation error of 1° with 0.5° SD for the true detector parameters and N=100 points along the cubchoric cube semi-edge, corresponding to a mean sampling step size of 1.4°. Increasing the number of sampling points decreases the error. The mean orientation error is insensitive to an error in the detector parameters up to about 5% error in detector parameter.