Pierre-Loup Schaefer

Grégory Chagnon

Alexandre Moreau-Gaudry

Advanced Sensors Placement for Accurate 3D Needle Shape Reconstruction

Needles are tools widely used in minimally invasive surgery. During such procedure the localization of the needle and its tip is a challenging situation because of the needle deformations due to the interactions with tissues. To tackle this problem, instrumented needles with sensors are currently developed to allow needle reconstruction and tip localization. In conventional surgery this difficulty is overcome by medical imaging. The interest of using an instrumented needle resides in the possible dispense of medical imaging. This papers develops new methods to reconstruct needles in three dimensions and to find the locations of sensors which minimizes the error of reconstruction of the needle. A notable feature of our method is that input data are based on real needle data, that should assure a better representativity of our results. Reconstructions simulated with 22 gauges 200mm long needles data show that the localization of the needle tip is more accurate by 18% to 52% with optimal sensors positions compared to equidistant sensors positions.

I. INTRODUCTION

Needles are more and more used in percutaneous medical procedures, such as, for instance, during interventional Radiology Procedures. During these procedures, the knowledge of the needle location is of primary importance with stakes being no others than the success or the failure of the procedures. Nevertheless, accurate location of the needle may be challenging because of the bending phenomenon which might appear when the needle is inserted. This phenomenon is caused by interactions between needle and tissue. It is also influenced by a lot of factors such as the size of the needle, its length or the type of tissues. Insertions in phantom using 18 gauges 20 cm have shown that tip deflection can be up to 2.8mm for a 6cm insertion [START_REF] Wan | Brachytherapy needle deflection evaluation and correction[END_REF] and up to 12 mm for a 10 cm insertion [START_REF] Podder | Effects of tip geometry of surgical needle an assessment of force and deflection[END_REF]. Furthermore deflections of 22 gauges 20 cm needles used for insertions in pig tissues can be as high as 25 mm [START_REF] Adeline | A generic three-dimensional static force distribution basis for a medical needle inserted into soft tissue[END_REF]. As a matter of fact, while using needles during minimally-invasive medical procedures, with no assisting tools, incertainty on the position remains. Furthermore current computer assisted interventional radiology systems make the assumption of non-deformability of the needle [START_REF] Carrat | Device and method for determining the position of an instrument in relation to medical images[END_REF], [START_REF] Leroy | System and method for imaging and locating punctures under prostatic echography[END_REF] and then don't provide direct access to the real shape of the needle, still needing additional medical imaging.

This article presents a method to reconstruct the 3D shape of a needle from the data retrieved from its sensors and a procedure to find the best locations of the sensors to minimize the tip localization error. The paper is organized as follows: Section II explains the three-dimensional reconstruction of the needle shape from sensors data. Section III describes the *This work was supported by the French ANR GAME-D 1 Univ. Grenoble Alpes / CNRS / TIMC-IMAG UMR 5525, France. Email: FirstName.Lastname@imag.fr method for evaluating best sensors position by processing real needle data. Section IV present and analyze the results, and is followed by Section V on future improvements.

II. 3D NEEDLE SHAPE RECONSTRUCTION A. Sensors informations

This section presents the process to obtain geometrical informations of the needle from the sensors. The purpose of the sensors is to measure the strain applied on the needle during the insertion. Differents types of sensors include optical fiber Bragg gratings [START_REF] Abayazid | 3d flexible needle steering in soft-tissue phantoms using fiber bragg grating sensors[END_REF], [START_REF] Park | Realtime estimation of 3-d needle shape and deflection for mri-guided interventions[END_REF] and strain gauges [START_REF] Yang | Conception et intgration de microsystmes sur un cylindre pour la mesure de ses dformations : application un outil du domaine de la sant[END_REF], [START_REF] Yang | Modelling and characterization of an instrumented medical needle in sight of new microsensor design for its insertion guidance[END_REF]. From a sensor triplet it is possible to obtain two geometrical parameters at the sensor location: the curvature κ and the bending angle θ [START_REF] Abayazid | 3d flexible needle steering in soft-tissue phantoms using fiber bragg grating sensors[END_REF], [START_REF] Roesthuis | Three-dimensional needle shape reconstruction using an array of fiber bragg grating sensors[END_REF]. Let (s i ) i=1,..,n the value such as (M(s i )) i=1,..,n are points of the needle where the triplets are located. Then for the triplet located at s i we can retrieve the curvature κ i = κ(s i ) and the bending angle θ i = θ (s i ). The set (κ(s i )) i=1,..,n and (θ (s i )) i=1,..,n constitutes a discretization of the function κ and θ . The interpolations of both sets gives us κ est and θ est , the estimate of curvature and bending angle from the sensors data.

B. Shape reconstruction

This section presents the method to obtain 3D needle shape estimation from the curvature and bending angle estimate κ est and θ est . Due to the geometry of the needles and the mechanical interactions with the tissues we make the widely accepted hypothesis that the mechanical torsion of needles is negligible [START_REF] Abayazid | 3d flexible needle steering in soft-tissue phantoms using fiber bragg grating sensors[END_REF], [START_REF] Park | Realtime estimation of 3-d needle shape and deflection for mri-guided interventions[END_REF], [START_REF] Yi | Spatial shape reconstruction using orthogonal fiber bragg grating sensor array[END_REF]. Given this assumption our reconstruction will be based on Rotation Minimizing Frame. Rotation Minimizing Frame (RMF) belongs to the family of 3D moving frame. A 3D moving frame is a frame composed of an orthonormal ordered basis vector defined along a space curve. In the rest of the study we use the notation (T, N 1 , N 2 ) to refer to the RMF. We have T, N 1 , N 2 ∈ R 3 and T is tangent to the curve. Differents moving frames can be defined for the same curve [START_REF] Bishop | There is more than one way to frame a curve[END_REF] and RMF has the property of minimizing its twisting motion when framing a curve [START_REF] Wang | Computation of rotation minimizing frames[END_REF]. Combined with the negligible torsion hypothesis this property results in the RMF being a materially fixed basis of the beam [START_REF] Rida | Rational approximation schemes for rotation-minimizing frames on pythagorean-hodograph curves[END_REF]. Calculus of the Cartan matrix under such assumptions give the following relationship:

d ds   T(s) N 1 (s) N 2 (s)   = κ(s)   0 cosθ (s) -sinθ (s) -cosθ (s) 0 0 sinθ (s) 0 0     T(s) N 1 (s) N 2 (s)   (1)
The signification of Eq.1 is that for every point of the needle M(s) with s ∈ [0, L], where L is the length of the needle, there is a relation at this point between the material frame (T(s), N 1 (s), N 2 (s)) and its derivative involving κ(s) and θ (s). Assuming known κ est and θ est , the estimate of curvature and bending angle from the sensors, and under appropriate boundaries conditions, the resolution of differential equation system Eq.1 gives the estimate of the RMF of the needle (T est , N1 est , N2 est ). Finally the estimate of the needle shape M est is obtained by integrating T est :

M est (s) = M 0 + s s 0 T est (x)dx, s ∈ [0, L] (2)

III. SENSORS PLACEMENT

A. Method

Sensor placement is a relatively few documented subject with [START_REF] Park | Realtime estimation of 3-d needle shape and deflection for mri-guided interventions[END_REF] and [START_REF] Adeline | A generic three-dimensional static force distribution basis for a medical needle inserted into soft tissue[END_REF] being the main articles on the subject. Approches made in both articles are to consider the optimal sensor location problem as a minimization one. The optimal sensor location is then defined as the location which minimizes the error made when reconstructing shapes with sensors at a certain location. Process used for evaluating the error of reconstruction on a needle for a given sensor location has the following steps. Firstly is defined a set of reference needles as input data. In both cases, these needles are calculated with 2D beam theory using force loads. Loads used by Park and al. [START_REF] Park | Realtime estimation of 3-d needle shape and deflection for mri-guided interventions[END_REF] are constitued of equidistantly reparted 2D forces over the needle with arbitrary amplitudes, while loads used by Robert and al. [START_REF] Adeline | A generic three-dimensional static force distribution basis for a medical needle inserted into soft tissue[END_REF] are more representative as the 2x2D forces comes from a preliminary study conducted on real needles. Then a calculation is made for given sensors locations, representing what the sensors output would be if they were really present at that place on each of the needle of the reference set. Under the 2D assumption the output of the sensors take the form of simple curvature. Eventually these values are used to reconstruct a needle shape which is then compared to the reference needle. Quality measurement of the reconstruction is made by using the tip deflection error, a widely used medical criteria for needle precision when performing a surgical act.

The method of sensor location optimization presented in this article differs from these previous methods by the use of a reference set based on real needle data and by the three-dimensional procedures, resulting in a more exhaustive approach. Fig. 2

B. Study of a set of needle

The aim of this section is to generate the reference set of needles which will be used in sensor location optimization process. We retrieve the needle shape from real needle data by B-Spline smoothing and then we compute the RMF of the B-Spline. The RMF will be used to simulate the data of the sensors (Fig. 2: Real 3D needle RMF).

1) Data: The study was carried out on a set of 54 CT scans of needle insertion into pig shoulder [START_REF] Adeline | A generic three-dimensional static force distribution basis for a medical needle inserted into soft tissue[END_REF]. Experiments involved standart 200 mm long stainless steel 22 gauges needles, which are often used in interventional radiology procedures. Choice of pig tissues is due to its similar biomechanical properties with human tissues [START_REF] Hocking | A review of the benefits and pitfalls of phantoms in ultrasound-guided regional anesthesia[END_REF].

2) Segmentation: In order to retrieve the voxels representing the needle we performed a segmentation on each of the scan. The aim is to separate the needle from the surrounding tissues, bones, etc... A region growing algorithm is particularly indicated in that case because of the difference of density of the needle and the density of the organs and tissues. We implemented the algorithm in a plug-in of the medical visualisation software CamiTK [START_REF] Fouard | Camitk: a modular framework integrating visualization, image processing and biomechanical modeling[END_REF]. Once the segmentation is processed stays only the voxels of the needle. We then have 54 data sets of voxels.

3) Smoothing: The main objective of this part is to reconstruct the shape of the needle from its voxels. The mathematical representation of the shape is the three-dimensional univariate mathematical function which maps the length between the proximal extremity and one point of the needle to the three-dimensional coordinates of this point. Needle voxels are the results of a noisy spatial discretization of that function, consequently it is possible to approximate the shape function by performing smoothing on needle voxels [START_REF] Christian | Smoothing by spline functions[END_REF]. The smoothing used here will be B-Spline smoothing, one of the most commonly used methods to smooth data. It consists to perform a regularized regression on a B-Spline basis. 3D B-Splines are piecewise polynomials defined by their order n, their knots (t i ) i=1,..,m and their control points (P i ) i=1,..,n c ∈ R 3 [START_REF] De | A Practical Guide to Splines[END_REF]. Let S(t) be the B-Spline which is defined by:

S(t) = n c ∑ i=1 B i,n (t)P i (3)
with the following recursive definition:

B i,0 (t) = 1 i f t i ≤ t ≤ t i+1 0 otherwise (4) B i,n (t) = t -t i t i+n -t i B i,n-1 (t) + t i+n+1 -t t i+n+1 -t i+1 B i+1,n-1 (t) (5)
Let the data points (Y i ) i=1,..,n d be the voxels coordinates and (w i ) i=1,..,n d be their values. Let (t i ) i=1,..,n d be the parameterization associed with the data points. The smoothing system then takes the following form:

n d ∑ l=1 w l S( tl ) -Y l 2 + λ 1 0 S (t) 2 dt (6)
The smoothing B-Spline estimate is defined to be the minimizer of the Eq.6 over the control points (P i ) i=1,..,n c . The left part of Eq. 6 reflects accuracy to the data points while the right part express the roughness of the estimate. The trade-off between these two parts is the variable λ , called the smoothing parameter. Because the problem of smoothing noisy data with a spline is well posed, the system admits one and only one solution. Thus, as a result, we have a set of 54 3D B-Spline. 4) Rotation Minimizing Frame: Each B-Spline estimate is then processed to determine its RMF. Method used to compute the RMF is called Double reflection method [START_REF] Wang | Computation of rotation minimizing frames[END_REF]. Choice of this method in particular is justified by its stability and its global approximation error of high order. Finally we obtain a set of 54 RMF corresponding to the material frames of the needles scanned in the input data.

C. Minimization

Finding optimal sensors positions is a n-dimensional problem, with n being the number of sensors. It belongs to the Mixed Integer Nonlinear Programming (MINLP) family as it is a problem with a nonlinear objective function (the reconstruction error) and that we restrain the sensors positions to integer (precision on the positions of the sensors doesn't need to be submilimeter). As the convexity of our problem is undetermined we used global optimization solvers (Genetic Algorithm [START_REF]MATLAB and Global Optimization Toolbox. Release R2014a[END_REF] and NOMAD [START_REF] Sébastien | Algorithm 909: Nomad: Nonlinear optimization with the mads algorithm[END_REF]).

IV. RESULTS

This section presents results from our study on optimal sensor locations on 22 gauges 200 mm long needles set. I presents the result of the minimization problem presented in Section III. It shows the optimal location of the sensors for a number of sensor going from one to seven. The optimization solvers didn't converge for a number of sensors superior to seven. Fig. 3 shows the cumulated tip deflection error from the reconstruction of our set of needles with two sensors. The optimal sensors position is the one for which the minimum reconstruction error is obtained. We see that for two sensors the reconstruction error function is smooth and almost convex. Left: Tip deflection error for equidistant and optimal sensors positions depending to the number of sensors. Right: Percentage of precision gained by using optimal sensors positions instead of equidistant sensors position Fig. 4 compares the mean tip error of reconstruction of our set with optimal sensors placement sensors versus equidistant sensors placement. It indicates also the percentage of precision gain by using needle with optimal sensors placement instead of equidistant sensors placement. We can see that the percentage of precision gained is comprised between 18% (2 sensors) to 52% (5 sensors). this range of number of sensors. The question of optimal locations for a higher number of sensors is not relevant as the number of sensors on a needle is limited. Moreover we see on Fig. 4 that the gain between optimal sensor location and equidistant sensor location starts to decrease for a number of sensor superior to five. The reason is that when the number of sensors increase they tend to be placed in an equidistant way. Fig. 4 presents the mean tip deflection error after reconstruction of our reference set of needle. Comparison with other works shows that our mean tip error is correct compared to the length of our needles: for a 4 sensors needle we have a mean tip error of 3.2 mm for 200mm long needles while Abayazid et al. have 1.8 mm for 120mm needles [START_REF] Abayazid | 3d flexible needle steering in soft-tissue phantoms using fiber bragg grating sensors[END_REF].

Here the ratio of mean tip error over the length of the needle is approximately the same. The improvement of the quality of the reconstruction evaluated in the Fig. 4 indicates that the reconstruction error falls from 18% (2 sensors) to 52% (5 sensors) when optimal sensors locations are used compared to equidistant sensors locations. This results exhibits that the gain of precision coming from an optimal sensors placement is not negligible. It confirms also that for a number of sensors going from one to seven our method gives better locations than a standart equidistant placement. Fig. 5 expose the concrete case of one of the needle of our reference set which has a huge deformation (tip deflection superior to 20 mm). Thus the reconstruction with only two sensors placed equidistantly gives a tip error of 9.6 mm. Keeping the same number of sensors but with optimized locations divide the tip error by 2 (4.8 mm) which is a little bigger than the the tip error with five sensors (3.9 mm). This example shows here that placing optimally the sensors can be a relevant solution to increase accuracy and should be considered before adding new sensors.

The use of real needle data combined with three dimensional processing contributes to the relevancy of our method by using less restrictive hypothesis. We believe that these choices improved the accuracy of reconstruction and the reliability of the sensors positions results.

Fig. 1 .

 1 Fig. 1. Flowchart indicating the steps of the 3D Needle Shape Reconstruction using the sensors measures.

Fig. 3 .

 3 Fig. 3. Cumulated squared tip deflection error of reconstruction on our set of needles depending on differents sensors locations (logarithmic scale).

Fig. 4 .

 4 Fig. 4.Left: Tip deflection error for equidistant and optimal sensors positions depending to the number of sensors. Right: Percentage of precision gained by using optimal sensors positions instead of equidistant sensors position

Fig. 5 Fig. 5 .

 55 Fig.4compares the mean tip error of reconstruction of our set with optimal sensors placement sensors versus equidistant sensors placement. It indicates also the percentage of precision gain by using needle with optimal sensors placement instead of equidistant sensors placement. We can see that the percentage of precision gained is comprised between 18% (2 sensors) to 52% (5 sensors).Fig.5presents an exemple of reconstrusction of a needle of our data set with differents number of sensors and differents placement.V. DISCUSSIONS AND CONCLUSION Fig.Ipresents the positions for which the solvers converged towards the same solution for a number of sensors comprised between one and seven. As these sensors positions minimize the reconstruction error of our set of needles we can assert that we found the optimal sensors locations for

  presents an overview of the differents steps of our method.

		Real 3D
		needle RMF
	Sensors locations	Sensors measures emulation	Sensors measures	3D Needle Shape
	Reconstruction error	Shape Comparison	3D Needle shape estimate	Recon-struction
	Fig. 2. Flowchart illustrating the evaluation of the reconstruction error
	from sensors location over a set of real needles

TABLE I TIP

 I DEFLECTION ERROR DEPENDING TO THE NUMBER OF SENSORS

	Table


Mathieu Lagrange
email: mathieu.lagrange@cnrs.fr

Mathias Rossignol

Grégoire Lafay

Visualization of audio data using stacked graphs

INTRODUCTION

The visual display of quantitative information [START_REF] Tufte | The Visual Display of Quantitative Information[END_REF] is at the core of the growth of human knowledge as it allows human beings to go beyond the limitation of natural languages in terms of precision and scale.

Defining what is the essence of a good visual display of quantitative data is non trivial and usually domain specific. That said, in most scientific fields, such displays serve two majors goals: 1) the routine interaction of the researcher with the data or the physical phenomenon and 2) the need of the researcher to motivate its claim to its peers. Both tasks require the display to fulfill the simplicity rule both in terms of production and design. First, the display shall be computed and adapted according to the need of the researcher very efficiently in order to allow an effective exploration of the data. Second, the display shall be able to convey at the first glance an important qualitative aspect about the data.

This paper is about the visualization of audio data, and audio data is originally made to be listened to. Therefore, we shall keep in mind that "all visual projections of sounds are arbitrary and fictitious" [START_REF] Schafer | The soundscape: Our sonic environment and the tuning of the world. Destiny books[END_REF]. That said, even if recorded versions of sounds can now be played back c Mathieu Lagrange * , Mathias Rossignol, Grégoire Lafay * . Licensed under a Creative Commons Attribution 4.0 International License (CC BY 4.0). Attribution: Mathieu Lagrange * , Mathias Rossignol, Grégoire Lafay * . "Visualization of audio data using stacked graphs", 19th International Society for Music Information Retrieval Conference, Paris, France, 2018. at convenience, it is still useful to represent them graphically as listening depends on time. On contrary, the visual display allows the reader to grasp a global view of the waveform at a glance. Also, the eye is less subject to stimulation fatigue and the visual display is very powerful to convey evidence as we are still fully into the print culture that since the Gutenberg invention gives an "uncritical acceptance [to] visual metaphors and models" [START_REF] Mcluhan | The Gutenberg Galaxy[END_REF].

We propose in this paper a display of audio data that is intuitive and gives information about the main dimensions of sound in a compact manner using stacked graphs [START_REF] Byron | Stacked Graphs-Geometry & Aesthetics[END_REF]. The display can be computed easily and efficiently 1 . In order to put this display into context, an overview of the routinely used type of displays is given, respectively from the perspective of the musician composer in Section 2 and the physicist in Section 3. We shall argue that the proposed display fully described in Section 4 can be thought of as the physicist's counterpart to a notational system introduced by Schafer [START_REF] Schafer | The soundscape: Our sonic environment and the tuning of the world. Destiny books[END_REF].

The display is then evaluated and compared to the commonly used waveform and spectrogram displays with two perceptual experiments. In the first experiment, the subjects have to distinguish between tones of different musical instruments by listening to the sounds or by considering the visual displays under evaluation. The protocol and the results for this experiment are presented in Section 5. In the second experiment, the subjects are asked to distinguish between saxophone performances of different level of instrumental expertise by listening to the sound and by considering the displays under evaluation. The protocol and the results for this experiment are presented in Section 6.

ABOUT NOTATION

From the phonetic alphabet for speech to the musical score for music, notation consists in putting together on a one or two-dimensional space symbols describing specific sound events. In a manner probably inherited from writing, time sequencing is usually depicted from left to right in the Western musical culture. Specific to the musical score is the use of the vertical axis to depict the pitch. A musical tone is therefore solely described in terms of time of appearance, duration, pitch and sometimes intensity. As such, the score is largely prescriptive and gives a tremendous amount of freedom to the musical performer in terms of interpretation.

In an intent to provide a more descriptive notation of musical objects, Schaeffer [START_REF] Schaeffer | Traité des objets musicaux[END_REF] designed a "solfège des object musicaux" that extensively apprehend the description of any kind of sound object. Perhaps because of its complexity this notation is hardly used today. In an effort to simplify this notation, Schafer proposed a notational system that can be considered for describing any kind of sound, be it a unique event or any kind of compound. The main rationale is to split the temporal axis from left to right into 3 parts corresponding to the attack, sustain and decay. For each part, its duration, frequency (related to the notion of mass as introduced by Schaeffer), fluctuations (related to the notion of grain as introduced by Schaeffer) and dynamics are displayed from top to bottom. Except for the frequency content that is depicted as a rough spectrogram contour, the other dimensions are described according to a specific alphabet of a few symbols. An example taken from [START_REF] Schafer | The soundscape: Our sonic environment and the tuning of the world. Destiny books[END_REF] of such annotation is given on Figure 1 for the sound of a church bell.

ABOUT MEASURE

When dealing with sound as a physicist, one wants to quantify mechanical properties and display them precisely. As in notation, the main aspect that is commonly looked for is the distribution of energy across frequency and time. The distribution of energy as a function of the modulation rate and the frequency scale of observations are less considered in the signal processing literature [START_REF] Anden | Multiscale Scattering for Audio Classification[END_REF][START_REF] Chi | Multiresolution spectrotemporal analysis of complex sounds[END_REF] but are shown to be perceptually important [START_REF] Dau | Modeling auditory processing of amplitude modulation. i. detection and masking with narrow-band carriers[END_REF][START_REF] Yang | Auditory representations of acoustic signals[END_REF].

Therefore, in order to display a sound on a twodimensional plane, one has to resort to a choice or a compromise. Either timing is emphasized and frequency neglected as in the waveform display 2a or frequency is emphasized and timing neglected as in the display of the Fourier spectrum 2b. A compromise is made by considering time and frequency respectively as horizontal and vertical axes of the two-dimensional plane as with the popular spectrogram magnitude of the short term Fourier transform, see Figure 2c. In such display, the use of a color code conveys information about energy.

That said, we believe that the spectrogram display still favors frequency over time. Spectral structure can be analyzed precisely, for example harmonicity, modulations, etc. Conversely, temporal dynamics and structure are harder to appreciate, as the way energy fluctuates in each sub bands has to be reconstructed from the color code.

The spectrogram is a display that is thus in our opinion very powerful for close inspection of a sound event that is active over a short period of time. Indeed, enlarging the time resolution quickly blurs the frequency resolution and may lead to a completely non informative display.

VISUALIZING SPECTRAL CONTENT USING STACKED GRAPH

With those limitations in mind, we propose in this paper to take a compromise that conversely favors time over frequency. In such display, the plane is therefore organized with time and energy as the horizontal and vertical axes respectively. The frequency is displayed as stacked layers displaying the level of energy across frequency sub bands of growing frequency range. Those layers can have colors assigned.

We seek a display that depicts information that is per- ceptually meaningful. Therefore, we consider spectral data projected on a Mel-scale [START_REF] Ss Stevens | A scale for the measurement of the psychological magnitude pitch[END_REF].

In order to improve legibility, colors are assigned to frequency layers according to their ranges with a color code ranging from blue (low frequency) to yellow (high frequency). The blue color is often associated with large phenomena, with the following adjectives: celestial, calm, deep, whereas the yellow color is often associated with transient phenomena that are highly energetic. Kandinsky in [START_REF] Kandinsky | Concerning the spiritual in art[END_REF] states that "Blue is comparable to low pitched organ sounds. Yellow becomes high pitched and can not be very deep". The color code is then chosen to be a linear gradient from blue (low frequency range) through green (middle frequency range) to yellow (high frequency range). In this paper, the gradient follows the LCH color model specified by the Commission Internationale de l' Éclairage (CIE) so that the perceived brightness appears to change uniformly across the gradient while maintaining the color saturation.

We argue that this display, termed spectral stack (spack), convey useful information about the sound. In particular, it conveys nicely, aside of fine details, the important dimensions retained by Schafer, see Figure 1.

To compute the spack display, a mel-scaled magnitude spectrogram is computed from the audio, see Figure 3. To each mel spectral band is assigned a given color code from dark blue (low frequency) to yellow (high frequency). At each time frame, the spack display is a stacking of the magnitude values of each mel frequency band, see Figure 4.

TASK 1: IDENTIFYING THE MUSICAL INSTRUMENT

The identification of the musical instrument used to play a tone rely largely on 2 factors, the spectral envelope and the attack [START_REF] Trevor R Agus | Fast recognition of musical sounds based on timbre[END_REF][START_REF] John | Multidimensional perceptual scaling of musical timbres[END_REF]. The spack display shall be able to conveniently display those factors. Indeed, the spectral envelope, i.e. the distribution of the energy across frequency is encoded using the stacking axis and color code. The attack is also well displayed as the spack focuses on the display of energy through time.

Protocol

Several tones played by four musical instruments: piano, violin, trumpet, and flute are considered as stimuli. Each instrument is played mezzo forte at 5 different pitches: C, D, E, F and G. For each sound, three visual representations are evaluated: waveform (W), spectrogram (Spe) and spack (Spa). For reference, the sound (S) is also considered 2 .

The test is a forced-choice categorization task. The sounds are displayed by gray dots on a 2 dimensional plane displayed on a computer screen. The dots can be moved freely within this plane and colored using 4 different colors, each corresponding to a given instrument. The correspondence is given to the subjects at the beginning of the experiment by the instructor: piano (black), violin (red), trumpet (magenta), and flute (green). If the sound modality is tested, the sound is played when the dot is clicked. If a visual modality is tested, the corresponding display is shown when the dot is clicked using the mouse.

Eight subjects, studying at the Engineering school "Ecole Centrale de Nantes", aged from 24 to 26 years, performed the test. Each subject reported normal hearing. They performed the test at the same time in a quiet environment using headphones. The sound level was set to a comfortable level before the experiment. A short introduction was given by the instructor for each display with a focus on the meaning of the axes and the color code. The subjects performed the evaluation using the sound modality first. The order of the three remaining modalities are ordered randomly among subjects to reduce the impact of precedence. The test is over when the subjects have assigned a color to each dot, this for all the evaluated modalities.

Results

Classification performance is evaluated as the number of couple of sounds played by the same instrument that have been assigned the same color divided by the number of couples. As can be seen on Figure 5, the task is trivial when listening to the sound, as the subjects achieve perfect classification. On overall, the classification is quite good for each of the graphical displays with a higher average performance for the spack display. Subjects verbally reported ease of use for the spack display.

TASK 2: ASSESSING THE LEVEL OF A SAXOPHONE PERFORMANCE

The control of the breath while playing the saxophone is crucial and can be monitored to assess the technical level of a saxophone player [START_REF] Robine | Evaluation of the technical leval of saxophone performers by considering the evolution of spectral parameters of the sound[END_REF]. For example, playing a single tone with sharp attack and constant amplitude during the steady state is non trivial and requires years of practice. Professional players typically practice such exercises on a daily basis as warm-ups and perform them with a trainer to get criticisms in order to improve their skills. Using graphical displays of their performance could be useful for them to spot during or after the performance. In order to be efficient, such display shall be intuitive with a few degrees of freedom in order to be easy to understand.

The validation of the spack display for such pedagogical needs is out of the scope of this paper. Nonetheless, we designed here a task that can demonstrate how several meaningful characteristics of the saxophone performance can be identified only by considering the graphical displays under evaluation.

In this kind of training, it could be useful for the trainer to have some kind of display of its performance. As the crucial part is to be able to control the air flow while playing in order to keep a stable amplitude and timbre, we hy- pothesize that the spack display may be a good candidate for such a task.

Protocol

The stimuli considered in this experiment are recorded performances of four saxophone players with a technical level assumed to be high or low (2 low, 2 high). Each player played several tones at pitch B and G. They were asked to play each note in three different ways: piano, forte and crescendo decrescendo 3 .

The test follows a XXY structure, where three performances are shown to the subject, one is at a given level (high or low) and the other two of the other level (low or high). The subject is then asked, based solely on the modality at hand, to select the one that is different from the two others. 24 triplets are randomly selected from the valid combinations of the above described stimuli. 16 subjects, studying at the Engineering school "Ecole Centrale de Nantes", aged from 24 to 28 years, performed the test in two sessions, 9 for the first session, and 7 for the second session. Each subject reported normal hearing. For each session, they performed the test at the same time in a quiet environment using headphones. The sound level was set to a comfortable level before the experiment. A short introduction was given by the instructor for each display with a focus on the meaning of the axis and the color code. The subjects performed the evaluation using the sound modality first. The order of the three remaining modalities are ordered randomly among subjects to reduce the impact of precedence. The test is over when the subjects have examined the 24 triplets for the 4 evaluated modalities. For each modality, the number of correct selection is averaged among the 24 triplets and then averaged among subjects. As can be seen on Figure 8, the task is more complex than task 1, as the score achieved using the sound modality is lower than task 1. This might be due to the fact that the task is less explicit than task 1. For the visual displays, the same ranking as task 1 is observed with a larger difference between each modality. A repeated measure ANOVA is used to test the potential significance of the type of display on the differentiation performance. A mauchly test reveals that the default of sphericity is not significant, thus no correction of the degrees of freedom of the Fisher test is needed. Table 1 presents the results of the Fisher test showing that the effect of the representation is significant p = 0.003. In addition, a multiple comparison test shows that the only significant differences are between Waveform and Spack p = 0.03 and Waveform and Sound p = 0.003. No significant difference is found between the remaining modalities: the Sound, the Spectrogram and the Spack displays.

Thus, if considering the graphical displays solely, only the Spack displays significantly improves upon the Waveform display. As can be seen on Figure 8, The spectrogram display have the largest dispersion of correct answer rate, i.e. the ratio of correct responses over the number of possible responses, termed p(c) in the following. Considering the distribution of p(c) for the spectrogram display shown on Figure 9, two modes can be observed contrary to the one of the spack display. Even though each sub-jects have been given the same introduction to each of the graphical displays, their familiarity with the standard displays may vary since some subjects had previous training in signal processing courses. This may explain the higher mode in the distribution of the spectrogram display. Even if this observation shall be considered with care due to the rather low number of subjects, this can lead us to conjecture about the influence of the familiarity of the subjects with the spectrogram display on the reported performance. The spack display does not exhibit the same distribution profile and prior familiarity cannot be assumed as the display was equally new to all subjects.

CONCLUSIONS

In this paper, we proposed a display based on the stacking of the envelopes of logarithmically spaced band pass filters. We have shown qualitatively that this kind of display may have some potential as it conveys nicely the distribution of the energy across time and frequency in a way that is an alternative to the one taken when considering the spectrogram.

When considering two evaluation tasks: 1) identifying the type of instrument played, and 2) identifying at which skill level a saxophone tone is played, the spack display compares favorably to more conventional displays, such as the waveform and spectrogram displays. Subjects reported ease of understanding and quick access to important aspects of the sounds.

Future work will focus on the design of validation tasks for the spack display using a wider range of audio data, namely speech and environmental data.

As the spack display is both compact and intuitive, it can be considered as an inspection tool while practicing a musical instrument in order to monitor the control of the nuance and the timbre while playing. Evaluation of the spack display in such a training use case would thus be of interest.

Figure 1 :

 1 Figure 1: Annotation of a church bell from Schafer [11].

Figure 2 :

 2 Figure 2: Standard displays of the sound of a church bell.

Figure 3 :

 3 Figure 3: Processing chain of the spack display.

Figure 4 :

 4 Figure 4: Spectral stack display (spack) of the sound of a church bell. The color code conveys nicely the modulation within each frequency band and the overall disappearance of the high frequency range.

Figure 5 :

 5 Figure 5: Classification performance of the different displays on Task 1 (identifying the musical instrument): sound (S) waveform (W), spectrogram (Spe) and spack (Spa). The star shows the average performance and the length of the vertical line is twice the standard deviation.

Figure 6 :

 6 Figure 6: Graphical displays of forte B tone. Several performance issues can be observed: lack of airflow control at the attack, change of pitch and loudness at 3 seconds and lack of steady airflow during the whole performance.

Figure 7 :

 7 Figure 7: Graphical displays of another forte B tone. Several performance issues can be observed, for example: lack of sharpness at the attack, change of timbre and loudness at 5 seconds.

Figure 8 :

 8 Figure 8: Boxplot display of the differentiation performance of the different displays on Task 2 (detecting the level of the saxophone player): sound (S), waveform (W), spectrogram (Spe) and spack (Spa).

Figure 9 :

 9 Figure 9: Histogram of the classification performance for the spectrogram (Spe) and the spack (Spa) displays. Only the spectrogram display exhibit two modes, suggesting different levels of expertise of the subjects.

Table 1 :

 1 Results of the repeated measure ANOVA evaluating the effect of the type of display on the performance.

		sum sq. df mean sq.	F	p-value
	Type	0.13	3	0.045	5.3	0.003
	Error	0.38	45	0.008
	6.2 Results

A reference implementation as well as all the data discussed in this paper is available at https://bitbucket.org/mlagrange/ paperaudiostackgraph

The sounds and the visual displays are available on the companion website

The sounds and the visual displays are available on the companion website

ACKNOWLEDGMENTS

The authors would like to acknowledge support for this project from ANR project Houle (grant ANR-11-JS03-005-01) and ANR project Cense (grant ANR-16-CE22-0012).


Mira Rizkallah

Francesca De Simone

Thomas Maugey

Christine Guillemot

Pascal Frossard

Rate Distortion Optimized Graph Partitioning for Omnidirectional Image Coding

Omnidirectional images are spherical signals captured by cameras with 360-degree field of view. In order to be compressed using existing encoders, these signals are mapped to planar domain. A commonly used planar representation is the equirectangular one, which corresponds to a non uniform sampling pattern on the spherical surface. This particularity is not explored in traditional image compression schemes, which treat the input signal as a classical perspective image. In this work, we build a graph-based coder adapted to the spherical surface. We build a graph directly on the sphere. Then, to have computationally feasible graph transforms, we propose a ratedistortion optimized graph partitioning algorithm to achieve an effective trade-off between the distortion of the reconstructed signals, the smoothness of the signal on each subgraph, and the cost of coding the graph partitioning description. Experimental results demonstrate that our method outperforms JPEG coding of planar equirectangular images.

I. INTRODUCTION

Nowadays, omnidirectional images are widely used for popular applications such as virtual reality and immersive communications. Omnidirectional images are spherical signals captured by cameras with 360-degree field of view. In order to use existing image and video processing algorithms, these signals are usually mapped to planar domain and stored as rectangular lattices. A commonly used planar representation for omnidirectional content is the so-called equirectangular representation [START_REF] Chen | Recent advances in omnidirectional video coding for virtual reality: Projection and evaluation[END_REF] (Figure 1). This representation presents strong warping distortions around the polar areas and corresponds to an equi-angular sample distribution on the spherical surface, which is non-uniform (Figure 2). An equirectangular image can be fed as input to existing state of the art encoders, but the equirectangular signal statistics differ from those of classical perspective images. Thus, using existing compression algorithms is sub-optimal [START_REF] Simone | Geometry-driven quantization for omnidirectional image coding[END_REF]. In this paper, we propose a graph-based representation for omnidirectional images, which takes into account the spherical geometry and provides a flexible way to efficiently store and compress the visual data. Specifically, we propose to represent an omnidirectional image by a graph, where the graph vertices correspond to the image pixels defined on the spherical surface. The edge weights capture the sampling grid on which the signal is defined. Such a flexible representation permits to go beyond traditional transform coding by moving from classical fixed transforms such as the Discrete Cosine Transform (DCT) to graph-based transforms that are adapted to the actual signal support, such as the Graph Fourier Transform (GFT) [START_REF] Fracastoro | Graph transform learning for image compression[END_REF] [START_REF] Shuman | The emerging field of signal processing on graphs: Extending high-dimensional data analysis to networks and other irregular domains[END_REF]. Due to the high spatial resolution typical of omnidirectional images [START_REF] Chen | Recent advances in omnidirectional video coding for virtual reality: Projection and evaluation[END_REF] , the graph that we propose to build would have a huge number of vertices (>500K). Consequently, the GFT computation in an actual coding pipeline would be unfeasible. Indeed, one would think about using sampling techniques on the sphere. However such coding schemes involve interpolation which makes the distortion control very complex. Instead, we propose an efficient graph partitioning strategy, which takes into account the geometrical information in order to optimize the smoothness of the signals on the subgraphs while keeping a small overhead to code the description of the partition. Finally, we propose a complete GFT-based lossy compression scheme using this partitioning and compare its performance to the classical DCT-based JPEG coding [START_REF]Digital compression and coding of continuoustone still images[END_REF].

Experimental results show that the partitioning provides an effective tradeoff between the smoothness of signals on the subgraphs and the cost of coding the partition. Moreover, the proposed coding scheme outperforms JPEG coding of planar equirectangular images, in terms of Rate-Distortion (RD) analysis using multiple quality metrics.

II. 360-DEGREE IMAGE AS SIGNAL ON A GRAPH

A 360-degree image I can be represented by a signal x ∈ R N defined on an undirected, 4-connected, weighted global graph G = {V, E, W} (See Figure 2), which consists of a finite set of vertices V defined on the image surface, with |V| = N , a set of edges E, and a weighted adjacency matrix W. For i = 1, ..., N , the signal value x i corresponds to the pixel color value at vertex i ∈ V. If there is an edge e = (i, j) connecting vertices i and j, the entry W i,j represents the weight of the edge, otherwise, W i,j = 0. We define the weight of an edge connecting adjacent vertices i and j via a Gaussian kernel weighting function:

W i,j = exp - d geo (i, j) 2 2θ 2 geo (1)
for some parameters θ geo , where d geo (i, j) represents the geodesic distance between vertices i and j capturing the sampling grid on which the vertices are defined (Figure 2). One known operator defined on the graph is the unnormalized graph Laplacian defined as L = D -W where D is a diagonal matrix whose i-th diagonal element D i,i is equal to the sum of the weights of all the edges incident to vertex i. L is a real symmetric matrix, thus, it has a complete set of orthonormal eigenvectors u l = {u l } associated to real nonnegative eigenvalues λ l = {λ l }, with l = 0, 1, ..., N -1, satisfying:

Lu l = λ l u l . (2)
We assume that the graph Laplacian eigenvalues are ordered, i.e., 0 = λ 0 ≤ λ 1 ... ≤ λ N -1 . The graph spectrum is defined as σ(L) = {λ 0 , λ 1 , ...λ N -1 }. Once we have a graph and a signal defined on its vertices, the eigenvectors of the Laplacian are used to define the graph Fourier transform (GFT) [START_REF] Shuman | The emerging field of signal processing on graphs: Extending high-dimensional data analysis to networks and other irregular domains[END_REF] of the signal x. Particularly, the GFT x of any function x : V → R defined on the vertices of G is the expansion of x in terms of the eigenvectors of the graph Laplacian:

x(λ l ) = x, u l = N i=1 x(i)u * l (i). (3)
The inverse graph Fourier transform is then given by:

x(i) = N -1 l=0 x(λ l )u l (i). (4)
The graph Laplacian eigenvalues and eigenvectors carry a notion of frequency. The eigenvector u 0 associated with the eigenvalue λ 0 is constant and equal to 1/

√

N at each vertex. The eigenvectors associated with small eigenvalues vary slowly across the graph, i.e., if two vertices are connected by an edge with a large weight, the values of the eigenvector at those locations are likely to be similar. The eigenvectors associated with larger eigenvalues oscillate more rapidly and are more likely to have dissimilar values on vertices connected by an edge with high weight.

A signal x is considered to be smooth on G if strongly connected vertices have similar signal values [START_REF] Zhou | A regularization framework for learning from graph data[END_REF]. This is usually quantified in terms of the laplacian quadratic form:

S(x) = x T Lx. (5)
In general, graph-based image compression methods use a graph representation as defined above, and perform a GFT to capture the main characteristics of the signal. The coefficients are then encoded instead of original values. The smoother the signal on a graph (smaller S(x)), the more its energy is concentrated in the low frequency GFT coefficients and the more it is easily compressible.

III. PROBLEM FORMULATION

In our case, the common computational limitation of the global graph representation is the maximum acceptable number of vertices in the graph for GFT computation, which limits the resolution of the visual signal that can be supported. In order to cope with the feasibility of the graph-based transform of the signal in high resolution omnidirectional images, the global graph G = {V, E, W} should be separated into several connected components, e.g. M subgraphs {G 1 , ..., G i , ..., G M } by pruning some unreliable edges. The i-th subgraph is G i = {V i , E i , W i } where V i are the vertices in the subgraph, with |V i | = N i < N , E i are their edges, and W i is the weights matrix. x i is the signal defined on the i-th subgraph. The signals on each of the subgraphs are then independently processed, and transformed separately using their respective local Laplacian L i . If the topology and weights of the global graph are fixed, in order to obtain a good compression performance, the graph partition should be chosen such that it leads to smooth representations of the signals inside different subgraphs. On the other hand, it should also be easy to encode, since it has to be transmitted to the decoder for signal reconstruction. Our problem is therefore how to split the fixed global graph into connected components, so that we achieve optimal RD performance and such that all the connected components contain less than N max nodes.. We first pose the problem as a rate-distortion optimization problem defined as:

min G={Gi} D( G) + γR C ( G) + βR B ( G) subject to N i < N max , ∀i (6)
G = {G i } is the global graph based on the geometry defined a priory using Equation [START_REF] Chen | Recent advances in omnidirectional video coding for virtual reality: Projection and evaluation[END_REF] , where some edges are removed. D( G) is the distortion between the original image and the reconstructed one, R C ( G) is the rate cost of the transform coefficients, and R B ( G) is the rate cost of the boundaries for the graph partitioning description. Each of these terms possibly depend on the chosen partition of the graph and of the coding scheme envisioned. We detail each one of them in the next section.

IV. R-D OPTIMIZED GRAPH PARTITIONING A. Distortion Estimation

Since the GFT is orthonormal and independant in each subgraph, the distortion term in the above problem D( G) is equal to the sum of distortions on all subgraphs:

D( G) = M i=1 D(G i ) = M i=1 x i -xi 2 = M i=1 xi -xiq 2 ,
where x i and xi are the original signal and decoded signal in the i th subgraph respectively. xi and xiq are the original and quantized signal GFT coefficients in the i th subgraph.

If we consider a uniform scalar quantizer with small quantization step q for all N coefficients, D( G) can be approximated by:

D( G) = q 2 N 12 (7)
and is thus independent from G. Therefore, the optimization problem ( 6) is reduced to minimizing the rate terms.

B. Rate approximation of transform coefficients

We can evaluate the rate of the GFT coefficients R C (G i ) in a subgraph i using the approximation in [START_REF] Hu | Multiresolution graph fourier transform for compression of piecewise smooth images[END_REF]:

R C (G i ) = S(x i ) = x T i L i x i = l λ l x2 i,l . (8)
The parameters λ l and x i,l are the eigenvalues of the local Laplacian, and the corresponding GFT coefficients of the signal x respectively. Hence, It is an eigenvalue-weighted sum of squared transform coefficients which depends on the underlying local graph L i . Such proxy assumes that the bitrate of the transform coefficients increases when the smoothness of a signal on the graph decreases. While the bitrate needed to code the DC component is not captured by this approximation, we assume that it is only dependent on the number of subgraphs, which can be tuned using the N max constraint in our optimization problem. The higher the N max , the lower the bitrate needed to code the DC coefficients.

C. Rate approximation of the subgraphs boundaries

In fact, in our problem we impose that the pixels of the same subgraph form a connected component. Thus, a common way to code the subgraph membership is to code the boundaries. In order to approximate the coding rate of a boundary B ij between two adjacent subgraphs G i and G j , we use the 4-directional differential freeman chaincodes (DCC) [START_REF] Freeman | On the encoding of arbitrary geometric configurations[END_REF] and estimate the coding rate of the boundary as its entropy computed as follows:

C B (ij) = -# l 4 k=1 p k log 2 p k , (9)
where # l is the number of chaincodes of the boundary and p k , k = 1 : 4 are the probabilities of each of the 4 directions.

D. Minimization of the total coding rate

Using (7) ( 8) and ( 9), the optimization problem in (6) becomes:

min G={Gi} M i=1 x T i L i x i + α 1 2 M i=1 j∈Ni C B (ij), s.t. N (G i ) < N max , ∀i (10)
where N i is the neighborhood of the subgraph G i . The second term is divided by 2 since we only have to code the boundary between any two neighboring regions once. Finding this optimal partition is in general a combinatorial task, so we solve it using traditional agglomerative approximation. To initialize our optimization process, we use the Normalized Cut [START_REF] Shi | Normalized cuts and image segmentation[END_REF] which is well known for favoring the highest smoothness inside partitions. For that, we build a new graph G NC with the same connectivity as G but with weights taking into account both the geodesic distance on the sphere and the euclidean distance in the Y space as:

w(i, j) = exp - d geo (i, j) 2 2θ 2 geo exp - d Y (i, j) 2 2θ 2 Y . (11)
Note that this graph G NC is only used in the Normalized Cut algorithm, and do not serve as a support for GFT thus will not be transmitted. To limit the computation time, the segmentation is performed with recursive 2-way cut algorithm: at each iteration, only the first 2 eigenvectors are computed exploiting the sparse nature of the laplacians. At the output of the initialization, we have an over-segmentation with nonoverlapping subgraphs R = {G 1 , ..., G i , ..., G K }. To model their spatial locality, we construct a subgraph neighborhood matrix E where E(i, j) = 1 indicates that the subgraphs G i and G j are adjacent in the image. In fact, merging any two adjacent subgraphs G i and G j implies re-considering the connections between adjacent pixels on the boundary between them (from the global graph), hence removing the boundary itself. At each iteration of the merging process, we find the two adjacent subgraphs G * i and G * j , which if merged, bring the most significant decrease of the criterion in [START_REF] Shi | Normalized cuts and image segmentation[END_REF] while not exceeding N max nodes in the merged region. In other words,

{G * i , G * j } = max Gi,Gj ∈R ∆R C (G i , G j ) + λ∆R B (G i , G j ), s.t. E(i, j) = 1, N (G i ∪ G j ) < N max (12)
where

∆R C (G i , G j ) = x T i L i x i + x T j L j x j - x T i x T j L i + D ij W ij W ji L j + D ji x i x j ∆R B (G i , G j ) = C B (ij).
∆R C (G i , G j ) and ∆R B (G i , G j ) essentially capture the difference in the rate needed to code the coefficients and the rate to code the boundaries between the two regions before and after merging, respectively.

If ∆R C (G * i , G * j ) + λ∆R G (G * j , G * j ) > 0, we merge G *
i and G * j into one subgraph, and repeat the process until the total rate cannot be further reduced. In the previous formulation and subgraph merging process, we assume that all subgraphs are having the same contribution to the global rate of the whole scheme. However, in our omnidirectional image application, we are interested in giving ideally more rate to the most useful part of the signal, allowing more rate to the subgraphs occupying the biggest surface in the sphere, favoring the merging on smaller surfaces in the spherical domain. Hence, we modify the initial RD gain of Equation ( 12) adding a normalizing factor equal to the area occupied by the merged region on the sphere:

{G * i , G * j } = max Gi,Gj ∈R ∆R C (G i , G j ) + λ∆R B (G i , G j ) A ij , ( 13
)
where A ij is the area on the sphere that the merged region occupies. Such normalization gives more priority to merging in the poles than in the equator. The final algorithm of the partitioning is detailed in Algorithm 1.

Algorithm 1: Rate-distortion optimized Graph partitioning for omnidirectional image coding Data: NcutLabels, Maximum number of Nodes in a Subgraph:

G N cut , N max Result: Labels after merging G f inal Initialization: G = {G i } = G N cut ; Construct the region neighborhood matrix E; Compute ∆R C (G i , G j ) + λ∆R B (G i , G j ) for all i, j
where

E(i, j) = 1 and N (G i ∪ G j ) < N max ; Repeat Find {G * i , G * j }, Eq. (13) G = G \ {G * i , G * j } ∪ {G * i ∪ G * j }; Update E based on the newly merged region Until max (∆R C (G i , G j ) + λ∆R B (G i , G j )) < 0; G f inal = G; V. EXPERIMENTAL VALIDATION
We now move to describe how we use the above graph partitioning algorithm in an omnidirectional image compression scheme. As pointed out in the previous sections, once we solve our optimization problem, we have two kind of information to transmit to the decoder side: the GFT coefficients of the signals in all subgraphs and the description of the partitioning. The transform coefficients are quantized using a uniform quantizer with a fixed step size q for all the bands, then coded with a simple entropy coder. In order to code the partition map, we use the arithmetic edge coder (AEC) proposed in [START_REF] Daribo | Arithmetic edge coding for arbitrarily shaped sub-block motion prediction in depth video compression[END_REF]. The contours are first represented by differential chaincode (DCC) [START_REF] Freeman | On the encoding of arbitrary geometric configurations[END_REF] and divided into segments. Then, to efficiently encode a sequence of symbols in a segment, AEC uses a linear regression model to estimate probabilities, which will be subsequently used in the arithmetic coder.

A. Validation of our rate proxy

Although we do not explore the true rate needed to code boundaries using AEC in our graph partitioning, we can show the accuracy of our proxy. During the optimization process, we compare the rate needed to code a boundary using AEC to our rate proxy using the entropy. Results are shown in Figure 3 in which the x-axis corresponds to the rate needed using AEC, and the y-axis corresponds to the our proposed rate proxy of R B . Although our rate proxy of the subgraphs' boundaries has a very small computation time with respect to AEC, the positive linear trend observed in the plot shows that it is a good approximation.

B. Coding results

We test our method on four grayscale omnidirectional images, namely Metro, Pool, Farm and Hotel shown in Figure 1. Each omnidirectional image is of size (512 × 1024). We test two versions of our scheme that we call WithoutGeometry and WithGeometry. In the first version, the geodesic distance is not taken into account which comes down to set θ geo = ∞ in the construction of graphs G and G NC . Moreover, the merging is done as explained in Equation [START_REF] Yu | Content adaptive representations of omnidirectional videos for cinematic virtual reality[END_REF]. On the other hand, the second version corresponds to our detailed scheme of the previous section taking into account the geometrical information in all stages: normalized cut, merging with equation ( 13) and transform coding. To evaluate the compression performance, we compute the PSNR in two different domains: the equirectangular domain between original and decoded omnidirectional images, and in the spherical domain after performing a uniform sampling of the spherical surface [START_REF] Yu | Content adaptive representations of omnidirectional videos for cinematic virtual reality[END_REF]. Figure 4 shows an example of subgraphs obtained using WithGeometry scheme after fixing λ to 600 for the Pool image, in the equirectangular and the spherical domains. It is clear that subgraphs are adhering to the objects borders in both domains and larger subgraphs are formed around the poles.

Results in Figure 5 show that WithGeometry leads to a better rate-distortion performance in both domains, compared to WithoutGeometry. There is two major explanations of this behavior. First, the global graph in the first case is more adapted to the omnidirectional signal: more specifically, the signal values in horizontally adjacent pixels around the poles are assumed to be more correlated than those which are horizontally adjacent in the equator. This is the case of most of omnidirectional images where poles usually consist of the floor or the sky. In practice, some images like Metro do not totally follow this assumption which explains the comparable performance observed for the two schemes. A second explanation is that in the WithGeometry case, the total rate is allocated more carefully taking into account the area occupied on the sphere. Furthermore, the obtained results show that our proposed schemes outperform classical DCT transform coding scheme in JPEG especially in the low bitrate range, although they can be further improved by optimizing the coding step namely in the quantization and arithmetic coding parts.

VI. CONCLUSION

In this paper, we have proposed a new graph-based framework for omnidirectional image compression. We introduced a new R-D optimized graph partitioning to cope with the feasibility of graph Fourier Transform on global graphs defined on high resolution images. The partition obtained provides an effective tradeoff between the smoothness of signals inside subgraphs and the cost of coding the partition description. Also, we showed that our methods outperform traditional DCT coding schemes at low bitrates. As future work, we investigate the use of different forms of laplacians and focus on the adpatation of the quantization and other coding tools, which may lead to further improvements to the coding performance. Comparison with traditional coding of planar representations other than the equirectangular one [START_REF] Chen | Recent advances in omnidirectional video coding for virtual reality: Projection and evaluation[END_REF], as well as analysis of RD performance on higher resolution test material, will also be performed.

Fig. 1 :

 1 Fig.1: Omnidirectional images (in equirectangular format) used in our experiments: two outdoor images (a and b) and two indoor images (c and d)[START_REF] Xiao | Recognizing scene viewpoint using panoramic place representation[END_REF]

Fig. 2 :

 2 Fig. 2: Equi-angular (i.e., non uniform) sampling on the sphere corresponding to the planar equirectangular representation. A graph is drawn on the spherical surface with vertices corresponding to pixels and edges connecting each pixel with its four closest neighbors (three closest neighbors at the poles).

Fig. 3 :Fig. 4 :

 34 Fig.3: Accuracy of the our rate proxy. x-axis: rate needed to code a boundary using EAC. y-axis: rate computed using our proxy

  Rate-distortion comparison. 1 s row: performance in the equirectangular domain. 2 nd row : performance in the Spherical domain

	PSNR(dB) in the equi-rect domain	25 30 35	0.1	0.2	0.3 Farm	0.4 JPEG WithGeometry 0.5 WithoutGeometry	PSNR(dB) in the equi-rect domain	26 28 30 32 34	0.1	0.2	0.3 Pool	0.4 JPEG WithGeometry WithoutGeometry 0.5	PSNR(dB) in the equi-rect domain	24 26 28 30 32 34	0.1	0.2	0.3 Hotel	0.4 WithGeometry JPEG WithoutGeometry	PSNR(dB) in the equi-rect domain	24 26 28 30 32	0.1	0.2 Metro	0.3 JPEG WithoutGeometry 0.4 WithGeometry
				bitrate(bpp)					bitrate(bpp)					bitrate(bpp)				bitrate(bpp)
	PSNR (dB) in the spherical domain	24 26 28 30 32 34	0.1	0.2	0.3 Farm	0.4 JPEG WithGeometry 0.5 WithoutGeometry	PSNR (dB) in the spherical domain	24 26 28 30 32 34	0.1	0.2	0.3 Pool	0.4 JPEG WithGeometry 0.5 WithoutGeometry	PSNR(dB) in spherical domain	24 26 28 30 32	0.1	0.2	0.3 Hotel	0.4 JPEG WithGeometry WithoutGeometry 0.5	PSNR(dB) in the spherical domain	24 26 28 30	0.1	0.2 Metro	0.3 JPEG WithoutGeometry 0.4 WithGeometry
				Bitrate(bpp)					Bitrate (bpp)				Bitrate (bpp)				Bitrate (bpp)
			Fig. 5: