Martin Rais
email: martin.rais@cmla.ens-cachan.fr

Jean-Michel Morel

Carole Thiebaud

Jean-Marc Delvit

Gabriele Facciolo

Improving wavefront sensing with a Shack-Hartmann device

In order to achieve higher resolutions, current earthobservation satellites use larger lightweight primary mirrors that can deform over time, impacting on image quality. We evaluated the possibility of compensating these deformations directly in the satellite by combining a deformable mirror with a Shack-Hartman wavefront sensor (SHWFS). The performance of the SHWFS depends entirely on the accuracy of the shift estimation algorithm employed, which should be computationally cheap to be executed on-board. We analyzed the problem of fast accurate shift estimation in this context and propose a new algorithm, based on a global optical flow method that estimates the shifts in linear time. In our experiments, our method proved to be more accurate and stable, as well as less sensitive to noise than all current state-of-theart methods, permitting a more precise on-board wavefront estimation.

Introduction

Adaptive optics (AO) is a well-known technology to sense and correct wavefront distortions. This technology is used in astronomy to produce sharper images from heavily aberrated wavefronts originated by atmospheric turbulence. This correction is usually performed through a deformable mirror which adapts to the measured wavefront correcting the distortion [START_REF] Wizinowich | First light adaptive optics images from the Keck II telescope: a new era of high angular resolution imagery[END_REF]. Since AO also allows to improve the performance of aberrated optical systems, it is widely used in several other contexts such as ophthalmology [START_REF] Wirth | Ophthalmic instrument having hartmann wavefront sensor with extended source[END_REF], microscopy [START_REF] Martin | Adaptive optics in microscopy[END_REF] and freespace laser communication systems [START_REF] Robert K Tyson | Bit-error rate for free-space adaptive optics laser communications[END_REF] among others [START_REF] Tyson | Principles of adaptive optics[END_REF].

A key component of an adaptive optics system is the wavefront sensing mechanism, i.e., the device used to precisely measure the distortion. A Shack-Hartmann wavefront sensor (SHWFS) is one such device. It uses an array of lenslets to measure the deformation of the incoming wavefront. The shift of each lenslet focal plane image is proportional to the mean slope of the wavefront in the subaperture onto this lenslet. It yields a discrete local approximation of the slope of the wavefront (Fig. 1). This deformation is usually measured by imaging a point source such as a star, and measuring the relative displacement between a reference image and all other subimages to compute the local gradient of the wavefront.

Recently, the community evaluated the possibility of cor-Figure 1: A Shack-Hartmann Wavefront Sensor measures the wavefront by computing the local shifts between the detected spots (in green) and the reference crosses (in black), which would occur if no deformation were present.

recting wavefront deformations on earth-observation satellites [START_REF] Poyneer | Scene-based shack-hartmann wavefront sensing: analysis and simulation[END_REF][START_REF] Escolle | Adapting large lightweight primary mirror to space active optics capabilities[END_REF][START_REF] Bonnefois | Comparative theoretical and experimental study of a Shack-Hartmann and a Phase Diversity Sensor, for high-precision wavefront sensing dedicated to Space Active Optics[END_REF] caused by the deformation of the primary mirror. In this setting, the problem of atmospheric turbulence is negligible. Indeed, in astronomical observations from the earth, the angle of view is extremely narrow. As a result, the light wavefront crosses a narrow solid angle of atmosphere and its perturbations due to turbulence have a great impact on image quality. In earth observation from a satellite, however, the angle of view is much larger, so the perturbations due to turbulence are relatively much smaller. However, the correction of optical aberration is becoming more and more important for high resolution earthobservation satellites. Indeed, in order to increase the spatial resolution of satellite images, i.e., its ground sample distance, a larger primary mirror is required to gain a higher angular resolution. Large mirrors must be thicker to avoid deformations which increases dramatically their weight and renders launching costs prohibitive. For this reason large yet lighter mirrors have to be considered. Their drawback is that time-varying deformations due to thermal effects and vibration severely deteriorate the image quality [START_REF] Peter C Chen | Advances in very lightweight composite mirror technology[END_REF]. To correct these deformations, a SHWFS device could be used to measure them by observing the earth, together with a deformable mirror to compensate for these deformations. As opposed to adaptive optics where the distortion has to be compensated in real-time due to the fast changes in the atmosphere, in our case the compensation for the mirror deformation is less frequent, roughly on timescales of seconds.

Coming from control theory, two schemes exist for performing wavefront correction. While in an open loop adaptive optics system the wavefront error is measured before it has been corrected by the wavefront corrector, in a closed loop system the measured wavefront deformation is the residual error after the correction of the previous estimation has been performed. The difference between both operational modes is important because in a closed loop system, the wavefront aberrations measured will be small, permitting to assume a maximum shift between images of up to a few pixels. In this article we assume a closed loop system.

Another important factor that affects the accuracy of wavefront sensing is the phenomenon of scintillation and phase anisoplanatism, which results in more complex patterns than simple global shifts between subapertures. The influence of this phenomenon was widely studied in the context of adaptive optics [START_REF] Poyneer | Experimental results for correlationbased wavefront sensing[END_REF][START_REF] Robert | Scintillation and phase anisoplanatism in shack-hartmann wavefront sensing[END_REF][START_REF] Védrenne | Shack-Hartmann wavefront estimation with extended sources: anisoplanatism influence[END_REF]. However, in remote sensing, because the phase aberration is produced at the telescope pupil, all parts of the image are affected by it in the same way, neglecting its incidence [START_REF] Poyneer | Experimental results for correlationbased wavefront sensing[END_REF].

Finally, once the shift estimation is performed, several methods can be used to reconstruct the wavefront from its local gradient estimations, namely the iterative zonal method [START_REF] William | Wave-front estimation from wavefront slope measurements[END_REF], the vector-matrix-multiply (VMM) method [START_REF] Herrmann | Least-squares wave front errors of minimum norm[END_REF] and the Fourier Transform Reconstruction (FTR) method [START_REF] Poyneer | Fast wave-front reconstruction in large adaptive optics systems with use of the fourier transform[END_REF]. The latter is recommended when the number of actuators is high, however, the VMM method obtains more accurate results under a SHWFS configuration using less than 12 × 12 subapertures [START_REF] Roopashree | Performance analysis of Fourier and Vector Matrix Multiply methods for phase reconstruction from slope measurements[END_REF].

As mentioned before, wavefront sensing in astronomy is usually performed using the stars. When observing the earth from space this task becomes more challenging. Extended scene vs point source observation. For an earth observation satellite, the SHWFS is used on extended scenes instead of point sources such as stars. This setup is called scene-based wave-front sensing (SBWFS), or extendedscene wave-front sensing (ESWFS). Because the scene is extended, a field stop has to be installed in front of the SHWFS, as seen in Fig. 1, so that the images given by the lenslet array do not overlap in the lenslet focal plane [START_REF] Michau | Shack-hartmann wavefront sensing with extended sources[END_REF]. This yields a grid of images, each one corresponding to one sub-aperture, which are shifted versions of the same scene. Accurately measuring these shifts permits to estimate the gradient of the wavefront. As a drawback, since we are dealing with landscapes larger than the captured image, achieving high accuracy on the shift estimation task gets challenging. Worse still, in most wavefront sensors, the extent of the source object normally reduces the contrast of the signal, thwarting accurate shift measurements [START_REF] Rousset | Wave-front sensors[END_REF]. Difference in subimage SNRs due to pupil occlusions. For long focal length telescopes commonly used on earthobservation satellites, the Korsch concept is the most common. The pupil of a Korsch telescope is generally occluded in the center by a secondary mirror. The arms used to hold this mirror also occlude the pupil. In these regions, the lenslets suffer a loss in the incoming signal, proportional to their percentage of occlusion. This configuration is depicted in Fig. 2a. An example of a SHWFS output in the CRT sensor is shown in Fig. 2b. From this figure, it can be seen how the SNR on the partially occluded lenslets is significantly lower. Limited on-board computational capacity. An important distinction when performing wavefront sensing from earth-observation satellites is its limited computational ca- pacity. Due to this constraint, several shift estimation methods proposed for SHWFS are not suitable on-board due to their high complexity. Unusable observations. Another difference when observing extended scenes, as opposed to the use of SHWFS with point sources, is the need to predict if the current scene permits to accurately estimate the wavefront aberration. Scenes such as clouds, sea or any textureless landscape can thwart all shift measurement methods, leading to poor wavefront estimation.

Contributions. Based on these differences, we present a new shift estimation method in the context of SHWFS used on extended scenes. Our contribution is threefold. We shall start with a review of the state-of-the-art on wavefront correction using SHWFS on extended scenes. Second, we propose the use of an iterative global optical flow method for shift estimation which presents several advantages over the conventional correlation methods. Third, we propose a fast and effective method for scene preselection that adds almost no further computational cost to the overall estimation using the proposed algorithm.

This article is organized as follows. In section 2 we review the state-of-the-art on shift estimation applied for SHWFS on extended scenes. In section 3 we present our optical flow method and its usage for wavefront correction from earthobservation satellites. In section 4 we focus on scene preselection methods. In section 5 we study the influence of the parameters on the proposed method, and thoroughly compare its performance with the state-of-the-art. We involve ground truth simulations provided by CNES and our own simulator. Section 6 is a conclusion.

State-of-the-art review

Since the Shack-Hartmann wavefront sensor was introduced in the late 1960s [START_REF] Ben | History and principles of shackhartmann wavefront sensing[END_REF], several algorithms were proposed to estimate the shifts using point sources such as stars. However, only few authors have studied the problem when the source is extended, as when observing the earth from space. Broadly, they can be categorized in correlation methods working in the spatial [START_REF] Michau | Wavefront sensing from extended sources[END_REF][START_REF] Michau | Shack-hartmann wavefront sensing with extended sources[END_REF][START_REF] Mats | Evaluation of image-shift measurement algorithms for solar shack-hartmann wavefront sensors[END_REF] or in the frequency domain [START_REF] Poyneer | Scene-based shack-hartmann wavefront sensing: analysis and simulation[END_REF][START_REF] Sidick | An Adaptive cross-correlation algorithm for extended scene Shack-Hartmann wavefront sensing[END_REF][START_REF] Sidick | Adaptive periodic-correlation algorithm for extended scene shack-hartmann wavefront sensing[END_REF][START_REF] Per Knutsson | Extended object wavefront sensing based on the correlation spectrum phase[END_REF][START_REF] Mats | Evaluation of image-shift measurement algorithms for solar shack-hartmann wavefront sensors[END_REF], phase correlation methods that estimate the shift directly in the frequency domain [START_REF] Per Knutsson | Extended object wavefront sensing based on the correlation spectrum phase[END_REF][START_REF] Stone | A fast direct fourier-based algorithm for subpixel registration of images[END_REF], iterative methods [START_REF] Sidick | An Adaptive cross-correlation algorithm for extended scene Shack-Hartmann wavefront sensing[END_REF][START_REF] Sidick | Adaptive periodic-correlation algorithm for extended scene shack-hartmann wavefront sensing[END_REF] that improve on previous approaches by progressively estimating the shift, and maximum likelihood (ML) approaches [START_REF] Gratadour | Sub-pixel image registration with a maximum likelihood estimator[END_REF] that incorporate a specific noise model to the problem and compute ML estimate as the solution of an optimization problem.

Correlation-based Methods

These methods compute a correlation score on a grid C(i, j) and interpolate it to determine the subpixel location of the peak. The methods mainly differ in the choices of correlation score and interpolation strategy. Spatial domain correlation-based methods. Michau et al. [START_REF] Michau | Wavefront sensing from extended sources[END_REF] were among the first to propose an experimental implementation for using a SHWFS on extended sources larger than the wavefront sensor field of view. Their method computes the discrete cross-correlation between the images and a reference subimage chosen from the central region of the wavefront sensor image pattern. To estimate the subpixel shifts, the correlation peak location is computed as the centroid of the pixels with a correlation higher than half the maximal observed correlation.

Löfdahl [START_REF] Mats | Evaluation of image-shift measurement algorithms for solar shack-hartmann wavefront sensors[END_REF] tested several shift estimation algorithms for Shack-Hartmann wavefront sensors observing the sun. By testing with several possible sources of errors such as noise, blur and bias mismatch, he evaluated five different correlation-based methods to obtain the correlation score C(i, j). Among the evaluated correlation algorithms, the best all-around performer proved to be the classical least squares approach or squared difference function (SDF). Another proposed correlation score is the Covariance Function in the Image domain (CFI) using trend-corrected versions of both images. Indeed, Smithson & Tarbell [START_REF] Tarbell | A Simple Image Motion Compensation System for Solar Observations[END_REF] showed that a linear trend in intensity shifts the covariance peak from the correct position, so a fitted plane has to be subtracted from both images beforehand. For the SHWFS case, this reduces to simply subtracting the mean value for each image. They also tested other two methods based on the absolute difference between both images (ADF), and its square (ADF 2 ) in order to better locate the minimum at the sub-pixel level.

To achieve subpixel precision, they evaluated four interpolation strategies to look for the minimum value on the correlation grid C(i min , j min ). The four algorithms they evaluated can be described as fitting a conic section to the 3×3-element submatrix s centered in the sample minimum C(i min , j min ). The evaluated algorithms differ on the number of pixels used and on whether the fitting is done in 2D or in each dimension separately. The 1D quadratic interpolation (1QI) fits a parabola in each dimension independently, while the 2D quadratic interpolation (2QI) fits the conic to the 3 × 3 neighborhood. They also evaluated a 1D Least Square (1LS) method on both dimensions, which averages the three values on one dimension and then applies least squares on the other, and a 2D Least Square (2LS) variant which includes the corner values in the procedure.

Frequency domain correlation-based methods.

Poyneer [START_REF] Poyneer | Scene-based shack-hartmann wavefront sensing: analysis and simulation[END_REF] studied the wavefront estimation problem using a SHWFS by observing the Earth from space using lightweight optics. By assuming periodicity on the input images, the author points out that minimizing the MSE between both images becomes equivalent to maximizing their periodic convolution, which is efficiently computed in the frequency domain using the cross-correlation theorem. To get subpixel precision the maximum at integer coordinates is refined by independently fitting a parabola on each dimension.

Löfdahl also evaluated a frequency-domain method [START_REF] Mats | Evaluation of image-shift measurement algorithms for solar shack-hartmann wavefront sensors[END_REF] which, similarly to [START_REF] Poyneer | Scene-based shack-hartmann wavefront sensing: analysis and simulation[END_REF], computes the covariance in the Fourier domain (CFF). However, both images are previously normalized to zero mean and windowed with a 2D Hamming window or a flat-top window to avoid ringing caused by the periodization. Again, the subpixel maximum is obtained by fitting a parabola to the grid.

Phase Correlation

The phase-correlation method was widely studied in the image processing domain [START_REF] Stone | A fast direct fourier-based algorithm for subpixel registration of images[END_REF][START_REF] Foroosh | Extension of phase correlation to subpixel registration[END_REF]. Let I(i, j) be an M ×N image, due to the Fourier shift theorem we know that

F I(i-∆ x , j -∆ y ) = Î(u, v)exp -j2π u∆ x M + v∆ y N (1)
then by computing the cross-power spectrum between both images and extracting the phase for each frequency, the matrix φ, called the phase correlation matrix, is given by

arg Ĉ (u, v) = φ(u, v) = 2π u∆ x M + v∆ y N . (2)
The shifts can therefore be computed directly in the frequency domain by fitting a plane passing through the origin of φ(u, v). Due to aliasing, some of these frequencies may be corrupted and distort the shift estimation. To avoid this problem Knutsson et al. [START_REF] Per Knutsson | Extended object wavefront sensing based on the correlation spectrum phase[END_REF] discard most of the corrupted frequencies from φ retaining only the two or four lowest frequencies. This estimate has the advantage of being the least sensitive to aliasing [START_REF] Stone | A fast direct fourier-based algorithm for subpixel registration of images[END_REF] hence the most reliable. Moreover assuming small shifts (smaller than half pixel) no phase unwrapping is required [START_REF] Dennis | Two-dimensional phase unwrapping: theory, algorithms, and software[END_REF]. The accuracy of these methods however, suffers considerably in low SNR situations [START_REF] Poyneer | Scene-based shack-hartmann wavefront sensing: analysis and simulation[END_REF].

Iterated Estimation

Since correlation-based shift estimation has a bias proportional to the shift magnitude [START_REF] Pham | Performance of optimal registration estimators[END_REF], then compensating the shift by resampling one of the images and iterating should progressively reduce this bias (see section 3). Sidick et al. [START_REF] Sidick | An Adaptive cross-correlation algorithm for extended scene Shack-Hartmann wavefront sensing[END_REF] proposed Adaptive Cross-Correlation (ACC) that estimates the shift using an approach similar to Knutsson [START_REF] Per Knutsson | Extended object wavefront sensing based on the correlation spectrum phase[END_REF] and uses this estimation to resample the second image in the frequency domain to iterate the procedure. Due to potential ringing artifacts on the image boundaries after the resampling procedure, the shift is estimated using the central part of both images with size N/2, involving one fourth of the image pixels in the computation. To make the shift estimation more accurate, it uses eight frequency components φ(n, m) with 0 ≤ m, n ≤ 2 to perform the least-squares fitting, excluding the center. This procedure is made iterative by accumulating both u and v shift estimates, until a predefined amount of iterations is reached or until ∆ s = [∆ 2

x + ∆ 2 y ] 1/2 , the increase of both shift estimates, is lower than a predefined tolerance value. This technique is then applied for every sub-aperture I k (i, j).

This algorithm was evaluated on simulated scenes by varying the shift between them [START_REF] Sidick | An Adaptive cross-correlation algorithm for extended scene Shack-Hartmann wavefront sensing[END_REF] and by using a real SHWFS configuration testbed [START_REF] Sidick | Adaptive Cross-Correlation Algorithm and Experiment of Extended Scene Shack-Hartmann Wavefront Sensing[END_REF] produced in the Jet Propulsion Laboratory at NASA [START_REF] Rhonda M Morgan | Testbed for extended-scene shack-hartmann and phase retrieval wavefront sensing[END_REF]. By relaxing the tolerance parameter of their iterative scheme, it achieves errors of the order of 0.05 pixels by using between 3 and 6 iterations.

This method presents two problems. Since it discards half of each image to perform the shift estimation, it usually fails when image contents appear away from the center. Also, as in [START_REF] Per Knutsson | Extended object wavefront sensing based on the correlation spectrum phase[END_REF], the shift estimation step of the method suffers severely in low SNR scenarios. Due to this reason, the authors updated the method in [START_REF] Sidick | Adaptive periodic-correlation algorithm for extended scene shack-hartmann wavefront sensing[END_REF] by replacing this step with the periodic-correlation technique of Poyneer [START_REF] Poyneer | Scene-based shack-hartmann wavefront sensing: analysis and simulation[END_REF]. This method, named Adaptive Periodic Correlation (APC), is more tolerant to noise than ACC.

Maximum Likelihood Estimator

Gratadour et al. [START_REF] Gratadour | Sub-pixel image registration with a maximum likelihood estimator[END_REF] studied the performance of the MLE within an adaptive optics context. Assuming that the two images I 0 and I 1 are unknown and contaminated with Gaussian noise their method reduces to minimizing

k (4σ 2 (k)) -1 |I 1 (k) -[I 0 (x -x 1 )] III (k)| 2 , (3)
with respect to the shift x 1 , where σ 2 (k) is the noise variance at pixel k and [] III a sampling operator that performs Fourier resampling, which is done by computing the inverse FFT of (1). To guarantee that the images can be resampled in the Fourier domain they are pre-filtered using a low-pass filter. This minimization is performed using an a conjugate gradient method. Although it improves over the typical crosscorrelation, it has high computational cost prohibiting its implementation on satellites [START_REF] Michau | Shack-hartmann wavefront sensing with extended sources[END_REF][START_REF] Poyneer | Scene-based shack-hartmann wavefront sensing: analysis and simulation[END_REF].

3 Accurate shift estimation using optical flow in the context of a SHWFS

In order to estimate the shift between two sub-apertures, we propose to use a gradient-based shift estimator (GBSE) based on the optical flow equation. The idea behind this methodology, as proposed originally by Lucas and Kanade [START_REF] Lucas | An iterative image registration technique with an application to stereo vision[END_REF], is to relate the difference between two successive frames to the spatial intensity gradient of the first image. Given the two images I 2 (x, y) = I 1 (x -v x , y -v y ) where v x and v y are the unknown shift coefficients, the first order Taylor approximation yields

I 1 (x, y) -I 2 (x, y) v x ∂I 1 (x, y) ∂x + v y ∂I 1 (x, y) ∂y . (4)
To estimate the global optical flow between I 1 and I 2 , the Lucas-Kanade algorithm assumes a constant flow for every pixel, which allows the construction of an overdetermined system of constraint equations Av = b, where A is composed of spatial intensity derivatives and b has temporal derivatives.

Emulating Simoncelli [START_REF] Simoncelli | Design of multi-dimensional derivative filters[END_REF], in order to increase the accuracy of the method by minimizing noise or aliasing influence, we looked for two kernel functions: an anti-symmetric kernel d to estimate the image gradients and a symmetric kernel k to prefilter the images. Using both kernels, matrix A and vector b become

A =    (d x * I 1 )(p 1 ) (d y * I 1 )(p 1 ) . . . . . . (d x * I 1 )(p n ) (d y * I 1 )(p n )    , b =    (k * (I 1 -I 2 ))(p 1 ) . . . (k * (I 1 -I 2 ))(p n )   
(5) where p i with i = 1 . . . n represents the ith pixel and n the image size. The shift is obtained by Moore-Penrose pseudoinversion

v x v y = I 2 x I x I y I x I y I 2 y A T A -1 I t I x I t I y A T b (6)
where I x , I y stands for d x * I 1 and d y * I 1 respectively, and

I t = k * (I 1 -I 2 )
is the derivative over time.

It is no coincidence that the results of the method depend on the inversion of this second moment matrix. As will be shown later, the determinant of this matrix is crucial for determining the limits on the estimation performance. This matrix will be used to discard ill-posed cases before actually performing the shift estimation. The classic rejection case is when the gradient is mostly oriented in a single direction. This unsolvable situation is known as the aperture problem. Its detection by using the Cramer-Rao lower bounds is detailed in section 4.

It is important to remark that by centering the Taylor development at 0 and taking up to the first order term, the method gets systematically biased and becomes less precise as the shift gets higher, i.e., estimated shifts larger than one pixel would not be correctly estimated. Furthermore, the noise in the input images is completely ignored by the algorithm, which also impacts on its performance. Both these bias sources are studied in detail in [START_REF] Robinson | Fundamental performance limits in image registration[END_REF][START_REF] Pham | Performance of optimal registration estimators[END_REF] and will be addressed next.

Iterative gradient-based shift estimator

Instead of dealing with the bias explicitly, it was shown in Pham et al. [START_REF] Pham | Performance of optimal registration estimators[END_REF] that both bias sources depend linearly on the shift magnitude. This justifies the use of an iterative method, which is able to significantly reduce the bias, provided an appropriate resampling method is used. Algorithm 1 performs k iterations, computing the shift v[i] by solving (6) (f indshif t), and iterates by reinterpolating the original image I 2 (Resample) using the total accumulated shift w. It is a variant of the Lucas-Kanade algorithm [START_REF] Lucas | An iterative image registration technique with an application to stereo vision[END_REF] known as Inverse Compositional Algorithm [START_REF] Baker | Lucas-kanade 20 years on: A unifying framework[END_REF], in which the gradients of the reference image are static throughout the iterations, so (A T A) -1 is only computed once.

In the next sections we will see that the gradient discretization used by f indshif t and the choice of interpolation method for image resampling are crucial for the precision of the GBSE method.

Algorithm 1 Iterative GBSE method.

1: procedure ILK(I 1 , I 2 )
Receives a pair of images

2: i ← 0; I 2 [0] ← I 2 ; w ← 0; 3: while i ≤ k do 4: v[i] ← f indshif t(I 1 , I 2 [i]) (6) 5:
w ← w + v[i] Accumulate total shift 6:

I 2 [i + 1] ← Resample(I 2 , -w) Use the original I 2 7: i ← i + 1 8:
return w Return accumulated shift

Gradient computation and image prefiltering

Gradient computation for optical flow methods has been thoroughly studied [START_REF] Simoncelli | Design of multi-dimensional derivative filters[END_REF][START_REF] Elad | Optimal filters for gradient-based motion estimation[END_REF][START_REF] Robinson | Bias minimizing filter design for gradient-based image registration[END_REF]. Since our shift estimation is focused on accuracy and the image size on SHWFS is typically small (usually smaller than 50 × 50 pixels), using a large kernel for computing image derivatives implies loosing valuable boundary values. This constrains the kernel to be compact yet precise and robust to noise. In fact, the impact on the accuracy and the robustness to noise of the gradient computation is a key factor to the final performance of the GBSE method.

Because the method should be computationally fast, only simple schemes could be afforded. A straightforward candidate is the centered differences kernel [1, 0, -1], however since the central pixel is ignored in the computation, its precision could be improved by taking contiguous pixels. For this reason, a backward difference kernel [1, -1] would seem more appropriate. However, the corresponding center of this derivative differs for each component as seen from the middle image of Fig. 3.

A more precise gradient estimation method, illustrated on the right of Fig. 3, computes the derivatives by performing convolution with d x and d y given by

d x = 1/2 -1/2 1/2 -1/2 d y = 1/2 1/2 -1/2 -1/2 , (7)
and computes the vector b in (5) by prefiltering both I 1 and I 2 using the half pixel bilinear shift kernel

k = 1/4 1/4 1/4 1/4 . (8)
This gradient estimation procedure called hypomode [START_REF] Rouge | Fine stereoscopic image matching and dedicated instrument having a low stereoscopic coefficient[END_REF], despite being simplistic usually improves the accuracy obtained by GBSE methods using finite difference gradient estimation. This is because it slightly blurs the input images, which alleviates both aliasing and noise. We also considered the normalized 2D Gaussian smoothing kernel k = g(x, y, σ g ) with standard deviation σ g , and its derivatives d x = -x σ 2 g g(x, y, σ g ) and d y = -y σ 2 g g(x, y, σ g ). The parameter σ g controls the amount of blur, thus using a small value would be less tolerant to noise, while a high value could remove textures useful for the shift estimation. In our experiments, we evaluated using σ g = {0.3, 0.6, 1} with supports of 3, 5 and 7 pixels respectively.

As seen from both Fig. 4a and Fig. 4b, the fastest and most accurate gradient estimator for small shifts (≈ 0.1) is the hypomode. Therefore it should be considered when the underlying deformations are small enough (which could happen on a closed loop system with frequent mirror corrections). However, if the maximum shift caused by the wavefront aberration is larger than half a pixel, computing the Gaussian derivatives with σ g = 0.6 offers the best balance between accuracy and tolerance against noise, as seen from figures 4c, 4d, 4e and 4f. In fact, using σ = 0.6 gives better results than using σ = 0.3. This is because a higher σ removes more noise from the derivatives, and blurring the images impacts less on the accuracy than neglecting the noise. Finally, it should be noted that although taking σ = 1 usually offers poor accuracy with small image sizes, its tolerance to noise becomes a decisive factor as the noise increases (as seen in figures 4e and 4f).

Interpolation methods for image resampling

In order to iterate the method, the second image needs to be shifted, as indicated in step 6 of Algorithm 1. To this end, five different interpolation methods were evaluated, namely bilinear, bicubic [START_REF] Keys | Cubic convolution interpolation for digital image processing[END_REF] and cubic spline [START_REF] De | A practical guide to splines[END_REF], together with resampling using the Fourier shift theorem [START_REF] Bernstein | Resampling Images in Fourier Domain[END_REF], which is evaluated with and without image periodization. Image periodization prevents the generation of ringing artifacts due to the discontinuities at the image boundaries. It amounts to resampling with DFT a mirrored version of the image which has no discontinuities when periodization is assumed (see Fig. 5). The choice of interpolation method is critical for the overall method accuracy. However precision comes at a higher processing cost. We evaluated the different interpolations by simulating a landscape with noise of σ n = 80, displacements between subapertures of up to 2 pixels, and computing gradients by using Gaussian derivatives with σ g = 0.6. For each method we computed the average error and standard deviation over 200 independent random realizations, using the GBSE method with 2 and 3 iterations. Results in Fig. 6 show a non-negligible performance difference between the methods. In terms of average error, both DFT-based methods achieve the best results, however no significant improvement is obtained by performing periodization, although a minimal reduction in the standard deviation is observed. The splines interpolation yields slightly less accurate results, however it should be considered if fast DFT hardware cannot be installed on board. Finally, both bicubic and bilinear interpolation perform significantly worse than the rest, although bicubic interpolation results had less variability with respect to other approaches. Lastly using three iterations does not improve considerably over using two, not justifying the increase in computational cost. However, although not shown in this experiment, this difference gets higher as the noise increases, which justifies its use.

Image intensities equalization

One requirement of the GBSE algorithm is to work on images with similar intensities. Since occluded sub-apertures receive less light, their intensities differ from the reference image. To this end, all subimages are equalized by normalizing their mean with the highest mean among all subapertures. When the point spread function (PSF) of each sub-aperture is known, another equalization strategy before shift estimation between a subaperture and the reference is to convolve both images with the PSF of the other one [START_REF] Bonnefois | Comparative theoretical and experimental study of a Shack-Hartmann and a Phase Diversity Sensor, for high-precision wavefront sensing dedicated to Space Active Optics[END_REF]. Our focus being on the shift estimation method, we chose to perform the equalization by the simple procedure explained above. It is computationally cheaper, yet obtains excellent results.

Multiscale implementation

The GBSE algorithm can be easily adapted to work in an open loop environment where aberrations are potentially larger, leading to larger shifts between subapertures. As mentioned before, if the shift is larger than one pixel, the GBSE method fails. However, by building a pyramid representation of the input images, (6) can be applied on each scale to estimate the shift between images, and this estimated shift can in turn be used to resample the second image on the following level of the pyramid. If more accuracy is desired, Algorithm 1 can be used to better estimate the shift at each scale. We compute the pyramid using an approximate dyadic Gaussian pyramid [START_REF] Peter | The laplacian pyramid as a compact image code[END_REF]. Starting from the coarse image at scale n > 1, the method is summarized in Algorithm 2.

Algorithm 2 Multiscale GBSE method.

1: procedure MSSE(I 1 , I 2 )
Receives a pair of images 2:

I 1...n 1 ← BuildP yramid(I 1 , n) Burt&Adelson's Gaussian
Pyramid [START_REF] Peter | The laplacian pyramid as a compact image code[END_REF] 3:

I 1...n 2 ← BuildP yramid(I 2 , n) i.e.,
while i > 1 do 6: v(i) ← ILK(I i 1 , I i 2 ) or f indshif t(I i 1 , I i 2 ) Alg. 1 or (6) 7:
w ← 2w + 2v(i) 8:

I i-1 2 ← Resample(I i-1
2 , -w)

9: i ← i -1 10: v(1) = f indshif t(I 1 1 , I 1
2 )

11:

return w + v(1)
Return accumulated shift

Scene preselection and robustness estimation

Before evaluating the deformation of the incoming wavefront, it is crucial to determine if the current landscape, i.e., what the satellite is observing at the moment, is suitable for performing wavefront sensing. This problem has been studied by Sidick et al. [START_REF] Sidick | Adaptive Cross-Correlation Algorithm and Experiment of Extended Scene Shack-Hartmann Wavefront Sensing[END_REF] and by Poyneer [START_REF] Poyneer | Scene-based shack-hartmann wavefront sensing: analysis and simulation[END_REF] who proposed several fail-safe criteria to determine if a SH image is acceptable.

In this context, we propose to perform two distinct validations to ensure an accurate wavefront correction. While the first one is based on calculating the lowest theoretical achievable error through the Cramer-Rao lower bound (CRLB) given by the noise and the image contents, the second validation, called the eigenratio test, measures the gradient distribution along the image through the ratio of the eigenvalues of the structure tensor. Both validations are complementary and should be verified together.

CRLB test

The best achievable accuracy, given the current landscape and the underlying noise, can be estimated by the Cramer-Rao lower bound. In the last decade, three main works have addressed the calculation of the error on the estimated shift caused by the noise. Robinson and Milanfar [31] used the Cramer-Rao lower bound to derive a performance limit for image registration. The CRLB imposes a lower limit on the mean squared error for any estimate of a deterministic parameter. Pham et al. [START_REF] Pham | Performance of optimal registration estimators[END_REF] continued on the same idea to derive a theoretical limit for image registration, followed by a study of the bias for the gradient based 2D shift estimation. Recently, Aguerrebere et al. [START_REF] Aguerrebere | Fundamental limits in multi-image alignment[END_REF] performed an in depth study of performance limits within a multi-image alignment context and derived several Cramer-Rao lower bounds depending on the conditioning of the problem.

By assuming that the noise is independent, homoscedastic, Gaussian and with the same variance σ 2 n on both images, and denoting by S the set of all pixels in the image, the Cramer-Rao bound for any unbiased 2D shift estimator is

var(v x ) ≥ σ 2 n S I 2 y Det var(v y ) ≥ σ 2 n S I 2 x Det (9)
where Det = S I 2

x S I 2 y -S I x I y 2 . We can therefore define a parameter ∆ CRLB , the maximum allowed error in pixels, that determines whether the current landscape is acceptable for performing wavefront estimation, by verifying

(var(v x ) + var(v y )) 1/2 ≤ ∆ CRLB . (10)
Yet, I x , I y represent the gradient obtained from the unknown noiseless image I. Thus this bound is only useful for a theoretical study. In practice, however, the required values can be approximated using the method of [START_REF] Pham | Bidirectinal bias correction for gradient-based shift estimation[END_REF], where they estimate the second moment matrix of image I 1 from noisy versions Ĩ1 and Ĩ2 . Another possibility would be to approximate these sums by computing the expected value using the derivative definition and a noise variance estimation. For example set I x (i, j) = I(i + 1, j) -I(i, j) and let Ĩ = I + n be the observed noisy image where n has distribution N (0, σ 2 n ). Since all n(i, j) are independent, by the law of large numbers we have

S Ĩ2 x = S (I x +n x ) 2 = S I 2 x + S I x n x + S n 2 x S I 2 x +2|S|σ 2 n .
It follows that we can estimate S I 2 x from the noisy image Ĩ by

S I 2 x S Ĩ2 x -2|S|σ 2 n . (11)
The other terms can be computed using similar calculations. The CRLB defined in [START_REF] De | A practical guide to splines[END_REF] expresses a direct relationship between the accuracy of the estimation and the SNR ratio, measured as a ratio of the noise to the square of the gradient. Furthermore, it follows from the form of the denominator of these expressions that if there is a strong correlation between I x and I y , the vertical and horizontal partial derivatives, then Det will be zero or very close to zero. The formulas show that this entails a high variance for the estimation of the translation. This is the so called Aperture Problem [START_REF] Shimojo | Occlusion and the solution to the aperture problem for motion[END_REF]. In that case the true motion is irrecoverable.

It must be noted that the Cramer-Rao bounds are derived for an unbiased 2D shift estimation. For the case of biased estimators, the bound is even higher. Also, the hypothesis made about the noise being white Gaussian and homoscedastic does not hold in a SHWFS context used from earth-observation satellites. However, the noise model for the sensor is usually known beforehand, therefore allowing its variance to be stabilized by applying a variance stabilization transform (VST) such as the well-known Anscombe transform [START_REF] Anscombe | The transformation of poisson, binomial and negative-binomial data[END_REF].

Finally, for the case of SHWFS, since shifts between images are independent, the limit on the registration accuracy is imposed by estimating the shift between two frames.

Eigenratio test As noted by [START_REF] Aguerrebere | Fundamental limits in multi-image alignment[END_REF], the CRLB is less precise where it is most needed, namely when the signal is dominated by noise. Indeed, in that case the numerator and denominator of (9) cannot be estimated reliably, being both the difference of equivalent terms. This is why we will be forced to recur to more robust safety computations.

We propose to rapidly discard a landscape for a shift estimation by involving the eigenratio score of the second moment matrix τ = A T A of (6). The matrix τ is also the structure tensor with unit weights associated to the image I 1 . The eigenvalues λ 1 , λ 2 of τ and their corresponding eigenvectors e 1 , e 2 summarize the distribution of the image gradient ∇I = (I x , I y ). Namely, if λ 1 > λ 2 , then e 1 (or -e 1 ) is the direction that is maximally aligned with the gradient within the image. More generally, the value of λ k , for k ∈ {1, 2}, is the average of the square of the directional derivative of I along e k . The relative discrepancy between the two eigenvalues is an indicator of the degree of anisotropy of the gradient in the image, measuring how strongly it is biased towards a particular direction (and its opposite). Then by calculating the ratio λ 2 /λ 1 , we obtain a number between 0 and 1, characterizing the dominance of a particular direction for the gradients of the image. Empirically it was verified that a ratio λ 2 /λ 1 < 0.2 effectively degrades the shift estimation task, and this metric could be used to discard badly conditioned situations.

To compute this ratio, since τ = A T A is positive definite and symmetric, then both of its eigenvalues λ 1 and λ 2 are real and non-negative and its computation is straightforward.

Results

Scene pre-selection evaluation

Figures 7a and7b display two examples of SH landscapes that do not validate the eigenratio test. In these images there is a strong dominant gradient direction that complicates the shift estimation. The aperture problem is visually obvious in both cases. Figures 7c and7d present a high CRLB. The first one corresponds to a project where the images are composed by only noise, while the second example has a higher SNR, however the noise is still considerable. Finally, two examples of valid landscapes are shown in figures 7e and 7f which were validated by both proposed tests.

Comparison with state-of-the-art methods

5.2.1 The CNES simulated database CNES (French space agency) provided realistic simulations of several SHWFS landscapes obtained from earth-observing satellites with their corresponding ground truth. The provided images were 37 × 37 pixels from a 12 × 12 grid following the occlusion schema of Fig. 2a. The simulated wavefront aberrations translated into displacements no larger than half a pixel. For each provided landscape, three different SNR settings were simulated. Fig. 8 displays some example landscapes.

Every described state-of-the-art algorithm and our proposed method using different interpolation/gradient estimation approaches was tested against the whole dataset. In table 1 the errors in pixels are displayed for the best performing methods. These are the method from Michau et al. [START_REF] Michau | Wavefront sensing from extended sources[END_REF], the SDF method performing 2D quadratic interpolation (SDF-2QI) from Lofdahl [START_REF] Mats | Evaluation of image-shift measurement algorithms for solar shack-hartmann wavefront sensors[END_REF], the periodic convolution approach from Poyneer [START_REF] Poyneer | Scene-based shack-hartmann wavefront sensing: analysis and simulation[END_REF] and the iterative APC method from Sidick et al. [START_REF] Sidick | Adaptive periodic-correlation algorithm for extended scene shack-hartmann wavefront sensing[END_REF]. For the proposed GBSE approach, three variants are considered with increasing computational cost. A fourth variant using the Gaussian derivative estimation method was included for comparison purposes.

The errors obtained for figures 8.1 and 8.2 are close to the theoretical limits predicted by their CRLB. On the other side, the error obtained with Fig. 8.3 and its CRLB implies that a strong improvement is possible. Line 4, corresponding to Fig. 8.4, is an example where every method fails to achieve accurate results, possibly due to the noise and the lack of texture, which is revealed by a lower bound considerably higher than for most other cases. Fig. 8.5 shows an example of a landscape where the Poyneer method [START_REF] Poyneer | Scene-based shack-hartmann wavefront sensing: analysis and simulation[END_REF] as well as PCC [START_REF] Sidick | Adaptive periodic-correlation algorithm for extended scene shack-hartmann wavefront sensing[END_REF] are clearly improved by the proposed approach. Fig. 8.6 proves that the aperture problem hinders most shift estimation methods, as evidenced by their resulting errors. a Rows in green represent landscapes considered valid (CRLB < 0.02 and ER > 0.2), while invalid landscapes are shown in red. Bold indicates lowest error. Averages correspond to the whole dataset, and its ER column displays the amount of landscapes processed. The Valid row represents the averages over all valid landscapes. The first column links each row with the landscapes shown in Fig. 8. Columns [START_REF] Mats | Evaluation of image-shift measurement algorithms for solar shack-hartmann wavefront sensors[END_REF], [START_REF] Michau | Wavefront sensing from extended sources[END_REF], [START_REF] Poyneer | Scene-based shack-hartmann wavefront sensing: analysis and simulation[END_REF] and [START_REF] Sidick | Adaptive periodic-correlation algorithm for extended scene shack-hartmann wavefront sensing[END_REF] refer to SDF-2QI, Michau et al., Poyneer and APC methods respectively. The variants of our approach 1it, 2it L, 3it S represent Alg. 1 using the hypomode derivative with 1, 2 (bilinear interpolation) and 3 (spline interpolation) iterations respectively. The last column stands for one iteration and Gaussian derivative with σ=0.6

The PCC algorithm particularly failed under the landscape present in Fig. 8.7 although the reason is not evident. Finally, line 9 (Fig. 8.9) shows a case where the SDF-2QI method was outperformed by all variants of the proposed approach.

Simulated experiments

Experimental Setup. We used a simulator to evaluate the performance of the proposed method. Given a 12-bit input image, a set of images was generated matching a SHWFS configuration provided by CNES. A 12 × 12 lenslet grid occluded as shown in Fig. 2a was assumed. For each lenslet a ground truth shift was randomly generated. Since we were simulating in a closed loop environment, the maximum displacement was kept to half a pixel. Each sub-aperture image was then obtained by DFT resampling the input image, extracting a 37 × 37 pixels subimage, multiplying by the occlusion factor, and lastly adding noise.

The error for each method was computed as the mean error for all valid sub-apertures. A sub-aperture was considered invalid and therefore discarded if it was occluded on more than 60% of its surface. The error for each sub-aperture was the Euclidean distance between the measured and the ground truth shift. To evaluate the robustness to image content, we used three subimages: a highly contrasted one that should not present any difficulty for shift estimation, a slightly more challenging one with its gradient mainly distributed on a single direction, and one from the sea with almost no signal (there are just a few pixels of land on the bottom). Different amounts of additive white Gaussian noise were simulated (with standard deviation σ n ∈ [1, ..., 150]), and all methods were evaluated on 100 noise realizations. Fig. 9 shows the subimages affected by three noise levels. Signal dependent noise was handled by application of a variance stabilization transform. Finally, the processing time of each method was measured using a non-optimized Matlab implementation on an Intel Xeon E5-2650 CPU.

We tested several variants of the proposed GBSE method. Each one was composed by an interpolation method for resampling selected among the five methods proposed (see section 3.3), a gradient estimation method chosen between the hypomode or one of the three Gaussian derivatives with σ g = {0.3, 0.6, 1} (section 3.2), and using up to three iterations per shift estimation. Also, different variants using the multiscale approach were evaluated (section 3.5) up to 3 scales, where for each scale, an iterative GBSE algorithm was used with up to 3 iterations. From these methods we retained the best performing non-iterative method (k = 1 and no resampling) which uses σ g = 0.6 for the gradient computation, and the best one with 3 iterations which uses DFT resampling with image periodization and σ g = 0.3. We compared them with all state-of-the-art algorithms presented in section 2. For practical reasons, we only display the best performing ones, namely, the ACC [START_REF] Sidick | An Adaptive cross-correlation algorithm for extended scene Shack-Hartmann wavefront sensing[END_REF] and the APC [START_REF] Sidick | Adaptive periodic-correlation algorithm for extended scene shack-hartmann wavefront sensing[END_REF] methods from Sidick presented in section 2.3, and the Poyneer [START_REF] Poyneer | Scene-based shack-hartmann wavefront sensing: analysis and simulation[END_REF] and Löfdahl SDF-2QI algorithm [START_REF] Mats | Evaluation of image-shift measurement algorithms for solar shack-hartmann wavefront sensors[END_REF] presented in section 2.1. The SDF-2QI algorithm computes the squared sum difference function and searches its maximum by using 2D quadratic interpolation. All other presented methods were also evaluated but discarded due to lower accuracies. Results in a closed loop system. We compared iterative and non-iterative algorithms separately. The average errors and the standard deviations for non-iterating methods on the first image type are shown in Fig. 10 (top left), together with the processing time. The non-iterative GBSE version using Gaussian derivatives with σ = 0.6 outperforms both Poyneer and SDF-2QI (which perform similarly in this image). GBSE is also more stable (less variability) and faster. In Fig. 10 (top right) both iterative methods ACC and APC were compared with GBSE using 3 iterations, DFT resampling and the hypomode derivative. Again the proposed method proves to be the most accurate and more stable. Even more, this method is the best also compared to non-iterated methods. Notice the high impact of noise on the ACC method. It is more precise than both non-iterating methods for low (σ ≤ 20). Nevertheless, it diverges for stronger noise. This problem is well-known for phase-correlation methods [START_REF] Poyneer | Scene-based shack-hartmann wavefront sensing: analysis and simulation[END_REF].

The results for the second image type (bottom of Fig. 10) are similar to the first one, although the average errors are higher. While the SDF-2QI performs slightly worse, the Poyneer method is considerably worse. This is due to discrepancies at the image boundaries resulting from the periodicity assumption of the periodic convolution used by Poyneer. This also explains why both ACC and PCC methods behave so poorly on this image. This also implies that computing the shift in the spatial domain is usually more stable than doing so in the frequency domain, unless some windowing is performed, which is prohibitive on such small images where the objects are sometimes close to the image boundaries. Yet again, the proposed methods are more precise, stable and with less variability than the state-of-the-art. The third image type did not pass the verification step and all the methods failed, so its results are omitted. Results under larger displacements. The same experiment was performed by simulating wavefront aberrations having displacements of up to 4 pixels, to evaluate robustness against high frequency aberrations. In Fig. 11 error together with the standard deviations are shown for the most representative methods. Due to the larger underlying displacements, multiscale approaches achieve better results as the noise increases, reproducing the results of [START_REF] Martin Rais | Iterative gradient-based shift estimation: To multiscale or not to multiscale?[END_REF]. In fact, a multiscale approach becomes mandatory if fewer than two iterations are performed, and the best performing method used three scales (where a single iteration was used in the coarsest scale, two iterations with FFTP interpolation in the second scale and three iterations with FFTP interpolation in the final finer scale). Although this method proved to be the most accurate in our tests, it also requires more computational resources as evidenced by its running time. However, using the iterative non-multiscale approach with three iterations and FFT interpolation showed comparable results while requiring fewer processing time. Nevertheless, the method proposed by Poyneer should be considered when lower processing times are required. Indeed, as expected, non-iterated versions of the GBSE algorithm failed when estimating the shift with large underlying displacements. Let us add the caveat that when measuring the aberration using the second slightly challenging landscape (green rectangle in Fig 9), the Poyneer method obtained inaccurate results even under high SNR conditions.

Conclusion

A new method for accurate sub-pixel shift-estimation, based on an iterative global optical flow method, was proposed in the context of a SHWFS on extended scenes for on-board aberration correction on earth-observation satellites. Using a telescope simulator developed by CNES, we showed that the proposed algorithm is more accurate, stable, robust to noise, and with lower variability than the current state-of-the-art, permitting a more precise wavefront estimation. By performing on-board real-time mirror correction, this opens the way to cheaper high resolution earth-observation satellites.

It should be noted that although this method was conceived for on-board earth-observation satellites, it could also be used in a ground-based wavefront sensing context. Indeed, for ground-based applications, the demanding processing time requirements existing in satellites are considerably reduced, in which case, due to the availability of more computational resources, the proposed method could still be applicable using more iterations and/or scales. Therefore, it is envisaged to test the proposed algorithm for the case of atmospheric turbulence correction on extended sources.

Funding. Centre National d'Etudes Spatiales (MISS Project); Office of Naval research (grant N00014-14-1-0023); Ministerio de Economia y Competividad, Spain (BES-2012-057113); Agence Nationale de la Recherche (ANR-12-ASTR-0035).

  (a) Occlusion schema in a 12 × 12 SHWFS. (b) Images obtained on an extended landscape.

Figure 2 :

 2 Figure 2: Example of a SHWFS occlusion scheme under a Korsch telescope. Both the secondary mirror and the arms used to hold it are clearly visible. For each lenslet, the decrease of the incoming signal is proportional to its occlusion.

Figure 3 :

 3 Figure 3: Used pixels (gray background) for fast gradient estimation methods and their exact localization (red dots). Left: Centered differences. Center: 1D backward difference for both ∂x and ∂y. Right: 2D backward difference.

  072s) σ=0.6 (0.054s) σ=0.3 (0.043s) (e) Max. shift: 2 pixels.

  072s) σ=0.6 (0.054s) σ=0.3 (0.043s) (f) Max. shift: 2 pixels.

Figure 4 :

 4 Figure 4: Gradient estimation methods comparison for three different maximum shifts (0.1, 0.5 and 2 pixels). Left: Highly contrasted image (simple case). Right: Gradients mainly distributed on a single direction (challenging case). Noise std. dev. according to 12-bit images.

Figure 5 :

 5 Figure 5: Example of FFT resampling with and without image symmetrization. Resampling with DFT produces ringing due to the discontinuities at the periodized boundaries. No visible ringing is observed after resampling with symmetrization.

Figure 6 :

 6 Figure 6: Comparison of evaluated interpolation methods using 2 and 3 iterations and Gaussian derivatives with σ g = 0.6. Image noise was fixed to σ n = 80 and the maximum displacement was set to 2 pixels.

Figure 7 :

 7 Figure 7: CRLB and eigenratios respectively for six examples of landscapes. Green represents valid and red implies failure.

Figure 8 :

 8 Figure 8: Examples of landscapes from the CNES database.

( 1 )

 1 0.014 0.66 0.019 0.121 0.026 0.054 0.030 0.017 0.016 0.022 (2) 0.006 0.44 0.029 0.106 0.038 0.158 0.014 0.008 0.008 0.013 (3) 0.006 0.19 0.048 0.096 0.105 0.092 0.035 0.036 0.036 0.040 (4) 0.041 0.59 0.127 0.113 0.090 0.794 0.070 0.064 0.069 0.078 (5) 0.011 0.77 0.039 0.113 0.069 0.070 0.026 0.016 0.016 0.022 (6) 0.030 0.01 0.081 0.099 0.083 0.081 0.053 0.043 0.040 0.042 (7) 0.019 0.69 0.023 0.111 0.048 0.134 0.034 0.021 0.021 0.025 (8) 0.026 0.63 0.045 0.114 0.039 0.083 0.048 0.032 0.031 0.037 (9) 0.010 0.34 0.043 0.113 0.029 0.042 0.028 0.015 0.014 0.021 (10) 0.027 0.61 0.065 0.112 0.051 0.645 0.046 0.037 0.043 0.047 Avg. 0.029 132 0.146 0.111 0.491 0.478 0.047 0.039 0.039 0.043 Valid 0.011 64 0.024 0.110 0.048 0.057 0.027 0.018 0.017 0.022 Time 0.296 0.025 0.084 0.221 0.034 0.165 0.311 0.052

Figure 9 :

 9 Figure 9: Top: Input image used for the simulations. Bottom: For each level of noise σ n = {1, 100, 200} (vertically separated), two different sub-apertures are shown: no occlusion and 57% occluded (horizontally separated). The dynamic ranges were stretched for viewing purposes.

Figure 10 :

 10 Figure 10: Performance for non-iterative and iterative methods compared to the best variant of the proposed GBSE approach. Top: First test image. Bottom: Second test image. Left: Non-iterative methods. Right: Iterative methods.

Figure 11 :

 11 Figure 11: Average error and standard deviations of selected methods and two variants of the proposed GBSE approach when simulating wavefront aberrations yielding displacements of up to 4 pixels. Left: First test image. Right: Second test image.

Table 1 :

 1 CRLB, EigenRatio (ER) and errors (in pixels) for selected sample landscapes of the CNES database. a

	CRLB ER [18] [21]	[26]	[39]

  the average

		0.1						0.15
		0.07 0.08 0.09			SDF-2QI (0.303s) POYNEER (0.0823s) APC0.001 (0.351s) 3 it. σ=0.6 (F) (0.277s) MS3(1,2(FS),3(FS)).HM (0.638s)	0.135 0.12 0.105			SDF-2QI (0.29s) POYNEER (0.0806s) APC0.001 (0.442s) 3 it. σ=0.6 (F) (0.27s) MS3(1,2(FS),3(FS)).HM (0.621s)
	Error (pixels)	0.04 0.05 0.06					Error (pixels)	0.06 0.075 0.09
		0.03						0.045
		0.02						0.03
		0.01						0.015
		0	0	20	40	60	80 100 120 140 160 180 200	0	0	20	40	60	80 100 120 140 160 180 200
							Noise std. dev.						Noise std. dev.


Simon Hampe
email: simon.hampe@googlemail.com

Michael Joswig
email: joswig@math.tu-berlin.de

Benjamin Schröter
email: schroeter@math.tu-berlin.de

Algorithms for Tight Spans and Tropical Linear Spaces

Keywords: Tropical linear spaces, Tight spans, Closure Systems, polymake 2010 MSC: 05B35, 05-04, 52B40, 14T05

We describe a new method for computing tropical linear spaces and more general duals of polyhedral subdivisions. This is based on an algorithm of Ganter for finite closure systems.

Introduction

Tropical linear spaces are among the most basic objects in tropical geometry [START_REF] Maclagan | Introduction to tropical geometry[END_REF]Chapter 4]. In combinatorial terms they form polyhedral complexes which are dual to regular matroid subdivisions of hypersimplices. Such subdivisions are characterized by the property that each cell is the convex hull of characteristic vectors of the bases of a matroid. Here the hypersimplices correspond to the uniform matroids. Research on matroid subdivisions and related objects goes back to Dress and Wenzel [START_REF] Dress | Valuated matroids[END_REF] and to Kapranov [START_REF] Mikhail | Chow quotients of Grassmannians. I. In I[END_REF]. Speyer instigated a systematic study in the context of tropical geometry [START_REF] Speyer | Tropical linear spaces[END_REF], while suitable algorithms have been developed and implemented by Rincón [START_REF] Rincón | Computing tropical linear spaces[END_REF].

Here we present a new combinatorial algorithm for computing tropical linear spaces, which are not necessarily realizable. This implemented in the software system polymake [START_REF] Gawrilow | polymake: a framework for analyzing convex polytopes[END_REF]. Moreover, we report on computational experiments. Our approach has two key ingredients. First, our method is completely polyhedral -in contrast with Rincón's algorithm [START_REF] Rincón | Computing tropical linear spaces[END_REF], which primarily rests on exploiting matroid data. Employing the polyhedral structure has the advantage that this procedure naturally lends itself to interesting generalizations and variations. In particular, this includes tropical linear spaces corresponding to non-trivially valuated matroids. Second, our method fundamentally relies on an algorithm of Ganter [START_REF] Ganter | Algorithmen zur formalen Begriffsanalyse[END_REF] for enumerating the closed sets in a finite closure system; cf. [GO16, §2.2] This procedure is a variant of breadth-first-search in the Hasse diagram of the poset of closed sets. As a consequence the computational costs grow only linearly with the number of edges in the Hasse diagram, i.e., the number of covering pairs among the closed sets. So this complexity is asymptotically optimal in the size of the output, and this is what makes our algorithm highly competitive in practice. The challenge is to implement the closure operator and to intertwine it with the search in such a way that it does not impair the output-sensitivity.

Kaibel and Pfetsch employed Ganter's algorithm for enumerating face lattices of convex polytopes [START_REF] Kaibel | Computing the face lattice of a polytope from its vertex-facet incidences[END_REF], and this was later extended to bounded subcomplexes of unbounded polyhedra [START_REF] Herrmann | Computing the bounded subcomplex of an unbounded polyhedron[END_REF]. Here this is generalized further to arbitrary regular subdivisions and their duals. Such a dual has been called tight span in [START_REF] Herrmann | Dressians, tropical Grassmannians and their rays[END_REF] as it generalizes the tights spans of finite metric spaces studied by Isbell [START_REF] Isbell | Six theorems about injective metric spaces[END_REF] and Dress [START_REF] Dress | Trees, tight extensions of metric spaces, and the cohomological dimension of certain groups: a note on combinatorial properties of metric spaces[END_REF]. The tight span of an arbitrary polytopal complex may be seen as a special case of the dual block complex of a cell complex; e.g., see [START_REF] Munkres | Elements of algebraic topology[END_REF]§64]. From a topological point of view subdivisions of point configurations are cell decompositions of balls, which, in turn, are special cases of manifolds with boundary. The duality of manifolds with boundary is classically known as Lefschetz duality (e.g., see [START_REF] Munkres | Elements of algebraic topology[END_REF]§70]), and this generalizes Poincaré duality as well as cone polarity. With an arbitrary polytopal subdivision, Σ, we associate a new object, called the extended tight span of Σ, which contains the tight span, but which additionally takes duals of certain boundary cells into account. In general, the extended tight span is only a partially ordered set. If, however, Σ is regular, then the extended tight span can be equipped with a natural polyhedral structure. We give an explicit coordinatization. In this way tropical linear spaces arise as the extended tight spans of matroid subdivisions with respect to those boundary cells which correspond to loop-free matroids. While a tropical linear space can be given several polyhedral structures, the structure as an extended tight span is the coarsest. Algorithmically, this has the advantage of being the sparsest, i.e., being the one which takes the least amount of memory. In this sense, this is the canonical polyhedral structure of a tropical linear space.

This paper is organized as follows. We start out with recalling basic facts about general closure systems with a special focus on Ganter's algorithm [START_REF] Ganter | Algorithmen zur formalen Begriffsanalyse[END_REF]. Next we introduce the extended tight spans, and this is subsequently specialized to tropical linear spaces. We compare the performance of Rincón's algorithm [START_REF] Rincón | Computing tropical linear spaces[END_REF] with our new method; see also Section 4.2 for further comments. To exhibit one application the paper closes with a case study on the f -vectors of tropical linear spaces.

Closure Systems, Lower Sets and Matroids

While we are mainly interested in applications to tropical geometry, it turns out that it is useful to start out with some fundamental combinatorics. This is the natural language for Ganter's procedure, which we list as Algorithm 1 below.

Definition 2.1. A closure operator on a set S is a function cl : P(S) → P(S) on the power set of S, which fulfills the following axioms for all subsets A, B ⊆ S:

(i) A ⊆ cl(A) (Extensiveness). (ii) If A ⊆ B then cl(A) ⊆ cl(B) (Monotonicity). (iii) cl(cl(A)) = cl(A) (Idempotency). A subset A of S is called closed, if cl(A) = A.
The set of all closed sets of S with respect to some closure operator is called a closure system.

The closed sets of a closure system form a meet-semilattice. In other words, the set S is closed and the intersection A ∩ B of closed sets A, B is closed. Conversely, each meet-semilattice arises in this way. See [START_REF] Caspard | The lattices of closure systems, closure operators, and implicational systems on a finite set: a survey[END_REF] for a survey on closure systems and lattices.

Classical examples include the following. If the set S carries a topology then the function which sends any subset A to the smallest closed set (defined as the complement of an open set) containing A is a closure operator, called the topological closure. If the set S is equipped with a group structure then the function which sends any subset A to the smallest subgroup containing A is a closure operator. Throughout the following we are particularly interested in the case where the set S = [n] := {1, . . . , n} is finite.

The closed sets of a closure system (S, cl) are partially ordered by inclusion. The resulting poset is the closure poset induced by (S, cl). The Hasse diagram of (S, cl) is the directed graph whose nodes are the closed sets and whose arcs correspond to the covering relations of the closure poset. We assume that all arcs are directed upward, i.e., toward the larger set. Ganter's algorithm (cf. Algorithm 1 below) computes the Hasse diagram of a finite closure system; cf. [START_REF] Ganter | Algorithmen zur formalen Begriffsanalyse[END_REF]. As its key property each closed set is pushed to the queue precisely once, and this entails that the running time is linear in the number of edges of the Hasse diagram, i.e., the algorithm is output-sensitive. [START_REF] Kaibel | Computing the face lattice of a polytope from its vertex-facet incidences[END_REF] proposed a method to compute the face lattice of a convex polytope P . This can be done in two different ways: A face of a polytope can either be identified by its set of vertices or by the set of facets it is contained in.

In the first case, the set S is the set of vertices and the closure of a set is the smallest face containing this set.

In the second case, the set S is the set of facets. Let F ⊆ S. The intersection of the facets in F is a face Q F of P . The closure of F is defined as the set of all facets which contain Q F . Note that, with this approach, Algorithm 1 actually computes the face lattice with inverted relations.

In both cases, the closed sets are exactly the faces of P and the closure operator is given in terms of the vertex-facet incidences.

Example 2.4. Instead of polytopes, one can also compute the face lattice of a polyhedral fan in much the same manner. The crucial problem is to define the closure of a set of rays which is not contained in any cone. The solution to this is to extend the set S to contain not only all rays, but also an additional artificial element, say ∞. Now the closure of a set F ⊆ S is either the smallest cone containing it, if it exists, or the full set S. In particular, this ensures that the length of a maximal chain in the face lattice of a kdimensional fan is always k + 1.

The following class of closure systems is ubiquitous in combinatorics and tropical geometry. The monographs by White [START_REF]Theory of matroids, volume 26 of Encyclopedia of Mathematics and its Applications[END_REF] and Oxley [START_REF] Oxley | Matroid theory[END_REF] provide introductions to the subject. Definition 2.5. Let S be a finite set equipped with a closure operator cl : P(S) → P(S). The pair (S, cl) is a matroid if the following holds in addition to the closure axioms:

(iv) If A ⊆ S and x ∈ S, and y ∈ cl(A ∪ {x}) \ cl(A), then x ∈ cl(A ∪ {y}) (MacLane-Steinitz Exchange).
This is one of many ways to define a matroid; cf. [Oxl11, Lemma 1.4.3] for explicit cryptomorphisms. The closed sets of a matroid are called flats.

Remark 2.6. For matroids it is not necessary to check for the minimality of the closed sets N i in Algorithm 1. In view of Axiom (iv) this is always satisfied. This application of the algorithm also demonstrates that, while the empty set is typically closed, this does not always need to be the case. In fact, the closure of the empty set in a matroid is the set of all loops.

For our applications it will be relevant to consider special closure systems which are derived from other closure systems in the following way. A lower set Λ of the closure system (S, cl) is a subset of the closed sets such that for all pairs of closed sets with A ⊆ B we have that B ∈ Λ implies A ∈ Λ.

Proposition 2.7. Let (S, cl) be a closure system with lower set Λ. Then the function cl Λ which is defined by

cl Λ (A) = cl(A) if cl(A) ∈ Λ , S otherwise (1)
is a closure operator on S.

Proof. Extensiveness and idempotency are obvious. We need to show that monotonicity holds. To this end consider two closed sets A ⊆ B ⊆ S. Suppose first that A lies in the lower set Λ. Then cl

Λ (A) = cl(A) ⊆ cl(B) ⊆ cl Λ (B). If, however, A ∈ Λ, then B ∈ Λ as Λ is a lower set. In this case we have cl Λ (A) = S = cl Λ (B).
Example 2.8. An unbounded convex polyhedron is pointed if it does not contain any affine line. In that case the polyhedron is projectively equivalent to a convex polytope, with a marked face, the face at infinity; see, e.g., [START_REF] Joswig | Polyhedral and algebraic methods in computational geometry[END_REF]Theorem 3.36]. So we arrive at the situation where we have a convex polytope P with a marked face F . Now the set of faces of P which intersect F trivially forms a lower set Λ in the closure system of faces of P . In this way combining Example 2.3 with Proposition 2.7 and Algorithm 1 gives a method to enumerate the bounded faces of an unbounded polyhedron. Ignoring the entire set S, which is closed with respect to cl Λ but not bounded, recovers the main result from [START_REF] Herrmann | Computing the bounded subcomplex of an unbounded polyhedron[END_REF].

Example 2.9. For a d-polytope P and k ≤ d the faces of dimension at most k form a lower set. This is the k-skeleton of P .

The closure operators from all examples in this section are implemented in polymake [START_REF] Gawrilow | polymake: a framework for analyzing convex polytopes[END_REF].

Extended Tight Spans

It is the goal of this section to describe duals of polytopal complexes in terms of closure systems. Via Algorithm 1 this gives means to deal with them effectively. For details on polyhedral subdivisions we refer to the monograph [START_REF] Jesús | Triangulations[END_REF].

Let P ⊂ R d be a finite point configuration, and let Σ be a polytopal subdivision of P. That is, Σ is a polytopal complex whose vertices lie in P and which covers the convex hull P = conv P. We call the elements of Σ cells; the set of maximal cells is denoted by Σ max and the maximal boundary facets (meaning those cells of Σ contained in the facets of P which are maximal with respect to inclusion) by ∆ Σ . Now we obtain a closure operator on the set

γ 1 γ 2 σ 1 σ 2 ∅ {σ 1 } {σ 2 } {σ 1 , γ 1 } {σ 1 , σ 2 } {σ 2 , γ 2 } {σ 1 , σ 2 , γ 1 , γ 2 }
S Σ := Σ max ∪ ∆ Σ by letting cl Σ (F ) := ∅ if F = ∅ , g ∈ S Σ σ∈F σ ⊆ g otherwise . (2)
for any F ⊆ S Σ . Note that cl Σ is basically the same as the dual operator in Example 2.3. In fact, the closed, non-empty sets in S Σ correspond to the cells of Σ, while the poset relation is the inverse containment relation. Now let Γ be a collection of boundary faces of Σ. This defines a lower set Λ Γ for the closure system (S Σ , cl Σ ), which consists of all sets F , such that σ∈F σ τ for any τ ∈ Γ. We will denote the corresponding closure operator by cl Σ Γ := cl Σ Λ Γ and we call the resulting closure system (S Σ , cl Σ Γ ) the extended tight span of Σ with respect to Γ.

If Γ = ∆ Σ , the closed sets are all cells of Σ which are not contained in the boundary. This is exactly the tight span of a polytopal subdivision defined in [START_REF] Herrmann | How to draw tropical planes[END_REF], which is dual to the interior cells. Note that this can also be obtained as the closure system of (Σ max , cl Σ ).

Example 3.1. Let P be {-1, 0, 1} with the convex hull [-1, 1] and Σ be its subdivision into intervals σ 1 = [-1, 0] and σ 2 = [0, 1], i.e., Σ max = {σ 1 , σ 2 }. The subdivision and its corresponding extended tight spans for Γ = ∅ and Γ = ∆(Σ) = {{-1}, {1}} = {γ 1 , γ 2 } can be seen in Figure 1. This example also demonstrates that we need to declare the closure of the empty set to be itself to ensure monotonicity.

If the subdivision is regular, i.e., induced by a height function h : P → R, we can actually coordinatize the extended tight span. Any regular subdivision with fixed height function is dual to a dual complex N Σ , which is a complete polyhedral complex in R d . More precisely, for every point x ∈ R d there is a cell of Σ, consisting of all points p ∈ P which minimize h(p) -p • x. Points which induce the same cell form an open polyhedral cell and the topological closures of these cells form the dual complex. It is a well-known fact that there is a bijective, inclusion-reversing relation between the cells of Σ and the cells of N Σ .

In particular, each maximal cell of Σ max is dual to a vertex and every boundary facet in ∆(Σ) is dual to a ray of N Σ . Hence, every closed set F of (S Σ , cl Σ Γ ) corresponds to a polyhedral cell ρ F and together these cells form a subcomplex of N Σ . More precisely, we denote by

T Σ,Γ := ρ F F S Σ closed w.r.t. cl Σ Γ (3)
the coordinatized extended tight span of Σ with respect to Γ. Its face lattice is by definition the poset of closed sets of (S Σ , cl Σ Γ ).

Tropical Linear Spaces

In this section we will finally investigate the objects that we are most interested in: valuated matroids and tropical linear spaces. We prefer to see the latter as special cases of extended tight spans. Valuated matroids were first studied by Dress and Wenzel [START_REF] Dress | Valuated matroids[END_REF]; cf. [MS15, Chapter 4] for their role in tropical geometry.

Let us introduce some notation. For a subset B of [n] of size r, let e B := i∈B e i . For a collection M ⊆ [n] r of such subsets we let

P M := conv{e B | B ∈ M } (4)
be the subpolytope of the hypersimplex ∆(r, n) which is spanned by those vertices which correspond to bases of M . The following fundamental characterization is known [START_REF] Edmonds | Submodular functions, matroids, and certain polyhedra[END_REF][START_REF] Izrail | Combinatorial geometries, convex polyhedra, and Schubert cells[END_REF].

Proposition 4.1. The set M comprises the bases of a matroid if and only if the vertex-edge graph of the polytope P M is a subgraph of the vertex-edge graph of ∆(r, n) or, equivalently, if every edge of P M is parallel to e i -e j for some i and j.

Throughout the following, let M be (the set of bases of) a matroid on n elements. In that case P M is the matroid polytope of M . The matroid M is said to be loop-free if B∈M B = [n]. The rank of M is r, the size of any basis. If P M is the full hypersimplex, then M = U r,n is a uniform matroid. The above description fits well with our geometric approach. Any function v : M → R gives rise to a regular subdivision on P M , which we denote by Σ M,v . The pair (M, v) is a valuated matroid if every cell of Σ M,v is again a matroid polytope. Then Σ M,v is called a matroid subdivision.

Example 4.2. The set M of subsets of {1, 2, 3, 4} with exactly two elements has cardinality six. Their characteristic vectors are the vertices of a regular octahedron embedded in 4-space. If we let v be the map which sends five vertices to 0 and the sixth one to 1, then (M, v) is a valuated matroid.

We will define tropical linear spaces as duals of valuated matroids. To this end let (M, v) be a valuated matroid of rank r on n elements. For a vector x ∈ R n , we define the set

M x := {B ∈ M | v(B) -e B • x is minimal} .
(5)

From the definition of the dual complex in Section 3 we see that the elements of M x correspond to a cell of Σ M,v and thus define a matroid. Note that for any λ ∈ R we clearly have M x = M x+λ1 .

Definition 4.3. The tropical linear space associated with the valuated matroid (M, v) is the set

B(M, v) := {x ∈ R n | M x is loop-free} /R1 . (6)
Note that our definition of a valuated matroid, as well as that of a tropical linear space are with respect to minimum as tropical addition. Tropical linear spaces generalize tropicalizations of linear varieties, i.e., the image of a linear space over a valuated field under its valuation map; cf. [MS15, §4.4].

The following is our main result. While it is easy to prove, it is relevant since it entails a new effective procedure for enumerating the cells of a tropical linear space via Algorithm 1.

Theorem 4.4. Let Γ be the set of boundary faces of Σ := Σ M,v which correspond to matroids with loops. Then

B(M, v) = T Σ,Γ /R1 , (7)
where T Σ,Γ is the coordinatized extended tight span defined in Section 3.

Proof. Let N Σ denote the dual complex of Σ. From our definition it is immediately clear that B(M, v) is a subcomplex of N Σ /R1. A cell in N Σ /R1 lies in B(M, v) if its dual cell in Σ is the matroid polytope of a loop-free matroid. Since any cell in Σ corresponds to a loop-free matroid, if and only if it is not contained in a boundary facet of a matroid with loops, the claim follows.

We call the resulting polyhedral structure of B(M, v) canonical.

Remark 4.5. Note that one can naturally replace Γ by the subset of maximal boundary faces corresponding to matroids with loops. These faces are defined by the equations

x i = 0 for i ∈ [n].
Example 4.6. If the valuation is constant then the matroid subdivision is trivial. It follows that the dual complex coincides with the normal fan of the matroid polytope P M . In this case B(M, v) is the Bergman fan of M , in its coarsest possible subdivision; cf. [START_REF] Hampe | Algorithmic aspects of tropical intersection theory[END_REF] for a proof.

The polyhedral complex B(M, v) reflects quite a lot of the combinatorics of the matroid M . For instance, the rank of M equals dim(B(M, v))+1. If L is the lineality space of B(M, v), then the number of connected components of M is dim(L) + 1; cf. [START_REF] Edmonds | Submodular functions, matroids, and certain polyhedra[END_REF][START_REF] Eva | Matroid polytopes, nested sets and bergman fans[END_REF].

Performance comparison

As mentioned in the introduction, there is an algorithm by Rincón [START_REF] Rincón | Computing tropical linear spaces[END_REF] for computing Bergman fans, i.e., tropical linear spaces with trivial valuation. An extension which can also deal with trivially valuated arbitrary matroids which may not be realizable has been implemented in polymake's bundled extension a-tint [START_REF] Hampe | a-tint: a polymake extension for algorithmic tropical intersection theory[END_REF]. It is this implementation we refer to in the following discussion. The original software TropLi by Rincón only takes realizable matroids as input. Rincón's and our algorithm are very difficult to compare for two reasons. First of all, a matroid can be encoded in numerous ways. For instance, in terms of closures, as in Definition 2.5, or in terms of bases, as in Proposition 4.1. Many further variants exist, and the conversion between these representations is often a non-trivial computational task. Below we will assume that all matroids are given in terms of their bases. The second problem is that the two algorithms essentially compute very different things. Our algorithm computes the full face lattice of the canonical polyhedral structure of a tropical linear space. On the other hand, Rincón's algorithm only computes the rays and the maximal cones of the Bergman fan, albeit in a finer subdivision. In this setup it is therefore to be expected that our approach is significantly slower. In particular, to identify the boundary cells (including the loopfree ones), we need to apply a convex hull algorithm to the matroid polytope before we can make use of our algorithm. Still, the discussion has merit when separating the timings for the different steps; cf. Table 1. We compute Bergman fans of all (isomorphism classes of) matroids of a given rank r on a given ground set [n] as provided at http://www-imai.is. s.u-tokyo.ac.jp/ ˜ymatsu/matroid); see also [START_REF] Matsumoto | Matroid enumeration for incidence geometry[END_REF]. Each matroid is given only in terms of its bases. We first apply Rincón's algorithm and then compute the Hasse diagram of the face lattice of the fan as described in Example 2.4. For our approach we split the computations into two steps: First we compute the convex hull of the matroid polytope, displayed under "CH" and then measure the running time of our closure algorithm "ETS" (extended tight span) separately. Times were measured on an AMD Phenom II X6 1090T with 3.6 GHz using a single thread and polymake version 3.1. We employed the double description method implemented in the Parma Polyhedral Library (via polymake's interface) for computing the convex hulls [START_REF] Bagnara | The Parma Polyhedra Library: toward a complete set of numerical abstractions for the analysis and verification of hardware and software systems[END_REF].

The results show that almost all of the time in our algorithm is spent computing the facets of the matroid polytope. On the other hand, if one aims at obtaining the same amount of information, i.e., the full face lattice, for Rincón's algorithm, this increases the computation time dramatically. This demonstrates that the finer subdivision produced by this algorithm is significantly worse in terms of complexity than the canonical subdivision.

We also like to point out that for non-trivial valuations our algorithm is, to the best of our knowledge, currently the only feasible method for computing tropical linear spaces.

Other ways to compute tropical linear spaces

For the sake of completeness we briefly comment on alternatives to computing with tropical linear spaces.

Let (M, v) be a valuated matroid of rank r on n elements. The associated tropical linear space B(M, v) consists of all points x in the (n-1)dimensional tropical projective space such that, for any (d+1)-subset τ of [n], the minimum of the numbers v(τ \ {τ i }) + x τ i , for i ∈ [d] is attained at least twice. The latter describes a circuit of (B, v). Dually, the cocircuits of that valuated matroid arise as follows. For any (d-1)-subset σ of [n] let v(σ * ) be the vector in (R ∪ {∞}) n whose jth coordinate is v(σ ∪ {j}) and ∞ if j ∈ σ. It now follows from a result of Yu and Yuster [YY07, Theorem 16] that B(M, v) is the tropical convex hull of all its cocircuits; see also [START_REF] Joswig | Affine buildings and tropical convexity[END_REF] for more details and an application.

This raises the question of how to compute the tropical convex hull of the cocircuits. First of all, there is more than one possible output format; cf. [START_REF] Joswig | Tropical convex hull computations[END_REF] for a survey. None of the known methods is very efficient; e.g., see [JSY07, Table 1] for some timings. Even worse, for most families of matroids, the number of cocircuits grows very quickly with r and n. That is to say, a tropical convex hull computation on such input is not computationally feasible under most circumstances.

Dually, the tropical linear space B(M, v) can be written as a tropical prevariety of linear forms, i.e., as the intersection of tropical hyperplanes coming from the circuits; cf. [MS15, (4.4.4)]. Assuming that there are about as many circuits as bases it is possible to apply a dual version of Algorithm 1 in this setting. Yet, compared with our approach, this is somewhat more involved since each tropical hyperplane needs to be represented as an ordinary polyhedral complex with Θ(n 2 ) many maximal cells.

A case study on f -vectors of tropical linear spaces

Throughout the rest of this paper we will restrict ourselves to valuations of uniform matroids. Equivalently, we study matroid subdivisions of hypersimplices (and their lifting functions). Speyer was the first to conduct a thorough study of the combinatorics of tropical linear spaces [START_REF] Speyer | Tropical linear spaces[END_REF]. He conjectured the following.

Conjecture 5.1 (Speyer's f -vector conjecture). Let n ≥ 1 and 1 ≤ r ≤ n. Let v be any valuation on U r,n . Then the number of (i -1)-dimensional bounded faces of

B(U r,n , v) is at most n-2i r-i n-i-1 i-1 .
To study this problem, one is naturally interested in some form of moduli space of all possible valuations on U := U r,n . This role is played by the Dressian Dr(r, n); see [START_REF] Herrmann | How to draw tropical planes[END_REF][START_REF] Herrmann | Dressians, tropical Grassmannians and their rays[END_REF]. It is a subfan of the secondary fan of P U = ∆(r, n), consisting of all cones which correspond to matroidal subdivisions. As a set it contains the tropical Grassmannian TGr p (r, n) for any characteristic p ≥ 0. This is the tropicalization of the ordinary Grassmannian over an algebraically closed field of characteristic p, and it consists of all cones of the secondary fan which correspond to realizable valuations on U , i.e., those which can be realized as valuated vector matroids in characteristic p; see [START_REF] Speyer | The tropical Grassmannian[END_REF] and [START_REF] Maclagan | Introduction to tropical geometry[END_REF]Chapter 4]. However, this inclusion is generally strict. In fact, the Dressian is not even pure in general.

Remark 5.2. For r = 2, the Dressian Dr(r, n) is equal to the tropical Grassmannian. Combinatorially, this is the space of phylogenetic trees; e.g., see [START_REF] Mikhail | Chow quotients of Grassmannians. I. In I[END_REF]§1.3] and [MS15, §4.3]. For r = 3 and 3 ≤ n ≤ 6, the equality Dr(r, n) = TGr p (r, n) still holds on the level of sets for each p ≥ 0. This is trivial for n = 3, 4, as there are no non-trivial subdivisions of P U in that case. For n = 5 it follows from duality and the statement for Dr(2, 5). The Dressian Dr(3, 6) was computed in [START_REF] Speyer | The tropical Grassmannian[END_REF]. Note that, while the Dressian and the Grassmannian may agree as sets, they can have different polyhedral structures. Understanding the precise relation between these fan structures is still an open problem for general parameters. The cases (3, 7) and (3, 8) are the first where the Dressian differs from the Grassmannian. The Dressian Dr(3, 7) was computed in [START_REF] Herrmann | How to draw tropical planes[END_REF]. In particular, the possible combinatorial types of the corresponding tropical planes (and thus, their possible f -vectors) were listed. The polyhedral fan Dr(3, 8) was computed in [START_REF] Herrmann | Dressians, tropical Grassmannians and their rays[END_REF].

The Dressian Dr(3, 8)

We wish to compute f -vectors of uniform tropical planes in R 8 /R1, i.e., tropical linear spaces corresponding to valuations on U 3,8 . To this end, we make use of the data obtained in [START_REF] Herrmann | Dressians, tropical Grassmannians and their rays[END_REF], which is available at http:// svenherrmann.net/DR38/dr38.html. There is a natural S 8 -symmetry on the Dressian and the web page provides representatives for each cone orbit.

We computed tropical linear spaces for each cone by choosing an interior point as valuation. For the sake of legibility, we only include results for the maximal cones of the Dressian. There are 14 maximal cones of dimension nine and 4734 maximal cones of dimension eight. The full data can be obtained at http://www3.math.tu-berlin.de/combi/dmg/ data/f_vectors_Dr38.

Convention. The f -vector of a tropical linear space L is the f -vector of its canonical polyhedral structure. The bounded f -vector of L is the f -vector of the bounded part of this structure. All counts are given modulo the S 8 -symmetry on the Dressian. There is only one bounded f -vector (n -2, n -3) for a tropical linear space that corresponds to a maximal cone in the Dressian Dr(2, n), since this linear space has the combinatorics of a binary tree with n labeled leaves. The generic tropical linear spaces in the Dressian Dr(3, 6) have a bounded f -vector which is either (5, 4, 0) or (6, 6, 1); see [START_REF] Herrmann | How to draw tropical planes[END_REF]. In the case of (3, 7) the (generic) bounded f -vectors read (7, 6, 0), (9, 10, 2) and (10, 12, 3).

Theorem 5.3. Every generic tropical plane in R 8 /R1 has one of four possible f -vectors:

If it corresponds to a nine-dimensional cone in the Dressian, its fvector is (13, 55, 63) and its bounded f -vector is (13, 15, 3). There are nine different combinatorial types of such planes; see Figure 2.

If it corresponds to one of the 4734 eight-dimensional maximal cones in the Dressian, there are three possibilities:

• There are 51 planes with f -vector (13, 56, 64) and bounded fvector (13, 16, 4). • There are 1079 planes with f -vector (14, 58, 65) and bounded fvector (14, 18, 5). • There are 3604 planes with f -vector (15, 60, 66) and bounded fvector (15, 20, 6).

There are 3013 different combinatorial types of such planes.

The maximal bounded f -vector (15, 20, 6) agrees with the upper bound in Conjecture 5.1.

Remark 5.4. Each of the nine different combinatorial types that correspond to a nine-dimensional cone contain a vertex (marked in white in Figure 2), which in turn corresponds to the matroid polytope of a parallel extension of the Fano matroid. This is a certificate that these tropical linear spaces are not realizable over any field of characteristic greater than two; see [Oxl11, Chapter 6 and Appendix]. Figure 3 illustrates the connected extensions of the Fano matroid; these are those that are loop-free.

Further computer experiments reveal the following details.

Proposition 5.5. Let p be 0, 3, 5 or 7. Then the intersection of the relative interior of a cone C in the Dressian Dr(3, 8) with the tropical Grassmannian TGr p (3, 8) is trivial if and only if a subdivision which is induced by a lifting in the relative interior of C contains the polytope of a Fano matroid extension as a cell.

Outlook

Higher Dressians

We have given an algorithm which computes tropical linear spaces for arbitrary valuations in reasonable time; computing all tropical linear spaces for Dr(3, 8) above only took a few hours on a standard desktop computer. This indicates that it is feasible to apply this algorithm more ambitiously, e.g., to Dressians with larger parameters. However, in these cases not much data is currently available. Computing higher Dressians is a challenging task in itself.

The next step would be to look at Dr(4, 8). While computing the full Dressian is, at the moment, beyond our means, we can consider the following construction by Speyer [START_REF] Speyer | Tropical linear spaces[END_REF]. Let M be a matroid of rank r on n elements. We define an associated valuation on U r,n by

v M (B) := r -rank M (B) , (8)
where B ∈ [n] r is a basis of the uniform matroid and rank M (B) = max

B ∈M {|B ∩ B |} (9)
is the rank of B in M . Speyer showed that the corank indeed defines a valuation and that the matroid polytope P M appears as a cell in the induced regular subdivision. There are 940 isomorphism classes of matroids of rank four on eight elements [START_REF] Matsumoto | Matroid enumeration for incidence geometry[END_REF]; our computation is based on the data from http:// www-imai.is.s.u-tokyo.ac.jp/ ˜ymatsu/matroid. For computing the tropical linear spaces given by the valuations defined above we employed the enriched version available at db.polymake.org. This is certainly not enough to provide a global view on the Dressian Dr(4, 8), but it gives us a first glimpse of relevant combinatorial features. There are 62 different bounded f -vectors of such tropical linear spaces, so we cannot list them all. Also, up to combinatorial isomorphism, there are 465 different subdivisions of the hypersimplex induced by these matroids. As an example, consider the matroid M := U ⊕4 1,2 ; see [Oxl11, Chapter 4.2] for more on direct sums of matroids. The bounded f -vector of the tropical linear space induced by v M is (14, 24, 12, 1). In particular, the last two entries already achieve the respective maxima conjectured by Speyer, which are (20, 30, 12, 1). Experiments suggests that this is generally true, i.e., if M = U ⊕d 1,2 , then the valuation on U d,2d gives a linear space whose bounded f -vector maximizes the last two entries. Among valuations of the form v M on U 4,8 , the maximal number of edges is in fact also 24. However, the maximal number of vertices is 15. This is achieved by the unique matroid with 56 bases and 14 hyperplanes, i.e., flats of rank three. For experts: This is a sparse paving matroid, which has the maximal number 16 of cyclic flats among all matroid of rank four on eight elements.

Further optimization

Many of the objects considered here, such as polytopes, fans and matroids, exhibit symmetries which are also visible in the corresponding closure systems. It seems desirable, therefore, to exploit this during the computation. For every orbit of a closed set, only one representative would be computed. In a first approach, this could be achieved by considering equivalent sets to be the same in Algorithm 1: Once, when collecting all minimal closures cl(N ∪ {i}) and again when checking if N i is already in the graph. One could then easily recover the list of all closed sets in the end, though reconstructing the full poset structure (i.e. without symmetry) would require significant computational work.

As mentioned in 4.1, the most expensive part in our computations is a convex hull algorithm for computing the subdivision and the facets of the matroid polytope. It is known that the facets can be described in terms of the combinatorics of the matroid; cf. [START_REF] Edmonds | Submodular functions, matroids, and certain polyhedra[END_REF][START_REF] Eva | Matroid polytopes, nested sets and bergman fans[END_REF]. It is unclear if such a description can be given for the regular subdivision.

Figure 1 :

 1 Figure1: A regular subdivision and its extended tight span for Γ = ∅ and Γ = ∆(Σ) = {γ1, γ2}, respectively. The latter is marked in gray.

Figure 2 :

 2 Figure2: The various combinatorial types of bounded parts of tropical linear spaces corresponding to nine-dimensional cones in the Dressian. Note that all of them share the same f -vector (13, 15, 3). The naming convention is P = pentagon, S = square, T = triangle. The star * indicates where the square has an additional edge attached.

Figure 3 :

 3 Figure 3: The three loop-free extensions of the Fano matroid.

Table 1 :

 1 Comparing running times for computing Bergman fans.

	(n, r)	# matroids	Rincón	Hasse	CH	ETS
	(6,2)	23	0.0	0.2	0.8	0.0
	(6,3)	38	0.0	0.4	1.6	0.0
	(7,2)	37	0.0	0.3	1.6	0.0
	(7,3)	108	0.0	1.5	5.8	0.2
	(8,2)	58	0.0	0.4	1.9	0.0
	(8,3)	325	0.3	6.0	21.4	0.8
	(8,4)	940	1.8	48.7	86.5	9.2

(Benjamin Schröter) 1 S. Hampe received support from Deutsche Forschungsgemeinschaft (DFG) (Priority Program 1489, Grant JO 366/3-2).

2 Research by M. Joswig and B. Schröter is carried out in the framework of Matheon supported by Einstein Foundation Berlin. Further support by DFG (SFB-TRR 109: "Discretization in Geometry and Dynamics" and SFB-TRR 195: "Symbolic Tools in Mathematics and their Application"). The authors would like to thank the Institute Mittag-Leffler for its hospitality during the program "Tropical Geometry, Amoebas and Polytopes".