Jean-Baptiste Caillau
email: jean-baptiste.caillau@univ-cotedazur.fr

Olivier Cots
email: olivier.cots@irit.fr

Pierre Martinon
email: pierre.martinon@inria.fr

ct: control toolbox -Numerical tools and examples in optimal control

Keywords: direct transcription method, shooting, bocop, nutopy, differential continuation, automatic differentiation

published or not. The documents may come

INTRODUCTION

Methods for solving optimal control problems with ordinary differential equations essentially fall into three main categories. Dynamic Programming (or Hamilton-Jacobi-Bellman method; see, e.g., ROC-HJ 1 for an implementation) computes the global optimum but may suffer from high computational costs, the so-called curse of dimensionality. Indirect methods based on Pontrjagin Maximum Principle are extremely fast and accurate but often require more work to be applied, in terms of mathematical analysis and a priori knowledge of the structure of the solution. Direct transcription methods offer a good tradeoff between robustness and accuracy and are widely used for industrial applications. There are excellent available direct solvers, usually making an intensive use of automatic differentiation (see, e.g., CasADI 2 ). For challenging problems, an effective strategy is to start with a direct method to find a first rough solution, then refine it through an indirect method. Such a combined approach is very commonly used; see, e.g., this paper [START_REF] Bonnard | Geometric and numerical methods in the contrast imaging problem in nuclear magnetic resonance[END_REF] for the optimization of contrast in medical imaging (MRI). This combination of direct and indirect methods has a lot of interest to solve optimal control problems that contain state or control constraints. In the example mentioned above, the interfacing between the two solvers Bocop 3 and HamPath 4 were done manually by ad hoc python or matlab layers. The aim of the current paper is to present a common project between several Inria and CNRS research teams from Sophia, Paris and Toulouse whose goal is to interoperate these solvers using a high level common interface. This project, coined ct: control toolbox,5 benefits from the support of Inria through a so-called Action de Développement Technologique. It is open source and welcomes external contributions.

Two examples are presented. First, a regulator example with constraints on the control; though very simple, the solution has a bang-singular structure that illustrates in an elementary fashion the interplay between direct (to capture the structure) and indirect (to compute efficiently and precisely) methods. It is solved using the new python interface of Bocop36 and the python package nutopy. 7The second example is taken from Riemannian geometry and more sophisticated. It leverages the use of differential continuation (homotopy methods) to compute conjugate and cut loci on an ellipsoid of revolution. In both examples, automatic differentiation is crucial (and used up to order three in the second one). The presented results are entirely reproducible online8 without any installation, using the Binder9 technology.

FIRST EXAMPLE: A SIMPLE REGULATOR PROBLEM

A simple regulator problem

We illustrate here the coupling of direct and indirect methods in optimal control, on a simple case of Linear Quadratic Regulator in dimension 2. More precisely, we

show how the direct method can provide the relevant information to initialize the indirect method accurately, namely with the correct control structure and with an estimate of junction times as well as state and costate values at the relevant times. We study the following problem, for which the optimal control structure consists in a bang arc followed by a singular arc.

             min 1 2 5 0 (x 2 1 (t) + x 2 2 (t)) dt ẋ1 (t) = x 2 (t), ẋ2 (t) = u(t), u(t) ∈ [-1, 1]
x(0) = (0, 1)

Direct method

We first solve the problem with the so-called direct transcription approach, that basically consists in solving the finite-dimensional nonlinear programming problem (NLP) resulting from applying a time discretization to the original optimal control problem (OCP). We reformulate the Lagrange cost as a Mayer cost by introducing a new state variable ẋ3 (t) = (x 2 1 (t) + x 2 2 (t))/2, x 3 (5) → min and use the software BOCOP for the numerical simulations. The definition of the problem consists in one C++ file for the problem functions, and one text file for all other parameters and settings. The file problem.cpp defines the final cost, dynamics and boundary conditions for the OCP. A separared definition file contains all other relevant information, such as problem dimensions, time discretization choices, bounds for the different variables and constraints, initial guess and numerical settings for the NLP solver Ipopt, etc.

We build the problem executable and solve the problem by launching the optimization.

import bocop problem_path = "." bocop.build(problem_path) bocop.run(problem_path)

Solution analysis

In addition to the state and control from the optimal trajectory, we can plot the Lagrange multipliers associated to the constraints for the discretized dynamics, that correspond to the costate variables from Pontrjagin Maximum Principle. Note that for the case of basic initial conditions x i (0) = x 0 i , the multipliers for these constraints will also match the costate variables at the initial time for the corresponding state variables. For this problem the third costate is constant and equal to -1, which is normal since the third state corresponds to the integral cost to be minimized. The solution exhibits a bang-singular structure, as can be clearly seen from Fig. 1 on the control variable in purple. The oscillations of the control over the singular arc are often encountered when using discrete transcription, however the averaged control usually corresponds to the correct singular control, meaning that the state dynamics should be close to the optimal trajectory. Here is the output of BOCOP:

Bocop returns status 0 with objective 0.377 and constraint violation 1.221e-12

Multipliers for initial conditions: [-0.94216793 -1.44190126 -1. ]

Indirect method

Using the information gathered on the solution from the direct method, we can now solve the problem accurately with an indirect method, with the knowledge of both the optimal control structure and an estimate of the junction time and state-costate values at this time. We use the package nutopy for the indirect shooting method. Applying Pontrjagin maximum principle, one has to consider the following Hamiltonian (in normal form since there are no terminal constraints):

H(x 1 , x 2 , p 1 , p 2 , u) = -(x 2 1 + x 2 2 )/2 + p 1 x 2 + p 2 u.
One checks that the control is either bang (equal to ±1) or singular. Singular arcs are actually of order one, given by the condition p 2 = 0 (singular arcs actually leave on the codimension two submanifold {p 2 = x 2p 1 = 0} of the cotangent bundle T R 2 R 2 ×R 2 , and the singular control (in feedback form) is u = x 1 . So we have a competition between three Hamiltonian flows, respectively associated with H ± := H(x, p, ±1), H s := H(x, p, x 1 ).

We define these three Hamiltonians and the corresponding flows in order to set up our shooting algorithm. (For this particular example, in view of the BOCOP run, we restrict to just -1 bang arcs.)

import nutopy as nt import numpy as np t0 = 0. tf = 5. x0 = np.array([0., 1.])

First, we define the (OCP) setting the running cost, that is the integrand of the Lagrange cost, and the dynamics.

# Running cost def df0fun(t, x, dx, u, du): df0 = x[0]*dx[0] + x[1]*dx[1] return df0 def d2f0fun(t, x, dx, d2x, u, du, d2u): d2f0 = dx[0]*d2x[0]+dx[1]*d2x[1] return d2f0
@nt.tools.

[0] = p1[1] # p2(t1) s[1] = x1[1]-p1[0] # x2(t1)-p1(t1) s[2] = pf[0] # p1(tf)

SECOND EXAMPLE: CUT AND CONJUGATE LOCI OF AN ELLIPSOID OF REVOLUTION

Preliminaries

We recall standard notions of Riemannian geometry that can be found in [START_REF] Berger | A panoramic view of Riemannian geometry[END_REF]; [START_REF] Do | Birkhäuser, Mathematics: Theory & applications[END_REF]. On a connected Riemannian manifold (M, g) of dimension n, the metric g induces a distance function d(q 0 , q 1 ) for any pair (q 0 , q 1 ) ∈ M 2 . The distance d(q 0 , q 1 ) being the infimum of lengths of continuously differentiable curves q joining the points q 0 and q 1 . If we denote by l(q) the length of a curve q, then on a connected and complete Riemannian manifold, for any pair (q 0 , q 1 ) ∈ M 2 , there exists a continuously differentiable curve q such that d(q 0 , q 1 ) = l(q). In this case, the solution curve q being necessarily chosen among the set of the so-called geodesic curves (or geodesics). The geodesics are so candidates as minimizers and we know that every geodesic is a minimizing curve at least on short distances.

Our main goal is to compute for a fixed q 0 ∈ M , the partial distance function d(q 0 , •) together with the set of minimizing geodesics, that is we want to find for any point q 1 ∈ M , the distance d(q 0 , q 1 ) and the geodesic q such that d(q 0 , q 1 ) = l(q). To do so, we compute the geodesics, which are projections of extremals given by the flow of the Hamiltonian vector field #-H associated to the Hamiltonian H given by the Legendre transform of the metric g.

The cut point of a geodesic is the point where it ceases to be minimizing. Fixing q 0 ∈ M , the set of cut points of all the geodesics starting from q 0 is the cut locus denoted Cut(q 0 ). Our main goal is thus equivalent to compute Cut(q 0 ).

We consider the academic case of an oblate ellipsoid of revolution to illustrate the numerical tools. The theoretical results related to this case may be found in [START_REF] Itoh | The cut loci and the conjugate loci on ellipsoids[END_REF]. In the Riemannian frame, the cut points may be of two kinds, either conjugate points or isochronous separating points. A conjugate point is a point where the geodesic q ceases to be optimal among the geodesics C 1close to q. The conjugate points are fold points of the geodesic flow and are due to the intrinsic curvature of the Riemannian manifold. This points may be computed by linearization of the flow of #-H, that is computing Jacobi fields. An isochronous separating point is a point where two distincts minimizing geodesics intersect, that is it is a self-intersection of a wavefront of the geodesic flow. The set of separating points is called the separating locus.

In the following, we present how to manipulate the nutopy package to compute the cut locus and related objects. See [START_REF] Bonnard | Geometric and numerical techniques to compute conjugate and cut loci on Riemannian surfaces[END_REF][START_REF] Facca | Computing the Cut Locus of a Riemannian Manifold via Optimal Transport[END_REF]; [START_REF] Itoh | Thaw: A Tool for Approximating Cut Loci on a Triangulation of a Surface[END_REF]; [START_REF] Sinclair | Loki: Software for Computing Cut Loci[END_REF] for others methods to compute cut loci in Riemannian geometry.

The Hamiltonian and flow for the oblate ellipsoid of revolution

The Cartesian equation of the normalized 2d ellipsoid of revolution may be written as

x 2 + y 2 + z 2 ε 2 = 1
The associated ellipsoid of revolution may be generated by a rotation of axis Oz of the parameterized curve y = cos u, z = ε sin u, with u ∈ [0, 2π] and where 0 < ε < 1 corresponds to the oblate (flattened) case while ε > 1 is the prolate (elongated) case. We parameterize the Cartesian coordinates by the azimuth θ ∈ [-π, π] and the latitude ϕ ∈ [-π/2, π/2], by the relations x = cos ϕ cos θ, y = cos ϕ sin θ and z = ε sin ϕ.

The restriction of the Euclidean metric of the ambient space R 3 gives the following induced metric on the ellipsoid: g = g 1 (ϕ) dθ 2 + g 2 (ϕ) dϕ 2 where g 1 (ϕ) = cos 2 ϕ and g 2 (ϕ) = sin 2 ϕ+ε 2 cos 2 ϕ. In the following, we consider the oblate case and we fix ε = 0.75 and q 0 = (0, 0). The associated Hamiltonian is

H = 1 2 p 2 θ g 1 (ϕ) + p 2 ϕ g 2 (ϕ)
that we implement in Fortran: From the Hamiltonian H, we define the Hamiltonian vector field defined on the cotangent bundle T * M (M being the ellipsoid) #-H(q, p)

def = ∂H ∂p (q, p), - ∂H ∂q (q, p) .
We finally define the Hamiltonian exponential map,

e t #- H : T * M → T * M , by e t #- H (q 0 , p 0 ) def = (q(t, q 0 , p 0 ), p(t, q 0 , p 0 )),
where the extremal z(•, q 0 , p 0 ) = (q(•, q 0 , p 0 ), p(•, q 0 , p 0 )) is defined as the maximal solution of the Cauchy problem ż(t) = #-H(z(t)), z(0) = (q 0 , p 0 ). Thus, the Hamiltonian exponential map gives us the flow of extremals, associated to the Hamiltonian system.

# Hamiltonian h = nt.ocp.Hamiltonian(hfun)

# Hamiltonian exponential map extremal = nt.ocp.Flow(h)

Geodesics

By homogeneity, the extremals may be parameterized fixing H = 1/2. This amounts to parameterize the geodesics by the arc length. Thus, given an initial point q 0 ∈ M , the initial covector p 0 must satisfy H(q 0 , p 0 ) = 1/2. One can write that p 0 ∈ S 1 = {p | ||p|| = 1}, defining the norm ||p|| 2 = 2H(q 0 , p) on T * q0 M . Hence, the initial covector may be parameterized by its angle α 0 ∈ [0, 2π). We write

p 0 = p(α 0 ) def = cos α 0 g 1 (ϕ 0 ), sin α 0 g 2 (ϕ 0 )
this parameterization, with p(α 0 ) ∈ T * q0 M . Finally, a geodesic starting from q 0 is the projection of an extremal parameterized by its initial angle α 0 and is given by the following classical exponential map:

exp q0 (t, α 0 ) def = π q e t #-
H (q 0 , p(α 0 )) ,

where π q (q, p) = q is the canonical projection on the state space.

# Initial covector def covector(q, alpha): g1, g2 = metric(q) p0 = [np.cos(alpha)*np.sqrt(g1), np.sin(alpha)*np.sqrt(g2)] return p0 # Geodesic def geodesic(t, alpha0): p0 = covector(q0, alpha0) q, p = extremal(t0, q0, p0, t) return q

Conjugate locus

A point q(t c ) on an extremal z = (q, p) is conjugate to q 0 = q(0) if there exists a Jacobi field δz = (δq, δp), solution of the linearized system along the extremal,

δ ż(t) = ∂ ∂z #- H(z(t)) • δz(t),
which is non-trivial (δq ≡ 0) and vertical at t = 0 and t c > 0 (called the conjugate time), that is δq(0) = δq(t c ) = 0. The conjugate locus is the set of such first points on extremals departing from q 0 . Conjugacy is classicaly related to local optimality of extremals in the relevant topologies. Let q : t → exp q0 (t, α 0 ) being a reference geodesic and introduce

F conj : R * + × [0, 2π) → R as F conj (t, α) def = det exp q0 (t, α) .
Then, in the Riemannian setting, q(t c ) is conjugate to q 0 if and only if F conj (t c , α 0 ) = 0. The conjugate locus from q 0 is thus given by Conj(q 0 ) = α0∈[0,2π) { exp q0 (t 1c , α 0 ) s.t. F conj (t 1c , α 0 ) = 0}, where t 1c has to be understood in the sense that it is the first conjugate time. To compute the conjugate locus, we need to compute the first conjugate times, that is to compute

α0∈[0,2π) {(t 1c , α 0 ) | F conj (t 1c , α 0 ) = 0}.
Algorithm. We have to compute a subset of F -1 conj ({0}). We use the numerical continuation method (or homotopy method) from [START_REF] Allgower | Introduction to numerical continuation methods[END_REF] implemented in the nutopy package. Under some regularity assumptions, the set F -1 conj ({0}) is a disjoint union of differential curves, each curve being called a path of zeros. To compute a path of zeros, we search a first point on the curve by fixing the homotopic parameter (here it will be α 0 ) to a certain value α * 0 and then calling a Newton method to solve F conj (•, α * 0 ) = 0. The Newton solver in the nutopy package is the hybrj code from the minpack library, see [START_REF] Moré | User Guide for MINPACK-1[END_REF].

When we have our initial point on the path of zeros, we use a Predictor-Corrector (PC) method with arc length parameterization to compute the differential curve. The nutopy package implements a PC method with a highorder Runge-Kutta scheme with adaptive step size for the prediction from Hairer et al. (1993). 10 The correction step is performed by a simplified Newton algorithm from (Hairer and Wanner, 1996, p. 119) with few iterations and where the Jacobian of the associated system to solve is not updated along the iterations.

For the prediction and the correction steps, we need to compute the Jacobian of F conj . It is computed by a combination of automatic differentiation (to get #-H and its derivatives up to order 2) performed with the tapenade software, see [START_REF] Hascoët | The Tapenade Automatic Differentiation tool: principles, model, and specification[END_REF], and variational equations to compute Jacobi fields and their derivatives with respect to the initial angle α 0 .

# Jacobi field @nt.tools.vectorize(vvars=(1,)) def jacobi(t, alpha0): p0, dp0 = covector(q0, (alpha0,1.)) (q,dq),(p,dp) = extremal(t0,q0,(p0,dp0),t) return (q, dq), (p, dp) Fig. 3. Geodesics and their envelope when the flow is folding, i.e. the conjugate locus (in red).

Cut locus

The cut locus Cut(q 0 ) is the union of the cut points of the geodesics departing from q 0 . A cut point is either a conjugate point or an isochronous separating point, see (Do Carmo, 1988, Proposition 2.2 p. 267), that is a point where two minimizing geodesics intersect with the same time (or length). A separating point (or splitting point) corresponds to a self-intersection of a wavefront. To compute the cut locus we need at this stage, to compute the separating locus (or splitting locus) and then compare for each geodesic its first conjugate time with its separating time. We introduce the following mapping:

F split (t, α 1 , q, α 2 ) def = qexp q0 (t, α 1 ), qexp q0 (t, α 2 ) .

The splitting locus is then given by solving F split = 0 since we have Split(q 0 ) ⊂ {q ∈ M | ∃ (t, α 1 , α 2 ) s.t. α 1 = α 2 and F split (t, α 1 , q, α 2 ) = 0}.

The splitting locus is also computed by homotopy.

# Equations to compute Split(q0) def Fsplit(y, alpha2):

t = y[0] alpha1 = y[1] q = y[2:4]
q1, _ = extremal(t0, q0, covector(q0, →alpha1), t) q2, _ = extremal(t0, q0, covector(q0, →alpha2), t) eq = np.zeros(4) eq[0:2] = q-q1 eq[2:4] = q-q2 return eq At the end, the cut locus of the oblate ellipsoid of revolution is the union of the splitting locus with the two conjugate points at the extremity of Split(q 0 ), that is Cut(q 0 ) = Split(q 0 ) ∪ {exp q0 (t * , 0), exp q0 (t * , π)}, where t * is the injectivity radius, the first positive time solution of F conj (t * , 0) = 0, see Fig. 4 for a view of the cut locus and see [START_REF] Itoh | The cut loci and the conjugate loci on ellipsoids[END_REF] for theoretical details.

Fig. 1 .

 1 Fig. 1. States, costates and control computed by BOCOP.The bang-singular structure of the control is clearly observed.

  arc x1, p1 = fbang(t0, x0, p0, t1) s[3:5] = x1 -x1_in # matching conditions s[5:7] = p1 -p1_in # matching conditions # singular arc xf, pf = fsing(t1, x1_in, p1_in, tf) s

#Fig. 2 .

 2 Fig. 2. States, costates and control computed by shooting (nutopy).

  subroutine hfun(x, p, e, h) double precision, intent(in) :: x(2), p(2), e double precision, intent(out) :

  0 of the costate at initial time, (x 1 , p 1 ) the state and costate values at the junction point between the (-1) bang arc and the singular one, and the time t 1 where this junction occurs. A standard solver, properly initialised thanks to the previous solution provided by BOCOP, is eventually called to solve the problem. As is clear from the obtained numerical results, a very precise solution is computed.

	# Shooting function
	def shoot(z):
	→tensorize(df0fun,d2f0fun,tvars=(2,3))
	def f0fun(t, x, u):
	f0 = 0.5 * (x[0]**2 + x[1]**2)
	return f0
	# Dynamics
	def dffun(t, x, dx, u, du):
	df = np.zeros(2)
	df[0] = dx[1]
	df[1] = du
	return df
	def d2ffun(t, x, dx, d2x, u, du, d2u):
	d2f = np.zeros(2)
	return d2f
	@nt.tools.tensorize(dffun,d2ffun,tvars=(2,3))
	def ffun(t, x, u):
	f = np.zeros(2)
	f[0] = x[1]
	f[1] = u
	return f
	o = nt.ocp.OCP(f0fun, ffun) # OCP
	Then, we define the Hamiltonians and flows. For the
	control laws we do not present the part of the code defining
	the derivatives.
	def ubang(t, x, p): # Bang control
	return -1.
	def using(t, x, p): # Singular control
	return x[0]
	# Hamiltonians and Flows
	hbang = nt.ocp.Hamiltonian.fromOCP(o, ubang)
	hsing = nt.ocp.Hamiltonian.fromOCP(o, using)

fbang = nt.ocp.Flow(hbang) fsing = nt.ocp.Flow(hsing) Assuming a bang-singular structure, it is then easy to define a shooting function whose unknowns are the value p

ct.gitlabpages.inria.fr/gallery

ct.gitlabpages.inria.fr/bocop3

ct.gitlabpages.inria.fr/nutopy

See the regulator and ellipsoid examples on the ct project gallery.

mybinder.org

The prediction part is the same as in the HamPath software, see[START_REF] Caillau | Differential continuation for regular optimal control problems[END_REF].

ACKNOWLEDGEMENTS

The authors thank the SED 11 team of Inria Sophia for their help and commitment on this project.

Supported by InriaHub /


Thomas Gasser

Philippe Ciais

Simon L Lewis

How the Glasgow Declaration on Forests can help keep alive the 1.5 °C target

Differing interpretations as to how countries

At last year's 26th UN Climate Change Conference of the Parties (COP26) in Glasgow, 141 countries committed to halt and reverse forest loss and land degradation by 2030 [START_REF]Glasgow leaders' declaration on forests and land use[END_REF]. It was part of one of several side deals designed to keep the objectives of the Paris agreement within reach. The UK government boasted that these nations had made a landmark pledge to end deforestation [START_REF] Gov | Over 100 leaders make landmark pledge to end deforestation at COP26[END_REF].

Yet, one crucial detail was left out: whether that deforestation will be gross or net. The distinction matters, because differing interpretations of how countries can "end deforestation" significantly impact future carbon dioxide emissions. Put simply, ending gross deforestation would be a major step forward for the climate. But considering only net deforestation could be anecdotal and even be detrimental to biodiversity.

Overall, human activities related to land use emit 4 gigatonnes (Gt) of carbon dioxide (CO 2 ) per year to the atmosphere (3). These are net emissions, comprising both carbon losses and carbon gains. On the loss side, the destruction of tree biomass during deforestation, for large-scale agriculture, cattle ranching, or shifting cultivation (cycles of cutting forest for agriculture, then abandoning to recover soil fertility, then returning), and the legacy emissions from harvested wood products, form the bulk of gross emissions, totaling 14 Gt CO 2 per year (3). On the gain side, secondary forest regrowth after reforestation, agricultural abandonment, or during shifting cultivation form the bulk of the gross carbon removal, totaling À10 Gt CO 2 per year (3). The 4 Gt headline emission figure for land use, then, comes from the simple imbalance between these two much larger gain and loss figures.

If countries interpret the Glasgow commitment as needing to halt gross forest area loss, in the way the UK government presented the pledge to the world, then the gross carbon source will rapidly reduce, whereas the carbon uptake from secondary forests will sustain for a few decades. Should this happen in all signatory countries, we estimate it would generate a net carbon uptake in 2030 of about 1 Gt CO 2 annually. This is a 5-Gt CO 2 saving from current policies, and it would contribute to closing the gap between the national pledges submitted under the Paris Agreement and the efforts needed to limit warming to 1.5 °C (4).

At the other extreme, if all countries aimed to halt net forest loss by expanding monoculture plantations to compensate for the areas lost to continued deforestation, the gross source will not be abated and the low-carbon plantations will deliver little removal, leading to carbon emissions in 2030 that are comparable with those of the 2010s.

It is therefore vital that countries specify how they intend to meet their commitment. To do so, they should produce separate estimates of their target areas of intact primary forests and of new secondary forests-the latter with a breakdown of areas of plantations and of naturally regenerating or restored forests. Countries could include this information in their new updated and enhanced climate pledges, known as nationally determined contributions (NDCs), to be announced before COP27 in Egypt in November 2022.

Key Contrasted Scenarios

To demonstrate the importance of the distinction between gross and net forest area loss, we developed a "business-as-usual" (BAU) scenario along with three illustrative scenarios that all fulfill the Glasgow commitment to halt forest area loss by 2030. We simulated global land use CO 2 emissions implied by these scenarios, using the compact Earth system model OSCAR and its bookkeeping module [START_REF] Gasser | The compact Earth system model OSCAR v2.2: Description and first results[END_REF][START_REF] Gasser | Historical CO2 emissions from land use and land cover change and their uncertainty[END_REF].

The BAU scenario linearly extrapolates land use activity trends from the last five years until 2030 and maintains that level afterwards. This leads to a decrease in forest area and a stabilization of land use CO 2 emissions around the last decade's level, at about 4 Gt CO 2 per year as of 2030, with slowly but steadily increasing gross emissions and gross removals (see Fig. 1).

In the "End Gross Forest Loss" scenario, gross deforestation in signatory countries is brought to a halt by 2030. This leads to a steep decrease of net emissions that switch to a net removal of CO 2 from the atmosphere slightly before 2030, reach about À2 Gt CO 2 per year after, and maintain this level of removal for about two decades. This source to sink transition is explained by the reduction in gross emissions (by about a factor two) from significantly stopping biomass destruction, whereas gross removals from forest regrowth remain mostly steady until they begin to slow as new forests grow older.

In the intermediate "End Net Forest Loss" scenario, gross deforestation in signatory countries is reduced, albeit only to the point of matching their BAU gross plantation/reforestation level. Despite reaching a net zero forest area change in 2030 similar to the previous scenario, the continuation of some deforestation (including shifting cultivation) does not reduce gross emissions enough to turn global forests into a carbon sink. In 2030, net emissions reduce by half compared with BAU to 2 Gt CO 2 per year, and they reach close to zero only in 2050.

In the "End Tree Cover Loss" scenario, the signatory countries follow their BAU gross deforestation level but increase their forest area by developing new plantations to balance these gross deforestation losses. This scenario contrasts with others because its net emissions are slightly higher than the BAU in 2030 (by about 10%), owing to increased gross emissions triggered by the conversion of land into plantations. Only later does the marginally increased gross removal catch up with gross emissions so that net emissions slowly go down, albeit not at a pace capable of compensating for much of the earlier emissions.

In summary, even though the three scenarios all appear to comply with the Glasgow commitment, they produce widely different net carbon gains. Our simulations demonstrate that the level of emissions reduction (if any) depends on whether gross or net deforestation is reduced to zero. These differences sum to globally significant amounts when considering land-use emissions through to 2050. Compared with the BAU, the End Gross Forest Loss scenario sequesters around 143 ± 38 Gt CO 2 by 2050, whereas the End Net Forest Loss scenario does half as much (68 ± 20 Gt CO 2 ), and the End Tree Cover Loss scenario produces no significant carbon sequestration (8 ± 16 Gt CO 2 ).

Should countries interpret the Glasgow pledge as a commitment to reduce net deforestation, the results could follow the unsatisfactory example of the Bonn Challenge [START_REF]The Bonn challenge[END_REF], a global call to restore 350 million hectares of forest by 2030. In this restoration challenge, 45% of the pledged area is converted to monoculture plantations [START_REF] Lewis | Restoring natural forests is the best way to remove atmospheric carbon[END_REF]. Although this is not what most people would think of as forest restoration, it does meet the UN Food and Agriculture Organization's definition of forests that encompasses pristine forests, restored forests, and monoculture plantations.

Similarly, although most would interpret the Glasgow Declaration on Forests as halting the loss of primary old-growth forests, it could be met by simply replacing high-carbon and high-biodiversity primary forests with monoculture plantations. This is why participating signatory countries must distinguish among all three types of forests-old-growth primary forests, regenerating and restored forests, and plantations-in the plans of how nations will comply with their commitment and in their progress reports. If they don't, the Glasgow Declaration risks resembling an intergovernmental greenwashing operation.

Fig. 1 .

 1 Fig. 1. Here we illustrate "business as usual" (BAU) along with three Glasgow-compliant scenarios and implied fluxes. Left panels illustrate the main land use change occurring in the scenarios. Right panels show the carbon fluxes to the atmosphere. Net emissions are made of gross emissions and gross removals. The BAU net emissions are recalled in other panels with a dashed line. Uncertainty bars are one standard deviation range. Supplementary Information includes further details of the scenario descriptions.

To whom correspondence may be addressed. Email: gasser@iiasa.ac.at. This article contains supporting information online at http://www.pnas.org/lookup/suppl/doi:10.1073/pnas.

2200519119/-/DCSupplemental. Published June 2, 2022. PNAS 2022 Vol. 119 No. 23 e2200519119 https://doi.org/10.1073/pnas.2200519119 1 of 3 OPINION Downloaded from https://www.pnas.org by 157.136.75.78 on June 24, 2022 from IP address 157.136.75.78.

of

https://doi.org/10.1073/pnas.2200519119 pnas.org Downloaded from https://www.pnas.org by 157.136.75.78 on June 24, 2022 from IP address 157.136.75.78.

ACKNOWLEDGMENTS. We thank R. A. Houghton and A. D. de Almeida Castanho for sharing their FRA2020-based land use and land cover change drivers. T.G. acknowledges support from the EU H2020 "CONSTRAIN" project (Grant No. 820829) and the FWF "ERM" project (grant P31796-N29).

Data and materials availability. The source code of OSCAR v3.1.2 is available at https://github.com/tgasser/OSCAR. The original GCB2021 outputs are available at https://www.globalcarbonproject.org/carbonbudget. Additional data and outputs are available on request.


Keywords: Zooplankton, Volume, Biomass, Copepod, Geometrical model, in situ imaging. 1

Accurate plankton biomass estimations are essential to study marine ecological processes and biogeochemical cycles. This is particularly true for copepods, which dominate mesozooplankton. Such estimations can efficiently be computed from organism volume estimated from images. However, imaging devices only provide 2D projections of 3D objects. The classical procedures to retrieve volumes, based on the Equivalent Spherical Diameter (ESD) or the best-fitting ellipse, are biased. Here, we present a method to correct these biases. First a new method aims to measure body area and fit an ellipse. Then, the body of copepods is modeled as an ellipsoid whose 2D silhouette is mathematically derived. Samples of copepod bodies are simulated with realistic shapes/sizes and random orientations. Their total volume is estimated from their silhouettes using the two classical methods and a correction factor is computed, relative to the known, total, volume. On real data, individual orientations and volumes are unknown but the correction factor still holds for the total volume of a large number of organisms. The correction is around -20% for the ESD method and +10% for the ellipse method. When applied to a database of ∼150,000 images of copepods captured by the Underwater Vision Profiler, the corrections decreased the gap between the two methods by a factor of 54. Additionally, the same procedure is used to evaluate the consequence of the bias in the estimation of individual volumes on the slopes of Normalised Biovolume Size Spectra and show that they are, fortunately, not sensitive to the bias. S1.6

Introduction

Plankton forms an extremely diverse community [START_REF] De Vargas | Ocean plankton. eukaryotic plankton diversity in the sunlit ocean[END_REF] and its members are keystone components of Earth's biosphere. First, photosynthetic plankton is responsible for about half of the fixation of carbon dioxyde from the atmosphere and therefore produces an equally large amount of dioxygen [START_REF] Behrenfeld | Biospheric primary production during an enso transition[END_REF]. The death and excretion of planktonic organisms result in massive amounts of carbon being sequestered to the seabed in the form of "marine snow"; this "biological carbon pump" is therefore an important contributor to the regulation of climate [START_REF] Volk | Ocean carbon pumps: Analysis of relative strengths and efficiencies in ocean-driven atmospheric co2 changes[END_REF]. Plankton is also a critical component of many marine food webs: it directly supports some of the largest fisheries on earth, off the coast of Chile for example [START_REF] Thiel | The humboldt current system of northern and central chile[END_REF], and some emblematic species such as corals. Finally, because plankton simply drifts, it cannot escape the conditions of the water mass it is embed-ded in. This makes planktonic organisms very sensitive to environmental change. Therefore, the contribution of plankton to the important processes described above will be influenced by the changes in Earth's climate [START_REF] Hays | Climate change and marine plankton[END_REF].

To quantify these roles of plankton and their possible change in time, accurate estimations of its biomass are needed. Major efforts have been made in that direction, to aggregate estimates from various techniques in panoceanic databases [START_REF] Buitenhuis | Maredat: towards a world atlas of marine ecosystem data[END_REF]Moriarty and O'Brien, 2013). These assessments are unfortunately fraught with uncertainties due, for example, to the difficulty of sampling plankton in open and deep seas, resulting in an uneven coverage, and to the differences between sampling techniques. Recently, the use of quantitative imaging instruments, i.e. automated cameras that take pictures at high frequency in a controlled manner that allows the computation of concentrations, have been highlighted as a one promising avenue for the estimation of plankton biomass and biodiversity at global scale [START_REF] Lombard | Globally consistent quantitative observations of planktonic ecosystems[END_REF]. However, most of these instruments take two-dimensional pictures, which reduce the organism to its projection on the imaged plane. The volume of organisms is estimated from this imperfect representation, mostly using two classic approaches detailed below (a spherical equivalent or an ellipsoidal approximation) and then converted to carbon biomass through a factor (their carbon density). For any organism that is not perfectly spherical, both approaches induce an error in the volume and therefore the biomass, that has been neglected so far.

Several of these quantitative imaging instruments target mesozooplankton (organisms between ∼ 200 µm and ∼ 2 mm in size), which are major contributors to biogeochemical fluxes in the ocean [START_REF] Buitenhuis | Biogeochemical fluxes through mesozooplankton[END_REF]. Within mesozooplankton, copepods represent about 85% of organisms [START_REF] Longhurst | Chapter 9 -the atlantic ocean[END_REF] and contribute largely to the carbon flux by grazing on smaller plankton, producing fecal pellets and dead carcasses that sink. The intensity of that flux is proportional to the biomass of organisms, particularly those present in the surface, mixed layer [START_REF] Buitenhuis | Biogeochemical fluxes through mesozooplankton[END_REF]. Because copepods are not spherical but rather ellipsoidal (Fig. 1), the estimation of their volume by imaging instruments is affected by the aforementioned projection error.

The purpose of the present work is to quantify the error in the estimation of the total volume of copepods from a collection of two dimensional images and to propose a correction of this estimation. This error can potentially change current estimations of carbon biomass in those organisms significantly, and its relationship with the carbon flux they generate. Both are important parameters of global biogeochemical models (Buitenhuis et al., (a) (b)

Figure 1: Examples images of copepods and theirs ellipse fits: (a) top or bottom view, (b) side view. On these projections, the bodies of copepods can be well approximated by ellipses, which we assume to be ellipsoids in 3 dimensions (see Appendix S1). Ellipse fits on copepod images are obtained by ZooProcess [START_REF] Gorsky | Digital zooplankton image analysis using the zooscan integrated system[END_REF] in red, affected by the antennae, and our method in green, that fits the prosome of the copepod better.

2006). We start by reviewing the two standard methods to estimate the volume of planktonic objects from their two-dimensional projection and detail the magnitude of errors they induce. We then present our correction method, followed by an application to a global data set of copepod images collected by the Underwater Vision Profiler [START_REF] Picheral | The underwater vision profiler 5: An advanced instrument for high spatial resolution studies of particle size spectra and zooplankton[END_REF]. Finally we conclude and present perspectives for this work.

Standard volume estimations

To safely infer a volume in three dimensions from a shape projected in two dimensions, the following assumptions are made:

• The distance between imaged objects and the camera is the same for all objects (or differences are negligible).

• One of the two following statements is true.

-The size of the object is negligible compared to its distance from the camera. Hence, even if the camera has a perspective acquisition geometry, it can be approximated well enough by a parallel one.

-The acquisition system follows a line scanner principle (then, its acquisition geometry is intrinsically parallel);

With these hypotheses, the imaging process can be schematically represented as in Fig. 2. All in situ plankton imagers presented in [START_REF] Lombard | Globally consistent quantitative observations of planktonic ecosystems[END_REF] satisfy these conditions. Then, an hypothesis on the three dimensional shape of the object has to be made. Two standard approaches exist: a spherical equivalent and an ellipsoidal hypothesis.

Using a spherical equivalent (M ESD )

There is a unique disk with the same area A as the organism's projected silhouette (the silhouette of the organism as observed on the image), and its diameter is ESD = 2 A π . The Equivalent Spherical Diameter estimation method (M ESD ) makes the assumption that the volume of the organism can be approximated by the volume of the sphere of diameter ESD, that is

V ESD = 4 3 π ESD 2 3 . (1)
C a m e r a

Figure 2: Representation of the geometrical setup of the imaging, with some notations: r 1 , r 2 , and r 3 the true semi-axes of an ellipsoidal object, with r 1 ≥ r 2 ≥ r 3 by convention; ρ 1 , ρ 2 the semi-axes of the projected ellipse, with ρ 1 ≥ ρ 2 by convention; A the area of the projected shape.

If the organisms were indeed spherical and the 3-D-to-2-D acquisition system performs a parallel projection, then V ESD would be the exact volume.

For ellipsoidal objects, like copepods, the projection silhouette is an ellipse of semi-axes ρ 1 and ρ 2 , with ρ 1 ≥ ρ 2 by convention (see Fig. 2). Its area is equal to πρ 1 ρ 2 . Therefore, the equivalent diameter is ESD = 2 √ ρ 1 ρ 2 . In that case, the error of the M ESD estimation can range from large underestimation to even larger overestimation, when the copepod is elongated in the direction of the optical axis (Fig. 3).

Using a best-fitting ellipse (M ELL )

A common alternative to M ESD is to fit an ellipse shape on the projection and construct an ellipsoid in three dimensions (M ELL ). It should be more appropriate for objects of ellipsoidal shape, such as copepods (assuming the antennae and urosome are thin/small enough for their influence on the volume to be negligible). It proceeds as follows: (i) an ellipse is fitted on the object silhouette, defining two semi-axes: ρ 1 and ρ 2 (Fig. 2), (ii) the smallest semi-axis of the fitted ellipse (ρ 2 ) is duplicated to form the triplet of semi-axes of an ellipsoid, (iii) the volume is computed as

V ELL = 4 3 πρ 1 ρ 2 2 .
(2)

Despite the fact that the estimation of the silhouette shape is more ap- propriate than with M ESD , it still leads to errors due to the projection from 3D to 2D (Fig. 3). We can remark that, within this ellipsoid model framework, V ELL is always lower than or equal to V ESD . Indeed, we assume ρ 1 ≥ ρ 2 , therefore

√ ρ 1 ρ 2 ≥ ρ 2 ⇔ √ ρ 1 ρ 2 3 ≥ ρ 1 ρ 2 2 ⇔V ESD ≥ V ELL
with equality when the projection silhouette is a circle.

Proposed method

Instead of proposing a novel volume estimation method, our approach was to study the errors made by the standard methods M ESD or M ELL in order to propose a procedure to compensate for these errors. Thus, the figures of past studies could be re-interpreted in light of the proposed corrections and marine ecologists could apply these corrections to future studies, sticking to their standard estimation method of choice. However, we noticed that copepod antennae, when visible, can affect, sometimes dramatically, the measurement of the projected area of their prosome (which constitutes the bulk of their volume) or the fitting of an ellipse to their silhouette (Fig. 1). Therefore, as a preliminary step, we propose area estimation and ellipse fitting approaches tailored to copepods, to eliminate that source of error. The area estimation is based on a form of mathematical morphology opening that first performs an erosion to discard the antennae and then a dilation to recover the area of the copepod prosome. The proposed ellipse fitting method relies on the same operations for discarding the antennae and ensuring a better fit of the ellipse (see Appendix S3 for details).

Principle for correction of total volume

To derive the correction to apply to a volume estimated with M ESD or M ELL , we model the body of a copepod by an ellipsoid (i.e. in three dimensions), study how an ellipsoid projects onto a plane, estimate the volume from its projection using M ESD and M ELL , and compare it with the known, groundtruth volume of the ellipsoid model. This process is performed on many ellipsoids, projected at many angles, and the correction is computed on the overall, total, volume because it is impossible to retrieve the projection angle and 3-D shape from only the projected silhouette of an object in a single imaging plane. This process amounts to (i) neglecting the volume of the antennae and urosome, which is an assumption implicitly made by the existing methods; (ii) not considering acquisition phenomena such as optical response and electronic noise, hence reducing the imaging process to the 3-D to 2-D projection onto the camera sensor plane, because those processes are not general (i.e. camera-dependent) and most often not characterised for plankton cameras; (iii) considering the acquisition projection parallel, which is most often the case in our applicative context (see section Section "Standard volume estimations") and simplifies the theoretical study (some supplementary materials still considers the more general perspective case, before applying it to the parallel case).

Silhouette of a projected ellipsoid

An ellipsoid centered on the origin is composed of the ensemble of 3-D points x verifying

x M x = 1

where M is a real, symmetric, positive definite, 3×3-matrix whose elements are denoted by m ij . The volume of the ellipsoid is defined by

V = 4 3 π det(M ) . ( 4
)
Matrix M encodes the overall size, shape (semi-axes ratios) and orientation of the ellipsoid. It can be written using a block matrix notation

M = M 11 M 21 M 21 m 33 (5)
where m 33 is a scalar (the dimensions of the other terms follow). If the ellipsoid is aligned on the axes of the coordinate system, then its form is

M =   1/r 2 1 0 0 0 1/r 2 2 0 0 0 1/r 2 3   (6)
where the r i 's are the semi-axes.

In Appendix S1, it is shown that the silhouette of the parallel projection of an ellipsoid is an ellipse, the two semi-axes of which are defined by

ρ i = m 33 λ i , i ∈ {1, 2} (7)
where

λ i = tr(P ) + σ i √ ∆ /2, (8)
P = m 33 M 11 -M 21 M 21 , (9)
|σ i | = 1 and σ 1 σ 2 = -1, (10)
∆ = tr(P ) 2 -4 det(P ) (11)
where tr(P ) is the trace of P , det(P ) is its determinant, and the σ i 's are chosen so that ρ 1 ≥ ρ 2 .

From these developments, we therefore have analytical (i.e. exact) definitions of an ellipsoid, its volume, its projection as an ellipse, the semi-axes and area of this projected ellipse.

Simulation of copepod bodies

By generating random ellipsoids that realistically represent copepod body shapes, observations of their projections on images can be simulated exactly, since the ρ i 's are functions of M (i.e. the semi-axes of the observed elliptical silhouette depends on the ellipsoid size, shape, and orientation).

The simplest choice to generate ellipsoids would be to draw the semiaxes from independent uniform distributions within appropriate ranges -valid ranges can be found in the literature (e.g. [START_REF] Conway | Identification of the copepodite developmental stages of twenty-six north atlantic copepods. occasional publication of the ma-rine biological association no. 21[END_REF]-and the orientation angle from an unconstrained uniform distribution. However, those parameters should ideally be adapted to the data set at hand, since copepod sizes, shapes, and acquisition viewpoints (or orientations) are related to environmental and imaging conditions. In the following, we propose a way to generate ellipsoids that follow a realistic, parametric copepod body model.

To generate random ellipsoids, one can directly generate random matrices M as defined in section Section "Silhouette of a projected ellipsoid". Alternatively, random axis-aligned ellipsoids can be generated, which are then randomly rotated. The advantage of this later procedure is that constraints are more easily imposed to create copepod bodies with realistic proportions.

Following Eq. ( 6), axis-aligned ellipsoids are defined by the three semiaxes r 1 , r 2 , and r 3 . Since ellipsoids will be randomly rotated next, these values can be chosen such that r 1 ≥ r 2 ≥ r 3 without loss of generality.

To generate random values of r 1 , r 2 , and r 3 , we need to define one Probability Density Function (PDF, or "distribution") per r i . These PDFs must be defined in accordance with the reality of copepod body shapes, either generically (e.g., according to the literature) or from the data at hand. Generating an ellipsoid then amounts to drawing one random semi-axis per PDF. If the semi-axes respect the ordering condition, then the ellipsoid is validated; otherwise, it is discarded, and a new round of random drawing must be performed. As such, the random process is of dimension three (r 1 , r 2 , r 3 ) with two conditions (r 1 ≥ r 2 and r 2 ≥ r 3 ).

Fortunately, the process can be simplified by noting that the error made by the M ESD or the M ELL method on the total volume estimation is (almost) invariant to scaling (see Appendix S2). In other words, the error computed from N ellipsoids defined by (r 1,n , r 2,n , r 3,n ), n ∈ [1..N ] is (almost) equal to the error computed by the same ellipsoids each scaled by a constant α n strictly positive (i.e. defined by (r 1,n /α n , r 2,n /α n , r 3,n /α n )). Choosing α n = r 1,n amounts to normalizing the ellipsoids so that their largest semiaxis is equal to one. Thus, the random process becomes two-dimensional (defined only by the axes ratios r 2 /r 1 and r 3 /r 1 ) with only one condition (r 2 /r 1 ≥ r 3 /r 1 ); the three, per-axis PDFs are replaced by two axes-ratio PDFs. This has two nice consequences: a statistical one and a practical one. Statistically speaking, to describe a random process through simulation, one needs "exponentially" more samples as the dimension increases (this is known as the curse of dimensionality). The number of samples, N , is limited by computational constraints; thus, reducing the dimension provides a higher quality description for the same N . Practically speaking, the shift from drawing semi-axes to drawing semi-axes ratios means that the proposed method only depends on the shape of copepods (prosome height over prosome length and prosome width over prosome length), not on their overall size, which can be considered more general (size varies across regions) and more stable.

The remaining question is how to define the PDFs of the semi-axes ratios r 2 /r 1 and r 3 /r 1 ? As for the per-axis PDFs, two reasonable options are literature-based and data-based. The literature may provide enough details to choose a PDF family (e.g., Gaussian) and set the parameters for each ratio (e.g., mean and variance for Gaussian). Alternatively, the ratios can be measured on physical samples or on images in which the copepods are seen from side (r 2 /r 1 ratio) and from the top or bottom (r 3 /r 1 ratio). Then, the required PDFs can be fitted on these measurements in a parametric (e.g., Gaussian, Beta, Gaussian mixture [START_REF] Redner | Mixture densities, maximum likelihood and the em algorithm[END_REF]) or non-parametric way (e.g., Kernel Density Estimation, KDE [START_REF] Parzen | On estimation of a probability density function and mode[END_REF][START_REF] Scott | On optimal and data-based histograms[END_REF]).

Finally, these ellipsoids, generated to match copepod body shapes, must be rotated to simulate a random acquisition viewpoint. In the absence of a strong a priori on the orientation of copepods relative to the camera, these rotations can simply be uniformly random, in all directions. But the procedure can easily be adapted to generate rotations favoring preferred orientations.

Corrected total volume

Once a set of N random ellipsoids with realistic proportions and various orientations is generated, their projection silhouettes are computed following Eq. ( 7), their volumes are estimated using M ESD (from the area of the silhouettes) and M ELL (from the semi-axes of the silhouettes), and the error between the total, true, volume of all ellipsoids and the sum of the estimated volumes is computed as

T * = N n=1 V n * N n=1 V n = W * W (12)
where the V n s are the true volumes of the generated ellipsoids and the V n * s are the corresponding estimated volumes by method " * " (ESD or ELL). Therefore, W is the true total volume of all ellipsoids in the simulation and W * is the estimated total volume.

Once T * has been estimated from simulated ellipsoids, it can be used to correct the total volume estimated from P actual images of copepods as follows

Ŵ * = P p=1 Ṽ p * T * = W * T * (13)
where Ṽ p * is the set of volumes estimated from the acquired images by method " * ', W * is their total, and Ŵ * is the corrected total estimated volume.

Note that the proposed correction formula in Eq. ( 13) provides no objective element to prefer the M ESD method over the M ELL method, or vice versa. Indeed, the respective correction factors theoretically allow to perfectly retrieve the true total volume. In practice though, the M ESD method might be a better option since the area measurement it relies on is more robust (i.e., less sensitive to acquisition noise and grayscale variations) than the ellipse fit performed in the M ELL method.

An algorithmic description of the proposed method is given in Appendix S4.

Extension to size spectra

Plankton images are often used to estimate the size of organisms and, in particular, to compute Normalised Biomass (or Biovolume) Size Spectra (NBSS). Indeed, the slope of these spectra is a proxy for the efficacy of the energy transfer from small to large organisms within an ecosystem [START_REF] Sprules | Surfing the biomass size spectrum: some remarks on history, theory, and application[END_REF]. An NBSS obviously depends on the measurement of the size of the organisms it encompasses and, therefore, its slope may be affected by the error in the size estimation due to the 3-D to 2-D projection. We will test this through a simulation procedure, similar to the one described above, and compare the slope value computed from the volumes of randomly generated ellipsoids (considered as the true value) to that derived from estimations of volume from their projections using M ESD and M ELL .

Since the NBSS depends on absolute sizes, the simulator based on standardised ellipsoids, defined by r 1 = 1 and two semi-axes ratios, must be reverted to use directly r 1 , r 2 , and r 3 , although this is a less favorable statistical context (see section Section "Simulation of copepod bodies"). As mentioned previously, the distributions of r 1 , r 2 , and r 3 can be defined from the literature or estimated from measurements. Then the constructed ellipsoids must be rotated, either uniformly or along a preferred orientation.

Experimental results

Data set

We applied the proposed method to a collection of 158,487 copepod images from the Underwater Vision Profiler 5 (UVP5) [START_REF] Picheral | The underwater vision profiler 5: An advanced instrument for high spatial resolution studies of particle size spectra and zooplankton[END_REF]. Images are grayscale, with a pixel size varying from 0.086 to 0.174 mm depending on the generation and configuration of the UVP5. These pixel sizes are used to rescale all measurements to millimeters before computation. For each image, the area measurement and the ellipse fit had been performed with the original image analysis software (ZooProcess). They were carried out again following the new image processing method described in Appendix S3. When the newly estimated volume from either M ESD or M ELL was below 0.1 mm 3 , the copepod silhouette detection was considered erroneous and the image was excluded from the data set. Images of partially cropped copepods (i.e. with a part of the body outside of the image) were also excluded, since they would bias the estimation of the area and the semi-axes. Around 2,500 images were excluded (<2% of the data set), leaving 155,936 copepod images for the following.

Total volume

As explained above, simulating copepod bodies for a correction of their overall volume requires PDFs for the two semi-axes ratios r 2 /r 1 and r 3 /r 1 and a rotation angle. We assumed uniform random orientation of copepods in the water column, since the data set is large, covers different depths and locations on the globe, and we do not know of any justification for a preferred orientation relative to the camera. The procedure for the axesratio is described below.

Estimating shape parameters

The axes ratios were measured on 295 images in which copepods were seen from the side and 265 images in which they were seen from the top or bottom. To build these samples, operators manually selected images in which the orientation of the copepod was clear, while constraining the selection to obtain distributions in latitude similar to those of the whole data set and distributions in length similar between the side and top/bottom samples (see Appendix S5). The constraint on latitude was meant to avoid biasing the samples towards a particular environment, since copepod morphology varies latitudinally. The constraint in prosome length should ideally have been checked against the whole data set, to avoid estimating the axes ratios on biases samples of the ∼150k images. However, the true prosome length, or 2 × r 1 in the 3-D ellipsoid, is unknown; only 2 × ρ 1 , the major axis of the projected ellipse, can be estimated. While r 1 ρ 1 in the side and top/bottom views, r 1 > ρ 1 in any other view, so only the distribution of ρ 1 in the side and top/bottom views can be compared.

The PDFs for r 2 /r 1 and r 3 /r 1 were then estimated from the measurements using a KDE with a Gaussian kernel of optimal width [START_REF] Scott | On optimal and data-based histograms[END_REF] (see Fig. 4).

Computing the correction factor

Once the axes ratio PDFs were fitted to the manually selected data, we used the simulation procedure presented above to compute the error in the estimation of volume T * , for both methods. We generated over 10 8 ellipsoids in order to cover the shape and orientation parameters space with a high enough resolution and, using Eq. ( 12), we obtained T ESD = 122% and Let us note that, as expected, T ESD is higher than T ELL , since we showed that the M ESD volume is always greater than or equal to the M ELL volume. Also, the correction factors are then simply computed as 1/T ESD and 1/T ELL .

T ELL = (a) (b) r 2 /r 1 r 3 /r 1
These correction factors were obtained under the assumption that they were invariant to ellipsoid scaling (Appendix S2). To verify that this is reasonable, we performed a similar simulation in the (less convenient) three parameters space (r 1 , r 2 , r 3 ), also using 10 8 ellipsoids. For that purpose, the distributions of r 1 , r 2 , and r 3 were needed. The generation of these distributions is detailed in the section regarding NBSS, below. In this framework, the errors in volume estimation were T 3D ESD = 129%, so +6% compared to T ESD , and T 3D ELL = 84%, -3% compared to T ELL . Part of this discrepancy may come from a coarser coverage of the parameters space in this three parameters case. It can be concluded that, indeed, the assumption of invariance to ellipsoid scaling made in the proposed method is reasonable.

Estimating and correcting the total volume

We estimated the total volume of the ∼150k copepods in our data set using various approaches: the M ESD or M ELL method, as computed originally by ZooProcess (ZP ) or using our improved method (Appendix S3), and with or without the corrections using the factors defined above. This is a total of eight estimations, presented in Table 1 and Fig. 5.

In theory (i.e. copepods are ellipsoids, the imaging system has perfect lenses and infinite resolution, surface measurements and ellipse fits are exact, the simulation parameters match the reality, an infinite number of samples are generated), we should obtain the same total volume estimations from M ESD and M ELL after correction. Naturally, this is not the case in practice but we can assess the effectiveness of the proposed volume correction by checking how the discrepancy between the M ESD or M ELL estimations decreases after applying the correction. This gap is divided by 4 when using the ZooProcess measurements and by 54 when using the improved versions. The fact that the corrected volumes ( ŴESD and ŴELL ) seem to converge is no proof that either one is the truth, but it at least suggests that the proposed correction method brings a significant improvement. The impact of the improved area and ellipse measurement can also be gauged in the same way. The discrepancy is divided by 2 when comparing W * ZP and W * , and by 22 when comparing the corrected versions, Ŵ * ZP and Ŵ * . This illustrates that some accuracy can also be gained in these image processing steps. Overall, if we compare the current state of the art (uncorrected total volumes obtained using classical image processing in ZooProcess) and the corrected total volumes obtained using improved image processing, the discrepancy is divided by 91, bringing the M ESD and M ELL estimations very close to each other.

Normalized Biovolume Size Spectrum (NBSS)

(a) (b) r 2 [mm] r 3 [mm] r 1 [mm] r 1 [mm]

Estimating shape parameters

As explained above, NBSS relies on absolute volumes and the distributions of the three semi-axes (r 1 , r 2 , and r 3 ) have to be estimated. They could have been estimated directly on the 295 + 265 samples used to estimate the distributions of ratios above. However, the distribution of r 1 for the 295 copepods viewed from the side shows a bias towards larger sizes, because it is difficult to tell whether a copepod is indeed viewed from its side when it is small (see Appendix S5). This bias has limited influence on the estimation of the distributions of the ratios, r 2 /r 1 and r 3 /r 1 , since the both relationships r 1 vs. r 2 and r 1 vs. r 3 are fairly linear (Fig. 6). Therefore, it was safely dismissed for the correction of the total volume, above, which relied solely on those ratios. This is not possible anymore. Instead, we build the distribution of r 1 from the ∼150k measurements of ρ 1 , its projection in the 2D plane, which can be seen as an underestimated version of r 1 . To correct for the underestimation, we multiply the values of ρ 1 by a constant, C r 1.12. This value is computed from the median r 1 /ρ 1 ratios observed in the previous simulation of the 10 8 ellipsoids projected as ellipses.

For a smoother distribution of r 1 , a fourth-order polynomial is fitted on the log-distribution in the range ∼ [1.5, 5.5] mm, and extrapolated out- side of this interval (Fig. 7). The lower limit (1.5 mm) corresponds to the start of the continuously decreasing section; before this, copepods are difficult to identify on UVP5 images and are undersampled. The higher limit (5.5 mm) is the point beyond which the density of samples is too low, making the distribution very noisy.

r 1 [mm]
The values of r 2 and r 3 are then deduced from the value of r 1 , following Fig. 6, were it can be seen that (i) r 2 and r 3 increase linearly with r 1 , (ii) the variance around the linear fit also increases with r 1 .

Finally, the values of r i s are generated as follows: a value of r 1 is drawn from its estimated distribution (Fig. 7), then values of r 2 and r 3 are drawn from a normal distribution with mean (µ) and standard deviation (σ) computed from r 1 : µ 2 = 0.3 × r 1 + 0.18, µ 3 = 0.21 × r 1 + 0.26, σ = 0.05 × r 1 (from the fits on Fig. 6). Only samples respecting the condition r 1 ≥ r 2 ≥ r 3 are kept.

Comparison of NBSS

Following the procedure described above, a set of 10 8 realistic ellipsoids are generated, with r 1 ∈ [1,6] mm. The actual ellipsoid volumes are used to compute the simulated ground-truth NBSS. Each ellipsoid is also projected, its volume is estimated with both M ESD and M ELL , and the corresponding NBSS are computed. In all three cases, the bin size is 0.1 mm 3 in log space. To assess whether the volume computation method influences the estimation of energetic transfer in the ecosystem, the slopes of a linear fit to each NBSS in the interval [1.7, 12.2] mm 3 (or [1.5, 2.9] mm) are computed (Fig. 8).

The ground-truth NBSS is between the M ELL and M ESD ones (in the linear part at least). More importantly, the three slopes are very close to each other: -1.2 for the ground-truth vs. -1.1 for M ESD and -1.2 for M ELL . To compare this with the range of natural variability in the data set, we com-puted the NBSS from images collected in polar (absolute value of latitude in [60 • ,90 • ]) and temperate (absolute value of the latitude in [20 • ,40 • ]) regions, between 0 and 150m depth. The slopes of these NBSS were -0.7 (polar) and -1.4 (temperate) with M ESD and -0.7 and -1.6 with M ELL . The amplitude of natural variability is therefore much larger than the variability induced by the volume estimation method. Therefore, despite the errors the estimation methods induce on individual volume (Fig. 3), both M ESD and M ELL seem to be valid approaches to compute NBSS and to infer the energetic transfer efficiency through a linear fit.

Discussion

We proposed a method to improve the standard estimations of the volume of copepods on images, which are based on a spherical equivalent or an ellipse fit. Our method relies on (i) more accurate measurements made on the projected silhouette and (ii) a correction factor computed from the error on the total volume due to the projection from 3D to 2D.

For the measurements, we proposed a method to estimate the projected area or to fit an ellipse on the copepod body only (Appendix S3), to avoid large biases caused by the antennae (and, more rarely, the urosome). Neglecting the volume of the antennae and urosome compared to that of the prosome seems appropriate in first order, and it is essentially what the standard ellipse fit does when it is not affected by antennae. Still, the validity of this assumption would need to be tested. The volumes of these different parts seem difficult to measure experimentally but could be assessed from detailed 3D scans of individuals, which we now have the technology for. A second aspect to test is whether our method indeed fits the projected silhouette better than the standard one. We noticed its does, on many images similar to Fig. 1 (see Appendix S3). The fact that we observed a significant reduction of the discrepancy between the total volume estimated with M ESD and M ELL when using this approach compared to the classic one also suggests a gain in accuracy. However, to assess its absolute performance, a ground-truth segmentation should be performed on a large number of images, by having human operators paint over the pixel of the image that belong to the prosome of the copepod. Then, a pixel-level match between this ground-truth and the two automated approaches (the classic one and ours) could be computed. This extremely labor-intensive effort was considered to be out of the scope of this study, especially since the proposed method is already a clear improvement over the classic one in the worst cases (Fig. 1 and Appendix S3). The core of this study is the correction of the error due to the 3D to 2D projection on the estimation of the total volume of copepods, and its application to the UVP5 data set. A potential weakness of this application is that the distributions of semi-axes ratios were estimated from a relatively small number of images (<300 for each). Since identifying copepods in a given orientation is very time consuming (and somehow subjective), a useful alternative would be to train a classifier that could automatically identify copepods seen from given angle within the total data set. Copepods were isolated from other organisms through a combination of machine learning and human classification, through the EcoTaxa application [START_REF] Picheral | Ecotaxa, a tool for the taxonomic classification of images[END_REF]. The same tools were used to train a custom classifier for side vs. top/bottom vs. other angle copepods. While it accelerated the collection of the examples in the sample sets, it did not achieve great accuracy, largely because of the overwhelming dominance of copepods seen from "other angles".

The other assumption was in the the choice of an uniform orientation distribution. While copepods in a given environment may orient themselves in a particular manner, vertically towards the surface for example [START_REF] Benfield | Estimating the in-situ orientation of calanus finmarchicus on georges bank using the video plankton recorder[END_REF], information on such behaviour is very scarce. It is, however, very likely to change with location, depth, time of day, organism age, condition, etc. Since our data set contains >150k organisms of various sizes, from different locations and times, it is currently not possible to assume an overall distribution different from a uniform one. In situ imaging instruments that do not disturb the water they image, particularly those that image the organisms from the side or in three dimensions (e.g. through holography [START_REF] Katz | Submersible holocamera for detection of particle characteristics and motions in the ocean[END_REF]), could yield more information on this orientation behaviour. From such data, a data-driven estimation of the orientation distribution could be performed. From 3D imagery in particular, or from realistic 3D models of copepods, 2D views can be generated and an image-to-orientation predictor can be trained. With such data, a last alternative could be to use an Augmented Auto-Encoder Convolutional Neural Network [START_REF] Sundermeyer | Implicit 3d orientation learning for 6d object detection from rgb images[END_REF] to generate a view of the organism with consistent orientation from inputs taken from various angles, hence isolating the orientation information in the innermost layer and allowing consistent measurements on the output image.

Finally, when working on either the semi-axes ratios (for the total volume) or the semi-axes themselves (for the NBSS), we did not take into account the joint links between the three parameters r 1 , r 2 , and r 3 . Indeed, the data at hand only allow to study the links between r 1 and r 2 (side views), and between r 1 and r 3 (top/bottom views). We did not have access to the triplet (r 1 , r 2 , r 3 ) for any given copepod. For this purpose, again, 3D im-agery from a multi-camera system is necessary.

Overall, we showed that estimating the volume of copepods with an Equivalent Spherical Diameter assumption leads to overestimating their total volume by about 20%, while using a ellipse fit on the projection underestimates it by about 10%. Given that mesozooplankton represent a carbon biomass of 0.19 PgC in the first 200 m globally and that copepods dominate this biomass (Moriarty and O'Brien, 2013), such variations would be far from negligible for global biomass estimates from images. On a data set of 150k copepod images from the UVP5, we showed that the proposed correction method can drastically reduce the discrepancy between the two standard estimations, especially when used in conjunction with the improved area and ellipse measurements. Despite these discrepancies in the total volume, a similar simulation approach showed that the same methods yield accurate NBSS slopes, likely because the cumulated error in each size bin is similar, leading to the same slope across bins. Both methods are therefore acceptable to estimate energy transfer from size spectra measured on 2D images.

S1 Projection of an ellipsoid

S1.1 Geometrical setup

A centered ellipsoid is defined by all 3-dimensional vectors verifying

x M x = 1 ( 14
)
where M is a positive definite1 3 × 3-matrix whose elements are denoted by m ij . Let (i, j, k) denote an orthonormal basis, and let O denote the origin. To study how this ellipsoid projects onto a plane using perspective projection, let us define

(i) an optical center e e =   0 0 -   (15)
where > 0 is such that e is outside the ellipsoid, and (ii) a projection plane Π described by its normal

n =   0 0 1   (16)
and its distance δ, > δ > 0, to the origin such that Π does not intersect the ellipsoid. The plane Π is equipped with the orthonormal basis (u, v) where u and v correspond to i and j respectively. Its origin O Π is located at the intersection between Π and the segment linking e to O. All these elements are illustrated on Fig. S1.1.

S1.2 Ellipsoid silhouette in 3-D

For some unit vector d, let x be defined as

x = e + τ d (17)
with τ > 0 and d • n > 0. The ellipsoid silhouette as seen from e is given by the set of vectors d such that the half-line described by x when τ varies is tangent to the ellipsoid2 . The point x is on the ellipsoid if and only if

(d M d)τ 2 + (2d M e)τ + (e M e -1) = 0, (18) e O

S1.2

which is of the form

ατ 2 + βτ + γ = 0. (19)
Therefore, the half-line described by x is tangent to the ellipsoid if and only if Eq. ( 19) has a unique solution 3 , that is if and only if β 2 -4αγ = 0, which is equivalent to

d Sd = 0 ( 20
)
where S is equal to

S = M e e M + (1 -e M e)M. (21)
The ellipsoid silhouette is defined by the solutions of Eq. ( 20) respecting, as mentioned earlier, the following two conditions: |d| = 1 and d • n > 0.

Let S be defined as

S = 1 2 S. (22)
Note that S can be used in Eq. ( 20) in place of S. To propose a more explicit form of S , let us use the following block matrix formulation

M = M 11 M 21 M 21 m 33 (23)
where

M 11 is a 2 × 2-matrix, M 21 is a 1 × 2-vector,
Therefore

S = M 21 M 21 m 33 M 21 m 33 M 21 m 2 33 -m 33 - 1 2 M (27) = M 21 M 21 -m 33 M 11 (m 33 -m 33 )M 21 (m 33 -m 33 )M 21 (m 33 -m 33 )m 33 (28)
3 for completeness, note that this unique solution is τ = -β 2α .

S1.3

where m 33 is defined as

m 33 = m 33 - 1 2 . ( 29
)
So finally

S = M 21 M 21 -m 33 M 11 (1/ 2 )M 21 (1/ 2 )M 21 (1/ 2 )m 33 . ( 30
)

S1.3 2-D silhouette

As mentioned earlier, Eq. ( 20) defines the silhouette of the ellipsoid in a perspective projection setup. On plane Π, this silhouette is defined by points p such that, for all solutions d to Eq. ( 20),

p = e + τ d (p -e) • n = -δ (31)
where τ is a scalar 4 . The first equation of ( 31) is equivalent to

d = 1 τ (p -e). (32)
Equation ( 20) can now be rewritten in terms of p (and S as noted earlier) as follows (p -e) S (p -e) = 0.

The point p can be written as

p = O Π + q (34)
where q is a vector whose third component is equal to 0. Using a block formulation, we have

q = q 1 0 (35)
and

S = S 11 S 21 S 21 s 33 . (36)
Then, Eq. ( 33) is equivalent to

q 1 P q 1 + Qq 1 = r (37)
4 combining the two equations of (31) together, one gets τ = ( -δ)/(d • n) S1.4

where5

P = S 11 , (38)
Q = 2( -δ)S 21 , (39)
r = -( -δ) 2 s 33 . (40)
One recognizes the equation of an ellipse which can be put into the following standard form

(q 1 -c) 1 r -Qc/2 P (q 1 -c) = 1 (41)
where c = -P -1 Q /2 is the center of the ellipse.

S1.4 Semi-axes for perspective projection

Let λ 1 and λ 2 be the two (positive) eigenvalues of P , λ 1 ≤ λ 2 . Then the semi-minor and semi-major axes of the ellipse defined by Eq. ( 41) are

ρ i = r -Qc/2 λ i , i ∈ {1, 2}. (42)
Gathering everything together, ρ i can be rewritten in terms of M as follows

r = -(1 -δ/ ) 2 m 33 , (43)
Q = 2 -δ 2 M 21 , (44)
P = M 21 M 21 -m 33 - 1 2 M 11 , (45)
c = -P -1 Q /2, (46)
λ i = (tr(P ) + σ i √ ∆)/2, ( 47
)
|σ i | = 1 and σ 1 σ 2 = -1, (48)
∆ = tr(P ) 2 -4 det(P ) (49)
where tr(P ) is the trace of P , det(P ) is its determinant, and the σ i 's are chosen so that ρ 1 ≥ ρ 2 .

S2 Invariance of volume estimation error to scaling S2.1 Individual volume

This section has no direct practical application. It only presents developments useful in Appendix S2.2. For an axis-aligned ellipsoid, M is diagonal. The diagonal components are related to the semi-axes r i as follows

m ii = r -2 i , i ∈ {1, 2, 3}. (52)
The ellipsoid volume is then classically given by

V = 4 3 πr 1 r 2 r 3 (53) = 4 3 π √ m 11 m 22 m 33 . ( 54
)
For a general ellipsoid (i.e., any orientation), the volume is

V = 4 3 π det(M ) (55)
where det(M ) is the determinant of M . Let V * denote an estimation of the true volume V , where * is ESD or ELL here. The relative error in volume estimation is defined as

E * = V * V . (56)
For writing Eq. ( 56) for the M ESD , it should be reminded that, since the projection silhouette is an ellipse of area πρ 1 ρ 2 , the equivalent radius is equal to √ ρ 1 ρ 2 . Then, the relative errors of the M ESD or M ELL methods are

E ESD = (ρ 1 ρ 2 ) 3/2 det(M ) See Eq. (1) (57)
E ELL = ρ 1 ρ 2 2 det(M ). See Eq. (2) ( 58
)
The following section demonstrates the invariance of these individual volume estimation errors to scaling, a mathematical result that strongly impact the practical procedure (see Section "Simulation of copepod bodies").

S2.2 Invariance of individual volume estimation error to scaling

S2.2.1 Common remarks

The purpose is to show that E ESD and E ELL do not depend on the absolute volume of the ellipsoid, which is a function of (r 1 , r 2 , r 3 ), but rather on the ellipsoid proportions (r 2 /r 1 , r 3 /r 1 ). One way to prove this statement is to show that E * (αM ) = E * (M ) for any α > 0. Indeed, if this holds, then choosing α equal to r 2 1 implies that αM is defined by the triplet (1, r 2 /r 1 , r 3 /r 1 ). Let ρ, resp. λ, be a generic notation for ρ 1 and ρ 2 , resp. λ 1 and λ 2 . The other useful reminders are

ρ = m 33 λ (59)
λ : eigenvalue of P (60)

P = m 33 M 11 -M 21 M 21 . (61)
Let us add a subscript α to these quantities to denote their expressions when M is replaced with αM . We have

m 33,α = αm 33 (62)
P α = α 2 P. (63)
It is also clear that if λ is an eigenvalue of P , then βλ is an eigenvalue of βP (P x = λx ⇒ βP x = βλx) for any β = 0. Therefore,

λ α = α 2 λ. (64)
Hence, it can be concluded that

ρ α = ρ √ α . (65)
Finally, note that we have the following property on the matrix determinant

det(αM ) = α 3 det(M ) (66) if M is a 3 × 3-matrix.

S2.2.2 M ESD method

As a reminder, the relative error in volume estimation of the M ESD method is

E ESD (M ) = E ESD (r 1 , r 2 , r 3 , ξ) (67) = (ρ 1 ρ 2 ) 3/2 det(M ). (68)
S2.

2

Then

E ESD (αM ) = (ρ 1,α ρ 2,α ) 3/2 det(αM ) (69) = (ρ 1 ρ 2 ) 3/2 √ α 3 α 3 det(M ) (70) = E ESD (M ). (71)

S2.2.3 M ELL method

As a reminder, the relative error in volume estimation of the M ELL method is

E ELL (M ) = E ELL (r 1 , r 2 , r 3 , ξ) (72) = ρ 1 ρ 2 2 det(M ). (73)
Then

E ELL (αM ) = ρ 1,α ρ 2 2,α det(αM ) (74) = ρ 1 ρ 2 2 α √ α α 3 det(M ) (75)
= E ELL (M ).

(76)

S2.3 Invariance of total volume estimation error to scaling

Let V i , i ∈ [1.
.n], be a set of (true) ellipsoid volumes, and let V i * and E i * be some corresponding estimated volumes by the method " * " and the associated individual volume estimation errors, respectively. The total volume estimation error is

T * = i V i * i V i = V * V ( 77
)
where X denotes the average of X. Now, suppose that each ellipsoid volume is scaled by a factor α i (U i = α i V i ), for example as a result of the normalization of the ellipsoids by dividing their semi-axes r i 1 , r i 2 , and r i 3 by r i 1 . How will the estimated volumes U i * vary with respect to V i * ? From Appendix S2.2, we know that E i * = V i * /V i is invariant to ellipsoid scaling when S2.3 " * " is ESD or ELL. Therefore, U i * /U i must still be equal to E i * , which implies that U i * = α i V i * . Hence, the total volume estimation error after scaling is

T * = i U i * i U i = i α i V i * i α i V i = αV * αV . ( 78
)
The covariance between two random variables, say a factor α and a volume V , can be written in terms of expected values as follows

Cov(α, V ) = E[αV ] -E[α] E[V ]. (79)
The situation of practical interest here is when α is related to the normalization of the ellipsoids, in which case V and α are not independent (which would otherwise immediately guarantee that Cov(α, V ) = 0). Indeed,

α = 1 r 3 1 (80) V = 4 3 πr 1 r 2 r 3 . (81)
However, they are not correlated (i.e., there relationship is not linear). Consequently, their covariance is equal to zero and

E[αV ] = E[α] E[V ]. (82)
The same conclusion can be drawn regarding the estimated volume V * . If the number of samples n is large enough, these results can be safely extended to average values so that

αV ᾱ V (83) αV * ᾱ V * . (84)
In the end,

T * ᾱ V * ᾱ V = T * . ( 85
)
This concludes the proof that, if the number of involved volumes is high enough, the total volume estimation error is largely invariant to scaling. As a consequence, the randomly generated ellipsoids used to determine the total volume correction factor (see Sections "Simulation of copepod bodies" and "Corrected total volume", and, in particular, Eq. ( 12)) can be safely normalized, for example by dividing their semi-axes by their largest one.

S2.4

S3 Improved surface estimation and ellipse fitting S3.1 Common steps

The general idea of the improved methods proposed in this appendix is to get rid of the antennae (and tail) before measuring the copepod silhouette surface or fitting an ellipse onto it. It is assumed that the binary mask of the copepod has been determined previously. We propose to compute the Inner Distance Map (IDM) of this mask and to erode it using a threshold. We fix the threshold to max(IDM)/2.1 in our experiments. This step allows to get rid of the antennae. Note that the binary mask could have been eroded directly using mathematical morphology. However, it would make use of a discrete so-called structural element (typically a discretized disk), which would lead to a coarser eroded shape. Given the small size of a copepod in our images, this could have a negative impact on the subsequent steps.

Next, to recover the original copepod body size, the outer distance map of the eroded mask is computed and thresholded using the same threshold as the one used for erosion. This amounts to dilate the eroded mask, but again in a finer way than if using mathematical morphology. The various steps are illustrated in Fig. S3.1.

S3.2 Surface estimation

The copepod surface estimation is performed by counting the number of pixels of its binary mask. The improved version simply counts the pixels of the mask obtained in Appendix S3.1 as opposed to counting the pixels in the original binary mask which includes the antennae.

S3.3 Ellipse fitting

When an object is described by a binary mask of pixels, the most classical ellipse fitting method interprets the pixels as the samples of a point cloud.

The covariance matrix of the cloud is computed. Its eigenvectors represent the best fitting ellipse orientation while its eigenvalues represent the semiaxes of the ellipse. A simple improvement of this method (or any other ellipse fitting method as a matter of fact) consists in rescaling the ellipse so that its area matches the object area. This is implemented by the software ImageJ that ZooProcess (ZP), the software used to process UVP images, uses. However, if this improvement allows to correct the fitted ellipse sur-S3.1

S3.2

face (which can be enough for some applications), it does not help that much for copepod volume estimation. Indeed, the precision of the small semi-axis is crucial, and it is not improved by the surface adjustment. As a reminder, the M ELL estimation of the volume is:

V ELL = 4 3 πρ 1 ρ 2 2 = 4 3 πρ 1 ρ 2 Surface ρ 2 Minor semi-axis . (86)
Alternatively, the ellipse could be fitted on the grayscale version of the object, that is using the pixel intensities as sample weights when computing the covariance matrix. However, we found that this alternative does not work well on the copepod images of our data set. Whatever the ellipse fitting method is, the starting point is the copepod mask. The fitting methods get distracted by the antennae, which can result in very bad ellipses (see the red ellipses in Fig. S3.2). Therefore, we proposed to fit an ellipse on the mask obtained in Appendix S3.1 instead of the original binary mask which includes the antennae (see the green ellipses in Fig. S3.2). green and for the new proposed method in red. We see that when antennae are not visible, the result is almost the same, but when they are visible, the classic method is not appropriate.

S3.3

S3.4

S4 The proposed method, step-by-step

This section gathers the results of the different sections into a step-by-step procedure for estimating the total volume of copepods given a data set of 2-D views. It is composed of two stages: a learning stage which has to be performed once for all, or whenever the expert thinks the proposed simulation procedure must be adapted to the data, and a "usage" stage which can be applied at will.

S4.1 Learning stage

1. Generate random ellipsoid samples that realistically represent a generic population of copepods, or a population following some characteristics inferred from the data set. The randomness must be constrained by the expert knowledge in the form of specific simulation parameters.

2. Compute the total volume of the ellipsoid samples. This represents the true total volume. See Eqs. ( 4) and ( 12).

3. For each ellipsoid sample, compute the projection ellipse (see Eq. ( 7)) and the estimated volume using either the M ESD (see Eq. ( 1)) or the M ELL method (see Eq. ( 2)).

4. Sum all the estimated volumes to get the estimated total volume.

5. Compute the total volume estimation error T * from the true and estimated total volumes (see Eq. ( 12)). This is the final product of the learning stage.

S4.2 "Usage" stage

1. For each copepod image of a data set, determine the copepod silhouette using an image segmentation method. On UVP images, a simple binarization using a fixed threshold is enough.

1.a. For the M ESD method, compute the silhouette area A (see Appendix S3.2) and the corresponding estimated volume (see Eq. ( 1)).

1.b. For the M ELL method, fit an ellipse onto the silhouette (see Appendix S3.3). Let ρ 1 and ρ 2 be the semi-major and semi-minor axes, respectively. Then compute the corresponding estimated volume (see Eq. ( 2)).

2. Sum all the estimated volumes to get the estimated total volume W * where * is either ESD or ELL.

S4.1 3. Compute the corrected total volume estimation Ŵ * by dividing W * with T * from the learning stage (see Eq. ( 13)).

S4.3 Algorithmic details

S4.3.1 Uniformly random rotations

This section defines the rotation matrices used to simulate random orientations of ellipsoids.

In order to generate an ellipsoid with a uniformly random orientation, we generate a random rotation matrix R and rotate an axis-aligned ellipsoid with it. The generation of an axis-aligned ellipsoid is described in Section "Simulation of copepod bodies". If the axis-aligned ellipsoid is represented by a matrix M (see Section "Principle for correction of total volume"), then the rotated ellipsoid is represented by the matrix

M rot = RM R . (87)
A general rotation matrix can be defined using three elementary rotation matrices

R = R z (Φ)R y (Θ)R x (Ψ) (88)
with R i (α), i ∈ {x, y, z}, defines the rotation by angle α around axis i. To generate a random rotation matrix, one has to randomly choose the angle triplet (Ψ, Θ, Φ). In order to guarantee the uniformity of the ellipsoid orientations, the angles Ψ, Θ, and Φ must be distributed adequately, that is

Ψ = U [0, 2π[ (89)
Θ = arccos(1 -2U [0, 1[) - π 2 (90) Φ = U [0, 2π[ (91)
where U [a, b[ is the uniform distribution between a (included) and b (excluded).

S4.2

S5 Distribution of selected sample images

To define the real-world distribution of the semi axes of the ellipsoid representing the body of copepods (r 1 , r 2 , and r 3 ) as well as a the ratios between them, defining the shape of the ellipsoid (r 2 /r 1 and r 3 /r 1 ), 295 copepods seen from the side (on which r 1 and r 2 are measurable) and 265 copepods seen from the top or bottom (on which r 1 and r 3 are measurable) were manually curated from a collection of >150k images.

To make sure that these small samples were representative of the whole data set, we check their latitudinal and size (i.e. r 1 ) distributions. The shape of the latitudinal distribution of the side and top/bottom views matches well that of the total data set (Fig. S5.1). The side views shows and excess at high latitude, likely linked with a bias in the size distribution (see below; copepods are larger at high latitudes), and a linked under-representation elsewhere. The pattern is opposite for the top-bottom views. However, no region is completely missed in the samples and even some details of the distribution (such as the two peaks around -40 • ) are captured. Therefore, we consider them representative enough. The theoretical expectation for the length distribution is an exponential decay [START_REF] Sprules | Surfing the biomass size spectrum: some remarks on history, theory, and application[END_REF], i.e. a linear decrease, in log-scale. This is approximately true once the lower detection limit of the camera is passed, after ∼ 1 mm (Fig. S5.2). However, the distribution of side views shows S5.1 an excess in the size range 2 to 3.5 mm. This is likely due to the fact that telling that a copepod is viewed from the top/bottom can be determined from the geometry of its antennae relative to its body, no matter its size; making sure that a copepod is viewed from the side requires additional details, which are easier to asses on larger individuals, inducing a bias in the manual selection of images. As explained in the main text, this has little consequence on the estimation of the distribution of the semi-axes ratios (r 2 /r 1 and r 3 /r 1 ) but does now allow the estimation of the distribution of r 1 from these samples only. S5.2

Figure 3 :

 3 Figure 3: Examples of volume estimations and errors made by M ESD andM ELL for an ellipsoid E with (r 1 , r 2 , r 3 ) = (12, 5, 3) (see Fig.2for the definition of r i ). The first row displays the simulated ellipse from various viewing angles. The table gives, in black, the rounded values of the volume for each method and, in colour, the percentage of under/over estimation compared to the true value, which is computed from the r i s. The individual volume estimation error is invariant to scaling (see Appendix S2), therefore any ellipsoid respecting r 2 /r 1 = 5/12 and r 3 /r 1 = 3/12 would lead to the same error percentages. The lower right plot shows the volume computed with each method for viewpoints regularly sampled along an arc turning around the ellipsoid; note that the M ESD estimation is always greater than or equal to the M ELL one.

Figure 4 :

 4 Figure 4: Distributions of r 2 /r 1 (a) and r 3 /r 1 (b) fitted on our data set. The markers are the normalized histograms of the samples and the solid lines are the Gaussian Kernel Density Estimates.

Figure 5 :

 5 Figure 5: Total volume estimated by M ESD (blue) or M ELL (red) estimated with ZooProcess (ZP ) or our improved measures, uncorrected ( W ) or corrected ( Ŵ ). The discrepancy between methods is highlighted by the dashed line. The proposed correction does diminish this discrepancy, although there is no way to tell which result would be more realistic between M ESD or M ELL before correction.

Figure 6 :

 6 Figure 6: (a) Relationship between r 2 and r 1 from 254 copepods seen from the side, for r 1 >1 mm. (b) Relationship between r 3 and r 1 from the 173 copepods seen from the top or bottom, for r 1 >1 mm. The colored dashed line are linear regressions fits, significant in both cases (p<0.01, R 2 =88% for (a) and R 2 =75% for (b)). The color scale of points represents the density of samples.

Figure 7 :

 7 Figure 7: Estimation of the distribution of r 1 as C r × ρ 1 . The blue and orange points are the (normalised) histogram values. The black curve is the fourth-order polynomial fit in log-space on the blue points only (size range within which the UVP5 images copepods in a quantitative manner).

Figure 8 :

 8 Figure 8: Simulated Normalized Biovolume Size Spectra (NBSS): simulated ground-truth (GT) in black; estimation from M ELL (resp. M ESD ) in red (resp. blue). The x-axis is given in volume (mm 3 ) but also in Equivalent Spherical Diameter (ESD, mm) for comparability with other work. The dashed lines show the linear fits (which are offset vertically to improved readability.

Figure S1. 1 :

 1 Figure S1.1: Geometrical setup of the ellipsoid model. The camera is represented by its optical center e, the sensor plane Π, and the sensor orthonormal coordinate system (O Π , u, v). The global orthonormal coordinates system is represented by (O, i, j, k). The axes i and u are parallel, and so are j and v. The ellipsoid center is at distance δ from Π and from e. Without loss of generality (for our problem), the ellipsoid center is at O, and e is on the axis k with the optical axis aligned with k. Consequently, O Π is also on the axis k.

  and m 33 is a scalar. Using such a block formulation, we have denotes a matrix of zeros matching the dimension of M ij . Then M e e M = 2 M 21 M 21 m 33 M 21 m 33 M 21 m = 2 m 33 .

(a )

 ) Figure S3.1: Copepod body mask computation as a common preliminary step for surface estimation and ellipse fit. Reading the figure in lexicographical order, each image is the result of the processing of the previous one. They are: (a) the input grayscale image, (b) the binary mask obtained by thresholding, (c) the inner distance map, (d) the eroded mask obtained by thresholding, (e) the outer distance map, and (f) the dilated mask obtained by thresholding (same threshold).

Figure S3. 2 :

 2 Figure S3.2: Multiple examples of ellipse fit based on the original mask in

Figure S5. 1 :

 1 Figure S5.1: Kernel density estimate of the latitudinal distribution of the images of all copepods and of the side or top/bottom views.

  Figure S5.2: Distribution of the length of the semi-major axis of the ellipse fitted in the two views of the copepods. The vertical axis is the number of observations, in log 10 scale. The horizontal axis is the semi-major axis r 1 , which is equal to ρ 1 in these viewpoints and approximates the half of the prosome length, in millimeters.

Table 1 :

 1 Table of the estimated volumes. The "Gap" column represents the absolute value of the difference between M ESD or M ELL volumes. W : estimated volume without correction; Ŵ : volume corrected with the proposed method; ZP subscript: the surfaces and ellipse fits were performed using classical methods, as used in ZooProcess. Without ZP subscript: the surfaces and ellipse fits were performed using the proposed improved approaches (see Appendix S3).

		M ESD M ELL Gap
		Unit: ×10 5 mm 3
	W * ZP 7.51 4.79 2.72
	Ŵ * ZP 6.16 5.51 0.65
	W *	5.57 3.95 1.62
	Ŵ *	4.57 4.54 0.03

positive definite matrices are, by definition, symmetric

for such a vector d, there is a corresponding value for τ S1.1

in case P is definite negative, P , Q and r must be replaced with their opposite S1.5

This section has no direct practical application. It only presents developments useful in Appendix S2.3. S2.1

Acknowledgments

The authors acknowledge the projects ARTIFACTZ (joint Université C ôte d'Azur and Université Laval) and WWWPIC (ANR-18-BELM-0003-01 and ANR-18-BELM-0003-02) for financial support. They also thank A McDonnell, A Rogge, A Waite, L Coppola, L Stemmann, L Karp-Boss, L Guidi, F Lombard, M Babin, R Kiko, T Biard who supervised the collection of UVP5 images through many cruises, L Drago for her help in sorting the side vs top/bottom copepod images, F Lombard for his advice on the NBSS work. The software development for this work was made in Python using the following open source libraries: Matplotlib (Hunter, 2007), NumPy (Van Der Walt et al., 2011), OpenCV (Bradski, 2000), Pandas (McKinney et al., 2010), Scikit-Image ( Van der Walt et al., 2014), Scikit-Learn (Pedregosa et al., 2011), SciPy (Jones et al., 2001).

S1.5 Semi-axes for parallel projection

From the previous section, it follows that, for a parallel projection (i.e.

= ∞), the semi-minor and semi-major axes have the following simpler expression

and

while λ i , σ i , and ∆ are unchanged. Note that δ no longer appears in the equations.