Ewelina Rupnik

Marc Pierrot Deseilligny

TOWARDS STRUCTURELESS BUNDLE ADJUSTMENT WITH TWO-AND THREE-VIEW STRUCTURE APPROXIMATION

Keywords: Global SfM, Bundle adjustement, Structure approximation

Towards structureless bundle adjustment with two-and three-view structure approximation

INTRODUCTION

Photogrammetry and computer vision Structure from Motion (SfM) algorithms underwent a remarkable evolution during the past two decades. Much of its evolution was driven by the advances in sensor technology, growing computer capabilities and democratisation of the fields through different software solutions [START_REF] Snavely | Photo tourism: exploring photo collections in 3D[END_REF], [START_REF] Pierrot Deseilligny | Apero, an open source bundle adjusment software for automatic calibration and orientation of set of images[END_REF], [START_REF] Moulon | Open-MVG: Open multiple view geometry[END_REF], [START_REF] Schonberger | Structure-from-motion revisited[END_REF]. Existing SfM approaches for pose estimation become computationally inefficient as the number of images increases. Two solutions exist: sequential and global methods. Sequential methods [START_REF] Schonberger | Structure-from-motion revisited[END_REF][START_REF] Snavely | Photo tourism: exploring photo collections in 3D[END_REF] generate one or multiple (e.g. hierarchical methods [START_REF] Klopschitz | Robust incremental structure from motion[END_REF][START_REF] Toldo | Hierarchical structure-and-motion recovery from uncalibrated images[END_REF]) seed images and sequentially concatenate overlapping images. Non-linear least squares bundle adjustment is then used to systematically eliminate accumulated errors [START_REF] Triggs | Bundle adjustmenta modern synthesis[END_REF][START_REF] Förstner | Photogrammetric computer vision[END_REF]. On the other hand, global methods [START_REF] Govindu | Combining two-view constraints for motion estimation[END_REF][START_REF] Moulon | Global fusion of relative motions for robust, accurate and scalable structure from motion[END_REF][START_REF] Wilson | Robust global translations with 1DSFM[END_REF] simultaneously and independently find relationships between pairs or triplets of images. These relationships are encoded in an epipolar graph which forms the basis of inferring the global rotations and translations. Global methods require the knowledge of image correspondences only in the first phase when computing the two-view or three-view geometries. This makes them faster and more computationally efficient. However, identifying outliers becomes challenging because the computation is bound to a pair or a triplet. Global solutions to SfM problems provide only initialisation of camera poses. Ordinarily, a final bundle adjustment (BA) including poses and point correspondences is run to calculate the optimal solution [START_REF] Triggs | Bundle adjustmenta modern synthesis[END_REF], which is a costly processing. The objective of this work is to go towards truly structureless bundle adjustment by refining the initial poses with very few points. We propose a structure (i.e. 3D points) approximation algorithm that relies on pairs and triplets of images. We firstly introduce a robust method for the relative motion calculation. Then, we redefine the per-relative-motion structure through the intermediary of an ellipsoid. The goal here is to generate a minimal set of 3D points that are sufficient to infer the local structure. This set is later used to refine the global solution in a bundle adjustment routine. Please note that we do not present a global motion approach in this work. Instead, to produce an initialization to the BA, we use the global approach of SfmInit [START_REF] Wilson | Robust global translations with 1DSFM[END_REF] or MicMac [START_REF] Pierrot Deseilligny | Apero, an open source bundle adjusment software for automatic calibration and orientation of set of images[END_REF][START_REF] Rupnik | MicMac-a free, open-source solution for photogrammetry[END_REF]. We evaluate our approach on the Strecha benchmark [START_REF] Strecha | On benchmarking camera calibration and multi-view stereo for high resolution imagery[END_REF], as well as sequential acquisitions by a mobile mapping platform and a drone. All methods are implemented and available from: github.com/micmacIGN . The sequential datasets and their respective ground truth are also available upon request.

RELATED WORK

On global motion methods Determination of global motions is performed separately for the rotational and translation parts [START_REF] Govindu | Combining two-view constraints for motion estimation[END_REF]. To calculate the rotations, the existing solutions: ignore the orthogonality constraint and employ linear methods [START_REF] Martinec | Robust rotation and translation estimation in multiview reconstruction[END_REF][START_REF] Arie-Nachimson | Global motion estimation from point matches[END_REF][START_REF] Moulon | Global fusion of relative motions for robust, accurate and scalable structure from motion[END_REF]; or optimise with robust 1 as Hartley et al. [START_REF] Hartley | L1 rotation averaging using the weiszfeld algorithm[END_REF] or as Chatterjee et al. [START_REF] Chatterjee | Efficient and robust largescale rotation averaging[END_REF] who base their algorithm on the Lie Group. A review of rotation averaging algorithms was published by [START_REF] Hartley | Rotation averaging[END_REF]. The translation component is always computed in the second place. Unlike for the rotational component, there is no direct way to compute scaleconsistent translations from a set of available pairwise constraints. Some methods exploit only pairwise constraints [START_REF] Govindu | Combining two-view constraints for motion estimation[END_REF], [START_REF] Martinec | Robust rotation and translation estimation in multiview reconstruction[END_REF], [START_REF] Arie-Nachimson | Global motion estimation from point matches[END_REF], [START_REF] Wilson | Robust global translations with 1DSFM[END_REF] but are susceptible to a degenerate solution when all cameras are aligned. Other methods formulate the problem using trifocal dependencies [START_REF] Courchay | Exploiting loops in the graph of trifocal tensors for calibrating a network of cameras[END_REF][START_REF] Zach | Disambiguating visual relations using loop constraints[END_REF][START_REF] Jiang | A global linear method for camera pose registration[END_REF][START_REF] Moulon | Global fusion of relative motions for robust, accurate and scalable structure from motion[END_REF]. Moreover, approaches based on linear solvers [START_REF] Arie-Nachimson | Global motion estimation from point matches[END_REF], [START_REF] Govindu | Combining two-view constraints for motion estimation[END_REF] and ∞ norm optimisation [START_REF] Sim | Removing outliers using the l\infty norm[END_REF][START_REF] Kahl | Multiple-View Geometry Under the L∞ -Norm[END_REF][START_REF] Moulon | Global fusion of relative motions for robust, accurate and scalable structure from motion[END_REF] can be distinguished. Methods formulating the BA estimation problem in terms of epipolar constraints [START_REF] Steffen | Relative bundle adjustment based on trifocal constraints[END_REF], [START_REF] Rodriguez | GEA optimization for live structureless motion estimation[END_REF], [START_REF] Indelman | Incremental light bundle adjustment[END_REF], [START_REF] Cefalu | Structureless bundle adjustement with self-calibration using accumulated constraints[END_REF] are not discussed herewith.

On outliers in relative motions Relative motions are the main ingredient in global pose estimation methods. This optimization scheme is particularly sensitive to outliers due to a limited number of motion observations. Since the motions are derived from automated algorithms for sparse correspondence search (e.g. SIFT [START_REF] Lowe | Distinctive image features from scaleinvariant keypoints[END_REF]), they inherently contain outliers, especially across ambiguous scenes with repetitive patterns, low texture, moving objects, etc. Existing approaches cope with the outliers by either eliminating them prior to optimising for global rotations and translations (i.e. filtering), or by adopting robust optimisation techniques. Examples of the latter are robust 1 optimization used for rotation averaging [START_REF] Hartley | L1 rotation averaging using the weiszfeld algorithm[END_REF][START_REF] Chatterjee | Efficient and robust largescale rotation averaging[END_REF] and translations solving [START_REF] Dalalyan | l1-Penalized robust estimation for a class of inverse problems arising in multiview geometry[END_REF][START_REF] Olsson | Outlier removal using duality[END_REF][START_REF] Zhu | Parallel structure from motion from local increment to global averaging[END_REF]. Some recent works have also addressed the evident complementarity of sequential and global methods. In Zhou et al. [START_REF] Zhu | Parallel structure from motion from local increment to global averaging[END_REF][START_REF] Zhu | Very large-scale global sfm by distributed motion averaging[END_REF] many local incremental SfMs help to discard mismatches and invalid motions. Among the solutions based on filters, the pioneering works [START_REF] Govindu | Robustness in motion averaging[END_REF], [START_REF] Sim | Removing outliers using the l\infty norm[END_REF], [START_REF] Martinec | Robust rotation and translation estimation in multiview reconstruction[END_REF] focus on pairwise constraints and remove the outliers by: RAN-SAC epipolar sampling [START_REF] Govindu | Robustness in motion averaging[END_REF]; using iterative removal techniques [START_REF] Sim | Removing outliers using the l\infty norm[END_REF]; or identifying false matches via partial reconstructions [START_REF] Martinec | Robust rotation and translation estimation in multiview reconstruction[END_REF]. Later works exploit the redundancy in the epipolar graph by looking at triplets of images. [START_REF] Courchay | Exploiting loops in the graph of trifocal tensors for calibrating a network of cameras[END_REF] introduce a trifocal graph parametrization and estimate robust global poses using linear programming with imposed loop constraints and RANSAC. Similarly to this work, Zach et al. [START_REF] Zach | Disambiguating visual relations using loop constraints[END_REF] uses loopy belief propagation to identify erroneous geometries in loops. [START_REF] Jiang | A global linear method for camera pose registration[END_REF] establishes a trifocal coplanarity constraint that minimizes a geometric error. Work in [START_REF] Moulon | Global fusion of relative motions for robust, accurate and scalable structure from motion[END_REF] fuses the Bayesian inference of Zach et al. [START_REF] Zach | Disambiguating visual relations using loop constraints[END_REF] and weighted graph filtering of Enqvist et al. [START_REF] Enqvist | Non-sequential structure from motion[END_REF] to remove outliers. Sweeney et al. [START_REF] Sweeney | Optimizing the viewing graph for structure-from-motion[END_REF] updates the fundamental matrices by optimizing a cost function defined over the trifocal point transfer. Finally, [START_REF] Wilson | Robust global translations with 1DSFM[END_REF] finds a purely pairwise method adapted to very large image collection where pairwise translations are projected to 1D subspace and a global ordering fitting the pairwise constraints is sought.

CONTRIBUTION AND OVERVIEW

To increase the robustness of the relative motions, our method combines several hypotheses, and tests each with a suitable direct algorithm (Section 3.1). To infer the triplets we take a similar perspective as [START_REF] Jiang | A global linear method for camera pose registration[END_REF] and [START_REF] Sweeney | Optimizing the viewing graph for structure-from-motion[END_REF]. We present a linear algorithm that rapidly calculates the triplet motions (Quick triplet algorithm). The robust triplet algorithm (Section 3.2) tests two hypotheses to see whether the camera in question has a short, normal, or a long focal length. As the long focal lengths approach the orthogonal projection, we replace the perspective camera model with the one proposed by [START_REF] Tomasi | Shape and motion from image streams under orthography: a factorization method[END_REF]. Otherwise, the pose of the "third" image of the triplet is found with the spatial resection [START_REF] Grunert | Das pothenotsche problem, in erweiterter gestalt, nebst bemerkungen über seine anwendung in der Geodäsie[END_REF]. We intentionally avoided the trifocal tensor parametrization knowing its unstable behavior on flat surfaces [START_REF] Julià | A critical review of the trifocal tensor estimation[END_REF]. Our structure reduction approach (Section 4) approximates the original points with an ellipsoid, and generates a small set of virtual points within its volume. This reduction allows for a lighter BA, without having to resort to PCG-based optimiza-tions [START_REF] Agarwal | Bundle adjustment in the large[END_REF][START_REF] Jeong | Pushing the envelope of modern methods for bundle adjustment[END_REF][START_REF] Kushal | Visibility based preconditioning for bundle adjustment[END_REF]. Therefore, we keep the access to the reduced camera system's covariance matrix. This is advantageous for many metrological applications that necessitate a degree of confidence associated with their measurements.

BUILDING THE RELATIVE MOTIONS

The computation of relative motion consistent within triplets of images is divided into three stages. Firstly, given a set of image correspondences, we create the epipolar hypergraph H(V, E), where each vertex Vi ∈ V represents an image, and the edge Eij connects a pair of images (Vi, Vj) (Section 3.1). The connection implies that relative motion {Rji, trji} is known. We refer to H(V, E) as a hypergraph because we have access to a list of connected vertices V k for a given Eij. The connected vertices are then explored to construct triplets complying with the respective edges (Section 3.2). We assume that cameras are calibrated, i.e. at least their focal lengths and the principal points are known. The local coordinate frame within an edge Eij or a triplet (Vi,Vj,V k ) is always associated with the first camera (i.e. Vi). In the text we often refer to the "third" image of a triplet by which we mean the image associated with V k .

Pairwise relative motions

Several direct algorithms are tested to compute robust relative motions between pairs of images. We begin with the hypothesis that the corresponding points contain a many outliers, and slowly relax the constraint by allowing more points in the computation. The full set of correspondences entering the method amounts to 500 and was chosen randomly. In total we test 6 calculation variants, see Table 1. The difference between variant T1 and T2 is in the way the points a drawn out the set of 500 points. While T2 is purely random, in T1 we bias the selection by forcing that the points are well distributed across the image plane. This way we avoid sampling points that are clustered in one zone. Variants T1-T4 are solved with 8-point algorithm [START_REF] Hartley | Multiple view geometry in computer vision[END_REF] throughout several RAN-SAC samples. After each draw the current solution is refined with a 2 solver. In T5 we employ the full set of points, and perform 1 and 2 estimations. The T6 variant tests the planarity of the scene by computing a homography and decomposing it to a relative motion [START_REF] Faugeras | Motion and structure from motion in a piecewise planar environment[END_REF]. Finally, the variant with the smallest re-projection error is refined in a coplanarity-based bundle adjustment.

Triplet motions

The triplet motions are determined ideally between "relevant" triplets of images. In order to select a subset of suitable triplets, MEss is the 8-point essential matrix algorithm, H is the homography decomposed to a relative motion [START_REF] Faugeras | Motion and structure from motion in a piecewise planar environment[END_REF]. (*) indicates that the points were drawn with a bias forcing a homogenous distribution within the image.

Variant NP ts NRANSAC 1 / 2 Model T1
we first determine the triplet motions between all possible candidates (see Quick triplet algorithm and Algorithm 1), and use them to define a per-triplet quality index (see Triplet selection).

The final triplet motions, calculated on a pre-selected set, result from the Robust triplet algorithm.

On straight lines intersection For the sake of clarity, we lay out the line intersection algorithm that is the building block of the Quick triplet algorithm. Let us define two straight lines as (cf. Figure 1 (a)):

P = P0 +p(P1 -P0), Q = Q0 +q(Q1 -Q0).
The intersecting 3D point belongs to the line that minimizes the distance between the lines P and Q, hence, it will be found on a line perpendicular to them. We define it as (P -Q), and impose the following:

(P -Q)•(P1-P0) = 0, (P -Q)•(Q1-Q0) = 0.
Using the above, parameters p and q are inferred. The final 3D point is equal to 1 2 (PInt p + PInt q ), where

PInt p = P0 +(P1 - P0) • p, PInt q = Q0 + (Q1 -Q0) • q.
Quick triplet algorithm The rotation of the "third" image in the triplet (Vi,Vj,V k ) is calculated twice: directly as R 3 = R31, and indirectly as R 3 = R21R32. The definite rotation being

R3 = 1 2 (R 3 + R 3 ).
The orthogonality of R3 is enforced by mapping it to its nearest orthogonal rotation with Singular Value Decomposition. The computation of the perspective center proceeds in two steps as presented in Figure 1 and Algorithm 1. Our objective is to force the perspective center to lie close to both the translation direction tr31 calculated in the preceding step, and the directions associated with the position of 3D points. To achieve this, our first step is to determine the 3D position of two corresponding points using image measurements in (Vi, Vj) (cf. Figure 1 (b)). Then, the new perspective center of V k results from intersecting the translation direction tr31 and the direction tr3P Int , defined by the image measurement in V k and attached at the position of the 3D point. For N corresponding points, one obtains N estimates of the perspective center (or N estimates of p and q, see Figure 1 (a)). The terminal value is then C3 = tr31 • p, where p = {pi}, p is the median of p and i ∈ [1, N ]. By definition, the proposed algorithm calculates consistent triplets exclusively from image correspondences of manifold 3 (i.e., points visible in 3 images). Statistically speaking, higher manifold points are less susceptible to outliers, therefore they are expected to produce more stable results. It is also required that a minimum of 8 3-manifold points exists. In practice, many more points are available and we decimate them by a factor of 50 when the number exceeds 500. Note that in this formulation, the computation of the perspective center is not degenerate if all three perspective centers are located on a straight line.

Algorithm 1 The Quick triplet algorithm

1: Input: a triplet motion {[R21|t21] , [R31|t31] , [R32|t32]} VP1, VP2, VP3 vector of points' correspondences N = VP1.size() = VP2.size() = VP3.size() 2: C1 = [0, 0, 0], R1 = 1 3: C2 = t21, R2 = R21 4: for i = 1 to N do 5: U1 = VP1 [i] 6: U2 = R2 • VP2 [i] 7: U3 = R31 • VP3 [i] 8: PInt = Intersect(C1 + U1, C2 + U2) 9:
{p, q} = Intersect pq(t31, PInt + U3) 10:

Vp.push(p) 11: end for Triplet selection Out of a number of initialized triplets, we want to select the K best ones. The notion of a best triplet is translated to a quality index Q and calculated for each triplet in H. The Q favours wide baselines between images (i.e. base-toheight ratios or bh) as presented here:

12: C3 = t31 • median(Vp) 13: R3 = N earestRot (0.5 • (R31 + R21 • R32))
R = bh bh + T L , (1)
Q = min(RV i V k , RV j V k ), (2)
where T L is a rough bh limit; RV i V k is computed between the first and the second image in a triplet, and RV j V k between the second and the third, correspondingly. To rapidly calculate the indices we exploit the triplet motions obtained with Quick triplet algorithm. We performed all the experiments with K = 1 as with the growing K we didn't get much improvement in accuracy but worsened the running times. The T L wes set to 0.15.

Robust triplet algorithm At this stage every edge Eij in H contains a list of at most K "third" images that were selected in the preceding step. We begin by re-estimating the poses of the "third" images for each Eij in H. This estimation is embedded in a RANSAC framework, and uses two direct estimation models: the spatial resection algorithm [START_REF] Grunert | Das pothenotsche problem, in erweiterter gestalt, nebst bemerkungen über seine anwendung in der Geodäsie[END_REF], and the Tomasi et al.approach to orientating images with long focal lengths [START_REF] Tomasi | Shape and motion from image streams under orthography: a factorization method[END_REF]. We always test both algorithms and choose the better result by comparing their reprojection errors. We set the number of RANSAC draws to 100, and take a random subset of 500 and 30 corresponding points for spatial resection and long focal length algorithms, respectively. Finally, a per-triplet bundle adjustment is run on the randomly reduced points.

STRUCTURE APPROXIMATION

Figure 2 illustrates a toy example of the structure approximation algorithm. Given a set of initial 3D points resulting from per-pair or per-triplet intersection (see On straight lines intersection), we calculate an ellipsoid inscribed in the points. The ellipsoid is represented by its eigenvectors, eigenvalues as well as its center of gravity. These parameters serve to generate a set of new, fictitious 3D points and their corresponding image measurements. In the same vein, Mayer [START_REF] Mayer | Efficient hierarchical triplet merging for camera pose estimation[END_REF] suggested a point reduction approach to speed-up the merging process in hierarchical SfM.

The ellipsoids are established for both triplets and pairs of images. To begin with, the center of gravity µ P and the covariance matrix K P are determined:

µ P = 1 N i=1 P dsP i • N i=1 Pi • P dsP i , (3)
K P = 1 N i=1 P dsP i • N i=1 P ds 2 P i • PiP T i -µ P • µ T P , (4)
where P = {Pi} are the point correspondences, P ds is the weight function that penalizes large reprojection errors (er) and small base-to-height ratios (bh), P ds = 1/ 1 + er α•bh 2 . In our experiments α was fixed to 10. Then, the eigenvalues and eigenvectors of the K matrix are retrieved with the Jacobi method [START_REF] Press | Numerical recipes in C++[END_REF]. We choose to approximate the structure with 5 symmetrically distributed points (see Figure 2 (d)) because it is the minimal number of points sufficient for direct pose estimation algorithms such as the 5-point algorithm (essential matrix), spatial resection as well as for a similarity transform that brings two stereo reconstructions to a common coordinate frame (e.g. in hierarchical SfM). Given the eigenvalues e and eigenvectors Ei, we generate a new point Hi from:

Hi = µ + e1E1 • hi 1 + e2E2 • hi 2 + e3E3 • hi 3 (5)
where µ is the ellipsoid's center of gravity; E1, E2, E3 are first, second and third unit eigenvectors and e1, e2, e3 are their corresponding eigenvalues; [hi 1 , hi 2 , hi 3 ] is the position of a fictitious point in the coordinate frame of the eigenvectors. In our 5-point configuration from Figure 2 (d), the following positions were used:

h1 = [-1, 1, 1], h2 = [1, -1, 1], h3 = [1, 1, -1], h4 = [0, 0, 0], h5 = [-1, -1, -1].
After the first generation of the new structure, we verify that its eigenvectors and eigenvalues were preserved. If they deviate from their initial values, we recalculate the new structure using an appropriate correction factor. Inspired by Mayer [START_REF] Mayer | Efficient hierarchical triplet merging for camera pose estimation[END_REF], we additionally extended our experiments to generating 5 randomly distributed points. The reported results are always a mean value over 100 or 1000 repeated random calculations. In the bundle adjustment, the reduced points are weighted according to the number of points that contributed to calculating the ellipsoid using: [START_REF] Lowe | Distinctive image features from scaleinvariant keypoints[END_REF], Ep, Ep+tri and Etri are respectively the fictitious image points derived from structure approximated on pairs, pairs and triplets, as well as triplets only.

P ds = 1 -N b M ax N b+N b M ax ,

RESULTS

We run experiments on three groups of datasets: the Strecha benchmark [START_REF] Strecha | On benchmarking camera calibration and multi-view stereo for high resolution imagery[END_REF], a sequential dataset (CAR) acquired by a mobile mapping system, and a drone dataset (UAV).

The number of images, number of corresponding points and number of triplet pairs per dataset is shown in Tables 2 and3.

In Tables 4,6 and Figure 3 the results are evaluated with respect to ground truth data, and with respect to other software solutions. To bring the results to the same coordinate frame, we perform a 7-parameter transformation on the cameras' perspective centers. The Average positional error refers to respective perspective centers' differences once the transformation has been applied. Strecha images are provided with their ground truth poses. The CAR dataset contains 647 images that amount to an approximately 2km trajectory without loops (see Figure 4). The camera ground truth poses were calculated using highly accurate ground control points measured in the field with classical surveying techniques (point positional accuracy in the range of few mm). To avoid having to measure many ground control points across the entire trajectory length, the vehicle moved in circles around a block of buildings, and during processing the image correspondences were extracted only between immediate cameras. The SfmInit [START_REF] Wilson | Robust global translations with 1DSFM[END_REF] solution did not suceed in orientating the CAR image set, therefore, we needed to resort to an alternative global SfM available in MicMac [START_REF] Pierrot Deseilligny | Apero, an open source bundle adjusment software for automatic calibration and orientation of set of images[END_REF]. A possible explanation is that the CAR dataset is a linear acquisition, with potentially many cameras located on the same 3D line. Such acquisitions form a degenerate case for algorithms based only on pairwise constraints. For the same reason the BA refinement on structure approximated with pairs of images (Ep) did not converge.

The UAV dataset contains 73 images over an area of 150x160m as illustrated in Figure 5. Analogously to the CAR dataset, the ground truth poses were calculated using ground control points measured in the field. Here as well, SfmInit did not manage to provide an initialisation to the final BA, hence, we resorted to MicMac [START_REF] Pierrot Deseilligny | Apero, an open source bundle adjusment software for automatic calibration and orientation of set of images[END_REF].

The followed general pipeline includes:

1. Extraction of correspondences with SIFT [START_REF] Lowe | Distinctive image features from scaleinvariant keypoints[END_REF];

2. Building of relative motions (Section 3);

3. Approximating the structure (Section 4), variants based on:

(a) pairs, Ep;

(b) triplets, Etri;

(c) pairs and triplets, Ep+tri;

(d) pairs and triplets with random point distribution, E rnd p+tri ;

4. Computation of global motions with SfmInit [START_REF] Wilson | Robust global translations with 1DSFM[END_REF] 3).

Running times Tables 5,6 and Figure 3 report on the running times. All our experiments were carried out on a machine with Intel Core i7-8550U, 1.8GHz x 8 processor. The Ours approaches concern time spent on the final BA. To deduce the integral processing time, one should add it to the time spent on calculating the initial solution. For a just evaluation, in Table 5, we look at the ratio Sf mIBA Ours , which compares the BA processing times with SIFT to our approximated structure approach. We are faster in all instances by a few factors.

Structure approximation variants

Approximating the structure with triplets only (Etri), or triplets and pairs (Ep+tri) of images turns out to be the best solution throughout all processed acquisitions. We can assume that for well connected images, which is the case for the three acquisitions tested, using triplets only is accurate enough and fast. For sparser connections, adding pairs may be indispensable. We have also observed a mediocre performance when the structure was approximated only with the pairs (Ep). In general, the refinement step succeeds in improving accuracy with respect to initial poses, and the variants with structure approximated by triplets of images perform at least as good as the sequential SfMs solution.

Deterministic or random fictitious points The results of the three experiments using random points do not stand out from the 5-point configuration. On the Strecha benchmark as well as the UAV dataset, the random distribution is only slightly worse than the best solutions. On the other hand, in the CAR dataset, the random selection sometimes outperforms the other three variants. We perceive that what discriminates the acquisitions is the scene geometry. In the CAR dataset, the scene is characterised by very large depth variations, which allows it to fully model the local structure even of randomly sampled.

Limitations

The structure approximation approach is not meant for image sets with unstructured connectivity graphs (e.g. Internet photo collections). The overwhelming redundancy of images does not reduce the information content across the scene, hence, does not accelerate the processing.

CONCLUSION

There are two leading contributions in this publication. First, we introduce a combinatorial approach to estimating relative motions which enhances the robustness to identified degenerate cases. Second, for a given relative motion, we propose a way of abstracting its 3D structure, such that the motion's estimation properties are preserved. The abstracted structure injected into a bundle adjustment effectively refines the initial solution at very low computational cost. The concept is validated with respect to ground truth data as well as other software solutions. Future work will concentrate on the purely structureless bundle adjustment, with exclusively view-dependent constraints and a thorough propagation of the minimal structure via relative and absolute motion covariance matrices. [START_REF] Pierrot Deseilligny | Apero, an open source bundle adjusment software for automatic calibration and orientation of set of images[END_REF], [START_REF] Schonberger | Structure-from-motion revisited[END_REF] SfMs, as well as our approach (Ours). SfmI is the initial solution found with [START_REF] Wilson | Robust global translations with 1DSFM[END_REF] that we relied upon in our BA. The E rnd p+tri corresponds the solution with structure approximated by pairs an triplets with randomly distributed points; a mean value over 1000 repeated calculations is reported; Ep, Etri, Ep+tri are the solutions where the structure is approximated with 5 points (see Figure 2 [START_REF] Pierrot Deseilligny | Apero, an open source bundle adjusment software for automatic calibration and orientation of set of images[END_REF], and in red continuous like it is Colmap's SfM [START_REF] Schonberger | Structure-from-motion revisited[END_REF]. MicMac and Colmap results were computed with with SIFT [START_REF] Lowe | Distinctive image features from scaleinvariant keypoints[END_REF]. In variant E rnd p+tri a mean value over 100 repeated calculations is reported. The high frequency drift is due to the cyclic nature of the acquistion geometry.

Figure 1 .

 1 Figure 1. Computation of the perspective center of the "third" image (C3) in the Quick triplet algorithm. C1 and C2 are perspective centers of images 1 and 2, respectively. Relative orientations between images 1 and 2 as well as 1 and 3 are known and related by a 7-parameter transformation. (a) A toy example illustrating the intersection of two straight lines defined by vectors ---→ P0P1 and ---→ Q0Q1 (cf. On straight lines intersection); (b) intersection of two lines defined by their image observations and respective perspective centers, C1 and C2; (c) the sought C3 lies at the intersection of the vectors tr31 and tr3P Int .

Figure 2 .

 2 Figure 2. A toy example of structure approximation in 2D and 3D, where (a) shows the initial 2D structure, (b) shows the fitted ellipse, (c) the fictitious 2D structure, and (d) is the equivalent in 3D.

Figure 3 .

 3 Figure 3. CAR dataset. Average position error (err) along the trajectory w.r.t. ground truth for several structure approximation variants: E rnd p+tri in black dashed line, Ep+tri in continuous magenta line, Etri in dotted blue line. MM (green dashed line) and MMInit (green dash-dotted line) are the sequential and global SfMs in MicMac [Pierrot Deseilligny, Cléry, 2011], and in red continuous like it is Colmap's SfM[START_REF] Schonberger | Structure-from-motion revisited[END_REF]. MicMac and Colmap results were computed with with SIFT[START_REF] Lowe | Distinctive image features from scaleinvariant keypoints[END_REF]. In variant E rnd p+tri a mean value over 100 repeated calculations is reported. The high frequency drift is due to the cyclic nature of the acquistion geometry.

Figure 4 .

 4 Figure 4. Left: excerpt from the CAR dataset. Right: (a) and (c) is the ground truth in top and side view; (b) and (d) is the equivalent Ep+tri result. Note the multiple facade effect in (b) and (d) due to error accumulation, and leads to a low frequency drift, as plotted in Figure 3.

Figure 5 .

 5 Figure 5. Left: excerpt from the UAV dataset. Right: (a) and (c) is the ground truth in oblique and side view; (b) and (d) is the equivalent SfmInit result. Note the erroneous initialisation of SfmInit in (b) and (d) where exclusively pairwise constraints are used.

Table 2 .

 2 where N b are the initial points, and N bMax was set to 100. Number of images (Nimages) and triplets (Ntri) in respective datasets.

		Dataset Nimages Ntri
		F-P11	11	52
		HJ-P8	8	18
		HJ-P25	25	228
		C-P10	10	44
		C-P19	19	68
		C-P30	30	214
		CAR	647	4167
		UAV	73	440
		Image correspondences		Ratio
	Dataset	Init.	Ep+tri	Etri	Ep	Init p+tri
	F-P11	31929	1796	1432	364	18
	HJ-P8	17084	736	540	196	23
	HJ-P25	96155	8158	6638	1520	12
	C-P10	25704	1628	1296	332	16
	C-P19	40008	2396	1850	546	16
	C-P30	148695	7726	6212	1514	19
	CAR	2175980 157878 128688	-	13
	UAV	170824	16388	13046	3342	10

Table 3 .

 3 Number of extracted features in respective datasets. Init are the features extracted with SIFT

Table 4 .

 4 Strecha benchmark. Average position error (err) w.r.t ground truth for different sequential [Pierrot Deseilligny,

					Accuracy err [mm]
		SfmI	Colmap MM	SfmI	Ours	Ours Ours	Ours
	Dataset				After BA	E rnd p+tri	Ep	Etri Ep+tri
	F-P11	163.8	1.4	1.8	1.3	1.3	33.6	1.0	1.0
	HJ-P8	70.5	2.8	5.2	2.7	2.8	2.0	2.6	2.5
	HJ-P25	95.3	3.9	5.8	4.3	3.8	4.8	3.5	3.4
	C-P10	105.5	5.2	3.3	5.5	5.7	44.4	5.9	5.5
	C-P19	1630.8	31.3	16.3	20.0	55.0	350.1 21.3	56.2
	C-P30	1346.3	14.7	11.2	13.3	16.7	263.5 12.7	12.6

Table 5 .

 5 (d)), and structure is approximated with pairs, triplets and the fusion of both, respectively. Strecha benchmark. Running times τ in conjuction with the results presented in Table4. Note that the times for Ours concern only the BA step. Times markes with ( * ) correspond to the total incremental SfM.

						Running time τ [s]
			SfmI Colmap	MM	SfmI	Ours	Ours Ours	Ours
		Dataset				BA	E rnd p+tri	Ep	Etri Ep+tri
		F-P11 HJ-P8 HJ-P25 C-P10 C-P19 C-P30	3.0 2.9 3.1 3.0 3.2 3.3	4.3 * 2.5 * 18.2 * 4.0 * 8.2 * 24.3 *	11.2 * 6.8 * 40.2 * 9.4 * 14.5 * 195.0 *	1.9 1 8.9 1.7 8.3 9.7	0.7 0.4 2.5 0.9 1.7 4.5	0.5 0.3 1.5 0.5 0.7 1.4	0.6 0.4 1.6 0.6 0.8 2.4	0.7 0.4 2.1 0.8 1.5 3.8
		Accuracy	Running time
		err [mm]	τ [s]
	Colmap		17.5	25.0
	MM		27.2	39.1
	MMInit		41.0	50.0
	MMInit BA		30.0	61.0
	Ours	E rnd p+tri	21.0	9.8
	Ours	Ep	27.4	3.8
	Ours	Etri	18.2	5.3
	Ours	Ep+tri	17.0	7.2

Table 6 .

 6 UAV dataset. Average position error (err) and running times (τ ). MMInit and MM are the global and sequential SfMs in MicMac. MMInit BA corresponds to the BA refinement on the initial solution. Note that the times for Ours concern only the BA step. In E rnd p+tri a mean value over 1000 repeated calculations is reported.


Inge Dox
email: inge.dox@uantwerpen.be

Jožica Gričar

Lorène J Marchand

Sebastien Leys

Paolo Zuccarini

Charly Geron

Peter Prislan

Bertold Mariën

Patrick Fonti

Holger Lange

Josep Peñuelas

Jan Van Den Bulcke

Matteo Campioli

Timeline of autumn phenology in temperate deciduous trees

Keywords: autumn phenology, xylem formation, foliar senescence, sink limitation, chlorophyll, wood growth

Cessation of xylem formation or wood growth (CWG) and onset of foliar senescence (OFS) are key autumn phenological events in temperate deciduous trees. Their timing is fundamental for development and survival of trees, ecosystem nutrient cycling, the seasonal exchange of matter and energy between the biosphere and atmosphere and affect the impact and feedback of forests to global change. A large-scale experimental effort and improved observational methods have allowed us to compare the timing of CWG and OFS for different deciduous tree species in Western Europe, in particularly silver birch, a pioneer species, and European beech, a late-succession species, at stands of different latitudes, of different levels of site fertility, and for two years with contrasting meteorological and drought conditions i.e. the low-moderately dry 2017 and the extremely dry 2018. Specifically, we tested whether foliar senescence started before, after or concurrently with CWG. OFS and CWG occurred generally between late September and early November, with larger differences across species and sites for OFS. Foliar senescence started concurrently with CWG in most cases, except for the drier 2018 and, for beech, at the coldest site, where OFS occurred significantly later than CWG. Behavior of beech in Spain, the southern edge of its European distribution, was unclear, with no CWG, but very low wood growth at the time of OFS. Our study suggests that OFS is generally triggered by the same drivers of CWG or when wood growth decreases in late summer, indicating an overarching mechanism of sink limitation as a possible regulator of the timing of foliar senescence.

Introduction

Phenology, the study of cyclic and seasonal plant and animal life cycle events, is important for understanding the functioning of temperate deciduous forest ecosystems, especially in the context of global change [START_REF] Lui | Delayed autumn phenology in the Northern Hemisphere is related to change in both climate and spring phenology[END_REF][START_REF] Chen | What"s hot about mercury? Examining the influence of climate on mercury levels in Ontario top predator fishes[END_REF]. Most phenological studies of deciduous forest ecosystems have investigated the onset of the growing season, mainly by characterizing the timing of canopy development and cambium reactivation. Spring phenological events are relatively simple to observe and have clear responses to environmental change [START_REF] Peñuelas | Phenology -Responses to a warming world[END_REF][START_REF] Guada | Response patterns of xylem and leaf phenology to temperature at the southwestern distribution boundary of Quercus robur: A multi-spatial study[END_REF], which facilitates these studies. In contrast, autumn phenology has been largely overlooked A c c e p t e d M a n u s c r i p t [START_REF] Gallinat | Autumn, the neglected season in climate change research[END_REF], even though the timing of autumn phenophases affects not only tree development (i.e. winter survival) and growth potential for the next year [START_REF] Estiarte | Alteration of the phenology of leaf senescence and fall in winter deciduous species by climate change: effects on nutrient proficiency[END_REF], but also fundamental ecosystem dynamics such as nutrient cycling, biomass production, CO 2 and H 2 O fluxes between forest and atmosphere, albedo, seasonal weather and species interactions [START_REF] Peñuelas | Phenology Feedbacks on Climate Change[END_REF][START_REF] Richardson | Influence of spring and autumn phenological transitions on forest ecosystem productivity[END_REF][START_REF] Richardson | Climate change, phenology, and phenological control of vegetation feedbacks to the climate system[END_REF][START_REF] Gallinat | Autumn, the neglected season in climate change research[END_REF].

The paucity of attention and studies of autumn phenological events is not due to a failure of ecologists to recognize their importance, but rather due to assessment difficulties. Onset of foliar senescence (OFS) and cessation of xylem formation, i.e. the cessation of wood growth (CWG), are two key processes of deciduous trees in autumn. Leaves actively degrade their components during senescence, including chlorophyll and the photosynthetic apparatus. This degradation releases resources and relocates them into overwintering tree organs before foliar abscission. Average efficiencies of N, P and K resorption can be as high as 62, 65 and 70%, respectively [START_REF] Estiarte | Alteration of the phenology of leaf senescence and fall in winter deciduous species by climate change: effects on nutrient proficiency[END_REF]. Timing of OFS is currently determined from ecophysiological, visual or remote sensing observations of canopy coloration, chlorophyll content or loss of foliar biomass [START_REF] Mariën | Detecting the onset of autumn leaf senescence in deciduous forest trees of the temperate zone[END_REF]. However, a standard method for identifying this crucial autumn phenophase is not yet available.

CWG represents the end of the annual wood formation. More specifically, wood growth ends when the production and maturation of xylem cells is completed [START_REF] Plomion | Wood formation in trees[END_REF]. Cell maturation includes the deposition of a secondary cell wall and cell wall lignification, followed by a programmed cell death. In contrast to onset of wood formation in spring (cell production and expansion), which can be studied by changes in stem diameter, observation of cell maturation requires laborious microscopic anatomical analyses [START_REF] Deslauriers | Ecophysiology and plasticity of wood and phloem formation[END_REF].

The few studies on CWG, foliar senescence and related processes (specifically, chlorophyll degradation, canopy coloration and leaf fall) have not established a clear timeline for the onset and progress of these processes in temperate deciduous trees. These studies either did not measure CWG in detail [START_REF] Fracheboud | The Control of Autumn Senescence in European Aspen[END_REF]Weih, 2009) or included foliar senescence only peripherally (González-González, 2013; Prislan, et al., 2013;[START_REF] Gričar | Intra-annual leaf phenology, radial growth and structure of xylem and phloem in different tree parts of Quercus pubescens[END_REF]. The time course and causal relationship between these processes remains unclear. A better understanding of the timeline between A c c e p t e d M a n u s c r i p t knowing their timelines would be valuable to understand the functional priorities of trees at the end of the growing season. Moreover, future projections of forest production and forest climate services can be improved by the elucidation of the autumn timeline in deciduous trees.

We tested three plausible hypotheses to account for the causal sequence between OFS and CWG. The first hypothesis assumes that foliar senescence starts before CWG, because photosynthesis may be more limited than wood growth in autumn due to accumulated damage to foliar cells throughout the season (Günthardt-Goerg & Vollenweider, 2007). This classical approach (Gallinat et al., 2015) has been widely applied in forest models and in the land surface component of the Earth systems models, where tree growth is often modeled to depend on the available photosynthates [START_REF] Fatichi | Moving beyond photosynthesis: from carbon source to sink-driven vegetation modeling[END_REF]. The second hypothesis assumes that foliar senescence starts after CWG as growth in autumn is more sensitive to environmental constraints than photosynthesis [START_REF] Fatichi | Moving beyond photosynthesis: from carbon source to sink-driven vegetation modeling[END_REF]Zuidema et al., 2018). The third hypothesis assumes that OFS and CWG occur concurrently, possibly because both phenophases respond to the same environmental cues, e.g. photoperiod. This hypothesis assumes that the triggering cues would not be perceived with different sensitivities by different tree organs. Three research questions were linked to these hypotheses: i) Do species follow the same hypothesis over their distributional range? ii) Do populations growing on sites with a different level of fertility or on the same site, but during years with contrasting meteorological conditions, follow the same hypothesis? iii) Do species with different traits respect the same hypothesis? To answer these research questions we tested the validity of the hypotheses for i) deciduous tree species of contrasting types, i.e. European beech (Fagus sylvatica L.), a late successional species that produces leaves mainly in spring, and silver birch (Betula pendula Roth.), a pioneer species that produces leaves throughout the whole growing season, for European populations at three locations (southern Norway, Belgium and northern Spain) of their distribution, spanning 20 degrees of latitude, ii) populations of these two species at the same location (Belgium), but with different soil fertilities and for two years (2017 and 2018) with different meteorological conditions and drought intensity, and iii) populations of other key European species, i.e. pendunculate oak (Quercus robur L., population in Belgium), which can be either a pioneer or a late-successional species, with A c c e p t e d M a n u s c r i p t leaf flush mainly in spring, and common aspen (Populus tremula L., population in Norway), a typical pioneer species, but with leaf flush mainly in spring for adult individuals.

The main analyses were based on detailed phenological observations of the timing of CWG and OFS for 64 mature trees during the autumn of 2017 and 16 trees during the autumn of 2018. CWG was determined through an improved method for an accurate determination of the end of xylem maturation based on microscopic analyses of wood cross-sections [START_REF] Gričar | Intra-annual leaf phenology, radial growth and structure of xylem and phloem in different tree parts of Quercus pubescens[END_REF]. OFS was determined with a breakpoint analysis applied to standard observations of canopy coloration (determined by changes in foliar color and loss of color as a result of leaf fall; [START_REF] Vitasse | Altidudinal differentiation in growth and phenology among populations of temperate-zone tree species growing in a common garden. The interaction between freezing tolerance and phenology in temperate deciduous trees[END_REF]; Methods) and leaf chlorophyll degradation [START_REF] Fracheboud | The Control of Autumn Senescence in European Aspen[END_REF]Methods). For 2018, an additional comparison between CWG and the end of xylem enlargement was done. We also compared our estimates of the OFS from canopy coloration and chlorophyll degradation to the timing of 50% canopy coloration and the timing when 50% of the leaves had fallen, which are commonly used proxies of OFS (Fracheboud et

Materials and Methods

Study species and sites

Stands of European beech (Fagus sylvatica L.) and silver birch (Betula pendula Roth.) growing in moderately fertile conditions ("fertile stands") were selected at three locations within their distributions: Ås, Norway (59°40′N, 10°46′E; 86 m a.s.l.; northern distribution), Brasschaat, Belgium (51°20′N, 4°30′E; 15 m a.s.l.; central distribution) and a mountainous area in Girona province, Spain (41°45′-42°11′N, 2°20′-2°28′E; 1075-1126 m a.s.l.; southern distribution; Fig. 1). Additional stands of the European beech and the silver birch in Brasschaat were selected at sites with low fertility ("infertile stands"). We also studied one fertile stand of common aspen (Populus tremula L.) in Norway and one fertile stand of pedunculate oak (Quercus robur L.) in Belgium.

We selected four (co)-dominant healthy individuals in each of the ten stands where we monitored OFS and CWG in the year 2017. One birch tree in Norway was excluded from the data set, because it was a suppressed understory and was erroneously selected. In the Belgian fertile stands of A c c e p t e d M a n u s c r i p t beech, birch and oak, four additional (co)-dominant healthy individuals that neighbored the trees selected for the study of OFS and CWG were selected for the study of foliar senescence and related processes (chlorophyll degradation, canopy coloration and leaf fall). For the fertile beech stand and the infertile birch stand, these 8 trees were studied during the year 2018.

Average daily temperature in Belgium was very similar over 2017 and 2018, with respectively 11.1 °C and 11.8 °C. Both years were warmer than the long-term average of 10.1°C (en.climatedata.org). In Belgium, 846 mm of precipitation fell in 2017 and 723 mm in 2018 compared to the longterm average of about 920 mm. In 2018, a severe spring/summer drought took place in Belgium with only 112 mm of precipitation during May, June and July (compared to long-term average of 212 mm).

The precipitation in May 2018 and the number of days with precipitation in July 2018 were "extraordinarily low" (recurrence time of 30 years). Both variables, precipitation and number of precipitation days, were extraordinarily low in June 2018 (meteo.be). Average temperature was respectively 7 and 11-13°C for the stands in Norway and Spain compared to long-term averages of respectively 5 and 11.1-15.5°C. Total precipitation in 2017 was 980 mm in Norway (long-term average: 785 mm) and 820-850 mm in Spain (long-term average: 880 mm). Additional information about the stands, trees and meteorological data for late summer and autumn 2017 and 2018 is provided in Fig. 2 and Suppl. Table 1.

Canopy coloration

Canopy coloration at day t (X t ) in late summer and early autumn was estimated directly as a percentage. Later in the season, when the process became more intense, X t was estimated by a combined rating of the percentage of leaves that had changed color (α t ) and the percentage of leaves that had fallen (β t ), as described by [START_REF] Vitasse | Assessing the effects of climate change on the phenology of European temperate trees[END_REF]:

𝑥 𝑡 = 𝛼 𝑡 × (100 -𝛽 𝑡 ) 100 + 𝛽 𝑡
Observations were recorded weekly from mid-September until 100% X t . Small proportions of the birch and aspen leaves in Norway were already colored (ca. 10%) at the start of the monitoring, so the canopy by tree climbers. For beech and oak, five sun-and five shade-leaves were sampled, whereas for birch only five sun-leaves were sampled as there was no clear difference between sun-and shade-leaves. Canopy coloration of the trees sampled for chlorophyll content was assessed as described above.

Determination of the onset of foliar senescence

The OFS was defined as the date when (i) seasonal canopy coloration sharply increased due to the start of natural foliar senescence (canopy coloration before this date had already slightly increased during the season due to stress damage to leaves) or (ii) chlorophyll content starts to decrease sharply.

The two phenophases are closely correlated and their measurements provide similar dates of OFS, except in case of severe drought (see above). We determined OFS by a breakpoint analysis of the time series of mean canopy coloration or chlorophyll content (see above). The breakpoint analysis was performed with the R package "segmented" [START_REF] Muggeo | segmented: An R Package to Fit Regression Models with Broken-Line Relationships[END_REF]. This package calculates multiple linear regressions for data with dependent variables that can be expressed by two or more straight lines with different slopes linked at a breakpoint and calculates these breakpoints (Fig. 3).

Wood growth cessation

Seasonal wood growth was monitored weekly from late August to late November in 2017 and, for two stands in Belgium, also from mid-March to late November in 2018. We measured the percentage of cells of the currently forming annual ring that were still in the wall thickening phase (and thus still developing and operating as a sink for assimilates) through microscopic analysis of micro-cores of stem wood (small samples 15-20 mm long and 2 mm in diameter). In 2018 this analysis was extended to the cells in the enlarging phase (preceding thickening) and in the mature phase (after thickening).

Sampling. The stem micro-cores were collected at breast height with a Trephor corer [START_REF] Rossi | Trephor: a new tool for sampling microcores from tree stems[END_REF]). An upward spiral pattern of sampling was followed throughout the season to minimize wound reactions. The samples were stored in 1. Microscopic analysis. The slides were examined under a light microscope (Leica DM 4000 B/M, Wetzlar, Germany) and a Leica LAS image-analysis system to determine the end of the formation of wood, defined as the point when the cells of late wood are lignified and have thus completed cell wall thickening [START_REF] Gričar | Intra-annual leaf phenology, radial growth and structure of xylem and phloem in different tree parts of Quercus pubescens[END_REF]. This stage is recognizable by a completely red cell wall through the safranin/Astra blue staining, whereas cells that have not completed wall thickening are colored blue (i.e. non-lignified) on the inner parts of the cell wall and red (i.e. lignified) on the outer part (Fig. 4).

The percentage of wall-thickening cells (WTC) was estimated weekly by measuring the width of the current year ring and the width of the layer of WTC. The width of the layer of the WTC can also be used directly as an index of growth (Suppl. Fig. 1), but we preferred to use the percentage of the WTC because of its value relative to annual growth. We defined CWG as the start of (at least) a 3-

A c c e p t e d M a n u s c r i p t U N C O R R E C T E D M A N U S C R I P T 9
week period with the percentage of WTC < 0.50%. This definition was adopted for multiple reasons.

First, the 0.50% threshold is an indication of severely reduced allocation of carbon to wood growth.

Second, using a 3-week period avoids associating the CWG with a sporadic short period when growth is limited due to a temporarily environmental constraint (e.g. drought). Third, this definition reduces the impact on the analysis of outliers and noise caused by minor heterogeneous growth patterns of the sampled tree stems. CWG was determined at the tree level, but then averaged (n=4-8) to calculate the stand value. For 2018, we performed measurements also for the whole growing season, recording the width of the current forming ring section with enlarging cells (enlarging zone), thickening cells (wallthickening zone) and mature cells (mature zone).

Statistical analysis

CWG and OFS. Significant differences in the timing between CWG and OFS were determined by identifying overlaps in confidence intervals based on the uncertainty method. The 95% confidence interval of OFS was derived for each stand from the mean (x) and standard error (SE) obtained from the "segmented" package (x-1.96SE, x+1.96SE), assuming a normal sampling distribution for the onset of senescence. The confidence interval for CWG was also derived from the mean and standard error of the mean (x-1.96SE, x+1.96SE). The difference between CWG and OFS was significant at p < 0.05 when the two confidence intervals did not overlap. Significant differences among OFS dates for different stands and years were determined in the same way considering the confidence interval of each OFS value. On the other hand, significant differences in CWG dates among stands in 2017 were obtained by one-way Anova followed by a post-hoc test with Tukey corrected p-values. For the Belgian stands, to check for CWG differences between the years 2017 and 2018 repeated measurements Anova were done.

Measurements methods of foliar senescence. Repeated measurements ANOVA"s, with an HC3 correction for non-homogeneity of variance ("car" package), followed by a Tukey"s post-hoc tests were was also performed to identify significant differences between OFS from canopy coloration data, OFS from the start of rapid chlorophyll degradation, 50% canopy coloration and 50% leaf fall, for fertile birch, beech and oak in Belgium.

Results

Onset of foliar senescence

In 2017, for birch, OFS occurred in late September in Norway and Spain but in mid-October in Belgium (p < 0.05; Fig. 5; Suppl. Table 6). For beech, OFS occurred in the same week in mid-October in Norway and Belgium and a week later in Spain (no sign. diff.). OFS in Belgium occurred at about the same time at the infertile and fertile stands for both species (no sign. diff.). Finally, OFS for oak in Belgium occurred in late October (2.5 weeks later than the local beech and birch, p < 0.05) and in late September for aspen in Norway (same week as the local birch; no sign. diff.). In 2018, OFS for the infertile birch stand and the fertile beech stand in Belgium occurred in late October (no sign. diff.). For these stands, difference in OFS between years was also not significant.

Cessation of wood growth CWG for both birch and beech occurred from late September to early October at all stands (no sign. diff.; Fig. 5; Suppl. Table 6). The results were thus more similar among stands than for OFS. Cell maturation (i.e. thickening of the cell wall and lignification) for beech at the Spanish stand, however, continued throughout the whole monitoring period, which lasted until early December. CWG also A c c e p t e d M a n u s c r i p t 6).

U N C O R R E C T E D M A N U S C

Seasonal pattern of wood growth

In 2018, the seasonal progress of the xylem enlargement zone in the infertile birch stand and fertile beech stand in Belgium showed an increasing width until May-June and a declining width afterwards, until complete enlargement in early August (Fig. 6). The thickening zone showed a similar trend but shifted, with an increasing width until late June -mid July, declining width until early September (Fig. 6) and full cessation of thickening in late September (Fig. 5). In more details for the thickening phase in autumn and for all stands, the percentage of WTC generally decreased rapidly to <0.50% in about 1-2 weeks after a period with a larger, generally stable, percentage of WTC from late August to mid-September. The percentage of WTC was 2-3-fold larger for beech than birch across locations before decreasing in mid-September. For each sampling date, no significant difference was observed in percentage of WTC between sites of different fertility for both birch and beech with the exception of Table 4). Similarly, no difference in percentage of WTC was observed between 2017 and 2018 in the monitored birch and the beech stands, with only a marginal difference in the first sampling date for birch (p = 0.070; Suppl. Table 5) The seasonal pattern (and rates) of wood growth in late summer and autumn did not change when the width of WTC was analyzed instead of their percentage (Suppl. Fig. 1).

Measurement methods of foliar senescence

In 2017, our detailed foliar measurements of autumn processes in Belgium indicated that OFS from data of chlorophyll degradation matched the estimates of OFS from data of canopy coloration (difference of 0.5-2.5 weeks; Fig. 7, Suppl. Table 2a; Suppl. Table 2b). On the other hand, our analyses indicated significant differences between OFS from canopy coloration and 50% threshold in canopy coloration for birch (p = 0.005) and between the former and 50% threshold in leaf fall for all species (p < 0.02). Concerning the estimates of OFS from chlorophyll degradation, they were found to be different from the 50% threshold in canopy coloration for beech and oak (p < 0.04) and the 50% threshold in leaf fall for all species (p < 0.01). In 2018, OFS from canopy coloration was recorded significantly earlier than in 2017 (p < 0.001; Suppl. Table 3; Suppl. Table 2b), presenting a significant difference with the OFS estimates from chlorophyll degradation (p < 0.001). However, early coloration in 2018 was related to drought stress and not to the natural process of leaf senescence [START_REF] Mariën | Detecting the onset of autumn leaf senescence in deciduous forest trees of the temperate zone[END_REF]. Therefore, OFS was derived from measurements of canopy coloration for the probable that CWG occurred shortly after our last sampling (or even shortly before it, as two sampling points were missed in late November 2017, Fig. 5). Formally, the pattern observed for beech in Spain follows Hypothesis 1. Practically, the wood growth at OFS was very low, as the percentage of WTC (1-3.5%), was substantially lower than the corresponding percentage in early September (13-15%).

The recorded pattern is inconsistent with previous findings from a study on European beech in Moncayo Park, Spain (del Castillo et al., 2016) where CWG occurred already in early August.

However, this stand is located on a drier location than our site (710 vs. 820 mm annual precipitation respectively), where wood growth may be constrained by drought. Also, as lignification is sensitive to A c c e p t e d M a n u s c r i p t

U N C O R R E C T E D M A N U S C R I P T 14
Site fertility and meteorological conditions. In Belgium, birch and beech had no different autumn phenologies at stands with different fertility, following hypothesis 3 in both fertile and infertile conditions. Similar phenology under different fertility levels was unexpected. In fact, as supported by fertilization experiments, growth is limited, and foliar senescence is expected earlier, under suboptimal nutrient conditions (Weih, 2009; [START_REF] Sigurdsson | Elevated [CO2] and nutrient status modified leaf phenology and growth rhythm of young Populus trichocarpa trees in a 3-year field study[END_REF]. However, as shown for other growth processes [START_REF] Campioli | Biomass production efficiency controlled by management in temperate and boreal ecosystems[END_REF], long-term adaptations might create differences between the influence of natural fertility and fertilization on plant growth. For the Belgian stands studied both in 2017 and 2018, drought did not have an effect on CWG and OFS. On the other hand, the difference between CWG and OFS became significant in 2018, following Hypothesis 2. As expected, the seasonal pattern of xylem production and enlargement was shifted compared to the seasonal pattern of xylem maturation, with the latter showing a later start and later end [START_REF] Rossi | A metaanalysis of cambium phenology and growth: linear and non-linear patterns in conifers of the northern hemisphere[END_REF]. However, while it started two weeks earlier, enlargement ended 3-4 weeks before the large majority of the maturation zone was completed for both the infertile birch stand and fertile beech stand in Belgium. In other words, the phase of enlargement was shorter than the maturation phase. This aspect, together with the fact that the enlargement zone showed a declining width after June (for birch even during June), confirms that maturation is less sensitive to drought than enlargement [START_REF] Arend | Seasonal change in the drought response of wood cell development in poplar[END_REF].

Concerning OFS, overall, our findings are comparable to the general range of values reported for the study species in multi-year field observations in areas comparable to our locations (e.g. [START_REF] Kint | Radial growth change of temperate tree species in response to altered regional climate and air quality in the period 1901-2008[END_REF][START_REF] Vitasse | Assessing the effects of climate change on the phenology of European temperate trees[END_REF].

Traits. Aspen had the same autumn phenology as birch. In comparison to beech, autumn phenology of oak was substantially later, with OFS and CWG occurring respectively 2.5 and 4 weeks later. Nevertheless, in Belgium, oak also followed Hypothesis 3, as did beech. We speculate that the low sensitivity to low temperature of the oak populations in Belgium [START_REF] Kint | Radial growth change of temperate tree species in response to altered regional climate and air quality in the period 1901-2008[END_REF] allows a longer period of wood growth for this species. Deciduous oak species also showed a longer cambial activity than conifers under stress conditions (drought) (Fernández-de-Uña et al., 2017, 2018).

Addressing our three research questions demonstrated that temperate tree species along their distributions, in general, have common dynamics of autumn phenology. For example, across the latitudinal gradient considered, OFS of beech took place within 10 days in mid-late October, OFS of birch in late September -mid October, while CWG for both birch and beech occurred from late September to early October at all stands. Notably we found that CWG and OFS were mostly concurrent (Hypothesis 3), or OFS occurred after CWG (Hypothesis 2) in harsher conditions e.g.

during the very dry 2018 (for both birch and beech) or at the coldest location (for beech only).

However, exceptions were observed, as the low but continued growth recorded at the time of OFS at the beech stand in Spain. This case is particularly interesting as, due to the elevation, the seasonal temperature regime of this site in 2017 was not substantially different than the one of northern locations e.g. in Belgium (Fig. 2). On the other hand, it is possible that southern genotypes might have more flexible CWG timing, as CWG may be more related to temperature rather than to photoperiod at southern locations. In fact, for leaf dynamics, [START_REF] Gill | Changes in autumn senescence in northern hemisphere deciduous trees: a metaanalysis of autumn phenology studies[END_REF] found that autumn phenology is more related to photoperiod at northern latitude and more to temperature at lower latitude. Also, we should note that only two years were studied and longer-term datasets (e.g. comprising a post-drought year, a year with an early cold spell) might provide more insights. Nevertheless, we are not aware of previous studies in which the final phase of seasonal wood growth, and its link to leaf senescence has been analyzed with so much detail as in this analysis.

Even if Hypotheses 2 and 3 are logical alternatives to each other, both may rely on the wellknown ecophysiological mechanism of sink limitation [START_REF] Fatichi | Moving beyond photosynthesis: from carbon source to sink-driven vegetation modeling[END_REF][START_REF] Fu | Variation in leaf flushing date influences autumnal senescence and next year"s flushing date in two temperate tree species[END_REF]. This principle states that photosynthesis is sustained by a tree only if the newly produced assimilates are consumed by plant organs or other sinks, such as woody tissues, which are important carbon sinks in late summer and autumn [START_REF] Cuny | Woody biomass production lags stem-girth increase by over one month in coniferous forests[END_REF]. The capacity of the carbon sink thus decreases when wood growth decreases, leading to the down-regulation of foliar photosynthesis. Application of the sink-limitation mechanism to autumn dynamics implies that foliar senescence is triggered when growth decreases in late summer (supporting Hypothesis 2 or 3). Our results are consistent with the mechanism of sink limitation (caused by CWG) for the timing of foliar senescence. Following a conservative interpretation, the anomalous case of the beech trees at the Spanish stand would represent an important exception to this pattern. However, the fact that wood growth at OFS in late October was substantially reduced than in early September, might not exclude a sink limitation even in this case. An experiment under completely different settings performed on birch and beech saplings in Belgium (Suppl. Fig. 2) further supports the sink-limitation mechanism. Monitoring OFS and CWG in branches of saplings grown in well-irrigated and fertilized pots indicated that foliar senescence started in early to mid-September after the rates of wood growth had decreased significantly, even though (slow) growth continued until mid-October.

Our methodology allowed us to detect subtleties in the autumn phenologies of temperate deciduous trees and to observe the elusive temporal relationship between CWG and OFS. Moreover, the overview of autumn phenology for birch, beech, aspen and oak provides the first detailed comparison of temperate deciduous trees for the timing of CWG and OFS. A better resolution of the timing of autumn processes is important for understanding tree functioning and potential causal relationships among autumn processes and for identifying their environmental drivers. These advances are key to improving simulations of tree growth in both forest models of biomass production and global models providing projections of climate change. In particular, modelling approaches simulating CWG after OFS should be amended. Improved modelling approaches should consider these phenophases as concurrent, with earlier CWG in case of extreme conditions. A c c e p t e d M a n u s c r i p t

A c c e p t e d M a n u s c r i p t we added a

  typical summer value (early August) measured from another year (2018) to the seasonal series of 2017. Chlorophyll degradation Chlorophyll content was estimated indirectly by measuring the relative chlorophyll content index (CCI) using a CCM-200 plus leaf absorption meter (Fracheboud et al., 2009; ADC Bioscientific Ltd., Hoddesdon, UK) every 2-3 weeks between 10:00 and 16:00 from late July to late November in 2017 and 2018. The CCI value was obtained by averaging measurements from leaves collected from

A c c e p t e d M a n u s c r i p t

  5-ml Eppendorf micro-tubes containing a 70:30 EtOH:H 2 O solution for conservation.Slide preparation. The micro-cores were cut to contain only the most recently created xylem rings (at least two), the cambium and the inner phloem. These segments were then dehydrated in a series of rising ethanol concentration, infiltrated with a clearing agent (UltraClear, J.T. Baker, Avantor Performance Materials, Center Valley, USA) and paraffin (Paraplast plus, ROTH, Karlsruhe, Germany) and embedded in paraffin blocks. The paraffin blocks were cut into 5-10 µm slices with a microtome (Leica Microsystems, Wetzlar, Germany). Next, the paraffin was removed with a clearing agent (UltraClear) and ethanol. The sections were then stained in an aqueous solution of safranin (Merck, Darmstadt, Germany) and Astra blue (Sigma-Aldrich, Steinheim, Germany). The slices were fixed by mounting in Euparal (Waldeck, Munster, Germany)[START_REF] Gričar | Intra-annual leaf phenology, radial growth and structure of xylem and phloem in different tree parts of Quercus pubescens[END_REF].

10 Seasonal

 10 pattern of xylem cells still in the wall thickening phase. A pairwise t-test with HC3correction was used for the seasonal change in percentage of WTC. Repeated measurements Anova, sometimes done after a log-transformation to ensure normality of the residuals and homogeneity of variance, were performed to compare the percentage of WTC for each sampling date. Similarly, t-tests were done to compare the percentage of WTC between fertile and infertile sites in Belgium. When necessary a Welch t-test was used or a log-transformation was performed to meet the assumptions.The normality of residuals and homogeneity of variances were tested with a Shapiro-Wilk test and a ncvTest or LeveneTest ("car" package), respectively. Significance was set at p < 0.05. All statistical analyses and graph construction used Rstudio version 3.2.2 (R Core Team 2016).

  difference in the first sampling event (DOY 233) at the birch site (p = 0.094; Suppl.

  2017 inter-population comparison and from measurements of chlorophyll degradation for the 2017-2018 inter-annual comparison.DiscussionDistributional range. Based on the comparison between OFS and CWG for all fertile stands of birch and beech, Hypothesis 3 is respected in four cases out of six, and always for birch. Beech in Norway follows Hypothesis 2, likely because wood growth is more sensitive to cold climate than leaf et al., 2014). The contrasting pattern of beech and birch in Norway might indicate that beech benefits more than birch from a prolonged period of photosynthesis, with later OFS of ca. 3.5 weeks. In fact, beech needs more return for its higher investment in thicker and longer-living leaves (Martìnez-Garza et al., 2005) and possibly a larger carbohydrate reserve to refoliate the entire canopy the next spring (Hoch et al., 2003; Michelot et al., 2012). Later OFS in beech than birch was also observed in Spain, but not in Belgium. The anomaly we observed at the beech stand in Spain (WTC at least until mid-November, with WTC still around 0.5% in early December) suggests a climatic influence due to the location in the Mediterranean region, in agreement with Sass-Klaasen & Eckstein (1995), who suggested that the formation of wood in beech during the second part of the growing season is mainly influenced by external factors. However, Sampling in spring of 2018, showed that the thickening phase of the 2017 ring was not fully completed by end March but only by early-mid May (data not shown). This indicates that the development of the latest formed cells in the previous ring was completed in the period of cambial reactivation in the following spring i.e. end of April / early May. In other words, there are no indications of a slow-growth activity in winter, which is more likely in warmer Mediterranean regions (Vieria et al., 2017; Cherubini et al., 2003), and it is

  temperature(Björklund et al. 2020) and cold spells[START_REF] Piermattei | The "blue ring'': anatomy and formation hypothesis of a new tree-ring anomaly in conifers[END_REF], significant inter-annual differences are possible, making comparisons between sites measured during different years difficult to interpret. At the same sites, CWG of Quercus robur and Quercus pyrenaica varied on average 20 days between two years (Pérez-de-Lis et al., 2017).

A c c e p t e d M a n u s c r i p t

  M. (2009) Genetic and environmental variation in spring and autumn phenology of biomass willows (Salix spp.): effects on shoot growth and nitrogen economy. Tree Physiology. 29(12):1479-1490. 49. Zuidema PA, Poulter B, Frank DC. (2018) A Wood Biology Agenda to Support Global Vegetation Modelling. Trends in Plant Science. 23(11):1006-1015.

Figure 1 :

 1 Figures references

A c c e p t e d M a n u s c r i p t 23

 23

Figure 2 :

 2 Figure 2: (a, d) Daily precipitation, (b, e) average 10-day mean temperature and (c, f) average 10-day

A c c e p t e d M a n u s c r i p t 24

 24

Figure 3 :

 3 Figure 3: Example of the breakpoint analysis to determine the day of year (DOY) when onset of foliar senescence (OFS) started according to a) change in canopy coloration and b) chlorophyll degradation for the infertile silver birch (Betula pendula) stand in Belgium in 2017.

Figure 4 :

 4 Figure 4: Example of a microscopic section prepared from a silver birch tree (Betula pendula) near the end of the growing season, with a few cells still in the wall thickening phase stained with Astra blue.

A c c e p t e d M a n u s c r i p t 25 Figure 5 :

 255 Figure 5: Percentage of stem cells in the wall thickening phenophase in late summer and autumn and

u s c r i p t 28 Figure 6 :

 286 Figure 6: Seasonal development of wood formation at (a) an infertile site of silver birch (Betula

A c c e p t e d M a n u s c r i p t 29 Figure 7 :

 297 Figure 7: Timeline of key foliar autumnal processes of the silver birch (Betula pendula), European

  al., 2009; Possen et al., 2014; Vitasse et al., 2011; Methods).

Table 6 )

 6 R I P T 11 occurred in late September for aspen in Norway (not significantly different from local birch), but only in early November for oak in Belgium, which was considerably later than for local beech and birch (sign. diff. with all other stands, p < 0.001). CWG in 2018 for the infertile birch stand in Belgium was

	at the end of September, with no significant difference compared to 2017. Similarly, no significant
	difference between years was found for CWG at the fertile beech stand in Belgium, also occurring at
	the end of September.
	Comparison OFS and CWG
	OFS occurred concurrently with CWG for all fertile birch stands. On the other hand, beech
	populations on fertile sites presented a less consistent pattern, with concurrent OFS and CWG in
	Belgium, earlier CWG than OFS in Norway (2 weeks) and no detectable CWG in Spain (Fig 6, Suppl.

. OFS and CWG were also concurrent at the infertile birch and infertile beech stand in Belgium (Fig 6, Suppl. Table 6). For the two Belgian stands (infertile birch and fertile beech) studied in both 2017 and 2018, no difference was recorded between OFS and CWG in 2017. On the other hand, the difference between OFS and CWG became significant in 2018 at both stands, with OFS taking place 4-5 weeks after CWG (Fig 6, Suppl. Table

Acknowledgements

This research was funded by the ERC Starting Grant LEAF-FALL (714916). M.C. is a Postdoctoral Fellow of the Research Foundation-Flanders (FWO). We would like to thank the institutions that gave permission to conduct research in the study areas: Agency for Nature and Forest of the Flemish Government (ANB), the Military Defense of Belgium, City of Brasschaat (Belgium), NMBU Norwegian University of Life Sciences (Norway), Natural Park and Biosphere Reserve of Montseny and Municipality of Vallfogona de Ripollès (Spain). Special thanks are given to Dirk Leyssens (ANB) and Daniel Guinart Sureda (Natural Park and Biosphere Reserve of Montseny). The meteorological data for the sites have kindly been provided by the Research Institute for Nature and Forest (INBO, Belgium) and meteorological stations at Ås (Norway), Montseny (Spain) and Sant Pau de Segúries (Spain).