A spherical harmonic transform approach to the indexing of electron back-scattered diffraction patterns

Ultramicroscopy 207 (2019) 112841

Contents lists available at ScienceDirect

Ultramicroscopy journal homepage: www.elsevier.com/locate/ultramic

A spherical harmonic transform approach to the indexing of electron backscattered diffraction patterns W.C. Lenthea, S. Singhb, M. De Graef a b

T

⁎,a

Department of Materials Science and Engineering, Carnegie Mellon University, 5000 Forbes Avenue, Pittsburgh, PA 15213, USA Lawrence Livermore National Laboratory, 7000 East Ave, Livermore, CA 94550, USA

A R T I C LE I N FO

A B S T R A C T

Keywords: EBSD Indexing Spherical harmonic transform

A new approach is proposed for the indexing of electron back-scattered diffraction (EBSD) patterns. The algorithm employs a spherical master EBSD pattern and computes its cross-correlation with a back-projected experimental pattern using the spherical harmonic transform (SHT). This approach is significantly faster than the recent dictionary indexing algorithm, but shares the latter’s robustness against noise. The underlying theory is presented, followed by example applications, one on a series of Ni EBSD data sets recorded with decreasing signal-to-noise ratio, the other on a large shot-peened Al data set. The dependence of indexing speed and memory usage on the SHT bandwidth is explored. The speed gains of the new algorithm are achieved by executing real-valued Fast Fourier Transforms, explicitly incorporating crystallographic symmetry in the crosscorrelation computation, and using efficient loop ordering to improve the caching behavior. The algorithm produces a cross-correlation array in the zyz Euler space; an orientation refinement procedure is proposed based on analytical derivatives of the Wigner d functions. The new approach can be applied to any diffraction modality for which the scattered intensity can be represented on a spherical surface.

1. Introduction Electron back-scattered diffraction (EBSD) patterns consist of a set of Kikuchi bands superimposed on a continuous background of backscattered electrons. These patterns have traditionally been indexed by extracting the positions and orientations of the bands by means of a Hough transform [1–3]. If the detector geometry has been calibrated, the normals to the lattice planes that gave rise to the bands can be determined, and a comparison of the angles between pairs of normals (or zone axes) and a set of pre-computed angles based on the known crystal structure then leads to the orientation of the lattice at the illuminated point on the sample. This analysis can be performed efficiently, and modern detector systems can be operated at a rate of several thousand patterns per second. When the pattern quality is poor, however, or when overlapping patterns occur, the indexing success rate of the Hough-based approach typically decreases [4], and alternative indexing methods must be used. One such approach is the dictionary indexing (DI) method [5], which relies on a comparison, using a similarity metric, of pre-computed EBSD patterns for a uniform sampling of orientation space with experimental patterns; the similarity metric is typically the dot product between the normalized EBSD patterns. This approach has been shown

⁎

to be robust against noise [4] and can handle overlapping EBSD patterns that arise near grain or phase boundaries. While the DI approach has been successfully applied to a number of material systems [6–9], the main drawback of this technique is the fact that the computational requirements increase with decreasing crystal symmetry. This is due to the fact that the size of the pattern dictionary is proportional to the volume of the fundamental zone (FZ) of orientations; for cubic symmetry, this zone has a volume that is 1/24th of the total volume of orientation space (π2 in the homochoric orientation representation [10]), but in lower symmetry systems this volume, and thus the number of dictionary patterns to be computed, increases significantly. For instance, for an average step size of 1.4∘ in orientation space, the number of dictionary patterns is 333,227 in the cubic 432 crystal point group, and 2,000,037 for orthorhombic 222 symmetry [11]. Hence, at the same angular step size, an indexing run for an orthorhombic material will take approximately six times as long as for a cubic material. Other researchers have attempted to mitigate the long execution times by preprocessing of patterns into clusters [12]. Nolze et al. [13] applied a similar pattern matching technique to pseudo-symmetry issues, and Tanaka and Wilkinson [14] analyzed the accuracy and precision of the use of pattern matching for pattern center, orientation, and strain determination.

Corresponding author. E-mail addresses: [email protected] (W.C. Lenthe), [email protected] (S. Singh), [email protected] (M.D. Graef).

https://doi.org/10.1016/j.ultramic.2019.112841 Received 2 May 2019; Received in revised form 17 July 2019; Accepted 2 September 2019 Available online 03 September 2019 0304-3991/ © 2019 Elsevier B.V. All rights reserved.

Ultramicroscopy 207 (2019) 112841

W.C. Lenthe, et al.

distributions of the back-scattered electrons (BSEs) for a given crystal structure, sample tilt, and microscope accelerating voltage; (2) a dynamical scattering simulation using the Bloch wave formalism combined with the Monte Carlo information; and (3) a geometrical model for the detector-sample geometry that allows for accurate interpolation of detector pixel intensities from the K-sphere. For details, we refer the interested reader to Callahan and De Graef, [15]. The anisotropic background is typically removed from the simulated and experimental patterns prior to their use in indexing computations by means of a combination of high-pass Fourier filtering and adaptive histogram equalization [23]. The Cartesian reference frame at the center of the Ksphere has its x axis aligned along the a Bravais basis vector, and z along the reciprocal c* direction, i.e., normal to the a–b plane.

The DI approach is based on a forward model for EBSD patterns that consists of a Monte Carlo simulation to obtain an estimate of the directional, energy, and depth distributions of the back-scattered electrons (BSE), and a dynamical scattering computation of the BSE yield on the Kikuchi sphere [15]. This spherical intensity data (see also [16]), which we refer to as the Master Pattern (MP), is then projected onto a square Lambert grid using an equal-area projection [17]; individual EBSD patterns can be generated using a bilinear interpolation on this square grid. Since the MP consists of data points on a sphere, an alternative approach would involve projecting the experimental patterns onto the sphere, and determining the 3D rotation that maximizes the correlation between the master pattern and the projected experimental pattern. Such a spherical correlation is widely used in other fields, e.g., in the astronomical community [18], for image processing [19,20], and shape registration [21], and its use for fast EBSD pattern indexing forms the subject of the present contribution. The use of spherical harmonic analysis (SHA) for pattern indexing was first proposed in [22]. While our approach is based on the same SHA theory, there are several differences between our algorithm and that used in [22]. The structure of the present paper is as follows: first we provide a brief description of the underlying spherical harmonic correlation theory. Then, in Section 2, we apply this approach to the EBSD indexing case, and we introduce a modified form of the EBSD master pattern that is better suited to the Spherical Harmonic Transform (SHT) in Section 2.6. We analyze the appropriate bandwidth to be used for the correlation analysis, as well as the peak finding and refinement steps. Finally, in Section 3, we provide a number of example indexing results and a comparison of the indexing speed of the new algorithm with that of the dictionary indexing approach. We conclude this paper with a summary and a discussion of the availability of the algorithm source code.

2.2. Spherical harmonic transforms The continuous Spherical Harmonic Transform (SHT) of a function f on the sphere 2 is defined in terms of the spherical angles (θ, ϕ) by: ∞

f (θ , ϕ) =

+ℓ

ℓ f^m Ymℓ (θ , ϕ),

∑ ∑

(1)

ℓ = 0 m =−ℓ

where ℓ f^m =

π

2π

∫0 dθ sin θ ∫0 dϕ f (θ, ϕ) Ymℓ (θ, ϕ),

(2)

where the over-bar indicates the complex conjugate and the spherical harmonics Ymℓ (θ , ϕ) have their usual normalized definition:

Ymℓ (θ , ϕ) =

(2ℓ + 1) ! (ℓ − 1)! ℓ Pm (cos(θ))eimϕ , 4(ℓ + m)!

(3)

the associate Legendre polynomial of degree ℓ and order m. with The coefficients for the discrete SHT can be defined as follows:

Pmℓ (x )

2. Theory

ℓ f^m ≈

∑ fk wk Ymℓ (θk , ϕk ),

(4)

k

2.1. EBSD master pattern

with fk the function value at grid point k, θk, ϕk the spherical coordinates, and wk the weighted solid angle; this expression is exact if the correct weight factors are used. If all the grid points on the sphere can be organized into constant latitude rings with equal angular spacing along each ring, then the discrete SHT can be computed efficiently via Fast Fourier Transforms (FFT) [24]. The discrete SHT is now given by:

The spatial distribution of back-scattered electrons can be represented on the Kikuchi sphere (K-sphere), as illustrated in Fig. 1 for Nickel at 20 kV. The plane tangent to the sphere represents the EBSD detector. The pattern on the detector is a gnomonic projection of the intensity distribution on the K-sphere; as the crystal orientation changes, the sphere rotates about its center while the detector remains stationary. The point of tangency between the sphere and the detector plane represents the pattern center, and the detector-sample distance uniformly scales the size of the detector pixels. The master pattern is computed using a physics-based forward model consisting of: (1) a Monte Carlo simulation that produces the spatial, energy, and depth

ℓ f^m =

Ny

∑ Gm,y Ymℓ (θy ),

(5)

y

where Ny

Gm, y ≡ wy exp(−imϕ0, y )

2π imx ⎞ ; Nx ⎠ ⎝

∑ fx,y exp ⎛− ⎜

k

⎟

(6)

In this expression, the data is organized into Ny rings with latitudes θy with the first point in ring y at longitude ϕ0,y. All Gm,y for a single y (i.e., for all m values) can be computed simultaneously via a single 1-D FFT. If a single weight factor is used for each ring (as opposed to each point in a least squares approach) and the rings are symmetric about the equator, then the ring weights can be computed as the solution to the following linear equation [25]:

1 ⎛ cos(2 θ1) ⎜ ⎜ cos(4θ1) ⎜ ⋮ ⎜ cos(2N θ ) y 1 ⎝

Fig. 1. The Kikuchi sphere is shown (left) for Nickel at 20 kV. On the right, a tangent plane representing the detector is drawn along with the gnomonic projection from the center of the sphere onto the plane; the point of tangency is referred to as the pattern center.

⎛ −11 ⎞ 1 ⎞ w1 − 1 ⎟ ⎛ w2 ⎞ ⎜ 3 ⎟ ⎜ −1 ⎟ 1 ⎟⎜ ⋮ ⎟ = ⋯ cos(4θ Ny ) ⎜ 15 ⎟ w N y ⋱ ⋮ ⋮ ⎟⎜ ⎟ ⎜ ⋮ ⎟ ⋯ cos(2Ny θ Ny ) cos( ±1) ⎟ ⎝ w Ny + 1⎠ ⎜ −2 1 ⎟ ⎠ ⎝ 4Ny − 1 ⎠ ⋯ ⋯

1 cos(2θ Ny )

(7)

The maximum bandwidth (highest ℓ) that can be computed is Ny/2. In 2

Ultramicroscopy 207 (2019) 112841

W.C. Lenthe, et al.

2.4. Cross correlation

the special case that the ring latitudes are zeros of the Legendre polynomial of degree Ny, the maximum bandwidth that can be computed is Ny . If the ring latitudes are symmetrically distributed about the equator, one can leverage the fact that the odd harmonics are anti-symmetric across the equator and the even harmonics symmetric to reduce the complexity of Eq. (5) by a factor of 2. Finally, if the function is real valued then the spherical harmonic coefficients are conjugate symmetric: ℓ ℓ f^−m = (−1)mf^m .

It is well known that the spherical harmonic functions can be rotated on the sphere by means of the Wigner D function; this operation is analogous to carrying out a translation by means of a phase shift in the Fourier domain. Consider a 3-D rotation R(α, β, γ), decomposed into three Euler angles (α, β, γ) by consecutive rotations around the zyz axes. Denoting the spherical coordinates by the symbol ω, the rotated spherical harmonic function can be written as

Ymℓ (ω) =

∑ Ynℓ (R−1ω) Dmℓ ,n (R),

(8)

n

where the symbol defined as:

Any set of grid points distributed into Ny rings with equally spaced latitudes and symmetry across the equator is suitable for implementing a fast discrete SHT, but rings coincident with roots of the Legendre polynomial are often preferred to achieve a higher bandwidth with the same number of points [25], making a Gauss-Legendre grid a common choice [24,26,27]. Many SHT implementations include a polar optimization since the magnitude of higher order harmonics decreases exponentially near the poles [28]. To avoid over-sampling, our discrete SHT implements two grid choices with fewer points near the poles: a square Lambert projection providing equal area [29] and a new ‘square Legendre’ grid which modifies the square Lambert projection by moving ring latitudes to Legendre polynomial roots. These grids implicitly implement polar optimization by effectively truncating the FFT of rings near the poles and are shown in Fig. 2. It may be advantageous to center the back projected pattern around regions of higher grid point density, but spatial resolution is also limited by the harmonics themselves. In general the square Legendre grid is preferred in this work to balance uniform grid density with bandwidth per point, but any grid can be used with the discrete SHT which accounts for less than 5% of the overall calculation time for reasonable bandwidths.

Dmℓ , n (R)

(9)

represents the Wigner D function which is

Dkℓ, m (α, β, γ ) = dkℓ, m (β )eimα eikγ ,

(10)

with the Wigner d function defined as:

β k+m β k−m (ℓ + k ) ! (ℓ − k )! cos ⎛ ⎞ sin ⎛ ⎞ Pℓk−−km, k + m (cos(β )); (ℓ + m) ! (ℓ − m)! ⎝2⎠ ⎝2⎠

dkℓ, m (β ) =

(11) in this expression, Pℓk, m (x ) is a Jacobi polynomial. Recursion relations can be used to efficiently compute the Wigner d function for all ℓ, k, and m at a single β [30,31]. The spherical cross correlation (SCC) between two spherical functions f and g can be computed using the Wigner D function

(f ★g )(R) =

∑

ℓ ℓ f^m g^n Dmℓ , n (R).

ℓ, m, n

(12)

Leveraging the definition of the Wigner D function and decomposing R into two rotations allows for the cross correlation to be computed efficiently via a 3D FFT [21]:

(f ★g )(α, β , γ ) = F−1

2.3. Crystallographic symmetry and the discrete SHT As described in Section 2.1, the EBSD master pattern inherits the symmetry of the underlying crystal structure projected onto the Ksphere surface. Due to the underlying crystallographic symmetry, the SHT of the K-sphere will inherit all the symmetry properties from the crystal lattice. In particular, this means that f mℓ = 0 when one or more of the following relations is valid, depending on the crystal point group:

ℓ ⎧ ∑ f^ g^ ℓ d ℓ ⎛ π ⎞ d ℓ ⎛ π ⎞ ⎫⎛α + π , β + π , γ + π ⎞. ⎨ ℓ m n m, k ⎝ 2 ⎠ k, n ⎝ 2 ⎠ ⎬⎝ 2 2⎠ ⎭ ⎩

(13)

The normalized spherical cross correlation over a mask function M can be computed using 3 un-normalized cross correlations and the integral of the mask function [32]:

(f ★g )(R) =

((f − fw ) ★ (g − g¯))(R) (f 2 ★M )(R) − 2fw (f ★M )(R) + (fw 2 ∫ M )(R)

(14)

where g¯ is the mean of g over the window and fw is the mean of f over the rotated window:

• Mirror plane in the x–y plane: mod (ℓ + m, 2) ≠ 0; • Inversion symmetry: mod (ℓ, 2) ≠ 0; • N-fold rotational symmetry about the z axis: mod (m, N ) ≠ 0.

fw (R) =

Other common crystallographic symmetry operators, e.g. rotational symmetry about other axes, do not create an easily leveraged structure in the spherical harmonic coefficients because they do not occur regularly in the individual spherical harmonics (basis functions). Fig. 3 shows schematic plots of the (ℓ, m) arrays for different point group symmetries; white pixels indicate zero SHT coefficients, black pixels indicate potentially non-zero coefficients. The plots are important for the cross-correlation computation described in the next section, because they allow the summations and Fourier transforms to be carried out more efficiently. A mirror plane at the equator reduces the computation of the cross correlation sums by a factor of 2. An N-fold rotational symmetry axis at the north pole reduces both the summation and the Fourier transform times by a factor of N. Note that this explicit inclusion of crystallographic symmetry relations allows for an efficient compression of the master pattern data into a relatively small set of spherical harmonic coefficients.

(f ★M )(R) . ∫M

(15)

As stated in the previous section, a large number of redundant multiplications can be avoided by leveraging the systematic zeros and the conjugate symmetry of real valued functions, leading to significant speed-up factors compared to the full computation. The systematic zeros from an equatorial mirror plane and/or inversion symmetry reduce the time to compute the value of a single point sum in Eq. (13). Systematic zeros from N-fold rotational symmetry about the z axis in f make entire m planes zero and manifest as periodicity in the real space cross correlation grid with (f ★g )(α, β , γ ) = (f ★g )(α + 2π / N , β , γ ) . This periodicity effectively reduces the 3D FFT size from (2ℓmax − 1)3 to (2ℓmax − 1) (2ℓmax − 1)2 . Rotational symmetry in g produces an equivalent N effect for n planes and γ. Furthermore, the glide plane symmetry of Euler space described in Section 2 allows for an additional speed-up, since there is no extra information in the cross-correlation volume for β > π. The conversion from the (α, β, γ) Euler angles to the more conventional (φ1, Φ, φ2) Bunge Euler angle triplet used in materials science is readily carried out by means of the following relation: 3

Ultramicroscopy 207 (2019) 112841

W.C. Lenthe, et al.

Fig. 2. The northern hemisphere of a square Lambert or square Legendre grid is shown with rings of constant latitude indicated by dashed lines (top left). The grid points for each constant latitude ring are shown schematically below with an FFT computed across each row during a SHT. Gauss-Legendre type grids have an extreme disparity between grid point size near the poles and equator (top right). The square Lambert grid has constant pixel size (blue line) while the square Legendre grid balances Legendre root latitudes with moderately uniform grid point sizes. Legend names indicate the number of constant latitude rings (equal to the side length for odd sized square grids). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

(φ1, Φ, φ2) = (α −

π π , β , γ + ). 2 2

following equation for the angular updates (Δα, Δβ, Δγ):

(16)

⎛ X , αα X , αβ X , αγ ⎞ ⎛ Δα ⎞ ⎛ X , α ⎞ ⎜ X , βα X , ββ X , βγ ⎟ ⎜ Δβ ⎟ = ⎜ X , β ⎟ ⎜ ⎟ ⎜X X , γβ X , γγ ⎟ ⎝ Δγ ⎠ ⎜⎝ X , γ ⎟⎠ ⎝ , γα ⎠

These Euler angles describe the rotation from f to g. If f and g correspond to the master and experimental patterns respectively, the rotation must be inverted to describe orientations in the crystal to sample frame convention.

(17)

where X is the cross correlation from Eq. (12), and the notation X,i indicates a partial derivative with respect to the variable i. Newton’s method achieves quadratic convergence, providing a significant advantage over derivative free techniques. In testing our implementation against perfectly rotated harmonics, convergence required 3 or fewer iterations with a typical error of ≈ 10−6 ∘. These expressions require summation over partial derivatives of the Wigner D function; the Wigner d function is numerically challenging to compute with recursive methods generally employed. Efficient expressions for the first and second derivatives of the Wigner d function have been derived and are listed in Appendix 1.

2.5. Orientation refinement The FFT method of computing the spherical cross correlation (Eq. (13)) yields the cross correlation values on a discrete grid of zyz Euler angles. However, since an analytical expression for the cross correlation as a function of rotation angles exists (Eq. (12)), the maximum cross correlation rotation can be refined efficiently via Newton’s method if an estimate sufficiently close to the peak is known. To do so, requires the computation of the Hessian and the Jacobian to solve the 4

Ultramicroscopy 207 (2019) 112841

W.C. Lenthe, et al.

Fig. 3. Point group symmetry (left) introduces systematic zeros in the spherical harmonic transforms of master patterns (right). Potential nonzero spherical harmonic 0 transform coefficients are indicated by black pixels and systemic zeros in white. The DC value f^ is in the top left corner of the images with ℓ increasing down the 0

image and m increasing to the right. Only non-negative m values are considered due to conjugate symmetry.

combined energy distribution for all pixels on the detector. Although this approach fails to capture variations in energy distributions across the detector, they are minimal in practice [33]. The experimental EBSD pattern is projected onto the sphere by inverting the expressions from [15] that project from the sphere to the detector plane. The intensities are rescaled to have a zero mean intensity, and the remainder of the sphere is then padded with zeroes. The reason for putting the mean intensity at zero is that the cross correlation of an arbitrary pattern will be higher near the brighter regions of the master pattern (i.e., mostly the zone axis locations); since the pattern window integrates over a higher valued region of the master pattern, the cross correlation would always be higher in the brighter regions of the master pattern. Zeroing the mean pixel value inside the window and rescaling so the standard deviation equals unity approximates the normalized cross correlation, resolving the brightness issue. The true normalized cross correlation additionally requires division by a denominator that is constant for a given master pattern and detector geometry (mask function). The SHT of the back projected experimental pattern is computed and used in the cross-correlation expression (Eq. (13)) to obtain an

2.6. Application to EBSD pattern indexing The point of tangency between the sphere and the detector plane represents the pattern center, and the detector-sample distance uniformly scales the size of the detector pixels. As discussed before, the master pattern is computed using a combination of a Monte Carlo simulation and a dynamical scattering simulation [15]. The anisotropic background is typically removed from the simulated and experimental patterns prior to their use in indexing computations by means of a combination of high-pass Fourier filtering and adaptive histogram equalization [23]. The Cartesian reference frame at the center of the Ksphere has its x axis aligned along the a Bravais basis vector, and z along the reciprocal c* direction, i.e., normal to the a–b plane. The physics-based forward model produces a distinct master pattern for each energy in a range of energy bins and an energy distribution for each detector pixel. To simulate a single EBSD pattern for dictionary indexing the energy distribution of each pixel is used to compute a weighted average of the master patterns for each energy. Spherical indexing requires cross correlation against a single function. A combined energy weighted master pattern is computed using a single

5

Ultramicroscopy 207 (2019) 112841

W.C. Lenthe, et al.

Fig. 4. A flow chart is shown for the indexing algorithm as implemented. The theoretical scaling for memory and computational complexity are shown as a function of bandwidth b in the bottom left and right corners of each step.

array of intensities in the (α, β, γ) Euler space. Finding the maximum intensity and refining its location using the approach described in the previous section then results in a final orientation for the experimental pattern. In multi-phase systems, the same approach can be used by repeating the spherical cross correlation computation for each master pattern; the highest peak in the set of auto-correlation volumes for all phases then identifies both the phase and the orientation. A schematic of the EBSD indexing algorithm is shown in Fig. 4.

3. Applications 3.1. Indexing examples In this section we apply the spherical indexing approach to two different experimental data sets. The first data set consists of 10 EBSD scans of the same region in a Ni sample using different detector gain and exposure time settings, keeping the total exposure constant [4, section 2.2]. The region of interest has 186 × 151 sampling points in a square grid with a step size of 1.5 µm; individual 8-bit patterns were binned to 60 × 60 pixels, and the pattern center coordinates are (x *, y *, z *) = (0.50726, 0.73792, 0.55849) in the EDAX/TSL convention.

6

Ultramicroscopy 207 (2019) 112841

W.C. Lenthe, et al.

Fig. 5. Spherical harmonic indexing for two bandwidths (53 and 63) is compared to traditional Hough transform indexing (HI) and dictionary indexing (DI) for multiple scans of the same sample with a range of camera gains and exposure times, resulting in approximately constant total exposure (see also [4]).

The sample tilt was 75.7∘. The Hough parameters were optimized to give the best possible indexing success rate for the optimal camera gain setting [34]. The second data set was acquired from a shot-peened Aluminum 7075-T651 alloy at the base of a shot peen crater. An area of 56 µm × 52 µm was scanned at 10 kV using a scan step size of 25 nm, resulting in 1931 × 1784 (about 3.4 million) patterns of size 160 × 120 pixels; for details about the sample preparation and pattern acquisition conditions we refer to Singh et al. [8]. It should be noted that the indexing success rate of all three indexing approaches (Houghbased, DI and SI) depends on accurate knowledge of experimental parameters, in particular the pattern center coordinates.

based indexing comparing favorably to the traditional Hough indexing in all cases for relatively low bandwidths. The left-most column shows a representative EBSD pattern at each detector gain setting, with a decreasing signal-to-noise ratio from top to bottom. Hough indexing (HI) and dictionary indexing (DI) results are shown in the second and third columns, respectively, as [001] inverse pole figure maps. The final two columns show the spherical indexing results for two different bandwidths. To achieve the robustness of dictionary indexing at the extreme noise level in the bottom row of Fig. 5 requires significantly higher bandwidths, and thus a longer computational time, as is illustrated in Fig. 6. Note that a bandwidth of 158 produces an IPF map that is comparable to the DI map. Pre-processing of the experimental patterns using a non-local means approach [35] has been reported to further improve the DI indexing results for the lowest signal-to-noise ratio, and we anticipate a similar improvement for the new spherical indexing algorithm. To validate the absolute accuracy of indexing, simulated EBSD patterns were indexed for a range of noise levels using the same

3.1.1. Ni data set results The set of 10 scans for different detector gains represents an opportunity to evaluate the robustness of the spherical indexing approach with respect to pattern noise. Indexing results for a range of bandwidths are shown in Fig. 5 with the noise tolerance of spherical harmonic 7

Ultramicroscopy 207 (2019) 112841

W.C. Lenthe, et al.

Fig. 6. At the extreme noise levels increased bandwidths are required to maintain the robustness of dictionary indexing (DI).

corresponding to 333,227 patterns [11]. This time does not include the subsequent refinement of the orientations, which took another 30 h. For the shot-peened Al data set, the complete SI run, including orientation refinement, for a bandwidth of 41 took 39 min, or 1466 patterns per second (i.e., about 125 times faster than the DI indexing run), making the approach comparable to real-time indexing by means of the standard Hough-based algorithm. The rates listed in Table 3 are for the Ni dataset but speeds for the Al dataset are comparable (within 3%) since bandwidth dominates computation time making indexing relatively insensitive to pattern size. Increasing the speed of the original dictionary indexing algorithm to the point where real-time indexing becomes feasible can in principle be accomplished by extending the algorithm to address multiple GPUs and a large number of compute threads; this would likely be prohibitively expensive in terms of compute hardware, and only possible in locations that have high performance computing facilities near the scanning electron microscope. Spherical indexing, on the other hand, requires far less compute power; in its current implementation, either in C++ or fortran-90, the code employs standard multi-threading facilities to run on multiple cores. For the Nickel data set described in Section 3.1, an indexing speed of 1450 patterns per second was achieved on a 24-core Linux CentOS 7 system, with Intel Xeon CPU E5-2670 v3 cores running at 2.3 GHz. A list of indexing rates for a range of bandwidths is given in Table 3. The execution times reported in Table 3 include pattern refinement as well as the symmetry-based optimization of the auto-correlation step. An analysis of the dependence of execution speed on bandwidth reveals that the indexing time per pattern is proportional to b3log b3, with b the bandwidth; the memory usage scales as b3. Both complexity and memory scaling are consistent with the (2b − 1)3 FFT in Eq. (12) (the ‘Correlate’ step in Fig. 4) being the dominant computational step. Results for the shot-peened Al data set using 48 compute threads (Intel(R)

detector geometry as the Ni dataset. The average solid angle of detector pixels is 4.06 × 10 −4 sr corresponding to a cone angle of 1.30∘. A dictionary of orientations was created using a coarse Euler angle grid with (φ1, Φ, φ2) = (29∘i, 19∘j, 23∘k ) for i = [0, 12], j = [0, 9], and k = [0, 15]. Noise was added to the patterns for a range of values N by re-scaling the pixel range of all patterns to [0, N], adding Poisson noise to each pixel using the pixel value as the mean, and then dividing by N. The indexing error is shown in Fig. 7 for a range of bandwidths and noise levels demonstrating excellent noise tolerance, accuracy, and precision. The large flat region in the cumulative distribution of the highest noise level (N = 1) indicates systematic misindexing of some patterns, with only 2.3% of errors falling in [2∘, 20∘]. Results are summarized in Table 1 and Table 2 with the mean, median, maximum, and standard deviation of indexing error below the angular resolution of an average pixel for all bandwidths and all but the highest noise level. Table 2 also includes signal-to-noise ratio (SNR) and peak signalto-noise ratio (PSNR) measurements. Correctly indexed points retain high accuracy and precision at the most extreme noise level if systematically misindexed points are excluded with a 10∘ threshold. 3.1.2. Shot-peened Al results Fig. 8 shows a comparison of Hough indexing (HI), dictionary indexing (DI) and spherical indexing (SI) for six different bandwidths. Black points in the HI results indicate un-indexed patterns. The only differences between the DI and SI results are apparent near the base of the shot peen crater (top of the figure); some of the smaller dynamically recrystallized grains are not fully resolved at the lowest bandwidths. The DI and SI results are virtually indistinguishable starting from a bandwidth of 53. The SI algorithm is significantly faster than the DI approach; dictionary indexing of the complete data set using 24 cores and an NVidia GeForce GTX 1080 GPU took a total of 82.5 h or 11.6 patterns per second, using a dictionary with 1.4∘ angular step size, 8

Ultramicroscopy 207 (2019) 112841

W.C. Lenthe, et al.

Fig. 7. Cumulative distributions and kernel density estimates of indexing error against a dictionary of simulated EBSD patterns is shown for a range of bandwidths (top) and a range of noise levels for a bandwidth of 88 (middle). The vertical gray lines are the average pixel size (1.30∘). The kernel density estimate for the highest noise level (1) only includes points below an error of 10∘. Representative simulated patterns are shown (bottom) scaled to the dynamic range of the original image (top row) and the full dynamic range (bottom row).

indexing. The SI algorithm also scales well with the number of compute threads, as shown in Fig. 10(a); the scaling is nearly linear over a range of 1 to 24 threads on a Linux CentOS 7 system with 24 cores. Memory usage is linear in the number of threads, as shown in Fig. 10(b). For the top curve in Fig. 10(a), the FFTW library [36], which is used for all Fourier transforms, was compiled with the Single Instruction Multiple Data (SIMD) option turned on, which results in a speed up by a factor of about 3 × compared to runs without SIMD support. The core of the cross-correlation algorithm involves matrix and vector multiplications, which suggests that the algorithm may benefit, in terms of execution speed, from an implementation on Graphical Processing Units (GPU), provided a suitable FFT library is available; this is a topic of ongoing research.

Table 1 Error distribution measurements are tabulated for indexing simulated patterns at a range of bandwidths demonstrating excellent precision and accuracy. Bandwidth

53

68

88

113

158

Mean Median Std. Dev. Max

0.576∘ 0.637∘ 0.257∘ 1.114∘

0.369∘ 0.388∘ 0.188∘ 0.785∘

0.257∘ 0.250∘ 0.134∘ 0.655∘

0.205∘ 0.201∘ 0.109∘ 0.523∘

0.135∘ 0.133∘ 0.068∘ 0.338∘

Xeon(R) CPU E5-2699 v4 @2.20GHz) are shown in Fig. 9(a) for memory usage and (b) for indexing time per pattern. At the lowest bandwidth, the indexing time becomes limited by disk access speeds when the patterns are read from a data file. Replacing reading from disk with a constant flow of patterns from a detector system will circumvent this disk access bottle neck, and indexing rates well above 3,000 patterns per second should be feasible, making the SI approach competitive with, but significantly more robust than, standard Hough-based

4. Discussion and conclusions The spherical indexing (SI) algorithm presented can robustly handle 9

Ultramicroscopy 207 (2019) 112841

W.C. Lenthe, et al.

Table 2 Error distribution measurements are tabulated for indexing simulated patterns at a range of noise levels for a bandwidth of 88 demonstrating excellent precision and accuracy. The last 2 columns are computed for a noise level of 1 partitioned about 10 ∘. Noise (N)

None

64

16

4

1

1 ( < 10∘)

1 ( > 10∘)

Mean Median Std. Dev. Max SNR PSNR

0.257∘ 0.250∘ 0.134∘ 0.655∘

0.274∘ 0.267∘ 0.133∘ 0.727∘ 13.5 dB 22.6 dB

0.311∘ 0.309∘ 0.148∘ 0.844∘ 7.53 dB 16.6 dB

0.433∘ 0.421∘ 0.199∘ 1.183∘ 1.51 dB 10.6 dB

13.414∘ 0.904∘ 19.953∘ 61.410∘ -4.53 dB 4.53 dB

0.751∘ 0.681∘ 0.550∘

41.841∘ 43.187∘ 11.077∘

pace with many EBSD detectors on a laptop and all but the fastest on a modest workstation, enabling real time forward model based indexing. The major weakness of SI compared to DI is the inability to capture variation in energy distribution across the detector. Although these

patterns that are difficult or impossible to index with traditional Hough indexing techniques. SI achieves quality comparable to dictionary indexing (DI) with a significant improvement in execution speed. At low bandwidths the SI implementation presented is fast enough to keep

Fig. 8. The full 3.4 million pattern shot-peened Al dataset is represented as an [001] IPF map; each vertical band corresponds to different indexing parameters or algorithm. For the Spherical Indexing approach, results for six different bandwidths are shown. The field of view measures 56 µm × 52 µm. 10

Ultramicroscopy 207 (2019) 112841

W.C. Lenthe, et al.

Table 3 Indexing rates (in patterns per second, pps) on a 4 core laptop [Intel Core CPU i7-7820HQ] and a 24 core workstation [Intel Xeon CPU E5-2670 v3] are given for a range of bandwidths. Note that these rates include the orientation refinement step described in Section 2.5. Bandwidth

Laptop

Workstation

53 63 74 88 113 158

216 pps 113 pps 71.1 pps 36.8 pps 19.0 pps 5.61 pps

1450 pps 860 pps 506 pps 298 pps 134 pps 45.6 pps

Fig. 10. Scaling of the indexing rate (a) and the memory usage (b) at a bandwidth of 88. The two curves in (a) correspond to two different ways in which the FFTW library can be compiled.

difference between a triclinic and a cubic indexing run. Explicit incorporation of crystallographic symmetry in the cross-correlation algorithm does provide a significant symmetry-dependent speed up by eliminating consideration of master pattern spherical harmonic coefficients that are zero by symmetry. Thus, there are indexing speed differences between individual point group symmetries. As a rule, the indexing speed improves with the order of the principal rotation axis in the point group; thus, EBSD patterns for hexagonal structures will see the largest speed-up. With respect to a spherical correlation computation that does not take any symmetry into account, our algorithm includes the following speed up factors: (1) execution of a real-valued 3D FFT skipping β > π (2 × dot product, 3 × FFT); (2) symmetry of real-valued spherical harmonic coefficients (4 × dot product); (3) efficient simultaneous computation of complex products of the form (a + bi)(c + d i) and (a + bi)(c − d i) (2 × dot product); (4) symmetry speed up due to an nfold rotation axis at the master pattern North pole (n × for both the dot product and FFT); and (5) symmetry speed up for a mirror plane at the equator (2 × dot product). Efficient loop ordering to improve caching behavior balanced with a reduction of duplicate operations provides an additional speed up. The overall speed improvements over a naive implementation of the spherical cross correlation is at least a factor of 16 × for triclinic symmetry, and more than 64 × for tetragonal and higher symmetries (for the dot product; 3 × and 12 × respectively for the FFT). This sets the current implementation apart from applications in other fields (e.g., astronomy [18]) where these symmetry-based speed ups are not generally possible.

Fig. 9. Memory (a) and indexing rate (b) scaling for indexing at a range of bandwidths using 48 compute threads. Below bandwidths of 53 for the shotpeened Al data set, the indexing times become disk limited on the system used for these computations.

variations are small they may be important to capture for some applications, e.g. strain measurements. In these cases SI can be used to determine orientation followed by existing refinement techniques developed for DI [37]. While the DI approach produces high quality orientation maps even for noisy EBSD patterns, the drawback of the algorithm is the need to compute a dictionary of patterns; the size of the dictionary increases inversely proportional to the order of the crystallographic rotational point group for the given crystal structure, so that a triclinic DI run takes about 24 times longer than a cubic run. In the new SI approach, on the other hand, there is no need to compute individual patterns since experimental patterns are correlated directly against the EBSD master pattern on the Kikuchi sphere. To first order, there is thus no longer any 11

Ultramicroscopy 207 (2019) 112841

W.C. Lenthe, et al.

Diffraction (TKD). In the TEM, Precession Electron Diffraction can potentially be carried out using the spherical indexing algorithm. Both transmission and reflection Laue X-ray diffraction can be described on a sphere and are hence amenable to spherical indexing. We are currently exploring all these techniques to evaluate the applicability of spherical indexing. A major requirement is the availability of an accurate forward model for the computation of the scattered intensity on a spherical surface. The spherical indexing algorithm and forward models for many of these modalities will be made available in release 4.3 of the open source EMsoft package (fortran-90 version, [38]); the new SI indexing algorithm will also be made available for non-commercial use as a standalone open source C++ package via GitHub [39].

Another novel aspect of the new SI algorithm is the fact that it is fully compatible with the use of a square Lambert projection for the storage of the master pattern. The original computation of the master pattern employs the equal-area sphere-to-square projection described in [17]; the lattitudes on the sphere can be replaced readily by an appropriate set of Legendre lattitudes that are optimal for the implementation of the spherical harmonic transform. The use of special windowing functions to place the experimental patterns on the sphere can be avoided by setting the mean intensity of the pattern to zero before padding the sphere with zeroes. This simplifies the handling of experimental patterns and improves the crosscorrelations by eliminating (or at least reducing) the issue of bright regions of the master pattern skewing the solution towards zone axis orientations. The refinement procedure described in Section 2.5 simultaneously refines all three Euler angles (α, β, γ) by making use of the explicit derivative expressions for the Wigner d function described in the Appendix. This represents an improvement over refinement procedures available in the literature [19] that alternate over α + γ and β to converge on the correct solution. Finally, the spherical indexing technique as described in this paper can be applied to other diffraction modalities; whenever the diffracted intensity can be represented on a sphere, the algorithm can be applied without any significant modifications. For SEM-based techniques, this includes Electron Channeling Patterns (ECP) and Transmission Kikuchi

Acknowledgements The authors would like to thank Dr. S. Wright for making the Ni data sets available, and Drs. B. Winiarski and T. Burnett for the shotpeened Al data set. The authors gratefully acknowledge funding from a DoD Vannevar-Bush Faculty Fellowship (# N00014-16-1-2821) as well as the computational facilities of the Materials Characterization Facility at CMU under grant # MCF-677785. The work of SS was performed partially under the auspices of the U.S. Department of Energy at Lawrence Livermore National Laboratory under Contract DE-AC5207NA27344 (LLNL-JRNL-773006)

Appendix A. Derivatives of the Wigner D function The computation of the α and γ partial derivatives of the Wigner D function

Dkℓ, m (α,

β, γ ) = dkℓ, m (β )eimα eikγ

(A.1)

is straightforward; however, the β partials require computation of the derivatives of the Wigner d function:

dkℓ, m (β ) =

β k+m β k−m (ℓ + k ) ! (ℓ − k )! cos ⎛ ⎞ sin ⎛ ⎞ Pℓk−−km, k + m (cos(β )) (ℓ + m) ! (ℓ − m)! ⎝2⎠ ⎝2⎠

(A.2)

which is considerably more complex. To the authors’ knowledge, closed-form expressions for the β partial derivatives of reported previously. The derivatives of the Jacobi polynomial can be expressed recursively as

Γ(k + m + ℓ + 1 + n) k + j, m + j ∂n k , m Pℓ (x ) = n Pℓ − n (x ) 2 Γ(k + m + ℓ + 1) ∂x n

=

dkℓ, m (β )

have not been

(A.3)

(k + m + ℓ + n)! k + j, m + j Pℓ − n (x ) 2n (k + m + ℓ)!

(A.4)

For the present pattern indexing application, we will assume that k, m, ℓ, and n are always non-negative integers. Note that some applications of spherical harmonics use half integers which would enable a different simplification. Using the expression for the derivatives of the Jacobi polynomial, compact recursive expressions for the β partial derivatives of the Wigner d function can be derived:

∂ ℓ dk, m (β ) = x 0 dkℓ, m (β ) − x1 dkℓ,+m1 (β ) ∂β

(A.5)

x 0 = (cos(β ) k − m)csc(β )

(A.6)

x1 =

(A.7)

(ℓ − k )(ℓ + k + 1)

∂2 ℓ dk, m (β ) = y0 dkℓ, m (β ) − y1 dkℓ,+m1 (β ) + y2 dkℓ,+m2 (β ) ∂β 2

(A.8)

y0 = (k 2 cos2 (β ) + m (1 − 2k )cos(β ) + (m2 − k ))csc2 (β )

(A.9)

y1 =

(ℓ − k )(ℓ + k + 1) (cos(β )(1 + 2k ) − 2m)csc(β )

y2 =

(ℓ − k )(ℓ + k + 1) (ℓ − k − 1)(ℓ + k + 2)

(A.10) (A.11)

These forms are particularly advantageous numerically since the evaluation of cos (β) is already required evaluating another trigonometric function via the relation:

12

for dkj, m

and csc(β ) can be computed without

Ultramicroscopy 207 (2019) 112841

W.C. Lenthe, et al.

csc(β ) =

sign(β ) 1 − cos2 (β )

(A.12)

If derivatives are required for all k and m values, they can be computed for non-negative k and m only and symmetry can be used to compute the other coefficients.

∂ ℓ ∂ d−k, m (β )=(−1)ℓ + m + 1 dkℓ, m (π −β ) ∂β ∂β

(A.13)

∂ ℓ ∂ dk, −m (β )=(−1)ℓ + k + 1 dkℓ, m (π −β ) ∂β ∂β

(A.14)

∂ ℓ ∂ d−k, −m (β )=(−1) k + m dkℓ, m (β ) ∂β ∂β

(A.15)

∂2 ℓ ∂2 d−k, m (β )=(−1)ℓ + m 2 dkℓ, m (π −β ) ∂β 2 ∂β

(A.16)

∂2 ℓ ∂2 dk, −m (β )=(−1)ℓ + k 2 dkℓ, m (π −β ) ∂β 2 ∂β

(A.17)

∂2 ℓ ∂2 d−k, −m (β )=(−1) k + m 2 dkℓ, m (β ) ∂β 2 ∂β

(A.18)

Combined, these expressions enable efficient tabulation of all first and second partial derivatives of the Wigner D function for all ℓ, k, and m values with respect to α, β, and γ. Supplementary material Supplementary material associated with this article can be found, in the online version, at 10.1016/j.ultramic.2019.112841.

Vol.3 [20] L. Sorgi, K. Daniilidis, Template gradient matching in spherical images, Image Processing: Algorithms and Systems III, vol. 5298, International Society for Optics and Photonics, 2004, pp. 88–98. [21] B. Gutman, Y. Wang, T. Chan, P.M. Thompson, A.W. Toga, Shape registration with spherical cross correlation, 2nd MICCAI Workshop on Mathematical Foundations of Computational Anatomy, (2008), pp. 56–67. [22] R. Hielscher, F. Bartel, T.B. Britton, Gazing at crystal balls - electron backscatter diffraction indexing and cross correlation on a sphere, under review in Ultramicroscopy, arXiv e-prints (2018) arXiv:1810.03211. [23] S. Pizer, E. Amburn, J. Austin, R. Cromartie, A. Geselowitz, T. Greer, B. Romney, J. Zimmerman, K. Zuiderveld, Adaptive histogram equalization and its variation, Comput. Vis. Graph. Image Process. 39 (1987) 355–368. [24] M. Reinecke, Libpsht–algorithms for efficient spherical harmonic transforms, Astron. Astrophys. 526 (2011) A108. [25] N. Sneeuw, Global spherical harmonic analysis by least-squares and numerical quadrature methods in historical perspective, Geophys. J. Int. 118 (3) (1994) 707–716, https://doi. org/10.1111/j.1365-246X.1994.tb03995.x. [26] P.J. Kostelec, D.N. Rockmore, S2kit: A Lite Version of SpharmonicKit, (2004). (https:// www.cs.dartmouth.edu/geelong/sphere/ ) [27] J.R. Driscoll, D.M. Healy, Computing fourier transforms and convolutions on the 2sphere, Adv. Appl. Math. 15 (2) (1994) 202–250. [28] N. Schaeffer, Efficient spherical harmonic transforms aimed at pseudospectral numerical simulations, Geochem. Geophys. Geosyst. 14 (3) (2013) 751–758. [29] D. Roşca, New uniform grids on the sphere, Astron. Astrophys. 520 (2010) A63. [30] P.J. Kostelec, D.N. Rockmore, Ffts on the rotation group, J. Fourier Anal. Appl. 14 (2) (2008) 145–179. [31] T. Fukushima, Numerical Computation of Wigner’s D-Function of Arbitrary High Degree and Orders by Extending Exponent of Floating Point Numbers, Technical Report, National Astronomical Observatory of Japan, 2016, https://doi.org/10.13140/RG.2.2.31922. 20160. [32] B. Huhle, T. Schairer, W. Strasser, Normalized cross-correlation using soft, 2009 International Workshop on Local and Non-Local Approximation in Image Processing, (2009), pp. 82–86, https://doi.org/10.1109/LNLA.2009.5278398. [33] F. Ram, M. De Graef, Energy dependence of the spatial distribution of inelastically scattered electrons in backscatter electron diffraction, Phys. Rev. B 97 (13) (2018) 134104. [34] S. Wright, Private Communication, 2018. [35] P.T. Brewick, S.I. Wright, D.J. Rowenhorst, NLPAR: non-local smoothing for enhanced EBSD pattern indexing, Ultramicroscopy 200 (2019) 50–61. [36] M. Frigo, S.G. Johnson, FFTW: an adaptive software architecture for the FFT, Proceedings of the 1998 IEEE International Conference on Acoustics, Speech and Signal Processing, ICASSP’98 (Cat. No. 98CH36181), vol. 3, IEEE, 1998, pp. 1381–1384. [37] S. Singh, F. Ram, M. De Graef, Application of forward models to crystal orientation refinement, J. Appl. Crystallogr. 50 (6) (2017) 1664–1676. [38] M. De Graef, EMsoft, (2019). (http://github.com/EMsoft-org/EMsoft ) [39] W. Lenthe, EMSphInx, C++ Implementation, (2019). (http://github.com/EMsoft-org )

References [1] N. Krieger, Image processing procedures for analysis of electron backscattering patterns, Scanning Microsc. 6 (1992) 115–121. [2] B.L. Adams, S.I. Wright, K. Kunze, Orientation imaging: the emergence of a new microscopy, Metall. Trans. A 24 (1993) 819–831. [3] N. Krieger Lassen, Automated Determination of Crystal Orientations from Electron Backscattering Patterns, The Technical University of Denmark, 1994 Ph.D. Thesis. [4] S. Wright, M. Nowell, S. Lindeman, P. Camus, M. De Graef, M. Jackson, Introduction and comparison of new EBSD post-processing methodologies, Ultramicroscopy 159 (2015) 81–94. [5] Y. Chen, P. S.U., D. Wei, G. Newstadt, M. Jackson, J. Simmons, M. De Graef, A. Hero, A dictionary approach to EBSD indexing, Microsc. MicroAnal. 21 (2015) 739–752. [6] F. Ram, S. Singh, S. Wright, M. De Graef, Error analysis of crystal orientations obtained by the dictionary approach to EBSD indexing, Ultramicroscopy 181 (2017) 17–26. [7] K. Marquardt, M. De Graef, S. Singh, H. Marquardt, A. Rosenthal, T. Hiraga, Quantitative electron back-scatter diffraction (EBSD) data analyses using the dictionary indexing (DI) approach: overcoming indexing difficulties in geological materials, Am. Mineral. 102 (2017) 1843–1855. [8] S. Singh, Y. Guo, B. Winiarski, T. Burnett, P. Withers, M. De Graef, High resolution low kv EBSD of heavily deformed and nanocrystalline Aluminium by dictionary-based indexing, Nat. Sci. Rep. 8 (2018) 10991. [9] F. Ram, M. De Graef, Phase differentiation by electron backscatter diffraction using the dictionary indexing approach, Acta Mater. 144 (2018) 352–364. [10] D. Rowenhorst, A. Rollett, G. Roher, M. Groeber, M. Jackson, P. Konijnenberg, M. De Graef, Tutorial: consistent representations of and conversions between 3D rotations, Model. Simul. Mater. Sci.Eng. 23 (2015) 083501. [11] S. Singh, M. De Graef, Orientation sampling for dictionary-based diffraction pattern indexing methods, Model. Simul. Mater. Sci. Eng. 24 (2016) 085013. [12] A.J. Wilkinson, D.M. Collins, Y. Zayachuk, R. Korla, A. Vilalta-Clemente, Applications of multivariate statistical methods and simulation libraries to analysis of electron backscatter diffraction and transmission Kikuchi diffraction datasets, Ultramicroscopy 196 (2019) 88–98. [13] G. Nolze, A. Winkelmann, A.P. Boyle, Pattern matching approach to pseudosymmetry problems in electron backscatter diffraction, Ultramicroscopy 160 (2016) 146–154. [14] T. Tanaka, A.J. Wilkinson, Pattern matching analysis of electron backscatter diffraction patterns for pattern centre, crystal orientation and absolute elastic strain determination–accuracy and precision assessment, Ultramicroscopy 202 (2019) 87–99. [15] P. Callahan, M. De Graef, Dynamical EBSD patterns Part I: pattern simulations, Microsc. MicroAnal. 19 (2013) 1255–1265. [16] A. Day, Spherical EBSD, J. Microsc. 230 (2008) 472–486. [17] D. Roşca, New uniform grids on the sphere, Astron. Astrophys. 520 (2010) A63. [18] P. Lemos, A. Challinor, G. Efstathiou, The effect of Limber and flat-sky approximations on galaxy weak lensing, JCAP 5 (2017) 14. [19] A. Makadia, L. Sorgi, K. Daniilidis, Rotation estimation from spherical images, Proceedings of the 17th International Conference on Pattern Recognition, 2004. ICPR 2004. volume 3, (2004), pp. 590–593, https://doi.org/10.1109/ICPR.2004.1334598.

13