Spatiotemporal Calibration of Electron Microscopes Spatiotemporal Calibration of Electron Microscopes

This chapter presents a method of calibration for scanning electron microscopes (SEMs). The described calibration method is a twofold step. The ﬁrst step evaluates the dynamical drift parameters. The second step estimates the parameters of the geometric projection and the static distortion. Both steps process the calibration parameters across the range of magniﬁcation scales, thanks to a representation with partial differential equations. The chapter is provided with an example of calibration of the JEOL JSM 820 and an example of application in metrology. The presented method is not the unique way of calibrating an SEM and can be a good start to inspire other methods of calibration.


Introduction
SEM calibration is required in many applications: 3D reconstruction of microscale and nanoscale specimens [1][2][3], deformation measurement [4], nanomaterial tracking [3,5,6], positioning and handling [7], mobile robot positioning [8], etc. The image acquisition process of a standard SEM has usually three stages: (i) The conditioning of the electron beam from the electron gun to the area of the observed specimen. In this stage, the electron beam goes through condensation process by magnetic lenses and a final deflection coils that move the beam in a raster fashion over a rectangular area of the specimen. (ii) The interaction between the deflected beam and the observed object, usually called specimen. When the deflected beam hits the specimen, it results in an emission of energy that can be amplified and detected by specific electronic detectors. (iii) The conversion of the emitted energy to pixel intensity. In this process, each position of the deflected beam that hits the specimen corresponds to a pixel in the image. An appropriate transform between the detected energy and the gray-scale values is then applied to obtain the corresponding pixel intensity. In the case of a distortion-free imaging process, the obtained image would be the result of the projection of the 3D rastered area onto the image plane. According to an increasing magnification scale, this projection goes from the perspective for small magnification to the weak perspective for average magnification and ends up with orthographic projection with high magnification. However in reality, stages (i) and (ii) bring some distortion artifacts to these projection transforms. Indeed, the hysteresis effect of the magnetic coils distorts the deflected beam during the observation process. These distortions cause some drifts of the pixels from their original projected positions. For instance, these drifts can exceed tenth of pixels (hundreds of microns) in an hour at 10 k× magnification with a JEOL JSM 820. These drifts can be modeled as a combination of static (spatial) and dynamic (spatiotemporal) mappings which will be composed with the projection mapping to obtain the acquired image.
State-of-the-art studies have established that dynamical drift is time dependent and magnification dependent but space independent [9]. Among estimation methods, the drift was estimated in full field based on digital image correlation (DIC) in [9][10][11][12]. In [13], it was estimated in the frequency domain using FFT. The trajectory flow of the drift over time has been modeled using of B-spline curve fitting [9]. A second-order dynamical system embedded in a Kalman filter was tested and validated to model a thermal drift calibration in scanning probe microscopy [14]. The spatial calibration which involves both static distortion calibration and projection calibration was modeled similarly as with classic optical imaging systems [15,16]. Some state-of-the-art works assumed a non-radial behavior of the static distortion and address this problem using B-splines to fit the spatial evolution of distortion and then warp it to a 3D-to-2D projection. Other works considered a perspective projection in the case of low magnifications (up to 5 k×) and an orthographic projection in the case of high magnifications (more than 5 k×) [1,9,17,18].
In this chapter, we go through an empirical calibration method to fit a system of partial differential equations (PDEs). The partial derivatives of the calibration parameters are estimated with respect to time and magnification. These parameters include the amount of drift in pixel, the amount of static distortion, the focal length, and the principal point. This modeling provides a systematic and flexible solution to this calibration problem. It allows us to smoothly update calibration parameters across the variation of magnification scales and to continuously compensate pixel drifts during experiments.

Notations
Two-dimensional points in homogeneous coordinates are denoted by symbols in typewriter font (e.g., u = (u x , u y , 1) T ). Three-dimensional points are indicated by plain letters (e.g., C = (C x , C y , C z ) T ). Matrices are denoted by uppercase sans serif font (e.g., A). This notation is also adopted for n-dimensional vectors. However, vectors providing a direction in 3D are represented using plain lowercase topped by an arrow (e.g., l). For convenience, and given two 3 × 1 vectors l and m, the dot product is indicated either using < . , . > or using regular matrix/vector multiplication (e.g., < l, m > = l T m), and the cross product is carried either using the symbol × or using the skew symmetric matrix (e.g., l × m = [ l] × m). ||.|| 2 denotes the vector 2-norm in any real vector space R n of finite dimension n. The symbols µm and nm designate, respectively, micrometer and nanometer unit distances. The abbreviation "w.r.t." stands for "with respect to."

Formulation of the imaging model
The imaging model of a SEM is the composition of three mappings [19,20] which are depicted in Figure 1: 1. A pixel-drift mapping from drifted to non-drifted image. It is both magnification and time dependent. To experimentally quantify this drift, we acquire images of a usual specimen (see Figure 2 left) at different times and for different magnifications. Correlation in the frequency domain of successive pairs of images is used to estimate the drift [13]. The dynamic of the drift with respect to time t and magnification g is modeled by two partial differential equations (PDEs) whose coefficients are estimated by means of the principal differential analysis (PDA) approach [21].

2.
A spatial distortion mapping which is magnification dependent. It may concern both radial and tangential distortions [22]. In this chapter, we consider only radial distortions.
3. A 3D-to-2D projection mapping. The type of projection to be applied depends on the magnification. For low magnifications, it can be considered as central perspective projection. But for magnifications above 1000×, it should be interpreted as parallel projective [9]. A specific calibration specimen is used [20]; see Figure 2 right. It contains squares of various sizes, enabling the calibration over a wide range of magnification. Images of this specimen are acquired for various magnifications and poses. An image-centered radial and tangential model that expresses the distorted points w.r.t. the undistorted points is developed. PDEs with respect to magnification of the evolution of the static distortion parameters are established. The estimation of the projection matrix where the magnification factor is embedded is proposed. A bundle-adjustment optimization [23] of the reprojection error between 3D points and their corresponding image pixels allows us to refine the estimated parameters.
Mathematically speaking, the above description can be formalized as: where Q is a 3D point of the observed specimen.q is the corresponding acquired image pixel. T d t,g : R 2 → R 2 and T s g : R 2 → R 2 are two-dimensional mappings (images x-axis and y-axis) which respectively represent the dynamical drift and the static distortion. Π g is the 3D-to-2D projective mapping. To retrieve a 3D ray incident from a 3D point of scanned scene, the corresponding pixel is first corrected from the drift effect; then it is statically undistorted and finally back-projected. In the following, we show how to estimate these mappings. Also, we provide examples of application on a JEOL JSM 820 SEM in secondary electron (SE) imaging mode for time ranging from 0 to 30 min and magnification ranging from 100× to 10 k×.

Drift modeling
The dynamical drift of an image frame at a time t can be represented as a trajectory flow (x(t), y(t)) T of pixels with respect to time: whereÎ t is the intensity of the acquired image which is submitted to drift. I t is the ideal image without drift effects. (x 0 , y 0 ) T is the non-drifted pixel position and (x(t), y(t)) T is the amount of drift. As described previously, SEM images are produced pixel by pixel with 10.5772/61947 a rastering process. An exact identification of the drift should be to follow the amount of drift flow in a pixel-wise manner (x(x 0 , t), y(y 0 , t)) T . In [10,19], a full-field correction of the image drift is proposed. For SEMs where the acquisition time t F of one frame is faster than the pixel-wise dynamic of the drift, it can be assumed as being a global pixel drift between frames especially at low magnification (less than 5000×). This global drift displacement can be easily assessed using image cross-correlation computation between frames. In this case, the first acquired image is considered as the ideal one at t = 0. Taking into account this assumption, equation (2) can be reconsidered aŝ From now on, we drop the hat symbol on the notation of acquired image unless it is worth to mention it. If two successive image frames I t−1 and I t contain the same view at different pixel positions, the cross-correlation integral has a large value at the vector (δx(t), δy(t)) T which corresponds to the drift of the features.
(δx(t), δy(t)) = arg max x,y with where i + x, j + y, x, and y are pixel coordinates running over the image domain. According to [10], at magnification higher than 5000×, the dynamical drift varies even during the scanning process. It showed that the amount of drift variation within an image ranged from 0.3 to 0.9 pixels at 10 k×. In this chapter, the use of an EKF allows us to face this amount of single pixel drift by assuming a Gaussian noise of 1 pixel in the state model [19].
For a given magnification range from g 0 to g f with a step G and a time interval from t 0 to t f with an image acquisition at each sample time T, the image drifts are estimated as follows: for g = g 0 to g f with step G do for t = t 0 to t f with step T do acquire images of the specimen pattern; estimate drift elements δx and δy between frame t 0 and current frame; end end Algorithm 1: Processing of Drift Data.
At the end of this step, two data matrices ∆ x (t, g) and ∆ y (t, g) of rows and columns are obtained. The next step is to use this data drift and PDA to evaluate the dynamics of the drift trajectory functions ∆ x (t, g) and ∆ y (t, g).

Estimating PDEs of the dynamical drift
When describing PDA for linear PDE models, Ramsay and Silverman [24] view the system dynamics as a linear differential operator (LDO) acting upon the process variables. For example, let ∆ x (t, g i ) be the observed data drift which varies w.r.t. the time parameter t at sampled time is assumed to be square integrable. In this chapter, we consider the identification of a second-order LDO which determines the first (speed)-and second (acceleration)-order parameters of the dynamic of the drift [19]: that comes as close as possible to satisfy the homogeneous linear differential equation: In other words, if we estimate the first and second derivatives D∆ x and D 2 ∆ x of ∆ x (t, g i ) w.r.t. time using finite differences, we wish the operator L to annihilate the drift function ∆ x (t, g) as nearly as possible. Thus, we seek a linear differential equation model so that our data satisfies: to the best possible degree of approximation. To carry out PDA, we adopt a least squares approach to the fitting of the differential equation model. The fitting criterion is to minimize, over w x 0 (g i ) w x 1 (g i ) , the sum of squared norms: A is a matrix of t f −t 0 T rows and 2 columns and is always of rank 2. b is a vector of elements. The solution of such an overdetermined least square problem is given as Solving equation (12) for each g i = g 0 , g 0 + G, . . . , g f gives rise to the following weight matrix: W is a matrix of 2 rows and g f −g 0 G columns. The first row represents a discrete sampling of the function w x 0 (g) and the second row a discrete sampling of w x 1 (g). Once again, computing the first-and second-order derivatives of the first and the second row using finite differences on the variation of the magnification gives rise to two differential equations, driving w x 0 and w x 1 variation w.r.t. magnification: The two real valued vectors α 0i α 1i , i = 0, 1 are estimated as in equation (12) after constructing the corresponding matrix A and vector b by using matrix W, the first and second finite differences of each row w.r.t. g. Equation (14) becomes then a PDE of second order that can be easily solved.
Henceforth, the differential equation related to drift function ∆ x can be expressed w.r.t. smooth magnification-dependent weight functions: Finally, equation (15) can be solved as a second-order PDE giving rise to a smooth drift function ∆(t, g). The development conducted to represent ∆ x can be easily followed to represent ∆ y : In order to take into account the noise in the data estimation and the finite differences computation, we choose to embed the differential equations (15) and (16) related to ∆ x and ∆ y in an EKF using a state modeling of the drift function.

Embedding PDEs of the drift in an EKF
If we assume that x(g, t) and y(g, t) are corrupted by a zero mean Gaussian noise ω x (g, t) and ω y (g, t) with covariance Q x (g) and Q y (g), then the stochastic state model can be written as , v x (g, t) and y(g, t) = y(g, t), v y (g, t) are the state vectors with v x and v y the speeds of the drift. The observed variable being the displacement flow x(g, t) and y(g, t), they are assumed to be corrupted by zero mean Gaussian noises γ x (g, t) and γ y (g, t) of covariances R x (g) and R y (g) . This space model representation allows us to write the EKF time update and prediction equations: The same EKF equations can be stated for the drift flow y(g, t).

The multi-scale drift calibration algorithm
In summary, the multi-scale drift flow is characterized in both magnification axis and time axis. Therefore, using the PDA approach, the PDEs w.r.t. time and the PDEs w.r.t. magnification are assessed. Assume a calibration through the range [g 0 , g f ] of magnifications. A calibration pattern with random shapes is positioned upon the stage inside the chamber.
The scene is static, and a set of images are taken at each sample time T. The different steps of the dynamical drift calibration can be summarized in Algorithm 2.

An example of dynamical drift calibration of the JEOL JSM 820
In this paragraph, we show how to use the method developed so far to calibrate the dynamical drift of the JSM 820, a SEM manufactured by JEOL (see Figure 4 for an illustration). The electron gun is equipped with a tungsten filament that can support from 0.3 kV up to 30 kV of acceleration voltage. The acquired images have a size of 512 × 512 pixels. The acceleration voltage is 15 kV, the scan rate is 15 frames per second, and the number of scans average is 8. The calibration is done for magnifications from g 0 = 100× up to g f = 30 k×.

Data drift estimation
To assess the pixel displacement between two frames, a specimen of particles of gold deposed above a layer of carbon is used (Figure 2 left). The particles are randomly positioned and have random shapes of different sizes so that the maximum of cross-correlation images of such a sample can be calculated with less errors [10]. According to the cross-correlation theorem [25], the cross-correlation can be calculated using the Fourier transform. The widely used FFTW3 [26] library is applied for Fourier transform calculations. In order to improve the accuracy of the peak calculation, the cross-correlation is combined with the frequency filtering. A set of 55 images per magnification scale are taken every 30 s. The magnification scale is tuned from 100× until 10 k×. Some of these images are shown in Figure 5 with an illustration of the cross-correlation, resulting peaks between the first frame and the following frames. The pixel displacement vector of all images across time and magnification is depicted in Figure 6. It can be shown that at higher magnification, the drift is more important and can reach up to (20,

PDE estimation of the dynamical drift and EKF embedding
Using PDA, the estimated differential equations associated to w x 0 and w x 1 are Similarly, the weights w y 0 (g) and w y 1 (g) describing the ODE of ∆ y (t, g) are the solution of the two following differential equations: The PDA study of the collected data shows that the time dependence is more likely to be a second-order differential equation and so is for the magnification dependence. This order may not be accurate because of the noise in the data. Thus, to take into account this noise, the two dynamical models of equations (15) and (16) are embedded in an EKF. The Gaussian noise associated to the state model and observation of ∆ x and ∆ y has a zero mean and 0.25 pixels of standard deviation. The plots of the prediction error are shown in Figure 7. Less than 0.9 pixel error in the x-axis and less than 1 pixel error in the y-axis can be observed. Also, 0.28 pixel of global RMS error in the x-axis and 0.23 pixel of RMS error in the y-axis can be observed.

Multi-scale static calibration
The multi-scale static calibration concerns the estimation of the static distortion parameters as well as the projection parameters. They are independent of time, but they are magnification dependent.

Static distortion calibration
In contrast with the dynamical drift, at low magnification, the static distortion is much more significant than at high magnification. This is due to the fact that at high magnification, the scanned area is much smaller than at low ones.

The static distortion model
The most commonly used distortion model represents this physical phenomenon as a decentered distortion which has both a radial and tangential component [27]: where ξ r , ξ t are, respectively, the radial and tangential distortion parameters. e is the center of distortion and λ = 2/1 + 1 − 4 ξ r r 2 u is called the factor of distortion.

Estimation of the center of distortion
The estimation of the center of distortion requires the use of a geometrically structured calibration pattern as the one shown in Figure 8. It consists on a planar grid of vertices {q c i } i∈N . The positions of the points q c i are assumed to be known in a Euclidean coordinate frame attached to the grid. If q d i are the corresponding points in the distorted image, then each pair of points (q c i , q d i ) is linked by the epipolar relation proposed in [28] and stated as [e] x is the skew-symmetric 3 × 3 matrix representing the cross product. H is the homography between the planar grid and the image plane. The matrix F may be called the fundamental matrix for radial distortion. It may be estimated using state-of-the-art methods [15], and the center of radial distortion can be estimated as the left epipole. In the case of no presence of radial distortion, the estimation of the fundamental matrix is unstable, and the value of e is meaningless. This situation can be detected during the estimation of the fundamental matrix.
The estimation of the distortion parameters ξ r and ξ t is processed iteratively. It is first assumed that ξ t = 0 and estimate e. Then ξ t and ξ r are initialized to zero and estimated using bundle adjustment methods [23].

Projection model
The 3D-to-2D projection mapping can vary from perspective to orthographic. State-of-the-art works [18,29] use either a perspective model at low magnification or an orthographic projection for high magnification with projection switch at the magnification of transition which is experimentally determined (usually 5 k×). In this chapter, a magnification-dependent projection model which smoothly switches from a perspective projection to an orthographic projection is detailed [20]. A perspective camera model can be written as where C ∈ R 3 is the position of the projection center. R ∈ SO(3) is the orientation of the projection frame. r i is the i-th row of R, and K is the matrix of intrinsics of the form f , a, and s are, respectively, the focal length, the ratio factor, and the skew parameter.
(p x , p y ) T are the coordinates of the principal point e 0 . The principal ray of the imaging system is in the direction of the vector r 3 , and the value d 0 = − r 3 T C is the distance of the planar grid origin from the camera center in the direction of the principal ray.
Multiplying the magnification by a scale factor g acts similarly as moving the camera center backward along the principal ray. The center of the camera is then moved to C − g r 3 .
Replacing C by C − g r 3 in equation (27) gives the projection matrix at magnification g: where the terms r i T r 3 are zeros for i = 1, 2 because R is a rotation matrix. The scalar d g = − r 3 T C + g is the depth of the world origin with respect to the imaging system center in the direction of the principal ray r 3 . The effect of zooming by a factor g is to move an image point q u on a line radiating from the principal point e 0 to the point q u = gq u + (1 − g)e 0 . From similar triangles, we obtain that The resulting projection matrix at a magnification g is When g → ∞, the projection mapping tends to an orthographic projection.

The static calibration method
Let us consider a calibration at the interval of magnification [g 0 , g t ]. This range is uniformly discretized at a sampling step of δg. The different steps of the static calibration can be outlined by Algorithm 3.
To estimate f (g), a(g), and s(g), we use second-order PDEs and we drive the same reasoning as the computation of drift PDEs in Section 3.4. Since f , a, and s are time independent, only one PDE of second-order is estimated for each parameter.
for g = g 0 to g f with step G do Data: A set of N image with grid points calibrate the projection matrix assuming ξ r (g i ) = 0 and ξ t (g i ) = 0 [15]; comptue f (g i ), s(g i ) and a(g i ). estimate distortion center e(g i ) as explained in Section 4.1.2; if e(g i ) is not degenerate then initialize ξ t and ξ r to zero and estimate them using bundle adjustment method [23]; else set e(g i ) to the principal point and ξ r = 0; initialize ξ t to zero and estimate it using bundle adjustment method [23]; end end estimate f (g), a(g) and s(g) using PDA; Algorithm 3: Static calibration

A small toy application: Cantilever deformation measurement
As a simple application to conclude this chapter, we propose to quantify the deformations of a cantilever that we deform with a micromanipulator. The cantilever is 35 µm long, 3.5 µm wide, and 300 nm thick. The micromanipulator is the Kleindiek MM3A-EM with a rigid tool mounted on the tip. The cantilever is fixed within a holder and is deformed by the moving of the tool tip; see Figure 11. Such an experiment may have several applications in the mechanical characterization of cantilevers [30], of biological deformable objects [31], etc. To check the accuracy of the calibration process, a simple test can be made. We estimate the repeatability of the measures for three different magnifications: 190×, 230×, and 260× [20]. The deformation measures which evaluate our calibration method use the magnification-smooth calibration parameters described in the previous section. It is worth to notice that these parameters were not estimated directly at these three magnification scales, but the estimated magnification-smooth functions allow us to find the calibration parameters at any magnification factor in the range [100× , 10 k×]. A set of 11 configurations are taken for each magnification factor (11 × 3 images). Initially the cantilever is straight and free from any contact with the planar surface. Then, it comes close to the tip of the cantilever which is progressively pushed forward by the MM3A-EM. After 7 × 3 acquired images at different configurations of the deformation, the cantilever is progressively dragged backward to the initial contact-free configuration; see Figure 12. Through the whole experiment, a time tracking frame acquisition is automatically processed using the PC processor's clock trigger. This time acquisition is important to retrieve the amount of pixel by which the acquired images have drifted. The acquired images are undrifted and undistorted with the estimated calibration parameters. We assume that during the deformation, the cantilever sweeps a virtual plane. The affine homography between image pixels and this deformation plane is estimated [16] by taking into account the estimated parameters of the projection model. The Euclidean stratification is done using the length of the cantilever provided by the manufacturer. The stratified homography gives us the mapping between distances in the image and their corresponding metric values in the deformation plane. The amount of deformation is measured as the distance between the tip of the cantilever at rest and its position after deformation. The measured deformation reaches a maximum of 250 nm. After drift and distortion correction, the standard deviation of the error among the three scales is of about 10 nm which is an acceptable amount of error at this scale of magnification.  Deformation measures for the three magnification scales with drift and static distortion correction. The repeatability 10 nm. It can be seen that the drift and distortion correction improve considerably the repeatability of the measurements for the three scales. 10.5772/61947

Conclusion
In this chapter, a spatiotemporal calibration model for SEM imaging systems was presented. This model has the advantage of being smooth with respect to magnification scales. Both dynamic (temporal) and static (spatial) mappings are treated. The evolution over time and magnification of the pixel drifts and of the spatial distortion and projection matrices are modeled by mean of PDA.