Bassel Haydar
email: bassel.haydar@femto-st.fr

Stéphane Chrétien
email: stephane.chretien@univ-lyon2.fr

Adrien Bartoli
email: adrien.bartoli@gmail.com

Brahim Tamadazte
email: brahim.tamadazte@femto-st.fr

3-dimensional OCT Compressed Sensing using the shearlet transform under continuous trajectories sampling

Keywords:

Background : Optical Coherence Tomography (OCT) is an emerging medical imaging technology. It is well suited to various medical applications requiring tissue imaging with micrometer resolution and millimeter penetration depth such as in ophthalmology and dermatology. Despite its numerous advantages, OCT has a long acquisition time for high-resolution images or volumes. This paper deals with the development of a Compressed Sensing (CS) paradigm for faster 3-dimensional OCT image acquisition. Methods : The proposed framework includes three main steps: 1) defining a random-like and parameterizable and continuous scanning trajectories that must be compatible with a smooth mechanical scan, 2) rasterizing the scanning trajectory to make it achievable by a physical system (i.e., galvanometer mirrors), and 3) incorporating a high sparsifying data technique so-called 3D shearlet transform into the compressed sensing scheme. Actually, shearlet transform is mathematically optimal for multidimensional data decomposition and has been proven more efficient than classical ones such as those obtained by wavelet or curvelet transforms. Actually, shearlet system provides a very efficient tool for encoding anisotropic features (such as edges in images) in multivariate problem classes. Results : Numerical simulations and ex vivo experiments were carried out. The obtained results showed the ability of the proposed method to recover OCT images and volumes with high fidelity for different subsampling rates and scanning schemes, demonstrating the relevance of the proposed approach.

Introduction 1.Overview and motivations

Optical Coherence Tomography is a powerful biomedical imaging technology that uses low coherent light to capture micrometer-scale resolution data from within optical scattering media (e.g., biological tissue). Actually, OCT data, sometimes called optical biopsy, can image a tissue at or near the resolution of the well-known histopathology procedure without the need for excisional biopsy. Consequently, the physician can take unlimited biopsies and reduce the need for physical biopsy that are more invasive and imply a delay in the interpretation, namely based on histopathological examination. An OCT imaging device allows acquiring cross-sections in a non-destructive and contact-less manner. Its operating principle is based on the use of low-coherence interferometry employing near infrared light [START_REF] Fujimoto | Optical coherence tomography: an emerging technology for biomedical imaging and optical biopsy[END_REF]. Initially, the OCT technology was developed for the ophthalmology and optometry fields where it can be used to obtain detailed images from within the retina, because of the translucent nature, the minimal scattering, the good light-tissue interaction, and the high-light penetration characteristics of the retina. Recently, it has also begun to be used in interventional cardiology to help diagnose coronary artery disease as well as in dermatology [START_REF] Zysk | Optical coherence tomography: a review of clinical development from bench to bedside[END_REF].

There are two main families of OCT imaging techniques: Fourier-Domain (FD) OCT and Time-Domain (TD)-OCT. The FD-OCT method is associated with high-speed scanning mechanism and better resolution as compared to the TD-OCT. In fact, the scan speed of TD-OCT systems depends on the mechanical cycle time of the moving reference mirror driver, whereas in FD-OCT, the reference mirror is fixed, which assists in sampling multiple points from the tissue structures, simultaneously [START_REF] Wang | Comparison of fourier-domain and time-domain optical coherence tomography in the measurement of thinnest corneal thickness in keratoconus[END_REF]. As results, FD-OCT can provide higher resolution, and frame-rate compared to the TD-OCT Chen, Cense, Pierce, Nassif, Park, Yun, White, Bouma, Tearney and de Boer [2005]. However, despite the superiority of FD-OCT systems over TD-OCT ones, it nevertheless has limitations for effective use in several medical applications. Effectively, FD-OCT offers a high frequency acquisition of 1D (optical cores) and 2D (cross-sections), however the required time for 3D data (OCT volume), their processing, transfer, and storage remains highly problematic. For instance, acquiring a volume of 10×10×3.5mm 3 using a standard FD-OCT system such as the Telesto II from Thorlabs®requires more than 1 minute and 10Go of RAM (Random Access Memory) for preview limiting real-time diagnosis or tissue monitoring. Additionally, in FD-OCT, to discretize and digitize the spectral interferograms at high-resolution and frame-rate, it is necessary to use high-performance array detectors such as CCD (Charge-Coupled Device) sensors [START_REF] Choma | Sensitivity advantage of swept source and fourier domain optical coherence tomography[END_REF]. The use of these kind of sensors increases significantly the cost of an FD-OCT system and then limits their deployment in hospitals, particularly in developing countries.

To overcome certain limitations of using OCT imaging modality, namely when the use of C-scan acquisition is required, we explore the potential of developing of Compressed Sensing (CS) scheme in aim to:

• increase the frame-rate of OCT volume acquisition, namely for applications in assisted surgical guidance and intervention

• reduce the amount of data required and subsequent processing for high-resolution image reconstruction [START_REF] Liu | Compressive sd-oct: the application of compressed sensing in spectral domain optical coherence tomography[END_REF];

• make possible to use smaller and less expensive array detectors, while guaranteeing spatial and axial resolutions that are faithful to those obtained by the expensive sensors usually used;

• reduce artifacts in OCT volumes mainly due to the low acquisition frame-time and possible physiological movements of the patient when it comes to in vivo examination.

Compressed Sensing was developed simultaneously by Donoho and Candès, Romberg and Tao Donoho [2006]; [START_REF] Candès | Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information[END_REF] who introduced randomized acquisition schemes and provided the strong mathematical underpinnings of CS theory. It consists of a paradigm which illustrates the possibility to acquire at sub-Nyquist rate and compress signals (measurements) all at the same time under sparsity assumptions, at the expense of potentially more involved computations for the recovery [START_REF] Foucart | A mathematical introduction to compressive sensing[END_REF]. Sparsity, which is the opposite of density, is an essential feature shared by many wavelet-type decompositions, that is leveraged in CS theory [START_REF] Davenport | Introduction to compressed sensing. chapter 1 of compressed sensing: Theory and applications[END_REF]. Most medical imaging systems, such as magnetic resonance imaging (MRI), computed-tomography (CT), OCT, ... are known to admit a high-level of sparsity in other decomposition bases such as wavelets, curvelets, shearlets, etc. This especially because that medical images are generally composed of large homogeneous areas with a low dynamic range in terms of contrast compared to conventional images [START_REF] Fei | Medical image fusion based on feature extraction and sparse representation[END_REF]. In addition, it is admitted to consider that the different medical imaging modalities (e.g., OCT) share the common feature of acquiring, totally or at least partially, the signals through spectral measurements of highly sparse signals. However, in many practical settings, natural sparsity is sometimes insufficient for accurate signal recovery and other tools have to be put to work in order to achieve better reconstruction in inverse problems of the CS type. One of the necessary conditions for the feasibility of compressed acquisition is the ability to decompose efficiency the signal/image in a sparsifying functions bases such as wavelet, curvelet, shearlet transforms, etc. The latter have recently been introduced as a new decomposition that is superior, in many respects, to wavelet multiple scale decomposition [START_REF] Mallat | A theory for multiresolution signal decomposition: the wavelet representation[END_REF]. Shearlet system Labate, Lim, [START_REF] Labate | Sparse multidimensional representation using shearlets[END_REF] provides an efficient compactly supported decomposition which extends the wavelet decomposition and was proved to outperform the wavelet decomposition in many applications such as sparsifying method. A notable feature of shearlet is that it can encode curvilinear singularities and others anisotropic features in a much better fashion than wavelet functions can do. As a result, shearlet provides a more precise separation of the morphologically distinct features of points and curves. Whereas wavelet systems are obtained using dyadic scaling and translations of a single function Ψ, shearlet systems can be generated using a single operator including parabolic scaling, shearing, and translation operators [START_REF] Kutyniok | Introduction to shearlets[END_REF]. From the point-of-view of CS, shearlet systems provide a high-sparsity representation of the signal, which is practically often optimal, a considerable advantage over other techniques such as Fourier and wavelet transforms.

Contributions

Most of the existing CS methods requires the use of incoherent pseudo-random or random subsampling of -sparse data [START_REF] Candes | Stable signal recovery from incomplete and inaccurate measurements[END_REF]; [START_REF] Candes | The restricted isometry property and its implications for compressed sensing[END_REF]. These subsampling schemes (also called masks or sensing matrices) are often of little relevance since they cannot be implemented on physical acquisition systems [START_REF] Bigot | An analysis of block sampling strategies in compressed sensing[END_REF]. In OCT context, the acquisition is processed through scanning scheme performed by a twodegreesof-freedom (dof) galvanometer mirror device and random sampling which leads to the fact that is inconsistent with the physical constraints that preclude the use of non-smooth trajectories. Consequently, in case of a physical implementation, acquiring few measurements in a random sampling method can result in a longer computation time in comparison to the acquisition of all data in raster mode.

In this paper, we investigate novel scanning schemes using continuous trajectories (e.g., spiral, rosette, and Lissajous), that consider the limitations imposed by the OCT mirrors in terms of scanning speed, scan area, and kinematics. We also develop a rasterization strategy which ensures that pixel coordinates are optimally chosen to follow the predefined trajectory (i.e., positions) on the sample to be scanned. Our implementation bridges our new scanning schemes with the existing digital shearlet transform toolbox provided by [START_REF] Kutyniok | Shearlab 3d: Faithful digital shearlet transforms based on compactly supported shearlets[END_REF]. Our scanning schemes are based on spiral, rosette or Lissajous shaped trajectories that offer rapid and efficient coverage of the -space. The proposed rasterization algorithm is embedded into our package and is designed to permit handy tuning of the subsampling rate, the scan area, and the variable density sampling along the trajectory. Finally, the developed methodology is shown to outperform current state of the art in the two following scenarios: 1) using a simulation based numerical validation framework in both 2D cross-sectional OCT image and 3D OCT volumes cases, and 2) using an experimental set-up equipped with a Fourier-Domain OCT device.

The rest of this paper is organized as follows. In Section 2, we present the compressed sensing framework, the sparsity representation as well as the shearlet transform. Also, the acquisition schemes, the proposed subsampling continuous trajectories, and the rasterization algorithm, both the sparsity and the incoherence are detailed in this section. Section 3 deals with the numerical validation when Section 4 discusses the experimental validations of the proposed methods and materials. Finally, a discussion on the obtained results with regard to the state-of-the-art is also provided.

Materials and Methods

Compressed sensing

Compressed sensing was discovered in the breakthrough papers by E. Candés, J. Romberg, T. Tao and D. [START_REF] Donoho | Compressed sensing[END_REF]; [START_REF] Candès | Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information[END_REF] that created via a mathematical tour de force, a new paradigm for joint signal acquisition and compression and explored its tight relationships with modern optimization and random matrix theory. It has subsequently triggered an extensive research effort exploiting sparse representation of signals and images via fast orthogonal and even non-orthogonal decomposition such as wavelet transforms. CS initially attracted the signal processing community's attention after Candès and Donoho were able to show how to use random sampling in order to break the Nyquist barrier that had set the what used to be thought of as the sampling frequency limit to any data acquisition procedure. Then, the theory developed into richly ramified research field [START_REF] Foucart | A mathematical introduction to compressive sensing[END_REF] with a versatile set of mathematical and algorithmic tools for efficient sampling of inherently sparse objects.

In mathematical terms, the problem can be stated as follows. Let be an object (a vector, a matrix or a tensor) in a Euclidean space which admits a -sparse representation in a dictionary Γ, i.e.,

= ∑ =1 Γ (1)
where ∈ ℝ is a -sparse vector, which means that a vector with no more than nonzero components. The observations are simply given by linear measurements of the form

= ⟨ , ⟩
where ⟨ , ⋅⟩ is a functional on the for = 1, … , , and we obtain a linear system

= (2)
where ∈ ℝ × and the rows of are given by

,∶ = Γ ; = 1, … , .
In a CS type of problem, our goal is to recover the sparse vector of components with as few observations (i.e., measures) as possible and therefore will be thought as small compared to . The main challenge is then to construct an observation matrix , which allow recovering with as small as possible for given values of and .

The CS problem can be solved unambiguously if there is no sparser solution to the linear system (2) than . Thereby, recovery is obtained by simply finding the sparsest solution to (2). If for any in ℝ we denote by ‖ ‖ 0 the 0 -norm of , i.e. the cardinal of the set of indices of nonzero components of , the CS problem is equivalent to

min ∈ℝ ‖ ‖ 0 s.t.
= .

(3)

Let us denote by Δ 0 ( ), the solution of problem (3) and Δ 0 ( ) is called a decoder1 . Thus, the CS problem may be viewed as a combinatorial optimization problem. Moreover, the following lemma is well known. [START_REF] Cohen | Compressed sensing and best -term approximation[END_REF]) If is any × matrix and 2 ≤ , then the following properties are equivalent:

Lemma 1. (see for instance

• The decoder Δ 0 satisfies Δ 0 ( ) = , for all ∈ Σ ,
• For any set of indices with # = 2 , the matrix has rank 2 where stands for the submatrix of composed of the columns indexed by only.

The 1 relaxation is given by min ∈ℝ ‖ ‖ 1 s.t.

= .

(4)

In the following, we will denote by Δ 1 ( ) the solution of the 1 -relaxation (4). From the computational viewpoint, this relaxation is of great interest since it can be solved in polynomial time. Indeed, (4) is equivalent to the linear program

min ∈ℝ ∑ =1 s.t. -≤ ≤ , and = .
The main subsequent problem induced by this choice of relaxation is to obtain easy-to-verify sufficient conditions on for the relaxation to be exact, i.e. to produce the sparsest solution to the underdetermined system (2). An algebraic condition was given by Candès, Romberg and Tao [START_REF] Candès | Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information[END_REF], called Restricted Isometry Property (RIP). We say that a matrix satisfies the , if, for every index subset with cardinality | | ≤ , and every ∈ ℝ

(1 -)‖ ‖ 2 ≤ ‖ ‖ 2 ≤ (1 + )‖ ‖ 2 . ( 5
)
The smallest value of in the previous inequalities is denoted by . Up to now, this condition could only be proved to hold with great probability in the case where is a sub Gaussian random matrix. Several algorithmic approaches have also been recently proposed in order to guaranty exactness of the 1 relaxation such as in [START_REF] Juditsky | On verifiable sufficient conditions for sparse signal recovery via 1 minimization[END_REF] and d [START_REF] Aspremont | Testing the nullspace property using semidefinite programming[END_REF]. A different approach to the study of decoder Δ 1 is based on the notion of incoherence. The coherence ( ) of the observation matrix is the largest absolute value among all scalar products of different columns of , i.e.

( ) = max ≠ ′ |⟨ , ′ ⟩| ‖ ‖ 2 ‖ ′ ‖ 2 . ( 6
)
The matrix is said to be incoherent when its coherence ( ) is on the order of 1∕ log( ). Several results have been obtained that guarantee exact recovery under the low coherence assumption.

In the noisy setting, the observation is given by

= ⟨ , ⟩ + (7)
and the problem can be addressed by solving

min ∈ℝ 1 2 ‖ -‖ 2 2 + ‖ ‖ 1 (8)
for specific values of ∈ ℝ + . Exact recovery cannot hold in this setting, but using incoherence, exact recovery of the support of was proved in [START_REF] Candès | Near-ideal model selection by 1 minimization[END_REF] in the case where the variance of the noise is known beforehand, and in [START_REF] Chrétien | Sparse recovery with unknown variance: a lasso-type approach[END_REF]

≥ • log + 1 . ( 9
)
Then the probability that the estimator ̂ defined by (8) with satisfying2

2 = var ‖ ‖ ‖ -̂ ‖ ‖ ‖ 2 2 log , ( 10
)
where a relevant range for var is given by var ∈

(1 -) 2 20(1 + ) spar ‖ ‖ 2 ;

(1 -

) 2 2(1 + ) spar ‖ ‖ 2 , ( 11
)
exactly recovers the support and sign pattern of is greater than 1 -228∕ .

Another approach for fast reconstruction makes use of the Iterative Hard Thresholding algorithm [START_REF] Blumensath | Iterative hard thresholding for compressed sensing[END_REF]. This algorithm is defined by the iterations

( +1) = ( ) + ( -( ) ) (12)
where is the threshold operator which sets all components to zero except for the largest among them. Furthermore, the main result from [START_REF] Blumensath | Iterative hard thresholding for compressed sensing[END_REF] is the following theorem.

Theorem 2.2. Let ♯ denotes the best -term approximation of . By assuming that satisfies the RIP property with 3 < 1∕ √ 32, then, for all , we have

‖ ( +1) -‖ 2 ≤ 2 - ‖ ♯ ‖ 2 + ̃ ( 13
)
where

̃ = ‖ -♯ ‖ 2 + 1 √ ‖ -♯ ‖ 1 + ‖ ‖ 2 . ( 14
)
The Iterative Hard Thresholding algorithm can also be studied from the point of view of incoherence. A breakthrough result of [START_REF] Maleki | Coherence analysis of iterative thresholding algorithms[END_REF] is the following Theorem. Finally, choosing between different estimators is always a matter of finding the good balance between several technological constraints. In the rest of this paper, we will use the Iterative Hard Thresholding method for estimating the sparse vector . The reason for this choice is that Iterative Hard Thresholding is lighter and faster than solving the LASSO problem (8).

Sparse representations

Sparsity is essential to the CS approach to reconstruction with few samples. In order to ensure the sparsity condition, it is necessary to represent the image in a suitable basis such as wavelet, curvelet or shearlet systems. In this work, we opted for the use of the shearlet transform, denoted by , for the sparse representation of the OCT images, because shearlet transform was observed to achieve the best sparsity among these three options. To demonstrate this, we applied successive decompositions of the "Lena" photography having the resolution of 512×512 pixels, using wavelet, curvelet, shearlet transforms, respectively. As can be highlighted in Fig. 1, the shearlet method outperforms both the wavelet and the curvelets ones as summarized in the following:

• wavelet transform: ≈ 15×10 4 null-coefficients;

• curvelet transform: ≈ 5×10 5 null-coefficients;

• shearlet transform: ≈ 2×10 6 null-coefficients.

In a recent work, we also demonstrated experimentally that shearlet decomposition outperforms both the wavelet and the curvelet Duflot, Krupa, Tamadazte and Andreff [2016].

Shearlet system

Shearlet theory provides an efficient mathematical tool for sparse image representation including geometry and multiscale analysis. It is considered as an extension of the wavelet transform, achieved by increasing their directional sensitivity in order to be more adapted for anisotropic image objects (e.g., edges and key points). Indeed the shearlet coefficients can be obtained by applying three successive anisotropic operations: dilation E and Ẽ, shearing S and translation (shift) g on a finite number of generating functions , ̃ ∈ 2 (ℝ 2 ). The anisotropic operators are defined using the following matrices.

The anisotropic dilation and shearing matrices are defined by: • translation: g ℎ(p) = ℎ(pg), (p, g) ∈ (ℝ 2 ) 2 , where g is translation vector, ℎ is a function defined in 2 (ℝ 2 ), and p represents the image-point coordinates.

• dilation: E = diag(2, √ 2) and Ẽ = diag( √ 2, 2) • shearing : S = 1 0 1 , ∈ ℤ
Note that, the operator S is used instead of the rotation parameter, namely used in a curvelet transform [START_REF] Candes | Curvelets: A surprisingly effective nonadaptive representation for objects with edges[END_REF], which can be considered as a significant advantage for discretization, as the integer lattice is invariant under the shear operator for any ∈ ℤ. In other words, similarly to wavelet support, it unified the treatment of the continuous and digital shearlet theory, which leads to a fast implementation of the shearlet system. These operators are used in the following definition of a discrete shearlet system Kittipoom, [START_REF] Kittipoom | Construction of compactly supported shearlet frames[END_REF].

Note that for a 2-dimensional scaling function ∈ 2 (ℝ 2 ), and generating shearlet , ̃ ∈ 2 (ℝ 2 ), a so-called cone-adapted shearlet system can be defined in order to cover the whole frequency plane S. and G. [2014]. Therefore, it is important to choose the generating shearlet functions and ̃ such that they are supported in different frequency cones as depicted in Fig. 2.

In the following, we will describe the mathematical background of the cone-adapted shearlet system that inspired our work. Note that in contrast to the standard discrete transform, it has a direction bias which allows separating the low-frequency region with a square centered around the origin [START_REF] Kutyniok | Introduction to shearlets[END_REF] as can be seen in Fig. 2.

By considering the scaling function , the generating functions , ̃ ) and the sampling constants b = ( 1 , 2 ) ∈ ℝ + × ℝ + , the shearlet system is defined as

( , , ̃ , b) = Φ( , 1 ) ∪ Ψ( , b) ∪ Ψ( ̃ , b), (15)
where

Φ( , 1 ) = { 1 g , g ∈ ℤ 2 } (16) Ψ( , b) = { , ,g ∶ ≥ 0, | | ≤ ⌈2 2 ⌉, g ∈ ℤ 2 } (17) Ψ( ̃ , b) = { ̃ , ,g ∶ ≥ 0, | | ≤ ⌈2 2 ⌉, g ∈ ℤ 2 } (18)
and

, ,g = E S Bg (19) ̃ , ,g = Ẽ S Bg ̃ (20) with B = diag( 1 , 2 ), B = diag( 2 , 1 )
and ⌈⋅⌉ denoting the ceiling function. Also, | | ≤ ⌈2 2 ⌉ is an important condition which allows varying the orientation of the shearlet support up to 4 .

It can be observed that the low frequency region is associated to Φ( , 1 ), when both the horizontal and the vertical cones correspond to Ψ( , b) with Ψ( ̃ , b) (Fig. 2).

Furthermore, the straightforward generalization of 2D shearlet is the 3D decomposition framework defined below. First, let us consider the 3-dimensional scaling function ∈ 2 (ℝ 3 ), generating shearlet , ̃ , ̂ ∈ 2 (ℝ 3 ), and the sampling constants b = ( 1 , 2 , 2 ) ∈ ℝ + × ℝ + × ℝ + for which a so-called pyramid-adapted discrete shearlet system is defined as

( , , ̃ , ̂ , b) = Φ( , 1 ) ∪ Ψ( , b) ∪ Ψ( ̃ , b) ∪ Ψ( ̂ , b)
where the sets

Φ( , 1 ) = { 1 g , g ∈ ℤ 3 } (21) Ψ( , b) = { , ,g ∶ ≥ 0, | 1 |, | 2 | ≤ 2 2 , g ∈ ℤ 3 } (22) Ψ( ̃ , b) = { ̃ , ,g ∶ ≥ 0, | 1 |, | 2 | ≤ 2 2 , g ∈ ℤ 3 } (23) Ψ( ̂ , b) = { ̂ , ,g ∶ ≥ 0, | 1 |, | 2 | ≤ 2 2 , g ∈ ℤ 3 } (24)
and

, ,g = E S Bg ̃ , ,g = Ẽ S Bg ̃ ̂ , ,g = Ê Ŝ Bg ̂ (25) with B = diag( 1 , 2 , 2 ), B = diag( 2 , 1 , 2 )
, so the scaling matrices are obtained as follows

E = diag(2, √ 2, √ 2), Ẽ = diag( √ 2, 2, √ 2), Ê = diag( √ 2, √ 2, 2), (26)
and shearing matrices

S = ⎛ ⎜ ⎜ ⎝ 1 1 2 0 1 0 0 0 1 ⎞ ⎟ ⎟ ⎠ , S = ⎛ ⎜ ⎜ ⎝ 1 0 0 1 1 2 0 0 1 ⎞ ⎟ ⎟ ⎠ , Ŝ = ⎛ ⎜ ⎜ ⎝ 1 0 0 0 1 0 1 2 1 ⎞ ⎟ ⎟ ⎠ , 1 , 2 ∈ ℤ (27)
As results, the shearlet transform allows decomposing an OCT image (respectively, a volume) in a sparse basis as coefficients.

Scanning trajectories

As mentioned in the introduction, usually, the literature related to the compressed sensing methods in imaging, almost all have one characteristic in common: the use of pseudo-random or random subsampling scheme to select the -sparse data that will reconstructed [START_REF] Rauhut | Compressive sensing and structured random matrices[END_REF]; [START_REF] Candes | A probabilistic and ripless theory of compressed sensing[END_REF]. This methodology is suitable for image sensors with an electronic data acquisition process such as CCD/CMOS cameras. Unfortunately, when it concerns imaging systems equipped with a mechanical device (scanning mirrors, magnetic coils, etc.), random scanning methods may be irrelevant owing to acquisition constraints.

To overcome this limitation, the notion of continuous trajectories with variable density samplers is discussed in this section. Continuous and smooth trajectories-based samplers are crucial to extend CS results to a physical implementation, especially for 3D OCT data acquisition and processing. The proposed trajectories are based on spiral, rosette and Lissajous scanning schemes which meet certain criteria such as:

• the design scheme should be reconstruction-guidance-based not observation-based;

• the trajectory must be continuous and easily achievable by the scanning system, i.e., considering the kinematics aspects;

• the scanning trajectory must be tunable in terms of sampling rate, length, sampling step, executing time, etc.

In other words, this results in creating proper 2D and 3D subsampling masks (i.e., sensing matrices) M for 2D or 3D OCT images, respectively.

Parametrizable scanning curves

As mentioned above, in this work, we designed three different continuous scanning trajectories (spiral, rosette, and Lissajous) (Fig. 3) deemed relevant to ensure an effective CS system while respecting the criteria cited.

There are different ways to implement the chosen scanning continuous curves. Easy tuning of the sampling parameters (sampling rate, step time, length, etc.), can be achieved using the following expressions:

• Spiral:

( ) = cos( ) y( ) = sin( ) (28)
with ∈ ℕ allows tuning the distance between two successive spires.

• Rosette:

( ) = sin( 2 +2 ) cos( ) ( ) = max cos( 2 +2 ) sin( ) (29)
where is the curvilinear abscissa, max is the curvature, is the number of rosette's petals,

• Lissajous:

( ) = 1 sin( 2 ( -1) ) y( ) = 2 sin(2 ) ( 30
)
where the parameters 1 and 2 represent the covered surface size, determines the base frequency of Lissajous curve.

Rasterization

To be able to successfully apply the designed sensing matrices (i.e., designed masks) and to match the coordinates of the selected measurements to the corresponding ones in the physical OCT device, we implemented a rasterization technique. The latter allows converting continues geometric curves of the vector format of pixels, points, lines, ... into an image and vice-versa. Among the existing rasterization algorithms, one can cite "Digital Differential Analyzer" (DDA). To meet perfectly the requirements of the proposed methods, we revisited DDA method in order to:

• meet the drawing speed (on-line rasterization) of complex scanning curves;

• convert the drawn trajectory defined in the image (in pixels) to the physical scanning curve (in metric) to be achieved into the sample by the galvanometer mirrors. Recall that the designed sensing matrices M is actually a binary matrix:

M( , ) = 1, acquire an A-scan at coordinates ( , ) 0, do nothing (31)
At the same time, this mask has to be generated from a continuous trajectory. Therefore, it is important to project parametric curve in the discrete basis of mask matrix. Also, the coordinates of M are used in the rasterization algorithm to map the pixels positions (in the image frame) to the metric locations (in the OCT frame) (Fig. 4). Moreover, with the proposed rasterization algorithm it is possible to control the global sampling rate, i.e. the percentage of the measures that will be acquired:

* = 100 × ∑ ℎ =1 ∑ =1 M( , ) Area of the trajectory (32)

Coherence evaluation

As discussed in Section 2, the notion of incoherence ( ) (6) of the observation matrix of the proposed continuous trajectories sampling (i.e., spiral, rosette, and Lissajous) was studied and compared to that of the traditional random sampling matrix. As can be seen in Fig. 5, the numerical values ( ) are very close, with a slight difference for the spiral scanning trajectory, this shows that the proposed sensing schemes enjoy sufficiently small coherence and can therefore be considered as relevant for compressed sensing-based recovery.

Results and Discussion

Numerical validation using realistic biological data

The proposed methods and materials were firstly validated numerically (in simulation) using ground-truth data. The first validation tests consist of using 2D OCT images of 512×512 pixels, while the second uses OCT volumes of 281×281×199 pixels directly acquired by the Telesto II system.

Evaluation criteria

In order to quantitatively assess the results obtained with different scenarios and OCT images, we implemented two criteria: 1) the peak signal noise ratio (PSNR), and 2) the structural similarity index (SSIM), which are generally used by the compressed sensing and more widely the image processing communities [START_REF] Miao | Quantitative image quality evaluation of mr images using perceptual difference models[END_REF]. The PSNR score, given in , is expressed by

= 10 log 10 2 (33)
where is the maximal pixel value in the initial OCT image and the is obtained by

= 1 ℎ ℎ ∑ =1 ∑ =1 I( , ) -R( , ) 2 (34)
with I representing an initial full-scanned 2D OCT slice (selected from the OCT volume), R is the reconstructed one. Furthermore, the SSIM score is based on the computation of three values, namely, the brightness , the contrast and the structural aspect . It is given by

= (R, I) (R, I) (R, I) (35)
where,

R, I = 2 R,I + 3 R I + 3 , R, I = 2 R I + 1 2 R + 2 I + 1
, and

R, I = 2 R I + 2 2 R + 2 I + 2 (36)
with R , I , R , I , and R,I are the local means, standard deviations, and cross-covariance for images R, I. The variables 1 , 2 and 3 are small numbers used to stabilize the division with weak denominator.

Numerical validation: 2D images

The numerical validation was performed following several scenarios. First of all, as inputs in the proposed algorithm we used three OCT volumes (only 2D slice images are depicted for a better visualization) shown in Fig. 6. Note that the first image (Fig. 6(a)) is the "Shepp-Logan" phantom, which is generally used in the literature dealing with medical image processing and CS. In addition, two OCT images of biological samples (i.e., part of grape and retina of a fish eye) were also used in this validation.

First, for each image test, we applied the proposed CS method using various subsampling rate ranging from 10% to 70% with a step of 10% defined using the constructed masks based the developed continuous trajectories introduced in Section 2.4. The obtained results are presented and discussed in the following.

Shepp-Logan phantom (30% of measurements)

The first numerical validation consisted of the reconstruction of the Shepp-Logan phantom image using 30% of measurements. The latter were obtained using the continuous subsampling trajectories. In Fig. 7, are compared the original (ground-truth) image (Fig. 7(a)) and the reconstructed one (Fig. 7(b)). As can be noticed, the recovered image (using the spiral subsampling method) is faithful to the ground-truth one. This is confirmed by the zoom-in thumbnail image. The second remark is that in the reconstructed image, the edges appear smoother than to those of the original image. Finally, the obtained results are very similar using the other two subsampling methods.

The qualitative study of the obtained results are discussed in Section 3.2 using the similarity scores, i.e., the PNSR and the SSIM.

Grape (20% of measurements)

Similarly, the CS method was validated using OCT images of a part of a grape. In this test, the subsampling rate is of 20%. The reconstructed images are depicted in Fig. 8: from left to right, using the spiral, the rosette, and the Lissajous trajectories subsampling methods, respectively. Again, one can highlighted the recovered OCT images are very similar to ground-truth ones. Note that, even the recovered images still similar from one subsampling method to another, Lissajous-based sensing matrix appears more interesting because it covers the entire OCT field-of-vision compared to both the spiral and rosette methods.

Numerical validation: volumes

The proposed methods were also validated using directly the OCT volume instead of individual 2D slices (B-Scans) as input of the proposed algorithm. One can compare the reconstructed OCT volume with the ground-trust one as can be seen in Fig. 9. It can be concluded that the recovered volume (using 30% of measurements) is very similar to the original one.

Quantitative analysis

To assess the performances of the proposed compressed sensing algorithm, we conducted numerous validation tests using both the different types of masks (i.e., spiral, rosette, and Lissajous) and various subsampling rates: 10%, 20%, 30%, 50%, and 70%. For each test, we computed the values of both the PSNR and SSIM (introduced at the beginning of this Section). Then, Table 1 summarizes the obtained numerical values of each test using both the retina of a fish eye and "Shepp-Logan" phantom images as inputs on the CS algorithm. Note that, generally, if the PSNR values are typical between 30 db and 50 dB (for 8-bits encoded images), then reconstructed image is considered faithful to the original one. Also, when the PSNR is equal or greater to 40 db, as consequence the quality of the reconstruction is qualified as similar to the original image.

As can be underlined the obtained PSNR values vary from approximately 40 for 10% of data to 45 db for 70% of data. The first remark is that even for only 10% of samples, the quality of the reconstructed image/volume is very interesting. This is confirmed by the second similarity score (i.e., SSIM) with numerical values varying from 0.97 to 0.99 for 10% to 70% of samples, respectively. The SIMM values are also interesting even only 10% of measurements are used during the recovery task. Remember that, generally, SSIM numerical values are in the range of [-1, 1], where value 1 (respectively, -1) is reachable when the similarity, between the compared images, is "perfect". Furthermore, the rosette and spiral sampling patterns only acquire 4 region of -space due to their circularspace support, as compared to Lissajous method. This factor is considered in the computation of the subsampling rate for each trajectory. The same performances assessment was achieved using the OCT image acquired on the fish eye. The PSNR and SSIM numerical values are reported in Table II in which one can highlighted that the results are very similar to those obtained for the "Sheep-Logan" phantom image. Again, the spiral-scanning trajectory gives better results comparing to rosette and Lissajous.

Experimental Validation

Experimental setup

The proposed materials and methods were validated in both simulation and experimentally. To carry out the experimental implementation of the CS algorithm, we designed a robotic setup which acts as a positioning platform combined with an OCT device, the Telesto II from Thorlabs® (Fig. 10(a)). The latter consists of a FD-OCT system which offer three types of acquisition modes: A-scan (1D optical core), B-scan (2D cross-sectional image) and C-scan ( × × ) volume). In fact, the Telesto II system is based on the principle of low-coherence interferometry because the velocity of light is extremely high. The low-coherence i.e., high-bandwidth light beam is directed to the target tissue and the scattered back-reflected light is combined with a second beam (reference beam), which was split off from the original light beam. The resulting interference patterns are used to reconstruct the A-scan. The latter represents the axial resolution in the tissue, of about 1-15 µm [START_REF] Drexler | State-of-the-art retinal optical coherence tomography[END_REF]. In addition, the light beam is able to move along the tissue in directions. This results in a compilation of A-scans with each of them having a different incidence point. From these series of A-scans, a two-dimensional cross-sectional image of the target tissue can be reconstructed which known as a B-scan. Moreover, a C-scan acquisition consists of the concatenation of parallel B-scans as summarized in Fig. 10(b-c).

Validation using a physical system

The experimental validation scenario is performed as follows: instead of using prerecorded images or volumes with further off-line subsampling them (using sensing matrices) to create the sparse data, we directly (on-line) controlled the galvanometer mirrors of the OCT probe to acquire sequentially only a limited number of measurements (A-Scans) following a predefined continuous trajectory. The coordinates (pixels at 1 in the sensing matrix) of rasterized trajectory are used as the inputs control of the galvanometer mirrors. As claimed previously, it is possible to tune the scanning trajectories parameters in aim to manage the acquisition rate (expressed as a percentage of the entire volume).

The first experiment was conducted by acquiring 30% of samples on the fish eye. The latter was placed under the OCT probe. Therefore, instead of acquiring the entire data, we controlled the galvanometer mirrors of the OCT device to acquire only 30% of measurements under a continuous trajectory (here a spiral). To do this, we used a spiralscanning trajectory. Actually, each position = ( , ) of the rasterized trajectory is considered as the control input of the galvanometer mirrors. In Fig. 11 is depicted simultaneously, the ground-truth OCT volume (acquiring using 100% of data), the 30% acquired data, and the reconstructed OCT volume. One can remark that the reconstructed volume is faithful to the ground-truth one.

In the second experimental test, we changed the retina sample by a part of a grape (less translucent sample compared

Discussion

The obtained results are satisfactory according to the metric evaluation criteria. It has been shown that it is possible to design non-random sampling matrices that are compatible with the physical implementation as widely studied in the related literature [START_REF] Donoho | Compressed sensing[END_REF]; [START_REF] Candès | Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information[END_REF]; [START_REF] Rauhut | Compressive sensing and structured random matrices[END_REF]. The studied sensing matrix offer similar features in terms of sparsity and incoherence comparing to the well-established random or pseudo-random methods.

In addition, bi-dimensional shearlet [START_REF] Duflot | Shearlet transform: A good candidate for compressed sensing in optical coherence tomography[END_REF] and wavelet [START_REF] Chitchian | Wavelet denoising during optical coherence tomography of the prostate nerves using the complex wavelet transform[END_REF] decomposition were studied the sparsity tool in CS paradigms applied to improve 2D OCT acquisition. However, to the best of our knowledge, it was the first time that tri-dimensional shearlet decomposition was implemented in a physical OCT system. Additionally, continuous trajectories combined with rasterization technique were already reported in the literature, however, they were designed to operate off-line using recorded image data [START_REF] Chauffert | Variable density sampling with continuous trajectories[END_REF].

Furthermore, the developed CS algorithm is not specific to optical coherence tomography modality but can also be used in different 3D imaging methods such as MRI or CT-scan. As a consequence, we evaluated the CS algorithm with 3D sensing matrices more compatible with the acquisition scheme of this type of imaging devices (i.e., MRI and CT).

Although the proposed methods and materials are intended to improve the acquisition frame-rate of an OCT imaging system, there are still improvements/optimizations to be considered. For instance, some parts of the CS algorithm consume much computation time, in particular those relating to the Fourier transform and its inverse. These parts can be implemented directly on a GPU [START_REF] Blanchard | Gpu accelerated greedy algorithms for compressed sensing[END_REF] to save up to 10 times more computing time.

Conclusion

In this paper, in order to overcome the limitations of an FD-OCT acquisition system, we proposed a new compressed sensing paradigm. The developed methodology consisted of four main parts: 1) design of different continuous scanning trajectories (off-line), 2) rasterize the subsampling trajectories to meet the physical conditions (off-line), 3) transform

Figure 1 :

 1 Figure 1: Sparsity comparison between wavelet, curvelet, and shearlet supports using Lena photography as benchmark.Note that, the histogram obtained by plotting the computed coefficients, sorted by their decreasing absolute values. One can reminder that the more the decomposition offer a large number of null-coefficients (or close to zero), the more the method allows a higher sparsity.

Figure 2 :

 2 Figure 2: Illustration of the fundamental working of a cone-adapted continuous shearlet system. The latter is able efficiently cover whole ℝ 2 and giving equal treatment of all directions.

Figure 3 :

 3 Figure 3: Illustration of the designed subsampling masks (black dots represent positions of the acquired measurements (20%)).

Figure 4 :

 4 Figure 4: Illustration of the rasterization algorithm operating in case of a spiral-like scanning trajectory.

Figure 5 :

 5 Figure 5: Comparison of the notion of incoherence between different subsampling methods. The numerical values are obtained for sensing matrices of 105×105 pixels.

Figure 6 :

 6 Figure 6: OCT volumes (here only a 2D slices are depicted) used to access the performances of the proposed algorithm: a) Shepp-Logan phantom, b) -OCT slice of a part of grape, and c) -OCT slice of the retina of a fish eye.

Figure 7 :

 7 Figure 7: [Sheep-Logan phantom] Comparison between (a) the original image and (b) the reconstructed one using 30% of measurements.

Figure 8 :

 8 Figure 8: [Grape] Reconstruction of OCT image using 20% of measurement and different subsampling techniques.

Figure 9 :

 9 Figure 9: Visual comparison between the original and the reconstructed OCT volumes.

Figure 10 :

 10 Figure 10: a) the global view of the OCT device (b-c) depicts the both the B-scan and C-scan acquisition modes, respectively.

Figure 11 :

 11 Figure 11: Experimental validation using a fish eye as biological sample (30% of measurements).

  to the retina). The obtained results are shown in Fig.12. Again, the reconstructed OCT volume (Fig.12(c)) is of good accuracy compared to the ground-truth one (Fig.12(a)).

Table 1

 1 Quality evaluation of the validation tests on both the retina of a fish eye and the Shepp-Logan phantom images using different subsampling rates and scanning trajectories.

					image type
	subsampling rate (%) continuous trajectory type	retina of a fish eye	Shepp-Logan phantom image
			PSNR (db) SIMM	PSNR (db) SSIM
		rosette	38.944	0.834	40.771	0.973
	10%	spiral	39.059	0.833	41.440	0.974
		Lissajous	37.825	0.859	40.908	0.973
		rosette	39.902	0.849	41.730	0.976
	20%	spiral	39.961	0.847	42.171	0.978
		Lissajous	38.911	0.882	41.790	0.977
		rosette	41.048	0.868	42.472	0.980
	30%	spiral	41.129	0.870	42.827	0.982
		Lissajous	40.103	0.910	42.420	0.981
		rosette	43.361	0.898	44.043	0.987
	50%	spiral	43.767	0.903	44.891	0.989
		Lissajous	42.006	0.945	44.004	0.988
		rosette	45.921	0.917	45.278	0.991
	70%	spiral	45.893	0.917	45.436	0.991
		Lissajous	44.202	0.969	45.073	0.991

In the general case where is not the unique sparsest solution of (3) using this approach for recovery is of course possibly not relevant. Moreover, in such a case, this problem has several solutions with equal 0 -"norm" and one may rather define Δ 0 ( ) as an arbitrary element of the solution set.

note that this eqnarray is implicit since depends on , but on the other hand, good algorithms for tuning exist as shown in[START_REF] Chrétien | Sparse recovery with unknown variance: a lasso-type approach[END_REF]

Acknowledgment

This work was supported by the French "Agence Nationale de la Recherche (ANR)" through the NEMRO project [ANR-14-CE17-0013-001] and the INSERM ROBOT Project: ITMO Cancer no 17CP068-00.

the OCT measurements in the shearlet supports for sparsifying, and 4) recover the entire OCT data using only few k-sparse data (on-line).

The proposed compressed sensing algorithm was successfully validated numerically and experimentally. First, this algorithm was evaluated: 1) quantitatively using various subsampling rate (from 10% to 70%) obtained thanks to the designed continuous scanning trajectories (e.g., spiral, rosette, and Lissajous), and 2) qualitatively using similarity scores such as PSNR and SSIM computed for each scenario (OCT image, subsampling rate, scanning trajectory type, etc.).

As claimed, the CS algorithm was also implemented directly on a commercial FD-OCT system. This means that, the galvanometer mirrors that equip the OCT system were controlled directly (using the coordinates of the scanning trajectories as inputs on the control loop) to acquire a limited number of A-scans. This data was then used to recover the whole OCT volume.

The obtained results in both simulation and experiments showed promising performance in terms of quality of the reconstructed OCT images and volumes. This fact is also confirmed by comparing the obtained reconstruction with the ground-truth data.

Future work will consist of the optimization of the CS code to improve its execution time. Actually, several parts of the algorithms were scripted without optimization in a MatLab framework. The ideal would be to implement the current scripts in C++ or in GPU (Graphic Processing Unit) for the direct and inverse shearlet transform.

Conflicts of interest

The authors declare that they have no conflicts of interest.


Andrey V Kudryavtsev

Valerian Guelpa

Patrick Rougeot

Olivier Lehmann

Sounkalo Dembele

Peter Sturm

Nadine Piat

Valérian Guelpa

Sounkalo Dembélé

Nadine Le Fort-Piat

Autocalibration method for scanning electron microscope using affine camera model

Keywords: Autocalibration, SEM, Affine camera, Global optimization

Introduction

The use of vision in robotics currently requires camera calibration, i.e. the knowledge of its motion (extrinsic parameters) and model (intrinsic parameters). In the present work, the aim is to achieve calibration of a Scanning Electron Microscope (SEM) directly from images of a specimen (like pollen grains presented in Fig. 1).

For SEM, the perspective effects can be neglected and the affine camera model is considered [START_REF] Mundy | Geometric invariance in computer vision[END_REF][START_REF] Wang | Guide to Three Dimensional Structure and Motion Factorization[END_REF]. But even if this subject is well studied, the calibration can be very complex Address(es) of author(s) should be given due to the following reasons. First, in most cases it requires a calibration object. It often means a special step of fabrication, which can be very expensive and time consuming especially when working with SEM as, at microscale, the calibration object is of micrometric size. Moreover, it is very difficult to guarantee the quality of its fabrication, which has a profound impact on the accuracy of further image processing. Secondly, the classic calibration [START_REF] Zhang | Flexible camera calibration by viewing a plane from unknown orientations[END_REF] needs to be done offline, which can be very restrictive in some applications where the calibration object can't be placed in front of the camera once the operation started. Thirdly, which includes partially the second point, there is a problem of maintainability of calibration parameters. In order to re-calibrate a camera, the main operation task has to be stopped.

It is important to mention that computation of intrinsic camera parameters is inseparable from camera motion, i.e. extrinsic camera parameters. In fact, autocalibration is a method of calibration, allowing to recover both the intrinsic and extrinsic camera parameters, and is carried out using the same images required for performing the visual task [START_REF] Olivier D Faugeras | Camera self-calibration: Theory and experiments[END_REF]. Most autocalibration methods use a projective calibration as a starting point: it is a calibration from point correspondances only, without any knowledge about camera and scene (scale factor, aspect ratio, skew between image axis, scene orientation, position of principal point, etc.), and it is different from the metric, i.e. true, one up to a 3D projective transformation. Affine calibration and euclidean calibration differs from metric one up to an affine transformation and a similarity transformation, respectively. In the literature, Euclidean calibration is often referenced as similarity or metric calibration [START_REF] Richard | Projective reconstruction and invariants from multiple images[END_REF]. Furthermore, in this work, the term metric will be referred only to a true calibration with the known scale. In our case, as there is no particular information about the object size only Euclidean calibration is achievable, i.e., an up-to-scale calibration. There is a variety of methods for projective calibration computation and the most common ones are factorization-based methods [START_REF] Sturm | A factorization based algorithm for multi-image projective structure and motion[END_REF][START_REF] Dai | Element-wise factorization for n-view projective reconstruction[END_REF]. Then, once the projective calibration is obtained, the goal of autocalibration algorithms is to determine a rectifying homography H from autocalibration constraints and then to transform the calibration to an Euclidean one. However, these methods are often blamed for instability [START_REF] Oliensis | A critique of structure-from-motion algorithms[END_REF].

In case of affine camera, by using a factorization technique, it is possible to directly obtain an affine calibration [START_REF] Tomasi | Shape and motion from image streams under orthography: a factorization method[END_REF]. It means that the plane at infinity is already in its canonical position and the goal is hence to determine the intrinsic parameters of the camera in order to upgrade the calibration from affine to Euclidean. A variety of methods are based on the use of some calibration constraints (zero skew, known aspect ratio, etc.) or motion constraints (pure rotation, planar motion, etc.) [START_REF] Pollefeys | Stratified self-calibration with the modulus constraint[END_REF][START_REF] Triggs | Autocalibration from planar scenes[END_REF][START_REF] Sturm | Critical motion sequences for monocular selfcalibration and uncalibrated euclidean reconstruction[END_REF]. One of the most significant article in the field of autocalibration of affine camera is the work of Long Quan [START_REF] Quan | Self-calibration of an affine camera from multiple views[END_REF]. In his paper, the author proposes a method allowing to upgrade the affine calibration (in this case obtained with a factorization algorithm) to the Euclidean one by using an optimization algorithm. The goal is to find a non-singular 3 × 3 homography matrix and the criterion of optimization is based on Euclidean motion constraints. One of the steps of estimation contains a Cholesky decomposition, which accepts only positive-definite matrices. A priori, this condition is not satisfied in the presence of noise in real images [START_REF] Hartley | Multiple view geometry in computer vision[END_REF]. Thus, the author proposed an elegant way to avoid this problem by imposing the constraints on the matrix to decompose and to assure by this that it is positive-definite. However, the experiments that will be presented further have shown that the method of Long Quan (further referenced as LQ) can not guarantee the respect of the metric constraints, such as aspect ratio close to one, at least in case of SEM.

As many other selfcalibration algorithms, LQ uses local optimization that aims to find a local minimum of the objective function, i.e. the minimum that is the closest to the initial solution. However, it is not sufficient for the application presented here, because the objective function contains non-linearities which results in non-convex cost function. It means that one should use more complex techniques allowing to expand the search space. More recently, this problem was addressed by the methods of global optimization. Fusiello et al. in [START_REF] Fusiello | Globally convergent autocalibration using interval analysis[END_REF] addressed this problem by an interval branch-and-bound method employed for numerical min-imization based on constraint on the fundamental matrix. In [START_REF] Chandraker | Autocalibration via rank-constrained estimation of the absolute quadric[END_REF], the global optimization was used to compute the dual image of the absolute conic to find a rectifying homography. However, even if these global optimization algorithms can guarantee a theoretical global optimality, their lock is a computational time which turns out to be a critical issue for some applications. According to Heinrich et al., the fundamental limitation of these algorithms is that the minimized objective function has no particular geometric meaning which results in instability of these methods [START_REF] Heinrich | Maximum likelihood autocalibration[END_REF]. In return, they proposed a method based on a maximum likelihood objective function.

This paper presents a new method of autocalibration of SEM, considered as an affine camera, which allows to directly obtain Euclidean set of camera matrices by means of global optimization, without passing by affine calibration. It is worth mentioning that it is considered that camera intrinsic parameters are constant and the scene is rigid, i.e. the change of 2D projections across the images is due only to the camera motion. The crucial topics contributing to the success of the algorithm are discussed. Section 2 introduces the specificities of Scanning Electron Microscope. Section 3 presents the representation of SEM considered as an affine camera. The formulation of optimization criteria that include error function, bound constraints and regularization is presented in Section 4. Section 5 is about the global optimization algorithms. Next, the robustness and convergence range of the auto-calibration is tested on synthetic images and compared to classical methods in Section 6. Finally the auto-calibration is tested on real SEM images in Section 7.

Background of SEM

From a physical point of view, electronic imaging is based on the interaction between a radiation of high energy electrons and the atoms of matter. Kinetic energy of an electron accelerated with the voltage E is eE, where e is the charge of the electron. The relation of De Broglie shows that the wavelength λ of the electron is inversely proportional to its level of energy, and then to its accelerating voltage E [START_REF] Amelinckx | Electron microscopy: principles and fundamentals[END_REF]. This leads to very low values of wavelength compared to light (radiation of photons), in the range of X-rays radiation. According to Abbe formula, the resolving power of a radiation (the minimal distance required to resolve two objects) is about λ /2, then with electrons the resolving power can be very low, i.e. the imaging resolution very high: for 1 V accelerating voltage, the resolving power is about 1.2 nm that is about 200 times smaller than optical imaging resolving power; for 2 kV accelerating voltage, the resolving power is about 0.02 nm that is about 40 000 times smaller than optical imaging resolving power [START_REF] Amelinckx | Electron microscopy: principles and fundamentals[END_REF] [START_REF] Joseph I Goldstein | Scanning electron microscopy and X-ray microanalysis[END_REF]. Finally electron imaging is a very high resolution imaging, up to atomic scale. Nowadays, SEM is a central observation instrument for material and life sciences, and micro-nano-engineering [START_REF] Ribeiro | Computer vision for nanoscale imaging[END_REF].

SEM has the following working principle. High energy electrons, that collide with the atoms of a sample, induce two types of scattering: elastic and inelastic [START_REF] Goldstein | Scanning electron microscopy and x-ray microanalysis[END_REF]. Elastic scattering comes with elastic shock, usually with kernel: there is no energy loss, instead, electron is diffracted and atomic bonds are broken. After the shock, the electron becomes a Back-Scattered Electron (BSE). Inelastic scattering comes with inelastic shock, usually with an electron of the cloud: there is a loss of energy. After the collision, the atom is ionized and emits a Secondary Electron or SE. The scanning of the sample surface and the collection of back-scattered or secondary electrons leads to a 2D image of that surface, reflecting the atomic number and topography, respectively. Fig. 1 shows some examples of pollen grains images from secondary electrons detector: the topography is very well presented.

A Scanning Electron Microscope comprises at least the following components (Fig. 2):

1. a gun (source) to generate the beam of electrons, 2. a condenser lens to make the beam cylindrical as possible, 3. an aperture and astigmatism lenses to reduce the diameter of the beam, 4. an objective lens to focus the beam onto the specimen, 5. some scanning lenses to scan the surface of the sample by the beam, 6. a detector to collect and amplify the emitted electrons and to generate pixels.

All these components are located inside a vacuum environment in order to avoid every possible interaction between gas molecules and electrons (both emitted and detected).

For the experiments, in this project, we used a Carl Zeiss Auriga FE-SEM along with its own computer, and a remote computer to run user applications (Fig. 3). The SEM, that is equipped with a six degrees-of-freedom micro-manipulator for object manipulation, features an electron column based on Schottky field emission gun and three detectors (SE detector in the chamber, SE detector in the column, BSE detector in the column). We used the SE detector in the chamber. Some specifications of the SEM are: magnification ranging from ×10 to ×10 6 , working distance from 0 to 20 mm, maximum frame rate of 20 Hz.

Geometrical modeling of SEM

While SEM has no matrix sensor, it still can be modeled as a camera. According to the literature [START_REF] Cornille | Automated 3-d reconstruction using a scanning electron microscope[END_REF][START_REF] Sinram | Calibration of a sem, using a nanopositioning tilting table and a microscopic calibration pyramid[END_REF][START_REF] Cui | Calibration of scanning electron microscope using a multi-image non-linear minimization process[END_REF], in case of SEM, the perspective effects can be neglected for magnification values bigger than ×1000 and then an affine model may be used.

Affine camera model assumes that all projection rays are parallel to each other, i.e. the camera centre is at infinity. That is why affine camera model is often referred as parallel projection model. With regard to SEM, consider the following situation which is quite common: the object with a size of 10 µm is visualized using SEM tuned to the working distance of 10 mm. The ratio between the distance cameraobject and the object size is equal to 1000. Evidently, in such conditions, the projection rays are very close to be parallel. This example gives the intuition behind the usage of affine camera model for SEM. Getting back to geometry, the pro- cess of image formation with an affine camera is represented schematically in Fig. 4.

The camera matrix P in this case is different from perspective projection by a parallel projection matrix Π:

P = KΠ R t ( 1
)
where K is the matrix of intrinsic parameters, R t is the motion matrix (R rotation matrix, t translation vector) and

Π =   1 0 0 0 0 1 0 0 0 0 0 1   .
As a result, the affine projection matrix is of the form:

P =   * * * * * * * * 0 0 0 *   (2)
Affine camera does not have principal point as all projection rays are parallel (camera center is infinitely far from the image plane). It has a direct impact on the form of intrinsic matrix K that has the following form in affine case:

K =   α f s 0 0 f 0 0 0 1   = diag( f , f , 1)   α s 0 0 1 0 0 0 1   (3)
with the three following parameters:

-Aspect ratio α that represents the ratio between height and width of a pixel. In ideal situation its value is equal or at least close to one. It is also true for SEM images which was confirmed in [START_REF] Cui | Calibration of scanning electron microscope using a multi-image non-linear minimization process[END_REF]. -Skew parameter s reflects the orthogonality level of x and y axis in image frame. Its value should be close to zero. -Overall scale factor f represents pixel m ratio. In addition to the special form of the intrinsic parameters matrix, parallel projection imposes some other interesting properties. First, image is invariant to the object displacement along the direction of projection (see Fig. 5.a). In other words, image is invariant to the distance between camera and object, which is the case in SEM. For instance, moving the sample closer to the electron beam or moving it away will not change the resulting 2D projection. As a result, the depth coordinate is completely lost in the process of image formation in case of parallel projection.

Secondly, the projection of the object is independent of translations in x and y directions of image frame if the relative coordinates are used both in 3D object frame and in camera frame (see [START_REF] Quan | Self-calibration of an affine camera from multiple views[END_REF]). If such translation is performed, only the position of the object in the image changes, but not the relative disposition of its feature points (see Fig. 5.b). A reference point is often chosen as the centroid of 2D projection which is the centroid of corresponding 3D points at the same time. The exact same points may be obtained by using relative coordinates and translating the centroid to (0, 0) .

In the case of SEM the displacement of 2D projections is due to the object motion (the motion of the robotic stage). However, for the remainder of the manuscript we will consider that the camera is moving and the object is fixed. Both situations are geometrically equivalent, yet, the latter one simplifies a lot the presentation.

Concerning the subject of motion estimation, it is known that motion parameters as well as the camera intrinsic parameters are contained in the camera matrix which has a special form in affine case (see Eq. 2) [START_REF] Hartley | Multiple view geometry in computer vision[END_REF]. It can be decomposed as:

P =   α f s 0 0 f 0 0 0 1   K   1 0 0 0 0 1 0 0 0 0 0 1   R t 0 0 0 1 (4) = A 2×2 0 2×1 0 1 K   r 11 r 12 r 13 t x r 21 r 22 r 23 t y 0 0 0 1   (5)
Thus, the 2D projection q of a 3D point Q is written as (in homogeneous coordinates):

q = PQ (6)
which is equivalent to

q x q y = α f s 0 f r 11 r 12 r 13 r 21 r 22 r 23   Q x Q y Q z   + α f s 0 f t x t y (7)
Finally, taking into account all the described properties and equations, the total number of parameters for an affine camera is:

3 intrinsic + 3 rotations + 2 translations = 8.
In this way, for N im images corresponding to as many camera poses, the total number of parameters is:

3 intrinsic + 3 rotations .N im + 2 translations .N im = 3 + 5.N im .

Problem formulation

Now the challenge is to find a method to compute the intrinsic and extrinsic parameters of affine camera. Taking as a basis a measurement matrix, we will define a principled cost function usable in an global search optimization method.

Measurement matrix

Mathematically, a set of feature correspondences is represented as a matrix W that will be further referred as measurement matrix. Consider that N im images of the same object were taken from different view points by the same but moving affine camera. It is possible to extract a measurement matrix W that contains the projections of N pt 3D points in different images:

W =      W 1 W 2 . . . W N im      =       q 1 1 q 1 2 . . . q 1 N pt q 2 1 q 2 2 . . . q 2 N pt . . . . . . . . . . . . q N im 1 q N im 2 . . . q N im N pt       (8)
with W i the set of features extracted from the i-th image:

W i = q i 1 q i 2 . . . q i N pt (9)
Several remarks are to be done about the measurement matrix:

the number of columns is equal to the number N pt of extracted features, one column corresponds to N im projections of one 3D point, the number of rows is equal to 3N im in homogeneous coordinated, all features are viewed in all images.

Cost function formulation

The problem of autocalibration can be formulated as the following optimization problem:

arg min ξ ∈R n f (ξ ) (10)
where ξ is a vector of camera's parameters (both intrinsic and extrinsic) and f (ξ ) is the cost or objective function representing the error. Obviously, the formulation of this function is a decisive point to obtain good properties for autocalibration such as robustness and accuracy. Knowing that every element of W can be obtained by the multiplication of 3D point coordinates and the matrix of camera in which it is projected, the expression for the estimated measurement matrix ( Ŵ) can be written as follows:

Ŵ =      Ŵ1 Ŵ2 . . . ŴN im      =      P 1 P 2 . . . P N im      Q 1 Q 2 • • • Q N pt = PQ ( 11
)
where P is a stack of camera matrices (3N im × 4 matrix containing camera matrices for all views and Q is the set of 3D points in homogeneous coordinates with the size 4 × N pt . Every camera matrix (P i ) has 8 degrees-of-freedom. Furthermore, it will be considered that images are taken with the same camera, thus only 5 extrinsic parameters have to be estimated separately for every image. Thus, we need to find such P and Q that would minimize the difference between W and Ŵ. In this formulation the total number of parameters is equal to:

length(ξ ) = 3 intrinsic + 3N im rotation + 2N im translation + 3N pt 3D points (12)
The first item in this equation stands for three intrinsic parameters (scale factor, aspect ratio and skew) of the camera;

the second and third ones represent motion parameters that are different for each view; the last component is the 3D coordinates of object points. Next, we undertake several steps allowing to reduce the number of parameters. Actually, as general optimization problems involve more variables to be optimized, it is harder to make them well-constrained. If they are not well-constrained, the optimization will proceed in a way that minimizes the global mathematical error, but that does not correspond to a solution of a real problem, i.e. has no physical meaning. Thus, it is very important to reduce the number of parameters as much as possible.

Step 1. On the basis of only images taken, only Euclidean calibration could be achieve, and not metric one. In order to make the upgrade to metric calibration, one can use the information about the pixel size for a given magnification which is generally given by SEM manufacturer. Here the scale factor f will be considered constant and equal to one, because the metric calibration is not needed. In this way, only two intrinsic parameters have to be founded.

Step 2. As the position of the world frame is unknown, we are free to fix its orientation equal to the frame of the first camera, so that: R 1 = 1 0 0 0 1 0 which allows to exclude three parameters:

length(ξ ) = 2 + 3(N im -1) + 2N im + 3N pt (13)
Step 3. The biggest term in [START_REF] Sturm | Critical motion sequences for monocular selfcalibration and uncalibrated euclidean reconstruction[END_REF] corresponds to 3D points. One of the possible ways to eliminate it is to use the pseudoinverse of the matrix P. This solution reduce drastically the number of parameters, in return for more calculations. Indeed, the stack of 3D points Q may be replaced by:

Q = P + W (14)
where P + is the pseudo-inverse of the matrix P. Here, we obtain it using SVD (Singular Value Decomposition) algorithm. Therefore, [START_REF] Triggs | Autocalibration from planar scenes[END_REF] transforms into:

Ŵ = PP + W (15)
and the number of parameters is reduced to:

length(ξ ) = 2 + 3(N im -1) + 2N im (16)
Step 4. The number of parameters can be further reduced by using relative coordinates and eliminating translations. Indeed, as stated in section 3 and Fig. 5.b, affine model is invariant to translations in image plan. The use of relatives coordinates will allow to eliminate t i and the number of parameters to estimate can be reduced to:

length(ξ ) = 2 + 3(N im -1) (17)
and the equation ( 15) takes the following form:

Ŵ = MM + W r (18)
In this expression, W r is the measurement matrix in nonhomogenious relative coordinates and M is:

M =      M 1 M 2 . . . M i      =      AR 1 AR 2 . . . AR i      ( 19
)
where M is a 2 × 3 upper-left components of affine camera matrix P.

Recall that

A = f α s 0 1 ( 20
)
and R i is the first two rows of the rotate matrix, according to [START_REF] Kudryavtsev | Full 3d rotation estimation in scanning electron microscope[END_REF]:

R i = R z ( i θ i-1 )R y (ρ i,i-1 )R z ( i-1 θ i ) (21)
where the angles θ are in-plane rotations and ρ is outof-plane rotation.

In-plane rotation can be estimated from a pair of affine images, but there is an issue with out-of-plane rotation estimation because of bas-relief ambiguity: from two affine views, the scene (motion and structure) may be recovered only up to the combination of depth variation and out-ofplane rotation [START_REF] He | Modified fast climbing search autofocus algorithm with adaptive step size searching technique for digital camera[END_REF].

With three or more views, this ambiguity disappears in general. In our implementation, we use the self-calibration approach of [START_REF] Kudryavtsev | Full 3d rotation estimation in scanning electron microscope[END_REF] which computes directly, from the affine factorization of three views, the scene (structure and pose of the three views), up to an Euclidean transformation, as well as the intrinsic parameters. We use this to initialize a global non-linear optimization of all variables, as described in the following. Note that using 3 views for this initialization step makes it straightforward and efficient to embed it in a robust estimation scheme such as RANSAC.

Ideally, we are looking for such set of parameters ξ that ensure Ŵ = W r . However, due to the presence of noise in the estimation of image features, this equality is never satisfied exactly, which means that a way of comparison needs to be found. As in common, a geometric distance d between estimated and measured points will be used in present work. Hartley and Zisserman in [START_REF] Hartley | Multiple view geometry in computer vision[END_REF] specify that the minimization of geometric error between measured and estimated 2D points is equivalent to finding such a Ŵ as close as possible to W r in Frobenius norm. Thus, the autocalibration task converges to an optimization problem with the following cost function:

f (ξ ) = N im ∑ i=1 N pt ∑ j=1 d( qi j , M i (ξ )M + i (ξ ) qi j ) 2 (22) = W r -Ŵ 2 F (23) = W r -M(ξ )M + (ξ )W r 2 F (24) = W r -MM + W r 2 F ( 25
)
The goal is to find a global minimum of this function subject to a set of constraints that will be developed in the next sections.

It is important to add that, in Euclidean calibration, the matrix M i is the result of the product of two matrices: an upper-triangular matrix of intrinsic parameters A i and a rotation matrix. Such matrices can be obtained using RQdecomposition of M i matrix. It is important to notice that the result of such decomposition is unique if rank(M i ) = 2. Thus, by computing M i as the product of upper-triangular matrix and a rotation matrix, we ensure that the obtained calibration is an Euclidean one.

Regularization

Many optimization algorithms encounter the following problem: they do not take into account the physical meaning of the optimization parameters, nor the stochastic nature of noise. Thus, due to the presence of noise, they find a solution that actually makes the global error smaller, but the values of parameters do not correspond to a realistic model. Moreover, the autocalibration problem may be considered ill-posed [START_REF] Hadamard | Sur les problèmes aux dérivées partielles et leur signification physique[END_REF] and strongly non-linear due to presence of multiplication of sine and cosine functions. The method allowing to compensate the issue is called regularization and more specifically Tikhonov regularization [START_REF] Nikolaevich Tikhonov | Solutions of ill-posed problems[END_REF]. The idea behind it consists in introducing more knowledges about the system, by adding weight coefficients to each of the optimized parameter. It allows to improve the likelihood that the minimization converges to the desired solution.

In current cost function [START_REF] Kudryavtsev | Full 3d rotation estimation in scanning electron microscope[END_REF], all parameters have the same weight. However, we would like to tell the algorithm that the values of intrinsic parameters should be close to the initial ones and that the biggest impact should be made on the values of out-of-plane rotations. In other words, during the optimization, the error should be minimized mainly by adapting motion parameters and not the intrinsic ones.

Mathematically, it translates in adding a standard regularization term of the following form:

r = Γ ξ 2 F ( 26
)
where Γ is a regularization diagonal square matrix that has the same number of rows as ξ . Each value on the diagonal represents the desirable impact of corresponding parameter, i.e. the weight. The bigger the value, the less the algorithm would tend to change the given parameter. Note that the definition of this regularization term is currently based on the assumption of a normal distribution of parameters, but this model will be refined in futur works.

In this work, we defined all weights as powers of 10:

-100, for intrinsic parameters -0.1, for slope angles θ 1 and θ 2 -0.01, for out-of-plane rotation ρ These values were found well suited experimentally.

Finally, the cost function with regularization term is written as:

f r (ξ ) = W r -M(ξ )M + (ξ )W r 2 F + Γ ξ 2 F (27)

Bound constraints

Finding a minimum of the objective function obtained previously is the subject of nonlinear optimization, which is typically a challenge without any additional information and thorough understanding of the nature of parameters. In presented formulation of the objective function all of the parameters ξ have an actual geometric meaning which allows to largely reduce the search space of the solution by implementing the bound constraints. It is all the more pertinent for global optimization, where starting points are actually generated inside predefined bounds. It means that tighter bounds lead to faster convergence and higher probability of finding the solution.

The constraints for the elements of M i matrix can be defined as follows. First, as regards the matrix of intrinsic parameters, as the common scale factor f is factored out, the constraints can be easily imposed: the value of aspect ratio should be close to one, because the pixels are generally squared, and the skew factor should be close to zero, which would denote that x and y axis of camera are perpendicular. Furthermore, as the metric calibration can not be obtained (the only possible is the Euclidean "up-to-scale" calibration), in case of constant focal length, the value of f can be fixed to any positive real value, e.g., to one. Secondly, the rotation matrix is decomposed into a sequence of three elemental rotations in a spherical coordinate system with angles θ 1 , ρ and θ 2 . Hereafter, we speak about intrinsic rotations which means that they don't occur about the axes of the fixed coordinate system, but about the axes of the rotating coordinate system, which changes its orientation after each elemental rotation. It results in a following constraints for the angles: θ 1 and θ 2 can vary in a range of (-90 • ; 90 • ). The angle ρ, the out-of-plane rotation may vary in range that depends on the acquisition process. In typical situation, this angle is equal to several degrees. Here it will be fixed in the

Initial value Bounds

Aspect ratio α 1 [0.5; 1.5] Skew s 0 [-0.5; 0.5]

θ 1 0 [-90 • ; 90 • ] ρ 0 [-90 • ; 90 • ] θ 2 0 [-90 • ; 90 • ] Total number of varying parameters 2 + 3(N im -1)
range (-90 • ; 90 • ). Thus, the constraint for all elements of M i matrix are combined in Table 1.

Here, constraints values are chosen larger than what is needed, in order to demonstrate the high convergence range of this method. However, applications could be done with reduced bounds if preliminary data are known.

Problem solution based on global search

At this point, two cost functions defined for autocalibration, f (ξ ) and f r (ξ ), the latter one containing the regularization term. Generally, local optimization is often used in computer vision and autocalibration in particular, which works well for convex problems. However, for strongly non-linear cost functions, in order to find the minimum, it is preferable to use global optimization algorithms.

Global optimization algorithms can be subdivided in two main classes [START_REF] Weise | Global optimization algorithms -theory and application[END_REF]: deterministic and probabilistic. In case of deterministic algorithms, the search space can be subdivided into multiple pieces, similarly to the divide and conquer strategy. This step is then followed by the exploration of each smaller region by a local solver. After that, the results are combined and the best one is taken as a global minimum. Deterministic methods can provide a certain level of assurance that the global optimum will be located. However, even if they can guarantee that the solution found is the global one, no algorithm can do it in a finite time [START_REF] Moles | Parameter estimation in biochemical pathways: a comparison of global optimization methods[END_REF]. It comes from the fact that the smaller the regions on which the space search is subdivided, the more times the local solver should be launched. At the same time, it becomes evident that if the sizes of the regions are close to zero, the likelihood of finding a global minimum increases. In return, the probabilistic algorithms generate the solution based on random variables. They include algorithms like simulated annealing and evolutionary algorithms. In this case, the solution can be found with relatively high time efficiency, however, the global optimality can not be guaranteed. In present work, we decided to apply both types of algorithms: Scatter Search with a local solver, presented by Ugray et al. in [START_REF] Ugray | Scatter search and local NLP solvers: A multistart framework for global optimization[END_REF], and the Genetic Algorithm [START_REF] Goldberg | Genetic Algorithms in Search, Optimization, and Machine Learning[END_REF]. For both algorithms, their MatLab implementations were used.

The Scatter Search algorithm (further referenced as GS for Global Search) is a population based meta-heuristic algorithm devised to intelligently perform a search on the problem domain algorithm. It consists in a generation of multiple trial points within finite bounds, which are candidate starting points for a local solver. These are then filtered to provide a smaller subset from which the solver attempt to find a local optimum. Then it evaluates remaining possible candidates and start a local solver. The best solution is retained as the global minimum of the problem. More detailed description of the scatter search algorithm is given in [START_REF] Glover | A template for scatter search and path relinking[END_REF].

The Genetic Algorithm (GA) is an heuristic algorithm of search, which is used to solve different optimization problems using random search, mutation and variation of parameters. It belongs to the class of evolutionary algorithms, which generate solutions by the techniques inspired by natural evolution. GA differs from other evolutionary algorithms by putting an accent on the use of cross-over operator, which conducts the operation of recombination of solution candidates.

6 Tests on virtual images

Accuracy on nominal case

In order to find the best solution, we tested four combinations:

-GSR, scatter search optimization with regularization -GS, scatter search optimization without regularization -GAR, genetic algorithm optimization with regularization -GA, genetic algorithm optimization without regularization

As stated previously, MatLab implementations of these algorithm were used. Their performance were compared to the state-of-the-art optimization algorithm of Long Quan [START_REF] Quan | Self-calibration of an affine camera from multiple views[END_REF], refereed as LQ in the text. The factorization method presented by Tomasi and Kanade [START_REF] Tomasi | Shape and motion from image streams under orthography: a factorization method[END_REF], refereed as TK, will be also used as reference. It consists in inferring scene geometry and camera motion from a stream of images, recovering shape and motion under orthography, without computing depth. This method mainly uses the singular value decomposition of the measurement matrix Ŵ to find the rotation matrix and the shape matrix. This method stays a standard to analyze the displacements of an object from multiple views.

A sequence of four synthetic images of a diamond (see Fig. 6) was used. The 3D structure and the cameras parameters are known, so the measurement matrix contains perfect noise-free correspondences. The faces between vertices are given only in a purpose of better visualization of the object. As it was mentioned previously, it is not possible to obtain the metric calibration without any additional information on the object structure or camera magnification. Thus, in order to analyze the results, only the properties that are preserved under similarity transformation might be used. Among them are ratio of lengths and angles. For the analysis four vertices of the object were chosen: a, b, c, d (see Fig. 7). The real values of the ratios of lengths and the angles between vertices (α ab , α bc , α cd ) as well as the estimated ones are presented in Table 2. It can be noticed that all global methods give a similar results, the estimation error for both ratios and angles is lower than thousandth of percent. This result is not highly significant, but confirm the viability of the presented approach and allows to proceed to more tests on noisy or real images. It is also interesting to look at the execution time (Table 3). It was measured for MatLab R2017b implementation of global optimization on the computer with the following parameters: Windows10 x64, 3.80 GHz Intel Xeon CPU, 32 GB of RAM. One can remark, that for the same output results, the genetic algorithm takes 100 times more time than scatter search method. Thus, genetic algorithm will no longer be used for tests.

Robustness to noise

We have several methods (LQ, GS, GSR,TK) that perform equally well for noise free images. However, this situation is far from reality as the set of feature coordinates always contain noisy measurements. Hence, it is important to test the performance of the algorithms at different noise levels. The procedure is exactly the same as above, except for a Gaussian noise with different variance is added to 2D image coordinates. For every noise level, the autocalibration was repeated 100 times. The results are represented in Fig. 8.

At this point we remark that GSR and TK outperforms both the LQ and GS methods. Even for high noise level, the error for estimation of length ratios stays inferior to 2%. In contrast, it may reach 4.5% for LQ method and 3.5% for the scatter search optimization without regularization (GS). The situation with angles is similar, but the difference between methods is lower. These results prove at the same time the efficiency of presented method and the importance of regularization term.

Convergence range

As it was already shown, global optimization is a time consuming technique. Typically, for a sequence of four images with 22 points, the autocalibration takes about 0.5 seconds. So it is important to verify if it is worth using global optimization instead of local one. To test that, we conducted the following experiment. Once again, four images are obtained in a way that the out-of-plane rotation between the first image and the second one is 5 degrees. Then, we run autocalibration with different initial conditions using GSR method and the same method but with Levenberg-Marquardt (LM) local optimizer (LMR, for Levenberg-Marquardt with regularization). The results are summarized in Table 4. As expected, with LMR algorithm, optimization suffers from multiple local minima of the cost function, so is needed a very good initial estimate to use this algorithm. In contrast, global optimization allows indeed the arbitrary choice of starting point as it is capable to compensate the error of 45 degrees.

Tests on real images

Five image datasets were used:

affine camera images:

-hotel: 10 images from the hotel image sequence [START_REF] Conrad | A paraperspective factorization method for shape and motion recovery[END_REF]. -SEM images from [START_REF] Ahmad P Tafti | 3dsem: A dataset for 3d sem surface reconstruction[END_REF]:

-brassica: 4 images of pollen grain of white turnip plant from a Hitachi S-4800 FE-SEM. -grid: 5 images of TEM copper grid from a Hitachi S-4800 FE-SEM. -SEM images acquired at FEMTO-ST Institute, taken using a SEM Zeiss AURIGA 60:

-pot: 7 images of pollen grain of aquatic, mostly freshwater, plant of the family Potamogetonaceae. -pot2: 6 images of another Potamogetonaceae pollen grain.

The hotel sequence is the image dataset used by Long Quan to validate his algorithm. In our implementation, 196 matches were extracted and the results of autocalibration, as expected, were exactly the same for all four methods (TK, LQ, GS, GSR). An example of camera poses as well as the sparse 3D calibration is shown in Fig. 9.

Back to SEM images, three datasets were analyzed. As all image datasets were taken using SEM, there is no available information on the location of the camera. Thus, we will use two criteria allowing to judge on the algorithm performance. First, the respect of metric constraints and, secondly, the reprojection error (Fig. 10). It is important to note that we use a leave-one-out method to do the reprojection: for a set of N images, we extract each image one time to estimate the 3D-points, then reproject these points according to the estimated pose for this image. In this way, the method's robustness has to be better than if the 3D-points are defined from the N images.

We consider that metric constraints are respected if the aspect ratio is within the interval (0.9; 1.1) although we know that it has to be close to one. For LQ and GS method, the constraints are never respected. The TK cannot estimate the metric constraints. GSR algorithm gives the desired result for all datasets as the constraints on aspect ratio and skew are already included in cost function formulation.

In terms of reprojection error, the algorithms gave good results for TK and GSR. This result is important because it proves that even if we prevent the change of intrinsic parameters, the GSR algorithm succeeds in keeping the same level of accuracy as TK.

Obviously, these two criteria do not allow full validation of the autocalibration: one can consider them as mandatory but not sufficient. Unfortunately, when working at such scales and with small objects, there is no measurement device capable to provide the ground truth. However, for some datasets, we know the actual angle of tilt: 3 degrees for brassica, 7 degrees for grid, and 3 degrees for pot. For other datasets, the rotation between images were performed as a combination of two articular robot movement which was not calibrated, which means that the angle of out-of-plane rotation cannot be extracted. To compare the rotation angles, the following procedure was performed. Once the sequence passed through autocalibration, the rotation matrices for every frame were extracted and then transformed into axis-angle representation. Therefore, if the rotation were performed in one movement, these angles should match the angles given by robot sensor. The results are presented in Table 5. The difference do not exceed 0.21 degrees on all tries.

Conclusion

In this paper, a new method of autocalibration for affine camera was presented. All its components have been computed: model as well as motion parameters. The presented method, being based on global optimization, has the following advantages comparing to the state-of-the-art techniques:

-All metric constraints are imposed directly on the optimized parameters. Moreover, with regularization, it is possible to guide the optimization towards the desired result, i.e. add a penalty score for excessive change of such parameters as aspect ratio and skew. -All optimization parameters have a physical meaning. It means that all known elements about camera location can be easily imposed. For example, one can impose the equality of out-of-plane rotation if the angle of tilt is known to be constant for all image sequence. -Thanks to the deep analysis of camera properties, the final number of optimized parameters is low and good initial estimate is provided from the elements of funda- With regard to microscopy community, we presented the first autocalibration algorithm for Scanning Electron Microscope.

It is worth adding that the crucial step in all autocalibration algorithms is the feature matching across the images. In case of affine camera less points are needed to constraint camera parameters, because of the reduced number of degrees of freedom. As it was shown before, we used the number of features from 20 to 200. If it is possible to increase the number of matches without loss of their quality, it will improve the final results, however, will slightly increase the execution time. Though, the quality of the correspondences

Fig. 1

 1 Fig. 1 Examples of SEM images with 100 µm pollen grains.

Fig. 2

 2 Fig. 2 Basic SEM internal structure.

Fig. 3

 3 Fig. 3 Carl Zeiss Auriga SEM used in experiments.

Fig. 4

 4 Fig. 4 Image formation under assumption of affine camera model.

Fig. 5

 5 Fig. 5 Illustration of affine camera properties: a) invariance to depth variation; b) invariance to in-plane translations of camera or object.

Fig. 6 Fig. 7

 67 Fig. 6 Four images of diamond sequence with 3D view. Cameras are represented by local coordinate axes.

Fig. 8

 8 Fig.8Robustness of different autocalibration algorithms to noisy measurements. TK and GSR obtain best results, mainly overlaid on curves.

Fig. 9

 9 Fig. 9 Result of autocalibration for hotel image dataset.

1 TKFig. 10

 110 Fig.10Analysis of reprojection error for different autocalibration algorithms: TK, LQ, GS, GSR. Red background means that metric constraints are not respected: typically, aspect ratio is greater than 1.1 or lower than 0.9. Grey background means that metric constraints are not calculated. For points inside circles, the reprojection error is lower than: 1 pixel for blue circle, 0.5 pixel for red circle. Numbers reflect the percentage of points inside the corresponding circle.

Table 1

 1 Examples of initial values and bound constraints for autocalibration. The rotation of the first camera is fixed and equal to I 2×3 .

Table 2

 2 Comparison of performance of different algorithms on a noise free sequence of synthetic images.

			Ratios of lengths			Angles
		a/b	b/c	c/d	α ab	α bc	α cd
	Real value	0.761	12.661	0.036	94.028	114.747	136.544
	LQ	0.761	12.661	0.036	94.028	114.747	136.544
	GSR	0.761	12.664	0.036	94.030	114.746	136.533
	GS	0.761	12.664	0.036	94.021	114.743	136.534
	GAR	0.761	12.664	0.036	94.030	114.746	136.533
	GA	0.761	12.664	0.036	94.021	114.743	136.534
	TK	0.761	12.661	0.036	94.028	114.747	136.544
	Errors:
	ε LQ , %	6.56e-13	4.21e-14	1.34e-13	4.53e-13	5.70e-13	6.24e-14
	ε GSR , %	5.91e-03	2.97e-02	1.91e-02	2.69e-03	1.71e-03	8.15e-03
	ε GS , %	2.75e-03	3.19e-02	1.70e-02	7.07e-03	3.99e-03	7.08e-03
	ε GAR , %	5.91e-03	2.97e-02	1.91e-02	2.68e-03	1.71e-03	8.14e-03
	ε GA , %	2.75e-03	3.19e-02	1.70e-02	7.07e-03	3.99e-03	7.07e-03
	ε T K , %	1.75e-13	4.91e-13	4.48e-13	4.08e-13	3.96e-13	1.25e-13

Table 3

 3 Mean execution time of different autocalibration methods for 4 images and 22 points.

	Method	LQ	GSR	GS	GAR	GA	TK
	Time, s	0.0156 0.408 0.842 109.118 101.074 0.0011

Table 4

 4 Comparison of convergence range between local optimizer (Levenberg-Marquardt, LMR) and global optimizer (GSR). Tests were run on noise-free images.

	Initial ρ 12 , in degrees	-20	0	1	2	3	4	5	6	7	8	20	45
	LMR	5.23 0.00 1.69 4.07 4.

Table 5

 5 Comparison between rotation angles in degrees given by robot sensors (true), and autocalibration algorithm (estimated).

			brassica		grid		pot
		true estimated true estimated true estimated
	ρ 12	3	2.9425	7	6.9103	3	3.0312
	ρ 23	3	3.0182	7	7.0851	3	3.0965
	ρ 34	3	2.8270	7	6.7993	3	3.0723
	ρ 45	-	-	7	6.8596	3	3.1131
	ρ 56	-	-	7	-	3	2.9628
	ρ 67	-	-	-	-	3	3.0676

is much more important than their quantity: even if a feature is an inlier, its quality may vary. Outliers are rejected at the step of fundamental matrix estimation.