Corinne Fournier
email: corinne.fournier@univ-st-etienne.fr

Frédéric Jolivet

Loïc Denis

Nicolas Verrier

Eric Thiebaut

Cédric Allier

T Fournel

Pixel super-resolution in digital holography by regularized reconstruction

Keywords:

teaching and research institutions in France or abroad, or from public or private research centers.

INTRODUCTION

Since in-line holographic imaging was proposed by Dennis Gabor [START_REF] Gabor | A new microscopic principle[END_REF], followed by digital holographic (DH) imaging [START_REF] Royer | Holographic velocimetry of submicron particles[END_REF][START_REF] Schnars | Digital recording and reconstruction of holograms in hologram interferometry and shearography[END_REF], it has become a method of choice for a wide variety of applications including fluid mechanics, mechanical inspection, and biomedical imaging. The method, also known as lensless imaging [START_REF] Repetto | Lensless digital holographic microscope with light-emitting diode illumination[END_REF], relies on a relatively inexpensive and easy to implement configuration in which no collection optics is needed to perform imaging. The role of the collection optics is advantageously replaced by data processing approaches aimed at simulating light back-propagation from the sensor plane to the plane where the objects to be investigated are located [START_REF] Schnars | Digital recording and numerical reconstruction of holograms[END_REF][START_REF] Kreis | Handbook of Holographic Interferometry, Optical and Digital Methods[END_REF][START_REF] Verrier | Off-axis digital hologram reconstruction: some practical considerations[END_REF], thus enabling 3D imaging from the recording of one 2D hologram. Its cost-effectiveness associated with the democratization of high resolution, and high definition imaging sensors made it possible to develop on-chip wide field holographic microscopes [START_REF] Tseng | Lensfree microscopy on a cellphone[END_REF] suited for the detection of bacteria and viruses [START_REF] Allier | Bacteria detection with thin wetting film lensless imaging[END_REF][START_REF] Mudanyali | Wide-field optical detection of nanoparticles using on-chip microscopy and self-assembled nanolenses[END_REF], cytometry [START_REF] Seo | Lensfree holographic imaging for on-chip cytometry and diagnostics[END_REF], or the characterization of protein aggregates [START_REF] Wang | Holographic characterization of protein aggregates[END_REF]. In this case, the resolution of the lensless microscope is driven by the pixel pitch of the chosen sensor, but can be enhanced using pixel super-resolution strategies [START_REF] Bishara | Lensfree on-chip microscopy over a wide field-of-view using pixel super-resolution[END_REF]. Simulation of light back-propagation is still often used for hologram reconstruction because it is simple, although it is prone to border effects and twin image noise, which dramatically reduce the reconstruction signal to noise ratio (SNR) and consequently the accuracy of the reconstruction. Holographic reconstruction is also possible based on another paradigm. Instead of transforming the acquired hologram, the aim is to find the object transmittance that best matches the measurements. This inverse problem (IP) approach extracts more information from the hologram and has been shown to solve two major problems in digital holography: (i) the need to improve the accuracy of the reconstruction, (ii) the need to enlarge the field-of-view beyond the physical limit of the sensor [START_REF] Sotthivirat | Penalized-likelihood image reconstruction for digital holography[END_REF][START_REF] Soulez | Inverse-problem approach for particle digital holography: accurate location based on local optimization[END_REF][START_REF] Soulez | Inverse problem approach in particle digital holography: out-of-field particle detection made possible[END_REF][START_REF] Denis | Inline hologram reconstruction with sparsity constraints[END_REF][START_REF] Bourquard | A practical inverse-problem approach to digital holographic reconstruction[END_REF]. It also leads to almost unsupervised algorithms (only a few tuning parameters are required). These approaches are sometimes referred to as compressive sensing methods [START_REF] Lim | Sampling and processing for compressive holography [invited[END_REF][START_REF] Brady | Compressive holography[END_REF][START_REF] Rivenson | Compressive fresnel holography[END_REF].

To overcome the limitation caused by the pixel pitch of the hologram, in 2010, Ozcan's team introduced pixel superresolution (SR) [START_REF] Bishara | Lensfree on-chip microscopy over a wide field-of-view using pixel super-resolution[END_REF], which was already used in other modalities [START_REF] Park | Super-resolution image reconstruction: a technical overview[END_REF][START_REF] Vandewalle | A Frequency Domain Approach to Registration of Aliased Images with Application to Super-Resolution[END_REF]. These authors suggested it was possible to recover a higher resolution hologram using a stack of rigidly shifted lower resolution holograms and to reconstruct them by back-propagation. Recently, it was shown that coupling IP with pixel SR can improve the accuracy of hologram reconstruction even more [START_REF] Verrier | Digital holography super-resolution for accurate three-dimensional reconstruction of particle holograms[END_REF][START_REF] Verrier | 3d tracking the brownian motion of colloidal particles using digital holographic microscopy and joint reconstruction[END_REF]. In this IP approach, a high resolution hologram is no longer computed, instead the reconstruction is performed directly from the hologram stack. However, despite their accuracy, these new DH-SR algorithms are based on a parametric reconstruction limited to the reconstruction of parametric objects (i.e. relying on a few parameters). For arbitrary shaped objects, the object 3D space can be sampled and a transmittance value estimated at each voxel. The image-hologram formation model can then be inverted using an a priori on the objects by means of a prior on the 3D object. This maximum posteriori (MAP) reconstruction has been shown to reduce twin images and to extend the usual field-of-view under sparse prior [START_REF] Denis | Inline hologram reconstruction with sparsity constraints[END_REF] or total variation prior [START_REF] Sotthivirat | Penalized-likelihood image reconstruction for digital holography[END_REF][START_REF] Bourquard | A practical inverse-problem approach to digital holographic reconstruction[END_REF][START_REF] Brady | Compressive holography[END_REF][START_REF] Marim | Off-axis compressed holographic microscopy in low-light conditions[END_REF][START_REF] Marim | Compressed sensing with off-axis frequency-shifting holography[END_REF]. We propose to combine pixel SR with IP approaches to perform optimal high resolution hologram reconstruction of arbitrary sparse objects from a stack of shifted holograms. This approach has already been used for conventional imaging [START_REF] Vandewalle | A Frequency Domain Approach to Registration of Aliased Images with Application to Super-Resolution[END_REF], we apply it to DH. The specificity of DH makes the results unique: high accuracy registration, out-of-field reconstruction. In the first part of the article, the pixel super-resolution issue is presented from an IP point of view. Then, an algorithm that alternates registration steps and transmittance reconstruction is presented. Finally, the approach is validated with both synthetic and experimental data resulting in high SNR transmittance reconstruction over an enlarged field-of-view. Here the proof of concept is limited to the reconstruction of a 2D plane, but can be extended to 3D sparse objects.

SR RECONSTRUCTION PROBLEM AS AN INVERSE PROBLEM

Pixel SR algorithms in DH attempt to extract a high-resolution opacity distribution (real or complex valued) from a sequence of low-resolution holograms. This issue can be solved using an IP because the imaging model (direct model) is well-known. In this section, we first describe the direct model and then mathematically formalize the reconstruction problem.

A. Direct model

The direct model gives a mathematical expression of the hologram intensity as a function of the object opacity. Let us consider the simple case of sparse semi-opaque objects located in a plane and illuminated by a collimated coherent beam of wavelength λ. This object is described by its opacity distribution ϑ, which is assumed to be real. A sensor positioned at a distance z from the object plane (Fig. 1) records the intensity of the wave diffracted by the object.

Assuming the holographed object l is narrow in width and satisfies πl 2 /(4λz) 1, the hologram intensity I(x, y) can be approximated as an incoherent summation of diffraction patterns created by each point of the object plane characterized by its opacity ϑ.

I(x, y) ∝ [1 -2.ϑ * h z ](x, y) (1)
where h z stands for the real part of the free space propagation impulse response, * for the 2D spatial convolution operator and ∝ is the proportionality sign.

In the Fresnel diffraction regime [START_REF] Joseph | Introduction to fourier optics[END_REF] (i.e. z 3 π l 4 /(64λ)), h z is given by:

h z (x, y) = 1 λz sin π(x 2 + y 2 ) λz . (2)
Accounting for pixel integration [START_REF] Kreis | Frequency analysis of digital holography with reconstruction by convolution[END_REF], a transversal shift of the sensor Δ i = (Δx i , Δy i ) and the sensor sampling, the intensity model of the i th hologram can also be expressed in a linear form: where (x p , y p ) is the center of the p th pixel, Π is the 2D rectangular function representing the photo-sensitive area of a pixel, δ Δ i is the Dirac distribution centered on (Δx i , Δy i ).

I i (x p , y p ) ∝ 1 -2.ϑ * h z * Π * δ Δ i (x p , y p ), (3)
This expression can be rewritten as

I i (x p , y p ) ∝ 1 -ϑ * h Π z,Δ i (x p , y p ), (4)
with h Π z,Δ i = 2.h z * Π * δ Δ i . (5)
In a pre-process step the holograms are divided by a background image (image with no object in the field), so that the proportionality factor is removed.

In the following, we use matrix notation as a simpler mathematical expression of the direct model. Bold notations are used to indicate a vector or a matrix. The image is represented in vector form by considering the 2D image as a 1D vector by stacking all columns of the image. The convolution of the kernel h Π z,Δ i with the sampled opacity ϑ is thus expressed as the multiplication of a large discrete convolution matrix H Δ i with the vector corresponding to the object opacity ϑ. The matrix H Δ i is made of pixel-shifted replications of h Π z,Δ i put into a vector form. The direct model m i representing the intensity of the i th hologram is in matrix notation:

m i = 1 -H Δ i ϑ (6)
where 1 represents the offset and is a vector of n "ones", H Δ i models diffraction and takes into account the relative displacement of the object with respect to the camera Δ i and ϑ is the unknown opacity distribution of the object, common to all holograms. Considering a centered model mi and a centered diffraction matrix H, the expression of the model is simplified [START_REF] Denis | Inline hologram reconstruction with sparsity constraints[END_REF]

to: mi = HΔ i ϑ ( 7
)
with mi = m i -1 n 11 t m i and H = -H + 1 n 11 t H. From the point of view of implementation, to evaluate this model, the convolutions are computed in the Fourier space using Fast Fourier Transform. The sampling of the model is thus given by the opacity distribution sampling. A super-resolved opacity distribution leads to a super-resolved model:

mi,SR = HΔ i ϑ SR (8)

B. SR Reconstruction

In an IP framework, the cost function enforces a certain fidelity of the final solution to the measured data. For the SR issue, the goal is to jointly reconstruct the super-resolved opacity from the whole hologram sequence. The fidelity term is thus composed of a sum of individual cost functions with a common super-resolved object opacity. Modeling the error between the modeled intensity and the measured hologram as white and Gaussian, the negative log-likelihood is proportional to:

1 N N ∑ i=1 HΔ i ϑ SR - di,SR 2 w (9)
where:

• N is the number of holograms,

• di is the i th centered hologram.

• the exponent ( SR ) indicates that data sampling is the same as opacity sampling ϑ SR ; in practice zeros have been inserted in di,SR at the position of missing pixels (e.g. for a super-resolution factor of 2, one column out of two are considered as missing (Fig. 3)),

• w is the weight vector that accounts for missing data (for pixel that are missing the weight is 0, otherwise it is 1).

More generally, it can also account for non-uniform noise variance [START_REF] Soulez | Inverse-problem approach for particle digital holography: accurate location based on local optimization[END_REF], [START_REF] Fournier | Digital hologram processing in on-axis holography[END_REF],

• . 2 w is the weighted L 2 norm defined as:

u 2 w = u, u w = ( n ∑ p=1 w p u 2 p )/( n ∑ p=1 w p ).
This inverse problem is ill-posed. In order to stabilize the inversion process, a regularization penalty function ρ, that enforces an a priori on the object to reconstruct, is added. This regularization is applied to the reconstructed superresolved object opacity ϑ SR .

The expression of the cost function becomes:

C ϑ SR , {Δ i } = 1 N N ∑ i=1 HΔ i ϑ SR - di,SR 2 w + τρ(ϑ SR ) (10)
where the coefficient τ indicates the weight of the regularization term. Reconstructing a super-resolved opacity using an IP framework is equivalent to solving the minimization problem:

θSR = arg min 0 ϑ SR 1 min Δ i C ϑ SR , {Δ i } ( 11
)
where bound constraints are imposed on the opacity ϑ to account for the physics of absorption: the opacity of semiopaque objects being between 0 (transparent) and 1 (opaque).

This IP formulation provides a rigorous mathematical framework for the SR problem. Indeed, the shifted holograms of the sequence are modeled using the image formation model and consequently, no arbitrary interpolation function is used. The interpolation is implicitly made by the image formation model. Furthermore, the penalty function (based on an a priori on the objects) acts as a denoising step in the reconstruction. In the following section, we detail how the minimization in equation ( 11) is achieved.

ALTERNATING OPTIMIZATION ALGORITHM

The joint minimization of Eq. 11, consists in optimizing two different quantities ϑ SR and {Δ i } i=1..N . It can be performed iteratively by alternating minimization [START_REF] Hardie | Joint map registration and high-resolution image estimation using a sequence of undersampled images[END_REF]. The cost function is minimized with respect to the high resolution object opacity ϑ SR and the shift parameters {Δ i } i=1..N alternatively until convergence of the algorithm. Figure 2 represents the k th step of the alternating optimization. The two minimization steps are detailed below.

A. Registration step

It is widely accepted that the accuracy of the shift estimation in the image sequence is of paramount importance for the efficiency of SR algorithms [START_REF] Farsiu | Advances and challenges in super-resolution[END_REF]. The shifts are usually estimated relative to a reference image (e.g., the first one of the sequence). The estimation is then performed by minimizing the square difference between the first image and a shifted version of the i th image. As the accuracy of the estimate is subpixel, interpolation functions have to be used. Whatever the choice of the interpolation function (linear, cubic spline function, etc.), it is arbitrary. The use of the IP approach fixes this issue, as interpolation is implicitly performed by the image formation model. The transverse coordinate system (O, x, y) of the first hologram is considered as reference. The opacity distribution is reconstructed relatively to this coordinate system. The shift Δ i of the i th hologram relative to the first one is also the shift between the i th hologram and the image formation model Hϑ SR . Thus, the problem is to minimize the residual between the shifted model and the i th hologram:

Δ i = arg min Δ i HΔ i ϑ SR - di,SR 2 w ( 12
)
An illustration is given in Fig. 3. The support of the model is larger than the data because of the diffraction and outof-field reconstruction of the opacity. Provided the opacity is well reconstructed, the shift estimate can be expected to be more accurate than by performing the registration on noisy holograms. Indeed, in our approach, the problem is equivalent to registering a large noiseless model with a hologram (which consists of a noisy shifted model), whereas in state-of-the-art approach, holograms are registered 2 by 2. The noise variance in our registration step is theoretically divided by two. It theoretically decreases the standard error on the shifts by a factor √ 2 the standard error on the shifts. Furthermore there are no cropping effects, because the model is much larger than the holograms.

Direct minimization of this cost function is computationally too expensive because it involves the computation of a large model for each model evaluation (e.g. six times the width of the data if the SR factor is 2 and the field-of-view extension ratio is 3) and 6 FFT to compute the weighted centered correlation. In practice, the shift estimation can be performed in a sampled space of the shifts. Equation ( 12) can be rewritten as a maximization of the weighted normalized centered cross correlation [START_REF] Soulez | Inverse problem approach in particle digital holography: out-of-field particle detection made possible[END_REF]:

Δ i = arg max Δ i HΔ i ϑ SR , di,SR w HΔ i ϑ SR 2 w di,SR 2 w ( 13
)
To avoid calculating multiple subpixel shifted models, the correlation is calculated on a map corresponding to integer shifts. The maximum of the correlation map is then positioned with subpixel accuracy by fitting parabolas (vertical and horizontal) on the correlation peak. The first guess concerning the shifts is obtained by correlation-based registration of the holograms.

B. Opacity reconstruction step

In the second optimization, the shifts are fixed and the opacity is estimated using a penalized maximum likelihood. Regularization may differ depending on the a priori on the object (e.g. total variation for objects with piecewise-constant opacity and sharp edges). Here, we assume a sparse semi-opaque object. A natural constraint is thus the L 1 norm of the opacity: ρ(ϑ SR ) = ϑ SR 1 . ϑ SR is also assumed to be positive and less than 1. The reconstructed opacity then corresponds to the solution of the following minimization problem:

ϑ SR = arg min 0≤ϑ SR ≤1 1 N N ∑ i=1 HΔ i ϑ SR - di,SR 2 w + τ ϑ SR 1 ( 14
)
Under positivity constraints, the L 1 norm corresponds to the sum of the opacities: ϑ SR 1 = ∑ j ϑ SR j . The minimization problem ( 14) is a smooth optimization problem under bound constraints that can be solved efficiently using a limited-memory quasi-Newton method with bound constraints, such as VMLMB [START_REF] Thiébaut | Optimization issues in blind deconvolution algorithms[END_REF]. Proximal algorithms are also widely used for L1 minimization. Use of the FISTA algorithm [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF] (an accelerated proximal method) to the problem (14) leads to Algorithm 1. We found that it converges slightly more slowly than VMLMB, but is useful because of its straightforward implementation, and could easily be accelerated by implementation on the GPU [START_REF] Endo | GPU-accelerated compressive holography[END_REF].

η 1 ← ϑ SR 0 t 1 ← 1 for k = 1, 2, . . . ϑ SR k ← S + τ,1 η k + 2 N N ∑ i=1 H * Δ i • diag(w) • di,SR -HΔ i η k (15) t k+1 ← 1 + 1 + 4 t 2 k 2 (16
)

η k+1 ← ϑ SR k + t k -1 t k+1 ϑ SR k -ϑ SR k-1 ( 17
)
where the positive soft-thresholding operator is applied coordinate-wise:

S + τ,α [u] = ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ u - τ 2 if α + τ 2 > u > τ 2 0 i fu ≤ τ 2 α if u ≥ α + τ 2 (18)
and H * Δ i represents the adjoint operator HΔ i . The soft-thresholding operator is simply generalized by considering a spatially varying threshold τ . This is useful to account for the difference between the norms of each line of H. To avoid over-penalizing elements with a small norm, we define τ i = τ ∑ j w j H2 i,j / ∑ j w j . This normalization can be negligible in the field-of-view but is mandatory for out of the field reconstruction.

EXPERIMENTS

Compared to reconstructions based on a single low-resolution hologram, the MAP-SR is expected to recover higher frequency information in the hologram space, thereby leading to i) better resolved reconstructions, and ii) an enlargement of the reconstructed field-of-view. Demonstrating the two improvements requires different experimental configurations.

To reconstruct an enlarged field-of-view (ii), a large recording distance is needed to capture the fringes coming from objects that are outside the field-of-view (due to pixel integration, the out-of-field signal is filtered out when the recording distance is too short). At the same time, large recording distances imply that high-frequency fringes from in-the-field objects fall outside the sensor area. The performance of the present method is demonstrated on both simulated and real holograms. The improvement in resolution (i) is demonstrated in the field-of-view of simulated holograms. The accuracy of shift estimation and the quality of the reconstruction are quantitatively compared with state-of-the-art approaches. The second expected improvement (ii) is shown on a stack of experimental holograms with a deliberately very low SNR in order to demonstrate that the proposed reconstruction is also efficient with noisy data. The quality of the reconstruction is compared with state-of-the-art reconstructions.

A. Synthetic data

Knowledge of the opacity distribution of the synthetic object and of the virtual relative shifts of the object with respect to the camera is used to quantify the errors in the reconstruction and in shift estimates.

A.1. Simulated holograms

In order to compute synthetic holograms, a plane transmittance object is simulated (Fig. 5). The wave amplitude on the sensor is computed by propagating the amplitude wave of the transmittance plane to the sensor and considering a plane reference wave. The complex Fresnel free space propagation model is used. Finally, the intensity is calculated by the square modulus of the complex amplitude. This simulation, in contrast to the direct model used in the reconstruction (see section A), is closer to the optical model but is non-linear (which is why it was not used to derive the algorithm). It consists of computing a high-resolution complex amplitude of the propagated wave on the sensor and then computing its intensity (square modulus of the complex amplitude). The pixel integration effect is simulated by convolving the intensity by a 2D rectangle function of the final pixel size. To simulate a shifted hologram, the Fresnel kernel, which is calculated from its analytical expression, is shifted with random shifts in x and y directions. To model a low resolution sensor, an under-sampling step by the same factor f SR is applied (see Fig. 4). Lastly, white Gaussian noise is added to simulate the electronic noise and a 12 bit quantization is applied on the holograms.

To simulate the hologram stacks, a 512 × 512 sampled transmittance distribution with a dot pitch of 3.7μm is used (see Fig. 5). The wavelength is set to 662nm and the object-sensor distance to z = 16.5mm. Holograms of 256 × 256 pixel with a pitch of 7.4μm and a fill-factor (i.e., active area over the total area of the pixel) of 0.5 are simulated. Stacks of N = 50 holograms are generated with shifts uniformly distributed in the range [-2, +2] pixels. The SNR (ratio of the magnitude of the signal to the standard deviation of the noise) is set to 5 in order to be close to an experimental value. The hyperparameter τ is set by minimizing the mean square deviation between the reconstruction of simulated holograms and the true transmittance.

A.2. Accuracy of the hologram stack registration

The first advantage of the MAP-SR method is a more accurate estimation of the shifts. In state-of-the-art approaches, a hologram is registered to another hologram [START_REF] Bishara | Lensfree on-chip microscopy over a wide field-of-view using pixel super-resolution[END_REF][START_REF] Vandewalle | A Frequency Domain Approach to Registration of Aliased Images with Application to Super-Resolution[END_REF], whereas in the proposed approach, registration is between a noiseless model and the noisy data, as discussed in section A. Assuming a perfect reconstruction and ignoring boundary effects, a theoretical gain of √ 2 in precision for the estimated shifts is expected. In contrast to conventional methods which only consider the overlapping area between two holograms to perform registration, our approach exploits all the pixels of the holograms thus avoiding boundary effects. Nevertheless, incorrectly estimated opacities in our model may cause registration errors.

To check the improvement of the shift estimations in the reconstruction, simulations were run on 100 stacks of 50 holograms with the same random shifts and SNR as described in section A.1. The first shift guesses are obtained by correlation based registration of the i th hologram with the first hologram. The errors in the shifts are shown in figure 6. They decrease during the course of the opacity reconstruction process. At the fifth iteration, the accuracy gain is 7.0 in y direction and 10.6 in x direction, which is a significant improvement in shift estimation.

A.3. Quality of the reconstruction

Four reconstructions were performed (Fig. 78): (a) back-propagation of a super-resolved hologram obtained with a state-of-the-art algorithm [START_REF] Kreis | Handbook of Holographic Interferometry, Optical and Digital Methods[END_REF][START_REF] Vandewalle | A Frequency Domain Approach to Registration of Aliased Images with Application to Super-Resolution[END_REF], (b) MAP reconstruction using one image [START_REF] Denis | Inline hologram reconstruction with sparsity constraints[END_REF] and (c) MAP reconstruction using the same super-resolved hologram as in (a), (d) the proposed IP-MAP SR reconstruction.

Figure 7a shows that the back-propagation reconstruction suffers from twin image noise. Figure 7b shows that the reconstruction using a single hologram is worse than the reconstructions using a stack of holograms displayed in figures 7c and d. To better visualize the difference between figure 7c and figure 7d, zoomed versions are provided in figure 8. The SR MAP reconstruction (d) is the best reconstruction. For a quantitative comparison, the peak signal to noise ratio (PSNR) was computed for each reconstruction : (a) 16.2 dB, (b) 24.1 dB, (c) 28.0 dB, (d) 39.2 dB. It can be seen that the PSNR was improved by a factor of 2.4 from reconstruction (a) to (d).

B. Experimental data

Using IP-MAP SR reconstruction, the simulations showed an improvement in the reconstruction in the field-of-view. For experimental holograms, we demonstrate the improvement in the reconstruction outside the usual field-of-view using experimental holograms. To check the robustness of the proposed approach versus noise, the SNR of experimental holograms was deliberately reduced by spraying water on the light path. The setup is shown in figure 9. Sprayed water droplets create a realistic disturbance that decreases the SNR of the object signal. Between each of the N recordings, the object was randomly shifted transversely (versus x and y). The test object is a metric crossed reticle (Edmund Optics) which has a sparse binary transmittance. The camera used for the experiments was a Prosilica GE 4900 camera. Its pixel pitch is 7.4μm and its fill-factor is around 1. The recorded holograms are 1024 × 1024 pixels. The wavelength of the laser is 662nm. The camera was positioned at z = 283mm from the object plane to record fringes of objects located outside the field-of-view.

Reconstruction was performed using a SR factor of 2 (the dot pitch of the opacity plane is 3.7μm). The hyperparameter τ was set empirically. The width of the field-of-view increased by a factor of 3 leading to a 2.3cm wide field-of-view and a surface of 5.2cm 2 .

In figure 10, the same four reconstructions, as in the synthetic data case, are displayed for the purpose of comparison. It can be seen that the proposed reconstruction in figure 10d makes it possible to reconstruct an object 3 times wider than the sensor, whereas using the back-propagation reconstruction method, it is difficult to reconstruct beyond the sensor field. This out-of-field reconstruction not only enables enlargement of the field-of-view but also renders shift estimation more accurate. Figure 11, shows that the in-field reconstruction is less noisy using MAP-SR reconstruction (d). The supplementary material (Visualization 1..4) makes it possible to zoom in to check the reconstruction of the border of the reticle. It should be noted that because the holograms are better registered with the MAP-SR approach (d), the hyperparameter is smaller in (d) than in (c) and hence easier to tune. In the example, the reconstruction (d) was obtained with a hyperparameter two times smaller than in (c).

CONCLUSION

An inverse problem approach was successfully applied to the super-resolved reconstruction of holograms. Starting from a collection of unregistered holograms, the proposed method alternates between a step of sub-pixel registration of the holograms and a step of joint inversion of all holograms. This method differs from recent super-resolution methods applied in digital holography and lensless imaging in two major aspects: (i) sub-pixel hologram registration is performed by matching holograms with a (noiseless) model obtained from the previous reconstruction, rather than by matching holograms to a reference (noisy) hologram; (ii) the reconstruction method inverts the hologram formation model by including prior knowledge (bound constraints, sparsity), rather than by light back-propagation. Those two modifications result in better reconstructions with reduced artefacts (twin-images, border effects), an extended field-ofview, a better signal-to-noise ratio, and improved spatial resolution. These improvements come at the cost of increased computational complexity compared to existing pixel super-resolution approaches. We believe the improvement in the reconstructions outweighs the drawback of this higher computational cost. In applications where reconstruction time must be kept short, special effort should be invested in accelerating the Fresnel propagator, using GPU implementations, and refining the minimization strategy, for example using continuation schemes.

Fig. 1 .

 1 Fig. 1. Illustration of the in-line hologram formation model.

Fig. 2 .

 2 Fig. 2. Illustration of the k th step of the alternating optimization algorithm.

Fig. 3 .

 3 Fig. 3. Illustration of the shift estimation step for a SR factor of 2. The unshifted and super-resolved hologram model ( HΔ i ϑ SR ) is shown in the background. The shifted and low-resolution captured holograms d i are superimposed The gray pixels corresponds to missing data in the super-resolved space. The white pixels corresponds to the measured gray levels. The pixel size and the shift lengths are enlarged for the purpose of illustration.

Algorithm 1 :

 1 Reconstruction of a sparse distribution of opacities input: holograms di,SR and shifts Δ i output: reconstructed opacity distribution ϑ SR initialization: initial opacity distribution ϑ SR 0 (can be set to zero)

Fig. 4 .

 4 Fig. 4. Illustration of the holograms simulation pipeline to generate the stack.

Fig. 5 .

 5 Fig. 5. Transmittance used to simulate holograms.

Fig. 6 .

 6 Fig. 6. Illustration of the refinement of x and y shift estimation throughout the opacity reconstruction process. The size of error bars represents 2 times the standard deviation on the shifts.

Fig. 7 .Fig. 8 .Fig. 9 .

 789 Fig. 7. Qualitative comparison of the reconstructions: (a) Back-propagation reconstruction using a super-resolved hologram [6, 13], (b) MAP reconstruction using 1 hologram [17], (c) MAP reconstruction using the same super-resolved hologram as in (a), (d) IP-MAP SR reconstruction using a stack of 50 holograms. The super-resolved hologram used for (a) and (c) is computed from the stack of 50 holograms with a state-of-the-art pixel super-resolution algorithm.

Fig. 10 .

 10 Fig. 10. Qualitative comparison of experimental hologram reconstructions: (a) Back-propagation reconstruction using a super-resolved hologram, (b) MAP reconstruction using 1 hologram [17], (c) MAP reconstruction using the same super-resolved hologram as in (a), (d) IP-MAP SR reconstruction using a stack of 15 holograms (see Visualization 1, 2, 3 and 4). The super-resolved hologram used for (a) and (c) was computed from the stack of 15 holograms with a stateof-the-art pixel super-resolution algorithm. The blue square in the center of the image represents the sensor borders. The green rectangle shows the zoomed areas displayed in figure 11.

Fig. 11 .

 11 Fig. 11. Zoom on the reconstructions displayed on Figure.10.a-d

Fundings

This project is funded by the Agence Nationale de la Recherche (ANR) (ANR-1 1-IDEX-0007, ANR-11-LABX-0063); DETECTION-CNRS DEFI IMAGIn 2015. Funding for this project was provided by a grant from la Région Auvergne-Rhône-Alpes.


MAHLER MEASURE ON GALOIS EXTENSIONS Francesco Amoroso

est destinée au dépôt

Introduction

Let α be a non-zero algebraic number of degree d. We let as usual M (α) the Mahler measure of α. Thus

M (α) = |lc(f )| d j=1 max(|α j ), 1)
where lc(f ) is the leading coefficient of a minimal equation of α over Z and where α 1 , . . . , α d are the conjugates of α.

Assume that α is not a root of unity (otherwise M (α) = 1). By an old result of Kronecker M (α) > 1 and Lehmer [START_REF] Lehmer | Factorization of certain cyclotomic functions[END_REF] asked if we could replace 1 by a real number > 1 which does not depend on α. This problem is still open, the best known result is a theorem of Dobrowolski (see [START_REF] Dobrowolski | On a question of Lehmer and the number of irreducible factors of a polynomial[END_REF]) which proves the lower bound

M (α) ≥ 1 + c(ε)d -ε
for all ε > 0, with c(ε) > 0 depending only on ε.

Recently, a construction of Smyth gives a renewed interest in lower bounds for the Mahler measure of a generator of a Galois extension (a problem first considered in [START_REF] Amoroso | Le problème de Lehmer en dimension supérieure[END_REF]). In this special case we can considerably improve the above lower bound. Let α be a non-zero algebraic number of degree d, not a root unity, such that Q(α)/Q is Galois. Then (see [START_REF] Amoroso | Lower bounds for the height in Galois extensions[END_REF])

M (α) 1/d ≥ 1 + c(ε)d -ε
for all ε > 0, with c(ε) > 0 depending only on ε.

The result of [START_REF] Amoroso | Lower bounds for the height in Galois extensions[END_REF] was partially motivated by a problem posed by Smyth during a recent BIRS workshop (see [3, problem 21, p. 17]). Let n ≥ 2 and β = β 1 , . . . , β n be the roots of z n -z -1, known to be irreducible for all n, and to have Galois group the full symmetric group S n . Put

α = β 1 1 β 2 2 . . . β n-1 n-1 .
Date: 20 august 2017.

Then (by an easy consequence of [START_REF] Smyth | Additive and Multiplicative Relations Connecting Conjugate Algebraic Numbers[END_REF] Lemma 1) the Galois closure of The quoted result of [START_REF] Amoroso | Lower bounds for the height in Galois extensions[END_REF] is related to the first question. One of the aim of this paper is to give an answer to the second question. More generally, in section 3 we show:

Q(β) is Q(α) of degree d = n! over Q. Smyth computed with Maple the first values of M (α) 1/d n d = n! M (α)
Theorem 1.

= a j -1 n i a i . Then M (β) cn|y| 1 ≤ M (α) 1/n! ≤ M (β) |y| 1 with |y| 1 = 1 n n j=1 |y j | and where c n ∼ 2 πn . 3) If α is a generator of Q(β 1 , . . . , β n )/Q, M (α) 1/n! ≥ M (β) (1+o(1)) √ n 8π .
The Mahler measure of a root of x n -x -1 is ≥ θ = 1.32..., the smallest Pisot's number (the root > 1 of x 3 -x -1). Thus Theorem 1.1 implies that Smyth's sequence tends to ∞.

A more classical way to construct a generator for the Galois closure of Q(β) is given by the proof of the Primitive Element Theorem, thus taking a general linear combination of the conjugates of β. In section 4 we give a proof, based on a simple discriminant argument, of the following partial analogous of Theorem 1.1.

Theorem 1.2. Let β be an algebraic integer of degree n, with algebraic conjugates

β 1 , . . . , β n . Assume that Q(β 1 , . . . , β n )/Q has Galois group S n . Let a 1 , . . . , a n ∈ Z and put α = a 1 β 1 + a 2 β 2 + • • • + a n β n . Then: 1) α is a generator of Q(β 1 , . . . , β n )/Q if and only if a 1 , . . . , a n are pairwise distinct. 2) Let V (a) be the Vandermonde 1≤i<j≤n (a j -a i ). Then, if n ≥ 5, 2 -1 24 |V (a)| • |disc(β)| 1/2 1 12n(n-1) ≤ M (α) 1/n! ≤ |a| 1 • M (β).
3

) If α is a generator of Q(β 1 , . . . , β n )/Q, M (α) 1/n! ≥ (n/4) 1/48 .
Theorems 1.1 and 1.2 may suggest some speculations on the behavior of M (α) for α a generator of a Galois extension.

Conjecture 1.3. Let α ∈ Q be a generator of of a Galois extension of degree d = n! with Galois group the full symmetric group S n . Then,

M (α) 1/d ≥ C(d) with C(d) → +∞ for d → +∞. Equivalently, in terms of the Weil height h(α) = 1 d log M (α), h(α) ≥ c(d)
with c(d) growing to infinity with d. This would be the first result of the kind "h → +∞".

Remark that the assumption on the Galois group is necessary, as a couple of exemples of [START_REF] Amoroso | Lower bounds for the height in Galois extensions[END_REF] show: take ζ e a e-root of unity and put α = 1 + ζ e or α = 2 1/e + ζ e , both of uniformly bounded height.

Aknowledgements. We are indebted to Eric Ricard for a convexity argument which is the key of the proof of Lemma 2.2. We warmly thank Tanguy Rivoal for equality (2.3).

Auxiliary results

Let H

n := {x ∈ R n | x 1 + • • • + x n = 0} and, for x ∈ H n , |x| 1 = 1 n n j=1 |x j | .
We remark that S n acts on H n (by σ(x) = (x σ(1) , . . . , x σ(n) )).

For x, y ∈ H n we set

s n (x, y) = 1 n! σ∈Sn 1 n n j=1
x σ(j) y j .

Thus s n is symmetric,

s n (x, y) = s n (y, x).
Moreover, for each y ∈ H n , the function x → s n (x, y) is a seminorm 1 , stable by the action on S n .

For h = 1, . . . , n -1 we define a vector z (n,h) ∈ H n by

z (n,h) j = n 2h if j = 1, . . . , h -n 2(n-h) if j = h + 1, . . . , n (thus |z (n,h) | 1 = 1). Let also c n = min 0<h,k<n s n (z (n,h) , z (n,k) ).
The aim of this section is to prove the following proposition.

Proposition 2.1. For x, y ∈ H n \{0} we have

(2.1) c n ≤ s n (x, y) |x| 1 • |y| 1 ≤ 1 . Moreover (2.2) c n ∼ 2 πn .
The above lower bound is clearly sharp, and the upper bound is almost sharp, since s n (z (n,1) , z (n,k) ) = n 2(n -1) (see Lemma 2.3 below).

Proposition 2.1 easily implies the lower bound for M (α) in assertion 2) of Theorem 1.1, as we shall see in the next session. In the rest of this section we present the proof of the proposition, which follows from the three lemmas below. Lemma 2.2. Let x, y ∈ H n \{0}. Let h = h(x) be the cardinality of the set of j ∈ {1, . . . , n} such that x j ≥ 0 and similarly for k = k(y). Then

s n (x, y) |x| 1 • |y| 1 ≥ s n (z (n,h) , y) |y| 1 ≥ s n (z (n,h) , z (n,k) ).
Proof. Let us prove the first inequality. Since x → s n (x, y) is homogeneous, we can assume |x| 1 = 1. Let A be the set of j ∈ {1, . . . , n} such that x j ≥ 0. Since x → s n (x, y) is invariant by the action of S n , we can assume A = {1, . . . , h}. We now use a convexity argument suggested by E. Ricard. Let G be the subgroup of σ ∈ S n such that σ(A) = A. Then

|G| -1 σ∈G σ(x) = z (n,h)
and, by the convexity of x → s n (x, y),

s n (z (n,h) , y) ≤ s n (x, y).
1 and indeed a norm if y = 0, for instance as a consequence of Proposition 2.1 below.

The second inequality follows from the first one and from the symmetry of s n :

s n (z (n,h) , y) |y| 1 = s n (y, z (n,h) ) |y| 1 ≥ s n (z (n,k) , z (n,h) ) = s n (z (n,h) , z (n,k) ).
Lemma 2.3. Let y ∈ H n and h, k ∈ {1, . . . , n}. Then

s n (z (n,h) , y) = n 2h(n -h) n h -1 S⊆{1,...,n} Card(S)=h j∈S y j and s n (z (n,h) , z (n,k) ) = n 2 (h -[hk/n])(k -[hk/n]) 2hk(n -h)(n -k) n h -1 k [hk/n] n -k h -[hk/n] .
Proof. We have

s n (z (n,h) , y) = 1 n! σ∈Sn 1 n n j=1 z (n,h) j y σ(j) = 1 n! σ∈Sn 1 2h h j=1 y σ(j) - 1 2(n -h) n j=h+1 y σ(j) = 1 n! σ∈Sn 1 2h + 1 2(n -h) h j=1 y σ(j) = n 2h(n -h) 1 n! σ∈Sn h j=1 y σ(j) = n 2h(n -h) n h -1 S⊆{1,...,n} Card(S)=h j∈S y j .
We now prove the second equality. From the first,

s n (z (n,h) , z (n,k) ) = n 2h(n -h) n h -1 S⊆{1,...,n} Card(S)=h j∈S z (n,k) j . For S ⊆ {1, . . . , n} of cardinality h we have j∈S z (n,k) j = n 2k |S ∩ {1, . . . , k}| - n 2(n -k) |S ∩ {k + 1, . . . , n}| = n 2k + n 2(n -k) |S ∩ {1, . . . , k}| - n 2(n -k) |S| = n 2 2k(n -k) |S ∩ {1, . . . , k}| - hk n .
This gives

s n (z (n,h) , z (n,k) ) = n 3 4h(n -h)k(n -k) n h -1 Σ n,h,k with Σ n,h,k = S⊆{1,...,n} Card(S)=h |S ∩ {1, . . . , k}| - hk n = min(h,k) j=max(h+k-n,0) S 1 ⊆{1,...,k} Card(S 1 )=j 1 S 2 ⊆{k+1,...,n} Card(S 2 )=h-j 1 j - hk n = min(h,k) j=max(h+k-n,0) k j n -k h -j j - hk n .
We now quote an equality suggested by T. Rivoal. Let q ∈ Z with max(h + k -n, 0) ≤ q ≤ min(h, k).

Then (2.3) q j=max(h+k-n,0) k j n -k h -j hk n -j = 1 n (h-q)(k-q) k q n -k h -q .
This formula can be easily verified by induction on the parameter q. It shows that

min(h,k) j=max(h+k-n,0) k j n -k h -j hk n -j = 0 and that Σ n,h,k = 2 [hk/n] j=max(h+k-n,0) k j n -k h -j hk n -j = 2 n (h -[hk/n])(k -[hk/n]) k [hk/n] n -k h -[hk/n] . Thus s n (z (n,h) , z (n,k) ) = n 3 4h(n -h)k(n -k) n h -1 × 2 n (h -[hk/n])(k -[hk/n]) k [hk/n] n -k h -[hk/n] = n 2 (h -[hk/n])(k -[hk/n]) 2hk(n -h)(n -k) n h -1 k [hk/n] n -k h -[hk/n] .
In the next lemma we give an asymptotic estimate for

s n (z (n,h) , z (n,k) ) = n 2 (h -[hk/n])(k -[hk/n]) 2hk(n -h)(n -k) n h -1 k [hk/n] n -k h -[hk/n] . Lemma 2.4. Let (n m ) n∈N , (h m ) m∈N , (k m ) m∈N with 0 < h m , k m < n m . Assume lim m→+∞ n m = +∞, lim m→+∞ h m n m = u, lim m→+∞ k m n m = v for some u, v ∈ [0, 1]. Then lim m→+∞ 2 2πuv(1 -u)(1 -v)n m • s nm (z (nm,hm) , z (nm,km) ) = 1.
Proof. By a continuity argument, we can assume u, v ∈ (0, 1). Since

lim m→+∞ [h m k m /n m ] n m = uv,
we have

n 2 m (h m -[h m k m /n m ])(k m -[h m k m /n m ]) 2h m k m (n m -h m )(n m -k m ) ∼ (u -uv)(v -uv) 2uv(1 -u)(1 -v) = 1 2 .
For the other factors in s nm (z (nm,hm) , z (nm,km) ) we use Stirling's formula, in the following version.

Let (n m ) m∈N , (A m ) m∈N , (B m ) m∈N with 0 ≤ B m ≤ A m . Assume n m → +∞, A m /n m → a, B m /n m → b as m → +∞ with 0 < b < a . Then A m B m ∼ 1 √ 2πn m × a b(a -b) × a a b b (a -b) a-b nm .

Thus,

n m h m -1 ∼ √ 2πn m × u(1 -u) × u u (1 -u) 1-u nm and k m [h m k m /n m ] ∼ 1 √ 2πn m × v uv(v -uv) × v v (uv) uv (v -uv) v-uv nm = 1 √ 2πn m × 1 uv(1 -u) × 1 u uv (1 -u) v-uv nm ; n m -k m h m -[h m k m /n m ] ∼ 1 √ 2πn m × 1 -v (u -uv)(1 -v -u + uv) × (1 -v) 1-v (u -uv) u-uv (1 -v -u + uv) 1-v-u+uv nm = 1 √ 2πn m × 1 u(1 -u)(1 -v) × 1 u u-uv (1 -u) 1-v-u+uv nm . Hence s nm (z (nm,hm) , z (nm,km) ) = n 2 m (h m -[h m k m /n m ])(k m -[h m k m /n m ]) 2h m k m (n m -h m )(n m -k m ) × n m h m -1 k m [h m k m /n m ] n m -k m h m -[h m k m /n m ] ∼ 1 2 √ 2πn m × u(1 -u) × 1 uv(1 -u) × 1 u(1 -u)(1 -v) × u u (1 -u) 1-u × 1 u uv (1 -u) v-uv × 1 u u-uv (1 -u) 1-v-u+uv nm = 1 2 2πuv(1 -u)(1 -v)n m .
Proof of Proposition 2.1. The upper bound in (2.1) is easily proved:

1 n! σ∈Sn 1 n n j=1 x σ(j) y j ≤ 1 n! σ∈Sn 1 n n j=1 |x σ(j) | • |y j | = 1 n n j=1 1 n! σ∈Sn |x σ(j) | |y j | = 1 n n j=1 |x| 1 • y j = |x| 1 • |y| 1 .
The lower bound follows from Lemma 2.2 and from the definition of c n . The asymptotic estimate (2.2) for c n is an easy consequence of Lemma 2.4, since max 0≤u,v≤1

uv(1 -u)(1 -v) = 1 16 .

Proof of Theorem 1.1

Let us prove the first assertion of the theorem. If a i = a j for some i = j, then α is not a generator of Q(β 1 , . . . , β n )/Q, since τ α = α for τ = (i, j). Assume now a 1 , . . . , a n pairwise distinct. Let σ be a permutation of S n such that σα = α. Then β

a 1 -a σ(1) 1 • • • β an-a σ(n) n = 1,
which implies, by [6, Lemma 1], that j → a j -a σ(j) is constant, say = κ. For l = o(σ) we have a 1 = a σ l (1) = • • • = a 1 + lκ, thus κ = 0 and ∀j, a σ(j) = a j . Since a 1 , . . . , a n are pairwise distinct, σ is the identity.

To prove 2), let x j = log |β j | ∈ R and remark that j x j = 0 and j |x j | = 2 log M (β). Moreover

2 n! log M (α) = 1 n! σ∈Sn n j=1 a j x σ(j) = 1 n! σ∈Sn n j=1
x σ(j) y j .

Inequality (2.1) of Proposition 2.1 then gives

c n 2 log M (β) n |y| 1 ≤ 1 n • 2 n! log M (α) ≤ 2 log M (β) n |y| 1 with c n ∼ 2 πn . Assertion 2) follows.
To prove the last assertion, we need the following elementary lemma:

Lemma. Let y ∈ H n with y j+1 -y j ≥ 1 for j = 1, . . . , n -1. Then

|y| 1 ≥ n -2 4 .
Proof. Let k be such that y k ≤ 0 < y k+1 . Then, for j = 1, . . . , k We can now prove 3). The integers y 1 , . . . , y n are pairwise distinct, since α is a generator of Q(β 1 , . . . , β n )/Q. Thus we may assume y j+1 -y j ≥ 1 for j = 1, . . . , n -1. From the lower bound for M (α) in 2) and from the lemma above we get:

y j ≤ y k -(k -j) ≤ -(k -j) while y j ≥ y k+1 + (j -k -1) ≥ j -k -1 for j = k +
c n |y| 1 ∼ 2 πn |y| 1 ≥ 2 πn n -2 4 ∼ n 8π .

Proof of Theorem 1.2

The proof of the first assertion of Theorem 1.2 is similar to the proof of the corresponding assertion of Theorem 1.1. which shows the upper bound for M (α) in 2). We now prove the lower bound. Let T n be the set of transposition of S n and, for τ = (i, j) ∈ T n , let A τ be the set of permutation of S n with support {1, . . . , n}\{i, j} and with no orbits of length 2.

A computation shows that n h=1 h! ≥ n n(n-1)/4 . Thus, by 2) and since |disc(β)| ≥ 1,

M (α) 1/n! ≥ 2

  1, . . . , n. Thus n|y| 1 = -

For

  any σ ∈ Gal(Q(α)/Q) we have |σα| ≤ |a| 1 max{|β 1 |, . . . , |β n |} ≤ |a| 1 M (β)

- 1 24 1 12n(n- 1 )

 111 |V (a)| • |disc(β)| 1/2 ≥ (n/4) 1/48 .

  1. Let β be an algebraic unit of degree n, with algebraic conjugates β 1 , . . . , β n . Assume that Q(β 1 , . . . , β n )/Q has Galois group S n . Let a 1 , . . . , a n ∈ Z

	and put
	α = β a 1 1 β a 2 2 • • • β an n .
	Then:
	1) α is a generator of Q(β 1 , . . . , β n )/Q if and only if a 1 , . . . , a n are
	pairwise distinct.
	2) Put y j

Lemma. For n ≥ 3, let Λ n be the set of permutation of S n with no orbits of length 1 or 2. Then |Λ n | ≥ n!/12.

Proof. The inclusion-exclusion principle shows that the set of permutation of S n without fixed points has cardinality

Notice that 1/3 ≤ κ(n)/n! ≤ 1/2 for n ≥ 3. Again by the inclusion-exclusion principle

By the lemma above (recall that n ≥ 5) the sets

Let τ , τ ∈ T n and let σ ∈ A τ , σ ∈ A τ . Then τ σ = σ (since the supports of τ σ σ are not the same). If τ = τ and σ = σ then τ σ = τ σ and σ = σ . Moreover, if τ = τ then again τ σ = τ σ (since σ has no orbits of length 2) and σ = σ (since the supports of σ and σ are not the same). Thus

From the two last displayed equations and from the inequality 1) .

We finally prove 3). Notice that the result is trivial if n ≤ 4, thus we assume n ≥ 5. The integers a 1 , . . . , a n are pairwise distinct, since α is a generator of Q(β 1 , . . . , β n )/Q. Thus we may assume a j+1 -a j ≥ 1 for j = 1, . . .


J M Díaz

O A González-Estrada

C I López

Biomechanical analysis of a cranial Patient Specific Implant on the interface with the bone using the Finite Element Method

Keywords: Biomechanics, Finite Element Method, Patient Specific Implant, Cranial Implant

New advance technologies based on reverse engineering, design and additive manufacturing, have expanded design capabilities for biomedical applications to include Patient Specific Implants (PSI). This change in design paradigms needs advanced tools to assess the mechanical performance of the product, and simulate the impact on the patient. In this work, we perform a structural analysis on the interface of a cranial PSI under static loading conditions. Based on those simulations, we have identified the regions with high stress and strain and checked the failure criteria both in the implant and the skull. We evaluate the quality of the design of the implant and determine their response given different materials, in order to ensure optimality of the final product to be manufactured.

I. INTRODUCTION

Biomechanics studies the mechanical behaviour of biological systems. It helps to provide solutions to the injuries the body is not capable of healing itself, i.e., cases in which natural recovery is not possible. This could be accomplished by locating orthopaedic implants where the tissue did not heal to help with the recovery of a patient [START_REF] Schmitt | Trauma biomechanics: introduction to accidental injury[END_REF].

The study of these biomedical implantable devices involves areas such as materials, medicine and engineering [START_REF] Juanola-Feliu | Market challenges facing academic research in commercializing nano-enabled implantable devices for in-vivo biomedical analysis[END_REF]. These disciplines cannot be mutually exclusive if we want to ensure that the final product is in agreement with the expectations of the design. Moreover, it is very important to evaluate the quality of the designs before manufacturing them, in order to find any possible existing flaws, optimise the implant prototype and guarantee the patient is getting a product that is going to enhance his quality of life [START_REF] Sumida | Custom-made titanium devices as membranes for bone augmentation in implant treatment: Clinical application and the comparison with conventional titanium mesh[END_REF] [START_REF] Schmitt | Trauma biomechanics: introduction to accidental injury[END_REF].

This testing process can be done through non-destructive tests by means of numerical simulation, using finite element analysis (FEA) [START_REF] Lantada | Handbook on Advanced Design and Manufacturing Technologies for Biomedical Devices[END_REF] to evaluate the implant-bone interface behaviour [START_REF] Harrysson | Custom-designed orthopedic implants evaluated using finite element analysis of patient-specific computed tomography data: femoral-component case study[END_REF]. The model for the skull-implant assembly is considered, subjected to given boundary conditions. The structural analysis is performed in order to evaluate the stress and strain fields. Failure criteria for the materials are also considered. It is expected that critical zones are located in the interface between the skull and the implant, and that will be the main subject of this work.

II. METHODS

A. Model identification and geometry definition

The models of the study are 3D skull reconstructions obtained with a multislice spiral CT scanner, Toshiba Aquilon, in DICOM format. Using the software Invesalius 3, with semiautomatic segmentation, they are converted to STL format. Later, they receive CAD treatment with Rhino®. Based on the reconstruction, two different implants are designed and all models imported to ANSYS 15.0.

For optimisation purposes, we removed some of the bone tissue, leaving only the bone portion of the cranial vault near to the defect, thus simplifying the model.

B. Material models

Mechanical properties of the materials to use in the numerical models are defined. The software used in the CAD reconstruction did not allow displaying the entire Hounsfield scale [START_REF] López | Experimental study to establishing the insertion torque on dental implant[END_REF]; therefore, the values used to represent the mechanical properties are taken from the literature. In the case of bone, Galicer [START_REF] Sánchez | Simulación del Bruxismo y su Efecto en la Estructura Ósea de la Mandíbula[END_REF] suggests that the value of the bone matrix density is . The bone elastic constants, Young modulus and Poisson ratio, are defined by the following equations [START_REF] Taylor | Determination of orthotropic bone elastic constants using FEA and modal analysis[END_REF]:

𝐸=1763 3,2 if >1,2 (1)
𝜐=0,32 if >1,2 (2)
By replacing the bone matrix density in (1), we have 𝐸 = 16725 MPa.

Table 1 shows the Polyether ether ketone (PEEK) and titanium alloy Ti6Al4V properties, suggested by Safi [START_REF] Safi | Investigation of the failure behavior of a cranial implant-skull model under different load conditions using FEM[END_REF]:

C. Contact type definition

Many numerical models are available to define the contact conditions. For this problem we have chosen a linear contact condition, defined using Multi-Point Constraint (MPC). Linear contact models are faster to solve and less computational demanding. Extension of this work to include nonlinear models is expected. The MPC condition defines a state of no separation between the implant and the bone tissue. MPC are recommended for non-conforming contact surfaces, as it is the case for the implant-bone assembly due to features of the design introduced in the CAD model.

D. Meshing

Meshing is a fundamental part in the FEA, since a coarse mesh can generate results with low accuracy, whereas a fine mesh can generate a high computational cost and poor numerical conditioning. Safi [START_REF] Safi | Investigation of the failure behavior of a cranial implant-skull model under different load conditions using FEM[END_REF] recommends the use of quadratic tetrahedral elements (SOLID 187).

Two nonconforming meshes were used: one for the implant and one for the skull. Both were generated with the Patch Independent method to overcome issues related to the STL features imported from the CAD. This process is ideal for models with complex geometries, with very sharp angles, as the algorithm ensures that the mesh is refined where needed, but maintains larger elements in areas where possible. Table 2 shows the parameters used to generate the meshes for implant 1, with a total of 529113 elements. The resulting meshes are shown in Figure 1. Table 3 shows the values of the parameters used to generate the meshes for the implant 2, with a total of 681643 elements. Figure 2 shows the mesh for the implant and bone tissue.

E. Boundary conditions and applied loads

In the skull part, displacement in all directions is restricted along the boundary defined by the cutting planes. For the implants, displacement is restricted in the x and y directions of the surfaces representing the perforations aimed for the location of the screws, thus simulating the screw clamping.

The applied load is a pressure of 1 MPa which will be applied on an circular area with a radius of 15 mm in the central zone of the implant, this is equivalent to statically applying 70 kg, as used previously in [START_REF] Safi | Investigation of the failure behavior of a cranial implant-skull model under different load conditions using FEM[END_REF]. Additionally, to simulate the pressure exerted by the screws, pressure is applied around the holes aimed to hold the screws. The value is defined in terms of the torque recommended by ASTM F-543 [START_REF] López | Experimental study to establishing the insertion torque on dental implant[END_REF] for HA 2.0 screws:

(3)

Where is the applied pressure, is the torque, is diameter of the head of the screw, giving a value of .

III. RESULTS

The quality of the implant designs will be assessed taking into account the stresses obtained for the implant and the skull, which should be below the material limits.

A. Skull with implant 1 using PEEK

The obtained results are analysed under the von Mises criterion, with an applied load P = 1 MPa. In Figure 3, we observe the obtained stress distribution. The greatest stress is located in the implant, on the region in contact with the fracture, and has a value of 96.291 MPa. This may be because the PSI implant contour adjusts tightly to the shape of the fracture, and the interference is large. It is recommended to change the CAD implant so that when the two parts come into contact interference is smaller.

B. Skull with implant 1 using titanium alloy (Ti6Al4V).

The simulation is performed changing the material assigned to the implant part to titanium alloy Ti6Al4V. The maximum stress for this model is 125.83 MPa. Compared to the stress for the PEEK model, this one is slightly larger due to the higher elastic modulus. Using Ti6Al4V, it can withstand a pressure P = 5 MPa, with a maximum stress of 629.16 MPa, below the elastic limit of . The maximum displacement for the Ti6Al4V.is 0.01543 mm, which compared with the PEEK model is much lower. This behaviour is ideal as it is intended that the implant is deformed as little as possible once it has been placed on the patient. However, notice that the more rigid Ti6Al4V has a stronger stress shielding effect, bad for bone resorption. Table 4 contains the results of the analysis for both implants. C. Skull with implant 2 using PEEK.

Considering the previous results we generate an enhanced implant model. We use the same boundary conditions for the implant 2. Figure 5 shows the stress distribution obtained for this model using PEEK. The maximum stress is located on the face that comes into contact with the bone, with 32.575 MPa. This implant model has a smoother contour in the area, which contact the fracture. Comparing the results obtained for the implant 1 using PEEK, and given that the behaviour of the stress of this model is linear, it can resist three times more load. Figure 6 shows the strain and displacements obtained, which are much lower compared to the implant 1 using PEEK. The maximum displacement has a value of 0.1448 mm. cranial-claib2014-PAPER_170_VF.docx

D. Skull with implant 2 using titanium alloy (Ti6Al4V).

Simulation is performed using the titanium alloy Ti6Al4V. Figure 7 shows the stress distribution for this case. Maximum stress is located in the holes provided to secure the implant with screws and have a value of 47 MPa. Different behaviour from the previous model can be observed because the titanium alloy is more rigid than PEEK and the contour that gets in contact with the bone is practically not deformed so the holes act as stress concentrators. Comparing this value with the obtained in implant 1 using Ti6Al4V, implant 2 resists a load three times higher. Figure 8 shows the strain and displacements respectively, with a max. displacement of 0.005 mm. IV. CONCLUSIONS For this work it was assumed that the bone behaves as an isotropic linear elastic material, which is an approximation to simplify the analysis. Meshing is a key part of the analysis and requires great effort in the preprocessing stage. Manual interaction is still required to produce good quality meshes for such type of problems. A linear contact model was considered, with prestress induced by the screws.

In the analysis performed, implant 2 is superior to implant 1 model, primarily because implant 2 is a redesign of the implant 1, in which the geometry is redefined to improve the contour. PEEK exhibits lower values of stress that the more rigid Ti6Al4V, which is in agreement with the known stress shielding effect of Titanium alloys. We show the feasibility of virtual prototyping for modelling PSI implants.

Extension of this work could include nonlinear models for contact conditions, and more complex material models, including varying bone density obtained from the scans.

Figure 1 :

 1 Figure 1: Final mesh for the model of the skull with implant 1.

Figure 2 :

 2 Figure 2: Final mesh for the skull with implant 2.

Figure 3 :

 3 Figure 3: Stress distribution in the implant 1 using PEEK.

Figure 4

 4 Figure4shows the strain and displacements, indicating that there are not large deformations. Structural stability is an important condition to ensure the health of the patient.

Figure 4 :

 4 Figure 4: Strain (a) and displacements (b) in implant 1 using PEEK.

Figure 5 :

 5 Figure 5: Stress distribution for the skull with implant 2 using PEEK.

Figure 6 :

 6 Figure 6: Strain (a) and (b) displacements in implant 2 using PEEK.

Figure 7 :

 7 Figure 7: Stress distribution in the implant 2 using Ti6Al4V.

Figure 8 :

 8 Figure 8: Strain (a) and (b) displacements in implant 2 using Ti6Al4V.

Table

  1. Young modulus E and Poisson ratio 𝜐 for PEEK and Ti6Al4V.

	Material	E [MPa]	𝜐
	PEEK	4000	0,4
	Ti6Al4V	110300	0,36

cranial-claib2014-PAPER_170_VF.docx

Table 2 .

 2 Parameters used for the mesh of implant 1.

	Parameter	Implant	Skull
	Element's maximum size	2 mm	5 mm
	Element's minimum size	0,75 mm	1,5 mm
	Number of cells across a gap	8	3
	Curvature angle	7,5°
	Smooth transition	Yes	Yes
	Number of elements	180303	167450

Table 3 .

 3 Parameters used for the mesh of implant 2.

	Parameters	Implant	Skull
	Element's maximum size	2 mm	3 mm
	Element's minimum size	0,75 mm	1 mm
	Number of cells across a gap	3	3
	Smooth transition	Yes	No
	Number of elements	199685	250001

Table 4 :

 4 Results obtained for the model of the skull with implants 1 and 2.

	Material	Max. stress	Max. strain	Max. displacement
		[MPa]		[mm]
	Implant 1
	PEEK	96,291	0,0281	0,3701
	Ti6Al4V Alloy	125,83	0,0013	0,0154
	Implant 2
	PEEK	32,575	0,0092	0,1448
	Ti6Al4V Alloy	47,034	0,0004	0,0050

ACKNOWLEDGMENT This work was supported in part by VIE-UIS under Code No. 1361. V. REFERENCES


Octavio Andrés González-Estrada

M León Uribe

Sundararajan Natarajan
email: snatarajan@iitm.ac.in

Juan José Ródenas
email: jjrodena@mcm.upv.es

Stephane Bordas
email: stephane.bordas@uni.lu

Oscar M León Uribe

Stéphane P A Bordas

Adaptivity and Error Control for the Polygonal Finite Element Using Recovery-based Error Estimates

Keywords: error estimation, adaptivity, recovery, polygonal finite element method

published or not. The documents may come


Juan S León Becerra

Octavio Andrés González-Estrada

Alberto Pertuz

Sebastián León Becerra

Damage evolution in composite pipes using a continuum damage mechanics formulation

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.


Ivan D Breslavsky
email: ivan.breslavskyi@mcgill.ca

Marco Amabili

Mathias Legrand

Annie Ruimi

DYNAMICS OF STATICALLY PRE-LOADED HUMAN AORTA

The response of human aortic segment under static and dynamic pressure is investigated. The aortic segment is modelled as threelayered hyperelastic circular cylindrical shell. Both geometrical and physical nonlinearities are taken into account. The effect of internal fluid is also taken into account. The results are presented in the form of pressure-deflection curves and frequency responses of forced vibrations.

PROBLEM DESCRIPTION

Static and dynamic analysis of circular cylindrical shell that models a segment of human aorta is carried out in this study. The shell is assumed to be three-layered, with each layer being hyperelastic, fiber-reinforced, viscoelastic material. Following the paper [START_REF] Holzapfel | Determination of material models for arterial walls from uniaxial extension tests and histological structure[END_REF], we chose a material model in form of combination of Neo-Hookean and Fung models. The strain energy density for this case has the following form:

𝑊 = 𝑊 𝑁𝐻 + 𝑊 𝐹 , (1)
𝑊 𝑁𝐻 = 𝐸 6 (𝐼 1 -3), (2)
𝑊 𝐹 = 𝐶(𝐸𝑥𝑝(𝑐 11 𝜀 𝑥 2 + 𝑐 12 𝜀 𝑥 𝜀 𝜃 + 𝑐 22 𝜀 𝜃 2 ) -1), (3)
where E is the Young's modulus of the shell material; 𝐼 1 is the first invariant of the right Cauchy-Green deformation tensor 𝐂 = 𝟐𝐄 + 𝐈, and 𝐄 is Green-Lagrange strain tensor. The material (1-3) is anisotropic material due to Fung term, which is included to describe anisotropic response of the fibers and their sharp increase in stiffness. Neo-Hookean term describes the isotropic part of the tissue response. The geometrical nonlinearity is described with the higher-order shear deformation theory [START_REF] Amabili | A new non-linear higher-order shear deformation theory for large-amplitude vibrations of laminated doubly curved shells[END_REF].

The viscoelasticity is taken into account through employment of Kelvin-Voigt model [START_REF] Lakes | Viscoelastic Materials[END_REF]. The shell is filled with quiescent fluid [START_REF] Amabili | Nonlinear vibrations and stability of shells and plates[END_REF]. Since we deal with large deflections and large strains, the displacement-dependent expression for pressure is used [START_REF] Amabili | Displacement dependent pressure load for finite deflection of shells and plates[END_REF]. Such load is proportional to the area of the deformed surface and is orthogonal to it.

Strain energy density in the form (1-3) is not polynomial in strain components, which significantly complicates the investigation of the dynamics. The analysis is simplified by expansion of strain energy density into truncated after order 4 series in strain components. Such model is able to describe the behavior of the vessel only locally, in the vicinity of a given configuration of interest around which the expansion is made. To reach highly deformed configurations, a sequence of successive local models has to be built.

NUMERICAL EXAMPLE

We consider the shell with geometrical parameters that are similar to the parameters of the straight thoracic aortic segment of human aorta. The geometrical parameters are [START_REF] Amabili | A three-layer model for buckling of a human aortic segment under specific flow-pressure conditions[END_REF]: radius 𝑅 = 0.01575 m, length 𝐿 = 0.126 m. Thicknesses for three layers are: for intima, the internal layer, ℎ 𝑖 = 0.00033 m, for media, the medium layer, ℎ 𝑚 = 0.00132 m, for adventitia, the external layer, ℎ 𝑎 = 0.00096 m. The material parameters are partially determined based on experimental data from [START_REF] Holzapfel | Determination of material models for arterial walls from uniaxial extension tests and histological structure[END_REF], partially adopted from [START_REF] Holzapfel | Determination of material models for arterial walls from uniaxial extension tests and histological structure[END_REF]. It is assumed that all three layers have the same density 𝜌 = 1200 kg/m 3 , the contained fluid density is taken 𝜌 𝐹 = 1060 kg/m 3 . The shell is simply supported at the ends. The boundary conditions, including the natural nonlinear boundary condition of zero stress resultant, are identically satisfied by the trigonometric functions of special form as shown in [START_REF] Amabili | A new non-linear higher-order shear deformation theory for large-amplitude vibrations of laminated doubly curved shells[END_REF]. The convergence study has been carried out and it is concluded that the 29 degrees-of-freedom system provides good approximation of the solution.

Initially, the static problem of inflation of the shell under the pressure load was studied. The pressure-defection curve is shown in Fig. 1. As can be seen, the stiffness drastically increases after certain strain threshold is reached.

Afterwards we analyzed free and forced vibrations around pre-loaded state of the entire artery as a composite shell and also single-layered shell to understand the contribution of different layers. For the three-layered artery model the configuration in the middle of physiological range (12.5 kPa pressure) was chosen for the expansion of the strain energy density around it. Vibrations with dominant first axisymmetric bending mode and with dominant mode with two half-waves in circumferential direction of the empty and fluid-filled shells are studied. The frequency response for the empty shell made of adventitial material is shown in Fig. 2. We can see that the nonlinearity is weak enough for this case, however, in the case of fluid-filled shell it is stronger.

CONCLUSIONS

The method for the analysis of nonlinear vibrations of arteries is presented. The method is able to deal with hyperelastic, viscoelastic, multilayered, fiber-reinforced shells that model blood vessels. The method is tested on the human aortic segment model. Under the static pressure initially soft shell becomes much stiffer with the pressure growth, which is the common feature for soft biological tissues. The dynamic nonlinearity is not very strong, but still nonlinear effects, like change in stiffening type, are observed.

Figure 1 .

 1 Figure 1. Static pressure-deflection curve. The abscissa shows central deflection.

Figure 2 .

 2 Figure 2. Frequency response of empty shell vibrations around pre-loaded configuration (adventitial material). Principal bending mode amplitude vs. frequency. Frequency response (thin line) and backbone curve (thick line).

ACKNOWLEDGEMENTS

The authors acknowledge the financial support of the NSERC Discovery Grant, Canada Research Chair and the Qatar National Research Fund NPRP 7-032-2-016.