Alice Fontbonne
email: alice.fontbonne@protonmail.com

Hervé Sauer

François Goudail

Comparison of methods for end-to-end co-optimization of optical systems and image processing with commercial lens design software

We compare three different methods to co-optimize hybrid optical / digital imaging systems with a commercial lens design software: conventional optimization based on spot diagram minimization, optimization of a surrogate criterion based on a priori equalization of modulation transfer functions (MTFs), and minimization of the mean square error (MSE) between the ideal sharp image and the image restored by a unique deconvolution filter. To implement the latter method, we integrate -for the first time to our knowledge -MSE optimization to the software Synopsys CodeV. Taking as application example the design of a Cooke triplet having good image quality everywhere in the field of view (FoV), we show that it is possible, by leveraging deconvolution during the optimization process, to adapt the spatial distribution of imaging performance to a prescribed goal. We also demonstrate the superiority of MSE co-optimization over the other methods, both in terms of quantitative and visual image quality.

Introduction

Currently, almost all imaging systems integrate both a complex optical system and digital processing algorithms. The idea of modeling the entire imaging chain to enable its optimization is not new, but it took several years before the first co-design approaches really appeared. The first to co-design an optical system were Cathey & Dowsky [START_REF] Dowski | Extended depth of field through wave-front coding[END_REF][START_REF] Cathey | New paradigm for imaging systems[END_REF], focusing on the co-design of a single optical element, a phase mask. The next step is to take into account more complex optical systems and to co-optimize them with digital post-processing algorithms. These methods are referred to as "joint optical-digital design" [START_REF] Burcklen | Joint digital-optical design of complex lenses using a surrogate image quality criterion adapted to commercial optical design software[END_REF], "integral design strategy" [START_REF] Wang | An integral design strategy combining optical system and image processing to obtain high resolution images[END_REF], "end-to-end lens design" [START_REF] Sun | End-to-end complex lens design with differentiable ray tracing[END_REF] or "holistic opticaldigital hybrid-imaging design" [START_REF] Vettenburg | Holistic optical-digital hybridimaging design:wide-field reflective imaging[END_REF] in the literature. Recent approaches attempt to bypass the expertise of commercial optical designers by using deep learning techniques for co-design [START_REF] Vincent Sitzmann | End-toend optimization of optics and image processing for achromatic extended depth of field and super-resolution imaging[END_REF][START_REF] Shang | Two-step training deep learning framework for computational imaging without physics priors[END_REF][START_REF] Nehme | Learning optimal wavefront shaping for multi-channel imaging[END_REF][START_REF] Sun | End-to-end complex lens design with differentiable ray tracing[END_REF][START_REF] Liu | End-to-end computational optics with a singlet lens for large depthof-field imaging[END_REF][START_REF] Halé | End-to-end sensor and neural network design using differential ray tracing[END_REF]. However, currently, only commercial optical lens design software (like Zemax ® OpticStudio ® or Synopsys ® CodeV ® ) allows the design of complex optical systems that are truly manufacturable. This is why it is interesting to integrate directly a co-design approach into this type of software. Indeed, the optical designer expertise allow them to use many hints in order to steer the lens optimization process [START_REF] Fisher | Optical System Design[END_REF], which usually consists of a sequence of local optimizations towards evolving intermediate goals to reach a final realistic design with good final image quality. This idea of integrating a co-design metric taking into account the entire hybrid system, namely, the mean square error (MSE), is not new in Zemax. Stork & Robinson [START_REF] Stork | Theoretical foundations for joint digital-optical analysis of electro-optical imaging systems[END_REF] were the first to implement it, followed by Vettenburg & Harvey [START_REF] Vettenburg | Holistic optical-digital hybridimaging design:wide-field reflective imaging[END_REF] and Wang et al. [START_REF] Wang | An integral design strategy combining optical system and image processing to obtain high resolution images[END_REF]. However, to the best of our knowledge and up to now, no publication claims to have implemented it in CodeV. Moreover, the results obtained with a MSE-based criterion and a commercial lens design software have never been compared to other co-design methods. Indeed, there are several methods to optimize a hybrid system from end-to-end. In particular, another joint optical-digital design method has been developed by Burcklen et al. [START_REF] Burcklen | Joint digital-optical design of complex lenses using a surrogate image quality criterion adapted to commercial optical design software[END_REF] to facilitate the use of the co-design paradigm by professional optical designers in CodeV. It is based on a "surrogate" merit function that contains terms classically available in lens design software but used in a non-standard way to implicitly take into account image processing. Our purpose in this work is to implement the MSE-based optimization criterion in the CodeV software, and to compare the obtained performance with two other optimization criteria. The example of application chosen for illustration is the design of a Cooke triplet aimed at providing good image quality over a wide field of view (FoV) of 40 • (±20 • ). The basic elements of the considered hybrid optical-digital system are described in Sec. 2. In Sec. 3, we apply the "conventional" design method, which considers the existence of the digital processing only after the optimization of the optical system. We establish the performance of the obtained hybrid system by evaluating its modulation transfer functions (MTFs), its effective MTF (taking into account the whole optical system and post-processing) and the visual quality of the restored images. In Sec. [START_REF] Code | V Lens System Setup Reference Manual[END_REF], we investigate in the same way the method based on a surrogate criterion implemented in Code V [START_REF] Burcklen | Joint digital-optical design of complex lenses using a surrogate image quality criterion adapted to commercial optical design software[END_REF]. Then, in Sec. 5, we describe the implementation of a method for directly optimizing the MSE under CodeV, and compare the performance of the obtained hybrid system with those provided by the two other methods. We conclude in Sec. 6 on the performance of the three studied methods and on the most efficient way to take into account the existence of a digital post-processing during the lens design process.

Description of the hybrid system of interest

This section aims to describe the basic specifications of the optical/digital hybrid systems considered in this article. Their goal is to obtain the best possible image quality over the whole FoV. A hybrid system is made of two blocks: the imaging system (i.e. the optical system and a sensor) and a digital post-processing algorithm. For the purpose of this study, we have chosen to consider a classical lens architecture, a Cooke triplet (described in Sec. 2.1), and a fast deconvolution algorithm, the average Wiener filter (described in Sec. 2.2). The optimization goal and the different investigated methods to achieve it are described in Sec. 2.3.

Optical architecture

The original purpose of a Cooke Triplet was to serve as a photographic lens. Consisting of three optical elements and a diaphragm, it is recognized for its good correction of aberrations over a wide FoV [START_REF] Dennis | Optical designing as an art[END_REF]. Its variants have been the subject of several studies [START_REF] Sharma | Design of a new five-element cooke triplet derivative[END_REF][START_REF] Sharma | Four-element lens systems of the cooke triplet family: designs[END_REF][START_REF] Darko | Optimization of the cooke triplet with various evolution strategies and damped least squares[END_REF] and it is still used today for the improvement of optical systems [START_REF] Stangner | Cooke-triplet tweezers: more compact, robust, and efficient optical tweezers[END_REF]. The Cooke triplet is thus a basic "building block" of conventional optical design, which is interesting to study in the context of joint optical/digital design. We will thus use it in this article as a basis for the description of the various optimization methods. We choose to set its focal length to 50 mm, its aperture to F/4 and its half-FoV to 20 • (it is enlarged compared to the initial 14 • half-FoV of the CodeV example from which this triplet is derived). We will consider that this lens is used with a panchromatic sensor of pixel size about 5 µm, which sets the Nyquist frequency to f N = 100 lp.mm -1 . Three wavelengths of the visible spectrum are considered, and the constraints applied to the system during optimization are standard (see Supplemental 1, section 1).

Post-processing algorithm

The digital post-processing applied to the image acquired by this system consists of a deconvolution to compensate for the aberrations of the optical lens. We want the deconvolution algorithm to ensure reconstruction regardless of the position in the FoV and to be fast and simple. For this purpose, we choose the linear average Wiener filter:

wΨ (ν) = 1 K K k=1 hψ k (ν) ⋆ 1 K K k=1 | hψ k (ν)| 2 + Snn(ν) Soo(ν) , (1)
where Ψ = {ψ 1 , ψ 2 , ..., ψ K } is a set of K positions in the FoV describing the whole FoV, hψ k (ν) is the local optical transfer function (OTF) of the imaging system at the position ψ k depending on the spatial frequency ν, S nn (ν) is the power spectral density (PSD) of the noise and S oo (ν) is the statistical PSD of the scene. We use a generic ideal image model with a power-law PSD [START_REF] Daniel | Origins of scaling in natural images[END_REF][START_REF] Van Der Schaaf | Modelling the power spectra of natural images: Statistics and information[END_REF] S oo (ν) ∝ ν -2.5 , which well represents natural scenes. The white noise PSD is such that the signal-to-noise ratio (SNR) on the raw image is 34 dB, with SN R = 10 log 10 S oo (ν)dν/ S nn (ν)dν . Note that in order to obtain rotationally symmetrical average Wiener filters, it is better to use a large number of positions in the FoV. In the following, unless otherwise stated, K = 400. Furthermore, in view of the expression of the Wiener filter in Eq. 1, we can define the "average MTF" as:

1 K K k=1 | hψ k (ν)| 2 .
(2)

Optimization goal

Our objective is to design a hybrid system with good image quality over the whole FoV. To reach this goal, we will investigate three different design approaches.

• The "conventional triplet" is obtained with a conventional lens design criterion based on the minimization of the average size of the spot diagram over the FoV. This method does not take into account the digital postprocessing during the optimization process (Sec. 3).

• The "EMTF triplet" is obtained by optimizing the surrogate criterion based on conventional optical metrics introduced in [START_REF] Burcklen | Joint digital-optical design of complex lenses using a surrogate image quality criterion adapted to commercial optical design software[END_REF]. It implicitly takes into account the post-processing (Sec. 4).

• The "MMSE" triplets are obtained by minimizing the average MSE between the deconvolved image and an ideally sharp image of the scene. More precisely, the image produced on the sensor at one position ψ k of the FoV can be modeled by h ψ k (r) * O(r), where O(r) is the sampled ideal scene image (r represents the spatial coordinates) and * denotes the convolution operator. This acquired image is then deconvolved with the filter w Ψ (r), defined in Eq. 1 in the Fourier domain, to restore its sharpness.

The restored image at the output of the hybrid optical-digital system can thus be modeled as:

Ôψ k (r) = w Ψ (r) * [h ψ k (r) * O(r) + n(r)] (3)
where n(r) is the detection noise. The mean-squared error (MSE) is then defined as:

M SE(ψ k ) = E Ôψ k (r) -O(r) 2 dr , (4)
where E[•] represents the mathematical expectation over the noise n(r) and the scene image O(r), which are both assumed to be zero-mean, stationary random processes of PSD S nn (ν) and S oo (ν) respectively. The MSE Eq. 4 can in fact be generically expressed in the Fourier domain by the sole noise and statistical scene image PSDs, the OTF hψ k (ν) of the lens and the restoration filter wΨ (ν) [START_REF] Robinson | Joint design of lens systems and digital image processing[END_REF][START_REF] Diaz | Increase in depth of field taking into account deconvolution by optimization of pupil mask[END_REF]:

M SE(ψ k ) = | hψ k (ν) wΨ (ν) -1| 2 S oo (ν) + | wΨ (ν)| 2 S nn (ν) dν . ( 5
)
The lens optimization criterion is then defined as the MSE averaged over all the considered FoV positions

Ψ = {ψ 1 , ψ 2 , • • • , ψ K }: M SE mean (Ψ) = 1 K K k=1 M SE(ψ k ) (6)
Note that the average Wiener filter defined in Eq. 1 is the linear deconvolution filter that minimizes M SE mean for a given optical system [START_REF] Diaz | Increase in depth of field taking into account deconvolution by optimization of pupil mask[END_REF], which leads to overall optimizations in the M SE mean sense. As all (local) optimization methods, its result depends on the starting point. We will consider two different systems: the "MMSE 1 triplet" will be obtained by taking the "conventional triplet" as starting point, and the "MMSE 2 triplet" will take the EMTF triplet as a starting point (Sec. 5).

In the remainder of this paper, we will describe the optimization of hybrid imaging systems with these three different approaches and compare the obtained performance.

3 Conventional approach

Definition of the conventional approach

To limit the aberrations in an optical system for a selected number of positions in the FoV, the conventional optimization criterion in CodeV is to minimize the sum of the squares of the root mean square (RMS) diameters of the spot diagrams (calculated from a set of desired positions in the FoV). Other criteria exist natively in CodeV: it is possible to minimize the wavefront error (WFE) or to maximize the value of the modulation transfer function (MTF) at selected spatial frequencies. Whatever the chosen criterion, it is always necessary to impose constraints on some parameters during the optimization (for example, set the focal length to a particular value defined by the specifications, or constrain the thickness at the edges of the lenses to have strictly positive values). This can be done by adding penalization terms like p 2 × (f ′ -f ′ target ) 2 (where p is a positive constant, defining the "weight" of this term) to the basic criterion, or by strictly enforcing the constraint using Lagrange multiplier. In this study, we optimize the sum of the squares of the RMS diameters of the spot diagrams under common CodeV optical constraints (see Supplemental 1, section 1). The CodeV optical design software uses by default the Levenberg-Marquardt (or Damped least squares) algorithm to minimize this cost function [START_REF] Levenberg | A method for the solution of certain non-linear problems in least squares[END_REF][START_REF] Marquardt | An algorithm for least-squares estimation of nonlinear parameters[END_REF][START_REF] Donald | Pseudo-second-derivative matrix and its application to automatic lens design[END_REF][START_REF] William | Numerical Recipes in C (2nd Ed.): The Art of Scientific Computing[END_REF], but other algorithms can be used or implemented [START_REF] Rogers | A new method for the optimization of optical systems: Comparisons and discussions[END_REF][START_REF] Houllier | Comparing optimization algorithms for conventional and freeform optical design[END_REF][START_REF] Furkan | Open-source optimization algorithms for optical design[END_REF]. Since the Levenberg-Marquardt algorithm is a local minimization algorithm, convergence to a global minimum is not guaranteed. The expertise of a lens designer is therefore necessary, both to choose a suitable starting point for the optimization [START_REF] Côté | Deep learningenabled framework for automatic lens design starting point generation[END_REF] and to steer it towards an acceptable local minimum. Here, the starting point is part of the CodeV example library. The triplet is optimized in the CodeV conventional way by setting as variable all element surface curvatures, element center thicknesses, element distances and also indices (n) and Abbe numbers (V ) of glass materials, letting general constrains (see Supplemental 1 Section 1) warrant that the system is physically realistic. Once the optimization is done, true existing materials close to the continuously optimized (n, V ) values are chosen from the CodeV glass materials library, among the inexpensive and common glasses from the Schott and Ohara catalogs, and the system is re-optimized, resulting in the lens shown in Fig. 1(a) and precisely described in Supplemental 1 (section 2), referred to as the "conventional triplet". In order to evaluate the performance of a hybrid system based on this triplet, we will consider several evaluation metrics. In Sec. 3.2, we analyze the MTFs at different positions in the FoV. It is a purely optical metric, useful to compare the raw image qualities on the sensor plane given by the different Cooke triplets.

In order to take into account the digital post-processing, we also analyze the effective MTFs (a metric that takes into account the optical system and the postprocessing, which will be defined in the next paragraph). Finally, in Section 3.3, we evaluate the global system performance through image simulation.

Analysis of the MTFs and effective MTFs

Figure 2(a) represents the tangential and sagittal MTFs of the conventional triplet at different positions in the FoV. It has the standard characteristics of a lens of this type: its structure allows a good correction of aberrations, which makes possible to obtain, on the axis, a MTF taking fairly high values (higher than 0.3 at the Nyquist frequency). However, performance degrades as the FoV increases: at 14 • half-FoV, the sagittal MTF presents nullings at the spatial frequency 40 lp.mm -1 . This means that at (and around) this frequency, for this particular FoV position and orientation, the signal is drowned out in noise and cannot be recovered. The relative frequency-wise noise level (i.e. S nn (ν)/S oo (ν)) corresponding to a SNR of 34 dB is shown with a black dotted line in Fig. 2(a) . However, the MTFs alone are not representative of the performance of the whole hybrid optical/post-processing system, since they do not take into account the image restoration with post-processing. A more representative metric of the entire hybrid system is the effective MTF defined as:

| heff ψ k (ν)| = | hψ k (ν) × wΨ (ν)| . (7)
In an ideal hybrid system, the effective MTF should be uniformly equal to 1, meaning that the signal components at all spatial frequencies are perfectly restored by the deconvolution filter. Figure 2(b) represents the effective MTFs of the conventional triplet. We can observe that for the highest frequencies (between 85 lp.mm -1 and 100 lp.mm -1 ), the level of all the effective MTFs drop compared to the corresponding MTFs. This is due to the fact that the Wiener filter makes a compromise (the best in the sense of the MSE) between signal reconstruction and deleterious noise reinforcement. On the contrary, for spatial frequencies lower than 45 lp.mm -1 , the effective MTF on axis is enhanced and exceeds 1, which generates, as we will see later, visual over-contrast. This is due to the fact that the on-axis MTF is higher than the "average MTF" (Eq. 2) to which the average Wiener filter is adapted. In consequence, this filter "over-compensates" the on-axis MTF. Furthermore, the MTF nullings are still present on the effective MTFs, since information is totally lost at these spatial frequencies and cannot be recovered by the average Wiener filter.

Analysis of final image quality

Let us now verify the observations made on MTFs and effective MTFs through image simulations. Figure 3 shows the scene used for these simulations. It is a large image, of 5144 × 5144 pixels, chosen to obtain a FoV of 40 • (±20 • ) on the diagonal of the sensor, with a pixel size close to 5 µm. Several subparts are selected on the raw image, to check the reconstruction of the details at three different positions in the FoV: on axis, intermediary position and extremity. Simulation of the observation of this scene through the Cooke triplet is performed with the image simulation tool of the CodeV software, with 20 × 20 variable PSFs in the FoV, in order to take into account all the defects of the optical system. The image is then corrupted with a white Gaussian noise so that the SNR is 34 dB. For the conventional triplet used without deconvolution, the resulting subpart images are shown in the first row of Fig. 4. One observes that the images appear blurred, especially those corresponding to large fields. Noise appears quite low, thanks to the high value of the SNR (34 dB). Figure 4(b) represents the same subparts after post-processing with an average Wiener filter constructed from OTFs corresponding to thirteen different positions in the FoV (see Supplemental 1 section 3 for more details). Deconvolution raises the contrast reduced by the optical system, and this can be shown theoretically with the mean image quality (IQ), expressed in dB and defined as

IQ mean = 10 log 10 1 M SE mean (Ψ) . (8)
It is a decreasing function of M SE mean defined in Eq. 6. It thus increases when the performance is enhanced. Table 1 contains the values of IQ mean for all the Cooke triplets optimized in this article. For the conventional triplet, the IQ mean goes from 12.9 dB up to 13.5 dB thanks to deconvolution. Nevertheless, this increase is not homogeneous across the FoV: the on-axis subpart image (Fig. 4(d)) is more contrasted than it should be (according to the ideal image in Fig. 3). The deconvolved middle field image (Fig. 4(e)) appears to be well reconstructed, while the deconvolved extreme field image (Fig. 4(f)) is only slightly more contrasted than the unprocessed subpart (Fig. 4(c)). These observations match the conclusions drawn from the effective MTFs. Indeed, Fig. 2(b), an effective MTF higher than 1 for low spatial frequencies indicates over contrast, while an effective MTF close to 0 indicates under contrast (blur).

4 Surrogate criterion based on optical metrics

Optimization method

Since our optimization goal is to obtain the best possible image quality over the whole FoV, it is possible to define the "best system" as having equivalent performance at any point of the FoV, after deconvolution. Since we have decided to use a single deconvolution filter, this is equivalent to optimizing the optical system under the two following general constraints.

• The MTFs at all the positions in the FoV should be nearly equal to each other, so that deconvolution by a unique filter yields equivalent performance across the FoV. To reach this goal, one should implement specific optimization constraints forcing the MTFs at all FoV positions to be as close as possible to each other for a few chosen spatial frequencies.

• The MTF values should be maximized to obtain the best possible frequencywise SNR before deconvolution (and avoid any MTF nullings in the spatial frequency range of interest). This approach will be called "co-optimization by MTFs equalization". It should be noted that the maximization of the MTF values is a very important point of this optimization technique. In addition to corresponding well to the goal of obtaining the best possible image quality over the whole FoV, this co-design criterion can be seen as an alternative (a "surrogate") to the MSE. Indeed, several studies have shown that optimizing performance in terms of MSE over a set of parameters (e.g. depths of field) leads to MTFs that are close to each other [START_REF] Rafael Falcón | Performance limits of binary annular phase masks codesigned for depth-of-field extension[END_REF][START_REF] Burcklen | Joint digital-optical design of complex lenses using a surrogate image quality criterion adapted to commercial optical design software[END_REF][START_REF] Fontbonne | Theoretical and experimental analysis of co-designed binary phase masks for enhancing the depth of field of panchromatic cameras[END_REF]. This criterion can be implemented in various ways in commercial optical software [START_REF] Burcklen | Joint digital-optical design of complex lenses using a surrogate image quality criterion adapted to commercial optical design software[END_REF]. In this work, we have disabled the default CodeV error function based on the sum of squared spot diagram RMS diameters and only enforced the following user-defined constraints:

• Favoring equality between the MTFs at a set of some Y-meridional FoV positions and the on-axis MTF (it is a set of weighted constraints, to allow some latitude and balancing on these equalizations). In order to guarantee the similarity of the MTFs whatever the azimuth, one has to take into account the tangential, sagittal and 45 • MTFs.

• Lower bounding the value of the on-axis MTF (it is an inequality constraint). This constraint acts somewhat like the minimization of the spot diagram diameters in order to get good image quality.

In practice, it is sufficient to enforce similarity of the MTFs at only two or three well chosen spatial frequencies. This method is very simple to implement in an optimization routine and the modifications of the weights of each constraint can be done in an intuitive manner to drive the whole optimization process on its way.

The conventional triplet is used as a starting point. During the optimization process, we use the same variables and the same general constraints as in the former conventional optimization in Sec. 3 (See Supplemental 1, section 1). We also use the procedure previously explained to switch from fictitious (n,V) glasses to real catalog materials. This results in the triplet shown in Fig. 1(b), whose precise characteristics are given in Supplemental 1 (section 2). It will be referred to as "EMTF triplet". On a structural point of view, this new configuration has a shape similar to the conventional triplet. The main difference is that the central biconcave lens has become larger and its material has changed (NSF8 instead of NSF10). No MTF (in any field position, and neither in tangential nor in sagittal) presents nullings. On the other hand, they are, on average, lower than in the conventional case. For example, the on-axis MTF, which is quite high in the conventional case, takes values below 0.2 at a spatial frequency of 40 lp.mm -1 , whereas the Nyquist frequency is 100 lp.mm -1 . It goes below the relative noise curve at 60 lp.mm -1 . The performance on the axis is therefore lower than for the conventional system. This difference in optical performance is confirmed by the RMS diameter of the spot diagrams, which are much larger -by a factor of about 5 -in the codesigned case (Fig. 5(b)) than in the conventional case (Fig. 5(a)). The shape of the spot diagrams Fig. 5(b) is characteristic of the presence of strong aperture aberrations (i.e. spherical aberration), which means that the MTF equalization co-optimization method deliberately added aperture aberrations to the original system, most probably to be able to decrease the field-variant aberrations, but also in a special way that does not introduce early MTF nullings.

Analysis of MTFs and effective MTFs

The study of this EMTF triplet would not be complete without considering the digital processing by the average Wiener filter for which it was (implicitly) co-optimized. Fig. 2(d), the effective MTFs obtained after deconvolution are, whatever the field position, close to 1 up to the spatial frequency 40 lp.mm -1 . This means that the reconstruction will be possible for all patterns with frequency below 40 lp.mm -1 , regardless of their position within the FoV. This was not the case for the conventional triplet, where the effective MTFs (Fig. 2(b)) took very heterogeneous values, reaching zero at some positions in the field. The effective MTFs of the EMTF triplet smoothly decrease only after 40 lp.mm -1 , due to the necessary trade-off made by the average Wiener filter between signal and noise.

Analysis of final image Quality

To compare the performance of the hybrid systems based on optimized triplets, we will use the IQ mean (Eq. 8), Tab. 1. It is seen that without deconvolution, the EMTF triplet yields a IQ mean 2 dB lower than the conventional triplet. After deconvolution, the performance of both triplets is identical: the performance gain by deconvolution is thus very important for the EMTF triplet (2.6 dB) and exactly compensates for the drop observed before deconvolution. This large performance gain after deconvolution can be easily seen by comparing the first row (without deconvolution) to the second (with deconvolution) of Fig. 6. In the first case, the images appear extremely blurred, much more than in the conventional case (Fig. 6(a-b-c)). On the other hand, the details appear well restored after deconvolution (Fig. 6(d-e-f)).

In order to facilitate the comparison between the hybrid optical/processing systems, two details of the image are precisely analyzed Fig. 7. They belong to the on-axis subpart image and to the extreme field subpart image. The observation of these two details provides a good understanding of how the differences in MTFs play into the imaging performance. We have already observed that the conventional triplet shows an over contrast in the on-axis position (Fig. 7(c)) compared to the ideal scene image (Fig. 7(a)). This over contrast is due to the excessively high value of the on-axis effective MTF at low spatial frequencies.

On the other hand, the effective MTFs for larger fields are relatively low, which means that the reconstruction is poor. In particular Fig. 7(d), the pinnacle ornaments are barely visible. Conversely, the MTFs of the EMTF triplet are very close to each other, which explains that imaging performance is similar for the on-axis detail (Fig. 7(e)) and the extreme field detail (Fig. 7(f)): the contrast is less important on axis than for the conventional triplet, but the reconstruction is better at the extremity of the FoV, where the ornaments of the pinnacle are clearly visible. It can also be noted that with the EMTF triplet, the deconvolved image is slightly more noisy than with the conventional one: this can be seen with the slight "orange skin" effect in the sky, Fig. 7(f). This observation suggests an increased sensitivity to noise of the hybrid system based on the EMTF triplet. This was expected since the MTFs of the EMTF triplet are globally lower than those of the conventional one, so that the frequency-wise SNR before deconvolution is smaller. However, since the MTFs of the EMTF triplet remain globally higher than the noise level, deconvolution is still effective, even in the extreme field (Fig. 7(f)): the noise is increased (and at the origin of the "orange skin" effect), but the signal is clearly enhanced. On the contrary, for the hybrid system based on the conventional triplet (Fig. 7(d)), the noise is not increased, so the sky remains smooth, but the signal is not well reconstructed, as the effective MTF is too low.

To conclude, co-optimization of the Cooke triplet by MTF equalization leads to similar average performance as the conventional triplet, after deconvolution. However, the distribution of the performance over the FoV is different: while the conventional method favors the positions in the field close to the axis at the expense of the extreme fields, the EMTF method homogenizes the performance, as originally intended. 5 The MSE as a co-design criterion

Implementation and use of the MSE criterion

So far, we have used the MSE to evaluate the performance of hybrid systems (through the IQ criterion), but not as a direct co-optimization criterion. An implementation of this criterion has already been done under Zemax [START_REF] Stork | Theoretical foundations for joint digital-optical analysis of electro-optical imaging systems[END_REF][START_REF] Wang | An integral design strategy combining optical system and image processing to obtain high resolution images[END_REF][START_REF] Vettenburg | Holistic optical-digital hybridimaging design:wide-field reflective imaging[END_REF], but until recently, it seemed to us impossible to do the same under the other well-known CodeV commercial lens design software, since it was not possible to directly access PSFs (or 2D OTFs) in a user-defined error function during optimization. Using the new features of CodeV version ≥ 11.0 (2017) that remove these limitations, we have been able to implement -for the first time to our knowledge -MSE optimization directly into CodeV. The details of the implementation are given in Supplemental 1 (section 3). In this section, we cooptimize a Cooke triplet with this method and compare its performance with the two previously described triplets. This allows a fair comparison since all these triplets have been optimized with the the same lens design software. As for any optimization, the final solution, which corresponds to a local minimum of the criterion, may depend on the chosen starting point. We have thus considered two different starting points. The "MMSE 1 triplet", represented on Fig. 1(c), has been obtained by choosing the conventional triplet as a starting point. The "MMSE 2 triplet", displayed on Fig. 1(d), has been obtained by starting from the EMTF triplet. The detailed characteristics of these two triplets a given in Supplemental 1 (section 2). By comparing Fig. 1(c) and Fig. 1(a) on the one hand, and Fig. 1(b) and Fig. 1(d) on the other hand, it is clear that the MMSE triplets have kept the main structural features of their starting points. However, we will see in the following that its imaging performance is different.

Analysis of MTFs and effective MTFs

Figure 2(e) represents the MTFs of the MMSE 1 triplet. We can observe that the on-axis MTF collapses quickly and shows a nulling at 50 lp.mm -1 . On the other hand, the other MTFs no longer show any nulling for low spatial frequencies, which is a significant improvement over the conventional system. As a result, Fig. 2(f), the effective MTFs higher than 0.38 for spatial frequencies under 40 lp.mm -1 , which was not the case in the conventional case (Fig. 2(b)). However, the on-axis effective MTF presents a nulling (the same as for the onaxis MTF). The optimization of M SE mean (see Eq. 6) thus clearly gave priority to the peripheral positions in the FoV over the on-axis MTF, whereas the onaxis MTF is more often favored in conventional optical design because there is naturally less aberration around the optical axis. This evolution was possible because in the MSE-based criterion, the on-axis performance does not have a larger weight than any other field. Figure 2(g) displays the MTFs of the MMSE 2 triplet. The MTFs are close to each other, as for the EMTF triplet. Similarly, the effective MTFs are close to each other (Fig. 2(h)). The MMSE 2 triplet is the only system that has high effective MTFs over a wide range of spatial frequencies, with most of its effective MTFs higher than 0.4 at 60 lp.mm -1 . In terms of effective MTFs, it can thus be considered as the best solution.

Analysis of final image quality

Let us now analyze the obtained image quality. It is seen in Table 1 that the MMSE 1 triplet and the MMSE 2 triplet both yield enhanced image quality after deconvolution. This was expected since they are optimized to have a smaller M SE mean , but the important conclusion is that the performance gain is significant. Indeed, for the MMSE 1 triplet, the image quality after deconvolution is increased by 1.1 dB. This means that the optimization on the MSE, by tightening together the previously disparate MTFs (Fig. 2(a) and Fig. 2(e)), made the deconvolution by the average Wiener filter more efficient. Similarly, the MMSE 2 triplet yields a lower image quality than its starting point (the EMTF triplet) before deconvolution, but a much better one after deconvolution. It shows the best average performance after deconvolution, reaching IQ mean = 16 dB. The performance gain with respect to the starting point is even larger in this case, since it amounts to +2.5 dB. Such a large increase in IQ mean can be explained by the fact that the effective MTFs are close to 1 for a wider range of spatial frequencies. This increased performance of the MMSE 2 triplet is confirmed in the simulated images, Fig. 8. Contrary to the conventional triplet, no over contrast is visible on the restored on-axis subpart image (Fig. 8(d)). The middle field (Fig. 8(e)) and extreme field (Fig. 8(f)) subpart images also appear very well restored. It is possible to observe the reconstruction more in details on Fig. 7(g-h). Compared to the other hybrid systems, we notice that when using MMSE 2 triplet, the restored image is slightly more blurred on the axis, but much sharper in the field: it is the method that best reconstructs the pinnacle ornaments. This behavior was predictable from our analysis of the effective MTFs. These results lead to several interesting conclusions. The main one is that the hybrid system based on the MMSE 2 triplet has the best performance among the systems considered in this article, both in terms of quantitative and visual image quality. This clearly demonstrates the benefit of implementing MSEbased end-to-end optimization directly in a commercial lens design software. The second conclusion is that both MMSE 1 and MMSE 2 triplets have kept the characteristics of their starting points in terms of optical architecture and properties (heterogeneity or homogeneity of MTFs). Since the MMSE 2 triplet performs better, this shows that, to reach our optimization goal (good image quality at all FoV positions), it was preferable to choose as a starting point a lens having already similar MTFs over the FoV. However, we can also suspect that MSE optimization, as implemented, may have difficulties to "jump out" local minima. The sensitivity of lens design to starting point is a long-standing problem, and a lot of efforts have been done to bring solutions to it in the framework of conventional optical design. An important perspective to the present work is to devote the same efforts to enhance the capacity of MSEbased co-optimization to find globally optimal solutions.

Conclusion

In this article, we have compared three different methods to co-optimize a hybrid optical/digital imaging system with commercial lens design software: conventional optimization based on minimization of the RMS diameter of spot diagrams, a surrogate criterion based on near equality of the MTFs and a true MSE criterion taking explicitly into account the digital processing in the optimization process. To implement the latter method, we integrated -for the first time to our knowledge -MSE optimization to the CodeV software. These three methods have been illustrated and compared on a concrete application: the design of a Cooke triplet having good image quality everywhere in the FoV. The obtained results demonstrate the superiority of MSE co-optimization over the other methods, both in terms of quantitative and visual image quality. Analyzing in detail the different resulting lenses, we have shown that the conventional triplet tends to favor the center of the FoV at the expense of the extreme positions, whereas the hybrid system based on the MSE (the MMSE 2 triplet) yields much more homogeneous performance thanks to co-optimization with deconvolution. This homogeneity comes at the expense of the performance on axis, but enable much better performance at peripheral FoV positions. This shows that it is possible, by leveraging deconvolution during the optimization process, to adapt the spatial distribution of imaging performance to a prescribed goal-which was, here, to reach good performance over the whole FoV. This work has many perspectives. First, we assumed that the deconvolution algorithm was a unique linear filter common to all FoV positions. This postprocessing method is clearly not adapted to the conventional triplet, whose MTFs are quite inhomogeneous over the FoV. They are much more homogeneous with the systems resulting from the two other methods, but there is still residual inhomogeneity. To take into account the inhomogeneity of the MTFs over the FoV and improve further the performance, it would be interesting to consider spatially variant post-processing algorithms [START_REF] Denis | Fast approximations of shift-variant blur[END_REF] in the co-optimization process. Furthermore, recent works have proposed to use home-made differentiable ray tracing connected with a deep neural-network [START_REF] Vincent Sitzmann | End-toend optimization of optics and image processing for achromatic extended depth of field and super-resolution imaging[END_REF][START_REF] Shang | Two-step training deep learning framework for computational imaging without physics priors[END_REF][START_REF] Nehme | Learning optimal wavefront shaping for multi-channel imaging[END_REF][START_REF] Sun | End-to-end complex lens design with differentiable ray tracing[END_REF][START_REF] Liu | End-to-end computational optics with a singlet lens for large depthof-field imaging[END_REF][START_REF] Halé | End-to-end sensor and neural network design using differential ray tracing[END_REF].

Comparison with the approach proposed in this paper, that uses a powerful commercial software and a simpler (and faster) post-processing algorithm, would be very interesting. Moreover, it would be interesting to improve the MSE optimization method, by improving for example its speed of execution, its resilience to the choice of the starting point, and also to study the effect of the SNR parameter, or optimize, using the multi-configuration possibility of the lens design program, an optical system that should optimally work for different apertures (e.g. variable iris diaphragm) and therefore different SNRs.

A Supplemental 1

This supplemental document presents the optical constraints used for the cooptimization of the Cooke triplets presented in this article (Sec. A.1), the detailed characteristics of the Cooke triplets co-optimized with various methods (Sec. A.2), and the MSE optimization (Sec. A.3).

A.1 Constraints for optimization of the Cooke triplets

The general constraints applied to the four Cooke triplets are given in Table 2.

For the (fictitious) glasses, which can continuously evolve during optimization in the (n, V ) or (n, ∆n) map, where n is the refractive index, V the Abbe number and ∆n the blue-red dispersion with V = (n -1)/∆n, we use the default CodeV constraints that enforce the glass coordinates (n, ∆n) to stay inside a realistic glass map zone defined by a convex polygon whose vertices correspond to the Schott glasses N-FK5, N-SK16, N-LAF2, SF4 [START_REF] Code | CODE V Introductory User's Guide[END_REF][START_REF] Code | V Lens System Setup Reference Manual[END_REF][START_REF] Code | V Optimization Reference Manual[END_REF].

Applied constraints

CodeV

A.3 Implementation of MSE-based co-optimization in CodeV

CodeV lens design software offers the ability to write advanced optimization macros, that may use externally coded efficient compiled extension in the form of Dynamic Link Libraries (DLLs). It also offers, since version 11.0, the possibility to calculate the point spread function (PSF) of a given field during the iterations of its internal optimization algorithm, thanks to the predefined CodeV Macro-PLUS function PSF DATA 1FLD. Therefore, we create a DLL, based on a code written in C, taking as arguments the inputs (2D arrays describing the PSFs) and the output (the MSE) which is calculated according to the values of the PSFs. Note that computing a too large number of PSFs during an optimization routine would not be practicable due to large computing time. Only thirteen of them have been used during the optimization. The distribution of these thirteen fields, used for the construction of the average Wiener filter useful for the deconvolution of the image obtained with a Cooke triplet of half-field 20 • , is presented in Fig. 9. Their positions guarantee a nearly rotational symmetry of the average Wiener filter, so that it has no preferred direction. In order to compute the value of M SE mean (Ψ) (with here K = 13) using Eq.( 3) of the main paper that expresses it in the Fourier domain, it is useful to use the FFTW free library available in C, which allows to efficiently compute the Discrete Fourier Transform of the PSFs, and thus obtain the optical transfer functions hψ k . This implemented MSE does not take into account defects that do not directly affect the shape of PSFs such as relative intensity change and distortion (which can be otherwise constrained if necessary during optimization). For optimization purposes, these two defects will be considered negligible or easily corrected by post-processing, as they generally have little effect on the visual quality of an image. The written C code must be compiled as a 32-bits DLL (CodeV is only available for this software architecture) using a Synopsyssupplied Makefile, extended to allow FFTW use. This allows one to obtain a DLL (and associated files) that can easily be called in a homemade CodeV macro function. This function can then be called directly in the CodeV optimization routine, like other functions of interest to the optical designer. The optimization goal is then to make the obtained MSE result as close as possible to 0.

Figure 1 :

 1 Figure 1: Scheme of Cooke triplets. (a) "Conventional triplet". (b) "EMTF triplet". (c) "MMSE 1 triplet". (d) "MMSE 2 triplet".

Figure 2 :

 2 Figure 2: (First column) MTFs and (Second column) Effective MTFs of the (ab) conventional triplet, (c-d) EMTF triplet (e-f) MMSE 1 triplet (g-h) MMSE 2 triplet. The curve legends for the all first column plots are the same as the one in (a) and the legends for the all second column plots are the same as the one in (b).

Figure 3 :

 3 Figure 3: Image of the Gloucester Cathedral (public domain picture) used as a scene for the comparison of Cooke triplets optimized with different methods. Three images are selected for their details and their positions in the FoV (on axis, middle field, extreme field).

Figure 4 :

 4 Figure 4: Details of the cathedral image obtained with the conventional triplet for an SNR of 34 dB, (a-b-c) without deconvolution and (d-e-f) with deconvolution.IQ meanConventional triplet

Figure 2 (

 2 Figure2(c) shows the MTFs obtained with the EMTF triplet. Contrary to the very disparate MTFs of the conventional triplet (Fig.2(a)), these MTF curves are very close to each other. No MTF (in any field position, and neither in tangential nor in sagittal) presents nullings. On the other hand, they are, on average, lower than in the conventional case. For example, the on-axis MTF, which is quite high in the conventional case, takes values below 0.2 at a spatial frequency of 40 lp.mm -1 , whereas the Nyquist frequency is 100 lp.mm -1 . It goes below the relative noise curve at 60 lp.mm -1 . The performance on the axis is therefore lower than for the conventional system. This difference in optical performance is confirmed by the RMS diameter of the spot diagrams, which are much larger -by a factor of about 5 -in the codesigned case (Fig.5(b)) than in the conventional case (Fig.5(a)). The shape of the spot diagrams Fig.5(b) is characteristic of the presence of strong aperture aberrations (i.e. spherical aberration), which means that the MTF equalization co-optimization method deliberately added aperture aberrations to the original system, most probably to be able to decrease the field-variant aberrations, but also in a special way that does not introduce early MTF nullings. The study of this EMTF triplet would not be complete without considering the digital processing by the average Wiener filter for which it was (implicitly) co-optimized. Fig.2(d), the effective MTFs obtained after deconvolution are, whatever the field position, close to 1 up to the spatial frequency 40 lp.mm -1 . This means that the reconstruction will be possible for all patterns with frequency below 40 lp.mm -1 , regardless of their position within the FoV. This was not the case for the conventional triplet, where the effective MTFs (Fig.2(b)) took very heterogeneous values, reaching zero at some positions in the field. The effective MTFs of the EMTF triplet smoothly decrease only after 40 lp.mm -1 , due to the necessary trade-off made by the average Wiener filter between signal and noise.

Figure 5 :

 5 Figure 5: Spot diagrams of (a) the conventional triplet and (b) the EMTF triplet.

Figure 6 :

 6 Figure 6: Details of the cathedral image obtained with the EMTF triplet for an SNR of 34 dB, (a-b-c) without deconvolution and (d-e-f) with deconvolution.

Figure 7 :

 7 Figure 7: On axis (first column) and the extreme field (second column) details of the cathedral scene. (a-b) Ideal images. (c-d) Conventional triplet and deconvolution. (e-f) EMTF triplet and deconvolution. (g-h) MMSE 2 triplet and deconvolution.

Figure 8 :

 8 Figure 8: Details of the cathedral image obtained with the MMSE 2 triplet for an SNR of 34 dB, (a-b-c) without deconvolution and (d-e-f) with deconvolution.

Figure 9 :

 9 Figure 9: Representative scheme of the K = 13 positions in the field considered for the construction of the average Wiener filter useful for the deconvolution of the image obtained with an optical system of half-field 20 • .

Table 1 :

 1 IQ mean before and after deconvolution for the different Cooke triplets.

Table 2 :

 2 Constraints applied on Cooke triplets (Strict inequality constrains for all relevant lens elements enforced by the Lagrange multipliers method).

	mnemonic Length [mm]

Acknowledgments

The work reported in this study is supported in part by the Agence de l'Innovation de Défense (AID) that provides half of a PhD fellowship to Alice Fontbonne.

Disclosures The authors declare no conflicts of interest.

Data Availability Statement Data is available upon request.

A.2 Characteristics of the co-optimized Cooke triplets

We give here the characteristics of the different Cooke triples co-optimized in this article: the conventional triplet (Tab. 3), the EMTF triplet (Tab. 4), the MMSE 1 triplet (Tab.


Bukeikhan Omarali
email: b.omarali@qmul.ac.uk

Francesca Palermo

Kaspar Althoefer

Maurizio Valle

Ildar Farkhatdinov
email: i.farkhatdinov@qmul.ac.uk

Tactile Classification of Object Materials for Virtual Reality based Robot Teleoperation

This work presents a method for tactile classification of materials for virtual reality (VR) based robot teleoperation. In our system, a human-operator uses a remotely controlled robot-manipulator with an optical fibre-based tactile and proximity sensor to scan surfaces of objects in a remote environment. Tactile and proximity data and the robot's endeffector state feedback are used for the classification of objects' materials which are then visualized in the VR reconstruction of the remote environment for each object. Machine learning techniques such as random forest, convolutional neural and multimodal convolutional neural networks were used for material classification. The proposed system and methods were tested with five different materials and classification accuracy of 90% and more was achieved. The results of material classification were successfully exploited for visualising the remote scene in the VR interface to provide more information to the humanoperator.

I. INTRODUCTION

The last few years have seen rapid adoption of Virtual Reality (VR) technologies in robot teleoperation. The key advantage of VR-based teleoperation interfaces is the threedimensional (3D) reconstruction of a remotely operated robot and its' environment that improves a human-operator's (further referred to as "operator") sense of presence [START_REF] Whitney | ROS Reality: A Virtual Reality Framework Using Consumer-Grade Hardware for ROS-Enabled Robots[END_REF], [START_REF] Omarali | Virtual reality based telerobotics framework with depth cameras[END_REF], [START_REF] Kohn | Towards a Real-Time Environment Reconstruction for VR-Based Teleoperation Through Model Segmentation[END_REF], [START_REF] Omarali | Workspace scaling and rate mode control for virtual reality based robot teleoperation[END_REF]. However, this also leads to operators mainly relying on visual feedback. An important aspect of teleoperation is providing the operator with reliable information about the objects' materials in a remote environment. For example, the operator may be required to collect all metallic objects in the remote environment which might be difficult if the operator can only rely on the visual reconstruction of the remote environment, which is typically implemented with a pointcloud geometric visualisation based on the data acquired from RGB-D cameras. Often such cameras acquire noisy data that reduce the accuracy of VR reconstruction.

We suggest using tactile sensing and classification to identify materials of objects in the remote environment and visually communicate those materials to the operator in the VR interface. The majority of objects' materials classification techniques are based on machine learning applied to images of objects [START_REF] Cimpoi | Deep filter banks for texture recognition and segmentation[END_REF], [START_REF] Xue | Differential angular imaging for material recognition[END_REF]. A convolutional neural network (CNN) architecture was used to classify materials from patches extracted from photos of objects, followed by localising and segmenting the materials in the original images [START_REF] Bell | Material recognition in the wild with the materials in context database[END_REF]. In situations with reduced lighting when computer vision methods can fail (for example teleoperation in extreme and hazardous environments [START_REF] Vitanov | A suite of robotic solutions for nuclear waste decommissioning[END_REF]), tactile exploration can be used for material recognition [START_REF] Boissieu | Tactile texture recognition with a 3-axial force mems integrated artificial finger[END_REF], [START_REF] Drimus | Object texture recognition by dynamic tactile sensing using active exploration[END_REF]. A combination of proximity, tactile and force sensing can provide important information on the explored material such as texture, stiffness and shape [START_REF] Tomo | Design and characterization of a three-axis hall effectbased soft skin sensor[END_REF]. Compliance properties of different objects via supervised classifiers using a hybrid force and proximity finger-shaped sensor was investigated in [START_REF] Konstantinova | Object classification using hybrid fiber optical force/proximity sensor[END_REF]. A fibre opticbased sensor has been designed and used in [START_REF] Palermo | Implementing tactile and proximity sensing for crack detection[END_REF], [START_REF] Palermo | Multi-modal robotic visual-tactile localisation and detection of surface cracks[END_REF], [START_REF] Palermo | Automatic fracture characterization using tactile and proximity optical sensing[END_REF] to recognise and classify fractures on surfaces implementing a random forest classifier.

In contrast to machine learning algorithms which need engineered features extracted from the data, deep learning models achieve higher performance when trained with high dimensional input, such as audio in speech recognition [START_REF] Hinton | Deep neural networks for acoustic modeling in speech recognition: The shared views of four research groups[END_REF] and figures in image recognition and detection [START_REF] Krizhevsky | Imagenet classification with deep convolutional neural networks[END_REF]. Various machine learning classifiers (Gaussian, K-nearest neighbours, support vector machine) were compared in [START_REF] Baishya | Robust material classification with a tactile skin using deep learning[END_REF] against a CNN model to classify explored surfaces via tactile skin sensor attached to an iCub humanoid robot. It was found that the deep learning model for tactile material classification provides more robust performance in comparison to feature learning-based methods. The authors of [START_REF] Gao | Deep learning for tactile understanding from visual and haptic data[END_REF] demonstrated that implementing a multi-modal approach based on CNN with both visual and physical interaction signals achieved more accurate results than vision only as well as superior classification compared to hand-designed features methods. Authors of [START_REF] Alameh | Dcnn for tactile sensory data classification based on transfer learning[END_REF] proposed an algorithm for recognising the object touched via human interactions on an electronic skin. The 3D tactile data generated by the skin were converted into 2D images and used as input for a CNN which outperformed classical tactile data classification algorithms.

We propose a method for tactile classification of materials for VR-based robot teleoperation. In our system, the operator remotely controls a robotic arm with an optical fibre-based tactile and proximity sensor to scan surfaces of objects in a remote environment. Tactile and proximity data as well as the robot's end-effector state feedback are used for the classification of objects' materials. Classification results are then used to visualize objects' materials in VR. Machine learning techniques such as random forest, convolutional neural and multi-modal convolutional neural networks were used for material classification.

To the authors' knowledge, this work is the first attempt to demonstrate the integration of material classification with tactile sensing in VR-based robot teleoperation. Our work proposes a novel material classification technique based on data collected with a fibre optics-based tactile sensor, humansupervisory control and visualisation of objects' materials in VR.

II. SYSTEM OVERVIEW

Key components of the system are shown in Fig. 1: Franka Emika's Panda robot (a), a fibre optics tactile sensor (b), an Intel Realsense 435i RGBD camera, Oculus Rift S headset and Oculus Touch handheld controllers (d). Two desktop computers were connected in a wired network: the operator's computer running the VR interface (built with Unity 3D), and the robot's computer running ROS. The robot was controlled by the operator using a VR-based robot teleoperation interface described in our previous work [START_REF] Omarali | Virtual reality based telerobotics framework with depth cameras[END_REF]. The operator sets desired motion path by manually placing waypoints in a VR reconstruction of the remote environment. These waypoints are then used to plan the robot's motion, which is previewed and accepted or rejected by the operator. The remote environment was visualized in VR using a point cloud from an Intel Realsense 435i RGB-D camera mounted on the robot's end-effector.

We used an integrated force and proximity finger-shaped sensor described in [START_REF] Konstantinova | Object classification using hybrid fiber optical force/proximity sensor[END_REF], [START_REF] Konstantinova | Fingertip proximity sensor with realtime visual-based calibration[END_REF] attached to the robot's endeffector. The sensor consists of 3D printed rigid and soft components that allow the finger to bend during interaction with the environment, as shown in Fig. 1 use reflected lights' intensity to measure the deformation of the finger. The fourth pair of optical fibre cables (P) is used to measure the distance between the tip of the finger and nearby objects. Each pair of the sensor's fibre optic cables was attached to a Keyence FS-N11MN lightto-voltage transducers that communicate with the robot's control computer at 400 Hz.

III. METHODS Tactile Data Collection. Supervised teleoperation was used for tactile exploration of five objects made of metal, paper, silicon, hard styrofoam, soft foam. All objects were flat and placed horizontally in the robot's workspace. Fig. 1(a) shows sample objects made of metal, styrofoam and soft foam. An expert operator was asked to set the reference tactile scan pose, which determined the location and the orientation of initial touch between the end-effector mounted tactile sensor and the object as well as the subsequent sliding (scanning motion). Setting the initial conditions for tactile exploration manually allows efficient collection of highly variable data representative of real-world use cases suitable for the classifiers' training and validation.

Objects were scanned using horizontal linear sliding motion that was generated from a reference pose set by the operator. The operator set the reference poses by manually placing an axis 3D mesh object in the VR teleoperation interface, see Fig. 2(a). The reference pose was then used to generate the scan start and reverse positions (both generated at the same height as the reference pose) for the tactile sensor. The slide was always performed along the sensor's y-axis projection on the object's horizontal plane, as it results in more distinct D1 and D2 sensor outputs.

The orientation of the tactile sensor with respect to the scanned object was determined by the orientation of the 3D axis mesh set by the operator in VR. The robot has always approached to and moved away from the scanned object along the sensor's z-axis. Different orientations of the sensor and corresponding approaches have resulted in different sensor deformation behaviours. The contact force exerted by the robot/finger on the scanned object (and by extent the penetration depth on softer materials) depended on the reference pose height set by the operator.

After the operator set the reference pose, the robot planned and executed a trajectory necessary to approach and get in contact with the scanned object, see Fig. 2(b,c), performed a slide Fig. 2(d), reversed the slide direction, slid back Fig. 2(e) and retracted. The length of the sliding motion was set to 56.5 mm. The motion plans were generated using Moveit and the planned trajectory's scan duration was automatically re-scaled to 8.3 s.

We have recorded raw outputs of the tactile sensor for object classification. We recorded the robot's end-effector's average position error during the scan as the mean of differences between the end-effector's desired and actual positions. The robot's position error can be used as an indicator of objects' softness. For example, the tactile sensor can deform and penetrate into soft foam resulting in a small position error, which is not the case for metal. We also recorded the reference pose orientation. The tactile sensor scan output varied depending on the sensor's orientation with respect to the scanned objects, hence desired orientation was used as a feature in classifiers.

The operator set the reference pose and the robot performed the scan three separate times. Then the operator changed the reference pose by either changing the reference orientation, or reference position on the scanned object's surface, or the reference scan depth. We have recorded 150 samples (scans) per class.

Data pre-processing for classification. A sample was recorded from 2 seconds before the slide start (to include the initial contact between the sensor and the scanned object) to 2 seconds after the slide end (to include the retraction of the sensor from the scanned object). We used spectrograms of raw tactile data for objects' material classification. Spectrograms were generated with 52 time segments and 50 spectral bands (i.e. each channel's spectrogram is a 52×50 matrix). We compared three classifiers: random forest, convolutional neural network (CNN) and multi-modal convolutional neural network (M-CNN). Random forest was made using scikitlearn library, CNNs were made using Keras with Tensorflow. Spectrograms were standardized using the maximum amplitude present in the dataset. The average end-effector position error was standardized using the maximum absolute value present in the dataset. The reference orientation was represented as a unit quaternion. We generated additional 50 synthetic samples per class by randomly copying existing samples and adding corresponding ±5% standard deviation to each of the spectrograms' time segment and frequency band (to each cell of the spectrogram's 52×50 matrix).

The dataset was split into 60%-20%-20% training, validation and test sets, respectively. We distributed samples into sets based on the tactile sensor's orientation with respect to the scanned objects. This was done to allow the classifiers to train with the data collected at different tactile exploration conditions. Fig. 3 demonstrates sensor orientations in radial angles present in the dataset (note that many samples overlap), where γ is the angle between scanned object's z-axis (scanned surface normal) and the tactile sensor's z-axis and θ is the angle from scanned object's x-axis to y-axis (both are collinear to robot's base frame's x and y-axis). The training set contained samples in upper and lower 30% γ angles, validation and testing sets were randomly chosen from the rest. Hence classifiers were tested and validated on tactile exploration conditions that were not present in the training set, ensuring classifiers' robustness.

Classification with Random forest. Random Forest [START_REF] Breiman | Random forests[END_REF] is a machine learning algorithm used for classification and regression built on an ensemble of multiple learning trees. Insensitive to over-fitting, it can produce reasonable predictions with little tuning and provides an effective way of handling missing data. Random Forest has been vastly used in remote sensing [START_REF] Belgiu | Random forest in remote sensing: A review of applications and future directions[END_REF], [START_REF] Palermo | Automatic fracture characterization using tactile and proximity optical sensing[END_REF], [START_REF] Palermo | Implementing tactile and proximity sensing for crack detection[END_REF]. We implemented a Random Forest classifier with 1000 estimators. The number of estimators was determined by a grid search. The classifier was given concatenated spectrograms, average end-effector position error and reference orientation as inputs.

Classification with Convolutional Neural Network. CNNs are commonly used in image recognition due to their ability to learn cross-correlations between multiple channels (RGB in case of images) and shift invariancy. They can also be used to learn patterns between multiple sensor signals as demonstrated in [START_REF] Raghu | Eeg based multi-class seizure type classification using convolutional neural network and transfer learning[END_REF], [START_REF] Kwon | Electroencephalography based fusion two-dimensional (2d)-convolution neural networks (cnn) model for emotion recognition system[END_REF]. We suggest a similar approach to classify objects' materials -using tactile sensor's multiple sensing channels. The CNN classifier only takes spectrograms as an input and its' model architecture is shown in Fig. 4. There was a dropout with 30% probability between fully connected layers. The output of the last convolutional layer was batch normalized. The model was trained with an early stopping triggered by no improvement on validation loss. We used the Adam optimizer with an exponentially decaying learning rate. The model's and training hyperparameters were determined using grid search.

Classification with Multi-modal Convolutional neural network. We added an extra input to the CNN described above that included the robot's end-effector's average position error and the reference orientation. The model architecture is detailed in Fig. 4. The extra inputs were concatenated with the flattened output of the last convolution layer. The multi-modal CNN retained the 30% probability dropout between fully connected layers and batch normalization after the last convolutional layer. The M-CNN was trained with settings similar to the CNN.

Visualization of object's materials in VR. We propose two visualization methods that communicate scanned objects' predicted classes to the operator in the VR-based robot teleoperation interface. Both methods rely on naive Octomap-based segmentation introduced in our previous work [START_REF] Omarali | Virtual reality based telerobotics framework with depth cameras[END_REF], which allows the operator to segment pointclouds into separate objects. Octomap [START_REF] Hornung | OctoMap: An efficient probabilistic 3D mapping framework based on octrees[END_REF] is a probabilistic mapping method that splits the space (for example the remote environment) into occupied, vacant and unknown cubes (voxels) using pointcloud as an input. Our naive segmentation detects when such occupied cube has been pointed and clicked by the operator and segments all points (of the pointcloud) contained within it as well as within other connected cubes. All segmented points are then stored locally as static segmented pointclouds.

1) Objects' predicted classes as color-coded Octomaps: As the operator sets the reference pose using the 3D axis mesh, the axis mesh checks for collisions with the Octomap using Unity's collider system. These collisions do not have any physical meaning and are simply used to detect which part of the Octomap the scan will be performed on. Once a collision is detected the corresponding Octomap cube and connected cubes are segmented and copied from the "live" Octomap and stored locally. Once the classification is finished the predicted class is used to colour the segmented Octomap according to the color-code given to the operator in the VR interface.

2) Objects' predicted classes as color-coded pointclouds: This visualization method is similar to the one above except instead of segmenting and storing Octomaps, it segments and stores corresponding pointclouds. Once the classification is finished the predicted class is used to color the segmented pointcloud according to the labels' color-codes.

IV. RESULTS

Tactile data. Fig. 5 presents three samples of raw sensor output for each material recorded at different sensor orientations and heights. There is a noticeable variance between different samples of each material. The leftmost vertical lines indicate the start of the slide, the middle one indicates the slide direction reversal, and the last vertical line indicates the stop of the slide.

For all materials the proximity sensor (P) spikes up on the sensor approach and retraction (before the slide begins and after it finishes respectively). The proximity sensor reading depends on the distance to the object. The proximity sensor only detects in a limited range and once the object is too close (or the sensor touches the object) the proximity sensor value drops to the baseline. The amplitude of the spike is also determined by the object's surface roughness, color and reflectivity. In the case of rough and porous materials (soft foam, medium styrofoam), the proximity sensor generates a noisy output even during the touch.

There was a delay between the beginning of the slide and deformation sensors responses. The delay was the time necessary for the sensing finger to deform, (the sensor deformation can be seen in Fig 2d,e). Similarly, responses of the deformation sensors were delayed after slide direction reversal and sensor retraction.

Spectral analysis. Example spectrograms are shown in Fig. 6. Spectrograms were distinct per material, which was noticeable visually. Spectrograms differed most in the 0-15Hz frequency band, where there were large amplitude spikes during the sensor's initial contact with an object, during deformation and sensor retraction. There was also a noticeable high-frequency signal present in rougher/softer materials as well, (see D2 for soft foam).

End-effector position error. The average end-effector position errors occurred along the scanned surface's normal as the robot failed to push in as deep as the operator intended to. Table I presents means and standard deviations of the robot's end-effector's average position error along the scanned surface's normal per material. As expected the harder materials had larger errors: the sensor deformed and pushed in deep into the soft foam but could not do the same with metal. However, T-test showed no statistical significance.

Classification metrics. Accuracy, precision, recall, and f1score classification metrics are used to validate and compare the used classification models. The results for the analysis with the implemented classifiers are presented in Table II. Random Forest, which is robust to outliers and requires little parameter optimisation, achieves the best results. The MCNN, with the extra inputs of the robot's end-effector's average position error and reference orientation, achieves results comparable to the Random Forest classifier while the baseline CNN generates the worst outcome. Thus, even with the small size of the dataset of the network, the extra inputs make the model more robust than the CNN and able to generalise better. Confusion matrices for CNN, M-CNN, and Random forest classifiers are shown in Fig. 7. One of the most challenging materials to classify for all three models is the paperstack, which is frequently confused with the silicon class. This may be due to a similar pattern in the frequency domain.

Predicted class visualization in VR. The Fig. 8 shows the operator's VR view of classification results as colored Octomaps and colored pointclouds as well as the GUI description of the color-coding. The operator can toggle between visualization modes and "live" pointcloud using a button press on a handheld controller. In the current implementation, both methods use static segmented clones of pointcloud and Octomap respectively, which means that updates to "live" pointclouds or Octomap of objects (for example if objects move or RGB-D camera changes position with respect to objects) would not be reflected in classified clones. This presents a technical challenge that requires object tracking as well as continuous segmentation.

V. DISCUSSION AND CONCLUSION Summary. We presented a method for tactile classification of materials for VR-based robot teleoperation. We have collected a dataset using a teleoperated robot. The dataset contains five objects with a wide range of different sensorto-object relative orientations and interaction forces. We have compared three classifiers using spectrograms of raw tactile data and robot's end-effector feedback: random forest, convolutional neural network and multi-modal convolutional neural network. We also introduced two visualisation methods to display the classification results to a human-operator in the VR using color-coded pointclouds and corresponding bounding Octomaps.

The results demonstrated that the random forest classifier had the highest accuracy at 94.5% followed by the multimodal convolutional neural network with 92%. The results of material classification were successfully employed for visualising the remote scene in the VR interface to provide more information to a human-operator.

Limitations. The proposed method currently works only with flat surfaces. In the case of non-flat surfaces, a similar approach can be used where the operator sets multiple waypoints that can be used to generate a non-flat trajectory.

The dataset was collected using constant scan length and duration. It is possible that the random forest classifier trained on this dataset may not be able to generalize well to scans performed at higher or lower speeds. The dataset was collected using teleoperation to be representative of infield operations. This however prevented us from collecting a larger amount of samples that would improve classifier training. Although the dataset contains samples collected at various positions, orientations and depths, it was collected by a single expert operator. It is a small chance that classifier could be overfitting towards this operator's specific mannerisms of setting the desired scan position.

Future work The future work would benefit from larger sample sizes collected by multiple expert operators as well as from an extended object set. The latter would also allow moving from hard-coded material based labels to more general roughness/hardness estimation.

Due to the imprecise nature of pointclouds the operator may put the reference pose too deep into the object potentially damaging the sensor or not deep enough, resulting in poor classification. In future work, we will make better use of the proximity sensor and haptic feedback to place the sensor more accurately at the object's surface.

Proposed methods require a full scan to be completed before the object's class is determined. A potentially interesting follow-up topic is real-time objects' materials classification using Long Short-term Memory networks.

Fig. 1 .

 1 Fig. 1. Overview of the system. a) a remotely operated robotic-manipulator with a tactile sensor (b) is commanded by a human-operator (d) to explore the surfaces of different objects. Material classification results are used to visualise the objects differently in the VR reconstruction (c).

Fig. 2 .

 2 Fig.2. a) The operator manually sets the scan reference pose: the pose of the 3D blue axis mesh is used as a reference of scan location with the second reference at 56.5mm offset. b) the robot approaches to the scan start position at t = 0; c) the slide begins at t = 2; d) the slide direction reverses at t = 6.15; e) the slide is finished at t = 10.3.

Fig. 3 .

 3 Fig.3. Distribution of different initial tactile scanning angular orientations of the sensor. γ is the angle between the objects' surface normal and the tactile sensor's z axis. θ is the angle between the projection of the tactile sensor's z-axis to the scanned object's xy-plane and the object's x-axis (object's x and y-axis are collinear with the robot's with an origin and desired scan position).

Fig. 4 .

 4 Fig. 4.Convolutional neural network (CNN) and multi-modal CNN classifiers. Classifiers are similar except for the additional input layer that contains the robot's end-effector position error and scan orientation quaternion.

Fig. 5 .

 5 Fig. 5. The raw output of the tactile sensor for different materials. Reds: D1, blues: D2, greens: D3, purples: P. Leftmost vertical blue line represents the start of the sliding motion; the middle vertical blue line represents the reverse in slide direction; the rightmost vertical blue line represents the end of the slide.

Fig. 6 .

 6 Fig. 6. Sample spectrograms (not normalized) for metal and soft foam. There are large spikes in 0-15Hz frequency range for both materials (arrows indicate zoom-in regions) and smaller spikes in high frequency for softer/rougher soft foam.

Fig. 7 .

 7 Fig. 7. Confusion matrices of Random Forest, Convolutional Neural Network, Multi-modal Convolutional Neural Network classifiers.

Fig. 8 .

 8 Fig. 8. Classification results visualized in Virtual Reality interface: a) classification results as color-coded Octomaps overlaid on corresponding objects; b) classification results as static color-coded segmented pointcloud clones.

This work is supported by the QMUL-Genoa Univ. PhD Program and by the UK EPSRC grants NCNR EP/R02572X/1 and EP/T027746/1.


John Maxwell Campbell

Proof of a conjecture due to Sun concerning Catalan's constant

Keywords: Catalan's constant, Harmonic number, Dilogarithm, L-function 2010 MSC: 11Y60, 11B65

In 2014, Sun had conjectured that the following formula involving Catalan's constant G holds, letting

This had been discovered experimentally, in the context of the study of special values of L-functions. In this article, we prove this conjecture, using various identities concerning the dilogarithm function, and using an argument relying on a 3 F 2 1 4 -series for Gieseking's constant. It appears that this conjecture due to Sun has not previously been solved, based on the relevant literature.

Introduction

Sun had conjectured that the equality

∞ k=0 2k k (2k + 1)16 k 3H 2k+1 + 4 2k + 1 = 8G (1)
holds true in the 2015 article [START_REF] Sun | New series for some special values of L-functions[END_REF], and this same conjecture is given in Sun's 2014 preprint [START_REF] Sun | List of conjectural series for powers of π and other constants[END_REF] listing his conjectures on powers of π, letting G = ∞ k=0

(-1) k (2k+1)
L(s, χ) = ∞ k=1 χ(k) k 2 , with L (s, χ 4 ) = ∞ k=0 (-1) k (2k + 1) s ,
letting non-principal Dirichlet characters modulo n be denoted as per usual. We are to make use of values associated with

L(2, χ 6 ) = 1 36 ψ (1) 1 6 -ψ (1) 5 6
in our proof introduced in this article that Sun's conjectured formula in [START_REF] Sun | New series for some special values of L-functions[END_REF] holds true.

The formula in (1) is given as Equation (3.18) in [START_REF] Sun | New series for some special values of L-functions[END_REF], which is given as part of part (i) of Conjecture 3.4 in [START_REF] Sun | New series for some special values of L-functions[END_REF]. The formula in (1) is also reproduced as 5

Equation (1.27) in [START_REF] Sun | List of conjectural series for powers of π and other constants[END_REF], with Conjecture 3.4 from [START_REF] Sun | New series for some special values of L-functions[END_REF] being cited. Based on the extant literature citing [START_REF] Sun | New series for some special values of L-functions[END_REF], including [START_REF] Sun | List of conjectural series for powers of π and other constants[END_REF][START_REF] Sun | Supercongruences involving dual sequences[END_REF][START_REF] Sun | Open conjectures on congruences[END_REF][START_REF] Mao | Two congruences involving harmonic numbers with applications[END_REF][START_REF] Chu | Hypergeometric approach to Apéry-like series[END_REF]7,[START_REF] Wang | Alternating multiple zeta values, and explicit formulas of some Euler-Apéry-type series[END_REF]9,[START_REF] Sun | Supercongruences involving Apéry-like numbers and binomial coefficients[END_REF][START_REF] Lupu | Analytic Aspects of the Riemann Zeta and Multiple Zeta Values[END_REF]12,[START_REF] Wang | Explicit formulas of sums involving harmonic numbers and Stirling numbers[END_REF][START_REF] Chu | Alternating series of Apéry-type for the Riemann zeta function, 40[END_REF][START_REF] Chu | Further Apéry-like series for Riemann zeta function[END_REF][START_REF] Lai | Elementary proofs of zagier's formula for multiple zeta values and its odd variant[END_REF][START_REF] Xu | Apéry-type series with summation indices of mixed parities and colored multiple zeta values[END_REF][START_REF] Wang | Proof of some conjectural hypergeometric supercongruences via curious identities[END_REF]19,[START_REF] Wang | Evaluations of some Euler-Apéry-type series[END_REF][START_REF] Xu | Sun's three conjectures on apëry-like sums involving har-55 monic numbers[END_REF], along with [22, §18], it appears that the problem of proving [START_REF] Sun | New series for some special values of L-functions[END_REF] had been open, prior to our proof of (1) introduced in this article.

For the sake of brevity, we assume some basic familiarity with the special

function Li 2 (z) := ∞ k=1 z k k 2 (2)
known as the dilogarithm. The Euler-Mascheroni constant is such that γ = lim n→∞ H n -ln n , and we recall that the digamma function refers to the special function such that the following equalities are satisfied [START_REF] Rainville | Special functions[END_REF], §9):

ψ(z) = d dz ln Γ(z) = Γ ′ (z) Γ(z) = -γ + ∞ n=0 z -1 (n + 1)(n + z)
.

Writing ψ = ψ (0) , the polygamma function is such that

ψ (n) (z) = d n dz n ψ (0) (z).
Theorem 1. Sun's conjectured equality in (1) holds true.

Proof. We are to make use of the moment formula

1 0 (-2k -1)x 2k log(1 -x) dx = H 2k+1 (3)
together with the dominated convergence theorem, in the following manner. By replacing the summand factor H 2k+1 in the series

∞ k=0 2k k (2k + 1)16 k H 2k+1
with the definite integral in (3), and then rewriting the resultant summand so as to obtain

∞ k=0 1 0 (-2k -1) 2k k x 2k log(1 -x) (2k + 1)16 k dx,
we may reverse the order of the operators ∞ k=0 • and 1 0 • dx according to the dominated convergence theorem. So, by the generalized binomial theorem, we have that the equality

∞ k=0 2k k H 2k+1 (2k + 1)16 k = 1 0 - 2 log(1 -x) √ 4 -x 2 dx ( 4
)
holds true. We have determined that an antiderivative for

-2 log(1-x) √ 4-x 2
is as below, and this is easily seen by differentiating the following expression and simplifying:

-

2 iLi 2 - 2e i sin -1 ( x 2 ) -i + √ 3 + iLi 2 2e i sin -1 ( x 2 ) i + √ 3 + 1 2 i sin -1 x 2 2 - sin -1 x 2 log 1 - e i sin -1 ( x 2 ) - √ 3 2 + i 2 -sin -1 x 2 log 1 - e i sin -1 ( x 2 ) √ 3 2 + i 2 + log(1 -x) sin -1 x 2 .
Taking limits as x → 0 and x → 1, this gives us that the definite integral in ( 4) is equal to:

2 iLi 2 - 2 -i + √ 3 + iLi 2 2 i + √ 3 -2iLi 2 - 2e iπ 6 -i + √ 3 - 13iπ 2 36 + 1 3 π log i √ 3 + 1 3 π log 1 - e iπ 6 - √ 3 2 + i 2 .
Rewriting the argument of the first dilogarithmic expression shown above as

Li 2 e -5iπ 6 and then applying a series multisection to (2) according to the residue classes of the indices modulo 6, and similarly for the other dilogarithmic expressions in our evaluation of ( 4), this can be used to show that ( 4) is also equal to

8G 3 - 4iπ 2 27 + 1 9 - √ 3 + i ψ (1) 1 3 + 1 9 √ 3 + i ψ (1) 2 3 ,
and this may be confirmed with Mathematica's FunctionExpand applied to our dilogarithmic form for (4). So, we have that Sun's series in (1) is expressible as

8G - 4iπ 2 9 + 1 3 iψ (1) 1 3 - ψ (1) 1 3 √ 3 + 1 3 iψ (1) 2 3 + ψ (1) 2 3 √ 3 + 4 ∞ k=0 2k k (2k + 1) 2 16 k .
As it turns out, it is actually known1 that the series

∞ k=0 16 -k 2k k (1 + 2k) 2 = 1.0149416064... (5)
is equal to a well-known mathematical constant known as Gieseking's constant [START_REF] Adams | The newest inductee in the number hall of fame[END_REF], with reference to the OEIS entry A143298 and the references therein, and it is known that Gieseking's constant is equal to

9 -ψ (1) 2 3 + ψ (1) 4 3 4 √ 3 .
To show that the series in ( 5) is equal to Gieseking's constant, we may start with the Maclaurin series for the inverse since, so as to write

∞ k=0 2k k z k y 2k 2k + 1 = sin -1 (2y √ z) 2y √ z ,
and by indefinitely integrating with respect to y, we obtain

1 2 √ z sin -1 2y √ z log 1 -e 2i sin -1 (2y √ z) - 1 2 i sin -1 2y √ z 2 + Li 2 e 2i sin -1 (2y √ z)
, and this leads us to an expression equivalent to Gieseking's constant, using previously known dilogarithmic expressions for this constant. So, we have shown that Sun's series in ( 1) is equal to the following: Applying an index shift to rewrite ψ (1) 4 3 in terms of ψ (1) 1 3 , this gives us that Sun's series in (1) is reducible to 8G, as desired.

8G - 4iπ 2 9 + 1 3 iψ (1) 1 3 - ψ (1) 1 3

3 .

 3 Since Sun's series is real-valued, we may omit the complex terms in the above expansion, giving us that Sun'

See https://math.stackexchange.com/questions/3162014, for example.

[START_REF] Mao | Two congruences involving harmonic numbers with applications[END_REF] (1998) 341-349.