Liangsheng Ge Beibei Wang

Nicolas Holzschuch

Precomputed Multiple Scattering for Rapid Light Simulation in Participating Media

Keywords: Participating Media, Precomputation, Multiple Scattering. !

Rendering translucent materials is costly: light transport algorithms need to simulate a large number of scattering events inside the material before reaching convergence. The cost is especially high for materials with a large albedo or a small mean-free-path, where higher-order scattering effects dominate. We present a new method for fast computation of global illumination with participating media. Our method uses precomputed multiple scattering effects, stored in two compact tables. These precomputed multiple scattering tables are easy to integrate with any illumination simulation algorithm. We give examples for virtual ray lights (VRL), photon mapping with beams and paths (UPBP), Metropolis Light Transport with Manifold Exploration (MEMLT). The original algorithms are in charge of low-order scattering, combined with multiple scattering computed using our table. Our results show significant improvements in convergence speed and memory costs, with negligible impact on accuracy.

Introduction

Computing illumination simulation in scenes with participating media is still a costly process, even with algorithms specifically designed to handle them. To render high-albedo materials, such as milk or skin, we need to simulate a large number of scattering events before convergence. Worse, the overall appearance of these materials is often very smooth, almost uniform, meaning we used a lot of computational power for an almost constant appearance.

In recent years, several algorithms have been introduced for rendering participating media, such as Virtual Ray Lights (VRL) [START_REF] Novák | Progressive virtual beam lights[END_REF], several extensions to photon mapping culminating with UPBP (Unified Points, Beams and Paths) [START_REF] Křivánek | Unifying points, beams, and paths in volumetric light transport simulation[END_REF] and Manifold Exploration Metropolis Light Transport [START_REF] Jakob | Manifold exploration: A markov chain monte carlo technique for rendering scenes with difficult specular transport[END_REF]. All these methods have greatly improved simulation of participating media. They still encounter problems for simulation of materials with high albedo where multiple scattering dominates. Rendering these materials degrades the quality of the simulation for the entire scene (see Figure 1): a lot of CPU power is used for the translucent material, leaving little for the rest of the scene. This explains the large amount of noise in Figure 1(b).

Dipole methods [START_REF] Jensen | A practical model for subsurface light transport[END_REF], [START_REF] Jensen | A rapid hierarchical rendering technique for translucent materials[END_REF] allow fast simulation of translucent materials by approximating surface response and precomputing incoming radiance. This can be combined with single scattering computed from existing simulation algorithms. Complete integration with existing algorithms is difficult because dipole methods usually rely on precomputing incoming radiance on the surface of the object.

We present a method for faster simulation of light transport in translucent material. Our method is easy to integrate with existing illumination simulations algorithms. We provide examples for path tracing (VRL), photon-mapping (UPBP) and Metropolis Light Transport (MEMLT). Our algorithm uses precomputed tables for multiple scattering, and relies on the original algorithm for single and double scattering. More specifically, we shoot photons from a point light source and a segment source in an infinite medium and we accumulate radiance for samples distributed in two spatial dimension and two angular dimension, or two 4D tables (one for point light source and one for segment source). We then approximate the lobes in these two 4D tables using axial symmetry for point sources and using asymmetric lobes for segment sources, resulting in compact representations of multiple scattering for homogeneous medium. These representations can be used for fast multiple scattering computation in many rendering algorithms by connecting the point / segment from light side and point / segment from camera side directly. Our method results in faster convergence for all algorithms, greatly reducing noise in the pictures with minimal impact on accuracy. Reducing the computational cost for participating media makes more computational power available for the rest of the scene, improving image quality overall.

In the next Section, we review some of the previous work on rendering participating media. In Section 3, we present our algorithm for precomputing multiple scattering and use it in existing algorithms. We explain implementation details in Section 4. We present our results, compare with previous work and analyze performances in Section 5, and conclude in Section 6.

Previous Works

Photon Mapping: Chandrasekhar [START_REF] Chandrasekhar | Radiative transfer[END_REF] introduced the radiative transfer equation, describing radiation transport in participating media. Jensen and Christensen [START_REF] Jensen | Efficient simulation of light transport in scenes with participating media using photon maps[END_REF] apply this equation to light transport, and presented an algorithm based on Photon Mapping. Jarosz et al. [START_REF] Jarosz | The beam radiance estimate for volumetric photon mapping[END_REF], [START_REF] Jarosz | A comprehensive theory of volumetric radiance estimation using photon points and beams[END_REF], [START_REF] Jarosz | Progressive photon beams[END_REF] extend this algorithm using beams instead of photons for faster computations, with less noise. Křivánek et al. [START_REF] Křivánek | Unifying points, beams, and paths in volumetric light transport simulation[END_REF] automatically selects between beams, points and paths in participating media using multiple importance sampling. Bitterli and Jarosz [START_REF] Bitterli | Beyond points and beams: Higher-dimensional photon samples for volumetric light transport[END_REF] further extend the idea by tracing photon planes and volume.

Lighting Stage

Volume Event Surface Event

Rendering Stage

Precomputation Table

Precomputation Stage

Lobe Gather Contribution

Point Source Segment Source

Main Axis Main Axis Ratios

Fig. 2. Our algorithm: we begin by precomputing two multiple scattering tables. We shoot particles from the origin of the point source (red dot) or point samples (red dots) from the segment source. The photons propagate in the imagining medium, being absorbed or scattered. For scattering event, we update the contribution of the photon to each sampling direction to the table. After the shooting, we get two 4 dimension tables. Then we approximate the lobes in these two tables with different features of the lobe: using axially symmetric lobes for point source and with asymmetric lobes for segment source.

In the end, we get two much more compact tables. This step is required only once per material. We use these tables for multiple scattering in existing simulation algorithms (here, VRL). In the lighting stage of VRL, we shoot, scatter rays, and store the rays (orange line) with first scatter events. In the rendering stage of VRL, the precomputation tables are used to compute the multiple scattering contribution from the light rays (orange line) to the camera rays (red line).

Photon-mapping-based methods can provide high quality simulations of light transport in participating media, but they usually require a long time to converge. In contrast, our method can provide similar quality with much faster convergence for highly scattering media, as we connect the light point / segment and camera point / segment through the tables with already converged radiance.

Diffusion Theory: Jensen et al. [START_REF] Jensen | A practical model for subsurface light transport[END_REF] approximate surface response for high-albedo materials using the Dipole method. Jensen and Buhler [START_REF] Jensen | A rapid hierarchical rendering technique for translucent materials[END_REF] present a fast method for computing dipole response, relying on precomputed incoming radiance at the surface. D'Eon and Irving [START_REF] D'eon | A quantized-diffusion model for rendering translucent materials[END_REF] improved the accuracy using Quantized Diffusion method. Habel et al. [START_REF] Habel | Photon beam diffusion: A hybrid monte carlo method for subsurface scattering[END_REF] combined photon beams and the diffusion model by distributing extended virtual sources along the photon beams. Frisvad et al. [START_REF] Frisvad | Directional dipole model for subsurface scattering[END_REF] introduced the Directional Dipole, taking into account the orientation of incoming light relative the surface.

Diffusion-based methods are fast, and designed to work with materials with a high albedo, like our own algorithm. Diffusionbased methods tend to be faster than our method to compute multiple scattering effects. In comparison, our method can handle highly anisotropic media, while diffusion theory based methods use similarity theory for anisotropic media, which leads to inaccuracy. Even for isotropic media, our method provides better quality, by tracing the light inside the medium and combining with VRL, UPBP or MEMLT.

Virtual Point Lights: Keller [START_REF] Keller | Instant radiosity[END_REF] proposed virtual point light method for fast GI computation. Hašan et al. [START_REF] Hašan | Virtual spherical lights for many-light rendering of glossy scenes[END_REF] improved it by introducing spherical lights to avoid singularity. Novák et al. [START_REF] Novák | Virtual ray lights for rendering scenes with participating media[END_REF] used virtual point lights and virtual ray lights for light transport inside translucent materials. Novák et al. [START_REF] Novák | Progressive virtual beam lights[END_REF] replaced virtual ray lights with virtual beam lights to remove singularities. They used importance sampling for the transfer between camera rays and virtual light sources.

Virtual Point Lights and its variants can simulate multiple scattering with high quality, at the cost of large amount of virtual points or rays, resulting in slow convergence speed, especially for highly scattering media. Compared to them, combining our method with VRL, UPBP or MEMLT is much faster, while providing similar results.

Path tracing / Metropolis Light Transport: Georgiev et al. [START_REF] Georgiev | Joint importance sampling of low-order volumetric scattering[END_REF] used the product of geometry and scattering terms for importance sampling in rendering participating media. Compared to our method, the joint sampling method limits to low-order scattering, while our method focuses on high-order scattering. Veach and Guibas introduced Metropolis Light Transport (MLT) for illumination simulation, mutating existing light transport paths to find more efficient transport. Jakob and Marschner [START_REF] Jakob | Manifold exploration: A markov chain monte carlo technique for rendering scenes with difficult specular transport[END_REF] proposed a new mutation strategy for Metropolis Light Transport based on manifold exploration (MEMLT), improving the sampling for paths involving specular and highly glossy surfaces.

Tracing light in highly scattering media usually leads to very long paths. These long paths require a large amount of mutations to provide noise free result. When combining with our precomputed table, the paths are shortened, and much fewer mutations are needed.

Precomputation based Methods: Donner et al. [START_REF] Donner | An empirical bssrdf model[END_REF] precomputed material response on the surface for translucent materials, storing the result as a function of surface position and outgoing direction. They stored double and multiple scattering effects together. Our focus on multiple scattering allows for a more compact representation. We also focus on scattering inside the material, leaving surface diffraction to the original algorithm.

Moon et al. [START_REF] Moon | Rendering Discrete Random Media Using Precomputed Scattering Solutions[END_REF] also precompute material response, but store it on a set of concentric spheres or shells. We use instead a planar grid, resulting in fewer interpolation artefacts. Meng et al. [START_REF] Meng | Multi-scale modeling and rendering of granular materials[END_REF] and Müller et al. [START_REF] Müller | Efficient rendering of heterogeneous polydisperse granular media[END_REF] extend the approach for granular materials, while keeping a concentric shell approach.

Wang et al. [START_REF] Wang | Point-Based Light Transport for Participating Media with Refractive Boundaries[END_REF] and Wang and Holzschuch [START_REF] Wang | Point-based rendering for participating media with refractive boundaries[END_REF] proposed a precomputed table for multiple scattering of a light source in the infinite media. This table was used in a specifically designed extension of Point-Based Global Illumination. We extend their work, but show the approach is compatible with almost all illumination simulation algorithms. We also introduce multiple scattering from a segment source and provide a more compact representation.

Precomputed Multiple Scattering Simulation

Radiative Transfer Equation

We consider a scene containing objects with translucent material. Each of these is assumed to be made of an homogeneous material, with index of refraction η, scattering coefficient σ s , absorption coefficient σ a and phase function p(ω, ω t ). We note the meanfree path inside the material (mfp), with 1/ = σ t = σ s + σ a .

Radiance leaving point x in the direction ω is the sum of exitant radiance from the nearest surface along this direction and in-scattered radiance from the medium among the whole length of the ray [START_REF] Chandrasekhar | Radiative transfer[END_REF] (see Figure 3):

L(x, ω) = T r (x ↔ x s )L(x s , ω) + s 0 T r (x ↔ x t )σ s (x t )L i (x t , ω)dt, (1)
where T r is the transmittance, defined as:

T r (x ↔ x t ) = e -σt x-xt , (2)
s is the distance through the medium to the nearest surface at x s = xsω, and x t = xtω with t ∈ (0, s). L(x s , ω) is the exit radiance from the nearest surface, which is governed by the rendering equation [START_REF] Kajiya | The rendering equation[END_REF]. L i (x s , ω) is the in-scattering radiance at x t from all direction ω t over the sphere of directions Ω 4π using the phase function p, defined as:

L i (x t , ω) = Ω4π p(ω, ω t )L(x t , ω t )dω t . (3)

Notations and Context

We separate between single-, doubleand multiplescattering effects, depending on the number of volume scattering events inside the translucent material. Single scattering corresponds to a light paths with only one scattering event inside the material, double scattering to paths with two scattering events and multiple scattering to paths with more than two scattering events. We only count the number of scattering events, independently of the number of internal reflections on the specular surface. Our algorithm is designed to work with any rendering framework. We describe here integration with Virtual Ray Lights (VRL), UPBP and MEMLT. We focus on computing multiple scattering; the original algorithms compute low-order scattering.

Precomputing Multiple Scattering Effects

Precomputing Multiple Scattering Radiance

In the precomputation step, we assume an infinite, homogeneous, participating medium, with the same properties as the original medium. We run two types of precomputations: with a point directional light source and with a segment light source. A segment source is defined with (o, v, l), where o is the origin and set at the spatial coordinates origin (0, 0), v is the direction of the segment and is aligned with the z axis, l is its length, set at 2× mfp. This is a compromise, considering both the efficiency and the accuracy. A larger size leads to inaccuracy for short light rays, and smaller size requires more accumulation for a given light ray. We treat the segment as a combination of point sources during the precomputation. All spatial dimensions are normalized by the mean-free-path of the material, to reduce the number of parameters. In both cases, we shoot photons from the light source, let them travel in the medium, being scattered or first, we compute lobe shapes indexed by cylindrical coordinates (ρ, z) for position, using spherical coordinates (θ, ϕ) for direction (a). We approximate these lobes using axially symmetric lobes for point source and with asymmetric lobes for segment source, whose axis (red dashed line) is computed to minimize representation error (b). In (c) and (d), we visualize the photon density with varying color and the optimal axis with stream lines.

absorbed, and accumulate their contributions in a table, exploiting the symmetry of revolution of the problem around the direction of propagation. This table is a discretized representation of the multiple scattering radiance in the medium, sampling in both the spatial and directional domain. The core idea behind the table computation is similar to volumetric photon mapping. Each grid cell in the table works as a kernel, gathering all the particles scattered inside it. In each grid cell, we sample the direction and evaluate the radiance at this grid for each sampled direction, with unit watt per steradian per square metre. The radiance functions L o for a point source and a segment source at position (ρ, z) with direction (θ, ϕ) are defined:

r p (ρ, z, θ, ϕ) = 1 N × V ρ,z Mp i=1 Φ i σ s p(θ, ϕ, ω i ), (4)
r s (ρ, z, θ, ϕ) = l N × V ρ,z × K Ms i=1 Φ i σ s p(θ, ϕ, ω i ). ( 5
)
Where N is the total count of photons (we use 500 K in practice, see Table 1). V ρ,z is the volume of grid cell at position ρ, z, M p and M s are the number of particles who reach the cell at position (ρ, z) for a point source and a segment source respectively, Φ i is the flux of the photon i and p(θ, ϕ, ω i ) is the phase function of the medium with incoming direction ω i and outgoing direction (θ, ϕ). K is the count of samples for segment source, and we sample the segment source per grid length uniformly in our implementation. Both problems (point and segment source) have a symmetry of revolution around the direction of propagation. We index points by their cylindrical coordinates (ρ, z) around this axis, and store contributions as lobes in spherical coordinates (see Figure 6). This first representation is 4-dimensional, with a memory cost of 209 MB.

Compressing Multiple Scattering Representation

We then compress this representation, using axially symmetric lobes for point light sources and asymmetric lobes for segment sources (see Figure 4).

The radial plane, defined by ρ and z, is also a symmetry plane for each lobe. We refer to it as the sagittal plane (see Figure 5). The intersection between the lobe and the sagittal plane is a

Lobe intersection

Radial plane

Optimal Axis cross plane Fig. 6. The radial plane, optimal axis and sagittal cross section. The lobe is parametrized as (θ, ϕ), where θ is measured with respect to the optimal axis and ϕ is the azimuthal angle.

sagittal cross-section. Inside the sagittal plane, we find the axis that separates the sagittal section into similar parts by minimizing the difference between the cross-section and its reflection around the axis. This gives us the optimal axis. The plane perpendicular to the sagittal plane that includes the optimal axis is called the coronal plane. See Figures 5, 6 and7.

For point sources, we approximate the lobes using radially symmetric lobes around the optimal axis. For each lobe, we replace the precomputed values for (θ, ϕ) with a sampled function of θ.

For segment sources, we extend the approach, but take into account the asymmetry of the lobes. We begin by approximating the sagittal cross-section using sampled (R, G, B) values for one side. For the other side of the sagittal cross-section and for the coronal section, we store the ratio between computed values and the sampled function for the first side of the sagittal section.

This representation combines accuracy and compactness. Each lobe shape requires 3n samples for point sources, 5n samples for segment sources, where n is the number of sample points in θ. The sagittal section Fig. 7. Sagittal and coronal cross-section. The sagittal plane is a plane of symmetry for the lobe, the coronal plane is perpendicular to it along the optimal axis. For both sections, we visualize both the computed lobe (red) and the approximation (dashed blue). first half of the sagittal cross-section is stored using (R, G, B); the ratios are averaged across (R, G, B) and stored using floats.

Lobe direction is stored using a single float. We use the angle between the optimal axis and the line joining the origin and the current point. Total memory cost after this lobe compression is 9.1 MB.

We further compress our tables using quadtrees, resulting in a memory cost of 0.5 to 1.5 MB depending on material properties.

The strongest approximation we made was for asymmetric lobes for segment sources. We found experimentally that it has a limited impact on overall accuracy (see Figures 7 and8, as well as Section 5.5).

Querying Values From Our Precomputed Tables

Depending on the rendering algorithm, we need to compute multiple scattering from either a point source or a segment source inside the medium, to either a point destination or a segment destination, resulting in four possible combinations (see Figure 9): point-to-point, point-to-segment, segment-to-point and segmentto-segment. VRL and UPBP use segment-to-segment, MEMLT uses point-to-point. Other rendering algorithms could require the other configurations; for example Beam-Radiance-Estimate would require point-to-segment.

Point to Point Contribution. Computing contributions from a point source to a point destination is a direct reading of the precomputed tables. We first transform the point target and ray direction into the table local coordinates, then extract the outgoing radiance from multiple scattering for this point and direction:

z = (P -v) • d, (6)
ρ = (P -v) -zd , (7) mult.
(P, ω t ) = L o ρ , z , T (v,d) (ω t ) . (8)
where (v, d) is the position and orientation of the precomputed table, (ρ, z) are the cylindrical coordinates around the axis of propagation, and T (v,d) (ω t ) is the direction corresponding to ω t in the frame defined by (v, d). Point to Segment Contribution. Contributions from a point source to a segment target begin in the same way (Equations 6 and 7). First, we transform the segment into the coordinates of the point source. We then place point samples on the transformed segment, gather outgoing radiance for each sample using the point sampling algorithm. We multiply this radiance with the transmittance along the segment and sum them:

mult. (S, ω t ) = D c K c Pk L o ρ k , z k , T (v,d) (ω t ) e -σtdk , (9)
where D c is the length of the camera ray, K c is the count of point samples along the camera ray, P k is a point sample along the camera ray which is transformed in the local frame defined by (v, d), and d k is the distance of the sample along the segment target.

Segment to Point Contribution. Contributions from a segment source to a point target also begins in the same way, transforming the point or segment into the coordinates of the segment source. The precomputed segment has a fixed length (2× mfp), but the length of the queried light segment may have different length. We subdivide the querying light segment into several unit segments, each having the same length as the segment source. We assume each of these unit segments have the same contribution to the camera sample, except for attenuation along the light segment. We gather radiance for the transformed point as follows:

mult. (P, ω t ) = v j L o ρ , z , T (v j ,d) (ω t ) e -σtdv j , ( 10
)
where v j is the start point of each unit segment along the segment source and dv j is the distance of this unit segment from the segment source origin.

Segment to Segment Contribution. Contributions from a segment source to a segment target is similar to the contributions from a segment source to a point source, except that we sample the queried camera ray into point samples at first, gather contribution for each of these point samples and then sum them (see Figure 10

): mult. (S, ω t ) = D c K c Pk v j L o ρ k , z k , T (v j ,d) (ω t ) e -σtdv j e -σtdk . ( 11
)
Fig. 11. Application to MEMLT: For each generated seed path, we replace the sub-path with multiple events with one special edge. For example, here the subpath from V 1 to V 3 is replaced by the dashed edge. This special edge carries multiple scattering from one side of the edge to the other. When mutating the path, the first endpoint of the special path is not allowed to be a mutation source. The weight of the special edge is computed from the multiple scattering table.

...

Light Surface Event Volume Event Surface Event Camera

Fig. 12. Definition of subpaths and subspaces. A path x 0 , x 1 , ...x k is subdivided into three sub-paths: x1 , x2 and x3 . x1 starts from the light source and ends at the first volume event. x2 starts from (not include) the first volume event and ends at (not include) the last volume event. x3 starts from the last volume event and ends at the camera. All the sub-paths which satisfy these three conditions construct three sub-path space: Ω 1 , Ω 2 and Ω 3 .

Visibility in the media and Boundaries. We do not consider visibility between the source (segment or point) and the query (segment or point) or boundaries. This is not an issue where the medium has a convex shape, but it can be a source of inaccuracy for scenes with occlusion or air inside the medium. We illustrate this in Section 5.9 and Figures 32 and33.

Application to Existing Rendering Algorithms

Our precomputation tables can be used in combination with many existing global illumination algorithms. In this section we describe integration with VRL, UPBP and MEMLT. For each of them, we let the original algorithm handle single and double scattering, and combine the result with multiple scattering from our precomputed tables.

Combination with Virtual Ray Lights (VRL)

Virtual Ray Lights (VRL) [START_REF] Novák | Virtual ray lights for rendering scenes with participating media[END_REF] is a two-pass ray-tracing algorithm: First render path from the light source inside the medium, with multiple bounces and scattering events and store incoming radiance on these rays. In a second phase, render rays from the camera and accumulate contributions from the rays computed in the first phase.

We follow the structure of the original VRL algorithm, except we do not compute multiple scattering inside the medium. During the first phase, we only store single scattering events in the volume (see Figure 2), but still allow for multiple bounces.

During the rendering phase, we consider rays in the participating medium for each path from the camera. We compute double scattering for these rays using the importance sampling methods described in the original paper [START_REF] Novák | Virtual ray lights for rendering scenes with participating media[END_REF], and multiple scattering using our precomputed tables. We need segment-to-segment contributions. Note that VRL does not compute single scattering effects. We keep the first phase, shooting photons inside the scene, bouncing on surfaces and scattering inside the medium. We restrict scattering to at most two scattering events, and store singlescattered and double-scattered beams separately.

In the final gathering phase, double-scattered beams are only used to compute double scattering, using the original UPBP algorithm. Single-scattered beams are used both to compute single scattering (using the original algorithm) and to compute multiple scattering using our precomputed tables. For the first event on the camera ray inside a medium, we loop over all the beams stored inside this specific medium and gather their contributions to the camera ray. We need segment-to-segment contributions.

Combination with Manifold-Exploration Metropolis Light Transport (MEMLT)

Metropolis Light Transport (MLT) computes paths from the camera to the light source, then mutates these paths to look for more interesting connections. Manifold Exploration is an extension of MLT, where mutations involving path vertices on specular surfaces are constrained by the specularity.

We position our contribution in a modified path integral framework. We first review the original path integral framework. Image pixel intensity I is a integral over the path space Ω:

I = Ω f (x)dµ(x). ( 12
)
x represents a path: x 0 , x 1 , ...x k , which connects the light source and the camera, where x 0 is on the light source and x k is on the camera lens, and x 1 , ..., x k-1 represents vertices on the surface or scattering event in the media. dµx is the differential measurement of the path, corresponding to area integration for surface vertices or scattering events. The measurement contribution function f (x) is the product of the geometry throughput G(x), scattering throughput χ(x), path transmittance T r (x) and visibility V(x). Given path x, without loss of generality, we assume it includes two surface events only and at least three volume scattering event.

x 1 and x k-1 are surface events and x 2 , ..., x k-2 (k > 5) are volume scattering events. We divide this path into three subpaths: x1 , x2 and x3 , where x1 = x 0 , x 1 , x 2 , x3 = x k-2 , x k1 , x k , and x2 = x 2 , ..., x k2 . In subpath x2 , the endpoints are decided by the other two subpaths.

We further divide the path space into three subpath spaces: Ω 1 , Ω 2 and Ω 3 , where Ω 1 represents all the supaths from the light source to the first volume scattering event, Ω 3 represent all the subpaths from the last volume scattering event to the camera and Ω 2 represents the other scattering events in the volume (see Figure 12). Thus, we have

Ω 1 = {x 1 } , Ω 2 = {x 2 } , Ω 3 = {x 3 } . ( 13
)
We rewrite Equation 12with the new subpaths definitions:

I = Ω1 f (x 1 ) Ω3 f (x 3 ) Ω2 f (x 2 )dµ(x 2 )dµ(x 3 )dµ(x 1 ). ( 14
)
We replace the subpath x2 with a special edge (see Figure 11). This edge transports multiple scattering from one path vertex to the other, represented as g(x i , ω i , x o , ω o ). We merge g into Equation 14:

I = Ω1 f (x 1 ) Ω3 f (x 3 )g(x 2 , ω 1 , x k2 , ω k-2 )dµ(x 3 )dµ(x 1 ), (15)
where ω i represent unit vector pointing from x i to x i+1 . The precomputed table is exactly an approximation of function g. Our precomputed table includes G, χ and T r. We ignore the visibility along the subpath. We replace g with a point to point contribution in our table:

I ≈ Ω1 f (x 1 ) Ω3 f (x 3 )mul.(x 2 , ω 1 , x k-2 , ω k-2 )dµ(x 3 )dµ(x 1 ). ( 16
)
To solve the new integral Equation 16, we keep the MEMLT algorithm for all paths outside the medium, as well as for paths with less than two scattering events inside the medium. Both the proposal distribution and the acceptance probability have the same form as MEMLT. We use the multiple scattering evaluated from the table as the weight for x2 and the pdf is set as 1. Both endpoints of the special edge are considered as non-specular, and the end point from the light source can not be considered as the mutation source. The invertibility is not changed in the imagined converged subpaths, except when MEMLT itself has non-invertible cases.

Implementation Details

Precomputed Table Structure

We store our precomputed tables for ρ and z going from 0 to 32 mfp and up to 64 mfp for highly anisotropic materials with large g. Spatial coordinates are normalized by the material meanfree-path (mfp). We compress the spatial representation using a quadtree. The quadtree structure is computed in a preliminary step: first, we shoot a small number of photons to compute the spatial subdivision. We begin with a 32 × 32 regular grid, accumulate photons in each cell, and subdivide cells when the photon count reaches a certain threshold (see Figure 13). Once we have the tree structure, we shoot a large number of photons to compute the lobes inside each cell.

Lobe Approximation

For the lobe approximations, we sample θ every 10 • , resulting in 18 samples for each curve. Including the optimal axis position, the memory cost is 55 floats for each lobe for a point light source, and 91 floats for each lobe for a segment source.

The approximation is very close to what we computed originally (see Figures 7 and8). Figure 14 shows a comparison between the computed and approximated sagittal sections of several lobes in the radial plane, with their position.

Results

Implementation

We have implemented our algorithm inside the Mitsuba Renderer [START_REF] Jakob | Mitsuba renderer[END_REF] for VRL and MEMLT, and inside smallUPBP [START_REF] Křivánek | SmallUPBP[END_REF] for UPBP. All timings in this section are measured on a 2.20GHz Intel i7 (40 cores) with 32 GB of main memory. Unless otherwise specified, all timings correspond to pictures with 512 × 512 pixels. We report the full computation time, for both our algorithm and the environment. Reference solutions are computed using UPBP and MEMLT.

All materials in our scenes are homogenous materials, with Henyey-Greenstein phase functions and refractive boundaries. Material properties are from Křivánek et al. [START_REF] Křivánek | Unifying points, beams, and paths in volumetric light transport simulation[END_REF], Narasimhan et al. [START_REF] Narasimhan | Acquiring scattering properties of participating media by dilution[END_REF] and Holzschuch [START_REF] Holzschuch | Accurate computation of single scattering in participating media with refractive boundaries[END_REF] (see Table 1).

Comparison with High-Quality Methods

Our precomputed multiple scattering greatly speeds up computations of multiple scattering, with all rendering algorithms. Figure 15 shows a comparison between the original methods and our methods, for the three algorithms. We display only multiple scattering for better comparison. Our algorithm produces pictures that are virtually identical to the references, while being at least one order of magnitude faster. We used a point light source for VRL and UPBP, and had to use an area light source for MEMLT.

A strong advantage of our method is that we only change the part of the rendering algorithm dealing with multiple scattering. The rest of the simulation algorithm is unmodified, so we can easily render participating media in a complex setup. Figure 1 shows two such scenes, with multiple refractive interfaces around the participating media. Our algorithm provides noise-free pictures, along with all expected illumination effects, such as diffuse and specular reflections and caustics.

Light transport simulation is more difficult in participating media with refractive interface, due to the absence of a direct path to the light source. Our algorithm also works on index-matched participating media, without refraction at the interface. We have used our method for simulation in mist with varying density (see Figure 16). We observe that the gains are more important on dense mist (see Figure 16 (a)) than for thinner mist (see Figure 16 (b)). This is consistent with our observations on other media: precomputed multiple scattering is especially interesting when multiple scattering dominates.

Comparison with Dipole Methods

We also compare our method to fast computation of multiple scattering using the Dipole [START_REF] Jensen | A practical model for subsurface light transport[END_REF] and Normalized Dipole [START_REF] Christensen | An approximate reflectance profile for efficient subsurface scattering[END_REF] methods, both implemented in the Appleseed Renderer. Figures 17 and18 show a side-by-side comparison between our method (UPBP), UPBP, dipole and normalized dipole. We show all scattering levels (single, double and multiple). For Dipole-based methods, we compute single scattering using UPBP and add it to double and multiple scattering computed using the Dipole approximation.

Dipole methods are very fast. Computing single scattering reduces this timing advantage (600 s for Figure 17), but the method is still the fastest to reach a noise-free picture. It tends to overestimate absorption by thin features, such as the edge of the nose in Figure 17 or the helmet in Figure 18. Our algorithm is slightly slower, but still substantially faster than UPBP, and produces pictures that are indistinguishable from the reference.

Comparison with Diffusion Approximation

The diffusion approximation is based on the observation that the light distribution in highly scattering media tends to become isotropic. In this case, the radiance in the medium from an isotropic point light source may be approximated by a two term expansion involving the radiant fluence φ(x) and the vector irradiance E(x) [START_REF] Jensen | A practical model for subsurface light transport[END_REF]:

L(x, ω) = 1 4π φ(x) + 3 4π ω • E(x), (17)
φ(x) = Φ 4πD e -σtrr(x) r(x) , (18)
E(x) = 1 + r(x)σ tr 4φr(x) 2 e -σtrr(x) r, (19)
D = 1 3σ t , σ t = σ a + σ s (1 -g), ( 20
)
where r is the distance from the light source, σ tr = √ σ a D is the effective transport coefficient and r is the unit vector pointing from the source to x.

The main differences between diffusion approximation and our precomputed table with a point light source are: diffusion approximation is limited to highly scattering media, and approximated all Diffusion Approx. Our Table Source Segment Source Diffusion Approx. Our Extended Table Fig. 19. Comparison between extended precomputed table (including double scattering) and diffusion approximation for a point source in infinite medium and segment source in semi-infinite medium.

media with an isotropic version; diffusion approximation includes both double scattering and multiple scattering, while our method decouples them. Figure 19 shows a comparison between the values computed by diffusion approximation and the values stored in our precomputed table, to which we added double scattering. For an isotropic medium (g = 0), with low absorption (σ s = {2, 2, 2}, σ a = {0.01, 0.01, 0.01}), we get the same results. Figure 20 shows a comparison between the diffusion approximation and our precomputed table (with double scattering added) for illumination simulation on a simple scene, again with an isotropic medium. We get similar results. Figure 21 shows the same comparison, but with a highly anisotropic material (g = 0.9). For this material, the diffusion approximation provides incorrect results, which our table provides results very similar to the reference.

Precomputed Multiple Scattering

The precomputed tables are an important part of our algorithm. Table 1 gives the computation time and memory cost (after compression) for the materials in our test scenes. Memory costs are almost the same for both tables. Materials with a high albedo require longer computations, as particles are less likely to be absorbed at each event.

Moon et al. [START_REF] Moon | Rendering Discrete Random Media Using Precomputed Scattering Solutions[END_REF] store the probability density on a set of concentric spheres or shells, computed in a preprocess step. During the rendering phase, for each scatter event in the medium, they find the sphere with maximum radius with no intersection between the sphere and the geometry. This sphere is sampled to get the next scatter position and direction using the stored probability distribution function. This allows larger step tracing, instead of smaller steps, accelerating computations in participating media (see Figure 23).

Figure 22 shows a comparison between our method and precomputed multiple scattering stored using the concentric spheres [START_REF] Moon | Rendering Discrete Random Media Using Precomputed Scattering Solutions[END_REF], with equal memory cost. Our storage method has the same dimension in space (2), and one less dimension in angle. Quadtree-based compression allows for adaptive representation, resulting in better approximation, especially near the propagation axis. Our method also avoids the Mach banding effects of concen- tric spheres, caused by an insufficient number of samples. We implement Moon et al. [START_REF] Moon | Rendering Discrete Random Media Using Precomputed Scattering Solutions[END_REF] method to compare with our method (VRL) (see Figure 24). For equal-time comparison, our method results in a noise-free picture. Both methods save time in material traversal, using different approaches: Moon et al. [START_REF] Moon | Rendering Discrete Random Media Using Precomputed Scattering Solutions[END_REF] make a bigger step using the sphere radius and the stored probability distribution function, while our method connects the light ray and the camera ray using stored precomputed multiple scattering.

Moon et al. [START_REF] Moon | Rendering Discrete Random Media Using Precomputed Scattering Solutions[END_REF] method was designed for discrete media, such as sand. When used in a continuous medium, we find its efficiency is related to the ratio between the mean-free-path of the medium and the object size. Figure 25 shows a comparison between our method and Moon et al. [START_REF] Moon | Rendering Discrete Random Media Using Precomputed Scattering Solutions[END_REF], for varying mean- Fig. [START_REF] Kajiya | The rendering equation[END_REF]. Comparison between our method (VRL), Moon [START_REF] Moon | Rendering Discrete Random Media Using Precomputed Scattering Solutions[END_REF] and Path Tracing for varying mean free path with equal time (150 seconds). A bigger value of PSNR means better quality.

free-path. Their method works better for smaller mean-free-path (meaning more opaque materials), as more scatter events can be saved. Our method provides better quality over the entire mfp range.

We make a strong approximation by replacing the full lobes in spherical coordinates by asymmetric lobes. Figure 26 shows the result of this approximation on multiple scattering: the impact is not visible in pictures. There is a small decrease on the Peak Signal to Noise Ratio (PSNR). We think this limited impact is related to our focus on multiple scattering. Double scattering lobes have a more pronounced elliptic shape [START_REF] Donner | An empirical bssrdf model[END_REF], even with a point light source. Our Algorithm UPBP Fig. 27. Error comparison between our method (implemented inside UPBP) and unmodified UPBP as a function of rendering time.

Performance comparison

Table 2 reports the PSNR between reference pictures and pictures generated with the rendering algorithms, both unmodified and using our precomputed multiple scatter. This quantitative data confirms the qualitative visual impression: the pictures generated by our algorithm are very close to the reference images, and better than the pictures generated by the unmodified algorithms for small computation times (resulting in larger PSNR). The significant reduction in time spent for each interaction allows us to render more iterations within the required time frame (UPBP), to allocate more virtual ray lights (VRL) or more samples (MEMLT).

Performance analysis: UPBP

Figure 27 shows how PSNR evolves as a function of rendering time for UPBP, both unmodified and using our precomputed multiple scattering, on the Stilllife Scene in Figure 1. Our algorithm contains a strong approximation of multiple scattering effects; as computation time increases, the errors we introduced become dominant, imposing a cap on the PSNR. For long computation time (over 1 h), the unmodified algorithm becomes more accurate.

Our algorithm only applies to multiple scattering. Its advantage depend on the importance of multiple scattering effects in the scene. Figure 28 shows the effect of changing the material meanfree-path while keeping all other parameters constant, along with a density map showing the usage of our precomputed material response. Our algorithm is used often for smaller mean-free-paths, where multiple scattering effects dominate, and scarcely for larger mean-free-paths, where low-order scattering effects dominate.

Performance Analysis: MEMLT

Figure 30 shows how RMSE evolves as a function of the number of mutations for MEMLT, both unmodified and using our precomputed multiple scatter, using a log-log scale. For this test, we used a perfect sphere (see Figure 29) to focus on multiple scattering and remove boundary issues. We use the result rendered with MEMLT with 8196 mutation count as the ground truth. The error decreases as the mutation count increases for both methods. We observe the same behaviour as with UPBP: for a small number of mutations, our method has a smaller RMSE. If the number of mutations is large enough (after 4096), the unmodified algorithm provides better results, as our algorithm is biased. On the other hand, the cost for each mutation is much smaller for our algorithm than for unmodified MEMLT, as we only deal with surface mutations. Figure 31 shows how RMSE evolves as a function of time (in log-log scale). For the same computation time, our method consistently outperforms MEMLT, and provides high quality results in a very short time.

Limitations

Our method was designed for homogeneous, continuous, translucent materials with refractive boundaries. We assume that illumination can propagate freely inside the material. This is not an issue in most of our test scenes, as the translucent material blocks have mostly convex shapes. To illustrate the impact of this limitation, we designed two different scenes: in the first, we place a diffuse black ball inside the Candle scene (see Figure 32). In the second, we create a hollow candle, keeping only a shell of thickness 3 mfp. (see Figure 33). For both scenes, the accuracy of our method decreases as it does not account for all the effects. Our algorithm tends to overestimate scattering effects inside the medium.

Conclusion and Future Work

We present a new method to store precomputed multiple scattering in participating media. Our method is easy to integrate with existing rendering algorithms. We use the original algorithms for single and double scattering, and our precomputed data for multiple scattering. Our method reduces noise, with limited impact on accuracy. Reducing the cost for simulating participating media allows for faster convergence overall. This method is especially interesting for materials with large albedo and small mean-freepath, where higher-order scattering effects dominate. It has limited impact -but also limited cost -for more transparent materials with larger mean-free-path.

In this paper we used precomputed multiple scattering for all scattering events after the second, and used the original algorithms for single and double scattering. An easy extension is to use the unmodified algorithms for more lower-order scattering (up to the Ours (MEMLT) MEMLT Fig. 31. Error comparison (RMSE) as a function of rendering time, using a log-log scale, for both our method (implemented inside MEMLT) and unmodified MEMLT for a simple Sphere Scene. To keep both methods within the same time range the mutation count for our method is from 16 to 4096, and for unmodified MEMLT it is from 4 to 128.

third or fourth scattering event) and use our precomputed tables for subsequent scattering events.

We also want to consider importance sampling to guide longer path length.

Fig. 1 .

 1 Fig. 1. Our algorithm (UPBP) (a) produces realistic pictures of translucent materials. There is substantially less noise than simulation with UPBP in equal time (b); our picture is indistinguishable from reference (c). Our precomputed tables are used only for multiple scattering effects (d), so mostly for highalbedo, small mean-free-path materials. Reducing computation time for these has improved quality for the entire picture. The difference images between our algorithm (UPBP) / UPBP with equal time and the reference are shown in (e) and (f). The resolution of the pictures is 1024 × 512.

Surface

Fig. 3 .

 3 Fig. 3. The Radiative Transfer Equation.

Fig. 4 .

 4 Fig. 4. Multiple scattering from a point source and a segment source:

Fig. 5 .

 5 Fig. 5. Definition of the sagittal plane (blue) and coronal plane (green).

Fig. 8 .

 8 Fig.8. Comparison between computed lobe and resulting approximation, looking at cross-sections perpendicular to the optimal axis. Medium with anisotropy g = 0.9, sample position[START_REF] Jensen | A rapid hierarchical rendering technique for translucent materials[END_REF][START_REF] Jensen | A rapid hierarchical rendering technique for translucent materials[END_REF]. The curves are isolines with the numbers representing the radiance on the right. The solid lines represent the computed lobe and the dashed lines represent the approximated lobes.

  Fig. 9. In rendering, we encounter four different types of radiance estimator.

Fig. 10 .

 10 Fig. 10. Contributions from a segment source (orange line) to a segment target (red line) with length D c . The segment target is transformed into the coordinates of the segment source. We divide the segment source into several unit segments (red square), which have the same length as the precomputed segment source. Each unit segment has a start point v j , with distance dv j from the segment origin. The segment target is sampled regularly into camera samples (red dots), namely P k with a distance d k from the start point of the camera sample. The contributions from the first unit segment are evaluated by directly reading the table, called unit contribution. The contributions from the following unit segments are evaluated by multiplying the unit contribution with the attenuation along the distance.

3. 5 . 2

 52 Combination with Unified Point, Beams and Paths (UPBP) Unified Point, Beams and Paths (UPBP) is an extension of the Photon Mapping algorithm. In a first step, we send photons and beams inside the scene, including the participating medium, storing incoming radiance at each of them. In a final gathering step, we trace rays from the camera and accumulate contributions from neighbouring photons and beams.

Fig. 13 .

 13 Fig.[START_REF] Habel | Photon beam diffusion: A hybrid monte carlo method for subsurface scattering[END_REF]. Precomputed table quadtree structure for segment sources, for materials with anisotropy g = 0.9 and g = 0.3.

Fig. 14 .

 14 Fig. 14. Comparison between computed and approximated lobes, in the sagittal plane, with their position. Segment light source. Material: g = 0.8.

Fig. 15 .

 15 Fig. 15. Comparison (multiple scattering only) between our method and the original methods over two test scenes (bumpy sphere and wax). Our algorithm generates pictures almost identical to the fully converged solution, while being at least one order of magnitude faster. VRL and UPBP use a point light source, MEMLT use an area light source. The BRDFs of both the scenes are dielectric, with IOR 1.5 in Bumpy Sphere Scene and 1.446 in Candle Scene.

Fig. 16 .Fig. 17 .Fig. 18 .

 161718 Fig. 16. Comparison (full solution) between our method (with UPBP) and UPBP on the Bathroom Scene, with added mist. Our method provides a good improvement on dense mist (a,b and c), compared to the thinner mist (d, e and f). Material parameters: (g = 0.4, σ s = {0.016, 0.016, 0.016}, σ a = {0.001, 0.001, 0.001}).

Fig. 20 .Fig. 21 .

 2021 Fig. 20. Comparison between extended precomputed table (including double scattering) and diffusion approximation in UPBP and UPBP as a reference on the Athena Scene. Material: ( g = 0, σ s = {2, 2, 2}, σ a = {0.01, 0.01, 0.01}).

Fig. 22 .Fig. 23 .

 2223 Fig. 22. Comparison between precomputed multiple scattering stored using Concentric Spheres [20] (top) and our method (bottom). Concentric spheres have lower spatial resolution, resulting in poor representation on the axis and Mach banding effects.

Fig. 24 .

 24 Fig.[START_REF] Wang | Point-based rendering for participating media with refractive boundaries[END_REF]. Comparison between[START_REF] Moon | Rendering Discrete Random Media Using Precomputed Scattering Solutions[END_REF] and our method (VRL) with equal time. In this scene, we use a large area light to make path tracing converge faster. The scale of the mfp is 1.

4 Fig. 26 .

 426 Fig. 26. Approximating lobes as 1D functions has little effect on overall accuracy.

Fig. 28 .

 28 Fig. 28. Comparison between our method (in UPBP) and UPBP, with varying mean-free-paths, rendered with 600 seconds. The albedo of all the materials are set as α = (0.24, 0.48, 0.96). g is 0.3. Top row: usage density for our precomputed table.

Fig. 29 .

 29 Fig. 29. Comparison between our method (MEMLT) and MEMLT on a Sphere Scene. For this scene boundaries have little effect, allowing us to focus on multiple scattering convergence. Material: Wax.

Fig. 30 .

 30 Fig. 30. Error comparison (RMSE) as a function of mutation count, using a log-log scale, for both our method (implemented inside MEMLT) and unmodified MEMLT for a simple Sphere Scene. A smaller value of RMSE means better quality.

Fig. 32 .Fig. 33 .

 3233 Fig. 32. Comparison between our method (VRL) and the reference (UPBP) for Candle Scene with and without an opaque diffuse black ball placed inside the candle (multiple scattering only). Our method (VRL) provides brighter results than the reference, as it ignores the visibility issues between camera rays and light rays. Material: wax.

TABLE 1

 1 Parameters for the materials used in this paper. We use 0.5 M particles for the simulations. Memory costs are for both tables, after compression. Computation times are for segment sources (S) and point sources (P).

	Name	α					g mem. time (s)
	R	G	B	R	G	B	MB S P
	B.S. 0.955 0.677 0.457 4.55 3.23 2.17 0.9 1.88 18 2
	Wax 0.980 0.962 0.750 0.65 0.63 0.59 0.8 1.24 28 2
	Milk 0.999 0.999 0.999 0.84 0.75 0.68 0.7 1.32 28 2
	Skin 0.958 0.838 0.678 1.29 0.95 0.67 0.0 0.45 15 2
	Marble 0.999 0.998 0.997 0.46 0.38 0.33 0.0 0.46 27 3

TABLE 2

 2 Computation time, memory costs and error for our test scenes. The count of multiple beams for all the scenes are 10K. VRLs. PSNR VRLs. PSNR time. #iterations. PSNR #iterations. PSNR time. samples. PSNR samples. PSNR

			Ours (VRL)		VRL		Ours (UPBP)	UPBP			Ours (MEMLT)	MEMLT
	scene time. s	K		K		m					m			K
	BumpS. 53	3	34.5	1	23.21	1	269	35.8	80	19.48 1.4	1024	32.12	32	24.01
	Wax	180	40	37.4	1	12.28	3	722	37.3	85	21.01 1.9	1024	35.47	40	23.99
	Stilllife.	-	-	--	-	--	10	1160	31.83	128	21.08 -	--	--	--	--
	Head.	-	-	--	-	--	10	2044	40.59	440	36.06 -	--	--	--	--
	Athena	-	-	--	-	--	20	4093	37.0	600	32.3	-	--	--	--	--
	Our Method		Equal time				Reference

Acknowledgments. We thank the reviewers for the valuable comments. This work has been partially supported by the National Key R&D Program of China under grant No. 2017YFB0203000, the National Natural Science Foundation of China under grant No. 61802187, 61872223, the Natural Science Foundation of Jiangsu under grant No. BK20170857, the fundamental research funds for the central universities No. 30918011320, Open Project Program of the State Key Lab of CAD&CG Zhejiang University under grant No. A1804 and ANR project ANR-15-CE38-0005 "Materials".


Aurélien Thorette
email: aurelien.thorette@univ-rennes1.fr

Marco Romanelli

Steve Bouhier

Frederic Van Dijk

Marc Vallet

Mehdi Alouini

Frédéric Van Dijk

Hybrid opto-electronic oscillator for single-sideband microwave photonics

come

INTRODUCTION

Optically carried single-sideband microwave sources are a key requirement in many microwave photonics system, as they are not sensitive to dispersion fading in long fiber links [START_REF] Smith | Overcoming chromatic-dispersion effects in fiber-wireless systems incorporating external modulators[END_REF]. They can be generated by miscellaneous methods, such as, e.g., optical filtering [START_REF] Park | Elimination of the fibre chromatic dispersion penalty on 1550 nm millimetre-wave optical transmission[END_REF], "period-one" dynamical effects in lasers [START_REF] Zhuang | Tunable photonic microwave generation using optically injected semiconductor laser dynamics with optical feedback stabilization[END_REF] or optical heterodyning [START_REF] Yao | Microwave Photonics[END_REF]. This last technique by which two optical frequencies beat on a photodetector, allows arbitrary RF frequencies and a 100% modulation depth. But it also requires the frequency difference to be somehow stabilized, as each laser will experience different environmental fluctuations. Recently, we have proposed a method based on optical feedback applied to a monolithic DFB laser [START_REF] Vallet | Self-stabilized optoelectronic oscillator using frequency-shifted feedback and a delay line[END_REF]. It relies on the ability of a semiconductor laser to synchronize its phase with an injected field [START_REF] Wieczorek | The dynamical complexity of optically injected semiconductor lasers[END_REF]. It was shown that this property can be used in a configuration similar to optoelectronic oscillators (OEO) [START_REF] Maleki | The optoelectronic oscillator[END_REF], where the beatnote is locked on itself after propagation along an optical delay line. Yet, the long-term stability of this setup was not satisfactory, as it could randomly switch between different resonances of the closed optical delay line, on a typical timescale of a few minutes. Here we propose a more compact setup that allows to overcome the modehopping issue, and paves the way to integration into a single photonic chip.

EXPERIMENTAL SETUP

The two laser beams are provided by a custom-made InGaAsP chip containing two DFB lasers operating at wavelengths λ 1 and λ 2 around 1.5 µm. Such monolithic dual-DFBs, detailed in Ref. [START_REF] Van Dijk | Monolithic dual wavelength DFB lasers for narrow linewidth heterodyne beat-note generation[END_REF], have proven to be efficient for heterodyning, as the lasers are identical by design and experience the same frequency drifts caused by environmental fluctuations. The frequency difference can be controlled by the two pump currents from quasi-DC up to tens of GHz. We choose here to operate at a frequency difference δν = 10 GHz. Additionally, the pump current of one laser can be modulated with a RF signal, superimposed on its DC input by a bias tee (as shown on Fig. 1). The modulation bandwidth of the laser has been measured to be about 14 GHz at -3 dB cutoff frequency.

The lasers are combined by an on-chip -3 dB coupler, and the output is then injected into a L f = 1 m-long microlensed single-mode fiber, in which around 500 µW of power can be coupled. This fiber is connected to a flat FC/PC connector, leadint to a 4% optical feedback into the lasers. The rest of the signal goes through two optical paths of different lengths (L 1 = 5 km and L 2 = 1 km), and each of them ends on a Nortel PP-10G photodiode. The RF signals are recombined by a directional RF coupler. Thus, the summation of the two delayed signals acts as a RF interferometer, with very selective phase resonances. This technique has been proven to be very efficient for reducing spurious phase noise resonances in standard OEOs [START_REF] Yao | Dual-loop opto-electronic oscillator[END_REF]. It has also been used as an optomicrowave frequency discriminator in activly stabilized of dual-frequency lasers [START_REF] Pillet | Dual-Frequency Laser at 1.5 µm for Optical Distribution and Generation of High-Purity Microwave Signals[END_REF].

The resulting signal is then filtered by a 80 MHz BW bandpass filter (Lorch, model CF7), centered at 10 GHz. After an amplification stage, it is fed back into the RF input of one of the lasers, generating small modulation sidebands, i.e., lower than -20 dB around λ 2 . Owing to the optical feedback, these modulation sidebands lead to mutual optical injection between the two lasers, which causes their optical phases to lock, thus stabilizing their frequency difference. A typical output optical spectrum is reported in Fig. 2. The three main peaks correspond to the emission at λ 1 and λ 2 , and to one modulation sideband of DFB2 respectively (the other sideband being coincident with the optical frequency λ 1 of DFB). We point out that the residual peaks at level < -55 dBm are spurs caused by the detection setup (OSA Apex 2083A).

If the output signal is submitted to dispersion fading, for instance in a long fiber link, the presence of a third optical wavelength λ 3 will not deteriorate significantly the quality beatnote modulation, as the amplitude level of this wavelength is 15 db below that of λ 2 . Indeed, the beating between λ 2 and λ 3 may interfere destructively with the beating between λ 1 and λ 2 , but a straightforward calculation shows that the modulation depth reduction is expected to be less than -1.6 dB.

STABILITY

We measure the RF spectrum at the output of a Discovery 401HG photodiode. Without any stabilization, the beanote linewidth is 600 kHz [START_REF] Van Dijk | Monolithic dual wavelength DFB lasers for narrow linewidth heterodyne beat-note generation[END_REF]. When the optoelectronic loop is closed, the linewidth is dramatically reduced to less than 100 Hz. A look at the electrical spectrum on Fig. 3 shows a 40 dB rejection between the carrier and the surrounding spurs. Under our laboratory conditions, the beatnote remains stabilized for hours.

On the largest span (Fig. 3a), it can be noted that external cavity modes, linked to the delay introduced by the optical feedback, and expected to be seen at 1/τ f = c/2nL f = 75 MHz are completely filtered out by the commercially-available RF filter. This removes the possibility for mode hopping on another external cavity mode, and contributes to a greater frequency stability compared to a stabilization scheme with a longer optical feedback [START_REF] Vallet | Self-stabilized optoelectronic oscillator using frequency-shifted feedback and a delay line[END_REF]. The smaller span of Fig. 3b shows the filtering effect of the dual-loop setup, as resonant frequencies from the two RF delays, respectively at 40 kHz and 200 kHz contribute to residual spurs at a level lower than 40dB from carrier. Electrical phase noise of the RF signal was measured using a Rohde & Schwarz FSVP, and is presented on Fig. 4. It reaches less than -70 dBc/Hz at 1 kHz offset from the carrier. Above 20 kHz, sharp spurious resonances associated with the RF delays can be seen, but the two-coils RF interferometer scheme allows to keep the phase noise below -75 dBc/Hz. This is a strong improvement compared to the single-coil scheme, which features peaks up to -50 dBc/Hz. Cautious engineering, based on a precise characterization of the noise sources, could be done to tune the length of the fibers and reduce this value by a few dBc/Hz [START_REF] Lelievre | A Model for Designing Ultralow Noise Single-and Dual-Loop 10-GHz Optoelectronic Oscillators[END_REF]. At low frequencies, better phase noise performances would be reached by implementing techniques found in standard OEO, such as additional feedback loops and temperature control on the fibers.

CONCLUSION

We have presented an efficient way to self-stabilize the frequency difference between two DFB lasers, combining optically injected semiconductor lasers and techniques commonly found in OEO designs, such as the dual RFloop scheme.

At variance with [START_REF] Vallet | Self-stabilized optoelectronic oscillator using frequency-shifted feedback and a delay line[END_REF], we did not use any intensity modulator nor optical amplifier, thus reducing drastically the length of the optical feedback loop and preventing modehopping of the oscillator. We have shown that the use of current modulation does not have an important impact on the amplitude of the generated RF signal, inducing only negligible dispersion fading.

The ability to shorten the optical feedback path opens interesting perspectives towards an integrated selfstabilized heterodyne oscillator. Further miniaturization would require optical components with more integrated features. Preliminary work has already been done in this direction [START_REF] Primiani | Tunable optoelectronic oscillator based on an integrated heterodyne source[END_REF]. Finally, an all-integrated setup would benefit from much more compact optical delay lines, such as optical microresonators.

3 FIG. 2 .

 32 FIG.2. Output optical spectrum, monitored by an optical spectrum analyser (OSA). The output signal is the beatnote at frequency difference δν = 10 GHz.

FIG. 3 .

 3 FIG. 3. Experimental RF spectrum with two different spans. (a) RBW = 1 MHz, (b) RBW = 10 kHz

FIG. 4 .

 4 FIG. 4. Measured phase noise of the RF beatnote

This work benefits from partial funding of European Defense Agency (EDA) through the HIPPOMOS project.