Jean-Baptiste Bellet

MULTIRESOLUTION GREEDY ALGORITHM DEDICATED TO REFLECTIVE TOMOGRAPHY

Keywords: v6.6. 2010 Mathematics Subject Classication. 78A97, 94A12, 65B99, 65Y20. Key words and phrases. Computational optics, reconstruction, acceleration, complexity

Reective tomography recovers the surfaces of a scene to be imaged, from optical images: a tomographic algorithm computes a full volumic reconstruction and then the surfaces are extracted from this reconstruction. For better performance, we would like to avoid computing accurately the full reconstruction, and we want to focus computations on the sought surfaces. For that purpose we propose an iterative multiresolution process. The initialization computes a coarse reconstruction, and the iterations renes it. To identify the voxels to be rened, we take advantage of the asymptotic behaviour of the reconstruction, with respect to its cut-o frequency: it discriminates the surfaces to be extracted. By the way the proposed algorithm is greedy: each iteration maximizes the accumulated intensity of the selected voxels, with prescribed volume. The combination of the complexity analysis and the numerical results shows that this novel approach succeeds in reconstructing surfaces and is relatively ecient compared with the standard method. These works pave the way towards accelerated algorithms in reective tomography. They can be extended to a general class of problems concerning the determination of asymptotically discriminated sets, what is related to the computation of singular support of distributions.

1. Introduction 1.1. Reective Tomography. In optics, reective tomography reconstructs a scene from images of scattered intensity [START_REF] Berginc | Scattering models for 1-D2-D3-D laser imagery[END_REF]. The idea is to apply a tomographic algorithm on reective projections [START_REF] Bellet | Reective ltered backprojection[END_REF][START_REF] Berginc | Simulation of 3D laser systems[END_REF][START_REF] Henriksson | Optical reectance tomography using TCSPC laser radar[END_REF][START_REF] Knight | Tomographic Techniques Applied to Laser Radar Reective Measurements[END_REF]: a ltered backprojection. It oers several ways of representing the scene [START_REF] Berechet | Advanced algorithms for identifying targets from a three-dimensional reconstruction of sparse 3D Ladar data[END_REF][START_REF] Berginc | Optical 3D imaging and visualization of concealed objects[END_REF][START_REF] Berginc | Simulation of 3D laser imaging[END_REF][START_REF] Gering | Object modeling using tomography and photography[END_REF] in three-dimensional optical imaging. The most intense values of the reconstruction are located near the original surfaces, up to artifacts.

The surfaces are extracted using this information, for example by thresholding [START_REF] Berechet | Procédé de discrimination et d'identication par imagerie 3D d'objets d'une scène[END_REF]. 1.2. Goal. This article is part of a search for improvement of the algorithms of reective tomography. Here we would especially like to tackle the computational eorts. The observation is that the technique computes a ltered backprojection, and thus it computes a full volume, whereas we look for surfaces only. To gain performance, we would like to concentrate ne calculations close to surfaces, and limit calculations far from surfaces. The chicken and the egg problem: compute only the voxels close to surfaces that are unknown before computing all the voxels. In this article we design a new algorithm that helps to solve this problem. This is a rst attempt to accelerate the solver of reective tomography using a sequential approach.

1.3. Proposed strategy. To achieve this, we add a priori about the reconstruction. Surfacic extraction by thresholding assumes essentially that the intense values of the reconstructed volume constitute the surfaces.

We have got recently more quantitative results about the reconstruction [START_REF] Bellet | Analyse asymptotique et géométrique de la tomographie réective[END_REF]. Assuming essentially that the input is piecewise smooth, these results formulate the reconstruction as a function of the cut-o frequency of the ltering, when this frequency plays the role of an asymptotic parameter. We propose to take advantage of this dependence to compute a reconstruction which is rened near the surfaces, iteration after iteration. This approach will compute reconstructions for several cut-o frequencies, and thus several resolutions, trying to focus the highest resolutions on surfaces. Such a multiresolution renement has similarities with an algorithm derived in [START_REF] Soussen | 3D reconstruction of localized objects from radiographs and based on multiresolution and sparsity[END_REF]; but the framework and the core of the solver are dierent. 1.4. Main results. The most important result of this article is a new algorithm dedicated to reective tomography: Algorithm 2 page 10. It aims at answering eciently to the following question: compute a small set of voxels which represents the surfaces of the scene, and whose cardinality is set in advance. It is a multiresolution greedy algorithm, guided by the asymptotic behaviour of the reconstruction under the assumptions of [START_REF] Bellet | Analyse asymptotique et géométrique de la tomographie réective[END_REF]. We also derive complexity bounds: Table 1; they show that the method is hoped to be competitive, at least for advantageous occurrences. And we extend these two-dimensional results to a three-dimensional set-up: Algorithm 3 and Table 4. Numerical examples show that the method is relevant on synthesis images with specular reection and diusion, simulated by the Gouraud model [START_REF] Gouraud | Continuous shading of curved surfaces[END_REF], even if the images are corrupted with quite a strong speckle noise. By the way the principle of the multiresolution method concerns a much more general class of problems that we introduce: it deals with the ecient determination of a set from asymptotic discrimination. 1.5. Organization. This paper is organized as follows. We rst recall the principle of reective tomography; we formulate the algorithm that is the reference method in the eld. Then we design a new algorithm in a two dimensional set-up, and we analyze its complexity. We test this algorithm on cases of imaging with occlusions. To nish with, we extend these results in three dimensions, and we show the relevance of the new algorithm, by comparing it with the standard algorithm on simulated optical images, with and without noise.

2. Surface extraction in reflective tomography 2.1. Reective tomography. We recall what 2D reective tomography is.

2.1.1. Record: reectogram. We consider a scene whose objects are opaque; we collect reective projections similarly as [START_REF] Bellet | Reective ltered backprojection[END_REF]. We assume that we measure in this way a reectogram P : (θ, s) ∈ Θ × [-R, R] → P (θ, s), where Θ = [0, |Θ|] ⊂ [0, 2π] is an interval of angles, and R > 0 is the half-width of a screen. More precisely let ∂O = ∂W ∪ ∂S, where ∂W denotes the (boundary of the) wall and ∂S denotes the (boundary of the) objects of the scene. As in Figure 1, for every angle θ ∈ Θ, and for every radial value s ∈ [-R, R], we consider the ray x • θ = s, where θ = (cos θ, sin θ); the scene is projected along this ray, along the direction -θ ⊥ = (sin θ, -cos θ). We record a reective projection on a screen which is parallel to θ, and between the objects and the wall:

(1) P (θ, s) = f (θ, y(θ, s)), s ∈ [-R, R], y(θ, s) ∈ {x ∈ ∂O : x • θ = s}, f : (Θ, ∂O) → R, where y(θ, s) denotes the visible point, and f is the physics of the problem. The opacity assumption means that for the ray x • θ = s, the visible point is the last point of ∂O when we travel along the ray in the direction -θ ⊥ . The reectogram P is the collection of the reective projections: P ≡ (P (θ, •)) θ∈Θ . The function f can take dierent forms depending on the application. In optics, the measurement is the power per unit area incident on the sensor: the irradiance [START_REF] Horn | Robot vision[END_REF]. So P (θ, s) = f (θ, y(θ, s)) represents the amount of light that is reected from y(θ, s) to the sensor. Other example: if f | ∂S = 1 and f | ∂W = 0, we measure 1 if the ray x • θ = s intersects the scene ∂S, and 0 otherwise, and then P represents the silhouettes of ∂S [START_REF] Laurentini | The visual hull concept for silhouette-based image understanding[END_REF].

Silhouettes can be obtained by binarising optical data. For the numerical examples of this article, we will always assume that the wall is black: The acquisition geometry recalls the acquisition of the Radon transform in transmission tomography: the scene is projected along parallel rays. Thus the ltered backprojection, which is a Radon inversion, was introduced to solve this problem [START_REF] Bellet | Reective ltered backprojection[END_REF][START_REF] Knight | Tomographic Techniques Applied to Laser Radar Reective Measurements[END_REF]. The denomination reective tomography" comes from that point:

f | ∂W = 0.
the method applies a tomographic solver on reective data. More precisely for ψ(σ) = |σ| and Ω > 0, we dene a lter whose frequency window is contained in [-Ω, Ω]:

ψ Ω (s) := 1 2Ω Ω -Ω ψ( σ Ω )e iσs dσ = 1 0 σ cos(Ωσs)dσ.
We represent the scene ∂S using the ψ Ω -ltered backprojection of P : Denition 1. For every Ω > 0, the ψ Ω -ltered backprojection of P is the two-dimensional function:

R *

[P ψ Ω ] : x ∈ R 2 → R * [P ψ Ω ](x) = Θ P (θ, •) ψ Ω (x • θ)dθ,
where denotes the convolution with respect to the radial variable:

P (θ, •) ψ Ω (t) = R -R P (θ, s)ψ Ω (t -s)ds.
Of course, the reconstruction R * [P ψ Ω ] is only a heuristic representation of ∂S, because P is not a Radon transform. The half-width Ω of the frequency window plays the role of a resolution parameter. The reconstruction depends on that parameter. We will use this dependency later in the article.

2.2. Discretization. We recall a discrete version of the tomographic algorithm.

2.2.1. Records. We assume that the angular interval Θ is discretized on the grid of the θ j = (j -1)∆θ, 1 j m, with ∆θ = |Θ| m . We assume that the radial interval [-R, R] is discretized on the grid of the t j = -R + (j -1)∆t, 1 j n, with ∆t = 2R n . We assume that the reectogram P (θ, t) is only known on the grid of the (θ j , t k ); the record is thus a uniform discretization of the reectogram P . We approximate the ltered backprojection R * [P ψ Ω ](x), for x scanning the grid of the (t i , t j ), 1 i, j n.

∆ω = π R , α = -n 2 ∆ω, β = n 2 ∆ω; we assume [-Ω, Ω) ⊂ [α, β), i.e. Ω -α = n 2 ∆ω. Let ω k = α + (k -1)∆ω, 1 k n. Then for t = t l , 1 l n, P (θ, •) ψ Ω (t) = π Ω R -R P (θ, s) 1 2π β α ψ1 1 [-1,1) ( σ Ω )e iσ(t-s) dσds.
Notation. For every vector u, whose size is n, the forward, resp. inverse, discrete Fourier transform of u is:

FFT u := e -i 2π n (j-1)(k-1) 1 k,j n u, IFFT u := 1 n e i 2π n (j-1)(k-1) 1 j,k n u.
These notations emphasize that such operations are computed using fast implementations: Fast Fourier Transform and Inverse Fast Fourier Transform.

Notation. Let ˆ denote a product component by component, composed by a shift: for all u, v, of size n,

uˆ v := (u 1 v n 2 +1 , u 2 v n 2 +2 , . . . , u n 2 v n , u n 2 +1 v 1 , u n 2 +2 v 2 , . . . , u n v n 2 )
.

By the left Riemann sum, the ltering is approximated by a discrete Fourier inversion of a discrete spectral product:

(2)

[P (θ, •) ψ Ω (t l )] 1 l n ≈ π Ω IFFT FFT[P (θ, t j )] 1 j n ˆ ψ1 1 [-1,1) ( ω k Ω ) 1 k n .
This rst step approximates the PΩ (θ j , t l ) = P (θ j , •) ψ Ω (t l ), 1 j m, 1 l n, for Ω n 2 ∆ω. The filtering function, written hereafter in pseudo-code, realizes this step for Ω = n 2 π R .

2.2.3. Discretization of the backprojection. Then we discretize the backprojection. By the left Riemann sum:

R * [ PΩ ](x) ≈ ∆θ m j=1 PΩ (θ j , x • θ j ).
The PΩ (θ j , x • θ j ) must be interpolated on the grid of the t l . We assume that t 1 < x • θ j t n , 1 j m; this assumption is satised if |x| R -∆t. Then by linear interpolation:

PΩ (θ j , x • θ j ) ≈ PΩ (θ j , t l ) x • θ j -t l+1 t l -t l+1 + PΩ (θ j , t l+1 ) x • θ j -t l t l+1 -t l , 1 l := x•θj -a ∆t n -1.
Hence this second step computes:

(3)

R * [P ψ Ω ](x) ≈ ∆θ ∆t m j=1 -PΩ (θ j , t l )(x • θ j -t l+1 ) + PΩ (θ j , t l+1 )(x • θ j -t l ).
This is done using the backprojection function.

2.2.4. Reconstruction. Combining the ltering and the backprojection, we compute and we plot the reconstruction:

H(i, j) := 1 1 |x| R-∆t R * [P ψ Ω ](x), x = (t i , t j ), 1 i, j n, Ω = n 2 ∆ω.
This represents the scene under the form of an image H of n 2 pixels.

Function [ PΩ (θ j , t l )] = filtering([P (θ j , t l )], R) Inputs: P (θ j , t l ), 1 j m, 1 l n, with θ j = (j -1)∆θ,

t l = -R + (l -1)∆t, ∆t = 2R n . Ω = n 2 π R , ω k = ( -n 2 + (k -1)) π R , 1 k n, ψ(σ) = |σ|. For every 1 j m, PΩ (θ j , t • ) = π Ω IFFT FFT[P (θ j , t l )] 1 l n ˆ ψ1 1 [-1,1) ( ω k Ω ) 1 k n .
Output: ltered data PΩ (θ j , t l ), 1 j m, 1 l n, for Ω = n 2 π R .

Function H=backprojection(x, [ P (θ j , t l )], ∆θ, R, ∆t) Inputs: x ∈ R 2 and P (θ j , t l ), 1 j m, 1 l n, where θ j = (j -1)∆θ,

t l = -R + (l -1)∆t, ∆t = 2R n . If |x| > R -∆t, H = 0 ; else, H = ∆θ ∆t m j=1 -P (θ j , t l )(x • θ j -t l+1 ) + P (θ j , t l+1 )(x • θ j -t l ), l = x•θj +R ∆t
.

Output: backprojection of P , evaluated at x :

1 1 |x| R-∆t R * [ P ](x).
Algorithm 1 Reference of reconstruction-extraction in reective tomography.

Inputs: reective projections P (θ j , t l ), 1 j m, 1 l n, with θ j = (j -1)∆θ, t l = -R + (l -1)∆t, ∆t = 2R n ; rate α of wished pixels.

Filtering:

[ PΩ (θ j , t l )] = filtering([P (θ j , t l )], R) Backprojection: computation of R * PΩ .
For every 1 i, k n, x = (t i , t k ), H(i, k) =backprojection(x, [ P (θ j , t l )], ∆θ, R, ∆t).

Selection: sort the (i, j, H(i, j)) by decreasing |H(i, j)|, then select the rst αn 2 .

Output: the αn 2 most intense pixels of the reconstruction H (sorted by decreasing intensity).

Surface extraction.

The following problem formulates the question of surface extraction. Problem 2. We know the P (θ j , t k ). Let α ∈ (0, 1]. Represent the scene ∂S using αn 2 pixels.

We use as a priori that the most intense pixels of H represent the surfaces, up to artifacts. Hence among the n 2 pixels of H, we select the αn 2 most intense. A way is to sort the set of pixels, by decreasing intensity |H(i, j)|; then we keep the rst αn 2 pixels. In this article sorting will be realized by a fast algorithm: merge sort. At the end we get the Algorithm 1. We will consider that it is the reference method to solve Problem 2.

2.4. Costs. We give complexity indicators of the reference algorithm. 3. Multiresolution reflective tomography 3.1. Asymptotic result. An asymptotic expansion of the reconstruction R * [P ψ Ω ] is derived in [START_REF] Bellet | Analyse asymptotique et géométrique de la tomographie réective[END_REF], the asymptotic parameter being the frequency parameter: Ω → ∞. Assuming essentially that the record P is piecewise smooth, the reconstruction satises:

Ω 1.5 R * [P ψ Ω ](x) = O (1) , O Ω -0.5 log Ω , O Ω -0.5 .
The intensity of the reconstructed surfacic points is O (1), or O Ω -0.5 log Ω . The intensity of the artifacts is O Ω -0.5 log Ω , whereas O Ω -0. noise. Thus we select a set of pixels whose intensity is high. We divide them into four (sub-)pixels. For all of the new pixels, the normalized reconstruction Ω 1.5 R * [P ψ Ω ] is computed with an adequate value of Ω: from a pixel to a sub-pixel, Ω is doubled. Then we iterate. This process is multiresolution: the reconstruction is computed for several values of the resolution parameter Ω. The idea is to increase iteratively the resolution of the computed intense pixels.

Let us describe the sets of pixels that we rene. We can indeed obtain several variants of the method, depending on this choice. A way is to rene a single pixel at each iteration: the most intense. In this paper we will use the following rule. At the highest resolution, the sought thin pixels ll a surface whose area is A units area (proportional to αn 2 ). For every iteration, we rene the most intense pixels, among the pixels whose resolution is not maximal; we constrain their area: adding these pixels with the pixels at maximal resolution (already computed) must ll an area A. This process is illustrated on Figure 2, with A representing the area of the 12 smallest pixels.

3.2.2. Greedy algorithm. The proposed algorithm can be interpreted as a greedy algorithm for the following problem:

Problem 3. Find a surface S, whose area is A, divided into pixels, and whose accumulated intensity is maximal. Here, the accumulated intensity is dened as the sum of the absolute values of the normalized reconstructions, computed at the center of the pixels.

Each iteration of the algorithm solves a similar problem, from pixels whose resolution is not maximal, and then renes the pixels of the optimal surface S. The iterations are stopped when the pixels whose resolution is maximal ll an area A. Each iteration increases (or does not change) the number of computed pixels at maximal resolution. The area of the surface, that is to be determined such that the thin pixels ll an area A, is a decreasing sequence (of the number of iteration).

In comparison solving Problem 2 by the reference algorithm computes all the pixels at maximal resolution and then solves Problem 3 in a single step. It cannot be proved that the greedy algorithm yields the same result; but it is designed to give a similar result with a lower cost.

3.2.3. Comments. According to the asymptotic result, rening the most intense pixels should rst rene the O (1), and then the O Ω -0.5 log Ω . This is desired because we want to rst rene the surfaces of the scene.

Such renements should yield reconstructions of the same order near the surfaces, even if they are computed for dierent Ω. By the way, for every wrong alarm, i.e. selection of a non-surfacic pixel: the result is O [START_REF] Bellet | Analyse asymptotique et géométrique de la tomographie réective[END_REF] and thus the renement should decrease the intensity. This is requested in order to reduce the noise and also to avoid rening more and more the noise in future iterations. As a result selecting and rening a pixel produces a positive eect for both cases (true and wrong alarms).

Let us notice that it makes sense to compare reconstructions thanks to the normalization factor Ω 1.5 : even on dierent grids, O (1) values can be compared, and they are expected to be larger than O (1) values. On the contrary, without the normalization factor Ω 1.5 , the comments of the last paragraph may become false. In particular comparing reconstructions computed for several values Ω would not make sense. Knowing the asymptotic behaviour of the reconstruction has been crucial to adjust the proposed method. It was used to nd a reasonable normalization; and it is now used to understand the behaviour of the proposed algorithm.

Finally the method does not need to eliminate pixels during the process. The algorithm decides itself which zones must be rened, and it has a desired behaviour for noise (see the wrong alarms above). For every computed pixel, there may exist a future iteration which renes it. This has two advantages. This avoids eliminating prematurely pixels that should be preserved. Secondly, if the aimed area is the area of the full image, then the algorithm converges to the reference reconstruction.

3.2.4. Extension. In this paper, the multiresolution greedy principle is derived and is studied in the context of reective tomography. But the ideas can be immediately extended to a general class of problems, concerning the ecient computation of a set which is discriminated thanks to the high-frequency asymptotics of a band-limited function: Problem 4. Let I Ω : R d → R be a band-limited function, where the resolution parameter Ω > 0 is a bound over the frequencies: the Fourier transform

F[I Ω ] : ξ ∈ R d → R d I Ω (x)e ix•ξ dx is supported in the ball {ξ ∈ R d : |ξ| Ω}.
We assume the following asymptotics, when Ω → ∞:

I Ω (x) = O (1) , if x ∈ S, O (1) , otherwise,
where S ⊂ R d is a small set (small Lebesgue measure). We assume that we are able to compute, or to approximate:

I ω (x), x ∈ R d , 0 < ω Ω,
where Ω is large and xed. Question: compute S.

The brute force method computes I Ω on a full grid associated to the largest resolution Ω and then selects the highest values. The multiresolution method starts the computation at a lower resolution on a coarser grid, and then renes iteratively sets of meshes. If the cost of evaluation of I Ω is high, if d is large, if S is small, and if Ω is large, then the multiresolution greedy process may accelerate the computations.

This class of problems is linked to the numerical determination of the singular support (or support) of a distribution that would be known only through a low-pass lter. For example, up to the normalization,

I Ω (x 1 , . . . , x d ) = 1 i d sinc Ωx i is a low-pass version of the Dirac distribution δ 0 in R d , with S = {0}. Let d = 10
, and Ω = 10 6 , and let us imagine that we have a black box whose input is (x, ω), x ∈ R d , ω Ω, and whose output is I ω (x). It is clear that the multiresolution process would outperform the brute force for determining the set S.

3.3. Computations at several resolutions. We precise the computational formulas, as function of the resolution.

3.3.1. Normalized reconstructions. We assume that n = 2 p , where p 1 is an integer. Let

0 k p. At scale k, the square [-R -∆t 2 , R -∆t 2 ] 2 is divided into square meshes with sides of length ∆t k = 2 p-k ∆t,
and whose center is a point of the grid:

G k = {x k ij := (-R -∆t 2 , -R -∆t 2 ) + ∆t k (i -1 2 , j -1 2
), 1 i, j 2 k }. At the coarser scale, k = 0, we get a single mesh which is the full square. At the ner scale, k = p, we get n 2 meshes, whose centers are the x p ij = (t i , t j ). Let the unit area be the area of a mesh at scale k = p; thus the area of a mesh at scale k is ∆t 2 k = 4 p-k units area. For the reconstruction at scale k,

at x = x k ij ∈ G k , we use the frequency Ω k = 2 k-p n 2 ∆ω.
H k (i, j) := 1 1 |x k ij | R-∆t Ω 1.5 k R * [P ψ Ω k ](x k ij ).
We could use this reconstruction, combining (2) for the ltering and (3) for the backprojection. But in the sequel we will use a more ecient scheme.

3.3.2. Multiresolution ltering. We propose indeed to downsample the ltered data P ψ Ω k (θ j , t l ), 1 j m, 1 l n. We have in mind to speed up the computations for temporary scales k < p, and to keep precision only for the nal scale k = p.

Thus we keep only the l such that t l ∈ -R + ∆t k N and the j such that θ j ∈ 2 p-k ∆θN.

For the scale k 1, let m k = 1 + m-1 2 p-k , θ (k) j = (j -1)2 p-k ∆θ = θ 1+(j-1)2 p-k , and t (k) l = -R + (l -1)∆t k = t 1+(l-1)2 p-k .
Hence the ltered data are:

PΩ k (θ (k) j , t (k) l ), 1 j m k , 1 l 2 k . Notation. Let cut p [x h ] = cuts p [x h ] = [x h ],
and for every k < p:

cut k [x h ] 1 h n = [0, x n 2 -2 k-1 +2 , x n 2 -2 k-1 +3 , . . . , x n 2 +2 k-1 ], cuts k [f i ] 1 i n = [f 1 , . . . , f 2 k-1 , 0, f n-2 k-1 +2 , . . . , f n ].
Hence for every angle θ

(k) j : PΩ k (θ (k) j , t (k) • ) ≈ 2 k-p π Ω k IFFT cuts k FFT[P (θ (k) j , t l )] 1 l n ˆ cut k ψ1 1 [-1,1) ( ω h Ω k ) 1 h n .
In comparison with (2) we cut the frequency range. This is a way to downsample with respect to t, and also to reduce the computational cost. It is indeed sucient to keep the frequencies

ω h ∈ [-Ω k , Ω k ) due to the support of ψ1 1 [-1,1
H ↓ k (ı, ) := 1 1 |x k ı | R-∆t k Ω 1.5 k R * [P ψ Ω k ](x k ı ),
where the notation

• ↓ remembers that R * [P ψ Ω k ] is computed from a downsampling. In a similar manner as (3), for x = x k ı such that |x| R -∆t k : H ↓ k (ı, ) ≈ Ω 1.5 k ∆θ ∆t m k j=1 -PΩ k (θ (k) j , t (k) l )(x • θ (k) j -t (k) l+1 ) + PΩ k (θ (k) j , t (k) l+1 )(x • θ (k) j -t (k) l ), l = x•θ (k) j +R ∆t k .
These computations are now accelerated due to the angle downsampling. The backprojection function is still adequate for such computations.

To summarize, at a scale 1 k < p the reconstruction H ↓ k is considered on a grid G k coarser than G p , with a frequency Ω k lower than Ω p ; downsampling speeds up the computations, and reduces the precision. At the ner scale k = p, there is no downsampling, the full precision is kept; the reconstruction coincides with the reference reconstruction, up to the normalization factor: H ↓ p = Ω 1.5 p H.

PΩ k (θ (k) j , t (k) l ), k 0 k p, 1 j m k , 1 l 2 k . Then we compute the full reconstruction at scale k 0 : H ↓ k0 (i, j), for x k0 ij scanning the grid G k0 ; the (k 0 , i, j, H ↓ k0 (i, j)) are stored in cells: Denition 2. A cell (or a pixel) is a quadruplet (k, i, j, H), where H = H ↓ k (i, j) is the normalized recon- struction, computed at scale k at x k ij ∈ G k .
The intensity of the cell is |H|. The associated mesh is a square with side of length ∆t k and with centre x k ij . The area of the cell is the area of this mesh: 4 p-k units area.

The computed cells are stored in a list L. Then we sort L by decreasing intensity |H|. We dene a list L p , which contains the cells computed with scale p. At the moment L p is void. To nish the task we must insert A = αn 2 cells in L p ; this denes an area A to be lled.

3.4.2. Iteration. Then we iterate the following process. At the beginning of the iteration, L contains cells with scale < p, sorted by decreasing intensity, L p contains cells with scale p. We would like L p to contain αn 2 cells. A > 0 cells are missing; they represent an area to be lled. Hence we rene the most intense cells of L, such that their area ll (at least) the area A. We update L, L p , A. The iteration decreases the area A to be lled as much as it increases L p .

We rst dene an intermediate function: refine. This function renes the most intense cell of L, and returns its area. Since L is assumed to be sorted, the rened cell is the rst cell of L. Writing this cell as (k, i, j, H ↓ k (i, j)), the returned area is a = 4 p-k . We divide this cell into four cells with scale k + 1 p:

(k + 1, I, J, H ↓ k+1 (I, J)), I = 2i, 2i -1, J = 2j, 2j -1. If k + 1 = p,
then the resolution of the computed cells is maximal, and thus they are inserted into L p ; else, we insert these cells into a temporary list L. The original cell (k, i, j, H ↓ k (i, j)) is obsolete after renement; hence it is deleted from L.

Using this intermediate function, an iteration is dened as follows. We rene several cells using an inner loop. The accumulated area of the rened cells is stored in a counter a. The inner loop stops as soon as the accumulated area a reaches the aimed area A. After this loop, the list L does not contain anymore the rened cells, the list L contains the new cells with scale < p, and the list L p has been completed with the new ones with scale p. To update the sorted list L of cells with scale < p, L is sorted by decreasing intensity |H|; then L and L are merged into L. We want |L p | αn 2 ; the area to be lled for this purpose is now A = max(0, αn 2 -|L p |).

Notation. We will sometimes index objects, using the number of the current iteration. At the beginning of the iteration I 1, the thinnest cells (scale k = p) are stored in the list L (I) p , the list of cells with scale k < p is L (I) , and the area to be lled is

A (I) = αn 2 -|L (I) p | > 0; at the end, we get respectively L (I+1) p , L (I+1)
and A (I+1) .

Stopping criterion. The natural stopping criterion is

: A = 0. The area A = max(0, αn 2 -|L p |)
represents indeed the number of missing thin cells; thus we iterate while A > 0.

Remark. We could add a maximum number of iterations I max for safety purposes; we did not write it in Algorithm 2 for simplicity.

We now justify that the algorithm stops, and we count the nal number of cells with maximal resolution. Proposition 1. The Algorithm 2 stops after N iterations, where N is a nite integer.

Proof. At the beginning and at the end of every iteration, adding the area of the cells of L p with the area of the cells of L gives always the total area n 2 . Hence if we enter in the outer loop with a missing area Proposition 2. For the output of Algorithm 2, the number of cells with maximal resolution is:

A = αn 2 -|L p | > 0, then the area n 2 -|L p | of L satises: n 2 -|L p | A.
(5)

|L p | = 4 αn 2 4 .
Proof.

A (N ) = αn 2 -|L (N ) p | > 0 = A (N +1) αn 2 -|L (N +1) p |.
We rst prove that the scale of the last four computed cells is p. If we enter only once inside the inner loop, then the last four computed cells strictly increase |L p | because A (N +1) < A (N ) . Else we enter several times inside the inner loop. In that case, if the scale of the last four computed cells is < p, we consider the penultimate iteration of the inner loop. At the end of this iteration, the nal list L (N +1) p has already been obtained; thus the accumulated area a satises: a |L ) , and thus the inner loop must stop. This is a contradiction: the penultimate iteration is not the last one.

(N +1) p | -|L (N ) p |. Then a A (N ) -(αn 2 -|L (N +1) p |) A (N
We now prove that |L

[ PΩ k (θ (k) j , t (k) l )] k,j,l =filteringMultiResol([P (θ j , t l )], k 0 , R) Inputs: P (θ j , t l ), 1 j m, 1 l n = 2 p , 1 k 0 < p, where θ j = (j -1)∆θ, t l = -R + (l -1)∆t, ∆t = 2R n . ω h = ( -n 2 + (h -1)) π R , 1 h n, ψ(σ) = |σ|.
For every 1 j m, P (θ j , t . ) = FFT[P (θ j , t l )] 1 l n . For every k 0 k p,

m k = 1 + m-1 2 p-k , Ω k = 2 k-1 π R , (t (k) l = t 1+(l-1)2 p-k , 1 l 2 k ) for every 1 j m k , (θ (k) j = θ 1+(j-1)2 p-k ) PΩ k (θ (k) j , t (k) • ) = 2 k-p π Ω k IFFT cuts k P (θ (k) j , t . )ˆ cut k ψ1 1 [-1,1) ( ω h Ω k ) 1 h n . Output: downsampled ltered data PΩ k (θ (k) j , t (k) l ), k 0 k p, 1 j m k , 1 l 2 k . Function a =refine(L p , L, L, [ PΩ k (θ (k) j , t (k)
l )], ∆θ, R, ∆t)

Inputs: L p , L, L: lists of cells with scale p, < p, < p; L is sorted...

From the rst cell (k, i, j, H ↓ k (i, j)) of L to four cells with scale k + 1. Delete the rst cell of L. The area of the rened cell is a = 4 p-k .

K = k + 1, Ω K = 2 K-1 π R , ∆t K = 2 p-K ∆t, ∆θ K = 2 p-K ∆θ. For all (I, J) ∈ {2i, 2i -1} × {2j, 2j -1}, x = (-R -∆t 2 , -R -∆t 2 ) + ∆t K (I -1 2 , J -1 2 ), (= x K IJ ) H = Ω 1.5 K backprojection(x, [ PΩ K (θ (K) j , t (K) l )] j,l , ∆θ K , R, ∆t K ), (= H ↓ K (I, J)) insert the new cell (K, I, J, H) into L p if K = p, into L otherwise.
Output: area a of the rened cell (and updated lists).

Algorithm 2 Multiresolution greedy algorithm dedicated to reective tomography.

Inputs: reective projections P (θ j , t l ), 1 j m, 1 l n = 2 p , θ j = (j -1)∆θ, t l = -R + (l -1)∆t, ∆t = 2R n ; rate α ∈ (0, 1] of wished cells; initial resolution 1 k 0 < p.

Initialization

List of cells: the lists L and L p are void.

Filtering: computation of the downsampled ltered data PΩ k (θ (k)

j , t (k) l ). [ PΩ k (θ (k) j , t (k) l )] k,j,l =filteringMultiResol([P (θ j , t l )], k 0 , R) Backprojection: computation of a coarse reconstruction at scale k 0 , H ↓ k0 . k = k 0 , Ω k = 2 k-1 π R , ∆t k = 2 p-k ∆t, ∆θ k = 2 p-k ∆θ. For all 1 ı,  2 k , x = (-R -∆t 2 , -R -∆t 2 ) + ∆t k (ı -1 2 ,  -1 2 ), (= x k ı ) H = Ω 1.5 k backprojection(x, [ PΩ k (θ (k) j , t (k) l )] j,l , ∆θ k , R, ∆t k ), (= H ↓ k (ı, )) insert the cell (k, ı, , H) into L.
Sort: sort the list L by decreasing intensity |H|. Area to be lled: A = αn 2 .

Iterations

While the number of thin cells is too small, i.e. A > 0. a = 0, L is void.

While the accumulated area of rened cells is too small, i.e. a < A, rene the most intense cell of L:

a = a + refine(L p , L, L, [ PΩ k (θ (k) j , t (k)
l )], ∆θ, R, ∆t).

EndWhile

Sorted list: sort L by decreasing |H|, then merge L and L into L.

Area to be lled:

A = max(0, αn 2 -|L p |).

EndWhile

Output: multiresolution reconstruction L ∪ L p ; L p contains 4 αn 2 4 cells at full resolution.

3.5. Analysis of the algorithm. We study some properties of the multiresolution algorithm in order to establish guaranties and to estimate its performance. See Table 1 page 14 for a synthetic view.

C (1) k0 = n 2 4 p-k0 , C (1) k0+1 = • • • = C (1)
p = 0, and the number of cells with scale < p is: k .

The following schemes describe the evolution of the number of cells of the dierent scales, during the iterations: for every 1 I N ,

C (I+1) k0 -C (I) k0 = -S (I) k0 , C (I+1) k -C (I) k = -S (I) k + 4S (I) k-1 , k 0 < k < p, C (I+1) p -C (I) p = 4S (I) p-1 . (8)
Using ( 5), the nal number of thin cells is: k . During the iteration I the rened cells reach the area A (I) , but if we delete the last rened cell, whose scale is some k I ∈ [k 0 , p -1], the area does not reach A (I) :

A (I) p-1 k=k0 4 p-k S (I) k < A (I) + 4 p-k I , 1 I N.
To nish with, the area to be lled satises: A (1) = αn 2 , A (N +1) = 0, and a recurrence relation:

A (I+1) = A (I) -4S (I) p-1 , 1 I N -1.
m k = O m 2 p-k IFFT of vectors of size 2 k = n 2 p-k : O m 2 p-k n 2 p-k log n 2 p-k . The total cost is: O (mn log n) + O (mn log n) + O m 2 n 2 log n 2 + • • • + O m 2 p-k0 n 2 p-k0 log n 2 p-k0 = O 1 + 1 + 1 4 + • • • + 1 4 p-k0 mn log n = O
O (mn) + O m 2 n 2 + • • • + O m 2 p-k0 n 2 p-k0 = O 1 + 1 4 + • • • + 1 4 p-k0 mn = O 4 3 mn .
3.5.3. Backprojection. The backprojection cost is not easy to evaluate because it depends on the data. We will derive bounds, based on bounds for the number of computed cells: Proposition 4. Let C k be the number of computed cells at scale k during Algorithm 2.

(i) C p = 4 αn 2 4 and the number of cells computed during the iterations satises:

(11) 4 3 (1 -4 k0-p )C p C p + C p-1 + • • • + C k0+1 C p + 1 3 (1 -4 k0+1-p )n 2 .
(ii) The lower bound is reached if, and only if, each cell computed at a scale k 0 < k < p was rened. 4 p-k . Then, due to the renement strategy, for every k 0 < k p, C k is a multiple of 4 and C k-1

1 4 C k . So C k 1 4 p-k C p , k 0 < k p.
Hence the number of pixels computed during the iterations satises:

1 + 1 4 + • • • + 1 4 p-k0-1 C p C p + C p-1 + • • • + C k0+1 C p + 1 4 + • • • + 1 4 p-k0-1 n 2 .
This gives (i). Then (ii) and (iii) are deduced by studying the conditions of equality.

The following proposition estimate the cost of backprojection; see its proof for estimations on the number of computed backprojections. Proposition 5. The total cost of backprojection for the Algorithm 2 is We prove (i) by considering the case where the number of computed reconstructions would be minimal.

(i) at least O π 4 m 2 p-k 0 n 2 4 p-k 0 + O 8 7 (1 -8 k0-p )αmn 2 ; (ii) at most O α + π
The number of cells that we compute during the iterations is at least the lower bound of Proposition 4.

Reaching this bound would mean that the algorithm would converge without useless computations: every computed cell with scale k 0 < k < p would be utilized (rened) to get the nal result. In that case, for every scale k 0 < k p, we would compute

C k = Cp 4 p-k reconstructions of cost O m 2 p-k
operations, with C p = O αn 2 ; the total cost of backprojection during the iterations would be:

O (mC p ) + O m 2 C p 4 + • • • + O m 2 p-k0-1 C p 4 p-k0-1 = O 1 + 1 8 + • • • + 1 8 p-k0-1 αmn 2 .
We now show (ii) by maximizing the number of computed cells; this means reaching the upper bound of Proposition 4. The algorithm would compute C p cells at scale p, and the full reconstruction, at every intermediate scale k 0 < k < p. In that case, the number of backprojections would be:

O α + π 4 1 3 (1 -4 k0+1-p ) n 2 .
The corresponding cost would be:

O αmn 2 + m 2 π 4 n 2 4 + • • • + m 2 p-k 0 -1 π 4 n 2 4 p-k 0 -1 .
The following result establish guaranties concerning the use of memory.

Proposition 6. The total number of cells |L ∪ L p | satises:

|L ∪ L p | 1 4 n 2 + 3 αn 2 4 = O (1 + 3α) n 2 4 .
Proof. Proof. We easily see the cases of equality (ii) and (iii). We prove (i). For every iteration, the increment of the largest scale of the cells of L is at most one. Starting with cells with scale k 0 , we need at least p -k 0 iterations to obtain cells with scale p: N p -k 0 . By the way, we compute at least 4 cells at every iteration; thus 4N is lesser or equal to the number of computed cells during the iterations, which is itself dominated using [START_REF] Henriksson | Optical reectance tomography using TCSPC laser radar[END_REF].

O n 2 4 p-k0 log n 2 4 p-k0 + O n 2 3 (1 -4 k0+1-p ) log n 2 3 (1 -4 k0+1-p ) .
Proof. The initialization sorts a list of size ( n

2 p-k 0 ) 2 , with cost O n 2 4 p-k 0 log n 2 4 p-k 0 . If k 0 = p -1, no additional sort. Else, k 0 < p -1,
S (J) k n 2 3 (1 -4 k0+1-p ) log n 2 3 (1 -4 k0+1-p ).
Proposition 9. The merges of Algorithm 2 cost

O N n 2 4 γ ,
with γ ∈ [0, 1], and the number of iterations N satises [START_REF] Horn | Robot vision[END_REF]. (This cost is null if k 0 = p -1).

Proof. If k 0 = p -1, no additional merge to be considered. Else, k 0 < p -1, each iteration sorts L, and merges L with L. This additional merge costs O | L| + |L| (the cost of the merges of the sort is counted in the sort cost). Summing over 1 I N , the merge cost to be estimated is:

O N I=1 |L (I+1) | . Combining
the area conservation [START_REF] Gouraud | Continuous shading of curved surfaces[END_REF], and a mean value formula:

p-1 k=k0 4 p-k C (I+1) k = n 2 -C (I+1) p = 1 γ (I+1) 1 p-1 k=k0 C (I+1) k , with 1
γ (I+1) 1 ∈ [4, 4 p-k0 ].
We deduce

|L (I+1) | = γ (I+1) 1 (n 2 -C (I+1) p
) from the decomposition [START_REF] Berginc | Optical 3D imaging and visualization of concealed objects[END_REF]. Therefore:

N I=1 |L (I+1) | = N I=1 γ (I+1) 1 (n 2 -C (I+1) p ) = γ 2 N n 2 -γ 3 N I=1 C (I+1) p , with γ 2 , γ 3 ∈ [4 k0-p , 4 -1 ].
We apply 1 J N 1 I J on (8), with C

(1) p = 0. We get:

1 J N C (J+1) p = 4 1 I N (N -I + 1)S (I) p-1 = (1 + (N -1)γ 4 )4 1 I N S (I) p-1 ,
with γ 4 ∈ [0, 1]. We now obtain from ( 9):

N I=1 |L (I+1) | = γ 2 N n 2 -γ 3 (1 + (N -1)γ 4 )4 αn 2 4 ,
which proves the writing:

O N n 2 4 γ , γ ∈ [0, 1].
3.6. Discussion. We discuss the interest of the new method, in comparison with the reference method, concerning the computational eort. The discussion is based on Table 1.

Algorithm 1 Algorithm 2 Filtering

-operations

O (2mn log n) O 7 3 mn log n -storage O (mn) O 4 3 mn Backprojection -number O π 4 n 2 at least O π 4 n 2 4 p-k 0 + O αn 2 4 3 (1 -1 4 p-k 0 ) , at most O αn 2 + π 4 n 2 3 (1 -1 4 p-k 0 ) -operations O π 4 mn 2 at least O π 4 mn 2 8 p-k 0 + O αmn 2 8 7 (1 -1 8 p-k 0 ) , at most O αmn 2 + π 4 mn 2 7 (1 -1 8 p-k 0 ) -storage O n 2 O n 2 4 (1 + 3α) Selection -operations O n 2 log n 2 O n 2 4 p-k 0 log n 2 4 p-k 0 . . . +1 1 k0<p-1 O n 2 3 log n 2 3 + O N n 2 4 , with p -k 0 N αn 2 4 + n 2 4•3 (1 - 1 4 p-k 0 -1 )
Table 1. Costs: reference method (Algorithm 1) versus multiresolution greedy method (Algorithm 2).

3.6.1. Filtering. The additional cost due to downsampled ltering at several scales is at most 33%, for both computation and storage. This step does not determine what method is the most ecient.

3.6.2. Backprojection. We need to compute at least αn 2 backprojections, with cost O αmn 2 , in order to solve Problem 2 using tomography. We compare the two methods with these lower bounds.

The reference method computes 100( π 4α -1)% additional pixels; the corresponding extra cost is 100( π 4α -1)%. Concerning the new method, for the most advantageous cases (lower bound reached), with a small k 0 , the number of computed backprojections would be close to 4 3 αn 2 , i.e. 33% additional backprojections; the extra cost would be close to 100 7 ≈ 14%. For the most disadvantageous cases (upper bound reached) the extra cost would be close to 100 π 4•7α %. Concerning the backprojection cost, if α is small then the new method is always better than the reference one. But if α becomes larger than 6π 4•7 ≈ 0.67, then the new method becomes more expensive for the most disadvantageous cases; if α becomes larger than 7π 32 ≈ 0.69, this stands true even for the most favourable cases.

The multiresolution method is designed to diminish the backprojection cost when we want to compute only a small portion of the full volume at full resolution. This analysis shows that the proposed method achieves its goal. It also reveals that if we want to compute a large portion of the full volume (let us say more than 67%), it is better to directly use the reference method.

3.6.3. Selection. For the main costs of selection (merges and sorts), let us rst consider the case k 0 = p -1.

The cost is dominated by the cost of the sort of the coarse reconstruction: O n 2 2 log n 2 , which is at least four times smaller than the selection cost for the reference. If α is small, then the new method should reduce the total execution time.

It is not so obvious for k 0 < p -1. If the number of iterations is reasonable, i.e. N = O log n 2 , then the cost is reasonable; its order is roughly the order of the sort cost of the reference method. But if N becomes large, in comparison with log n 2 , then the main cost comes from the merges, and is O N n 2

4

. In particular, if N is close to 4m, then the merge cost is close to the global cost mn 2 of the reference method. In such a case the multiresolution algorithm would not really decrease the execution time. Worse still, if N is close to its upper bound, then the merge cost may be close to a n 4 term; this is prohibitive.

Unfortunately the number of iterations is not known in advance. From a pessimistic point of view, we do not know in advance if the multiresolution process will outperform the reference one (for k 0 < p -1). From a more optimistic point of view, the proposed analysis does not take into account the structure of the data set, nor the kind of solver. But in practice the algorithm is not used on arbitrary data, but on reectograms. And the solver (normalized ltered backprojection) has been adjusted for such data. According to the expected behaviour of the renement process, as discussed in the subsubsection 3.2.3, we hope that the costs will be in practice close to the lower bounds (and not the upper bounds).

3.6.4. Performance indicators. The multiresolution algorithm is designed to focus iteratively on the surfaces of the scene, from a coarse initial reconstruction, and to reach a set of thin pixels, whose number is set in advance. The eciency of the algorithm to achieve its goal appears under several forms in the complexity analysis.

Firstly, by the number of computed intermediate cells; this is related with the backprojection cost. For the most advantageous cases, the algorithm would focus on the parts to be extracted and would compute only the necessary cells. For the most disadvantageous ones, the algorithm would spread out the computations, by computing all the cells of all the intermediate scales. We dene an indicator which measures this phenomenon: Denition 3. For the Algorithm 2, the number of intermediate cells

S = C p-1 + • • • + C k0+1 is between the lower bound S 0 = 1 3 (1 -4 k0+1-p )4 αn 2 4
and the upper bound S 1 = 1 3 (1 -4 k0+1-p )n 2 . The focus F is dened by:

F := S 1 -S S 1 -S 0 ∈ [0, 1].
Thus: F close to 1 means that few useless cells were computed, and that the algorithm eciently focused on the parts to be extracted. On the contrary, F close to 0 means that many useless cells were computed. In practice S is computed by incrementing a counter during the iterations, and F is computed at the end of the execution. (Of course, F makes sense only if k 0 < p -1).

Secondly, by the number of iterations N ; this is related with the selection cost, which conditions the interest of the multiresolution algorithm. For the optimal case, N = p -k 0 . If N is a few multiples of 2p, the method is still relevant. But if N becomes very large, in comparison with p, then the new method is not interesting concerning the execution time.

In the numerical tests, we will check the eciency of the Algorithm 2 by checking a posteriori the values of F and N . 3.6.5. Modied algorithm. We can dene variations of the Algorithm 2, without changing the main principle: rening the most intense pixels. For example, downsampling the angles is not mandatory. This would need more computational eort, but would increase the precision of the intermediate computations. This may improve the coarse reconstruction, the selections, and the nal result.

Other example, if we want to be sure that the selection cost is not O n 4 , we can add a lower threshold βn 2 ∈ (0, αn 2 ], set in advance, for the area of the rened zone. The condition of the inner loop, a < A, is replaced by: a < A or a < βn 2 . After N iterations, the accumulation of the rened areas is at least N βn 2 . But this cannot exceed n 2 (p -k 0 + 1), and thus N p-k0+1 β

. The disadvantage of this version: we rene more cells than the strict necessary when A < βn 2 ; this may increase the backprojection cost.

Numerical results

4.1. Implementation. The multiresolution algorithm was implemented in Fortran. The computations are performed using double precision. The Fourier transforms are computed using the library fftw3. For the sorts, we implemented a merge sort, by decreasing absolute value.

Concerning the data structures, we dened a type cell for the cells. The lists of cells are stored in arrays which are dynamically allocated. L p is an array of C p cells. L and L are stored respectively at the beginning and at the end of a single array, whose size is large enough. To locate the beginning and the end of L p , L, L in these large arrays, we dene sentinels: the insertions or deletions, at the top or at the queue of the lists, update the sentinels (increment or decrement).

The code is compiled with gfortran-4.6. The execution takes place on a workstation HP Z820, processors Intel Xeon E5-2609, 2.40GHz. During the execution, we measure time dedicated to: initial ltering, initial backprojection, initial sort, and iterations. 4.2. Full example. We consider the silhouettes of a scene containing two circles: in model ( 1), f = 1 on the circles, and f = 0 on the wall: P (θ, s) = 1 if the ray x • θ = s intersects at least one of the circles, P (θ, s) = 0 otherwise. See Figure 3 for an image of the scene, and the associated reectogram P of size m × n = 805 × 256 (θ is the abscissa, and s is the ordinate). By the way we applied the multiresolution algorithm. On Figure 5, we represent the initial coarse reconstruction of size 32 2 (k 0 = 5). We represent the nal multiresolution reconstruction, and the 0. At the end, the extracted reconstruction is not as nice as the reference reconstruction, but it looks similar, and the computation time was divided by 10. We observe some convergence of the reconstructions to the reference reconstruction. We notice that for k 0 = p -1, p -2 (where p = log 2 n = 9), the results after extraction are relatively close to the reference.

For the smallest values of k 0 , i.e. for the coarsest initial grids, the result contains relevant information, but with a degradation (holes). We grant that k 0 cannot become too small: the algorithm is based on the high-resolution asymptotics of the reconstruction. Table 2. Performances as a function of the initial grid, for the reconstructions of Figure 7.

The performance indicators are in

is reasonable, in comparison with 2p = 18. The measured times are consistent with these indicators: the multiresolution method is clearly faster than the reference. Of course, for k 0 = p -1, the method converges after a single iteration.

To conclude this example, the case k 0 = 7 = p -2 is a good compromise between computation time and quality of the extracted reconstruction. We propose more generally to set k 0 such that p -k 0 is a small integer. in priority the circles, then the straight artifacts, then the noise. This is explained by the asymptotics:

Ω 1.5 R * [P ψ Ω ](x) is O (1)
on the circles, O Ω -0.5 log Ω on the lines, and O Ω -0.5 elsewhere. This sets the priority order when we select the most intense pixels. That was the idea of the algorithm. We also observe that α must be appropriate: if α is too small, the reconstructed part contains too few points; increasing α increases the number of reconstructed points, but it also increases the artifacts and the noise. In practice we could imagine several steps: start with a small α, then launch the iterations again to rene more cells, and so on.

In Table 3, we reported the performance indicators. The execution time increases lightly with α when α is very small; then it increases sublinearly (from β = -4). In particular we cannot hope vanishing this time, even if α is very small: the method contains incompressible costs (including the initialization).

The other important conclusions of the table concern the focus. We observe that the focus F decreases when α increases. The main decreases appear at the end of the table. We must link this with the Figure 8. The algorithm focuses on some parts of the scene, by taking advantage of the asymptotic behaviour when Ω β = log increases. As soon as the parts corresponding to the right behaviour are reconstructed at full resolution, the algorithm works on pixels which contains noise, and not on signals containing the sought behaviour. Thus the algorithm cannot focus anymore on well targeted regions, and so F decreases.

We have already mentioned that the focus F gives information on the eciency of the algorithm. We now see that it gives indications about the quantity of useful signal in the nal image. If F is very close to 1, increasing α (or relaunching iterations) may improve the quality of the reconstruction. On the contrary, if F is small, the computed reconstruction may be already corrupted by noise, and it may be preferable to keep only the most intense computed pixels; or we can alternately start again with a smaller α. These comments also suggest another variation of the algorithm, taking into account the focus in the stopping criterion: continue while the focus is above some xed threshold. 4.5. Complex scene. We add an example which is more dicult: see Figure 9. The size of the reectogram is m × n = 3217 × 1024. The scene contains: a small circle with f = 1, a large circle with f = 0.77, a polygon with f = 0.63, a croissant with f = 85 on the convex portion, and two parts for the concavity, f = 0.63 and f = 0.73. The reconstruction Ω 1.5 R * [P ψ Ω ](x) is essentially expected to be O (1) on convex portions, and O Ω -0.5 log Ω on lines and isolated points.

On Figure 10, we represent the full reference reconstruction, of size n 2 = 1024 2 , and the extraction of the 0.05 × 1024 2 most intense pixels. The reconstruction and the extraction took 286 seconds. On the second line of Figure 10, we represent the reconstructions of the multiresolution method: for the initial grid of size 128 2 (k 0 = 7), and the nal 0.05 × 1024 2 thinnest pixels. For this test the focus is F = 0.948, and the number of iterations is N = 33; the test took 25.9 seconds. By zooming we can see that the reference is slightly better; but it took 11 times more time.

Remark. For clarity, the images of Figure 10 do not represent H, but |H| sign(H). For the visualization, we project the thinnest voxels on three orthogonal planes, using the Maximum Intensity Projection (MIP -see for example [START_REF] Wallis | Three-Dimensional Display in Nuclear Medicine and Radiology[END_REF]): for each plane of projection, along every line orthogonal to the plane, we project the intensity of the most intense positive voxel (or 0 if all the voxels of the line are negative). For clarity we apply the square root at each image, before printing. 5.4.3. Reconstructions at 5%. On Figure 12, we represent the reference reconstruction, before extraction, and after extraction of α = 5% of the voxels. The computation took 8680 seconds. On Figure 14, we represent the greedy reconstruction, with initialization of size 128 2 • 512, and α = 5%. The computation needed 24 iterations, for a total time of 1060 seconds; this is 8 times smaller than the reference time. The focus is 0.962. To observe the dierences between the extracted reconstructions, we subtracted the MIPs of the greedy reconstruction from the MIPs of the reference: see Figure 13. The renderings of the two methods are nally similar. 5.4.5. Speckle noise. We would like to test the algorithm on a more challenging case with speckle noise. We simulate Gouraud images of the Stanford Bunny as before, but with a striped pattern on the surface: see Figure 16. After the rescaling P := 1 + P max(P ) , the range of the values of P is [START_REF] Bellet | Analyse asymptotique et géométrique de la tomographie réective[END_REF][START_REF] Bellet | Reective ltered backprojection[END_REF]. Then we add a speckle noise: each value P (θ, s, z) is replaced by P (θ, s, z) := P (θ, s, z)(1 + ε(θ, s, z)), where ε contains independent realizations of the gaussian N (0, 1). After these changes, the rst line of Figure 16 becomes the second one.

Reconstructing the surface with such a level of noise is quite challenging.

We compare the reference with the greedy method for α = 0.001, and initial size 64 2 • 512: see Figure 17. The computation took 8760 seconds for the reference, while it took 237 seconds for the greedy method (the ratio of time is 37). The greedy method converges in 96 iterations, with focus 0.996. The initial reconstruction of the greedy method is regularized thanks to the frequency cut; this helps the method to concentrate the computations on the surfaces. At the end, the cloud of the greedy reconstruction looks less diuse than the reference one. 6.1.1. Multiresolution greedy algorithm. This solution is a multiresolution greedy algorithm, derived from the high-frequency asymptotics of the reective ltered backprojection. The idea: start from a coarse reconstruction on a coarse grid, then rene iteratively sets of voxels. Each iteration selects and renes a set of voxels maximizing the accumulated normalized ltered backprojection, with a volume constraint.

In comparison with the reference, this approach does not compute the full volume at full resolution; the computations focus themselves on the regions of interest. The principles of the algorithm can be extended for a general class of problems, concerning the ecient determination of a set xed by asymptotics. 6.1.2. Performances. The expected consequence is a reduction of the computation time. Concerning the complexity, we derived bounds for the costs. The lower bounds show that the method can be competitive for the most advantageous occurrences. The upper bounds are prohibitive, but they do not take into account the special structures of the data and of the solver. We measure a posteriori the eciency of the algorithm, by looking at the number of iterations and a focus indicator. The observed values in the numerical tests reveal that the algorithm quickly converges and eciently focuses. 6.1.3. Settings. There is a price to be paid; two parameters must chosen in advance: initial resolution and rate of wished voxels. If they are too small, the reconstruction might be too sparse: pieces of surfaces not reconstructed. Based on empirical observations, a good compromise between computational eort and quality of the result: reconstruct a few percents of the full volume, from an initial resolution that is near the full resolution. In such a case, the reconstruction is expected to be close to the reference, but with a signicant reduction of computation time.

6.1.4. Accelerated solver. At the end, the proposed method is an alternative to the reference method, and especially for large volumes of data. It accelerates the reference solver and it oers new opportunities to develop an algorithm which is as fast as possible. 6.2. Perspectives. We suggest several options to be explored. 6.2.1. Complexity. The theoretical analysis of complexity must be deepen. The bounds that we derived do not explore the special structure of the data, nor the form of the solver. It would be interesting to rene the bounds, in a framework to be determined. This question is dicult. Another point to be studied is the average-case complexity. 6.2.2. Filtering. We used the Ram-Lak lter: the method stands on the asymptotics of the corresponding reective ltered backprojection. We can tune the process for other lters. If we want another lter, we must rst study the associated asymptotics to nd the right normalization, and to check that the algorithm will have the expected behaviour. 6.2.3. Parallel computing. By the way, in tomography, a standard way of gaining in performance is parallel computing. Here we introduced an acceleration principle in a sequential algorithm. It could be interesting to derive an algorithm which combines such a principle with parallel computing. 6.2.4. Extension. To nish with, the extended problem about the determination of an asymptotically discriminated set is relatively general, and thus it should be possible to apply the multiresolution greedy principles to other elds. It would be worth checking the relevancy of the concept on other problems that involve asymptotic methods such as stationary phase analysis.

Figure 1 .

 1 Figure 1. Reective projection for an angle θ ∈ Θ: the wall ∂W and the scene ∂S are projected in the direction -θ ⊥ on a screen, along the rays x • θ = s, s ∈ [-R, R]. The set of the projected visible points, y(θ, s), s ∈ [-R, R], is represented in bold. The dashed lines represent two rays of projection, from the visible points y(θ, s) and y(θ, s ), to the screen.

Figure 2 .

 2 Figure 2. Principle of the multiresolution algorithm: compute an image with coarse pixels, then rene sets of intense pixels, iteratively.

  .

3. 2 . 1 .

 21 Iterative rening. First we compute a normalized reconstruction Ω 1.5 R * [P ψ Ω ](x), for x scanning a coarse grid, and Ω being associated with the mesh size. Then we iteratively rene this reconstruction as follows. The intense values are expected to represent surfaces, and the low values are expected to represent

| - 4 < αn 2 .

 42 If the inner loop contains only one iteration, it is clear: |L (N +1) p | -4 = |L (N ) p | < αn 2 . Else, at the end of the penultimate iteration of the inner loop, if |L(N +1) p | -4 αn 2 , then a |L (N +1) p | -4 -|L (N ) p | A (N ) -(αn 2 -|L (N +1) p | + 4) A (N ), which is a contradiction as before.3.4.4. Solution to the main problem. The Algorithm 2 computes a set of cells: L ∪ L p . Displaying these cells, using their meshes and values, provides a multiresolution image. We can also display only the thinnest cells: L p . The obtained image is the new solution that we propose to solve the Problem 2.

Function

3. 5 . 1 .

 51 Evolution and conservation laws. Let 1 I N + 1. For every k 0 k p -1, let C (I) k denote the number of cells with scale k in L (I) , and let C (I) p the number of cells (with scale p) in L (I) p . In particular, after the initialization:

p- 1 ,

 1 of cells with scale k which are selected and rened. The inner loop builds L(I) , with L(I) = 4 |L (I+1) | = |L (I) | -S

4 .

 4 Furthermore the area of the cells of L (I) ∪ L

Proposition 3 . 7 3

 37 The ltering of Algorithm 2 costs O mn log n operations. The ltered data are stored in an array of size O4 3 mn . Proof. The ltering starts by m FFT of vectors of size n: P (θ j , •) = FFT[P (θ j , •)], 1 j m. This costs O (mn log n) operations. Then for each frequency Ω k , k 0 k p, the main cost comes from

7 3

 7 mn log n ; the constants of the O are of the same order than the constant of the cost of FFT.After the computations, we save: an array of sizem k × 2 k = O m 2 p-k n 2 p-k , for each k 0 k p.The storage cost is:

2 , thus C k n 2

 22 (iii) The upper bound is reached if, and only if, all the cells of all the intermediate scales k 0 < k < p were computed. Proof. At the end of the algorithm, for the scale k = p: C p = C (N +1) But more generally, for all k p and I N + 1, we have only C k C (I) k , due to deletions. The size of the full image at scale k is n 2 p-k

4 1 7 ( 1 -

 471 8 k0-p ) mn 2 . Proof. The initialization computes a reconstruction on the grid G k0 of size n 2 4 p-k 0 , from m k0 angles. This costs O π reconstruction. We must now add the costs of the iterations.

(ii) If each iteration 1 I p -k 0 renes αn 2 4 4 + n 2 4 • 3 ( 1 -

 124431 p-k 0 -I+1 cells with scale k 0 + I -1, then N = p -k 0 . (iii) If each iteration renes a single cell, and if the number of computed cells is maximal (upper bound of (11)), then N = αn 2 4 k0+1-p ).

3. 5 . 5 .

 55 Selection. We count the main costs related to the selection of intense values: merge sorts and (extra-)merges. Proposition 8. The merge sorts of Algorithm 2 cost:

  each iteration sorts the list L of rened cells with scale < p, with cost O | L| log | L| . Summing these costs over the iteration number 1 I N , the cost to be estimated is After the initialization |L (1) | = n 2 4 p-k 0 , and C (1) p = 0. And at the end, from (4), |L (N +1) | n 2 -C (N +1)

Figure 3 .

 3 Figure 3. Scene with white circle (left) and its reectogram 805 × 256 (right).

Figure 4 .

 4 Figure 4. Reference for the reectogram of Figure 3: full reconstruction of size 256 2 (left), and after extraction of 5% of the pixels (right).

Figure 5 .

 5 Figure 5. Multiresolution reconstruction for Figure 3: initialization of size 32 2 (left), multiresolution image (middle), and the 0.05 × 256 2 thinnest pixels (right). On Figure 4, we represent the full reconstruction of size n 2 = 256 2 , and the reconstruction after extraction of the most intense pixels, α = 5%. Computing the reconstruction and the extraction took 3.61 seconds.Qualitatively we nd what was already known: the most intense pixels correspond to the surfaces of the scene, up to artifacts (the artifacts being the lines that are tangent to both circles).

  05×256 2 thinnest pixels. This multiresolution computation took 0.364 second. Qualitatively we observe what was announced: the pixels are rened near the surfaces (and artifacts) and they stay coarse elsewhere. The algorithm succeeds in focusing the computations near the relevant places. The indicators show the eciency: we got F = 0.951 for the focus, and N = 25 = 1.5 • 2p + 1 iterations.

4. 3 .Figure 6 .

 36 Figure 6. Scene (left) and its reectogram 1609 × 512 (right).

Figure 7 .Figure 6 .Figure 8 .

 768 Figure 7. Reconstructions for the Figure 6: coarse initialization of size (2 k0 ) 2 (left), multiresolution reconstruction (middle), and the thinnest pixels (right). From top to bottom: k 0 = 4, 5, 6, 7, 8; on the last line: reference reconstruction.

Figure 9 .

 9 Figure 9. Scene (left) and its reectogram 3217 × 1024 (right).

Figure 10 .

 10 Figure 10. Reconstructions for the Figure 9. Top: reference before (left) and after extraction (right). Bottom: initialization (left) and nal reconstruction (right) by the multiresolution method.

5. 4 .

 4 Numerical tests. 5.4.1. Simulation of optical images. Using the triangulation of the Stanford Bunny [16], we simulated a sequence of m = 805 images of size n × n = 512 × 512, by turning around the bunny; the angular step is 2π m . These images are computed using the Gouraud model of Matlab, for orthographic projections (parallel rays). These synthesis images contain a component of scattering, and a component of specular reection.See Figure11for a few images of the sequence. The pixels values depend particularly on the light source position, the pixel location and the normal to the surface at the visible point.

Figure 11 .

 11 Figure 11. Simulated images of the Stanford Bunny, using the Gouraud model.

Figure 12 .

 12 Figure 12. Reference reconstruction: full volume of 512 3 (top), and extraction of 0.05 • 512 3 voxels (bottom).

Figure 13 .

 13 Figure 13. Greedy MIP subtracted from the reference MIP, for α = 5%.

Figure 14 .

 14 Figure 14. Greedy reconstruction: initialization of size 128 2 •512 (top), and reconstruction at 5% (bottom).

Figure 15 .

 15 Figure 15. Reference reconstruction (top), greedy reconstruction with initialization 128 2 • 512 (middle), and dierence between the two (bottom); α = 1%.

Figure 16 .

 16 Figure 16. Gouraud images of a striped bunny: without noise (top) and with speckle noise (bottom).

Figure 17 .

 17 Figure 17. Reconstructions for the noisy striped bunny. Reference reconstruction (top), greedy reconstruction with initialization 64 2 • 512 (bottom); α = 0.001.

1 .

 1 Synthesis. These works provide an original solution to the following problem: extract the surfaces of a scene by reective tomography.

  2.2.2. Discretization of the ltering. First the ltering is computed for every angle θ of the grid. Let

  2.4.1. Filtering. The size of the ltered array PΩ is the size of the data array P : mn. For every θ = θ j , 1 .4.3. Selection. Lastly the cost of surface extraction is predominated by the cost of a merge sort of n 2 elements: O n 2 log n 2 . 2.4.4. Synthesis. The costs are summarized in the rst column of Table1. We will choose m and n of the same order. Hence the total cost is dominated by the cost of backprojection: O m π 4 n 2 . That is the reason why we would like to tackle this cost in the sequel: by decreasing the number of backprojections, and by reducing the cost of some backprojections.

	j m, we compute FFT[P (θ j , •)] with cost O (n log n), we multiply in Fourier space with cost O (n), then we
	compute	π Ω IFFT[•] with cost O (n log n). In total, computing the array PΩ costs O (2mn log n) operations.
	2.4.2. Backprojection. The array H contains n 2 elements. Backprojections must be computed for the (t i , t j )
	inside the disk |x|	R -∆t. Dividing the area of the disk |x| < R by the area of the square [-R, R] 2 ,
	we get operations. In total, O m π π 4 ; thus the number of computed backprojections is about 4 n 2 operations are needed to compute H. π 4 n 2 . Every backprojection costs O (m)

2

  5 represents a noise. This result is consistent with the numerical observations: the highest values represent the surfaces of the original scene. Let us compute the normalized reconstruction Ω 1.5 R * [P ψ Ω ] for several values of Ω, large enough. If the asymptotics at x is O (1), the result is expected to decrease when Ω increases. But if the asymptotics at x is O (1), the result should not really change. We propose to take advantage of this model for better performance: injecting this a priori into the solver should increase the eciency of the reconstruction-extraction process.

	But this result is stronger: the heuristic reconstruction is quantied as a function of the resolution
	parameter Ω.

3.2. Principle. The main idea that we will develop is described below; see Figure

2

  ) . The cut is operated using cut k and cuts k . In order to obtain a hermitian spectral signal and a real signal, we cancelled the Fourier coecients associated with the frequency -Ω k , k < p (see the denition of cut k and cuts k ). Due to these cuts, the IFFT is applied on a vector whose size is 2 k ; at the end, the points of the spatial grid are the t

	(k)

l , 1 l 2 k . The filteringMultiResol function computes the downsampled ltering for the dierent resolutions. 3.3.3. Ecient normalized reconstructions. Then for all 1 ı,  2 k , let:

  Thus the inner loop will stop: it is possible to nd and rene a set of cells of L, whose area is at least A. By the way the area of L, n 2 -|L p |,

	is lled by at most	1

[START_REF] Berechet | Procédé de discrimination et d'identication par imagerie 3D d'objets d'une scène[END_REF]

(n 2 -|L p |) cells with scale < p. Thus:

(4) |L| 1 4 (n 2 -|L p |).

Concerning the outer loop, if the iterations do not stop, we consider the successive missing areas: A = αn 2 -|L p | > 0. This sequence decreases and takes its values in a nite set. Thus it is constant after a nite number of iterations; so is the number of thin cells |L p |. But the number of cells of L ∪ L p is |L| + |L p | and strictly increases. As a result |L| indenitely increases. This is in contradiction with the bound (4). This proves that the algorithm stops.

  The stopping criterion implies |L p | αn 2 . Since L p is built by adding 4 cells by 4 cells, |L p | is a multiple of 4. It remains to prove that |L p | < αn 2 + 4.Let us consider the last iteration, whose number is N :

  .5.4. Number of iterations. The following Proposition gives bounds for the number of iterations; if we do not use the structure of the data and of the solver, they are optimal:

	Proposition 7. (i) Algorithm 2 stops after N iterations, with
	(12)	p -k 0 N	αn 2 4	+	n 2 4 • 3	(1 -4 k0+1-p ).

By combining the inequality (4) and |L p | C p = 4 αn 2 4 , the number of cells in memory always satises: |L| + |L p | 1 4 n 2 + 3 4 C p (also true during the iterations).

3

Table 2 .

 2 The focus is always very high. The number of iterations N size 2 k0 focus F iterations N time (s)

	16	0.984	45	0.728
	32	0.981	43	0.752
	64	0.984	24	0.784
	128	0.988	13	1.12
	256	×	1	4.14
	512	×	0	31.8

Table 3 .

 3 Performances as a function of the rate of pixels α, for the reconstructions of Figure8.

J.-B. BELLET

Three-dimensional extension

We extend the multiresolution greedy algorithm to a problem in three dimensions, for an acquisition by slice.

5.1. Reconstruction by slices. We measure a three-dimensional reectogram, by slice: the collected reective projections are the P (θ, s, z), where for every altitude z, P (•, •, z) is a reectogram in the horizontal plane x 3 = z. The full reconstruction is obtained horizontal slice by horizontal slice: for every altitude z, the reconstruction, restricted to the plane x 3 = z, is (x 1 , x 2 ) → R * [P (•, •, z) ψ Ω ](x 1 , x 2 ). We assume that the parameters θ, s, x 1 , x 2 are discretized as before. Concerning the altitude, we assume that z scans the values z l , 1 l n, with constant step ∆z = z l+1 -z l .

Similarly as the two dimensional problem, the question is: represent the 3D scene using at least αn 3 thin voxels, where α ∈ (0, 1]. The reference method consists in computing a full reconstruction of size n 3 , by computing the reconstruction slice by slice, i.e. for every z l , 1 l n. Then the voxels are sorted by decreasing intensity, and we select the rst αn 3 voxels.

We would like to take advantage of the principle of Algorithm 2 in order to outperform this 3D process.

The easiest way is the following: we directly apply Algorithm 2, for every slice of altitude z l , 1 l n. Once Algorithm 2 is implemented, this is easy: we have only to add a loop over l. But this approach has a disadvantage: the 3D structure of the problem is not utilized. We will prefer operating selection on the whole volume, and not slice by slice. This is what we are going to do. 5.2. Multiresolution computations. We assume: n = 2 p with p 1 an integer. Let 0 k p. At scale k, we divide the volume

, and whose center is a point of the grid:

Thus a renement (increment of scale) is operated in the two horizontal directions; vertically, the discretization is always the thinnest possible. At the coarser scale, k = 0, the volume contains n slices with width ∆z.

At the thinner scale, k = p, the volume contains n 3 thin voxels, with (x p ij , z ) = (t i , t j , z ). The volume unit is chosen such that the volume of a parallelepiped with scale k is ∆t 2 k ∆z = 4 p-k volume units. The altitude z of every point of G k belongs to the altitudes where P is known. Thus, at scale k, at

n, is downsampled, as in 2D.

5.3. Multiresolution greedy algorithm.

5.3.1. Principle. We use the same principle as before: we start from a coarse reconstruction, and we iteratively rene the most intense voxels. At each iteration, we rene a set of voxels whose volume lls the volume to be lled. Essentially, in comparison with 2D, we must add an altitude coordinate to the objects.

5.3.2.

Filtering. The ltering is operated angle by angle, scale by scale, altitude by altitude, and with a downsampling in t and θ. Here we take care of memory usage, so we combine reading of the data and ltering. We avoid the simultaneous storage of several full volumes: see the function filteringMultiResol3D. For every angle θ j , the image with angle θ j is loaded and ltered in Fourier space; then, for every scale k such that the angle θ j is selected by the downsampling in θ, for every altitude z , we invert the cut Fourier data.

5.3.3.

Iterative rening of cells. The cells that we now dene contain also an altitude. And now, rening a voxel with scale k, whose volume is 4 p-k , produces 4 sub-voxels. We get at the end the Algorithm 3, using the function refine3D.

5.3.4. Costs. Similarly as the 2D analysis, we get Table 4 for the costs of Algorithm 3. Here again we will check a posteriori the value of N and we will compute the focus F , based on the number S of voxels computed at scales k 0 < k < p:

l , z )] k,j,l, =filteringMultiResol3D(k 0 , n, m)

Inputs: initial scale k 0 , sizes n = 2 p and m.

Read ∆t, ∆θ.

For every angle θ j , 1 j m, for every altitude z , 1 n, read [P (θ j , t l , z )] 1 l n , compute Pψ (θ j , t . , z ) = FFT[P (θ j , t l , z )]

invert the cut frequency data:

Output: parameters R, ∆t, ∆θ and downsampled ltered data PΩ k (θ (k)

Inputs: L p , L, L: lists of cells with scale p, < p, < p; L is sorted; ... From the rst cell (k, i, j, ,

Output: volume v of the rened voxel (and updated lists).

Algorithm 3 Multiresolution greedy algorithm for 3D reective tomography.

Inputs: Sizes n = 2 p , m; rate α ∈ (0, 1]; initial scale 1 k 0 < p.

Initialization

Lists of cells: the lists L and L p are void.

Reading and ltering: computation of the downsampled ltered data PΩ k (θ (k) j , t

Sort: sort the list L by decreasing intensity |H|. Volume to be lled: V = αn 3 .

Iterations

While the number of thin voxels is too small, i.e. V > 0. v = 0, L is void.

While the accumulated volume of the rened cells is too small, i.e. v < V , rene the most intense cell of L :

l , z )], ∆θ, R, ∆t).

EndWhile


Marianna Apidianaki

Xavier Tannier

Cécile Richart

Datasets for Aspect-Based Sentiment Analysis in French

Keywords: French Aspect-Based Sentiment Analysis, opinion mining, user reviews .ilsp.gr:8080

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

Aspect Based Sentiment Analysis (ABSA) is the task of mining and summarizing opinions from text about specific entities and their aspects [START_REF] Hu | Mining and Summarizing Customer Reviews[END_REF][START_REF] Ganu | Beyond the Stars: Improving Rating Predictions using Review Text Content[END_REF]. We have created and present two new datasets for the development and testing of ABSA systems for French which have been released in the framework of the SemEval-2016 Task 5 "Aspect-Based Sentiment Analysis". The first dataset comprises French restaurant user reviews annotated with relevant entities, aspects and polarity values. For this domain, both training and test data was provided. The second dataset is a smaller set of museum reviews for which no training data is available, dedicated to out-of-domain ABSA evaluation. The ABSA task was first introduced for English in the SemEval 2014 evaluation campaign [START_REF] Pontiki | SemEval-2014 Task 4: Aspect Based Sentiment Analysis[END_REF], where restaurant and laptop reviews annotated with aspect terms, categories and their polarity were provided for training and testing ABSA systems. The task was repeated in SemEval 2015 with a different, more unified, framework where aspect categories were defined as combinations of an entity type, an attribute type and a polarity value [START_REF] Pontiki | SemEval-2015 Task 12: Aspect Based Sentiment Analysis[END_REF]:

(1)

The fajitas were delicious, but expensive.

{FOOD#QUALITY, TARGET: fajitas}!POSITIVE {FOOD#PRICES, TARGET: fajitas}!NEGATIVE

(2) Great for a romantic evening.

{AMBIENCE#GENERAL, TARGET: NULL}! POSITIVE

An out-of-domain subtask was also proposed where annotated hotel reviews were provided for testing but no training data was released. In 2016, the SemEval ABSA task became multilingual [START_REF] Pontiki | SemEval-2016 Task 5: Aspect Based Sentiment Analysis[END_REF]. New datasets were released for English allowing systems to be tested on the same domains as in previous years (laptops, restaurants and hotels), but datasets were also developed in new languages using a common set of annotation guidelines. More precisely, the multilingual dataset created for this year's task comprises Restaurant reviews in six languages (English, French, Dutch, Russian, Spanish and Turkish), Hotel reviews in Arabic, Consumer Electronics reviews in three languages (English, Dutch and Chinese), Telecom reviews in Turkish and Museum reviews in French. All but the last dataset were released for Subtask 1 (Sentence-level ABSA) and part of the data was annotated at the text level for Subtask 2 (Textlevel ABSA). The French Museum reviews dataset was released for Subtask 3 "Out-of-domain ABSA" where only test and no training data was provided. 1 In what follows, we describe the data collection procedure and the annotation guidelines that were developed for the two domains addressed in French. The data and the annotation guidelines are publicly available under a noncommercial, no redistribution license 2 through META-SHARE [START_REF] Piperidis | The META-SHARE Language Resources Sharing Infrastructure: Principles, Challenges, Solutions[END_REF], 3 a repository devoted to the sharing and dissemination of language resources, and on the SemEval-2016 ABSA task website. 4

Datasets and Annotation

Data Collection

French datasets were developed for two of the SemEval-2016 ABSA subtasks. For in-domain sentence-level ABSA (Subtask 1) the dataset comprises annotated restaurant reviews while for out-of-domain ABSA (Subtask 3) annotated museum reviews were released. For the first subtask, the restaurant domain, both training and test data was provided. In Subtask 3, the participating teams had the opportunity to test their systems in a previously unseen The restaurant and museum user reviews were collected online and manually annotated. Participants in the French ABSA task were given a tool for downloading the unannotated reviews and were asked to assign the relevant annotations following the guidelines described in the next section. Each system's output was compared to the manually assigned ('gold') annotations during evaluation.

Annotation Schema and Guidelines for Restaurant Reviews

The French data was annotated by a native speaker linguistics student, then revised and corrected by the task organizers (a native and a non-native French speaker). Given a review text about a restaurant, the task of the annotator was to identify opinions expressed towards specific entities and their attributes and to assign the respective aspect category and polarity labels. The category values (entity-attribute pairs, E#A) had to be chosen from the inventory of entities and attributes provided for annotation of restaurant reviews in the SemEval-2015 task.

To increase the comparability of the multilingual datasets and the applicability of the systems, the same guidelines were applied for annotation in this domain in different languages. In the French guide, the English annotation examples were replaced with examples in French. In total, six entity labels are available in the restaurant domain and can be combined with specific attribute labels. The possible combinations are described in Table 2.

Annotation was done at the sentence level. If the sentence contained an explicit mention to the revised entity (e.g. a specific food or restaurant), the label was assigned to that specific occurrence of the expression in the text. Otherwise, if an entity was implicitly referred (for example, through pronouns) or inferred in a sentence, the label was assigned to the whole sentence and no specific target expression was marked in the text, as illustrated in the following examples.

(3) Carte des vins sympa. ) We have used the delivery service several times and have never been disappointed, but on-site service is horrible.

(5) Pourtant les plats sont bons et la déco est sympa.

{FOOD#QUALITY, TARGET: plats}!POSITIVE {AMBIENCE#GENERAL, TARGET: déco}!POSITIVE
) However, the food is good and the decoration is nice.

Cross-lingual Applicability of the Restaurant Guidelines

The ABSA-2016 task on Restaurant reviews involved six languages (English, French, Dutch, Russian, Spanish and Turkish). To ensure the comparability of the datasets and facilitate the cross-lingual applicability of the systems, the guidelines used for annotation were the same for all languages and corresponded to the guidelines provided for English in ABSA-2015 [START_REF] Pontiki | SemEval-2015 Task 12: Aspect Based Sentiment Analysis[END_REF]. This did not pose any particular problems in French at the level of category values, as the entities and attributes identified in the texts are similar in these closely related languages. We have however spotted the recurrent use of some expressions in the French reviews for which additional guidelines had to be given as these were not covered in the English guide. An example is the expression "rapport qualité/prix" in French (value for money), which assesses the quality of a restaurant or of the food/drinks served, compared to the cost. Sentences containing this expression are generally assigned two labels in the French dataset: one expressing the polarity regarding the prices (e.g. RESTAURANT#PRICES, FOOD#PRICES or DRINK#PRICES) and a second one describing the quality of the food/drinks, or the general quality of the restaurant.

To decide whether the entity addressed in a sentence is the restaurant or the food/drinks, it is often necessary to look at the preceding sentences in the text. In example 6 the entity referred to is the RESTAURANT, while in example 7 it is the FOOD because the preceding sentence in the text was: Perso j'avais pris le tournedos à 16 euros. ) Great value for money!

The MISCELLANEOUS label is assigned to attributes not covered by the other labels. In the French data, we have assigned this label to sentences describing the hygienic conditions in a restaurant, its style and the services provided (take-away, delivery, web reservations, etc.). (9) Déjà il y a un gros loupé niveau hygiène : le cuisinier sort dehors avec ses gants puis retourne en cuisine sans les changer...

{RESTAURANT#MISCELLANEOUS, TARGET: NULL}

!NEGATIVE

) There is a serious problem with hygiene: the cook comes out wearing his gloves and then goes back to the kitchen without changing them...

Although the applicability of the English category and attribute labels in French was quite straightforward, the identification of the entities referred to in the text turned out to be a not so easy task, especially in cases of multi-word expressions. This task is easier in English where dish names often correspond to nominal compounds, for example lobster ravioli, mashed potatoes or porcini mushroom pasta special. In French, dish names are often described by multi-word expressions containing prepositional phrases. The ambiguity of the French prepositions makes more complex the identification of the entity towards which an opinion is expressed in the text and the selection of the category attribute that is being addressed. In the following example, (10) Nous avons pris les pâtes au foie gras et cèpes, celles-ci baignaient dans de la crème et de la crème balsamique! {FOOD#QUALITY, TARGET: pâtes au foie gras et cèpes} !NEGATIVE

) We took the foie gras and porcini mushrooms pasta, which was floating in cream and balsamic cream! the prepositional phrase au foie gras et cèpes describes the ingredients of the recipe (foie gras and porcini mushrooms), so the whole expression pâtes au foie gras et cèpes is the target. In other cases, the prepositional phrase describes the preparation or serving mode and might be considered as a part of the target or not. For instance, poissons à la plancha (fish a la plancha) is a target of the category FOOD as it denotes a specific way of cooking fish with no polarity value, but in les frites maison sont à volonté (homemade french fries at will) the target is frites (french fries) and gets the label {FOOD#STYLE&OPTIONS ! POSITIVE} because the preparation (maison (homemade)) and serving mode (à volonté (at will)) are attributes with positive polarity.

Annotation Schema and Guidelines for Museum Reviews

The Museum domain had not been addressed in previous years or in other languages, so we had to compile a new set of guidelines for this domain. Annotated museum reviews were provided as the gold standard for evaluation in the outof-domain Subtask 3 where no training data was released. Six entity labels are available for this domain in the annotation guidelines and can be combined with specific attribute labels. Some of them are specific to the museum domain (e.g. TOUR GUIDING, COLLECTIONS) while others are also found in the restaurant domain (SERVICE, LOCATION). The possible entity-attribute combinations are described in Table 3.

The annotation for the out-of-domain substask was also done at the sentence level. Similarly to Subtask 1, if an explicit mention to the revised entity E (e.g. a museum) is contained in the text, the label is assigned to that specific occurrence of the expression. If an entity is implicitly referred or inferred, the label is assigned to the whole sentence, as shown in the following examples.

( ) Everything has been done to put forward all these treasures.

(13) béton brut partout aucun revêtement sur les murs... {MUSEUM#ARCHITECTURE, TARGET:

NULL} !NEGATIVE
) raw concrete everywhere, no coating on the walls...

Composition of the Datasets

In an effort to keep the training and test datasets as balanced as possible with respect to the polarity of the opinions expressed in the texts, we collected an equivalent number of user reviews with different rankings based on a 5-star rating system (excellent, very good, average, poor, terrible). As a consequence, the retained reviews are uniformly distributed over all possible ranks. Of course, a review may convey opposing sentiment for different aspect categories (i.e. E#A pairs). For example, a restaurant with good service might be considered too expensive or badly located, and a museum that houses interesting collections might provide poor facilities. Tables 4 and5 contain statistics by annotation type for each domain. We show the number of occurrences of an E#A pair with positive (POS), negative (NEG) or neutral (NEU) polarity in the data. For the Restaurant domain, we provide counts for the training and test sets separately and combined.

ABSA Evaluation at SemEval-2016

The SemEval-2016 ABSA evaluation is presented in detail in the task description paper [START_REF] Pontiki | SemEval-2016 Task 5: Aspect Based Sentiment Analysis[END_REF]. In this section, we provide the ranking of the systems that participated in the French ABSA task along with some information on the way that the gold review annotations were used for evaluation. More information about the systems can be found in the system description papers of the SemEval task.

In the sentence-level subtasks (1 & 3), given an unannotated review text, the participating systems had to identify tuples containing three types of information:

• Slot 1: the aspect category, i.e. an entity and attribute (E#A) pair towards which an opinion is expressed in the text

• Slot 2: an opinion target expression possibly used in the text to refer to the reviewed entity E of a pair E#A

P = {positive, negative, neutral}

In Slot 1, the evaluation assesses whether a system correctly identifies the aspect categories towards which an opinion is expressed in the texts. The categories returned by a system are compared to the corresponding gold annotations and precision (P), recall (R) and F-1 scores are calculated. In Slot 2, the evaluation assesses whether a system correctly identifies the expressions used in a sentence to refer to the reviewed entities. The returned list of target expressions is compared to the gold annotations to calculate P, R and F-1 scores. For Slot 3, the total accuracy score for polarity classification is retained. System performance for all slots is compared to a baseline score. The baseline system for Subtask 1 selects the categories and polarity values to be assigned to the test sentences using a Support Vector Machine (SVM) based on bag-of-words (BoW) features and trained on the manually annotated training data. Participants were free to decide the domain(s), subtask(s) and slot(s) they wished to participate in. Five systems participated in Subtask 1 for French. The out-of-domain Subtask 3 had no participants. However, the data developed for this subtask remains available and can be used for future research. The results in F-1 score (Slots 1 and 2) and the accuracy scores (Slot 3) obtained in Subtask 1 are reported in Table 6. The submissions were marked as constrained (C) or unconstrained (U) depending on whether the participants used only the provided training data or if they also exploited information from additional resources (e.g. lexicons or additional training data). We give in boldface the score of the system that performed best in the evaluation for each slot. The baseline systems and evaluation scripts are available for download from the SemEval-2106 ABSA website. More details on the evaluation procedure and the results are provided in the task description paper [START_REF] Pontiki | SemEval-2016 Task 5: Aspect Based Sentiment Analysis[END_REF].

Conclusion

We have presented two new datasets for Aspect-Based Sentiment Analysis in French. These are the first French datasets annotated with this type of information and were initially developed and released in the framework of the multilingual SemEval-2016 ABSA task. Two subtasks were addressed, namely in-domain sentence-level ABSA and out-of-domain ABSA (Subtasks 1 and 3) where usergenerated reviews in two domains were provided. Altough the datasets were released for training and evaluation in the framework of the SemEval ABSA task, they remain available and can be used for research purposes.

Table 2 :

 2 Possible entity-attribute pairs in the Restaurant domain.

	place est déplorable.

{RESTAURANT#MISC, TARGET: NULL}!POSITIVE {SERVICE#GENERAL, TARGET: service}!NEGATIVE

Table 3 :

 3 Possible entity-attribute (E#A) combinations in the Museum domain.

	appropriate label sometimes requires taking into account
	information at the text level.
	(6)	Un moment agréable avec des amis pour un rapport
		qualité/prix imbattable.
		{RESTAURANT#GENERAL,

TARGET: NULL}

!POSITIVE

{RESTAURANT#PRICES, TARGET: NULL}!POSITIVE ) A good time with friends and a great value for money.

(7) Perso j'avais pris le tournedos à 16 euros.

) Me, I took the tournedos at 16 euros.

Très bon rapport qualité prix ! {FOOD#PRICES, TARGET: NULL}!POSITIVE {FOOD#QUALITY, TARGET: NULL}!POSITIVE

Table 4 :

 4 Statistics by annotation type in the Restaurant domain.

	Entity	Attribute		Train					Test		Train
			POS NEG NEU Total POS NEG NEU Total + Test
		GENERAL	175	214	14	403	58	63	4	125	528
	RESTAURANT	PRICES	30	60	12	102	9	10	4	23	125
		MISCELLANEOUS	27	62	17	106	9	19	6	34	140
		QUALITY	332	340	43	715	137	128	22	287 1002
	FOOD	PRICES	31	45	10	86	6	23	6	35	121
		STYLE & OPTIONS 118	97	30	245	49	45	15	109	354
		QUALITY	28	20	3	51	15	5	1	21	72
	DRINK	PRICES	2	14	1	17	1	4	1	6	23
		STYLE & OPTIONS	12	9	1	22	4	6	3	13	35
	SERVICE	GENERAL	218	305	14	537	70	111	6	187	724
	LOCATION	GENERAL	34	3	3	40	20	1	4	25	65
	AMBIENCE	GENERAL	157	43	6	206	63	18	7	88	294
		Entity		Attribute		POS NEG NEU Total
			GENERAL		122	64	15	201
				PRICES		16	24	1	41
		MUSEUM	COMFORT ACTIVITIES		13 4	27 3	6 2	46 9
			ARCHITECTURE		47	11	4	62
			MISCELLANEOUS	10	24	7	41
		LOCATION	GENERAL		24	2	1	27
		SERVICE	GENERAL		11	25	0	36
		FACILITIES	GENERAL PRICES		14 1	16 6	4 0	34 7
		TOUR GUIDING	GENERAL		9	4	3	16
		COLLECTIONS	SETUP INTEREST		28 172	62 81	6 22	96 275

11) Très beau musée. {MUSEUM#GENERAL, TARGET: musée} !POSITIVE ) Very nice museum.

(12)

Tout est mis en oeuvre pour mettre l'ensemble de ces trésors bien en valeur.

Table 5 :

 5 Statistics by annotation type in the Museum domain.

	{COLLECTIONS#SET UP, TARGET: trésors} !POSITIVE {COLLECTIONS#INTEREST, TARGET: trésors} !POSITIVE

Table 6 :

 6 Results of the SemEval-2016 French ABSA Subtask 1.

	System	Slot 1 (F-1) Slot 2 (F-1) Slot 1&2 (F-1) Slot 3 (Acc.)
	XRCE/C	61.207	65.316	47.721	78.826
	IIT-TUDA/U	57.875	66.667	-	72.222
	IIT-TUDA/C	57.033	-	-	-
	INSIGHT-1/C	53.592	-	-	73.166
	UFAL/U	49.928	-	-	-
	UWB/C (1st run)	-	-	-	75.262
	UWB/C (2nd run)	-	-	-	74.319
	baseline/C	52.609	45.455	33.017	67.4

• Slot 3: the sentiment polarity to be assigned to each identified <category, target> pair, from the set