Robots et mères-porteuses : la confusion des personnes et des choses

  M.-A. Frison-Roche, La disparition de la distinction de jure entre la personne et les choses : gain fabuleux, gain catastrophique, D. 2017. 2386 Dans un article d'une lucidité et d'une densité remarquables, Marie-Anne Frison-Roche montre comment se prépare la liquidation de la distinction juridique des personnes et des choses. Deux exemples frappants le montrent : le premier robot parlant devenu citoyen a été présenté comme une personne et les mères porteuses sont présentées comme des choses (parfois comme des « fours »). Pour l'auteur, la nouveauté ne porte pas sur les faits (il y a toujours eu des automates et des mères porteuses) mais sur l'idée qu'on s'en fait (p. 2387). Cette volonté de changer les choses en personnes et inversement a pour raison l'argent, à savoir le mécanisme de marché qui crée de nouveaux désirs (y compris sexuels avec des robots de ce type). Ce serait alors au politique d'affirmer que la distinction des personnes et des choses ne relève pas du fait mais du droit pour ainsi éviter toute confusion. Certes, le droit n'est pas un objet naturel, il est une construction. Mais il vient doubler l'ordre naturel pour veiller à sa préservation. Au sens propre du terme est une écologie : il veille au maintien des équilibres dans la relation entre les êtres vivants et le milieu social dans lequel ils vivent. La preuve du caractère construit et instrumental du droit s'évince de la possibilité de le mettre au service de finalités abominables, c'est l'exemple archi-rebattu du droit nazi. Aujourd'hui, la question est bien, comme l'explicite Anne-Marie Frison-Roche, de savoir si droit est au service de la personne ou du marché. Sa priorité est-elle l'être ou l'avoir ? Cette question appelle une hiérarchie des valeurs qui ne peut être donnée dans le droit car celui-ci est purement volontaire et construit. Mais, a contario, cela ne veut pas dire qu'avant le droit, il n'y a aucune valeur et seulement des faits. Selon une idée remarquable de Julien Freund, le droit est précisément le médiateur entre la morale et la politique, le moyen de trouver leur équilibre (Freund, Droit et politique. Essai de définition du droit, APD 1971. 15). Le droit exprime une transaction entre les valeurs et l'opportunité : la morale et la politique en sont les conditions de possibilité. Dire cela c'est refuser simultanément d'être un positiviste absolu qui considère que tout le droit est politique ou d'être un jusnaturaliste absolu qui pense un droit universel dans les choses ou dans la raison humaine. Bien au contraire, nous soutenons que même en dehors du droit l'être et l'avoir sont irréductibles, qu'il y a une confusion impossible entre ces deux ordres. S'il revenait au droit, et seulement à lui, de dire qu'on ne peut confondre les personnes et les choses, que pourra-t-on lui opposer s'il déclare le contraire ? La question ne se réduit ni à la politique ni à l'éthique car elle est aussi ontologique. Elle engage l'identité des êtres et des choses qui peuplent notre monde. Le problème devient juridique lorsque cette identité permet de répondre à la question de savoir si nous pouvons traiter dignement un robot ou vendre des enfants à naître. Notre vision des êtres et des choses nourrit le droit. Il serait dangereux de penser qu'il puisse avoir le dernier mot car le droit est autant épée que balance.

	Frédéric Rouvière
	Professeur à l'Université d'Aix-Marseille
	Laboratoire de théorie du droit

Si nous souscrivons sans réserve au regard pénétrant que porte l'auteur sur ce délicat problème, nous sommes plus interrogatifs à propos de son argument central à savoir que « la distinction entre la personne et les choses n'est pas naturelle, elle est juridique » (p. 2387). C'est de cet argument que découle sa conclusion : c'est au politique de prendre l'épée du droit et de poser un socle fondateur comme on frappe la table (p. 2389).

L'opposition présentée entre le droit et le fait nous conduit à questionner les rapports qu'ils entretiennent. Le droit, libre de la construction de ses catégories, peut-il ignorer des différences naturelles ? La question s'est déjà posée à propos de l'empire des fictions juridiques. Les canonistes ont voulu cantonner celles-ci à une imitation de la nature (notamment en matière d'adoption, Y. Thomas, Fictio legis : l'empire de la fiction romaine et ses limites médiévales, Droits 1995. 17). De la même façon, si la liberté du droit est totale politiquement (le droit peut dire qu'un homme est une femme, qu'un robot est une personne, etc.), cette liberté se mesure toutefois à l'aune du naturel. La nature donne le modèle qui permet de penser l'artifice de la volonté politique, un artifice qui se cristallise dans le droit. La distinction des personnes et des choses est alors juridique parce qu'elle est naturelle, parce qu'il paraît difficile de nier que les pierres et les plantes sont assimilables à l'humain et, par extension, que les machines soient humaines ou que les humains soient des machines.


How can we elaborate the SIT4PLANET, a generic syllabus that enables university stakeholders to embed mindfulness practices in any European course and campus activities ?

Technical objective: Creation of the first Mindfulness Toolbox for Faculties and Students.

Research objective: European survey to understand the main challenges and concerns faced by students in 2022 and survey of the existing mindfulness practices in education.

-Open Pass for Climate OP4C is an Erasmus+ funded project, 250 000 €.

Problematic : How can we elaborate the SIT4PLANET, a generic syllabus that enables university stakeholders to embed mindfulness practices in any European course and campus activities ? The first SDG to be shared by English (17 nationalities), French and Lithuanian students.


Jawher Jerray
email: jerray@lipn.univ-paris13.fr

Adnane Saoud
email: adnane.saoud@centralesupelec.fr

Laurent Fribourg
email: fribourg@lsv.fr

Asymptotic Error in Euler's Method with a Constant Step Size

Keywords:

Introduction

Different techniques have been proposed in the literature to provide a deeper understanding of different optimization algorithms through the use of continuoustime differential equations. Indeed, the use of differential equations makes it possible to exploit the plethora of tools developed within dynamical systems analysis and control community in order to analyse stability and convergence properties. In this context, the authors in [START_REF] Su | A differential equation for modeling nesterov's accelerated gradient method: Theory and insights[END_REF] derived a continuous-time ver-sion of the Nesterov accelerated gradient and used a Lyapunov formulation to estimate the convergence rate. Finally, they transformed the continuous time Lyapunov function into a discrete-time version and provided a new proof of the convergence rate of the Nesterov accelerated gradient method. A similar Lyapunov-based approach has been used to provide proofs of the convergence rates for the Nesterov accelerated algorithm in [START_REF] Shi | Understanding the acceleration phenomenon via high-resolution differential equations[END_REF] and the accelerated mirror descent method in [START_REF] Krichene | Accelerated mirror descent in continuous and discrete time[END_REF]. However, all the aforementioned approaches fail to provide a general discretization procedure that generate a provably convergent optimization algorithm, except of the work in [START_REF] Zhang | Direct Runge-Kutta discretization achieves acceleration[END_REF], where the authors used a Runge-Kutta scheme for a particular case of an accelerated Nesterov algorithm. Indeed, the main difficulty comes from the fact that designing a Lyapunov function of a continuous-time ordinary differential equation (ODE) is the "art of the designer".

In this paper, we explore another strategy that does not require an a priori knowledge of a Lyapunov function, but rather relies on the contractivity properties of the continuous-time differential equation. More precisely, we give here conditions, as general as possible, on the flow of ODEs in order to guarantee the convergence to 0 of the error in Euler's method, making use of an analytic formula δptq bounding this error at time t [START_REF] Le Coënt | Control synthesis of nonlinear sampled switched systems using Euler's method[END_REF].

Related work.

Let us mention that the concept of contractivity has been previously used in the literature, where the authors in [START_REF] Cisneros-Velarde | A contraction theory approach to optimization algorithms from acceleration flows[END_REF] were able to generate a contractive continuous-time version of the accelerated Nesterov algorithm.

Let us also note that an important part of the literature studies the convergence of the discretization error as the step size h Ñ 0. In contrast here, the size h is constant while the number of steps k Ñ 8. Such a convergence at constant h was the subject of pioneering work by [START_REF] Dahlquist | Error analysis for a class of methods for stiff non-linear initial value problems[END_REF] who developed the notion of A-stability and studied the disk around the origin within which the error remains bounded. However, these studies concern a particular type of linear differential equation of the type 9

x " ´αx with α ą 0 (cf. [START_REF] Efimov | Realization and discretization of asymptotically stable homogeneous systems[END_REF] for extension to "homogeneous systems"), and not arbitrary differential equations as here. [START_REF] Viswanath | Global errors of numerical ode solvers and Lyapunov's theory of stability[END_REF] provides an, linear in time, upper bound on the discretization error under some stability-type conditions on the solutions to the differential equation, for general numerical methods, withoud showing the convergence of the discretization error to zero. The authors in [START_REF] Calvo | Global error estimation with adaptive explicit Runge-Kutta methods[END_REF] provides an upper bound on the discretization error while using an explicit Runge-Kutta method with an adaptive and variable step size. Similarly, [START_REF] Kraaijevanger | Contractivity of Runge-Kutta methods[END_REF] shows the convergence of the error dynamics while using Runge-Kutta method and for a fixed step size under the restrictive circle condition [START_REF] Dahlquist | Generalized disks of contractivity for explicit and implicit Runge-Kutta methods[END_REF]. In comparison with existing approaches in the literature, in this paper, we provide sufficient conditions in terms of strong monotonicity and co-coercivity to show the convergence of the discretization error to 0 for the explicit Euler's method, while using a constant step size.

Regarding the analysis of the gradient descent algorithm, it is well-known that the error (between the sequence generated by the algorithm and the optimal point) is non-decreasing when the step size satisfies h ă 2{L, where L is the Lipschitz constant of the gradient [START_REF] Dahlquist | Error analysis for a class of methods for stiff non-linear initial value problems[END_REF][START_REF] Sanz-Serna | Contractivity of Runge-Kutta methods for convex gradient systems[END_REF][START_REF] Nesterov | Introductory Lectures on Convex Optimization: A Basic Course[END_REF]. The same constraint is used here for the gradient descent (see Corollary 1), but, when coupled with an assumption of strong monotonicity, we show that the error tends to 0. We also make use of the property of co-coercivity (Theorem 3) whose role has already been observed by Zhu and Marcotte [START_REF] Zhu | Co-coercivity and its role in the convergence of iterative schemes for solving variational inequalities[END_REF] in the different context of iterative schemes for solving variational inequalities.

The remainder of this paper is organized as follows. Section 2 recalls some classical results on convexity, strong monotonicity and co-coercivity. Section 3

gives some sufficient conditions on the flow of the differential equation and the discretization time-step to make the error in Euler's method converge to 0. We give some final remarks in Section 4.

Preliminaries

Notation. The symbols N, N ą0 , and R and R ą0 denote the set of positive integers, non-negative integers, real and non-negative real numbers, respectively.

The Euclidean norm is denoted by }.}. For a differentiable map x Þ Ñ f pxq P R, ∇f denotes the gradient of f with respect to x. For an invertable matrix A P R mxˆmx , A J denotes its transpose. For x P R n and r ą 0, Bpx, rq " tz P R n : }x ´z} ď ru.

We recall some definitions and classical results on convexity, monotonicity and co-coercivity (see, e.g., [START_REF] Nesterov | Introductory Lectures on Convex Optimization: A Basic Course[END_REF][START_REF] Boyd | Convex Optimization[END_REF][START_REF] Ryu | A primer on monotone operator methods[END_REF]) that will be used through the paper.

Definition 1. Consider a continuously differentiable function f : R n Ñ R. The function f is convex if, for all x, y P R n :

f pyq ´f pxq ě p∇f pxqq J py ´xq.

If f is convex and bounded from below, then a minimum of f exists in some

x ˚P R n , and ∇f px ˚q " 0.

Definition 2. A function g : R n Ñ R n is • L-Lipschitz continuous if for all x, y P R n :
}gpxq ´gpyq} ď L}x ´y}.

• strongly monotone if there exists m ą 0 such that, for all x, y P R n : pgpxq ´gpyqq J px ´yq ě m}x ´y} 2 .

(1)

• co-coercive if there exists a positive constant a such that for all x, y P R n : pgpyq ´gpxqq J py ´xq ě a}gpyq ´gpxq} 2 .

Remark 1. Property (1) of strong monotonicity is equivalent to saying that g has a negative one-sided Lipschitz constant (OSL) λ (see [START_REF] Dahlquist | Error analysis for a class of methods for stiff non-linear initial value problems[END_REF][START_REF] Hairer | Solving Ordinary Differential Equations II[END_REF]; cf Definition 3 below), equal to ´m.

We consider a differential system of the form

9 xptq " gpxptqq (2)
with initial condition x 0 P R n , where g : R n Ñ R n is a differentiable Lipschitzian function of constant L ą 0. Without loss of understanding, we will denote by xpt; x 0 q, or more simply by xptq, the solution of this system at time t.

We denote by y k the (explicit) Euler discretization of 9

xptq " gpxptqq at time t " kh, for k P N, where h ą 0 is a constant step size. Given an initial point y 0 P R n , y k is defined, for k ě 1, by:

y k " y k´1 `hgpy k´1 q (3)
Proposition 1. On a domain where ´g is strongly monotone (i.e., g has an OSL constant λ ă 0), any two solutions xpt; x 0 q and xpt; y 0 q of (2) with initial conditions x 0 and y 0 respectively, converge to each other exponentially in time:

}xpt; x 0 q ´xpt; y 0 q} ď }x 0 ´y0 }e λt .

It follows that there is at most one stationary point x ˚(with gpx ˚q " 0).

This result states that strong monotonicity leads to a form of contraction (as studied in [START_REF] Cisneros-Velarde | A contraction theory approach to optimization algorithms from acceleration flows[END_REF][START_REF] Wensing | Cooperative adaptive control for cloud-based robotics[END_REF]).

The following result can be found in [START_REF] Chambolle | Continuous optimization, an introduction[END_REF] and allows to relate the concepts of convexity, Lipschitzness and co-coercivity.

Proposition 2. Consider a continuously differentiable function f : R n Ñ R.

If f is convex and ∇f L-Lipschitz, then ∇f is co-coercive of constant 1{L.

Proposition 3. Consider a function g : R n Ñ R n . If the following conditions are satisfied:

1. ´g is co-coercitive,
2. the step size h ă 2{L, 3. gpx ˚q " 0 for some x ˚P R n (existence of a stationary point).

Then }gpy k q} Ñ 0 with rate Op1{kq for the averaged iterates, where ty k u kPN is defined as in [START_REF] Krichene | Accelerated mirror descent in continuous and discrete time[END_REF].

A proof (adapted) from [START_REF] Gower | Convergence theorems for gradient descent[END_REF] is given in Appendix in order to be self-contained.

In the following, we show how the result of Proposition 3 applies to the case of the gradient descent algorithm. We first recall the link between the gradient flow and gradient descent algorithm.

Consider a function f : R n Ñ R, the gradient descent algorithm generates a sequence tx k u kPN described as:

x k`1 " x k ´h∇f px k q (4)
where h ą 0 is a constant step size. This algorithm is generally used to resolve optimization problems of the form min xPR n f pxq for a function f . The algorithm in (4) can be seen as the (explicit) Euler discretization of the gradient flow differential equation described by: 9

x " ´∇f pxq

for a continuously differentiable function f : R n Ñ R.

From Proposition 2 and Proposition 3, it follows:

Proposition 4. Consider a continuously differentiable function f : R n Ñ R.

If the following conditions are satisfied

1. f convex of minimizer x ˚,
2. ∇f L-Lipschitzian, 3. the step size satisfies h ă 2{L.

Then }∇f px k q} Ñ 0 as k Ñ 8, with rate Op1{kq for the averaged iterates, where tx k u kPN is defined as in (4).

Asymptotic Error in Euler's Method

Error bound in Euler's method

Our objective is to analyze the error between the solutions of the differential equation ( 2) and their Euler approximations in [START_REF] Krichene | Accelerated mirror descent in continuous and discrete time[END_REF]. For this reason, we define a function δ : R Ñ R in order to upperbound }y k ´xptq} for t " kh (see [START_REF] Le Coënt | Control synthesis of nonlinear sampled switched systems using Euler's method[END_REF] and Theorem 1 below).

Definition 3. Given a sequence tµ k u kě0 of nonnegative reals, we define δ µ k : r0, hs Ñ R ě0 as follows: for all t P r0, hs, if λ ă 0 :

δ µ k ptq " ˆµ2 k e λt `C2 k λ 2 ˆt2 `2t λ `2 λ 2 `1 ´eλt ˘˙˙1 2 if λ " 0 : δ µ k ptq " ´µ2 k e t `C2 k `´t 2 ´2t `2pe t ´1q ˘¯1 2 if λ ą 0 : δ µ k ptq " ˆµ2 k e 3λt `C2 k 3λ 2 ˆ´t 2 ´2t 3λ `2 9λ 2 `e3λt ´1˘˙˙1 2
where

• L denotes the Lipschitz constant for g,

• C k :" L}gpy k q},
• λ is the "one-sided Lipschitz constant (OSL)" (or "logarithmic Lipschitz constant" [START_REF] Aminzare | Contraction methods for nonlinear systems: A brief introduction and some open problems[END_REF]) associated to g, i. e., the minimal constant such that, for all x 1 , x 2 P R n :

pgpx 1 q ´gpx 2 qq J px 1 ´x2 q ď λ}x 1 ´x2 } 2 , (6)
Note that the constant C k depends on the value of yptq at time t " kh4 .

The constant λ can also be calculated "locally" for the zone occupied by xptq at time t P rkh, pk `1qhs, using a nonlinear optimization solver (see, e.g., [START_REF] Jerray | Robust optimal periodic control using guaranteed Euler's method[END_REF] for details).

Let us now consider the sequence tµ k u kě0 where µ k is defined recursively, for k ě 1 as

µ k " δ µ k´1 phq.
We can now define δ µ0 p2hq as δ µ1 phq with µ 1 " δ µ0 phq, and more generally, we can extend the definition of δ µ0 p¨q on r0, 8q as follows: for all k ě 0 and t P r0, hs, δ µ0 pkh `tq " δ µ k ptq.

In particular, it is easy to see that: δ µ0 pkhq " δ µ k p0q " µ k and δ µ0 ppk `1qhq "

δ µ k phq " µ k`1 .
Theorem 1.

[5] Consider the system 9 xptq " gpxptqq, with g : R n Ñ R n . Let xptq and py k q be defined as in (2) and (3) respectively, and µ 0 :" }y 0 ´x0 }.

Then, for all t " kh with k ě 0:

}y k ´xptq} ď δ µ0 ptq.
In the following, we give sufficient conditions on g and h ensuring that δ µ0 pkhq (hence }y k ´xpkhq}) converge to 0 as k Ñ 8. In the rest of the paper, we will always assume that g is Lipschitz continuous with constant L ą 0.

Strong monotonicity

Lemma 1. Consider the system 9 xptq " gpxptqq, with g : R n Ñ R n . Let xptq and py k q be defined as in (2) and (3) respectively, and µ 0 :" }y 0 ´x0 }. If the following conditions are satisfied:

1. }gpy k q} Ñ 0 as k Ñ 8 with a convergence rate r k , 2. g of OSL constant λ ă 0 (i.e., ´g strongly monotone).

Then δ µ0 pkhq Ñ 0 as k Ñ 8 with a convergence rate r k .

Proof. By Taylor-Lagrange's theorem, we find:

δ 2 µ0 ppk `1qhq ď δ 2 µ0 pkhqe λh `C2 k h 3 {p3|λ|q. ( 7
)
Since C k " L}gpy k q} Ñ 0 as k Ñ 8 with rate r k , and e λh ă 1, it follows that there exists ℓ P N such that: δ 2 µ0 ppk `1qhq ă δ 2 µ0 pkhq for all k ě ℓ. Hence, the sequence tδ µ0 pkhqu kěℓ is decreasing, and converges to a limit δ ˚ě 0 with rate r k . This limit δ ˚must satisfy (7), i.e.:

pδ ˚q2 ď pδ ˚q2 e λh `C2 k h 3 {p3|λ|q.
Since C k " L}gpy k q} Ñ 0 as k Ñ 8 and e λh ă 1, one must have δ ˚" 0.

Therefore δ µ0 pkhq Ñ 0 with rate r k .

Remark 2. Let us mention that the result in Lemma 1 can be relaxed to show that the sequence δ µ0 pkhq, k P N, is bounded if the sequence }gpy k q}, k P N is bounded.

From Lemma 1 and Theorem 1, it follows:

Theorem 2. Consider the system 9 xptq " gpxptqq, with g : R n Ñ R n L-Lipschitz continuous. Let xptq and py k q as defined in (2) and (3) respectively. Suppose:

1. }gpy k q} Ñ 0 with a convergence rate r k , 2. g of OSL constant λ ă 0 (i.e., ´g strongly monotone).

Then δ µ0 pkhq Ñ 0 and }y k ´xpkhq} Ñ 0 with a convergence rate r k as k Ñ 8.

Note that Theorem 2 applies even if }y k } Ñ 8, as shown in Example 1.

Remark 3. For the sake of simplicity of the statement, we assume in Theorem 2 (as well as Theorem 3 and Corollary 3 below), that the property of Lipschitz continuity and strong monotonicity (as well as, later on, co-coercivity and existence of stationary point) hold over all R n . In fact, it is enough that these properties hold on a subdomain D Ă R n only, provided that D be "invariant" in the following sense:

Dℓ P N @k ě ℓ : Bpy k , δ µ0 pkhqq Ď D 5 .
Example 1. Let gpxq " 1 2 ?

x , y 0 " 1, h " 0.5, µ 0 " 0.5, and let D " r1, `8q. We have y k`1 " y k `h 1 2 ? y k ą y k , therefore ty k u is an increasing sequence which converges to `8, and gpy k q " 1 2 ? y k converges to 0. Using the software ORBITADOR [START_REF] Jerray | Orbitador: A tool to analyze the stability of periodical dynamical systems[END_REF], one calculates that L ď 0.5 and λ ă 0 on D (see Figure 1). Since }gpy k q} Ñ 0 (Figure 2), it follows by Theorem 2: δ µ0 pkhq Ñ 0 and |y k ´xpkhq| Ñ 0 as k Ñ 8 (see Figures 3 and4).

Application to gradient descent

An immediate consequence of Theorem 2, using Lemma 1 and Proposition 4 with g " ´∇f , is the following result (which avoids mentioning }gpy k q} Ñ 0 as an explicit assumption).

Corollary 1. Consider a continuously differentiable function f : R n Ñ R.

Suppose:

5 Bpx, rq denotes the ball of centre x P R n and radius r ě 0. Example 2. 6 Let us consider the following function f : R 2 Ñ R defined for 6 adapted from t ě 0, by:

f pwq " 1 2 Σ 2 i"1 pb i ´pw 1 ¨ai `w2 qq 2 ,
with w " pw 1 , w 2 q, a 1 " 1, a 2 " 2, b 1 " 5 and b 2 " 7. Let us compute its minimum w ˚" pw 1 , w 2 q by gradient descent, with initial point w 0 " pw 0 1 , w 0 2 q " p9, 10q. Using the software ORBITADOR [START_REF] Jerray | Orbitador: A tool to analyze the stability of periodical dynamical systems[END_REF], we find L ă 7 and λ ď ´0.187 (see Figure 5). For h " 0.2 (which satisfies: h ă 2{L), it follows by Proposition 4 that }∇f pwpkhqq} Ñ 0 and C " L}∇f pwpkhqq} Ñ 0 respectively, as k Ñ 8

(see Figures 6 and7 respectively). Letting µ 0 " 0.1, it follows from Corollary 1 that δ µ0 ptq Ñ 0 (see Figure 8 for t " kh ď 95 and Figure 9, with a change of scale, for 95 ď t " kh ď 120), y k Ñ w ˚and xpkhq Ñ w ˚, with w ˚" p2, 3q (which is the stationary point of ∇f and the minimizer of f ).

Co-coercivity

The following result is a generalization of Corollary 1.

https://codingvision.net/gradient-descent-simply-explained-with-example 4. gpx ˚q " 0 for some x ˚P R n (existence of a stationary point).

Then we have:

• x ˚is the unique stationary point of R n ,

• y k Ñ x ˚and xpkhq Ñ x ˚as k Ñ 8 with rate Op1{kq for the averaged Example 3. Consider the differential equation 9

x " gpxq with gpxq " ´4x 3 6x BITADOR, one calculates L ď 12, where L is the Lipschitz constant of g. Let D " r1.25, 1.75s. For µ 0 " 0.1 and h " 0.12 ă 2{L, one can show, using ORBITADOR:

• λ ă 0 on D (see Figure 12),

• Bpy k , δ µ0 pkhqq Ď D for all k ě 12, and

• ´g co-coercive of constant 1{L on D.

Besides, x ˚" 1.5 P D is a stationary point (gpy ˚q " 0). From Theorem 3, it follows: δ µ0 pkhq Ñ 0, (see Figure 13), x ˚is the unique stationary point of D, y k Ñ x ˚and xpkhq Ñ x ˚as k Ñ 8 (see Figure 14). We can also check C " L}gpy k q} Ñ 0 (see Figure 15). Note that x ˚is here the minimizer of the non-convex function f pxq " x 4 ´2x 3 `2 (with ´∇f pxq " gpxq). See Figure 16.

Conclusion

In this paper, we have focused on the following problem: given an ODE 9

x " gpxq and the sequence py k q generated by its discretization with explicit Euler's method, what are the sufficient conditions on g to ensure that the discretization error }xpkhq ´yk } converges to 0 as k Ñ 8. In particular, under certain properties of g called "strong monotonicity" and "co-coercivity", we have shown the convergence of the discretization error to 0 (Theorem 3). This can shed new light on the relationship between the convergence of continuous differential according to a Gaussian distribution, then the gradient descent converges to a globally optimal solution. In future work, we will use the approach developed in this paper to show the convergence of the gradient descent algorithm for the class of neural networks considered in [START_REF] Du | Gradient descent provably optimizes over-parameterized neural networks[END_REF].

Figure 1 :

 1 Figure 1: Evolution of λ

Figure 2 :

 2 Figure 2: Evolution of C k " L}gpy k qq}

Figure 3 :

 3 Figure 3: Evolution of δµ 0

Figure 4 :

 4 Figure 4: Evolution of y k

Figure 5 :

 5 Figure 5: Evolution of λ computed locally (λptq ď ´0.187)

Figure 6 :

 6 Figure 6: Evolution of }∇f pwptqq}

Figure 7 :Theorem 3 .

 73 Figure 7: Evolution of Cptq " L}∇f ptq} varying from 918 (for t " 0) to 0.0272 (for t " 30)

Figure 8 :Figure 9 :

 89 Figure 8: Evolution of δptq for t P r0, 95s

Figure 10 :Figure 11 :

 1011 Figure 10: Evolution of w 1 ptq (which converges to w 1 " 2)

Figure 12 :

 12 Figure 12: Evolution of λ

Figure 13 :Figure 14 :

 1314 Figure 13: Evolution of δµ 0 (which converges to 0)

Figure 15 :

 15 Figure 15: Evolution of C " L}gpyqq}

Figure 16 :

 16 Figure 16: Graph py k , f py k qq

This slight modification of the definition of C, as originally given in[START_REF] Le Coënt | Control synthesis of nonlinear sampled switched systems using Euler's method[END_REF], is justified by a simple inspection of the proof of Theorem 1 in[START_REF] Le Coënt | Control synthesis of nonlinear sampled switched systems using Euler's method[END_REF].

, and its Euler discretization with y 0 " 0.25 and h " 0.12. Using OR-

Appendix: Proof of Proposition 3

Proof. This proof Let α :" Lh ą 0. We have:

Now, from the co-coercivity of ´g and gpx ˚q " 0, we have:

Therefore, from (8), we derive:

Therefore, by dividing by αp2 ´αq (which is positive since 2 ´α " 2 ´hL ą 0):

with β " αp2 ´αq ą 0.

Summing up from k " 1, . . . , K and dividing by K gives

It follows that }gpy k q} Ñ 0 as k Ñ 8 with rate Op1{kq for the averaged iterates.


Physically accurate spectral rendering generally brings further computational complexity to RGB rendering pipelines which precludes real-time rates. In this report, we show that by building on several physical a-prioris for the spectral phenomena in standard astronomical scenes, we achieve to obtain a model greatly reducing the complexity of the spectral aspect for volumetric rendering -by up to 83 factor in our examples. By projecting the total local contribution on a small optimal basis, thus providing scalability for the method, our rendering scheme then makes it possible to render a volumetric scene with any number and spectra of lights in real-time and no significant visual discrepancies.

Acknowledgement

I would like to give my sincere thanks to my supervisor Fabrice Neyret for his help throughout the project, and the opportunity to work on such interesting research question at Maverick. His constant availability for questions, discussions, and feedback truly enriched my understanding of the subject and the quality of my work.

I also would like to thank Antoine Richermoz for developing and sharing his Vulkan framework shaderV with which the shaders of this project were run, providing an efficient tool that helped me get results without having to build a whole GPU framework from scratch during these five months.

Résumé

Le rendu spectral physique implique généralement une augmentation de la complexité de calcul par rapport aux pipelines de rendu RGB, ce qui exclut les possibilité de rendu en temps réel. Dans ce rapport, nous montrons qu'en nous appuyant sur plusieurs a-prioris physiques pour les phénomènes spectraux dans les scènes astronomiques standards, nous parvenons à obtenir un modèle réduisant considérablement la complexité de l'aspect spectral pour le rendu volumétrique -d'un facteur allant jusqu'à 83 dans nos exemples. En projetant la contribution locale totale sur une petite base optimale, assurant ainsi l'extensibilité de la méthode, notre schéma de rendu permet alors de rendre une scène volumétrique avec n'importe quel nombre et spectre de lumières en temps réel et sans défauts visuels significatifs.

Introduction

This work is part of a long-term project involving the ANR Galaxy/veRTIGE between INRIA, RSA-Cosmos and Observatoire de Paris-Meudon [Anr]. The goal of this project is to render complex and physically accurate astronomical scenes in real-time for entertainment and educational purposes. We focus here in particular on the accurate spectral rendering of astronomical scenes.

In astrophotography, the process of obtaining a picture often involves stacking and compositing a set of pictures taken using specific filters, as illustrated in figure 1.1. Each of these filters puts emphasis on a specific (large or narrow) range of wavelengths in the scene during acquisition (ex: ultraviolets, infrared, x-ray), which corresponds to summing the intensity Figure 1.1: Stacking process for a picture of NGC 1512 combining the image acquisition obtained using 7 different filters Image credit: European Space Agency, HST contribution of all the ray lengths in the range of interest. This process is too complex to be accurately computed by traditional rendering pipelines using RGB approximation. This approximation assumes that all the light of the world and its interactions are only made of three pure colors (pure red, green and blue), while to be able to accurately render complex scenes involving non-linear light transformation (as in the case with nebulae) the range of wavelengths needs to be considered as a continuous spectrum. There is thus already a real challenge in trying to produce accurate renders of such astronomical objects without considering the real-time constraints. Moreover, we want to still be able to explore a significant variable space during render: the ability to freely move in the space around the nebulae is a prerequisite in order to get a compelling interactive exploration, but we would also want to be able to play with the filter spectra at run-time, similarly to how an operator would change a filter on the JWST in order to take a picture focusing more on an ultraviolet range, or an infrared one. This work's purpose is thus to produce a complete model for nebulae rendering, taking into account all of these previous requirements, all of that under real-time constraints. In section 2, we introduce in more detail the spectral rendering problem and describe our requirements for a relevant solution, and look at the literature approaches to this problem. Then, we present in section 3 an extensive mathematical description to model the spectral phenomenon (reflection, emission, attenuation, filtering) taking part in an arbitrary astronomical scene, basing our models on physical observations. We use this model in section 4 to propose an efficient method for the spectrally accurate computation of these scenes under real-time constraints. We also show and evaluate the results of our implementation of the method and its performance.

2

State of the art 2.1 Spectral volumetric rendering

Nebulae are astronomical objects made of dust, which can be seen as participating media, meaning that the medium by itself interacts directly with the light-rays passing through it. Due to this nature, nebulae cannot be accurately rendered using traditional surface-based rendering techniques.

Unlike those methods, volumetric rendering considers the properties (density, absorption...) of the medium filling the space. This allows for the accurate representation of phenomena such as fog, smoke, and other participating media that interact with light, like in our case nebulae.

The Volumetric Rendering Equation (VRE) [START_REF] Kajiya | Ray Tracing Volume Densities[END_REF][START_REF] Kajiya | The Rendering Equation[END_REF] describes the radiance received on a point x along a viewing direction w. Computing this equation on an array of pixels1 produces an image of the scene. In most applications, the VRE equation 2.2 is computed for each pixel over three color channels (RGB), associated with characteristic coefficients depicting the scattering, extinction, and transmission properties of the material over three pure wavelengths (red, green, blue) to approximate human vision.

L local,λ (t, w) = T (t, λ )[ σ s (t, λ )L s (t, w, λ ) + σ a (t, λ )L e (t, w, λ )]

T (t, λ ) = exp - t u=0 σ t (u, λ )du (2.1)
T is the light extinction from the camera to point t, σ a , σ s and σ t are respectively the light absorption, scattering, and extinction coefficients, L s is the external light arriving at a point and scattered in a particular viewing direction (its expression can be more or less complex, modelling multiple bounces through the medium or not) and L e if the light emitted at point t directly by the material.

L λ (x, w) = E t=0 L local,λ (t, w)dt (2.2)
The validity of the RGB model for rendering depends especially on two assumptions: the rendering pipeline should only use linear transformations on the spectra, and it should be possible to express any color as a sum of the three RGB base colors. However, these assumptions are often no longer true in most scenes rendered with volumetric rendering. As the generic extinction expression depends on the wavelength of the ray, results can indeed greatly differ between the RGB approximation and the expected render for the accurate simulation of arbitrary bodies. Two metamers cannot be visually differentiated without considering the spectrum as a continuous function over the wavelengths (the RGB representation not being able to differentiate "pure" yellow from a sum of red and green). Light refraction in glass or water is another example of a common physical phenomenon which cannot be simulated using a tri-ray model for the spectrum due to dispersion (visible with rainbows).

The relevance of taking the spectral aspect into account becomes even more important in the case of the study of astronomical objects, as is the case in this report. Astronomers often try to focus on specific wavelength ranges when acquiring astronomical photographs, in order to better discern the contribution of certain elements that have more impact at those wavelengths.

The Spectral Volumetric Rendering Equation (SVRE) written in equation 2.3 is a more general equation which takes the wavelength dependence into account, integrating the spectrum over the sensor sensibility.

L(x, w) = ∞ λ =0 S(λ )L λ (x, w)dλ (2.3)
In practice, the SVRE introduces an additional computational complexity based on the spectral integral. In real-time applications of the SVRE [START_REF] Bergner | Interactive Spectral Volume Rendering[END_REF], the incoming color is often numerically integrated (sampled) along the two dimensions (wavelength and ray length), which makes the outputs rely even more on numerical resolution than regular volumetric rendering, a value linked with the number of computations done for each integration step. Moreover, this discretization of the space can introduce color artefacts like spectral aliasing.

In the case of real-time rendering of astronomical objects, the need for speed and interactivity poses additional constraints on the rendering pipeline. Real-time rendering requires that images be generated in a matter of milliseconds: this often precludes the use of computationally expensive techniques such as spectral rendering, since ray-marching to compute the spatial integral of the VRE 2.2 is already a costly process.

Basis functions for analytical integration in literature

To drastically decrease the complexity drawback of the additional spectral integral, a relevant technique is to try to compute analytically the spectral aspect. This would completely nullify the drawbacks of having to keep a spectral discretization, and hopefully enable the computation of the integral between any two bounds in constant time. The most generic form of the SVRE generally isn't integrable analytically, but a common approach in that case is to introduce basis functions to model the different terms of the SVRE. By carefully defining a set of basis functions that can be integrated analytically over a "difficult" integral (and ensuring that their integrability will remain stable over their composition), the computational load of the integral can be significantly lowered. In the case of spectral rendering (but not necessarily volumetric), a few works have proposed several basis sets to model the involved spectra and ensure a fast and accurate integration (for example using polynomial [START_REF] Raso | A Piecewise Polynomial Approach to Shading Using Spectral Distributions[END_REF], Fourier series [START_REF] Peercy | Linear Color Representations for Full Speed Spectral Rendering[END_REF], or more generic spectrum decomposition using PCA [START_REF] Otsu | Reproducing Spectral Reflectances From Tristimulus Colours[END_REF][START_REF] Waddle | Real-Time Spectral Rendering[END_REF]). However, most of these sets are defined for surface rendering applications and do not transfer well to a volumetric rendering, where the spectrum is composed in different ways (products and sums). Moreover, the SVRE introduces terms depending on both the spatial coordinate and the spectral coordinate (generally in a non-separable way), which requires a different model than one depending only on the wavelength.

There are several characteristics that need to be taken into account in order to build such a basis set, the main requirement being that it should be integrable, i.e. have a simple analytical solution to its spectral integration. It should also be robust, meaning that compositions of the elements of the basis set should also remain integrable. This will enable huge computational gains and decrease the time needed to compute a picture, which will be measured at the end of the report. Finally, it should be generic enough but still able to accurately fit and model each term of the SVRE 2.3 that are to be spectrally integrated thanks to it. By building on a set of assumptions obtained from the study of experimental and theoretical data of nebulae and cosmic dust clouds, we will design such a basis set in section 3 that will ensure an analytical solution to the spectral integral part of the SVRE, enabling the possibility to greatly decrease the computational load of the volumetric integration while still providing accurate results in the case of astronomical bodies.

Three desirable characteristics for our optimal rendering pipeline

We list below three desirable characteristics of a good renderer:

• Fast: The renders should be computed in a reasonable time on reasonable hardware, with real-time rendering being the best goal

• Precise: Accurate computations, as close as possible to the simulated object

• Interactive: The parameters of the renderer should be editable (No long pre-computations should depend on these parameters)

A volumetric renderer able to reach a maximum score in the three characteristics is hard to design (as precision and speed can often be incompatible), and areas of focus can greatly vary regarding the requirements of a project (for example, a model needed to freely explore a parameter space will have a stronger focus on Speed and Interactivity).

The purpose of this work is to be able to obtain an interactive render method making it possible to explore the parameter space of most of the spectral functions of the SVRE, and still be able to compute the full integral in a fast and precise way. In effect, the idea would be to be able to emulate the rendering process of telescope imagery, being able to change the filters and explore a complex scene in real-time.

Astrophysical assumptions

In the following section, the typical astrophysical scene is described, as well as the different spectra that will contribute to the SVRE 2.3 and the intensity obtained on a camera sensor (in practice, a telescope sensor and a large bank of filters). Models are presented to accurately describe those spectra and justify the functions that are used in the contribution section 4.

Scene composition and nebulae

A typical scene in astrophysical rendering contains stars as light sources and a participating medium, in the case of this work an astronomical dust cloud or nebula. Nebulae can be of three types (illustrated in figure 3.1): • Reflection nebulae can reflect the incoming light of nearby stars.

• Emission nebulae are clouds of ionized gases which are also directly emitting light. The emitted light is restricted to certain wavelengths linked to the recombination of ionized atoms in the nebula cloud.

In effect, many nebulae are potentially a combination of these three optical phenomena. However, one of these phenomena often holds a greater magnitude in the final light contribution which eventually gives its type to the nebula. For example, we can simulate an absorption nebula by considering that the extinction term of the SVRE 2.3 is pre-dominant in front of the reflection and emission terms of the equation.

In other words, any nebula can be seen as made up of three kind of contributions that can be gradually added together to form a more complete and accurate render. The dark nebula is the simplest model by mainly being contributed by the extinction of light, the reflection nebula adds on top of that the reflection of nearby stars, and emission adds a local emission term. Each of these three terms can be seen in the VRE equation 2.2 and treated additively to account for a complete render.

Mathematical models

As the scope of this project is to efficiently compute the double spectral-spatial integral of the SVRE 2.3 by analytically integrating the spectral part, the first step is to propose a mathematical model for the terms of the integral depending on the wavelength λ and the spatial coordinate t. For each phenomenon involved, we present the classical physical equation or experimental data describing it, and then a mathematical function to model it. All the following models are summed up in table 3.1. Figure 3.2 illustrates the attenuation spectra of various astronomical dust clouds. For high wavelengths (λ > 0.25 µm), the spectra appear to be linearly correlated with the inverse of the wavelength. This range greatly coincides with the domain of interest in astrophysical imagery, as the filters on Hubble and JWST telescopes range from .1 to 30 µm.

Attenuation spectrum

Considering that domain of high wavelengths (low frequency), the following assumption is then made for the spectral attenuation, assuming a linear model: τ spectrum (λ ) ≈ C T λ , with C T a constant1 . For lower wavelength, this approximation doesn't hold any more, as the shape of the function of the inverse of λ appears to become quadratic. This isn't a problem for the scope of this project which focuses on the simulation of commonly studied wavelengths in astronomy, which coincides with the range of "high wavelengths" where the linear approximation for the attenuation still remains accurate.

Finally, the total local attenuation can be obtained by plugging the different contributions in T local (t, λ ) ≈ e -τ space (t)τ spectrum (λ ) with τ space the formula of the spatial optical depth, in effect obtained by sampling the scene along a ray and thus let free here.

External light emission spectrum

In the following two sections (3.2.2, 3.2.3), we separate the modelling of external lights (stars and objects illuminating the scene from remote positions, implicitly contained in the L s term of the SVRE 2.3) from internal lights (the light emitted by the participating medium, for example the ionized gas part of an emission nebula, implicitly contained in the L e term of the SVRE).

We assume the standard hypothesis that any external body emits light following Planck's radiation law (3.1), thus notably depending on the temperature of the body.

I Planck,T (λ ) = 2hc 2 λ 5 1 e hc λ kT -1 (3.1)
We then distinguish two different models for these emission spectra, that can both be used in the following procedure:

• For low wavelengths in the studied spectral domain, Wien's approximation of Plank's law [START_REF] Wien | XXX. On the division of energy in the emission-spectrum of a black body[END_REF] is an accurate approximation for emitting bodies, expressed as the product of an inverse polynomial in λ and an exponential in e -A/λ (as written in equation 3.2). This approximation remains accurate for λ ≪ hc k B T .

I Wien,T (λ ) = 2hc 2 λ 5 e -hc λ k B T = A λ 5 e -B λ
With A and B two constants in λ (3.2)

• For higher wavelengths where Wien's approximation doesn't hold true any more, we assume that Plank's law can be projected on a piece-wise polynomial function of λ of arbitrary degree.

Internal light emission spectrum

Light emitted by the participating medium follows the allowed emission rays of ionization of the material making the dust cloud (most often hydrogen, helium, and oxygen [AG11; FT79]). The most accurate model for these rays is thus a set of Dirac distributions, only allowing those specific wavelengths to emit light. The integrals of these Dirac aren't always equal to 1, as they are to be pondered by the intensity of the considered rays. Based on observation of the sensors used in the main telescopes (3.3, 3.4), we assume that the filter spectra can be approximated as low-degree polynomials or piece-wise-defined polynomial functions. In the case of trapezoid-looking sensor functions, the model can for example be separated into three polynomials modelling the left edge, top edge, and right edge. For smoother sensor spectra, a similar approach is used with polynomials of higher degrees.

Sensor sensibility or filter spectrum

We stress here the importance of the model of the filter spectra, both as a way of getting an accurate fitting for the actual spectra function, but also with the goal to ensure an explorable functional space in the end. The optimal renderer that is built should be able to explore the space of all the potential physically accurate filter functions so that an astronomical scene could be seen through the "eyes" of multiple filters and compared in real time. Thus, it makes sense to let freedom in the choice of the spectra models (such as using a family of polynomials of arbitrary degree) in order to explore this space.

Phase function

The phase function is an implicit term in the expression of the incoming light L s in the SVRE. It describes the intensity at which incoming light will be scattered in every direction. The standard angular phase function is often modelled as Rayleigh scattering, which shape is illustrated in figure 3.5. In this work, we assume that our participating medium's incoming light phase function is isotropic (independent of the angle). This hypothesis is one of the strongest in the model and is forced by the proposed method in section 4, which doesn't keep information about the star positions, required for the computation of the phase function. However, this restriction could be relaxed with more work and an extended method: a way to relax this hypothesis and let the phase function free is discussed in the conclusion.

Participating medium coefficients

We assume that the coefficients σ s , σ a are separable in the form: σ i (t, λ ) = σ i,1 (t)σ i,2 (λ ), and the spectral part can be approximated as a polynomial in λ :

σ i (t, λ ) = σ i,1 (t)P σ i (λ ).
Full freedom is left on the expression of the spatial separable terms σ i,1 (t).

Assumptions summary

I s (λ ) I Wien,T (λ ) λ ≪ hc k B T Light emission (external) I s (λ ) P L (λ ) NA Light emission (internal) I e (λ ) ∑ δ (λ ) NA Phase function Φ(w) C Φ NA Medium coefficients σ i (t, λ ) σ i,1 (t)P σ i (λ ) NA Table 3.1: Assumptions summary 2 P x are polynomials, C x are constants 4
Real-time spectral integration

Analytical integration of the SVRE

Given the a-prioris discussed in part 3 and summed up in table 3.1, we show below that the SVRE 2.3 can be expressed as a product of polynomials in λ and exponentials of the form e -a(t)/λ . By only focusing on the spectral aspect (with t a position on the ray), the general equation F(t) = λ f (t, λ )dλ can now be expressed in a more compact way1 :

F(t) = λ P(λ )e -g(t)/λ dλ = N ∑ n=0 P n λ
λ n e -g(t)/λ dλ with:

P = P S P L P σ N = deg(P) g(t) = (τ(t) + τ s (t))C T
We prove in the appendix (6.1) that this form of the SVRE achieves to possess the integrability trait that was discussed in the previous section: All functions of this model can be analytically computed over the wavelength. Instead of relying on spectral discretization for the numerical integration, it is now possible to compute the exact value of the spectral integral (for the chosen model of functions). Moreover, this form for the SVRE allows the expression of the definite integral as a recurrence relation, which is relevant to the following implementation part.

Algorithm and implementation

Using the recurrence defined in 6.1, one bound of the spectral integral of f (t, λ ) can be computed in O(N) operations, with N the degree of the spectra product polynomial P. Algorithm Algorithm 1 Analytical integration of the spectral integrand function SPECTRALINTEGRATION(P, A, λ 0 , λ 1 ) ▷ P polynomial, A = g(t) return INTEGRALBOUND(P, A, λ 1 ) -INTEGRALBOUND(P, A, λ 0 ) end function

function INTEGRALBOUND(P, A, λ i ) x ← exp(-A/λ i ) α n ← -E 1 (A/λ i ) ▷ Exponential integral 1 β n ← λ i • x res ← COMPUTETERM(α n , β n , A, 0, P[0]) for n ← 1...∥P∥ do NEXTALPHABETA(λ i , n, A, α n , β n ) res ← res +COMPUTETERM(α n , β n , A, n, P[n]) end for return res end function function COMPUTETERM(α n , β n , A, n, c) return ( α n A n+1 + β n ) • c end function function NEXTALPHABETA(λ i , n, A, α n , β n ) α n ← -β n -A n • α n β n ← β n • λ i • n n+1
end function 1 shows a pseudo-code for the integration in O(N). Thus, the full spectro-spatial integral can be computed in O(N • M) with M the spatial resolution for a light ray in the scene. The naive rendering scheme for a nebulae scene would be to compute a spatial numerical integral for each ray shot by the camera, sampling the scene at regular intervals. For each of these samples, the spectral integral can then be computed analytically. This nullifies the need of relying on a spectral discretization for doing a double numerical integration, but isn't by itself an efficient enough solution to the original performance problem. Indeed, efficient numerical integration methods such as Gaussian quadrature integration will be able to obtain results as precise as those obtained with our analytical integration in the same complexity, for a similar polynomial fitted function of fixed degree.

Acceleration by projection

The formula obtained in section 4.1 for the spectral integrand can be transformed in a way that will reduce the operational load on run-time for the computation of the pixel intensity. We start by defining mathematical tools and spaces that will be used in order to rearrange the spectral integral terms in a more convenient manner.

With D ⊂ R + a domain, let V be the vector space generated by the family F = { f a | a ∈ R + }, with f a defined as:

f a : D -→ R λ -→ e -a/λ
F is also a linearly independent family, as no exponential of arbitrary argument can be expressed as a linear combination of exponentials of different arguments. As it is also spanning over V, F is a basis of V.

There is an analogy between the spectro-spatial integral that we want to solve analytically and the vector space V. The numerical integration of the spatial integrand for each ray shot in the scene can actually be expressed as an element of V, and thus decomposed as a linear combination of elements of F. This is shown in equation 4.1, where h(t j , w) can be seen as a linear coefficient for each element f g(t j ) of F, the full sum giving an element of V.

λ ∈D P(λ ) E t=0 h(t, w)e -g(t)/λ dtdλ ≈ λ ∈D P(λ ) N ∑ j=0 h(t j , w)e -g(t j )/λ ∆ jdλ = λ ∈D P(λ ) N ∑ j=0 h(t j , w) f g(t j ) (λ ) ∆ jdλ (4.1)
To compute the spectral integral of an element of V, the naive method is to keep the decomposition on the basis of F in order to compute each of these spectral integrals that are known analytically. However, in the most general case, there will be as much element of F that will take part in the sum as the spatial resolution of the model (in the case of equation 4.1, N terms). In other words, a spectral integral must be computed for every spatial step in the scene.

In order to reduce the number of computations needed to compute the spectral integral, we explore an approach to reduce the cardinality of the basis B used to decompose the spatial integral of V, while also keeping the cardinality of each element of the basis as low as possible. The second part is crucial in order to reduce the number of computations. The family F presented before has a maximal cardinality (card(F) = card(R + )), but a minimal element-wise cardinality: each element of F is made of a single component f a . On the other hand, the result of a Principal-Component Analysis would output the complete opposite of this basis: the PCA basis would have minimal cardinality (given the precision that is to be reached, PCA will output a few vectors containing most of the information of V). However, the PCA basis will have maximal element-wise cardinality, as each element of the PCA basis will be a linear combination of all elements of F. The idea here is to try to find a middle ground between the two bases.

Following on equation 4.1, the goal is to find the best possible basis B = {v i |i ∈ {0...N v }} so that the spatial integral can be projected and resolved in the fewest operations. This generic process is shown in equation 4.2, which shows a change of basis, from B = F to an arbitrary orthonormal basis that will be determined in the following part2 .

λ 2 λ =λ 1 E t=0 f (λ )e g(t)/λ h(t, w)dtdλ ≈ λ 2 λ =λ 1 f (λ ) N ∑ j=0 f g(t j ) (λ ) h(t j , w)∆ jdλ ≈ λ 2 λ =λ 1 f (λ ) N v ∑ i v i (λ ) N ∑ j h(t j , w) < f g(t j ) , v i > ∆ j dλ (4.2)
By choosing a good enough basis, we can make another approximation by ordering it by the magnitude and truncating the linear combination of the basis vectors to project each vector on the approximate subspace. This is similar to the approximation of only keeping the first N v vectors of a PCA basis in order to cover up to N p percentage of information on the vector space.

Basis determination

The approach defined in the previous part is only relevant and usable if the model chosen for the basis has a good enough convergence in terms of coverage on the vector space.

In the most generic case, with an infinite family of non-co-linear functions (as in our case), the precision of the projection over a restricted basis can greatly vary, as full coverage of the information can only be assured by a basis of infinite cardinality. In the following part, we discuss a way to obtain the most precise basis possible for our specific vector space V, and how this choice for the function models makes it possible to still get accurate results with a low cardinality (truncated) basis.

Function analysis

We quickly analyse the function shape of elements of F which will be projected on the basis. Knowing the range of elements it can reach is linked with the error rate between projections of elements. The function that is considered here is, for all {a ∈ R + }:

f a (x) = e -a/x : {x ∈ R + } -→ [0, 1]
Two relevant function limits are:

f a (x) a→0 --→ 1 and f a (x) a→∞ ---→ 0 (4.3)
A good basis should be able to model these two limits that are easily attainable in our model: a → 0 corresponds to an absence of density in the optical path, and a → ∞ corresponds to an extremely high density.

Scalar product and induced norm

In the following part, the usual scalar product between functions is used over the same spectral domain D as the definition domain of each f a ∈ F. The norm used to compute distances and error analysis between functions is the norm L 2 . Any notation ∥a∥ in the next sections of the report refers to the second norm ∥a∥ 2 , the most usual norm for functions. Moreover, as the goal is to ensure that the integral of the projection is similar to the integral of the element projected, it makes more sense to look at the norm L 2 (in effect an integration over the spectral domain) instead of other norms (like the maximum or mean of a function). Mentions of the distance between two elements are expressed as the norm of the difference between those elements.

∀a, b ∈ V < a, b > = D a(x)b(x)dx ∥a∥ 2 = < a, a > 1/2 d(a, b) = ∥a -b∥ (4.4)

Obtaining a precise basis

To ensure the generation of a precise basis which can be truncated and still allow for an accurate approximation of the projected elements, an iterative approach can be used in order to choose the best next possible element to add to the basis. This constraint can be seen as choosing the element that will minimize the norm of the difference between an expected element of F and its projection, and this for every element of F. This can be expressed as minimizing the following equation by choosing the best possible basis B:

∑ f ∈F ∥ f -proj B ( f )∥ (4.5)

Application

This problem can generally be solved following a gradient descent to minimize the previous function and obtain the best possible basis. In the following application, this is approximated by a greedy approach in O(card(B) • card(F c )), with F c a discretized subset of the infinite family F. The basis B is constructed element by element, choosing the best possible element b at each step that will minimize

∑ f ∈F ∥ f -proj B ∪ b ( f )∥.
In the following part, the spectral domain is fixed to D = [0 µm, 1.0 µm]. A typical astrophysical filter has a range length between 0.1 and 1 µm, and the minimal wavelengths observed by Hubble and JWST are around .2 µm while the maximal values are around 30 µm (see figures 3.3 and 3.4). The "hardest" part to define our projection is around low wavelengths, as it is where small changes in the optical path will have greater repercussions on the spectral function. This is why for generality, the following application will focus on this "hard" region, while more effective results for higher wavelength are discussed in the subsection 4.2.2.

With this fixed spectral domain, results are shown below with a first step where the basis is initially generated with the spectral function associated with the minimal optical path value of our model (a = 0). This corresponds to the first limit of the spectral function shown in equation 4.3, f a a→0 --→ 1. This ensures that the lowest densities can be modelled by our basis, as well as the extremely high densities (the projection of f a (x) a→∞ ---→ 0 over a basis containing 1 can trivially be 0 by nullifying the projection coefficient). The graph in figure 4.1 shows the error function for four different bases built on arbitrary candidates for the second element (b 1 = 1, b 2 = {e -.32/λ , e -.78/λ , e -1.0/λ , e -1.3/λ }). The y-axis ranges from a distance of 0 (meaning that an element and its projection are the same considering our distance metric) to 1 (meaning that an element and its projection are the farthest they can be (here d(1, 0) = D (1 -0) 2 dx = 1)).

The best candidate for the basis is the one that will minimize the integral of the error function. A continuous plot of the areas of the error functions is thus shown in the first graph of 4.2. The minimization step for a second basis element gives a spectral function corresponding to an optical path of 0.45 (b 2 = e -.45/x ), as illustrated by the minimum in the function. The error function between all elements of F and their projection on this theoretically optimal basis B = {1, e -.45/x } is then plotted on the second figure of 4.2. As expected, there are two local minima at 0 on the two optical path values of our basis (0 and .45), and another one at the limit a → ∞ where the function converges to 0. These correspond to points that are fully expressible on this basis, while the non-zero values of the curve correspond to elements that will lose some of their information after projection. There are also two local maxima, one at low densities (≈ .1), and one at medium densities (≈ 1). This means that elements sampled in the 3D scene with an optical length close to these local maxima will lose part of their information when projected on the selected basis. On the other hand, elements with high optical length will lose almost no information, meaning that their spectral integral can be fully computed by their projection on B.

It can be noted that while the criterion used to decide the next best basis element 4.5 is an "objective" criterion, as it achieves to get the basis that will get a minimal distance between a plane and its projection, it is also possible to tweak it and focus on specific areas of the argument space. The second graph of figure 4.2 shows that the optimal 2-element basis has a local The argument for the second element of the basis being constructed that minimizes the error between all elements and their projection is obtained by minimizing the curve on graph 1, for a value of (a = .45).

maximum for elements of a very low argument. To better fit a basis in a scene containing a lot of elements with that low arguments, it would be possible to tweak the second basis element towards this local maxima. This would mean that the entire sum of the error function would increase, but that the low densities would be better modelled. We chose here to stay with an "optimal" basis without tweaking it to stay general.

At this second minimization step, we thus get a two-element basis where the maximum error between any element and its projection is at ≈ 0.015. The process can then be iterated multiple times from this extended basis B = {1, e -.45/x } until the expected precision is obtained, and when the projections get precise enough in comparison with the actual function vector.

Figure 4.3 illustrates the decrease of the error when the basis is extended, which follows a logarithmic decrease. For each element that is added to the basis, the error is divided by 10. It is to note that this goes in pair with a linear increase in the run-time complexity of the algorithm, as it will be discussed in the implementation part. Each element of the basis indeed needs to be integrated analytically once on the sensor.

Higher spectra domain -Telescope wavelengths

Considering the wavelengths of interest for Hubble and JWST, the wavelength ranges are around 0.1 µm to 30 µm, and the length of these ranges often does not exceed 1 micron. With these parameters, we get significantly better results than when using the previous wavelength range, as the most difficult part to integrate of the model function is located near low wavelengths. Likewise, a spectral domain spanning over a smaller range of wavelengths will enable a basis of similar cardinality to generate better results. Examples are illustrated in the figure 4.4, for a spectral domain of D = [1 µm, 10 µm]. The error is divided by 100 for each component added to the basis, instead of 10 for a range of lower wavelengths. The same precision can be attained in half the number of basis component f a with higher wavelengths, even though the The error between all elements and their projection is divided by 100 for each element added to the basis. This is 10 times better than the results obtained for lower wavelengths.

Implementation and error analysis

The SVRE is computed over a scene containing spectral information (that was modelled in the previous steps), but also information about the local density at each point of the scene (used for example to compute the optical depth). The simulation of an accurate density model is outside this project's scope and has been partially covered in previous works as part of the underlying ANR [Anr]. The density structure for our comparison test scenes were built using Fractional Brownian Motion (fBm) applied to Perlin noise, but could be replaced by any arbitrary model for a nebulae density. We chose the combination of fBm and Perlin noise, as by tweaking the set of parameters (amplitude, frequency, octaves) we can obtain multiple kinds of spatial density distribution to test our model. This is important, as the previous section showed that the results can vary based on the distribution of the density in the scene, with the higher relative errors located around low densities for small-cardinality basis-based projection.

Basis generation

Algorithm 2 Basis generation function GENERATEBASIS(N, a 0 , a 1 )

for i ← 0...N do b ← FINDBESTMONOME(B, a 0 , a 1 )▷ Gradient descent, find best element to add to B b ← GRAMSCHMIDT(B, b) ▷ Generate orthonormal vector to B B ← B ∪ b end for return B end function
Algorithm 2 illustrates a pseudo-code for the generation of an accurate basis, according to the steps discussed in section 4.2.1. This function is to be called only once at the start of the program, and the generated basis can then be passed to every other rendering procedure without needing to be re-computed at every frame. This step is relatively costly, as the implicit findBestMonome function iterates over all the possibles elements from the family F in order to select the element that, once added to the basis, will minimize the sum of errors for all elements in the vector family. The complexity of this part is of O(card(F c ) • N B ). N b is the cardinality of the basis that is to be generated, and card(F c ) is the cardinality of a discretized subset of the infinite family F. The greater the cardinality of B will be, the more accurate the results of the projections will be in the following steps, but the computational load will be linked with it. On the other hand, the cardinality of the subset F only contributes to this basis generation step. There is no run-time disadvantage to have as much spatial resolution in this discretized family for this basis generation step, as it will ensure greater precision on the basis generation.

Pre-computations

One benefit of the projection method is that we can pre-compute the shadows (external lights optical depth) and sample them in constant time during execution. Algorithm 3 illustrates the pre-computation of the shadows, while algorithm 4 shows how to use these values on run-time in order to compute the projected spatial-spectral integral. In practice, the contributions of each ▷ p is a pixel on the camera array r ← GETRAY(p) for j ← 0...N do ▷ Spatial numerical integration p ← GETPOINTALONGRAY(r, ▷ Get argument of basis element i c ← PROJECTIONCOEF(i, b, p) ▷ Compute projection coefficient for the ith term v ← v + c • SPECTRALINTEGRATION(P, a, λ 0 , λ 1 ) end for STORE(v, p) ▷ Storing intensity value in pixel p end function element of the discretized space (in 3D, a voxel of fixed length) are computed and stored in a data structure (in GLSL, a texture3D) that can be sampled at run-time. The size of this texture needs only to be the same as the input density texture, containing information about the local density number for each voxel of the discretized space.

Both algorithms presented are designed to fully take into account the benefits of GPU programming, by being completely parallelizable for each voxel in the scene. It is to note that a choice of implementation for the load function in algorithm 4 must be made. The approximation that the texture filled in the pre-computation step can be sampled continuously (and each value computed with trilinear interpolation) must be made in case one wants to use the full potential of most GPU implementations (the trilinear interpolation being optimized for most GPU cards). In the other case, this load function must be implemented to compute or approximate the value at any point in the scene or restrain the ray to follow the voxels of the discretized scene.

Error analysis

By comparing the results with a ground truth computed with a double numerical integration (spatial and spectral), the error stays below 6% for any pixel on a computed scene as shown in figure 4.5 for a projection method with an optimal basis of cardinality only 3. This means that only three spectral integrations are done for the projection method (with the goal of retrieving most of the information of the function vectors that were projected on the basis), while the full integration method would have had to do as many spectral integrations as the spatial resolution. We achieve to model in a visually correct way both the low and high densities in our nebulae scenes, with no significant color artefacts or discontinuities.

Complexity

We observe in effect a great increase in performance when using the projection method instead of the numerical one. The pre-computation is solved in a fixed time and helps reduce the operational load of the run-time integration. Thus, only card(B) analytical integrations are done at run-time. Part of the time complexity at run-time is converted to space complexity. As discussed previously, we fill a voxel texture with the pre-computations, in effect taking N space space for each ray. Table 4.1 references the time complexities for each method (separating the pre-computations and the run-time integration for the projection method). Table 4.2 illustrates the performance results when computing a particular scene and comparing the methods.

The projection scheme used in this section doesn't rely on the expression of the spectra polynomial (product of the emission and filter spectra models), which can thus remain free and modified freely at run-time. In effect, one can explore any scene in real time and slide or modify the filter functions used for the spectral integration, as it would be the case when using a real telescope.

Method

Multiple emission spectra in the scene: Light-Density projection

A typical nebulae scene can contain a huge number of light sources, as the dust cloud can cover multiple solar systems in its span, and thus multiple stars. The default hypothesis is to state that each of these stars can have a different emission spectrum, and thus interact differently with the participating medium.

In the naive integration scheme, the full rendering pipeline needs to be computed once for each emission spectrum (and thus each star), as the polynomial in terms of lambda will differ for each light source.

We discuss here a solution to address this problem, making it possible to compute a scene with multiple light sources and multiple spectra at a reduced cost.

Bi-dimensional projection

Section 4.2 presents a way to generate a basis of functions for expressing the optical paths of the elements of volume.

We can extend this method to a multidimensional basis, making the basis dependent on both the optical path and the spectrum of the light. Our light spectrum model is taken from the black body equations, which reduces the dimensionality of the possible emission spectra to one (the temperature of the star directly infers the emission spectra). Thus, the basis to be generated will depend on two variables: the optical length and the temperature of the star. Function 4.6 shows analogous vector spaces to the one presented in part 4.2.

With D ⊂ R a domain, let V be the vector space generated by the family F = { f a,T |(a, T ) ∈ (R + ) 2 }, with f a,T defined as:

f a,T : D -→ R λ -→ P T (λ )e -a/λ (4.6)
In effect, the generation of the basis follows a similar approach as the one described in the previous section. The only difference is that the variable space for the gradient descent is two-dimensional, but all the operators remain the same.

Application

Algorithm 5 Pre-computation of projections

function PRECOMPUTE(v, b) ▷ v is a voxel of the discretized 3D space, b a basis for L i ← 0...∥lights∥ do l ← OPTICALPATHTOLIGHTSOURCE(v, lights[L i ]) p ← p + PROJECT(l, lights[L i ], b) ▷ Projection on basis b end for STORE(v, p) ▷ Storing projection in voxel v end function
For verification purposes, the emission polynomials for the family F are generated as linear interpolations between two "edge" polynomials: P 0 (x) = (x+1)(x) 2 and P 1 (x) = (x-2)(x-1) 2 , P T (x) = (1 -T )P 0 (x) + T P 1 (x) with T ∈ [0, 1]. This range covers most of the possible aspects for emission polynomials, with functions associated with low temperature emitting most photons in the high wavelengths, and functions associated with high temperature emitting most photons in the low wavelengths. For the rest of the parameters, we follow the same decisions as in the previous section 4.2, with the extension for the two-dimensional basis when needed.

For the application, a basis with cardinality 4 is chosen in order to describe the full scene. The same algorithm as described in the previous section 4.2.2 to generate the best possible basis is applied, and outputs four basis vectors to project each element of the scene. The complexity of this first step is significantly greater than for a basis only relying on the optical path, as the function to minimize is now a 2D surface and no longer a 1D function. Still, a gradient descent approach can be applied and will give an optimal basis according to our criteria.

Figure 4.6 illustrates the results that are obtained when computing a scene with 6 lights of different spectrums. In effect, all of these light contributions are being projected onto the four basis elements which are spectrally integrated at the end of the pipeline. The algorithm applied is extremely similar to both algorithms 3 and 4 presented in the previous section, the only difference being the nature of the basis and the actual implementation of the project method, which now takes into account the spectrum of the emitting light sources.

The results obtained appear decent, with low visual differences between the full double integration method and the projection one. We note that for a basis of similar length (3 and 4), we get worse results with this spatial-spectral basis in comparison with the previous one. This can be explained because the basis is now covering a two-dimensional space, and information about both the optical depth and the spectral polynomial needs to be projected. It is still possible to decrease this error by augmenting the number of elements of the basis, but this also goes along with an augmentation of the run-time computational load.

Complexity

Table 4.3 illustrates the complexity of the algorithm in comparison with a double integration for a variable number of lights (L lights ). The space complexity is the same as in the previous scheme. In effect, while the pre-integration step takes more time as the projection of each light contribution must be computed, the run-time integration isn't impacted at all. This means that the scene can include any number of lights and spectra, and the complexity will remain the same.

Acceleration by pre-integration

In computer rendering, pre-integration [EKE01] is a method used to precompute parts of the results of the rendering equation (in this case the SVRE) before runtime, and to use these precomputed results during execution to reduce the real-time computational load by only having to sum small elements during the spatial integration part.

While the generic methods for pre-integration cannot be applied in the spectral case where the intervals have been inverted, it seems tempting to add a few approximations in our model to be able to apply the same pre-computations and reduce the operational load. A few approximations were made in the following part in order to explore this possibility, but no good results were reached with them in comparison with the projection 4.2 method. The previous method is thus to be preferred, as it gets better results in fewer computations and better scalability. We still chose to present this approximation below, as it seems to be the more "naive" way to approximate the literature pre-integration.

As in the previous part, we start from the spectral integral form of the SVRE 2.3, and we make several approximations in order to express it under a convenient form, as shown in equation 4.8:

λ 2 λ =λ 1 f (g(t), λ ) dλ = N ∑ k=0 P k λ 2 λ =λ 1 f k (g(t), λ )dλ (4.7) N ∑ k=0 P k λ 2 λ =λ 1 f k (g(t), λ )dλ (4.8a) ≈ N ∑ k=0 P k C k (λ 1 , λ 2 ) exp (-[g(t)] • ε(λ 1 , λ 2 )) (4.8b) = N ∑ k=0 P k C k (λ 1 , λ 2 ) exp (-[τ(t) + τ s (t)] • ε(λ 1 , λ 2 )) (4.8c) ≈ N ∑ k=0 P k C k (λ 1 , λ 2 )• (4.8d) exp -[ t ∑ j=0 σ t ( j)∆ j + τ s (t)] • ε(λ 1 , λ 2 ) (4.8e) = N ∑ k=0 P k C k (λ 1 , λ 2 ) exp (-[σ t (0)∆t] • ε(λ 1 , λ 2 )) 1 • (4.8f) t ∏ i=1 exp -σ t (i)∆ j • ε(λ 1 , λ 2 ) 2 exp(-τ s (t) • ε(λ 1 , λ 2 )) 3 (4.8g) with: C k (λ 1 , λ 2 ) = λ k 2 -λ k 1 k (4.9)
We discuss below the two approximations made in demonstration 4.8:

• 4.8b : The spectral integral is approximated as an exponential in g(t)ε(λ 1 , λ 2 , k) and pondered by a term C k (λ 1 , λ 2 ). This approximation makes sense by looking at arbitrary analytical results of the spectral integral, assuming that the last term in E 1 (g(t)/λ ) is negligible. The remaining terms can thus be put in this form, two significant results being that the two limits of the function (in g(t) -→ 0 and g(t) -→ ∞) are respected by this approximation.

• 4.8e: Discretization of the analytical spatial integration from numerical spatial integration (necessary in all cases for the methods presented in this work, in order to preserve full freedom for the spatial density function).

The final term 4.8g is expressed in a very convenient way for our problem. The first part of the term 4.8g (red) is exactly the result of the spectral integral over a small element of the domain where the spatial σ s coefficient is uniform (i.e. the first voxel from the camera). The second term (green) is a product of all spectral integrals between the start of the ray and the current point, without external illumination, and divided by the constant term ∑ N k=0 P k C k (λ 1 , λ 2 ). Finally, the third term (blue) is only the external light contribution at point t.

Knowing the exact function ε isn't relevant in this case, as it will be computed by the analytical resolution of the spectral integral. The purpose of setting this function is simply to become able to separate the exponential in a product.

Implementation and error analysis

This way of writing the spectral integral makes it straightforward to compute in a front-to-back algorithm using two precomputed values for each voxel, as depicted in algorithms 6 and 7:

Algorithm 6 Pre-integration approximation function PRECOMPUTE(v) ▷ v is a voxel of the discretized 3D space ρ ← LOCALDENSITY(v) p ← SPECTRALINTEGRATION(P, ρ, λ 0 , λ 1 ) STORE(v, p) ▷ Storing spectral integral in voxel v end function
The results obtained with this method do not really fall into the criteria that were decided at the beginning of the report (section 2). While the renders are indeed fast and explorable to some extent (full exploration of the space is possible as the pre-integration only computes the contributions for each small voxelized element, which are to be summed up at run-time), they aren't close enough to the ground truth to appear "visually" accurate, especially in comparison with the previous projection method (4.2). Example results are shown in the appendix (6.3). This distance between the ground truth and the actual renders from this approximation should come from the assumption that one term of the spectral integral (E 1 ) can be nullified in front of the other contributions. This should be an assumption that is too strong for accurate renders and doesn't retrieve all of the information that would make for a relevant render. Moreover, there is no obvious way to increase the accuracy of this scheme and get better results ▷ p is a pixel on the camera array

r ← GETRAY(p) p ← GETPOINTALONGRAY(r, j) C k ← COMPUTEPKCK(P, λ 1 , λ 2 ) prod ← 1 v ← 0 for j ← 0...N do ▷ Spatial numerical integration p ← GETPOINTALONGRAY(r, j) v local ← LOAD(p)
▷ Loading the spectral integral value at position p

if j! = 0 then v local ← v local /C k end if v ← v + σ s • prod • ∆ j prod ← prod • v local end for STORE(p, v)
▷ Storing intensity value v in pixel p end function (whereas for the projection scheme, increasing the cardinality of the basis will significantly increase the results). Thus, the previous method should be preferred.

Conclusion & future works

In this report, we presented an efficient technique to compute a physically accurate spectral volumetric render for astronomical scenes. We achieve to compute in real-time complex scenes taking multiple seconds per frame with the ground truth method, without significant visual changes, reaching an improvement of 83% on the double numerical integration performances. This technique makes it possible to dynamically explore a nebulae scene in real-time and dynamically change the filters used by the simulated telescope in order to retrieve multiple kinds of spectral information, mimicking the classical imagery done with telescopes such as Hubble and JWST. This technique is also easily scalable, as increasing the spectral precision relies on augmenting the cardinality of the considered basis, and increasing the spatial resolution relies on the resolution of the underlying voxel array.

While this basis integration method is general and doesn't necessarily rely on the spectral model that was presented in section 3, further research could try to generalize the model for any kind of scene (and not only astrophysical ones). The model for the attenuation in astronomical dust clouds seems to be linearly correlated with the inverse of the wavelength, but this relationship would no longer be true when considering any arbitrary volumetric objects. Different models would have to be chosen in order to apply a similar integration method. One of the "hard" assumptions made by our model in section 3 is also that the phase function for the participating medium must be isotropic, as in any case the information on the position of the emitting star would be lost during the projection part of the spectral pre-computations. However, it could be possible to relax this hypothesis by using another form of projection to keep the angle-dependant intensity information, for example on the basis of 3D spherical harmonics.

Furthermore, one area of improvement for our model would be to take into account multiple scattering. The integration scheme currently only models one light bounce, but the fact that the model presented in section 4.3 allows the projection of any number of light sources in the scene at constant cost appears to be a good indication for a potential extension of the model. It should be possible to compute all of the light bounces in a first pass as secondary lights and project them on the same bi-dimensional basis. While the cost for the pre-computation would greatly increase in order to compute all the possible N-bounce interactions, the cost at run-time would remain constant (and in our case, real-time).

To conclude, with precise a-priori knowledge of physics and a relevant model for complex phenomenons, it is possible to render in great detail, and in a few milliseconds, complex astronomical scenes which would otherwise be impossible to compute in real-time.

Addendum

Analytical integration

The following equations demonstrate a relevant recurrence relation for the computation of the analytical antiderivative of the product of a polynomial and an exponential (the applied basis function defined in section 4).

The only remaining term left to be computed is the E 1 1 function, which is only sampled twice per definite integration of the whole basis function set due to the recurrence relationship. It can be tabulated or approximated (for example with a Chebyshev approximation [START_REF] Cody | Chebyshev Approximations for the Exponential Integral Ei(x)[END_REF]).

λ n e -g(t)/λ dλ = λ n+1 n + 1 e -g(t)/λ + g(t) n + 1 e -g(t)u u n+1 du With u = 1 λ = β n (g(t), λ ) + g(t) n + 1 α n (g(t), λ ) (6.1)

α n (x, λ ) = -β n-1 (x, λ ) -x n α n-1 (x, λ ) α 0 (x, λ ) = -E 1 (x/λ ) (6.2)

β n (x, λ ) = λ n n + 1 β n-1 (x, λ )
β 0 (x, λ ) = λ e -x/λ (6.3)

Negative powers

For the sake of generality (allowing the numerical integral of models integrating Wien's approximation for emission spectra 3.2.2), we show below the antiderivative for negative powers of lambda.

λ -1 e -g(t)/λ dλ = E 1 (g(t)/λ )

1 Exponential integral function E 1 λ -2 e -g(t)/λ dλ = e -g(t)/λ g(t)

∀k ≥ 3 λ -k e -g(t)/λ dλ =u k-2 e -g(t)u du With u = 1 λ = -u k-2 e -g(t)u g(t) + k -2 g(t) u k-3 e -g(t)u du

We can then use a similar recursive relationship as in the previous section 6.1 to compute the remaining antiderivative of the form u k e -au du. .1 illustrates how the projection method presented in section 4 compares to other numerical methods. Two Newton-Cotes methods and Gaussian quadrature converge slower than the projection method for most functions that were observed (excluding degenerate cases for the polynomial coefficients and the exponential arguments). We believe these results not to be accurate enough in comparison with the double numerical integration, especially when compared with the projection method discussed in section 4.2. While the color of the main part of the nebulae appears similar for both methods, shadows, and external light attenuation doesn't appear clearly in the pre-integrated results, and multiple discrepancies appear in the low-density areas of the nebula cloud.

Comparison of convergences for projection

Results of the pre-integration method

Figure 3

 3 Figure 3.1: (a) Absorption nebula (black object) (b) Reflection nebula (blue object) (c) Emission nebula Image courtesy: nasa.gov

Figure 3

 3 Figure 3.2: Data courtesy: Cardelli & al 1989

Figure 3 . 3 :

 33 Figure 3.3: Some Hubble filters (61 in total) Image courtesy: HST documentation

Figure 3

 3 Figure 3.5: Angular aspect of Rayleigh scattering for an incoming ray coming from the left

Figure 4

 4 Figure 4.1: Visualization of various error functions for arbitrary choices of the second basis member The second basis element has a form: e -ARG_2/x

Figure 4 . 2 :

 42 Figure 4.2: Visualization of the first minimization stepThe argument for the second element of the basis being constructed that minimizes the error between all elements and their projection is obtained by minimizing the curve on graph 1, for a value of (a = .45).

Figure 4 . 3 :

 43 Figure 4.3: Evolution of error along basis extension The error between all elements and their projection is divided by 10 for each element added to the basis. For the first plot, the elements of the basis B are represented in their argument form (B[i] = e -b i /x )

Figure 4 . 4 :

 44 Figure 4.4: Increased projection precision when considering higher wavelength rangesThe error between all elements and their projection is divided by 100 for each element added to the basis. This is 10 times better than the results obtained for lower wavelengths.

  j) l ← LOAD(p) ▷ Loading light optical path at position p x ← PROJECT(ρ + l, b) ▷ Projecting the optical path ρ + l on basis b e ← PROPAGATE(p, j, x, e) ▷ Propagate the projected contribution on the basis through the light ray end for for i ← 0...|b| do ▷ Analytical integration of each member of b a ← ARG(b[i])

Figure 4 . 5 :

 45 Figure 4.5: Comparison between numerically integrated renders and renders obtained by a projection on a 3-element basis. (a) Double numerical integration (b) Projection approximation (c) Absolute error ×16

Figure 4 . 6 :

 46 Figure 4.6: Comparison between numerically integrated renders and renders obtained by spectral-spatial projection on a 4-element basis (a) Double numerical integral (b) Projection approximation (c) Absolute error ×6

Figure 6 . 1 :

 61 Figure 6.1: Comparison of convergences for projection and numerical methods on the integration of element f (x) = ((x -1)(x -.2)(x -.5)(x -.7)(x)100 + 1)e -.1/x

Figure 6

 6 Figure 6.2: Comparison between numerically integrated renders and renders obtained by pre-integration (a) Double numerical integral (b) Pre-integrated results approximation (c) Absolute error ×4

Table 3 .

 3 1 sums up the assumptions discussed in the previous parts.

	Name	Symbol	Model 2	Domain restriction
	Sensor sensibility, Filter	S(λ )	P s (λ )	NA
	Attenuation	T(t, λ )	e	-τ(t)C T λ	λ > 0.25 µm
	Light emission (external)

  Algorithm 3 Pre-computation function PRECOMPUTE(v, b) ▷ v is a voxel of the discretized 3D space ρ ← LOCALDENSITY(v) for L i ← 0...∥lights∥ do l ← l..OPTICALPATHTOLIGHTSOURCE(v, lights[L i ])

	end for
	STORE(v, l)	▷ Storing optical depth l in voxel v
	end function
	Algorithm 4 Run time integration
	global variables
	P, Spectral polynomial
	λ 0 , λ 1 , spectral integration bounds
	B, basis
	end global variables
	function COMPUTEINTENSITY(p)

  Complexity Double Netwon integration O(N space • (M wavelength + N space )) Naive analytical spectral integration O(N space • (Deg(P) + N space )) Pre-process O(N space • N space ) Run-time integration O(N space •Card basis + Deg(P) •Card basis )

	Table 4.1: Comparison of algorithm complexity between numerical integration of the SVRE
	and projection approximation
	Method	Avg Time per frame 3
	Newton	260ms
	Naive analytical computation	160ms
	Run-time projection	30ms
	Table 4.2: Comparison of performance (time to compute each frame) between numerical inte-
	gration of the SVRE and projection approximation on scene A (figure 4.5) for a similar preci-
	sion level

  Table 4.4 illustrates the program performances for computations on scene B (4.6). space • (M wavelength + N space ) • L lights ) Projection (pre-process) O(N space • (Deg(P) •Card basis + N space ) • L lights ) Run-time integration O(N space •Card basis + Deg(P) •Card basis )

	Method	Complexity
	Double Netwon integration O(N Table 4.3: Comparison of algorithm complexity between numerical integration of the SVRE
	and spectral-spatial projection approximation
	The run-time integration performance of the SVRE isn't impacted by the number of lights in
	the scene
	Method	Avg Time per frame 4
	Newton	2500ms
	Run-time projection	30ms
	Table 4.4: Comparison of performance (time to compute each frame) between numerical
	integration of the SVRE and projection approximation for scene B (figure 4.6).

For planar cameras, the position x varies on the camera plane, and w is the direction from the focal point to the pixel position.

In section 4, C T will be fixed to C T = 0.56 to model an average of the dust clouds attenuation plotted on figure 3.2.

For the sake of simplicity, we consider the case where the emission spectrum is modelled as a polynomial in λ , and not according to Wien's approximation as presented in section 3.2.2. Choosing to do so would only require extending the polynomial P to take into account the potential five inverse powers of lambdas introduced by Wien's approximation. This is completely doable (and the model is still integrable in a recurrence as presented in appendix 6.1), but it introduces over-complexification that is not needed in order to understand the gist of the model and integration scheme.

< f g(t j ) , v i > is the scalar product of e g(t j )/λ and v i .

N=100,M=20,C b =3,Deg(P)=5,GPU=Nvidia GTX 1080 Ti, scene A

L lights = 6, N=100,M=20,C b =4,Deg(P)=5, GPU= Nvidia GTX 1080 Ti, scene B