Microtexture region segmentation of eddy current testing data using a structural prior

Crystallographic orientation refers to the set of rotations that generates the orientation of a grain with respect to a global frame of reference. In this work, it is described in terms of three Euler angles $(\psi_1,\theta, \psi_2)$: ψ₁ is rotation about the z-axis, θ the subsequent rotation about the new x-axis, and ψ₂ the final rotation about the new z-axis [27]. Generally, $\psi_1 \in [0, 2\pi]$, $\theta \in [0, \pi]$, and $\psi_2 \in [0, 2\pi]$. However, symmetry of the HCP crystal permits rotations to be expressed in the ranges $\psi_1 \in [0, 2\pi]$, $\theta \in [0, \pi/2]$, and $\psi_2 \in [0, \pi/3]$. Note that the axis of the HCP crystal parallel to the z-axis is called the c-axis, and the plane perpendicular to the z-axis is called the basal plane. A diagram of the unit cell for grains with hexagonal symmetry depicted in the three rotations that define crystallographic orientation is shown in figure 2.

Zoom In Zoom Out Reset image size

Ultrasound measurements are performed at spatial locations, in a grid pattern with 50 µm spacing, across the sample with a single element ultrasonic transducer. Both the sample and transducer are immersed in water to ensure proper coupling of acoustic energy. The transducer is an Olympus V3330 with a nominal frequency of 50 MHz, 6.35 mm diameter element, and a 5.08 mm focal length. The transducer also includes an integrated fused silica delay line that allows the ultrasonic wave to form a uniform wave front prior to focusing. The transducer is de-focused to slightly shift the time-of-flight of the reflected longitudinal wave from that of the Rayleigh surface wave. This is required to facilitate data analysis. Time-of-flight refers to the time required for the generated sound wave to propagate to and interact with the sample and then reflect back to the transducer to be measured.

A digitized voltage, directly proportional to the acoustic pressure across the entire transducer element, is recorded at each spatial location as a function of time. This is commonly referred to as an A-scan. A time gate is defined—in this case, a time window within which we expect the Rayleigh surface wave to propagate back to the transducer element. This time gate is allowed to shift in time as we scan across the sample by tracking the time-of-flight of the higher amplitude reflected longitudinal wave in a separate tracking gate. This accounts for any subtle local variations in sample geometry.

The maximum value of the voltage within the time gate for the Rayleigh surface wave is taken from the A-scan captured at each spatial location and used to produce the image shown in figure 3(b). The color scale in this image shows the measured maximum voltage within the gate, which is directly proportional to the amplitude of the Rayleigh surface wave, mapped across the surface of the sample. This form of presentation is commonly referred to as a C-scan.

The amplitude of this back-scattered Rayleigh surface wave has been found previously to be sensitive to grain structure, though the exact mechanisms are not fully understood [28]. Furthermore, as previously mentioned, since this is a single element transducer with a spherically focused beam, we are simultaneously measuring the back-scattered response as the wave propagates along the surface radially in all directions from the point of focus and therefore we do not obtain useful orientation information.

All ECT data used in this paper are simulated. The forward model used to simulate ECT data is the approximate impedance integral (AII) model (see [7]), which is based on Maxwell's equations. We briefly review its formulation here. In differential form, Ampere's Law and Faraday's Law are given by

$$\begin{align} \nabla \times \mathbf{H} & = j \omega \epsilon \mathbf{E} + \overline{\sigma} \cdot \mathbf{E} \end{align} \tag{ 1 }$$

$$\begin{align} \nabla \times \mathbf{E} & = -j \omega \mu \mathbf{H}, \end{align} \tag{ 2 }$$

where ω is the angular frequency, µ is the permeability, ε is the scalar permittivity, and $\overline{\sigma}$ is the anisotropic conductivity tensor. We assume there are no external fields and that the only current sources are the induced currents in the conductive material and the displacement currents. From [29], (1) and (2) can be combined to yield the expression

$$\begin{equation} \oint_S \left(\mathbf{E} \times \mathbf{H}\right) \cdot \mathbf{n} dS = \int_V j \omega \mu \mathbf{H} \cdot \mathbf{H} + j \omega \epsilon \mathbf{E} \cdot \mathbf{E} + \mathbf{E} \cdot \left(\overline{\sigma} \cdot \mathbf{E}\right) ] dV. \end{equation} \tag{ 3 }$$

where S is the surface enclosing a flaw in the conductive medium and V is the volume of the flaw. For our application, S is the entire surface of the sample. From [30], this expression can be related to the impedance of an ECT coil to yield

$$\begin{align} \Delta Z & = \frac{1}{I^2} \oint_S \left(\mathbf{E}_b \times \mathbf{H}_a - \mathbf{E}_a \times \mathbf{H}_b\right) \cdot \mathbf{n} dS \nonumber \\ & = \frac{1}{I^2} \int_V \left[ j \omega \mu \left(\mathbf{H}_a \cdot \mathbf{H}_a - \mathbf{H}_b \cdot \mathbf{H}_b\right) + j \omega \epsilon \left(\mathbf{E}_a \cdot \mathbf{E}_a - \mathbf{E}_b \cdot \mathbf{E}_b\right) \right. \nonumber \\ &\left. \quad + \mathbf{E}_a \cdot \left(\sigma_a \mathbf{I}_{3 \times 3} \cdot \mathbf{E}_a\right) - \mathbf{E}_b \cdot \left(\overline{\sigma_b} \cdot \mathbf{E}_b\right) \right] dV, \end{align} \tag{ 4 }$$

where all fields with the subscript a are the result of exciting the ECT coil above an isotropic homogeneous material with conductivity σ_a , while all fields with the b subscript are the result of exciting the coil above an anisotropic, polycrystalline material. Due to the relatively low conductivity changes in a rotated grain, the Born approximation can be applied to (4). Under this approximation, we have that $\mathbf{H}_a \approx \mathbf{H}_b$ and $\mathbf{E}_a \approx \mathbf{E}_b$, in which case (4) becomes

$$\begin{equation} \Delta Z = \frac{1}{I^2} \int_V {\mathbf{E}_a \cdot \left(\sigma_a \mathbf{I}_{3 \times 3} - \overline{\sigma}_b\right) \cdot \mathbf{E}_a} \; dV. \end{equation} \tag{ 5 }$$

The rotated, anisotropic conductivity tensor $\overline{\sigma}_b$ is a function of the local crystallographic orientation and is given by

$$\begin{equation} \overline{\sigma}_b = \left[ \begin{array}{ccc} \sigma_{11} & \sigma_{12} & \sigma_{13}\\ \sigma_{12} & \sigma_{22} & \sigma_{23} \\ \sigma_{13} & \sigma_{23} & \sigma_{33} \end{array} \right], \end{equation} \tag{ 6 }$$

where

$$\begin{equation} \begin{array}{cc} \sigma_{11} = \sigma_{xx} + \left(\sigma_{zz} - \sigma_{xx}\right)s_1^2 s_2^2 & \sigma_{12} = s_1 c_1 s_2^2 \left(\sigma_{xx} - \sigma_{zz}\right)/2 \\ \sigma_{13} = -s_1 c_2 s_2 \left(\sigma_{xx} - \sigma_{zz}\right) & \sigma_{22} = \sigma_{xx} + \left(\sigma_{zz} - \sigma_{xx}\right)s_2^2 + \left(\sigma_{xx} - \sigma_{zz}\right)s_1^2 s_2^2 \\ \sigma_{23} = c_1 c_2 s_2 \left(\sigma_{xx} - \sigma_{zz}\right) & \sigma_{33} = \sigma_{zz} +\left(\sigma_{xx} - \sigma_{zz} \right)s_2^2 \end{array} \end{equation} \tag{ 7 }$$

and c_i and s_i are the sine and cosine of the ith Euler angle, respectively. Within α-titanium, conductivity differs along the c-axis but is isotropic in the basal plane. ECT is sensitive to changes in the first two Euler angles because they correspond to changes in the orientation of the c-axis. However, it is not sensitive to changes in the third Euler angle due to the fact that it is simply a final rotation in the basal plane. The fact that $\overline{\sigma}_b$ in the final model is dependent on only ψ₁ and θ is consistent with the fact that ECT is not sensitive to changes in the third Euler angle within α-titanium. For further details on this model, we refer to ([7]).

Figures 3(c) and (d) show the imaginary portion of the simulated ECT data of the large grain specimen for two rotations of an elliptical ECT coil with axes equal to 200 µm and 635 µm. The simulation was a two-step process. First, the reference electric field E_a from (5) was computed using COMSOL Multiphysics. Specifically, a linear FEM simulation was used to compute the reference electric field values at the surface of an isotropic, homogeneous half-space. Second-order triangular elements with size equal to one-third of the skin depth were used with Lagrange interpolation. Given the field E_a , the change in impedance for a fixed location of the coil was computed by solving (5) using numerical integration in MATLAB. The conductivity tensor $\overline{\sigma_b}$ was computed at each point in the specimen using the orientations provided by the PLM data shown in figure 3(a). The spatial resolution of the PLM data was 15 µm, which is a sufficiently small step size to ensure the numerical integration converges. Note that the integration is performed over the span of the coil field, which is a 2 mm square centered at the coil location. The full scan of the specimen is performed by moving the center of the simulated coil 100 µm, interpolating E_a at the new location, and recomputing (5) over the new span of the coil.

Because this model only requires a single solve of an FEM model to obtain E_a , it is computationally efficient. A full raster scan of the specimen can be performed without needing to solve an FEM model at each measurement point. On an Intel Core i9 desktop with a 2.4 GHz CPU the simulation took around 230 s. Note that the full PLM dataset is spatially larger than the portion shown in figure 3(a). As such, we did not need to be concerned about edge effects in our simulation.

We chose to use simulated ECT data due to the fact that the method requires the ECT data and SAM data to be registered to each other, which is non-trivial. In this paper, we registered the SAM data to the PLM data using a landmark registration algorithm built into MATLAB. By simulating ECT data from the PLM data, we could ensure that it would be registered to the SAM data. The development of an algorithm to register SAM and ECT data is a topic of ongoing work, but is not the topic of this paper. While we did choose to use simulated ECT data, we note that the AII model was validated against experimental ECT data in [7]. Thus, we feel it is appropriate to use here.

Lastly, we use an elliptical absolute ECT coil which exhibits directional sensitivity to the underlying crystallographic orientation. The elliptical coil is necessary to recover information about both the first and second Euler angle. This is opposed to using a circular coil, which is only sensitive to the second Euler angle. We simulated ECT data from two coil rotations, both of which will be used in the inversion. As will be demonstrated in section 5, we require data from multiple rotations of the ECT coil to obtain unique orientation information.

Matching component analysis (MCA) is a data analysis method which is used to determine common features between two data domains by mapping each to a lower dimensional common domain [31]. It was originally developed for transfer learning; one of its initial applications was to find a common domain between experimental and simulated synthetic aperture radar (SAR) data to improve classification. We use an extension of MCA known as covariance generalized matching component analysis (CGMCA) [32]. For this application, the two data domains of interest are (1) the underlying segmentation of a known microstructure into MTR and (2) the corresponding ECT data of the microstructure.

Let $X_1 \in \mathbb{R}^{d_1}$ and $X_2 \in \mathbb{R}^{d_2}$ be random variables modeling the two data domains. Then CGMCA seeks to find maps g₁ and g₂ that are solutions to the constrained optimization problem

$$\begin{align} \text{minimize}\quad \quad &\mathbf{E} \left[\left\| g_1\left( X_1 \right) - g_2\left( X_2 \right) \right|^2 \right] \end{align} \tag{ 8 }$$

$$\begin{align} \text{subject to} \quad \quad &g_i: \mathbb{R}^{d_i} \rightarrow \mathbb{R}^k \end{align} \tag{ 9 }$$

$$\begin{align} & \mathbf{E}\left[g_i\left( X_i \right) \right] = 0, \end{align} \tag{ 10 }$$

$$\begin{align} & \mathbf{E} \left[g_i \left(X_i \right) g_i \left(X_i \right)^T \right] = \Gamma_i, \quad i \in \left\{1,2\right\}, \end{align} \tag{ 11 }$$

where $k \unicode{x2A7D} \min(d_1, d_2)$ is the dimension of the common domain, and the symmetric positive definite matrix $\Gamma_i \in \mathbb{R}^{k \times k}$ is the covariance matrix prescribed to the map g_i . It was shown in [32] that if we require that

$$\begin{equation} g_i\left(x\right) = A_ix + b_i, \end{equation} \tag{ 12 }$$

for $A_i \in \mathbb{R}^{k \times d_i}$ and $b \in \mathbb{R}^k$, then a closed form solution exists to the CGMCA optimization problem. Note that if the choice of $\Gamma_i$ in (11) is restricted to the identity matrix, this reduces to the original MCA problem. The algorithm to compute the CGMCA maps is provided in appendix.

The basic procedure for segmentation using CGMCA was developed in [33]. The CGMCA maps g₁ and g₂ are trained on a set of known microstructures with known MTR segmentation matched to the corresponding ECT data of each microstructure, as shown in figure 4. Given a new set of ECT data y with unknown segmentation, we first map y to the common domain using g₂. We then recover the unknown segmentation by inverting g₁ to map back to the segmentation domain. Formally, we wish to find the unknown segmentation $s^*$ such that

$$\begin{equation} s^* = \arg\min \left\|g_1\left(s\right) - g_2\left(y\right) \right\|_2^2. \end{equation} \tag{ 13 }$$

Recall that g₁ and g₂ are affine linear transformations; thus we can find $s^*$ by solving the linear system

$$\begin{equation} A_1 s = A_2y + b_2 - b_1. \end{equation} \tag{ 14 }$$

Note that $A_1 \in \mathbb{R}^{k \times d_1}$, with $k \lt d_1$; as such, this system is underdetermined. A possible solution is given by

$$\begin{equation} s^* = A_1^\dagger \left( A_2 y + b_2 - b_1 \right). \end{equation} \tag{ 15 }$$

However, as was shown in [33], this solution will not behave as a segmentation should, i.e. it will not be spatially piecewise constant. Thus, additional constraints on the solution are needed.

Zoom In Zoom Out Reset image size

In [33], the piecewise constant requirement was enforced by utilizing level set inversion to find $s^*$, while in [34], a penalty term was added to the inversion. Specifically, in [34], a regularization method developed for image deblurring in [35] was applied, which we review briefly in the next subsection. This paper proposes a new method to impose the piecewise constant requirement using the SAM data as structural information in the inversion, as will be discussed later in the section.

Rather than estimate s directly, the method proposed in [35] estimates the increments (i.e. the difference between neighboring pixels) of s. Denote the vector containing the horizontal and vertical increments of s by α. This method assumes that the entries of α are mutually independent and follow a conditionally Gaussian prior. The variance of each increment is also an unknown and is estimated as part of the inverse problem. More specifically, the method estimates the reciprocal of the variance of each increment, which we denote by λ. In [35], an exponential hyperprior is applied to each entry of λ to encode the idea that the majority of the increments will be close to zero, while still allowing for jumps when needed.

We briefly review the implementation of the image deblurring prior here; for further detail, we direct the reader to [35]. The biggest challenge in implementing this prior is writing the forward model in terms of α. Let N_e denote the number of increments in s. Note that we assume in the following calculations that $s \in \mathbb{R}^{N \times N}$ has been stacked into a single column vector of length N². We can compute α from s,

$$\begin{equation} Bs = \alpha, \end{equation} \tag{ 16 }$$

where $B \in \mathbb{R}^{N_e \times N^2}$ is a difference operator that computes both the horizontal and vertical increments of s. Because $N_e \gt N^2$, we cannot invert B to obtain s from α. It can be shown that $\text{rank}(B) = N^2 - 1$. Denote by $B_2 \in \mathbb{R}^{N^2-1 \times N^2}$ the matrix that contains a set of linearly independent rows of B, and let $\alpha_2 \in \mathbb{R}^{N^2-1}$ denote the vector that contains the corresponding edges.

The fact that the number of linearly independent edges is equal to $N^2-1$ agrees with the fact that an image cannot be reconstructed from increments alone; therefore, we must apply another constraint. Let $v \in \mathbb{R}^{N^2}$ be a vector with non-zero entries and let ν be the projection of $s$ onto $v$, that is, $\nu = v^Ts$. Define the matrix $C \in \mathbb{R}^{N^2 \times N^2}$ as $C = [B_2^T \; \; v]^T$, and then we can write

$$\begin{equation} C s = \left[\begin{array}{c} \alpha_2 \\ \nu \end{array} \right]. \end{equation} \tag{ 17 }$$

Because C is full rank, it can be inverted to find $s$ given the increments α ₂ and ν. This allows us to write the forward model in terms of the increments.

Next, define the augmented vector $\xi \in \mathbb{R}^{N_e + 1}$ and the permutation matrix $P \in \mathbb{R}^{N^2 \times N_e + 1}$ as

$$\begin{equation} \xi = \left[ \begin{array}{c} \alpha \\ \nu \end{array} \right], \quad P\xi = \left[ \begin{array}{c} \alpha_2 \\ \nu \end{array} \right]. \end{equation} \tag{ 18 }$$

Now we can write an equation to recover $s$ from ξ , specifically

$$\begin{equation} s = C^{-1}P\xi. \end{equation} \tag{ 19 }$$

In the inverse problem, there is an additional constraint placed on α ; specifically, that the four increments surrounding each interior point in $s$ sum to 0. This allows us to solve for the remaining entries of α in terms of α ₂ . (See section 4 in [35] for more details.) This condition is enforced by the matrix $M \in \mathbb{R}^{(N-1)^2 \times N_e+1}$, so that we require $M\xi = 0$. Lastly, define the matrix D as

$$\begin{equation} D = \left[ \begin{array}{c} A_1 C^{-1} P \\ M \end{array} \right]. \end{equation} \tag{ 20 }$$

Our forward model in terms of ξ is given by

$$\begin{equation} D \xi = \left[ \begin{array}{c} A_2y + b_2 - b_1 \\0 \end{array} \right] + g, \quad g \sim N\left(0, \gamma^2 \mathbf{I}\right). \end{equation} \tag{ 21 }$$

We have a Gaussian error term in our forward model to signify that the segmentation will not map exactly to the ECT data in the common domain. Additionally, this term expresses the fact that we require only that $M\xi \approx 0$, instead of $M\xi = 0$. We could have two different error terms for our CGMCA forward model and the additional constraint; however, for the sake of simplicity, we combine it into one.

We assume that the entries of the unknown ξ are mutually independent and follow a conditionally Gaussian prior. Furthermore, each entry of the reciprocal of the variance λ follows an exponential hyperprior. The optimal values of ξ and λ are given by the maximum a posteriori (MAP) estimates, which are solutions to the minimization problem

$$\begin{equation} \left(\xi, \lambda\right) = \arg \min \displaystyle{\frac{1}{\gamma^2}} \left\|D \xi - \left[ \begin{array}{c}A_2y + b_2 - b_1 \\ 0 \end{array}\right] \right\|_2^2 + \left\| \Lambda \xi \right\|_2^2 + \beta\sum_{k = 1}^{N_e + 1} \lambda_k - \sum_{k = 1}^{N_e + 1} \log \lambda_k. \end{equation} \tag{ 22 }$$

The entries of the diagonal matrix $\Lambda \in \mathbb{R}^{N_e + 1 \times N_e+1}$ are given by $\sqrt{\lambda}$ and β is the reciprocal of the mean of the exponential hyperprior.

The procedure to estimate ξ and λ is iterative, and requires alternating between optimizing ξ and λ . The procedure is as follows:

• (i)  
  Let h = 0. Initialize λ^h .
• (ii)  
  Update ξ by minimizing (22) with respect to ξ ; specifically, find ξ ^h by solving

  $$\begin{equation} \left[\begin{array}{c} D \\ \gamma \Lambda^h \end{array}\right]\xi = \left[ \begin{array}{c}A_2 y + b_2 - b_1 \\ 0 \end{array} \right] \end{equation} \tag{ 23 }$$

  in the least squares sense. Note that we first multiplied (22) through by γ² in order to obtain (23).
• (iii)  
  Update λ by minimizing (22) with respect to λ . The update for λ is explicit, and is given by

  $$\begin{equation} \lambda_\ell^{h+1} = \displaystyle \frac{1}{\left(\xi_\ell^h \right)^2 + \beta}. \end{equation} \tag{ 24 }$$
• (iv)  
  Increase h to h + 1 and repeat from step (ii) until convergence.

In this work, we run this algorithm for a fixed number of iterations. Then $s^* = C^{-1}P\xi^*$, where $\xi^*$ is the optimal value of ξ as found by the iterative procedure.

In this work, we propose applying a penalty term that incorporates structural information from the SAM data. We discuss the implementation of our method here. Denote the $\ell$th entry of the unknown discretized segmentation $s \in \mathbb{R}^{N \times N}$ as $s_\ell$ for $1 \unicode{x2A7D} \ell \unicode{x2A7D} N^2$ and denote the set of indices of the neighboring pixels to $s_\ell$ as $n_\ell$. We assume that $s_\ell$ is roughly equal to the weighted average of its neighbors, that is,

$$\begin{equation} s_{\ell} \approx \displaystyle{\frac{1}{\sum_{k \in n_\ell} \omega_\ell^k}}\sum_{k \in n_\ell} \omega_\ell^k s_k, \end{equation} \tag{ 25 }$$

where $\omega_\ell^k$ is the kth entry of the weight vector $\omega_\ell \in \mathbb{R}^{1 \times N^2}$. Note that the weight vector differs for each entry $s_\ell$.

In our application, we extract the MTR boundaries from the SAM data. This edge information is contained in the matrix e, which is the same size as the unknown segmentation s. The $\ell$th entry of e is defined as follows,

$$\begin{equation} e_\ell = \left\{\begin{array}{ll} 1, & \;\; s_\ell \text{lies on an edge} \\ 0, & \;\; \text{else} \end{array} \right. . \end{equation} \tag{ 26 }$$

Given the edges, the weight vector is formed in the following way. If $k \notin n_\ell$, then $\omega_\ell^k = 0$. Otherwise, $\omega_\ell^k$ is given by

$$\begin{equation} \omega_\ell^k = \left\{\begin{array}{cc} 1, & e_k = e_\ell \\ \epsilon, & e_k \neq e_\ell \end{array} \right., \end{equation} \tag{ 27 }$$

where $0 \lt \epsilon \lt 1$ is the coupling constant. The smaller the value of ε, the more neighboring pixels that lie on an edge are allowed to differ from their non-edge neighbors. This construction results in a spatially piecewise constant solution; however, unlike our previous methods which attempted to directly estimate MTR boundaries from ECT data, this approach assumes that we have prior knowledge about the boundaries from the SAM data.

Define the $\ell th$ row of the matrix $R_\epsilon \in \mathbb{R}^{N^2 \times N^2}$ by

$$\begin{equation} R_\epsilon\left(\ell,:\right) = \displaystyle \frac{1}{\sum_{k = 1}^{N^2} \omega_\ell^k} \omega_\ell. \end{equation} \tag{ 28 }$$

Then the full regularization operator $L_\epsilon \in \mathbb{R}^{N^2 \times N^2}$ for the structural prior is given by

$$\begin{equation} L_\epsilon = \mathbf{I} - R_\epsilon, \end{equation} \tag{ 29 }$$

where I is the $N^2 \times N^2$ identity matrix. Using this regularization method, the value of $s^*$ is given by

$$\begin{equation} s^* = \arg \min \left\| A_1 s + b_1 - \left(A_2 y + b_2\right) \right\|_2^2 + \rho^2 \left\|L_\epsilon s \right\|_2^2, \end{equation} \tag{ 30 }$$

where ρ is the regularization parameter. This implies that $s^*$ can be found by solving the linear system

$$\begin{equation} \left[ \begin{array}{c} A_1 \\ \rho L_\epsilon \end{array} \right]s = \left[ \begin{array}{c} A_2y + b_2 - b_1 \\ 0 \end{array} \right] \end{equation} \tag{ 31 }$$

in the least squares sense. Note that unlike the method described in the previous subsection, this technique does not require an iterative process to find $s^*$.

We used a standard canny edge detection algorithm built into MATLAB to detect grain boundaries from the SAM data shown in figure 3(b). Figure 5 shows the detected edges overlaid on top of maps of the first and second Euler angles. Note that although SAM cannot provide unique orientation information, it is sensitive to changes in both the first and second Euler angles. The grain boundaries as determined by the SAM data match up very well with the true grain boundaries from the specimen. These edges are used to build our structural prior when we perform inversion on the ECT data.

Zoom In Zoom Out Reset image size

We generated a set of $n = 18\,000$ simulated microstructures and simulated ECT data sets to train the CGMCA maps. Each microstructure is 1 mm × 1 mm in size and consists of two regions, with each region defined by a single crystallographic orientation. A sample microstructure is shown in figure 6. Each pixel in the simulated microstructure represents an area of $25 \; \mu \text{m} \times 25 \; \mu\text{m}$. For each microstructure, we simulated ECT data using an elliptical coil with axes equal to 200 µm and 635 µm. The assumed liftoff of the probe was 90 µm. The simulated data was recorded with a step size of 100 µm. The ECT data collection was simulated for two rotations of the elliptical coil (0^∘ and 90^∘); the imaginary portion of simulated ECT data for a sample simulated microstructure for both rotations is also shown in figure 6.

Zoom In Zoom Out Reset image size

Our first test was to segment the data along the second Euler angle. Thus, the regions in each training segmentation were labeled based on the value of the second Euler angle for that region. Note that each training segmentation was downsampled from 25 µm step size to 50 µm in order to match the resolution of the SAM data.

We trained the CGMCA maps on the segmentations matched to the simulated ECT data sets. Let Z₁ and Z₂ denote the matrices that contain the training data for the segmentation and ECT data, respectively. The dimensions of a single training segmentation are N × N, where N = 21. Each column of $Z_1 \in \mathbb{R}^{N^2 \times n}$ contains one of the training segmentations stacked into a single vector. A single set of ECT measurements is of the size M × M, where M = 11. Each column of $Z_2 \in \mathbb{R}^{4M^2 \times n}$ contains the real and imaginary portion of the ECT signal from each rotation for the corresponding segmentation in Z₁ stacked into a single column vector.

We chose the matrices Γ₁ and Γ₂ for each CGMCA map in the same manner as in [33]. Specifically, we assume that $\Gamma_1 = \Gamma_2$ and assign a single covariance matrix $\Gamma \in \mathbb{R}^{k \times k}$ to the common domain. Let $\mu_1 \in \mathbb{R}^{N^2}$ be the mean vector of the columns of Z₁. Denote the matrix containing the centered segmentation data sets as S₁, where S₁ is given by

$$\begin{equation} S_1 = Z_1 - \mu_1. \end{equation} \tag{ 32 }$$

Then denote the SVD of S₁ by

$$\begin{equation} S_1 = U_{S_1} \Sigma_{S_1} V_{S_1}^T, \end{equation} \tag{ 33 }$$

and denote the projection of S₁ onto the first k columns of S₁ as $P_1 \in \mathbb{R}^{k \times n}$. The covariance matrix Γ is given by

$$\begin{equation} \Gamma = \frac{1}{n} \sum_{\ell = 1}^n \left(p_\ell - \overline{p}\right)\left(p_\ell - \overline{p}\right)^T, \end{equation} \tag{ 34 }$$

where $p_\ell$ is the $\ell$th column of P₁ and $\overline{p}$ is the mean vector of the columns of P₁. This choice of Γ was shown to work well for in the segmentation algorithm in [33].

The dimension of the common domain k is chosen based on the singular values of S₁, in a similar manner as in [33]. Specifically, k is chosen as the index where the singular values of S₁ begin to level off. Here, we let k = 50.

Note that our training ECT data sets have dimensions 1 mm × 1 mm. However, the ECT data of the large grain Ti specimen represents an area of 33.4 mm × 18 mm. Thus, we need a way to process a larger data set through the segmentation algorithm. We achieve this by processing the data in 1 mm × 1 mm squares. The process is outlined below.

• (i)  
  Extract a 1 mm × 1 mm square of data from the large ECT data set to form y.
• (ii)  
  Map y to the common domain, i.e. compute $A_2y + b_2$.
• (iii)  
  Find $s^*$ that satisfies

  $$\begin{equation} s^* = \arg \min \left\| A_1s + b_1 - \left(A_2y + b_2\right) \right\|_2^2 \end{equation} \tag{ 35 }$$

  while also satisfying the constraints imposed by either the image deblurring prior or the structural prior.
• (iv)  
  Repeat for another extracted data set that partially overlaps the previous one.
• (v)  
  Repeat steps (i)–(iv) until the entire data set is processed.

Once each of the ECT data subsets have been processed, we need to combine them into a single final segmentation, which we denote by $S$. Let $x_\ell$ denote the $\ell$th spatial location in S. After completing the process outlined above, suppose there are a total of $m_\ell$ smaller segmentations that cover $x_\ell$. Denote the kth segmentation that contains $x_\ell$ by $s^{\ell k}$. Then the value of the large segmentation S at the point $x_\ell$ is given by

$$\begin{equation} S\left(x_\ell\right) = \displaystyle{\frac{1}{m_\ell}}\sum_{k = 1}^{m_\ell}s^{\ell k}\left(x_\ell\right). \end{equation} \tag{ 36 }$$

Processing the large data set in smaller portions is beneficial because each portion can be run independently in parallel. Additionally, smaller data sets are computationally faster to run through the segmentation algorithm. The next section will discuss the application of the segmentation algorithm to the large grain ECT data.

As previously mentioned, We simulated ECT data sets of the large grain specimen for two rotations of the ECT coil (specifically 0^∘ and 90^∘). The imaginary portion of the noiseless simulated data sets is shown in figures 3(c) and (d). To both the real and imaginary portions of this data, we added zero-mean Gaussian noise with standard deviation equal to 5% of the maximum value of the signal. Note that although we only show the imaginary portions, both the real and imaginary portions of the ECT data were used to estimate the segmentation.

We applied the process outlined in section 4.3 to the noisy ECT data to find the segmentation of the specimen along changes in the second Euler angle. We estimated the segmentation assuming both the SAM structural prior and the image deblurring prior. For the structural prior, we used a value of ε = 0.005 for the coupling constant and ρ = 5 for the regularization parameter. For the image deblurring prior, we chose $\gamma = 2e-3$. The results are shown in figure 7. The image deblurring prior does well at estimating the orientation of the individual regions. However, it struggles to resolve some of the boundaries and smaller features. In contrast, the segmentation using the structural prior has sharp boundaries, and is able to resolve the small features that are missed with the image deblurring prior.

Zoom In Zoom Out Reset image size

To further illustrate this point, we also show a zoomed-in image of a subset of each result in figures 7(d)–(f). These images highlight the inability of ECT to resolve certain features on its own, such as the small teardrop region in the bottom middle of the image. In addition, it is obvious that the estimated boundary of the larger region is much clearer with the addition of the structural prior.

Another consideration in comparing these two methods is the computational time required to compute each segmentation. Both estimates were computed on a Intel Core-i9 with a 2.4 GHz CPU with 8 cores. The result using the image deblurring prior took around 30 min to compute, while the structural prior took approximately 12 min. This is due to the fact that the image deblurring prior requires an iterative process to find the unknown segmentation, while the structural prior only needs to solve a single linear system. Thus, the addition of structural information greatly decreases the computational time.

We also consider the effect of missing information in the structural prior. We removed some of the edges shown in figure 5 and reran the ECT inversion using the structural prior. The results are shown in figure 8. Note that we only show the result for the zoomed-in area where edges were removed, as the rest of the solution remains the same. For comparison, we also show the result with all edge information, as shown in figure 7.

Zoom In Zoom Out Reset image size

Without edge information, the algorithm reverts to a smoothness prior. While it will still detect an orientation change, boundary information is completely lost. Furthermore, the teardrop shaped feature is no longer visible. Ideally, we would like for the algorithm to revert to the image deblurring prior if no structural information is available. Although the boundary of the small feature is not well defined using the image deblurring prior (as seen in figure 7(d), the feature is still visible. Future work will modify the structural prior so that it will revert to the image deblurring prior if necessary, similar to what was done in [19], which utilized total variation regularization if structural information was not available.

As previously mentioned, we use ECT data from multiple rotations of the ECT coil in order to produce the images in figure 7. However, for illustrative purposes, we also estimated the segmentation using ECT data from a single rotation of the coil. In particular, we computed the segmentation using data from only the 0^∘ rotation and only the 90^∘ rotation using both the image deblurring prior and the structural prior. These results are shown in figure 9. When using the image deblurring prior, neither rotation on its own is sufficient to recover information about the grain edges. However, even with the edge information provided by the structural prior, multiple rotations of the coil are still needed to provide correct orientation information. This is due to the fact that two different values of θ can produce the same ECT response for a single rotation of the coil depending on the value of ψ₁. Adding a second rotation of the coil resolves this uncertainty, as is evident in figure 9.

Zoom In Zoom Out Reset image size

As shown in [33], we can also use this method to segment the ECT data along changes in the first Euler angle, with a caveat. Recall that for the HCP crystal, $\psi_1 \in [0, 2\pi]$. However, due to the symmetry of the crystal, it is more difficult to develop a labeling scheme segmentation along the first Euler angle. In particular, certain crystallographic orientations can be represented by multiple different values of ψ₁. For instance, if $\theta = \pi/2$, then ψ₁ = 0 or $\psi_1 = \pi$ describe the same orientation. This ambiguity will cause the algorithm to fail. To avoid this non-uniqueness issue, we enforce $\psi_1 \in [0, \pi/2]$ when running the segmentation algorithm. Developing a labeling system for ψ₁ is the topic of future work.

We relabeled the training segmentations based on the value of the first Euler angle and reran the process described in section 4.3 assuming both the image deblurring prior and the structural prior. The results are shown in figure 10. Similar to the results with the second Euler angle, the estimates of the orientation using either prior are very similar. However, as before, the image deblurring prior is unable to resolve many of the boundaries between grains, while the structural prior provides much sharper boundaries. We note that ECT is generally less sensitive to the value of ψ₁, particularly if the value of θ is close to zero. As such, the estimates of ψ₁ using either prior are much less accurate for values of θ close to zero.

Zoom In Zoom Out Reset image size

We can combine the segmentation results for each Euler angle into a single image. Figure 11 shows the combined results for the image deblurring method and the structural prior. Each point in the result is colored based on the estimated value of the first and second Euler angle at that point. Also shown in figure 11 is the PLM data of the specimen. The result using structural prior matches the PLM data very well. The advantage of the structural prior is twofold: it improves the resolution of the final image, and it speeds up computation time. The structural prior shows promise to take advantage of the benefits of each imaging modality.

Zoom In Zoom Out Reset image size