Accurate reconstruction of EBSD datasets by a multimodal data approach using an evolutionary algorithm

TL;DR

This paper introduces a novel multimodal data approach using an evolutionary algorithm to accurately correct distortions in EBSD datasets, improving phase differentiation without human annotation.

Contribution

The method combines segmented images and an evolutionary algorithm to correct EBSD distortions, enabling precise, automated phase identification in large, complex datasets.

Findings

Achieves distortion correction with accuracy on the order of pixel size.

Applicable to multiphase materials with feature sizes down to 1 μm.

Does not require human annotation or strong assumptions on distortion shape.

Abstract

A new method has been developed for the correction of the distortions and/or enhanced phase differentiation in Electron Backscatter Diffraction (EBSD) data. Using a multi-modal data approach, the method uses segmented images of the phase of interest (laths, precipitates, voids, inclusions) on images gathered by backscattered or secondary electrons of the same area as the EBSD map. The proposed approach then search for the best transformation to correct their relative distortions and recombines the data in a new EBSD file. Speckles of the features of interest are first segmented in both the EBSD and image data modes. The speckle extracted from the EBSD data is then meshed, and the Covariance Matrix Adaptation Evolution Strategy (CMA-ES) is implemented to distort the mesh until the speckles superimpose. The quality of the matching is quantified via a score that is linked to the number of overlapping pixels in the speckles. The locations of the points of the distorted mesh are compared to those of the initial positions to create pairs of matching points that are used to calculate the polynomial function that describes the distortion the best. This function is then applied to un-distort the EBSD data, and the phase information is inferred using the data of the segmented speckle. Fast and versatile, this method does not require any human annotation and can be applied to large datasets and wide areas. Besides, this method requires very few assumptions concerning the shape of the distortion function. It can be used for the single compensation of the distortions or combined with the phase differentiation. The accuracy of this method is of the order of the pixel size. Some application examples in multiphase materials with feature sizes down to 1 $\mu$m are presented, including Ti-6Al-4V Titanium alloy, Rene 65 and additive manufactured Inconel 718 Nickel-base superalloys.

Figures

Tables

Table 1. Table 1: Precision of the location of the phase boundaries versus grain boundaries in the profiles of Fig. 8 , in the Rene 65 dataset

Boundary	$\Delta x$ ($\mu m)$
a	0.0
b	0.0
c	0.0
d	0.0
f	0.1
g	0.2
j	0.0

Equations

$$\text{Similarity}({\mathcal{I}}^{EBSD},{\mathcal{I}}^{BSE})=2\frac{|{\mathcal{I}}_{seg}^{EBSD}\cap {\mathcal{I}}_{seg}^{BSE}|}{|{\mathcal{I}}_{seg}^{EBSD}|+|{\mathcal{I}}_{seg}^{BSE}|}$$

$$f^{\ast }$$

$$\displaystyle=\mathop{\mathrm{argmax}}_{f}\Big{(}2\frac{|f(\mathcal{I}^{EBSD})_{seg}\cap\mathcal{I}_{seg}^{BSE}|}{|f(\mathcal{I}^{EBSD})_{seg}|+|\mathcal{I}_{seg}^{BSE}|}\Big{)}$$

$$f(x,y)=\{\begin{array}{l}{x}^{\prime }={\sum }_{n=0}^P{\sum }_k^n{c}_{k,n-k}^x{x}^k{y}^{n-k}\\ y^{\prime }={\sum }_{n=0}^P{\sum }_{k=0}^n{c}_{k,n-k}^y{x}^k{y}^{n-k}\end{array}$$

$$f(x,y)=\{\begin{array}{l}{x}^{\prime }=c_{0,0}^x+c_{1,0}^{x}x+c_{0,1}^{x}y+c_{1,1}^{x}xy+c_{2,0}^x{x}^2+c_{0,2}^x{y}^2\\ y^{\prime }=c_{0,0}^y+c_{1,0}^{y}x+c_{0,1}^{y}y+c_{1,1}^{y}xy+c_{2,0}^y{x}^2+c_{0,2}^y{y}^2\end{array}$$

$$x_k^{(g+1)}\leftarrow m^{(g)}+{\sigma }^{(g)}\mathcal{N}(0,{\mathcal{C}}^{(g)}),fork=1,2,...,\lambda .$$

$$\Delta x=\Vert x_{Ph.B}-x_{GB}\Vert$$

Full text

††journal: Materials Characterization

1 Introduction

1.1 Limits of phase differentiation by EBSD

Electron Backscattered Diffraction (EBSD) is a powerful tool for gathering local crystallographic data in multiphase materials. This technique enables a fast mapping of microstructures with a good angular resolution (typically 0.5$\degree$), based on the collection of diffraction patterns and the analysis of Kikuchi bands on those patterns. Initially implemented in scanning electron microscopes (SEM), it has been adapted to transmission microscopy (TEM) under the name of Transmission Kikuchi Diffraction (TKD) [1, 2]. In the past years, tridimensional EBSD (3D EBSD) by serial sectioning has also emerged as a new tool to reconstruct and analyze materials in three dimensions [3, 4, 5, 6, 7, 8]. To index the collected patterns, most commercial software use algorithms based on Hough transformations to detect the characteristic bands on the Kikuchi patterns. This process has the advantage of speed but suffers from some limitations. Because of the diffraction patterns they produce, some phases are intrinsically challenging to index with this method. This is the case for the finely distributed $\beta$ phase in $\alpha-\beta$ Titanium alloys, for example. Another well-known limitation of this method happens in the presence of two phases of similar structures. This is the case of $\gamma$ and $\gamma^{\prime}$ phases in Nickel and Cobalt-base superalloys. The ordered $L1_{2}$ $\gamma^{\prime}$ phase produces extra very low-contrast bands compared to the Face Centered Cubic (FCC) $\gamma$, that are undetectable by Hough transforms. In the case of phases producing similar Kikuchi patterns but having significantly different chemistries, one can employ multi-modal data acquisition. Nowell and Wright [9] have proposed the simultaneous collection of X‐ray energy dispersive spectroscopy (XEDS) data, along with EBSD data. This method has been shown to enable phase differentiation in complex alloys, and reveal unique crystallographic relationships between phases [10, 11, 12]. However, XEDS-EBSD requires a significant difference in the chemistry of the phases of interest and a minimum number of counts on the XEDS spectra, which can lower the speed of data collection. The acceleration voltage also has to be high enough to enable the detection of the characteristic X-rays emission of the elements of interest. Increasing the acceleration voltage would also result in more counts on the EDS spectra and increase the indexing speed, but is detrimental to the spatial resolution. Indeed, the XEDS interaction volume is about one micrometer at 20 kV in Nickel alloys and increases greatly with the acceleration voltage.

Using the different chemical composition of the phases, other methods involving multi-modal data combination have been implemented. Phases of different composition exhibit different intensities on back-scattered electrons (BSE) images, due to their different Z number. Payton and Nolze [13] have proposed a method that involves the acquisition of BSE images on a tilted sample, using specific BSE detectors located above the phosphorous screen of the EBSD detector. Those authors have shown that this method enables a more accurate phase differentiation than XEDS-EBSD, since the contrast on the back-scattered images is very sensitive to the variations in chemical composition. They also demonstrate the capability of that method for segmenting small features precisely, with typical dimensions of the order of a few micrometers. The limit of this technique is that not all the EBSD detectors are equipped with such built-in diode sensors. The collection of back-scattered data is usually made by relatively large retractable BSE detectors that have several quadrants, at a 0$\degree$ tilt angle. On the other hand, EBSD data is collected at a 70$\degree$ tilt angle. This, among other reasons, leads to complex distortions in the EBSD versus the BSE data, which makes their precise recombination very challenging.

1.2 Distortions in EBSD data

EBSD data is known to suffer from distortions that arise from many instrument and detector artifacts [14]. The following types of distortions can be distinguished:

First order spatial distortions: those are induced when tilting the sample at 70$\degree$. An area that would be a square at 0$\degree$ tilt turns into a trapezoid if the rotation axis of the stage is not perfectly parallel to the surface of the sample. The borders of the trapeze are also misaligned with those of the non-tilted area. Those distortions can be modeled by a first order affine transformation, involving a translation and a rotation.

• 2.

Second order spatial distortions: the electron beam is deflected by a set of lenses that induce a barrel distortion on the final image. This distortion is mostly visible at low magnification. It can be modeled by a second order polynomial function.

• 3.

Space-time distortions: when a sample is being exposed for some time to the electron beam, charging effects occur. The accumulation of charges on the surface of the sample leads to the deflection of the beam during scanning. This effect occurs in any scanning process. It is usually the most pronounced in the beginning of the scan and becomes more stable as the scan progresses. A good example is shown in [15]. The amount of drift is related to the conductivity of the sample itself, the conductivity of the mounting system that is being used, the voltage and aperture, as well as the cleanliness of the scanned surface. The amount of drift is also more important at high than low magnifications [16]. There is no simple physical model to describe the pattern of the drift distortions.

Given the complexity of the net resulting distortion, it cannot be easily calculated.

2 Method

In this context, we propose a method that enables the compensation of the distortions in EBSD data and adds phase differentiation if wanted, by accurately matching it with BSE images collected at a 0$\degree$ tilt angle, at individual pixel precision and over broad areas. To do so, the Covariance Matrix Adaptation Evolution Strategy (CMA-ES) is used to calculate the relative distortion function between the EBSD data and the BSE image.

2.1 Working with speckles

With regard to the recombination of multi-modal data, one first intuitive approach is to use sets of matching points over the area of interest in complementary images. Some methods aiming at correcting the distortions in EBSD maps have been implemented using this principle [15]. They require the manual selection of a set of reference and corresponding points in the BSE image and the EBSD data. The higher the number of pairs of points, the greater the precision. However, this step can be time-consuming and does or made difficult by the lack of contrasted features in both EBSD and BSE data. Furthermore, it does not scale up to large or tridimensional datasets, as it may require human annotations for every slice in the dataset. In most 3D EBSD datasets, the process of data collection consists in successive movements of the stage to positions that are practical for machining, imaging, and collection of diffraction patterns. The placement of the stage from a slice to another is never exactly the same [17]. This modifies the set of reference points on every slice.

Here, we propose to use speckles of similar features without any human annotations, allowing the recombination of multi-modal data to scale up to large EBSD datasets. To do so, a high number of automatically generated reference points are combined by using microstructural features as speckles, such as pores, a second-phase, or precipitates. Since the goodness of the superimposition of the speckles can be quantified, it can be used as a parameter to maximize, enabling an automated matching process. Therefore, quantifying the goodness of the superposition of the speckles is critical.

Defining a similarity measure between the speckles

In order to superpose the EBSD and BSE speckles, a quantitative measurement of their similarity is needed. The Dice similarity metric [18] (also known as F1-score) describes the relative overlap between the segmented images. More precisely, this score measures the overlap between segmented pixels. Formally, given the binary EBSD speckle $\mathcal{I}^{EBSD}\in\{0,1\}^{I\times J}$ and the binary BSE image $\mathcal{I}^{BSE}\in\{0,1\}^{I\times J}$ where $1$ encodes the segmented object and [math] the background image, and ${I\times J}$ their domain of definition, the Dice similarity for is computed as follows:

[TABLE]

where $|.|$ encodes the size of the underlying set and $\mathcal{I}_{seg}^{.}$ only corresponds to the segmented pixels in the image. Thus, $|\mathcal{I}^{EBSD}_{seg}\cap\mathcal{I}^{BSE}_{seg}|$ corresponds to the number of matching pixels and $|\mathcal{I}^{EBSD}_{seg}|$ and $|\mathcal{I}^{BSE}_{seg}|$ respectively correspond to the total number of pixels in the EBSD and BSE images. This similarity measure counts the number of overlapping segmented pixels normalized by the sum of segmented pixels. It ranges from 0 (no overlapping segmented pixels) to 1 (perfect matching). The value of the ${Similarityfunction}(\mathcal{I}^{EBSD},\mathcal{I}^{BSE})$ is referred to as the ’score’ in the following. As the Dice similarity is both intuitive, robust and easy to compute, it has been extensively used in several image segmentation applications such as medicine [19].

2.2 Correction work-flow

The global process of the reconstruction is shown in figure 1. Starting from the initial EBSD data and the BSE image, two speckles are generated. In the first step (initial alignment), the BSE speckle is rescaled so its pixel size corresponds to that of the EBSD data. The rescaling is done using a nearest neighbor interpolation. The BSE speckle image is also rotated and translated so it is pre-aligned to the EBSD data. To do so, a grid-search over a set of affine transformations is performed, to distort the re-sized speckle over three parameters: two translations along the x and y directions, $t_{x},t_{y}$, and a rotation of angle $\alpha$. The user is free to choose the widths of the translation and the rotation angle, as well as their increment size. The minimum increment for the translations is one pixel. The step size for the angle can be refined as much as wanted. All the possible translations and rotations within the provided ranges are tested. The transformation leading to the best score is saved and applied to the BSE speckle. The latter is then used in the second step, as the reference speckle for the compensation of the distortions. In the second step (CMA-ES on fig. 1), the EBSD speckle is meshed on a regular grid, and the grid is distorted using the CMA-ES optimizer, in an iterative process. The initial and new locations of the points of the grid are used as pairs of control points, to estimate the distortion function. At each iteration, the score produced by the distorted EBSD speckle is calculated. The distortion function corresponding to the distorted speckle that exhibits the best score is used to determine the new location of the points in the new EBSD file and the segmented BSE image is used to fill in the phase data (if needed). The parametrization and calculation of the distortion function are explained in the following.

2.3 The distortion function

Formally, the goal is to find the distortion function $f$ that maximizes the similarity between the EBSD and the BSE speckles. This leads to the optimization procedure described in eq. 3:

[TABLE]

As described in section 1, many overlapping distortion phenomena occur while collecting the EBSD data, leading to a non-linear distortion. In the following, the nature of the distortion function is assumed to be polynomial. It is a very weak assumption as polynomial functions can approximate a large range of complex functions. Besides, the proposed method can be easily extended to other function approximators such as Gaussian processes, trees etc. Formally, a polynomial distortion of degree $P$ is defined as follows:

[TABLE]

For instance, given a polynomial fit of degree 2, the coordinates $(x,y)$ of a pixel are transformed into $(x^{\prime},y^{\prime})$:

[TABLE]

Given a polynomial order, the goal is to find the optimal set of weights $c^{*}$ that maximizes the similarity measure. This set of weights $c$ can be computed by two means: by directly looking into the space of parameters $c$ or by generating intermediate matching points that are used to perform a polynomial regression.

While the first approach may be more direct, the manifold of parameters $c$ is very difficult to explore as a small change of parameters can lead to drastically different distortions. Therefore, in the proposed method, a set of matching points is generated, that is used to regress the polynomial coefficients $c$. As a small change in the location of the matching points lead to a close distortion, this process allows a better granularity in exploring the space of the distortions. More precisely, a regular grid and a distorted one are generated, and then the coefficients $c$ are determined by a least-squares fitting method that maps the distorted grid to the regular one. This function is described below:

Generate a distorted mesh grid $M=((x_{0},y_{0}),\ldots,(x_{N},y_{N}))$

• 2.

Generate a regular mesh grid $M^{\prime}=((x_{0}^{\prime},y_{0}^{\prime}),\ldots,(x_{N}^{\prime},y_{N},))$

• 3.

Find the distortion $f$ parametrized by $c$ such as c=\mathop{\mathrm{argmin}}_{\hat{c}}\sum_{i}^{N}||\bigl{(}\begin{smallmatrix}x_{i}^{\prime}\\ y_{i}^{\prime}\end{smallmatrix}\bigr{)}-f_{\hat{c}}(x_{i},y_{i})||_{2}

where, $||.||$ is the euclidean norm. Note that the dimension of the mesh grid is tuned by the user via the number of points $N$ and/or the step-size between points. The influence of that parameter is discussed in section 4. As a result, finding the optimal distortion $f^{*}$ requires to find the optimal set of matching points $M^{*}$ that better describes the underlying distortion. These matching points are then used to estimate the distortion itself. Counter-intuitively, the distorted grid may not match physical reality as discussed further in section  4.4.

Finally, the polynomial distortion function is called from the python package Sklearn [20] and Skimage [21]. The next section describes how the pairs of matching points are created, using the CMA-ES optimizer.

2.4 Using CMA-ES as a randomized black-box optimization to correct the distortions

The goal of the optimization is to find the function $f$ that maximizes the $Similarity$ function between the speckles by finding the best distorted mesh grid $M$ (which is non-differentiable). To do so, a black-box optimizer method is used. Black-box optimization algorithms are efficient at solving non-linear, non-convex optimization problems, in continuous domains. They also overcome the deficiencies of the derivative-based methods in complex multidimensional landscapes that are rugged, noisy, have outliers or local optima, or when no error-gradient is available. All those features make the black-box optimizer relevant for the present problem, where the shape of the distortion is unknown and the process of superposing speckles intrinsically implies the $Similarity$ function to reach local maxima. The influence of noise present on the EBSD speckle is discussed in section 4.5.

The black-box algorithm that is used in the present case is CMA-ES, which stands for Covariance Matrix Adaptation Evolutionary Strategy. CMA-ES has become a standard tool for continuous optimization and has been applied in various fields of research, such as image recognition for biology [22, 23], energy [24], chemistry [25, 26]. Though, to the authors knowledge, it has not been applied in the field of electron microscopy. CMA-ES belongs to the family of evolutionary strategies (ESs), they are iterative algorithms, based on the principles of natural selection. CMA-ES involves a parametrized distribution (a multivariate normal distribution) that evolves throughout the iterations [27]. ESs follow four distinct steps: initialization, sampling, evaluation and update. The initialization step consists in initializing the probability distribution parameters and sampling an initial population of individuals accordingly. At each iteration, a new generation of individuals is created. Those individuals are candidate solutions to the optimization problem, which goodness can be evaluated on a fitness function. After evaluation of the fitness of each individual in that population, the statistics of the distribution are updated according to the algorithm at hand. A new population is re-sampled according the new distribution and the process is repeated until a termination criterion is reached. Those operations are described below:

Initialize the distribution parameters $\theta$

For each generation $g=0,1,...number\ of\ iterations$:

Sample $\lambda$ individuals $x_{1},x_{2},...,x_{\lambda}\in\mathbb{R}^{n}$ according the probability distribution $P_{\theta}(x)$

Evaluate them on the fitness function

Update $\theta\leftarrow F(\theta,x_{1},x_{2},...,x_{\lambda},Similarity(x_{1}),Similarity(x_{2}),...,Similarity(x_{\lambda}))$

where $P_{\theta}$ is a probability distribution that describes where the good solutions are believed to be and $F(.)$ is some update rule. In the CMA-ES, $P_{\theta}$ is a multivariate normal distribution [27]. Those steps are detailed in the following.

Initialization

The distribution $P_{\theta}$ encodes the distribution of the distorted spatial mesh grid. More precisely, the distribution is a multi-variate normal distribution where each dimension encodes either the x-coordinate or y-coordinate of a given point of the grid. For example, if the mesh is a 4x4 grid, the multi-variate normal distribution has 32 dimensions (16 points of 2 coordinates each). The initial distribution is defined by centering the mesh distribution on a regular grid. The standard deviation then encodes the initial acceptable moving distance of the matching points. Both the spacing of the points of the grid and the standard deviation are chosen by the user. In practise, CMA-ES slowly and iteratively distorts the regular grid, in order to retrieve a set of control (matching) points that encode the distortion. The concept of mesh grid here is to understand as an array of points that do not have any topological connectivity. In other words, the permutation of two points of the grid is allowed and does not affect the regression of the distortion function. The properties of the mesh are discussed in more detail in section 4.4.

Sampling

New distorted meshes are generated following the distribution $P_{\theta}$. The coordinates are rounded, and matching points outside the range of the speckle dimension are kept.

Evaluation

The fitness function is the $Similarity$ described in eq. 1. The coefficients $c$ of the distortion function are first regressed by matching the distorted mesh grid to the regular one. The similarity measure is then computed and returned as the fitness score.

Update

The basic CMA equation for sampling the search points at the generation $g+1$ is given in eq. 6.

[TABLE]

where $x_{k}^{(g+1)}\in\mathbb{R}^{n}$ is the $k-th$ offspring of the generation $g+1$. $m^{g}\in\mathbb{R}^{n}$ is the mean value of the distribution at the generation $g$. $\sigma^{(g)}\in\mathbb{R}^{+}$ is the step size at the generation $g$. $\mathcal{N}(0,\mathcal{C}^{(g)})$ is a multivariate normal distribution with a mean of zero and covariance matrix $\mathcal{C}^{g}$. $\mathcal{C}^{g}\in\mathbb{R}^{nxn}$ is the covariance matrix of the distribution at the generation $g$. It is symmetric definite positive and describes the geometrical shape of the distribution. The initial value of $\sigma$, $\sigma_{0}$, is picked by the user and $\mathcal{C}_{0}=\mathcal{I}$. They both evolve throughout the iterations, as the population evolves. The self-adaptation of those parameters is the key point for the rapid convergence of the optimization[28]. The equations governing the update of $\sigma^{(g)}$, $\mathcal{C}^{(g)}$ and $m^{g}$ are described in A. In the present case, this update shifts the spatial distribution of the points of the distorted grid in order to search for the optimal distortion function.

Termination

This process is repeated until the termination criterion is reached. In the present case, a criterion based on the maximum number of steps is used. More advanced methods can be used to detect when the optimization (and thus the fitness score) starts plateauing. If several speckles are known to have the same distortion, they also can be used as a validation criterion to avoid over-fitting. Once the CMA-ES algorithm is finished, the means of the distribution $P_{\theta}$ are used as the final distorted mesh grid.

2.5 Generation of the new EBSD map

Once the maximum number of calls has been reached, the polynomial function leading to the best superposition of the speckles is applied to the EBSD data. In this process, the original grid of the EBSD is kept and the Euler angles, confidence indexes and other data is interpolated using the polynomial function. The points of the grid that end up containing no data are filled with zero-values for the confidence indexes and the $(0,0,0)$ triplet of Euler angles. The rotated and translated segmented BSE image is used to fill-in the phase data of the new EBSD file. For that reason, and as already mentioned, the segmentation of the BSE image has to be done as accurately as possible. Some application examples are shown in section 3.

2.6 Source code and data availability

The full distortion pipeline was developed in Python. The code source, the hyperparameters and one data set are available at https://github.com/MLmicroscopy.

3 Results

All materials have been characterized in a FEI Versa 3D SEM, equipped with a TSL EDAX EBSD Hikari Plus Camera with an indexing speed of 300 to 400 frames per second and a 4x4 binning. The computer for the distortion correction was equipped with a Quad Core processor, 4.2 GHz and 64 GB of RAM. Most computations required about 10 to 15 minutes.

3.1 $\beta$ phase in $\alpha-\beta$ Titanium alloy

A sample of Ti-6Al-4V with an $\alpha-\beta$ structure has been embedded in bakelite and electro-chemically polished. It was then mounted on the SEM stage and copper tape was used to enable electrical conductivity from the sample holder to the polished surface. An area of 14 x 20 $\mu m$ has been characterized with a 20 kV acceleration voltage and a step size of 40 nm. One purpose of collecting this dataset is to better resolve the phase map, using the better quality of the segmentation of the BSE image, with fine details. This mounting system was purposely used in order to induce a large drift of the beam during the EBSD scan.

The orientation map for this sample is shown on fig. 2-a. The corresponding BSE image is shown on fig. 2-d, where the $\beta$ phase exhibits a lighter contrast than the $\alpha$ phase. Fig. 2-b shows the initial phase map with the $\beta$ phase in white. This phase map is superimposed to the Index Quality map on fig. 2-c, where the $\beta$ phase is colored in red. The comparison of fig. 2-b and -e shows that only a small fraction of the $\beta$ phase was indeed identified as such on the initial EBSD data. Its initial area fraction is 3.4%, versus 9.9% in the segmented BSE image. The CMA-ES optimizer has been applied on a 25 x 25 points mesh grid, with an initial standard deviation $\sigma_{0}=20$ pixels, using a polynomial order of 3 and 5,000 iterations. Fig. 2-f shows the final phase map, colored as fig. 2-c. It superimposes well with the Index Quality map. The final area fraction of $\beta$ phase is 9.9%, as in the BSE speckle. The contribution of the CMA-ES optimizer to compensate higher order distortions and match all the features of the speckles is discussed in section 4.1.

3.2 $\gamma$ and $\gamma^{\prime}$ phases in Rene 65 superalloy

Rene 65 is a polycrystalline Nickel-based alloy that has been designed for turbine disk applications. Its chemical composition (wt%) is Ni-16 Cr – 13 Co – 3.7 Ti – 2.1 Al – 4 Mo – 4 W – 1 Fe – 0.7 Nb – 0.05 Zr – 0.016 B [29]. The sample was extracted from a ring that was subjected to a sequence of forging operations followed by solution annealing and aging treatments, 1065$\degree$C for one hour and 760$\degree$C for 8 hours respectively. It was prepared using conventional metallography techniques, followed by vibratory polishing using a 0.04 $\mu m$ $Al_{2}O_{3}$ suspension. Its microstructure consists of fine $\gamma$ matrix grains, which average equivalent diameter is 10-12 $\mu m$. Spherical primary $\gamma^{\prime}$ precipitates are located on the grain boundaries, with an equivalent diameter in the range of 1-4 $\mu m$. EBSD maps and corresponding BSE images have been acquired using an acceleration voltage of 20 kV with a step size of 0.1 $\mu m$ over an area of 150 x 200 $\mu m$. Figure 3-a shows the initial EBSD map colored according to the orientation of the crystals projected along the normal to the polished surface. Figure 3-b shows the corresponding BSE image on which the precipitates exhibit a darker contrast. Figure 3-c and -d show the BSE and EBSD speckles, respectively, created by segmenting the $\gamma^{\prime}$ precipitates. Fig. 3-c has been created by segmenting the features which equivalent diameter is smaller than 5 $\mu m$, considering the annealing twin boundaries as internal defects. This criterion enables a rough segmentation of the precipitates. Note that some precipitates may have been omitted from the segmentation and some small grains may have been mistakenly segmented as precipitates. The influence of the quality of the EBSD speckle on the goodness of the compensation of the distortions is discussed in section 4.5. This segmented image was then used as the EBSD speckle for the CMA-ES optimizer, along with a segmented BSE image of the same area. The CMA-ES optimizer was used on a mesh grid with a spacing of 24 pixels between consecutive points, an initial standard deviation of 5 pixels and a polynomial order of 3. The initial grid is presented on figure 3-e, the final grid is presented on figure 3-f. The superimposed speckles after applying the CMA-ES optimizer are shown on fig. 3-g. The distortion correction leads to a final score of 0.62 after 1,800 iterations. Figure 3-h shows the final phase map with the precipitates colored in red, superimposed with the Image Quality map displayed in grey scale. The location of the precipitates corresponds well with the location of the grain boundaries on the Image Quality map, over most of the area of interest. The precision of the method is discussed in sections 4.6 and 4.4.

3.3 Pores in Additive Manufactured Inconel 718

The alloy Inconel 718 has been manufactured by the Electron Beam Melting (EBM) technique, as described in [30]. A sample was cut by Electrical Discharge Machining (EDM) and polished parallel the build direction, using conventional metallography processes followed by 0.05 $\mu m$ colloidal silica. Image and EBSD data have been acquired under a 20 kV acceleration voltage, using a 0.5 $\mu m$ step size over an area of 300 x 700 $\mu m$. Fig. 4-a shows the Secondary Electrons (SE) image that was used for the SE speckle. The beam contamination zone corresponding to the area of the EBSD map exhibits a trapezoid shape. A standard cleaning procedure involving slight cleaning of the EBSD data has been applied, consisting in removing isolated groups of less than 3 pixels having the same orientation, leading to a change of 1% of the data. Since the pores could not be indexed, the TSL software assigned random orientations to those pixels. The inverse pole figure map, projected along the build direction, is shown on fig. 4-b. Some of the grains have their $<001>$ axis preferentially aligned with the build direction (vertical), a common feature of additive manufactured face-centered cubic materials. A partial EBSD speckle of the pores has been created by segmenting the pixels of lowest Confidence Indexes (CI). On the other hand, the pores have been accurately segmented from a SE image of the same area, in order to create the SE speckle. The CMA-ES strategy has been applied, using a 65 x 65 points mesh grid, an initial standard deviation $\sigma_{0}=20$ pixels, a maximum sampling of 5,000 iterations and a polynomial order of 3 for the distortion function. Fig. 4-c shows the superimposed speckles after the CMA-ES strategy. The EBSD speckle consists of light blue points on a transparent background, and the SE speckle consists of dark red points on a light red background. Thus, as shown in the insert on the figure: light red pixels (labeled ”1” on the insert) correspond to the SE background, dark red pixels (labeled ”2”) correspond to pores of the SE speckle that did not match any pores of the EBSD speckle, light blue pixels (labeled ”3”) are pores of the EBSD speckle that did not match the SE speckle dark blue pixels (labeled ”4”) are pores that were successfully matched. 111using a method based on the CI to segment the pores can lead to inaccuracies in the segmentation: some pixels may have been segmented as pores even though they are not (typically those pixels correspond to grain boundaries). From this basis, the final EBSD map is shown on fig. 4-d. Only the points with a confidence index higher than 0.2 are plotted. The color code is the same as fig. 4-b and the pores are colored in black. Most of those pores are located on the grain and sub-grain boundaries. After applying the distortion, some data points ended up out of the initial grid, which leads to the black pixels on the edges of fig. 4-d. This dataset illustrates the limits of this method to compensate the distortions over broad areas when the features are unevenly distributed over the area.

4 Discussion

4.1 Versatility of the method

Several methods have been proposed in the literature, for correcting the distortions in EBSD data, using various algorithms to make up for the drift distortions. Zhang et al. [15] have shown that a thin plate spline function enables compensation of the distortions in EBSD data. Contrary to most methods, the use of the CMA-ES strategy does not assume any specific shape of the distortion function, nor constrain the order of the polynomial function. The use of speckles and a score to quantify the goodness of the superposition enables the matching of both images at the resolution of the pixel. Fig. 5 shows a comparison between the superposition of the speckles in the Ti6Al4V dataset, after the compensation of affine distortions only (a) and after the CMA-ES procedure (b). On a blue background, the EBSD speckle is colored in red and the BSE speckle in white. Compensating the affine distortions enables a good overlap of the speckles. The CMA-ES optimizer then compensates the finer distortions, even very local, and enables to achieve a better matching of the speckles despite an important drift of the beam. Three examples are presented in the inserts of Fig. 5, where the laths are accurately superimposed.

Figure 6 shows the maps of the vertical and horizontal components of the distortion function in the Inconel 718 material. The origin and x and y directions of the scan are indicated with arrows. The electron beam scanning pattern starts from the top left corner and goes line-by-line from top to bottom. A polynomial order of 3 has been used. Both maps show that the maximum error occurs on the first points of the scan. Indeed, drift phenomena are the highest in the beginning of the scan. The magnitude of both components then decrease along both x and y directions. Along the x direction, $f_{x}$ decreases relatively more than $f_{y}$, which means that the error on the x’ coordinate decreases along this direction, whereas the error on the y’ coordinate remains relatively constant. From a line to another (along the y direction), the magnitude of $f_{x}$ decreases, which means that the error on the location of the x’ coordinate decreases from a line to another. A similar observation can be made about the vertical component $f_{y}$. Those trends are consistent with most drift phenomena, where the first scanned pixels that are the most subjected to drift. Drift phenomena tend to decrease as the scan proceeds: along a line, and from a line to another. A quantitative analysis of those trends can be found in B.

4.2 Sensitivity to the choice of the initial parameters in CMA-ES

The influence of the CMA-ES parameters on the goodness of the matching has been studied in the Ti-6Al-4V dataset. Figure 7 shows the final scores obtained for different numbers of points per width and height of the image, and initial step sizes. The obtained surface is convex, with a global maximum reached for (35,21), pointed by the red arrow on the figure. For a fixed number of grid points, the score admits a maximum for a given initial step size $\sigma_{0}$. This behavior is visible on the four points on the far right of the figure, on the part of the envelope highlighted by a blue curve. For 55 grid points per length of the image, the score reaches a maximum of 0.5128 for an initial step size of 27 pixels. A similar behavior occurs for any given initial step size but varying numbers of points in the mesh grid. Note that a similar behavior is expected for all data sets, but for different values of the parameters.

4.3 Precision of the reconstruction

The accuracy of the reconstruction of the phases can be quantified by comparing the location of the phase boundaries and grain boundaries on the reconstructed maps. Such measurements have been performed on the Rene 65 dataset and displayed on fig. 8. A subset of the corrected EBSD map has been selected where three line profiles have been traced, crossing isolated or clusters of $\gamma^{\prime}$ precipitates. The point-to-point misorientation is plotted in black and superimposed to the phase data, which is plotted in red. The first profile crosses a single precipitate. The phase boundaries and grain boundaries ”a” and ”b” are indicated on the corresponding profile. The transition from a phase to the other corresponds well with the crossing of high-angle boundaries. The second profile crosses a cluster of two precipitates which were separated on the segmentation of the BSE image. The phase boundaries ”c”, ”d” and ”f” are well superposed with the grain boundaries. Note that the boundary ”d” consists of a single pixel. The profile also crosses an annealing twin boundary ($\Sigma 3$, 60$\degree$ rotation around a $<111>$ axis) present in one of the precipitates at the marker ”e”, characterized by a 60$\degree$ misorientation. The third profile also crosses two precipitates which were not separated on the segmented BSE image. The grain boundaries of the first precipitate are labeled ”g” and ”h” on the profile. The second precipitate contains an annealing twin boundary ”i”, that is crossed right before the grain boundary ”j”. This grain boundary is well superimposed with the phase boundary. The error on the location of the phase boundaries has been measured on all the profiles, by comparing their location $x_{Ph.B}$ to that of the grain boundaries $x_{GB}$ through the parameter $\Delta x$ of eq. 7. The results are displayed in table 1. The error is less than 0.2 $\mu m$ for all the boundaries (i.e. 2 pixels), reaching $\approx$ 0 $\mu m$ for most of them. The error on the location of the two phase boundaries that were not exactly correlated with the grain boundaries (”f” and ”g”) corresponds to a misplacement of 1 and 2 pixels respectively. This gives an indication of the precision achievable with this correction method, assuming a good segmentation of the BSE image is achieved. A segmentation in which the boundaries of the objects have been eroded or blurred impacts the quality of the reconstruction.

[TABLE]

4.4 Convergence, stability and repeatability of the CMA-ES optimization

During the CMA-ES optimization procedure, new distorted mesh grids are generated to regress the polynomial distortion. Counter-intuitively, these grids are not constrained to encode a realistic distortion; they are only optimized such that the final distortion is meaningful. Thus, one must not use the distorted mesh grids outside the polynomial regression step, as several different mesh grids can lead to the same distortion function. Indeed, there exists an infinite number of sets of points that extrapolate to the same final function. Those points do not have to always be on the curve itself. Even though phenomena of grid points swapping does not affect the calculation of the distortion function, the authors suggest that the user sets the $\sigma_{0}$ parameter as less than the distance between two consecutive grid points, if the user wants to avoid this phenomenon. To a lesser extent, several polynomial coefficients may also lead to a very similar resulting distortion function, only differing around the border of the speckle.

Despite this apparent limitation, CMA-ES turns out to be highly reproducible in practice. For example, fig. 9 shows the typical pattern of convergence of the CMA-ES algorithm: fig. 9-a shows the evolution of the score as a function of the number of iterations in the CMA optimizer and 9-b shows the evolution of the step size throughout the iterations. Starting from a score value of about 0.4, the CMA optimizer enables a quick convergence towards a steady-state value of about 0.7. This steady state is reached after less than 5000 iterations. As the score increases, the step-size $\sigma$ decreases with the number of iterations, indicating little variation around the mesh grid ground.

The CMA-ES strategy has been applied 100 times on the Rene 65 dataset. Each run consisted in 2000 iterations, with a step size of 75 points, an initial standard-deviation $\sigma_{0}$ of 5 pixels. A polynomial order of 3 was used for the polynomial function. The CMA optimizer produced an average score of 0.6587 with an associated standard deviation of 0.0015. The similarity score between distorted speckles over the 100 runs was 0.9353 in average, indicating a good precision and repeatability of the process. Figure 10-a shows the evolution of the score throughout the iterations: the dark blue curve corresponds to the mean score over 100 runs. The light blue colored area around the curve corresponds to the lower and upper bounds of the score (mean - standard deviation, and mean + standard deviation, respectively). Those values show that the optimization is stable and repeatable over several runs. Some steps are clearly visible on the first 500 iterations and correspond to the new bounds determined by the creation of new individuals at each generation in CMA-ES. Based on those statistics, a heat map has also been generated, displayed on fig. 10-b. This map shows the segmented BSE speckle colored according to the number fraction of times that a pixel was present at a given (x’, y’) location. In other words, a pixel that appears in white (value 1.0) was always assigned to the $\gamma^{\prime}$ phase. On the opposite, a pixel that appears in black was never assigned to this phase. The more consistent and precise the optimization, the sharper the contrast on this map. The shape of the precipitates appears clearly on fig. 10-b, which is consistent with the good repeatability suggested by fig. 10-a. This map also gives information about the precision of the reconstruction over the whole map. The center of the precipitates usually appears with a value of 1.0, however their boundaries are usually less precise. The deviation on the reconstruction of the boundaries of the precipitates is illustrated on two examples, labeled ”A” and ”B”. The precipitate ”A” is close to the border of the image, and the precipitate ”B” is in the center. The inserts on the side of the figure show that the boundaries of those two precipitates are not reconstructed with the same consistency over the 100 slices. There is only one layer of pixels having a value lower than 1.0 on the precipitate B, versus 2 to 3 rows on the boundary of the precipitate A. The boundaries of most of the precipitates located on the edges of the map have a similar coloring. This indicates that the CMA-ES optimizer leads to consistent results in the center of the map, with a deviation of about 1 pixel as the location of the phase boundaries are approached. At the borders of the map, the consistency on the location of the phase map is about 2 to 3 pixels. This corresponds, in this dataset, to an associated error of 0.2 to 0.3 $\mu m$ on the location of the phase boundaries. This error remains much smaller than the actual size of the features of interest. In other words, the reconstruction does not artificially create additional features (precipitates) nor assign the wrong phase to any feature.

Again, it should also be noted that different final grids can describe the same polynomial distortion function, since the calculation of the function is not affected by the permutation of the points of the mesh grid. Moreover, different polynomial functions can describe the same physical distortion. This phenomenon is enhanced if the points of the speckle are sparsely distributed or mostly located in the center of the region of interest.

4.5 Sensitivity to the accuracy of the EBSD speckle

As mentioned previously, the use of speckles eliminates the manual selection of reference points. Another advantage of using speckles and a score based on the number of superimposed pixels, is that the EBSD speckle can be only a small subset of the original image. This is the case in the Rene 65 dataset, where the BSE speckle contains all the $\gamma^{\prime}$ precipitates but the EBSD speckle contains the smallest features in the dataset. Some precipitates are missing on this speckle, and small grains are mistakenly segmented as precipitates. This does deteriorate the theoretical maximum score achievable, but does not prevent a good superposition of the speckles, as already shown on fig. 3-c. The influence of the quality of the EBSD speckle has been studied on the Inconel 718 dataset: some noise was added (wrong features and random noise) and some pores were also removed from the EBSD speckle on purpose, in comparison with the one used on 3.3 (fig. 4). The amount of features added or deleted is quantified by the $\Delta f$ parameter, which quantifies the percentage of pixels points added or deleted to the speckle containing the correct number of pores: $\Delta f=\frac{numberofpixelsinwrongspeckle}{numberofpixelsinrightspeckle}$. Figure 11 shows the best score obtained as a function of the percentage of pixels added or deleted in the EBSD speckle. Empty and filled markers correspond to the initial and final scores in CMA-ES, respectively. The speckle at $\Delta f=0\%$ corresponds to the EBSD pattern in which only the pores were segmented. It leads to the best score achievable. The deletion of points leads to a drastic decrease of the score. The addition of ”wrong” cavities, however, leads to a slight decrease of the score but stagnates to a roughly constant value, as the number of points added increases. This indicates that the calculation of the distortions is fairly tolerant to the presence of noise or wrong elements in the EBSD speckle. However, missing points are more detrimental.

4.6 The case of $\gamma-\gamma^{\prime}$ nickel-base superalloys: comparison with other mapping techniques

The $\gamma$ and $\gamma^{\prime}$ phases of nickel-based superalloys are well known to be challenging to differentiate, on the single basis of their diffraction patterns. Methods like conventional EBSD with an indexation of phases based on the Hough transform are not efficient. Other indexation techniques have been proposed, such as the dictionary approach by De Graef [31]. In this approach, the collected patterns are compared to a dictionary of synthetically generated patterns. This technique enables the indexation of patterns with a better angular resolution [32] as well as the differentiation of phases in complex alloys [33] and geological materials [34]. However, it can be computationally expensive and -this far- has not been able to differentiate $\gamma$ and $\gamma^{\prime}$ phases. Combined EDS-EBSD methods have been used in the last years but suffer from long indexing times (and so, increased beam drifting), post-processing times and limited resolution [12]. Using a completely different technique but a similar approach to dictionary indexing, with ion channeling imaging, the iCHORD method of Langlois et al. [35] has been shown to enable the differentiation of phases [36]. The channeling contrast is used to obtain information about the orientation. On the other hand, the secondary ion signals from $\gamma$ and $\gamma^{\prime}$ phases are sufficiently contrasted to enable the identification of phases. In a recent comparative study, Vernier et al. [37] have shown that this method enables a more accurate mapping than combined EDX-EBSD methods, since it enables to identify precipitates of size down to 150 nm. The scanning speed is also 3 times faster than EDS-EBSD methods. The time necessary for the acquisition and indexation of the patterns however, is typically similar than that of the combined EDS-EBSD method. With a post-processing time of about 15 to 30 minutes, the multimodal data recombination process using the CMA-ES optimizer appears to be a competitive alternative to the other methods proposed in the literature. The only constraint is the acquisition of a BSE image of sufficient quality (in terms of contrast), and its precise segmentation. Indeed, in this method, the size of the smallest objects resolved depends not only on the step size of the EBSD map but also directly on the quality of the segmentation of the BSE image. Assuming a good segmentation, the spatial resolution of this technique is about a few pixels, which is better than EDS-EBSD methods, and comparable to the iCHORD method. It can also be easily implemented in any SEM equipped with a conventional EBSD detector.

5 Conclusion and perspectives

A new method for the compensation of the distortions and improved phase differentiation in EBSD data has been developed. The principle consists in using an electron image of the same area than the EBSD data, taken at a 0$\degree$ tilt angle, and using it as a reference. Similar features are segmented out of the EBSD data and electron image. Then, the CMA Evolutionary Strategy is applied in order to match the speckles. The goodness of the superposition is measured by a score, which ranges from 0 to 1. Crystallographic and phase data are then recombined in a new EBSD file. This new file has the same grid than the initial EBSD file.

This method has been applied successfully to nickel-based superalloys and a titanium alloy. Assuming a precise segmentation of the phases on the electron image, this method can reach a precision of a couple of pixels over broad areas, despite important drift phenomena. It can be used on any EBSD dataset, as long as two speckles of similar features can be generated. It is relatively tolerant to inaccuracies and noise in the EBSD speckle. Full automation and short computation times make it a competitive post-processing method, when compared to other options proposed in the literature.

Acknowledgements

The authors gratefully acknowledge a Vannevar Bush Fellowship, ONR Grant N00014-18-1-3031.

We acknowledge the following agencies for research funding and computing support: CHISTERA IGLU and CPER Nord-Pas de Calais/FEDER DATA Advanced data science and technologies 2015-2020 Henry Proudhon from Mines ParisTech is acknowledged for prolific conversations and exchange of ideas. Olivier Pietquin, Jérémie Mary and Philippe Preux from Inria are acknowledged for insightful conversations. Xuedong Shang from INRIA is acknowledged for his various comments on CMA-ES. McLean Echlin from the University of Santa Barbara is acknowledged for providing the Ti-6Al-4V sample. Mickael Kirka from the Oak Ridge National Laboratory is acknowledged for providing the additive manufactured Inconel 718 sample. The Carlton Forge Works company (PCC Corporation) is acknowledged for providing the Rene 65 material.