# A generic topography reconstruction method based on multi-detector back scattered electron images

Jan Neggers, Eva Héripré, Marc Bonnet, Simon Hallais, Stéphane Roux

► To cite this version:

Jan Neggers, Eva Héripré, Marc Bonnet, Simon Hallais, Stéphane Roux. A generic topography reconstruction method based on multi-detector back scattered electron images. Strain, 2022, 58 (5), pp.1625-1646. 10.1111/str.12416. hal-04502013

HAL Id: hal-04502013

https://hal.science/hal-04502013v1

Submitted on 13 Mar 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
ARTICLE TYPE

# A generic topography reconstruction method based on multi-detector back scattered electron images

Jan Neggers*<fnref n="1" />,<fnref n="2" /> | Eva Héripré<fnref n="1" /> | Marc Bonnet<fnref n="2" /> | Simon Hallais<fnref n="3" /> | Stéphane Roux<fnref n="2" />

Correspondence

*Jan Neggers, CentraleSupélec, Email: jan.neggers@centralesupelec.fr

Abstract

Surface topographies can be reconstructed from backscattered electron (BSE) images captured from different detector orientations. This article presents a very general approach to this problem, in the spirit of photometric stereo methods, allowing for arbitrary BSE detector number (at least 3) and shapes. The general idea is to both determine the (non-linear) model parameters and compute the surface topography so that the modeled images match at best the acquired ones. Three samples are used for validation of the measured topography with respect to AFM measurements. RMS errors in the range of 10-35 nm, or 1-1.5% of total sample height, are obtained.

KEYWORDS:

BackScattered electron images; Topography; SEM;

1 | INTRODUCTION

The evolution of the 3D topography of a sample surface during mechanical loading is of great interest for understanding the mechanical behavior of its constituting material (e.g.,<fnref n="1" />,<fnref n="2" />). At small scales, this topography information can be even more precious as micro mechanical phenomena tend to influence surface roughness<fnref n="3" />. The natural tools that are available for studying the surface topography are Scanning Probe Microscopes (SPM), such as the Atomic Force Microscopy (AFM). In this article, a method is proposed to extract 3D topography data from multiple Back-Scattered Electron (BSE) detectors in a Scanning Electron Microscope (SEM). SEM acquisitions are common with in-situ experimentation and generally more widely applicable at a wide range of imaging length scales. This method does not compete with AFM topography measurements in terms of resolution. However, occasionally, no AFM is available, or the desired field of view is too large for an AFM. Moreover, successful acquisitions of AFM measurements in in-situ mechanical tests can be extremely challenging⁴ and time-consuming. In contrast, simultaneous acquisitions on different BSE detectors can be obtained without additional manipulation of the observed specimen, thus making time series very accessible.

Obtaining the topography using SEMs is nothing new, and in general SEM topography reconstruction methods can be classi-
fied into two groups; 1) Photogrammetry⁵,⁶ and 2) Photometric Stereo⁷,⁸,⁹,¹⁰. The first class of methods uses the parallax created
by observing the sample from multiple points of view to reconstruct the topography. For SEMs, this requires tilting the sample
to capture a second acquisition. Tilting the sample back and forth is not always feasible, especially with an in-situ tester mounted
on the tilt-stage. Moreover, this procedure augments the total acquisition time significantly making the method less useful for
time-sensitive tests. The second class of methods relies on how the sample surface changes intensity as a function of the detec-
tor location. This reduction in intensity can be due to two physical phenomena. First, the observed intensity can be reduced
due to partial occlusion of the detector from the point of incidence. Second, the intensity changes due to the surface orientation
directing more electrons to some directions. The phenomena is present for nearly all types of detectors found in SEMs¹¹,¹²,¹³.

<footnote n="1">¹ CentraleSupélec / CNRS / Université Paris-Saclay, MSSMat, F-91190, Gif-sur-Yvette, France</footnote>

<footnote n="2">² Université Paris-Saclay / ENS Paris-Saclay / CNRS, LMT - Laboratoire de Mécanique et Technologie, F-91190, Gif-sur-Yvette, France</footnote>

<footnote n="3">³ CNRS / École Polytechnique, LMS - Laboratoire de Mécanique et Technologie, F-91128, Palaiseau, France</footnote>
However, to faithfully reconstruct the surface, this gray-level intensity has to be converted to surface orientation. Consequently, it has to be predictable.

The method proposed in this article falls in the second class of methods. Specifically, it has been developed for multi-segment BSE detectors. Classically, four-quadrant BSE detectors have been used for this. Because the four segments in such a detector have a similar sensitivity, the difference of two diagonally opposed detectors is a good proxy for the surface tilt along that diagonal. Quite a number of interesting articles have been devoted to these four-quadrant setups<sup>14,15,16,17</sup>. However, these four-quadrant detectors are less common in today's SEMs. Currently, detectors with multiple annular segments where each ring can be segmented in 120° segments are also found, e.g., figure 1. These detectors still provide sufficient information to perform the reconstruction, however, they require a more general algorithm. In photometric stereo methods, it is common to rotate the sample to generate additional detector locations<sup>18</sup>. This is equally possible within the SEM and would naturally enhance the reconstruction<sup>19</sup>. However, for the presented work, it was chosen to not use sample motion because in *in-situ* settings, this is often not possible. Moreover, capturing all detectors with a single scan has the advantage that all images share the same scanning artifacts and image distortions typical in SEMs.

Recently, a simple algorithm was proposed that decomposes the set of multi-detector images into principal component images<sup>20</sup>. The principal component images represent a linearization of the problem with respect to its proper basis, in the same spirit as most methods, but could handle an arbitrary set of BSE detectors. This type of linearization can obviously be questioned.

The starting point of the proposed method can be seen as following the same spirit as photometric stereo. Indeed, it is natural to assume initially that the detected back-scattered electrons have undergone sufficiently many elastic scatterings that the influence of the incident direction is completely lost. It is known from experiments that this is not strictly true. For instance, crystalline materials exhibit special channeling and diffraction angles, and even in general, the scattering distribution tends to be more complex than pure Lambertian<sup>21</sup>. However, the needed integration of the response over the entire extent of the BSE detectors by itself transform the response into a more complex model than an “extented” Lambertian model, and the resulting approach would inherit from the fragile assumption of an ideal multiple scattering regime. Rather, the large solid angle of the detector, its position, and its offset and gain, together with the effective interaction of the incident electrons with the surface will all be assembled into an effective response, which is designed to be robust. Consequently, the method will be more accurate for amorphous materials with moderate surface slopes, even if the assumption of a perfectly Lambertian surface will *not* be assumed to be valid. The methodology proposed here is completely generic. However, no chemical contrast is considered here. A difference in chemical composition has a multiplicative effect on the observed intensity. This would require an additional degree of freedom per pixel. For this article, the choice was made to not include this in the algorithm. The consequence is that the presented method is only valid for homogeneous materials and where channeling contrast can be safely ignored (nanocrystalline or amorphous materials).

The method discussed in the paper is using an explicit modeling of the observed detector intensity. These models are parameterized equations that predict the observed intensity for a given topography for a given detector for each pixel. An inverse method is applied to obtain the set of parameters that describe the detector configuration while at the same time finding the topography that minimized the difference with the actual acquisitions.

Interestingly, some authors have successfully applied the photometric stereo method by considering the detectors as optical light sources and the incident beam as the point of view<sup>22,23</sup>. The optical metaphor allows the connection to existing photometric stereo algorithms. However, it relies on the Lambertian scattering assumption to work. Conversely, our proposed algorithm uses the Lambertian assumption only for the introduction of a first simple model, which is further extended to more general cases.

Section 2 introduces the construction of the modeling of the detector response based on very general considerations, and assumptions that are progressively relaxed. The algorithm used for the joint determination of the surface topography and of the model parameters is then detailed. In Section 4, the application of the proposed method to three samples is discussed together with a quantitative comparison with corresponding AFM measurements. Eventually, the main conclusions are summarized in Section 5.

## 2 | METHOD

### 2.1 | Image formation

As the electronic beam is focused on the surface at position **x**, detector *i* collects an intensity which is denoted $f^{(i)}(\mathbf{x})$. The target of this study is to reconstruct the surface topography $z = h(\mathbf{x})$. The first step is to formulate a generic model that describes
<figure><img src="image_2.png" /><figcaption>**FIGURE 1** An example of a typical BSE multi-detector sold by ThermoFisher. It has 8 segments which are grouped in 4 clusters, here the ABS variant is shown and used for the cases studied in this article which groups the segments angularly (instead of annularly for the CBS variant)</figcaption></figure>

the response $\hat{f}^{(i)}$ based on characteristics of the surface that has been hit, so that the model based images $\hat{f}^{(i)}$ match the data acquisition $f^{(i)}$.

Because detectors are wide, in general of non-simple shape (such as a ring segment), and may be affected by different acquisition artifacts (e.g., offset or gain being imperfectly controlled, or slightly different from one detector to another one) it is difficult if not impossible to resort to first principles to compute $\hat{f}^{(i)}$.

It is instead proposed to measure the response of each detector. This may appear as a kind of bootstrap problem since the final goal is to measure the topography of a sample, and the starting point is to measure the response of the detector with no more information than images of a surface of unknown topography. However, the following sections will show that such an approach may work.

The first step consists in formulating a generic response model, and hence a number of assumptions are first formulated to make the problem as simple as possible. Then this first framework can be later extended in different directions, relaxing the initial restrictions.

All detectors are supposed to obey the same generic model, although with potentially different model parameters. The detector label (i) written as a superscript $f^{(i)}$, is thus dropped in this section, where only one detector is considered first. Let us first assume that the response is (H1) local, i.e., solely dependent on the sample surface at the location where the electron beam has been focused on, and (H2) "translationally invariant" (over the range of the observed region). This latter property (H2) implies that the response $\hat{f}$ does not depend per se on the position $x$ which is targeted, but only on the local surface topography. This assumption obviously restricts this first analysis to single-phase materials where chemical and crystallographic contrast can be ignored.

As a further simplifying assumption, it is considered that (H3) the topography can be considered as a “small perturbation” about a flat surface, that may conventionally be written $z = 0$. These hypotheses lead to the consequence that the response can be described by its Taylor expansion to first order, namely

$$ \hat{f}(h) = p_0 + p_1 h + \underline{p_2} \cdot \underline{\nabla h} + \underline{p_3} \cdot \underline{\nabla \otimes \nabla h} + h.o.t. \quad (1) $$

Taking into account derivatives at different orders is necessary to define properly what is meant by small perturbations. Indeed, surface elevation $h$ may remain very small, while its slope $|\nabla h|$ may be large, and the same would hold for curvature. This question will be further addressed below when discussing the different parameters.

$p_0$ describes the response of the detector for the reference perfect plane $h = 0$. It is obviously an important parameter to include in the description of the detector, although it does not inform on the topography.

$p_1$ characterizes the change in response as the surface is translated normal to the beam. The range of variation for $h$ considered here (typically in the micrometer range at most) is very small as compared to other physical length scales relevant here such as working distance or detector width, (several mm) so that it is proposed to neglect this contribution, and set $p_1 = 0$.

$p_2$ describes the variation of the response with the orientation of the normal to the surface. It is well known that this orientation has a major influence on the observed BSE image gray level. Being the first non-zero parameter order (apart from the trivial $p_0$
which does not depend on topography), this 2D vector parameter is the most important characteristic of a detector that will be used in the sequel.

$p_3$ captures the influence of curvature on the detector response. Dimensional analysis indicates that a ratio of any $p_3$ component to $|\overline{p_2}|$ one provides a length scale $\lambda$ such that radius of curvature smaller than $\lambda$ are needed to make this term comparable to the previous. Curvature of the surface can only be felt at the scale of penetration $\xi$ of backscattered electrons in the material (say in the range 10-100 nm). Such curvatures would signal very rough surfaces which are not expected to be resolvable by such a technique. It is therefore natural to neglect those terms, and set them to 0, $p_3 = 0$. It is even more justified to neglect all other higher order terms for similar reasons. This is quite fortunate since, as higher and higher order $j$ terms are considered, the number of components of $p_j$ increases ($p_1$ is a scalar, $p_2$ a 2D vector, $p_3$ a symmetric tensor, ...), and hence the number of model parameters would increase very fast with maximum considered order. However, locally, high curvature may be encountered at specific locations over the surface, and hence it is possible to include such curvature corrections to the dominant signal which is controlled by the orientation. Because this is a correction, one may further introduce a hierarchy in the different components of $p_3$: the curvature tensor may be decomposed in a spherical component, the mean curvature, (given by the trace of the tensor), and a deviatoric part (trace-less). The former is isotropic, while the latter is turned into its opposite when the surface is rotated by $\pi/2$ about the z-axis. Integration of the signal over a wide detector naturally dampens the deviatoric contribution, whereas the mean curvature is unchanged. Thus a minimal refinement including curvature would involve only the Laplacian of the topography, $\nabla^2 h = \partial_x^2 h + \partial_y^2 h$, which is a scalar, and hence the conjugate parameter $p_3$ is also a scalar.

Therefore, the proposed generic linear model is written

$$ \hat{f}(h) = p_0 + \underline{p_2} \cdot \overline{\nabla h} + p_3 \overline{\nabla^2 h} \quad (2) $$

involving three parameters, $p_0$ and $\underline{p_2}$, and optionally a fourth one $p_3$.

Let us propose a simple geometrical picture of the proposed model: We have seen that basically only the surface normal $\underline{n}$ matters. This result holds for any point $M$ on the detector. Similarly, any such point is essentially at a fixed distance from any point $S$ of the sample surface (the field of view on the sample is very small as compared to the detector mean distance), or, introducing the origin $O$ of the sample reference plane, $|OS| \ll |OM|$, and hence $\overline{SM} \approx \overline{OM}$. However, the detector width is such that the unit orientation vector $\underline{m} = \overline{OM}/|OM|$ may explore a wide solid angle. Now if both surface and detector rotate about the z-axis, the response is expected not to change, and hence the response of the detector depends not on $\underline{n}$ but only on the relative angle between $\underline{m}$ and $\underline{n}$. Rather than discussing 3D orientations or points on the unit sphere, it is much more convenient to project $\underline{n}$ and $\underline{m}$ on the plane $z = 0$, defining respectively the two vectors $\underline{n}_2$ and $\underline{m}_2$. The former (in the small slope approximation) is identical to $\underline{n}_2 = \overline{\nabla h}$. The latter $\underline{m}_2$, as $M$ is varied to all positions within detector $i$ draws on the plane a domain $\Omega^{(i)}$ that looks like the shape of the detector. Let us stress the fact that the global detector response is the integral over the entire domain $\Omega^{(i)}$ of a elementary local response. Thus, it is obvious that the azimuth sensitivity with respect to $\underline{n}$, i.e., the polar angle of $\underline{n}_2$ is small, because of the aperture effect of a wide detector.

$\underline{p_2}$ is also a vector which can be shown in the $(x, y)$ plane. Its orientation is naturally expected to be directed towards the center of mass $\overline{M}^{(i)}$ of $\Omega^{(i)}$. It is therefore natural to introduce the polar angle $\theta^{(i)}$ of $\overline{M}^{(i)}$. It may appear convenient to rotate the surface topography gradient, $\overline{\nabla h}$, by an angle $\theta^{(i)}$, so that the resulting vector $R(\theta^{(i)})\overline{\nabla h}$ is written $(\partial_{\theta^{(i)}}h, \partial_{\theta^{(i)}}h)$, where the first component is the slope of the surface in the direction of the detector (on average), and the second is perpendicular to it. Introducing this rotation in the above model yields

$$ f(h) = p_0 + \tilde{p}_2 \cdot \overline{R(\theta)} \cdot \overline{\nabla h} + p_3 \overline{\nabla^2 h} \quad (3) $$

where $\tilde{p}_2 = \overline{R(-\theta)} \cdot \underline{p_2}$ defines a rotated vector. The reason why this expression is convenient, is that for a detector which is symmetric with respect to the plane containing the z axis and the point $\overline{M}^{(i)}$, $\partial_{\theta^{(i)}}h$ does not influence the response for this linear model, and hence $(\tilde{p}_2)_t = 0$. In the following, this rotation is considered, but the assumption $(\tilde{p}_2)_t = 0$ is only used for initialization. Then upon further calibration of the model $\tilde{p}_2$ is treated as a full vector. It is to be noted that this is mathematically completely equivalent to the initial writing Eq. 2, and rotation is not necessary. It allows however to correct for a poor initial determination of the angle $\theta^{(i)}$, or for detector with a complex (non-symmetric) geometry, it could even allow to define $\theta^{(i)}$, from the condition that $(\tilde{p}_2)_t = 0$. Or, if this direction is known (e.g. for having been calibrated earlier), then one may take advantage of the condition $(\tilde{p}_2)_t = 0$ to reduce the number of parameters of the model to no more than 2 (per detector).
## 2.2 | Generalization

Let us now re-discuss the assumptions that were listed above in the construction of our model to see whether one could extend or generalize the proposed modeling.

### (H1): Locality

Backscattered electrons stem from a small neighborhood (of order of $\xi$) of the focal point of the beam because of multiple scattering. This is the reason why the response may depend on slope and curvature of the surface. As above discussed, any topographical feature below this scale $\xi$, cannot be resolved, and when gradient and curvature are considered, they are to be computed over a low-pass filtered surface topography with a cut-off frequency of $1/\xi$. Apart from this caveat, locality may be questioned in the case of shadowing. For very steep topographies, some part of the surface may appear hidden from some surface elements of a detector, because say of a salient bump. In this case, it is not only the topographical characteristics of the local surface point that dictates the detector response. This difficulty may possibly be tackled with techniques proposed for the “photometric stereo” problem$^8$, but it goes much beyond the scope of the present paper.

### (H2): Translational invariance

For a wide field of view, the direction of the backscattered electrons from the sample, $\overline{SM}$ can no longer be confused with $\overline{OM}$. Thus a small correction arises in the determination of $\partial_n h$, proportional to the projection of $\overline{OS}$ on $e_\theta$. The proportionality constant is dependent on the orientation of the detector and hence cannot be quantified precisely a priori. Similarly a linear correction is expected when the studied point is at elevation $z = h$. This effect leads to

$$ \hat{f}(\underline{x}, h) = p_0 + \bar{p}_2 \cdot (R(\theta) \cdot \underline{\nabla} h + (\alpha \underline{x} \cdot e_{\theta} + \beta h)e_{\theta}) + p_3 \underline{\nabla}^2 h \quad (4) $$

where the parameters $\alpha^{(i)}$ and $\beta^{(i)}$ could be included in the parameter list. If this correction is not taken into account, a spurious curvature may arise in the reconstructed surface. For conciseness the correction is denoted as

$$ \underline{c}(\underline{x}, h) \equiv (\alpha \underline{x} \cdot e_{\theta} + \beta h)e_{\theta} \quad (5) $$

### (H3): Small perturbation

The most obvious limitation of the above linear model is the assumption of slight perturbations, especially for slopes, which may assume values of order 1. However, there is no way to assess a priori the range of validity of the linear description. However, considering the fact that $\hat{f}$ is mostly a function of $(\partial_n h + c, \partial_t h)$, one may propose higher order correction easily by expanding the first two terms of $\hat{f}$, denoted $\varphi(\partial_n h + c, \partial_t h)$ as a polynomial of higher order

$$ \varphi(u, v) = p_0 + p_{2n}u + p_{2t}v + q_{nn}u^2 + 2q_{nt}uv + q_{tt}v^2 + r_{nnn}u^3 + 3r_{nnt}u^2v + 3r_{ntt}uv^2 + r_{ttt}v^3 + \dots \quad (6) $$

or equivalently using a vector and tensor notation, with $\underline{u} = (u, v)^T$,

$$ \varphi(\underline{u}) = p_0 + \underline{p}_2 \cdot \underline{u} + \underline{q} \cdot (\underline{u} \otimes \underline{u}) + \underline{r} \cdots (\underline{u} \otimes \underline{u} \otimes \underline{u}) + \dots \quad (7) $$

Moreover, when one can trust the symmetry of the detector with respect to the plane going through its center and the beam axis, then the simplification $p_{2t} = 0$ can be further generalized to $q_{nt} = r_{nnt} = r_{ttt} = 0$ for all coefficients of odd terms in $t$. However as earlier mentioned, such symmetries will not be used in the following analysis. It is also worth noting that, although the proposed generalization opens the way to higher and higher slopes, this is at the expense of a significant inflation of parameters for the model. The reported application cases below will show that expansions up to polynomials of 4th order can be dealt with without instabilities that could have been feared.

To recapitulate, the proposed approach has lead us to the following form for each detector

$$ \hat{f}(h) = \varphi \left( \underline{R}(\theta) \cdot \underline{\nabla} h + \underline{c}(\underline{x}, h) \right) + p_3 \underline{\nabla}^2 h \quad (8) $$

where $\varphi$ is a bivariate polynomial of given order. When the order is larger than 1, this model is non-linear, and when seeking for a small correction $\delta h$ to the current determination of $h$, it is useful to compute the tangent model that describes the change in the BSE image, $\delta \hat{f}$

$$ \begin{aligned} \delta \hat{f} &= \underline{\nabla} \varphi \left( \underline{R}(\theta) \cdot \underline{\nabla} h + \underline{c} \right) \cdot \left( \underline{R}(\theta) \cdot \underline{\nabla} \delta h + \beta \delta h e_{\theta} \right) + p_3 \underline{\nabla}^2 \delta h \\ &\equiv \psi \left( \underline{R}(\theta) \cdot \underline{\nabla} \delta h + \beta \delta h e_{\theta}, h \right) + p_3 \underline{\nabla}^2 \delta h \end{aligned} \quad (9) $$

where the second expression is introduced here to formally underline that the tangent model has the same formal structure as the full model, with $\psi$ being a linear function of its first argument, but with a non trivial vector field (comparable to $p_2$) whose value
depends on the precise location of the pixel through the current determination of the topography (in contrast with $p_2$ which is constant.

Let us finally stress that with the proposed model (when parameters are unknown) there is no way to set an absolute scale. If the topography is scaled by an arbitrary factor $\mu$, and correspondingly the parameters of the polynomial $\varphi$, are scaled by $\mu^{-m}$ where $m$ is the order of the corresponding parameter, then $f$ remains unchanged. However, relative topographies can be trusted. In order to lift the above degeneracy, a convention is adopted. Namely, the mean over all detectors of $|p_2|$ is set to 1. All other parameters and the topography itself are rescaled according to their dimensional homogeneity.

## 2.3 | Connection with previous modeling

Let us assume that the detector obeys perfectly a fourfold rotation symmetry, as expected for a four-quadrant detector, and moreover that each quadrant has precisely the same physical characteristics, and acquisition parameters. Then $\theta^{(k)} = (k-1)\pi/2$ with $k = 1, ..., 4$. Hence the center of mass vectors of the detector point respectively towards $e_x, e_y, -e_x$ and $-e_y$. Moreover, because of the symmetry of each detector with respect to the direction $\overline{OM}^{(k)}$, only the component of $n$ along $^{(k)}$ will be felt, and hence, the linear model with no laplacian dependency can be written

$$ 
\begin{aligned}
f^{(1)}(h) &= p_0 + |p_2|\partial_x h & f^{(2)}(h) &= p_0 + |p_2|\partial_y h \\
f^{(3)}(h) &= p_0 - |p_2|\partial_x h & f^{(4)}(h) &= p_0 - |p_2|\partial_y h
\end{aligned}
\qquad (10) $$

This simple system can easily be inverted to provide

$$ 
\begin{aligned}
\partial_x h &= \frac{p_0 \frac{f^{(1)}(h) - f^{(3)}(h)}{|p_2| \frac{f^{(1)}(h) + f^{(3)}(h)}{|p_2| \frac{f^{(2)}(h) - f^{(4)}(h)}{|p_2| \frac{f^{(2)}(h) + f^{(4)}(h)}{|p_2|}}}}{\frac{f^{(1)}(h) + f^{(3)}(h)}{|p_2| \frac{f^{(2)}(h) + f^{(4)}(h)}{|p_2|}}}
\end{aligned}
\qquad (11) $$

These equations coincide with the analysis proposed by Lebiedzik$^{16}$, and used in Paluszyński and Slówko$^{13}$. Thus the proposed model can be seen as a generalization of this classical approach, to arbitrary detector number and geometry and curvature effects and non-linearities that may arise for steep slopes. Even without non-linearities, or curvature effect, let us underline that the proposed approach offers much more freedom and robustness to the analysis.

## 2.4 | Algorithm

After having described the derivation of the detector response model, the algorithmic aspects are discussed. The general philosophy is to start with a simplistic model as an initial guess. Then from the different detector images a first approximation of the surface topography is computed. This topography is then used to refine the model of each detector considered individually. And the procedure is repeated, alternatively refining topography and detector model, up to a fixed point.

Initialization is set to a trivial model where $p_0^{(i)} = \langle f^{(i)} \rangle$ (which somehow assumes that the surface is globally flat, but $p_0^{(i)}$ could equally well to set to the minimum or the maximum value of $f^{(i)}$, with no prejudice on convergence) and $p_2^{(i)} = e_{\varphi(i)}$, while all other parameters are set to 0. Let us note that the latter condition is consistent with the chosen conventional normalization $\langle |p_2^{(i)}| \rangle_i = 1$.

The algorithm 1 details the main steps of the global procedure. Some of these steps, labeled A, B or C, are further described below.

### 2.4.1 | (A) Computation of gradient and curvature

One detector only delivers a partial information which is obviously not sufficient to determine the entire sample topography. Thus, to evaluate topography gradient or curvature, it is needed to assemble the information coming from all detectors.

Let us consider the general case where the model is non-linear. Topography gradient and curvature are computed as corrections to the current determination of the topography. Therefore, from the current topography (including $h = 0$ in the first pass), residuals of BSE images are computed. These residuals $\delta \hat{f}(\mathbf{x})^{(i)}$ are the difference between the actual acquisition $f(\mathbf{x})^{(i)}$ and the one, $\hat{f}(\mathbf{x})^{(i)}$, computed from the response model of detector $i$. Then the gradient $\gamma(\mathbf{x})$ and curvature $\kappa(\mathbf{x})$ of $\delta h$, treated as three independent (unrelated) scalar fields, are computed from the tangent model, Eq. 9. Because the tangent model is linear in $\delta h$ — hence in $\gamma_x, \gamma_y$ (and $\kappa$) — these estimates are obtained from a redundant linear system, and solved in the least squares sense.
**Algorithm 1** Pseudocode for our algorithm. The three procedures labelled A, B and C are detailed in the main text

<table>
  <tbody>
    <tr>
      <td>procedure TOPOSEM(f<sup>(i)</sup>)</td>
      <td></td>
    </tr>
    <tr>
      <td>Initialize model parameters</td>
      <td></td>
    </tr>
    <tr>
      <td>Compute gradient (and laplacian) of h from model and BSE images f<sup>(i)</sup></td>
      <td>▷ (A1)</td>
    </tr>
    <tr>
      <td>Compute topography h</td>
      <td>▷ (B1)</td>
    </tr>
    <tr>
      <td>while ||δh|| ≥ ε do</td>
      <td></td>
    </tr>
    <tr>
      <td>Calibrate model parameters (from total model) p</td>
      <td>▷ (C)</td>
    </tr>
    <tr>
      <td>Normalize all parameters and topography so that &lt;|p<sub>2</sub>|&gt; = 1</td>
      <td></td>
    </tr>
    <tr>
      <td>Compute residual detector images δf<sup>(i)</sup></td>
      <td></td>
    </tr>
    <tr>
      <td>Compute gradient (and laplacian) from tangent model</td>
      <td>▷ (A2)</td>
    </tr>
    <tr>
      <td>Compute topography correction δh</td>
      <td>▷ (B2)</td>
    </tr>
    <tr>
      <td>Update topography h ← h + δh</td>
      <td></td>
    </tr>
    <tr>
      <td>end while</td>
      <td></td>
    </tr>
    <tr>
      <td>end procedure</td>
      <td></td>
    </tr>
  </tbody>
</table>

2.4.2 | (B) Computation of topography

As earlier mentioned, $\gamma_x$, $\gamma_y$ and $\kappa$ are treated as three independent fields, and hence computing the topography correction $\delta h$, consists in solving the system

$$
\begin{align}
\nabla \delta h(\mathbf{x}) &= \gamma(\mathbf{x}) \tag{12} \\
\nabla^2 \delta h(\mathbf{x}) &= \kappa(\mathbf{x}) \nonumber
\end{align}
$$

The difficulty of this problem is that the expected cross-derivative identity, $\partial_x \partial_y(\delta h) = \partial_y \partial_x(\delta h)$, is generally not satisfied with the predetermined $\gamma$, $\partial_x \gamma_y \neq \partial_y \gamma_x$. The same holds obviously for the laplacian which violates the expected relation $\gamma_{x,x}+\gamma_{y,y} \neq \kappa$.

Different integration strategies can be envisioned. However, the least sensitive to noise can easily be determined by formally considering a white noise being added to $\gamma$ and $\kappa$, computing it effect on $\delta h$ and minimizing its variance$^{24}$. This leads formally to an optimal solution. Writing the problem in Fourier space,

$$
\begin{equation}
\begin{aligned}
i \underline{k} \tilde{\delta} h(\underline{k}) &= \tilde{\gamma}(\underline{k}) \\
-k^2 \tilde{\delta} h(\underline{k}) &= \tilde{\kappa}(\underline{k})
\end{aligned}
\tag{13}
\end{equation}
$$

where the oversign $\tilde{\cdot}$ denotes the Fourier transform, and $\underline{k}$ the wavenumber vector, one sees the linearity of the problem which is solved in the least squares sense (optimal with respect to Gaussian noise) according to

$$
\tilde{\delta h} = - \frac{i \underline{k} \tilde{\gamma} + \omega^2 |\underline{k}|^2 \tilde{\kappa}}{(|\underline{k}|^2 + \omega^2 |\underline{k}|^4)} \quad (14)
$$

where $\omega^2$ is a weight (the ratio of variances of the noise affecting gradient and curvature). $\omega$ can also be seen as the cross-over inverse wavenumber below which $g$ sets the gradient, and above which the curvature dominates.

2.4.3 | (C) Calibration of model parameters

The procedure (c) consists in a recalibration of the entire model for each detector. As written, the model — even when “non-linear” — is linear with respect to the model parameters, $p_0$, $\underline{p}_2$, $p_3$, $q$, $r$, etc (globally denoted as $p$). Thus, assuming that the current topography is given, re-calibrating the detector response is again a (highly redundant) linear system, which again is solved in the least squares sense.

Let us note however that when the model is enriched to include the effect of the large field of view, Eq. 4, then treating model parameters including $\alpha$ altogether would imply a breakdown of linearity because of the coupling between $\alpha$ and $p_2$. However, the problem is linear in $\alpha$ and linear in the set of all $p$. Thus it is chosen to optimize these two groups of parameters sequentially.

2.5 | Pixel weight

Both (A) and (C) procedures allow a more tailored determination of either topography gradient or model parameters. Indeed, in both cases, a least square minimization is performed, assembling different pieces of information at each pixel and each detector.
In the above presentation, all pixels and detectors were considered with the same weight. However, some pixels may have a slope such that they are directly facing one detector and not other ones, yet two independent measurements are sufficient to solve for the slopes. Hence the detector with the worst orientation may be omitted for any pixel, or more generally a weight can be attributed to all pixel/detector pair to lower the influence of the least trust-able piece of information, leading to a mere weighted least-squares solution. This can be seen as a kind of “training”, since the more reliable information is used preferentially, and, in turn, allows to adjust the parameters of the low weight detector model. This feature is particularly at play for the non-linear models, where knowledge issued from the small slopes regions helps designing the appropriate response for unfavorable geometry.

Moreover, it will be seen in the following, that BSE images can easily be saturated for a few pixels $f^{(i)} = 0$ or $2^N - 1$ gray levels, for a $N$-bit deep encoding. This helps benefiting from a wide dynamics in each image, however a gray level recorded as 0 actually means $f^{(i)} \le 0$, and similarly for $2^N - 1$ which means $f^{(i)} \ge 2^N - 1$. Generally, the number of such pixels is low, yet, their meaning should not be trusted. Giving a very low weight (say 1%) to those values is a convenient way to disqualify them whenever their value can be provided by other detectors, still avoiding the risk of too many omissions that would make the problem ill-posed (for a null weight).

Some freedom can be exerted here. Assuming that a surface normal pointing toward the detector may be more reliable than for other orientations, we introduce arbitrarily a reference value $d$ for the normal component of the slope, and a range of confidence, $\sigma_d$ such that the following weight is considered

$$w = \exp\left(-\frac{(\partial_n h - d)^2 + (\partial_l h)^2}{2\sigma_d^2}\right) \qquad (15)$$

## 3 | EXPERIMENTAL TEST CASE

The proposed algorithm is tested on three different samples called; **Pyramid**, **Dome** and **Babel**, see figure 2.

<figure><img src="image_3.png" /><figcaption>FIGURE 2 3D shapes of the three samples discussed in this article obtained from AFM scanning</figcaption></figure>

The first two are part of a micro-topography calibration sample manufactured by M2C<fnref n="1" /> the manufacturing details of which are discussed in ²⁵. The provided sample is a single 6 mm square silicon chip mounted on a 12 mm square aluminum substrate. Four patterns are created on the chip using a combination of FIB deposition and milling of platinum. Three of the four geometries are similar stepped pyramids the forth is a dome shape. For this article only one pyramid and the dome are used, henceforth called **Pyramid** and **Dome**. The stepped pyramid has a square base with side of roughly 12 µm and has 3 steps, the first two about 600 nm in height and the third about 800 nm giving it a total height of about 2 µm. The dome is a spherical cap with a height of approximately 1 µm.

These approximate descriptions are provided to understand the shape. The true shape of the objects is provided by SPM. The manufacturer of the sample has described this SPM method as a highly accurate long range AFM measurement ²⁶. This SPM measurement and the accurate coordinates of the 308 landmarks that surround the shapes are provided together with the sample. These landmarks are annular divots with a diameter of 600 nm and a depth of 120 nm. These landmarks are very useful when

<footnote n="1">¹Microscopy Measurement & Calibration, http://www.m2c-calibration.com</footnote>
orienting or aligning different data sets. However, they have proven challenging in the reconstruction process. This is because their topography is sufficiently steep to partly occlude material points for the detectors. This point will be further discussed below.

The third sample **Babel** is one created in house using an FEI ThermoFisher dual beam microscope Helios Nanolab 660. A similar manufacturing strategy was adopted as for the calibration samples. First a homogeneous 3 µm thick layer of platinum was deposited using ion assisted deposition with a beam current of 430 pA and an acceleration of 30 keV. Next, the desired shape was milled using the FIB with the same beam settings. The base of the tower has a diameter of 5 µm with a total height of 1.5 µm.

The **Babel** shape was designed with a specific purpose: The proposed reconstruction algorithm measures the surface slopes and consequently integrates the slope to obtain the surface topography. In general, high slopes are more challenging for a plethora of reasons. In both the **Pyramid** and the **Dome** sample, there are abrupt changes in slope that form (several) complete loops. Consequently, there is no “high fidelity” path from the base to the summit. Any error on the reconstruction of the steep slopes will accumulate in the integration process. In contrast, the **Babel** shape has a gentle path to the summit while also having steep embankments. The objective is that such a sample can be used, through the parameter model adjustment, to “learn” the detectors response for high slopes, based on a topography that is more reliably computed from small-slope paths. This sample can also highlight the quality and weakness of the integration method.

<figure><img src="image_4.png" /><figcaption>**FIGURE 3** The two orthogonal AFM scans for the **Babel** sample and their difference after registering them to the same coordinate system.</figcaption></figure>

The **Babel** sample has been measured by AFM to provide a topography to compare the proposed method to. The 10×10 µm² areas are scanned with a Dimension Icon AFM from Bruker with a scan definition of 2048×2048 pixels. A sharp Si₃Ni₄ tip with a nominal 2 nm tip radius and 15° and 25° front and back angle respectively is applied in PeakForce Tapping mode. The sample has been scanned twice in two orthogonal scan directions. The two scans provide some insight in the uncertainty of the AFM measurements. Figure 3, shows the two scans and their difference after registration to the same coordinate system. The difference is of the order of 6 nm averaged (RMS) over the entire domain. The largest differences occur when the AFM has to ascend or descend a steep slope.

For the BSE acquisition part, the three geometries are captured with the same SEM with an acceleration of 5 keV, a beam current of 0.4 nA and a dwell time of 10 µs per pixel. For each case, the working distance was 4 mm, the image size was 3072 by 2048 pixels with a horizontal field width of [30, 30, 19.5] µm for the samples respectively (Figure 4). It is important to highlight that this SEM is capable of capturing the four detector signals during a single scan. The consequence is that any scanning artifacts will be equal for all detectors. Some SEMs do not allow this modality and require a new scan to obtain each detector signal and thus may require some image registration to correct any sample drift or scanning artifacts.
# 4 | RESULTS

## 4.1 | Reproduction of BSE images

The proposed method consists in reproducing the acquired image from a model that has as inputs the topography *h* and a set of parameters per detector {*p*}. These model reproductions are then compared to the actual acquisitions to compute a residual. This residual is iteratively reduced up to some remainder. The final residuals are interesting to analyze, they show precisely where the actual acquisitions differ from what the model can account for. However, the models at hand do not attempt to account perfectly for all the ways in how the observed intensity can be influenced, *i.e.*, some residual is expected.

<figure><img src="image_5.png" /></figure>

**FIGURE 4** For the three discussed samples (Pyramid, Dome and Babel), four outputs for detector A; Data $f^{(A)}$, Model $\hat{f}^{(A)}$, Confidence *w* and the Residual $\delta f^{(A)}$. The left three columns are plotted with gray values (GV) ranging from 0 to 255 with 0 colored red and 255 colored green to highlight saturation. The residual is plotted with values relative to the total dynamic range (255).

Figure 4 shows the Data $f^{(A)}$, Model $\hat{f}^{(A)}$, Confidence *w* and the Residual $\delta f^{(A)}$ for detector A for all three samples. All detectors provide very comparable results. The left two columns, where the acquisition and the model are shown, the extreme gray values (i.e., 0 and 255) are highlighted in red and green respectively. SEM images are often saturated because of a human preference for high contrast, but saturated pixels constitute a difficulty for such algorithms. The locations of saturation are mostly present at the extreme slopes, either facing towards the detector or away from the detector. For the Pyramid and Dome samples these features are most prominent inside the deep crevasse of the landmarks. These landmarks are not of major interest. Nevertheless, it will be interesting to see how the proposed method handles areas inside the image where the main assumptions are violated. The confidence images (third column) show that the areas on the far side of the detector have been assigned a lower weight (see section 2.5). Finally, the image residuals show that the most pronounced residuals are in the locations with the highest slopes. Clearly, there is some part of the real electron/sample interaction that is not fully taken into account by the proposed model. Nevertheless, errors on these locations remain localized and do not tend to contaminate the rest of the image.

## 4.2 | Detector response

After reconstruction, both the topography and the acquisitions are known. This allows the gathering of the observed intensity as a function of the surface orientation. Figure 5, shows the detector response, again for detector A, as a function of the surface
gradients $\partial_n h$ and $\partial_l h$, where the first is a tilt towards the detector and the second a tilt in a perpendicular direction. This response is shown for the Model $\hat{f}^{(A)}$ and the Data $f^{(A)}$ both as a color map and by the topography of the plotted surface. The 3D surface allows for a more natural estimation of the linearity of the surfaces. For the **Babel** sample, the Data response is nearly linear. For the **Pyramid** the response is linear with a visible cross-shaped signature while for the **Dome** the response is linear for small surface gradients with a steep non-linearity at a certain slope.

The first conclusion from these response figures is that a linear model should be able to represent the image formation reasonably. All surfaces seem slightly curved, so perhaps a low order polynomial representation does improve the model. However, accurately describing the fine details of these response maps would require an unreasonably high order description.

<figure><img src="image_6.png" /><figcaption>**FIGURE 5** The reflectance maps for the three discussed samples (Pyramid, Dome and Babel), i.e., the detector response as a function of surface orientation for the Data $f^{(A)}$ and the Model $\hat{f}^{(A)}$. The shown surface is representing the same value as the color map.</figcaption></figure>

Figure 6 shows the histogram and residual of the pixels in this surface orientation space. The histogram is showing counts, i.e., number of pixels for that orientation, using a logarithmic color map. About $6 \cdot 10^6$ pixels (about 95%, shown in red) are at a slope close to zero. Consequently, the majority of the parameter identification is weighted for these low angles. This makes it natural to expect the obtained parameter identification to perform better at low angles. This is confirmed by the residual maps. The areas with elevated residuals are typically for slopes that are sparsely populated.
<figure><img src="image_7.png" /><figcaption>**FIGURE 6** The RMS residual relative to the dynamic range and the histogram as a function of surface orientation for the three discussed samples (Pyramid, Dome and Babel), for detector A.</figcaption></figure>

**4.3 | Topography**

From an algorithmic point of view, the topography is a byproduct of minimizing the residuals. Nevertheless, it is the main goal for applying the algorithm. Figure 7 shows the three reconstructed topographies. On the surface, the algorithm seems to have been successful. Perhaps foremost, the algorithms have proven very robust. The landmarks of the **Pyramid** and **Dome** samples are quite challenging. Not only are the slopes steep, they are sufficiently narrow that the detector is partly occluded. This is a non-local phenomena that is not at all accounted for in the proposed method. Consequently, it is expected that the produced topography is of low quality in these areas. Remarkably, the algorithm remains well behaved in these areas, and these regions where the reconstruction cannot be successful do not contaminate their surrounding.

<figure><img src="image_8.png" /><figcaption>**FIGURE 7** The reconstructed surfaces for the three samples obtained with a single calibration parameter, the z-amplitude. Let us emphasize that the **Babel** sample is reconstructed with the BSE images over a wider field of view than the corresponding AFM shown in Figure 2.</figcaption></figure>

Next, the topographies computed from the BSE images are compared to their counterparts measured with an AFM. This allows for a quantitative evaluation of the method. Before such a comparison can be done, the two data sources have to be registered. First, the AFM data is corrected for average height and constant slope. Next, a Digital Image Correlation (DIC) routine is applied
that only allows for affine motion. Allowing the DIC routine to also correct for constant deformations is important here because both data sets rely on scanning methods that often are skewed slightly because of drift.

The proposed method requires one calibration step. The computed topography is scale free, and height has to be scaled to metric units with a reference. The conversion from the measured $h(x)$ to absolute height depends on many acquisition parameters that are often not controlled for, e.g., working distance, beam acceleration, analog to digital conversion, etc. To construct a quantitative topography, a height conversion is required. This factor can be found by having an object of known height in the field of view, or by measuring a specific point ex-situ. In the following, the scale factor is chosen so as to minimize the difference between the two topographies.

<figure><img src="image_9.png" /><figcaption>**FIGURE 8** Difference between topography maps (AFM minus BSE-based) in nanometers for the **Pyramid**, **Dome** and **Babel** samples from left to right</figcaption></figure>

Figure 8 shows the topography difference between the AFM measurements and the BSE reconstructions. Although the AFM measurements have some error, it is assumed that these are sufficiently small such that the major part of the shown error fields are due to errors in BSE reconstruction. A few general remarks can be made. (i) The **Pyramid** shape creates these concentric areas that are separated by high slopes. Consequently, the high uncertainty on these slopes tends to accumulate. This caused the first two steps to be overestimated while the final step underestimated. (ii) Similarly, the summit of the **Babel** shape has the steepest slopes causing the majority of the error to be concentrated there. (iii) All three fields show a smooth parabolic (or long wavelength) error. It is to be noted that this spurious overall curvature has only a minute effect on the surface normals, and the orientation of the latter is the most sensitive characteristic which is captured in the BSE images. Hence, the uncertainty on these long-wavelength components is the highest. In contrast small scale details are much more robustly reconstructed.

The full-field errors shown in Figure 8 are summarized in Table 1. The used BSE detector has four segments. However, only three of them contain directional information. The annular detector D has been proven of little impact. The proposed algorithm accepts the fourth image without any modifications. However, the overall error remains roughly the same. Additionally, the **Pyramid** and **Dome** cases were processed with the landmarks masked. For both the three and four detector variants, this improved the overall error slightly, but not significantly. Finally, the average (RMS) errors range from 14 to 32 nm, which is about 1% to 1.5% of the total height. Depending on the application it can be argued that this range of error is acceptable.

<table><caption>**TABLE 1** Topography rms errors, $(\delta h^2)^{1/2}$. (all corrections included)</caption><thead><tr><th rowspan="2">Sample</th><th rowspan="2">Height</th><th colspan="2">3 BSE images</th><th colspan="2">4 BSE images</th></tr><tr><th>Full image</th><th>Masking landmarks</th><th>Full image</th><th>Masking landmarks</th></tr></thead><tbody><tr><td><b>Pyramid</b></td><td>2.23 µm</td><td>32.2 nm</td><td>30.7 nm</td><td>32.3 nm</td><td>30.4 nm</td></tr><tr><td><b>Dome</b></td><td>1.32 µm</td><td>14.2 nm</td><td>10.3 nm</td><td>14.4 nm</td><td>9.8 nm</td></tr><tr><td><b>Babel</b></td><td>1.55 µm</td><td>19.0 nm</td><td>—</td><td>19.3 nm</td><td>—</td></tr></tbody></table>
5 | CONCLUSION

In this paper, a method is proposed to reconstruct the 3D topography from multi-segment BSE images. The method is constructed from the fundamentals of BSE image formation inside an SEM. From basic principles, a model is independently setup for each detector. The detector specifics such as its sensitivity, shape and location are captured with a few parameters which are automat-ically tuned in the proposed procedure. Both the topography and the detector parameters are obtained in a single optimization process that minimizes the difference between the actual acquisitions and the respective model images for each detector.

The topography is obtained from two orientation images which are integrated to obtain the topography. As a consequence, the number of unknowns is two per pixel, plus a couple of detector parameters. This requires a minimum of three detectors to be used to solve the problem. Consequently, the proposed method is highly general and works with many possible detector configurations. Any number (greater than 2) detector segments can be used which are not required to be symmetric, diagonally opposed, of the same size, etc.

The proposed method requires a single (optional) calibration parameter to convert the produced topography map from mean-ingless arbitrary values to a height in meters. This single parameter can be obtained if some feature in the view is of a known height, or if the sample is measured by other means before or after the experiment. The calibration is nevertheless optional. Often, 3D visualizations of the topography evolution can be highly revealing. For these applications, the non-calibrated topographies can be used directly.

The method was analyzed on three test samples for which AFM data were also obtained. Two of the three are part of a commercially available calibration target while the third was made in house. The novelty of this third “Babel” sample is that it combines high and low slope areas in such a way that a low slope path is preserved to the summit. This causes it to be less sensitive to integration errors and thus a possible improvement to calibration targets for topography in general.

In the comparisons between the BSE reconstructions and the AFM data, the RMS errors for the three samples were [14, 19, 32] nm for samples of height [1.32, 1.55, 2.23] µm respectively. The results in a relative error of about 1.5% for the worst case sample. It is quite legitimate to assume that there are many cases where such an error is acceptable.

Aside from the well established highly accurate family of scanning probe measurements, the proposed approach, purely based on SEM, provides a reasonable (relative) accuracy. The acquisition of BSE images is often easier to combine with in-situ measurements making them much more accessible. Consequently, the proposed method may offer a solution where SPM is challenging to apply or simply not available.

Acknowledgements

This work has benefited from the financial support of the LabeX LaSIPS (ANR-10-LABX-0040-LaSIPS) managed by the French National Research Agency under the "Investissements d'avenir" program (n°ANR-11-IDEX-0003-02). SEM images were acquired from an SEM microscope which is part of the MATMECA platform, supported by the "Investissements d'avenir" program (n°ANR-10-EQPX-37 MATMECA).

The images and code discussed in this work are available on request.

BIBLIOGRAPHY

References

1. Adibi Sara, Branicio Paulo S., Liontas Rachel, et al. Surface roughness imparts tensile ductility to nanoscale metallic glasses. *Extreme Mechanics Letters*. 2015;5:88–95.

2. Chern Andrew H., Nandwana Peeyush, McDaniels Robert, et al. Build orientation, surface roughness, and scan path influ-ence on the microstructure, mechanical properties, and flexural fatigue behavior of Ti–6Al–4V fabricated by electron beam melting. *Materials Science and Engineering A*. 2020;772(November 2019):138740.

3. Shibutani Yoji, Koyama Atsuhiro. Surface roughness effects on the displacement bursts observed in nanoindentation. *Journal of Materials Research*. 2004;19(1):183–188.

4. Giessibl Franz J.. Advances in atomic force microscopy. *Reviews of Modern Physics*. 2003;75(3):949–983.
5. Piazzesi G. Photogrammetry with the scanning electron microscope. *Journal of Physics E: Scientific Instruments*. 1973;6(4):392.

6. Samak D., Fischer A., Rittel D. 3D Reconstruction and visualization of microstructure surfaces from 2D images. *CIRP Annals - Manufacturing Technology*. 2007;56(1):149–152.

7. Woodham Robert J. *Photometric Method For Determining Surface Orientation From Multiple Images*. 1980.

8. Ikeuchi Katsushi, Horn Berthold K.P.. Numerical shape from shading and occluding boundaries. *Artificial Intelligence*. 1981;17(1-3):141–184.

9. Beil W., Carlsen I. C.. A combination of topographical contrast and stereoscopy for the reconstruction of surface topographies in SEM. *Journal of Microscopy*. 1990;157:127–133.

10. Tafti Ahmad P., Kirkpatrick Andrew B., Alavi Zahrasadat, Owen Heather A., Yu Zeyun. Recent advances in 3D SEM surface reconstruction. *Micron*. 2015;78:54–66.

11. Jablonski A.. Analytical applications of elastic electron backscattering from surfaces. *Progress in Surface Science*. 2003;74(1-8):357–374.

12. Paluszyński J., Slówko W.. Surface reconstruction with the photometric method in SEM. *Vacuum*. 2005;78(2-4):533–537.

13. Paluszyński J., Slówko W.. Measurements of the surface microroughness with the scanning electron microscope. *Journal of Microscopy*. 2009;233(1):10–17.

14. Mignot Christophe. Color (and 3D) for Scanning Electron Microscopy. *Microscopy Today*. 2018;26(3):12–17.

15. Paluszyński J., Slówko W.. Compensation of the shadowing error in three-dimensional imaging with a multiple detector scanning electron microscope. *Journal of Microscopy*. 2006;224(1):93–96.

16. Lebiedzik J. An Automatic Topographical Surface Reconstruction in the SEM. *Scanning*. 1979;2:230–237.

17. Carlsen I C. Reconstruction of True Surface-Topographies in Scanning Electron Microscopes Using Backscattered Electrons. *Scanning*. 1985;7:169–177.

18. Pintus Ruggero, Podda Simona, Vanzi Massimo. Image alignment for 3D reconstruction in a SEM. *Microelectronics Reliability*. 2005;45(9-11):1581–1584.

19. Shi Boxin, Wu Zhe, Mo Zhipeng, Duan Dinglong, Yeung Sai-Kit, Tan Ping. A Benchmark Dataset and Evaluation for Non-Lambertian and Uncalibrated Photometric Stereo. In: :3707–3716; 2016.

20. Neggers Jan, Héripré Eva, Bonnet Marc, et al. Principal image decomposition for multi-detector backscatter electron topography reconstruction. *Ultramicroscopy*. 2021;227(December 2020):113200.

21. Darliński A.. Measurements of Angular Distribution of the Backscattered Electrons in the Energy Range of 5 to 30 keV. *physica status solidi (a)*. 1981;63.2:663–668.

22. Pintus Ruggero, Podda Simona, Vanzi Massimo. An automatic alignment procedure for a four-source photometric stereo technique applied to scanning electron microscopy. *IEEE Transactions on Instrumentation and Measurement*. 2008;57(5):989–996.

23. Miyamoto Atsushi, Chen Deshan, Kaneko Shun’ichi. Bootstrapping de-shadowing and self-calibration for scanning electron microscope photometric stereo. *Measurement Science and Technology*. 2014;25(10).

24. Frankot Robert T., Chellappa Rama. Method for Enforcing Integrability in Shape From Shading Algorithms. *IEEE Transactions on Pattern Analysis and Machine Intelligence*. 1988;10(4):439–451.

25. Ritter Martin, Dziomba Thorsten, Kranzmann Axel, Koenders Ludger. A landmark-based 3D calibration strategy for SPM. *Measurement Science and Technology*. 2007;18(2):404–414.
26. Dai Gaoliang, Pohlenz Frank, Danzebrink Hans Ulrich, Xu Min, Hasche Klaus, Wilkening Guenter. Metrological large range scanning probe microscope. *Review of Scientific Instruments*. 2004;75(4):962–970.

<figure><img src="image_10.png" /></figure>