Adrien Maglo

Clément Courbet

Pierre Alliez

Céline Hudelot

Progressive Compression of Manifold Polygon Meshes

This paper presents a new algorithm for the progressive compression of manifold polygon meshes. The input surface is decimated by several traversals that generate successive levels of detail through a specific patch decimation operator which combines vertex removal and local remeshing. The mesh connectivity is encoded by two lists of Boolean error predictions based on the mesh geometry: one for the inserted edges and the other for the faces with a removed center vertex. The mesh geometry is encoded with a barycentric error prediction of the removed vertex coordinates and a local curvature prediction. We also include two methods that improve the rate-distortion performance: a wavelet formulation with a lifting scheme and an adaptive quantization technique. Experimental results demonstrate the effectiveness of our approach in terms of compression rates and rate-distortion performance.

Introduction

Surface meshes are of common use in a range of application domains such as computer-aided design, simulation, medical imaging, digital heritage and entertainment. The increasing needs for high precision models lead to the generation of complex meshes which must be stored and transmitted over heterogeneous networks. As data storage has a cost and network bandwidths do not grow as fast as the size of these data, solutions must be found to reduce the size of these models through mesh compression.

Among the mesh compression algorithms, the progressive ones allow during decompression to first obtain a coarse version of the mesh. This first level of details (LOD) is then progressively refined as more data is decompressed, until the input mesh is restored. The general goal is to achieve the best rate-distortion (R-D) performance in the sense that each LOD decoded must be as close as possible to the original mesh with the minimum amount of transmitted data.

Previous work on progressive mesh compression has focused on the compression of triangle surface meshes. However, a significant number of carefully designed meshes are composed of polygon faces. In addition, many recent work focused on quad mesh processing. Moreover, for applications such as remote scientific visualization, meshes can contain not only triangular faces and the decompression must restore the initial connectivity. While some approaches have been proposed for single-rate compression [START_REF] King | Connectivity compression for irregular quadrilateral meshes[END_REF][START_REF] Isenburg | Face fixer: compressing polygon meshes with properties[END_REF][START_REF] Kronrod | Efficient coding of non-triangular mesh connectivity[END_REF][START_REF] Isenburg | Compressing polygon mesh connectivity with degree duality prediction[END_REF][START_REF] Khodakovsky | Nearoptimal connectivity encoding of 2-manifold polygon meshes[END_REF] or randomaccessible compression [START_REF] Courbet | Random accessible hierarchical mesh compression for interactive visualization[END_REF], to our knowledge no approaches were proposed for their progressive compression. Our main contributions are summarized as follows:

1. We first propose a simple progressive mesh compression algorithm that compresses any 2-manifold mesh with arbitrary face degrees. 2. From this base algorithm, we describe a curvature prediction method and a connectivity prediction scheme to further reduce the size of geometry and connectivity. 3. We also include two complementary methods that improve the R-D performance. The first method consists in a wavelet formulation of the geometry compression. It contains a lifting step that slightly improves the R-D performance without increasing the final compression ratio. The second method is the adaptive quantization algorithm from [START_REF] Lee | Rate-distortion optimization for progressive compression of 3d mesh with color attributes[END_REF].

It further improves the R-D performance but increases the final compression rate.

Previous work on progressive mesh compression

Connectivity-based

Hoppe introduced the concept of progressive meshes (PM) [START_REF] Hoppe | Progressive meshes[END_REF]. The idea is to incrementally decimate a mesh using the edge collapse operator. The compressed representation consists of the base mesh followed by all parameters required for the incremental reverse operations, called vertex splits. The main advantage of this scheme is its high multi-resolution granularity, together with the possibility to perform selective refinement during decoding. Such granularity is achieved at the cost of low compression rates: in the order of 37 bits per vertex (bpv) with 10 bits quantization. Popović and Hoppe in [START_REF] Popović | Progressive simplicial complexes[END_REF] generalized the PM representation to arbitrary simplicial complexes. They introduced the generalized vertex split and its inverse, the vertex unification operations. With this representation, a model requires about 50bpv with 10 bits quantization.

In order to come closer to compression rates of single rates methods, some methods were proposed to encode the vertex split operations in batches. Taubin et al. [START_REF] Taubin | Progressive forest split compression[END_REF] build a progressive mesh compression scheme inspired by the single-rate topological surgery algorithm. Their progressive forest split representation encodes a manifold triangular mesh with a base mesh and a sequence of forest split operations. The forest split operation consists in cutting the mesh through several sets of connected edges, filling the generated holes with triangles and relocating the vertices. Pajarola and Rossignac [START_REF] Pajarola | Compressed progressive meshes[END_REF] improved the compression rates by imposing some restrictions for choosing the candidates to the edge collapse operations: the operations are grouped into batches during the traversal of a spanning tree. They also improve the geometry coding by using a butterfly predictor. Karni et al. [START_REF] Karni | Efficient compression and rendering of multi-resolution meshes[END_REF] designed a progressive compression scheme which enables the fast rendering of all the LODs. The first step of their algorithm is to create an efficient vertex rendering sequence composed of series of incident vertices. The mesh is then decimated by collapsing edges along this sequence. Better compression ratios are achieved with these approaches (about 30bpv [START_REF] Taubin | Progressive forest split compression[END_REF] and 22bpv [START_REF] Pajarola | Compressed progressive meshes[END_REF] for 10 bits quantization). Nevertheless the multiresolution granularity is impacted compared to the PM representation.

Other progressive compression schemes use vertex removals instead of edge collapses. Li and Kuo [START_REF] Li | Progressive coding of 3-d graphic models[END_REF] pioneered a method based on vertex removal followed by a local patch retriangulation. The connectivity is encoded with a local index which specifies the patch neighborhood pattern and a global index which locates this pattern in the whole mesh. The geometry data is encoded with a barycentric error prediction. The authors also pioneered the idea of adapting the vertex quantization along the transmission of the LOD. Cohen-Or et al. [START_REF] Cohen-Or | Progressive compression of arbitrary triangular meshes[END_REF] used the same decimation mechanism combined with patch coloring to encode the face locations with respect to each patch. They achieved compression rates competitive with single rate techniques (about 23bpv with 10 bits quantization). Alliez and Desbrun [START_REF] Alliez | Progressive compression for lossless transmission of triangle meshes[END_REF] proposed what could be seen as a progressive version of the Touma-Gotsman single-rate encoder [START_REF] Touma | Triangle mesh compression[END_REF]. At each iteration, the mesh is decimated by two deterministic patch traversals, the connectivity being encoded through the valence of the removed vertices. The geometry is encoded through the patch barycentric error prediction in a local Frenet frame. The obtained compression rates are about 13bpv with 10 bits quantization.

Valette et al. build a progressive mesh compression algorithm through a wavelet framework [START_REF] Valette | Wavelet-based progressive compression scheme for triangle meshes: Wavemesh[END_REF]. The initial mesh is progressively decimated with a subdivision scheme tailored to irregular meshes. The connectivity data is composed of all face subdivision operations. The geometry is encoded through a wavelet lifting scheme. The compression rates are slightly better than those from [START_REF] Alliez | Progressive compression for lossless transmission of triangle meshes[END_REF] (about 19bpv with 12 bits quantization).

Geometry-based

Observing that the compressed size of geometry is generally higher than the one of connectivity, Gandoin and Devillers [START_REF] Gandoin | Progressive lossless compression of arbitrary simplicial complexes[END_REF] designed a mesh compression algorithm driven by the geometry. In their scheme, the vertex positions are stored in a kD-tree and the vertex occurrences in each cell are entropy coded. The connectivity is encoded using vertex splits. Peng and Kuo also developed a geometry-driven progressive mesh compression algorithm [START_REF] Peng | Geometry-guided progressive lossless 3d mesh coding with octree (ot) decomposition[END_REF] based on an octree data structure. This coder predicts the connectivity of the mesh from the neighbor vertices geometry during the vertex splits thanks to pivot vertices. This algorithm compresses triangle meshes with about 15bpv with 12 bits quantization.

Improving the Rate-Distortion trade-off

Optimizing the R-D trade-off has been the main focus of recent research on progressive mesh compression. Lee et al. [START_REF] Lee | Rate-distortion optimization for progressive compression of 3d mesh with color attributes[END_REF] showed that the rate-distortion trade-off of the Alliez-Desbrun (AD) coder [START_REF] Alliez | Progressive compression for lossless transmission of triangle meshes[END_REF] is improved by using an adaptive quantization method. The idea consists in choosing between performing a decimation operation or a global quantization operation. They obtain both better R-D performances and compression rates (about 1bpv improvement with 12 bits quantization). Ahn et al. [START_REF] Ahn | R-d optimized progressive compression of 3d meshes using prioritized gate selection and curvature prediction[END_REF] proposed another improvement for the AD coder through an optimized mesh traversal to maximize the number of removed vertices per decimation step. A curvature prediction is also used for encoding the geometry. The latter shares the general idea of spectral methods (see. Section 2.4) because a topology-based Karhunen-Loève transform concentrates the distribution of geometry residuals. The residuals are entropy coded with a bit plane coder. Decimation conquests are interleaved with the transmission of bit planes to improve the R-D performance. They significantly improve the compression rates of the AD coder (about 4bpv reduction with 12 bits quantization).

Valette et al. cast the progressive mesh compression problem as a mesh generation problem [START_REF] Valette | Progressive lossless mesh compression via incremental parametric refinement[END_REF]. The algorithm starts from a coarse version of the initial mesh that is progressively refined using a Delaunay mesh generation approach. When all vertices of the original mesh have been decoded, the initial connectivity is restored by flipping edges. This algorithm is shown to compress efficiently (about 15bpv with 12 bits quantization) and provides good rate-distortion performances. The complete connectivity restoration process is however not guaranteed to succeed. The idea of computing the best decimated version of an initial mesh has been recently further investigated [START_REF] Peng | Feature oriented progressive lossless mesh coding[END_REF]. The algorithm starts from the initial mesh vertex set and recursively splits it into several child subsets. Each time a new vertex subset is generated, a representative vertex of this set is computed. In this hierarchy, the number of children of a set is entropy coded. The offsets between a representative and its parent representative are quantized and entropy encoded. The hierarchy is encoded with the connectivity information in a specific order to achieve the best R-D trade-off. This algorithm yields compression rates at about 16bpv with 12 bits quantization.

Laplacian operator-based

Other schemes propose to compress the vertex positions using the frequency domain through the use of the mesh Laplacian operator. Karni and Gotsman [START_REF] Karni | Spectral compression of mesh geometry[END_REF] compress the geometry of a mesh by computing its spectrum with the eigenvector decomposition of its Laplacian matrix. The spectral coefficients, after being quantized and entropy coded, are sufficient to decompress a good approximation of the initial mesh. The compression, albeit lossy, achieves excellent R-D performance with few coefficients. Mamou et al. [START_REF] Mamou | Shape approximation for efficient progressive mesh compression[END_REF] devised an algorithm that computes the Laplacian matrix of a mesh. The mesh is then approximated with a heat equation and a minimal set of control points. The vertex locations are encoded as residuals from the approximation. The connectivity is encoded by the TG single-rate encoder [START_REF] Touma | Triangle mesh compression[END_REF]. This scheme achieves an excellent compression ratio (about 10bpv with 12 bits quantization) despite a high complexity due to the time spent at solving the heat equation.

Wavelet-based

When the restoration of the initial mesh connectivity is not crucial, some authors resort to semi-regular remeshing. Khodakovsky et al. [START_REF] Khodakovsky | Progressive geometry compression[END_REF] propose an algorithm that remeshes the initial mesh, then compresses it using a wavelet transform based on the Loop filter and a zero-tree coder. This scheme was later improved [START_REF] Khodakovsky | Compression of normal meshes[END_REF] through a normal mesh representation. Payan and Antonini [START_REF] Payan | 3d mesh wavelet coding using efficient model-based bit allocation[END_REF][START_REF] Payan | An efficient bit allocation for compressing normal meshes with an error-driven quantization[END_REF] allocate the bits across the wavelet subbands for the standard and normal mesh representations. In general, the wavelet-based algorithms incorporating a remeshing step provide better compression ratio than pure lossless algorithms.

Handling polygon meshes

An indirect approach to deal with polygon meshes consists of first triangulating the polygon mesh before resorting to an existing method restricted to triangle meshes. This approach was already proposed by Taubin et al. [START_REF] Taubin | Progressive forest split compression[END_REF] to extend the progressive forest split algorithm. While being simple at first glance, this approach is not when caring about restoring the initial connectivity of the finer level. It requires encoding an extra information to remove the edges added during triangulation, and conceptually adds more connectivity by adding more edges than necessary. The theoretical result of Tutte's entropy [START_REF] Tutte | A census of planar maps[END_REF] states that the entropy of a planar graph is expressed in bits per edge. This result and our experiments (see Section 6.1) confirm the intuition that adding extra edges increases entropy and hence that the size of the compressed triangulated mesh is in general superior to that of the compressed original mesh. Li et al. state also that their compression scheme [START_REF] Li | Progressive coding of 3-d graphic models[END_REF] can be used on polygonal meshes. However, they provide no details about how the algorithm would be adapted. It can also be noticed that the experimental results presented in [START_REF] Taubin | Progressive forest split compression[END_REF] and [START_REF] Li | Progressive coding of 3-d graphic models[END_REF] do not compare favorably with recent state of the art techniques. In recent work [START_REF] Peng | Geometry-guided progressive lossless 3d mesh coding with octree (ot) decomposition[END_REF], Peng and Kuo discuss the progressive compression of polygon meshes by an octree coder. As their algorithm can compress arbitrary connectivity between vertices, it is possible to modify the face construction algorithm to reconstruct polygon faces. The mesh connectivity is encoded through vertex splits and efficient prediction of pivots vertices. By definition, pivot vertices are connected to the two vertices of the edge generated by a vertex split, but there is no pivot vertex when the adjacent faces of this edge are not triangles. This encoding scheme is therefore optimized to compress triangle meshes and is not best suited to polygon meshes.

Our algorithm compresses progressively and effectively manifold polygonal meshes with arbitrary face degrees while still be competitive with previous approaches specialized to triangle mesh.

Base algorithm

Compression

The proposed algorithm is based on mesh decimation. It is composed of four main procedures that are successively repeated until the initial mesh M n cannot be further simplified.

• The decimation step consists in applying the patch decimation operator to generate the mesh LOD M l-1 from M l . This operator removes vertices and adds new edges to the mesh. M l connectivity is transformed into face and edge Boolean flags. Face flags indicate if a face has a removed center vertex. Edge flags indicate if an edge was inserted by a remeshing operation. The prediction residuals of the positions of the removed vertices are also stored for later encoding.

• The patch encoding step builds the face symbol list S f and the residual list S r from the face flags and residuals. It consists of a deterministic face conquest of the mesh to encode the faces with a removed center vertex and the position of the removed vertices.

• The edge encoding step builds the edge symbol list S e from the edge flags. A deterministic edge conquest encodes which edges have or have not been inserted.

• The entropy coding step compresses the S f , S r and S e symbol lists.

Decimation step

The decimation step tries to simplify as much as possible a mesh LOD M l to generate M l-1 using the patch decimation operator. A patch is a set of faces with a common center vertex. The algorithm first attempts to form a patch on the mesh. The degree of its center vertex v j must greater than two (see Figure 3(a)). If it succeeds, it then creates a polygon f i around v j . To do so, it splits every faces around v j that is not a triangle by adding a new edge between the two vertices of this face that are connected to v j . We call this remeshing operation re-edging (see Figure 3(b)). This operation aims at generating faces with low degree. The inserted edges, are marked. v j is removed (see Figure 3(c)) and the residual

r f i = p v j -b f i , (1)
where p v j is the position of v j and b f i is the barycenter of f i vertices, is stored. The barycenter of the vertices of a face f i is simply defined as:

b f i = 1 |V f i | v k ∈V f i p v k ,
where V f i is the set of the f i vertices. r f i represents the geometry data. f i becomes a face of the mesh. It is marked as having a removed center vertex and can no longer be implied in a patch decimation operation in the current decimation step.

The patch decimation operation by definition always removes one vertex v j from the mesh. If v j has t incident triangles, the variation of the number of faces in the mesh caused by the operation is equal to 1t. This means that, if v j is only incident to non-triangle faces, t = 0 and the number of faces is incremented by one, while the degree of the neighbor faces is decremented by one (see Figure 3). Other patch decimations on neighbor vertices may make these faces progressively disappear. If v j is only incident to triangle faces, then no re-edging is needed and the number of faces in the mesh decreases (see Figure 4). If v j is a degree d vertex, f i generation can be seen as the merging of d triangular faces, as achieved by superface creation in [START_REF] Kalvin | Superfaces: Polygonal mesh simplification with bounded error[END_REF]. But in this previous work the process is driven by the mesh geometry. In our algorithm, the merging is driven by the connectivity because it must ensure for the decompression that all f i vertices were connected to v j . Besides, re-edging does not triangulate all mesh faces since it only happens where patch decimations are proceeded.

Patch decimation operations must preserve the manifold property of the mesh, so that if a patch border vertex shares an edge with more than 2 other patch border vertices, the patch is not decimated. The generated faces are not necessarily planar. This is not a problem since non-planar faces can be good local approximations. Concave faces however may be problematic to render and may lead to further deteriorations during decimation. In our scheme, a face is said concave if its edges, projected on a plane directed by the face normal, form a concave polygon. A triangle face is by definition always convex. The normal of a non-planar face is computed with Newell's method [START_REF] Tampieri | Newell's method for computing the plane equation of a polygon[END_REF]. If the generation of concave faces is not allowed, the first decimation steps generate only convex faces. When it is no longer possible, the next decimation steps further simplify the mesh by allowing the generation of concave faces. These last unconstrained decimations allow reducing the size of the base mesh that is not compressed in our current implementation. During the decompression, the first LOD displayed to the user is the first that contains only convex faces. The algorithm can also skip the decimation of important vertices to minimize the distortion through a volume-based metric as in [START_REF] Alliez | Progressive compression for lossless transmission of triangle meshes[END_REF].

Once the current patch decimation is over, the algorithm attempts to create other patches with unmarked faces of the mesh. The patch decimation order is not constrained. In our current implementation, the algorithm starts from a random seed vertex and progressively conquers the whole mesh by trying to generate new patches with adjacent vertices that do not belong to already marked faces. An example is depicted in Figure 2. A new mesh LOD M l-1 is obtained when no more patch decimation operations can be performed. The result of the decimation step is the simplified level M l-1 with a set of marked inserted edges and a set of marked faces with a removed center vertex (see Figure 5). About 30% of the vertices of a LOD are removed during a decimation step.

The decimation algorithm stops when no more decimation steps can be performed. The lowest LOD M 0 is called the base mesh. Its size depends on the complexity of the mesh but is generally less than 1% of the total file size. As shown by Figure 6, the decimation step preserves the regularity of the infinite regular structure it decimates. In practice however the regularity gets progressively worse during decimation, as the meshes are not specifically designed to preserve regularity during decimation.

Patch encoding step

The patch encoding step produces the list S f of Boolean face symbols s f i and the list S r of the residuals r f i . s f i codes if f i has a removed vertex or not. The algorithm uses a gate FIFO queue to perform a deterministic conquest of the mesh. A gate is an edge between a conquered and an unconquered face. The first gate of the patch encoding step is specified for the whole mesh compression and hence can not be removed. When a face f i is conquered, its symbol s f i is added to S f . If s f i is true, we also add the projection of r f i in the local Frenet frame to S r (see Figure 7). To avoid the post-quantization step mentioned in [START_REF] Alliez | Progressive compression for lossless transmission of triangle meshes[END_REF] and slightly reduce the entropy, we use the bijection proposed in [START_REF] Lee | Adaptive coarse-to-fine quantization for optimizing rate-distortion of progressive mesh compression[END_REF] to transform the coordinates. The gates of f i are then added to the queue in the counterclockwise order starting from the current gate. The next conquered face is the one pointed by the next gate in the queue. The same conquest order is followed by the decoder during the decompression. An example of a patch encoding step is depicted in Figure 8 (a).

Edge encoding step

To build the edge symbols list S e , a deterministic full conquest of the mesh edges is performed using an edge FIFO queue. When an unconquered edge e k is encountered, a symbol s e k is generated to code if e k was inserted. Note that if e k belongs to two faces that do not have a removed center vertex, it is inevitably original. Therefore, no symbol is generated in this case. The neighboring edges of one e k vertex are then added to the queue. The next edge to conquer is extracted from the queue. An example of an edge encoding step is given by Figure 8 (b). The previously described face and edge binary symbol lists S f and S e represent the connectivity data. The simplification process, depending on the LOD, may lead to biased binary distributions of their symbols (see Figure 11). Therefore, an entropy coder with one adaptive context per list is used to encode S f and S e . The geometry data, the S r list, is also entropy coded with two adaptive contexts: one for the tangential components of r f i and one for its normal component. In our current implementation, we used the range encoder from Schindler [START_REF] Schindler | s range coder[END_REF]. Figure 9 depicts the structure of the compressed data generated by our coder.

Entropy encoding

Decompression

The mesh decompression starts by reconstructing M 0 . Then, by applying the reverse operations of the decimation step, the successive LOD M l are progressively restored to finally obtain M n . Thus, for each LOD M l , a first full conquest is performed to decode the faces with a removed center vertex. These vertices are then inserted at their original positions from the encoded geometry data. A second full conquest decodes and removes the edges that were not present in M l .

Improving connectivity and geometry encoding

The scheme described above allows compressing and decompressing any 2-manifold polygon mesh. We describe next improvements to further reduce the size of the connectivity and geometry. They make our scheme competitive in term of compression ratio with triangle specialized approaches. The typical improvement is 0.3bpv for geometry and 0.7bpv for connectivity.

Predicting connectivity from geometry

As depicted in Figure 10, after one decimation conquest, the average area of the faces with a removed center vertex is greater than the average area of the other faces. This observation allows predicting connectivity from geometry. The prediction algorithm works as follows. During the patch encoding and decoding conquests, the average area of the two types of faces are progressively updated when new faces are conquered. The predicted type of a new conquered face is the one which has the closest average area value to the current face area. A binary symbol is generated to indicate if the prediction is verified. This symbol replace s f i and is also later entropy coded. The binary distribution of this new symbol is for most cases more biased than the simple coding distribution.

For some connectivity however, such as the decimation of a regular hexagonal mesh (see Figure 6(c)), the average area of the two types of faces can be equal. In this case, the type of a face can be predicted from the type of its neighbors. When the difference between the face type pourcentages is below a threshold (experimentally set at 10%), our algorithm predicts the type of a face as the inverse of the type that is the most represented among the adjacent already conquered faces. A similar algorithm is used for the prediction of the edge symbols s e k . If the polygons are well-shaped, then the inserted edges are in general longer than the original edges due to the re-edging process.

To predict the group a face or an edge belongs to, we use average values of face areas and edge lengths computed on the part of the mesh already conquered. We make the assumption that the mesh is uniformly sampled. When the current LOD regularity is bad or the last decimation step removed few vertices, the simple coding scheme is more effective. Therefore, at the be-ginning of each connectivity entropy encoding step, a binary symbol is generated to indicate if the connectivity prediction algorithm is used or not.

Curvature prediction for geometry encoding

As described in Section 3.1.1 the geometry data is composed of the removed vertex residuals r f i . To improve the geometry data encoding, we use a curvature prediction method. Instead of directly encoding the residual r f i , we encode:

g f i = r f i -αl f i , (2)
where l f i is the average Laplacian of the set V f i of the vertices of f i . We define the Laplacian of a vertex v j as:

l v j = 1 |V v j | v k ∈V v j p v k -p v j ,
with V v j the set of neighbor vertices of v j . Therefore we have:

l f i = 1 |V f i | v k ∈V f i l v k .
In our experiments, we set α = 0.5. For most meshes, g f i has a more biased distribution than the one of r f i and hence can be more effectively entropy coded.

Improving the Rate-Distortion

We now describe the inclusion to our codec of two methods that improve its rate-distortion performance. The first is based on a wavelet decomposition with a lifting step. It improves the R-D ratios at low rates by about 25% without impacting the final compression rate. The second is the adaptive global quantization method taken from [START_REF] Lee | Rate-distortion optimization for progressive compression of 3d mesh with color attributes[END_REF]. It improves the R-D ratios at low rates by about 35% but increases the final compression rate by 1 to 2 bpv.

Wavelet formulation of the geometry compression

We formulate here the mesh geometry compression as a wavelet decomposition using the lifting scheme [START_REF] Sweldens | The lifting scheme: A custom-design construction of biorthogonal wavelets[END_REF]. The idea of using a wavelet decomposition for the geometry compression of irregular meshes is not new [START_REF] Valette | Wavelet-based progressive compression scheme for triangle meshes: Wavemesh[END_REF]. However the wavelet decomposition we present is specific to our method as we use different mesh decimation operators and geometry encoding schemes.

When a new level of details is generated, two types of geometry data are computed during the wavelet decomposition:

C l-1 = A l .C l , (3)
D l-1 = B l .C l . (4)
C l-1 is the m × 3 global matrix of the coarse coefficients, the m vertex coordinates of M l-1 . D l-1 is the p × 3 global matrix of the detail coefficients, the p local r f i values that are encoded. The analysis filters A l and B l are defined by the decimation operations completed to generate M l-1 . A l is the matrix that extracts the vertex positions of M l-1 from the vertex positions of M l . B l is the matrix that computes all the details coefficients, as locally defined in (1) (see Figure 12 (a)). During the decompression, once the connectivity has been decoded, the coarse coefficients C l can be obtained from the coarse coefficients C l-1 because the vertices of M l-1 are a subset of the vertices of M l . The details coefficients D l-1 are decoded from the compressed data. From (1), we have:

p v j = r f i + b f i .
So, thanks to this local formula, it is possible to recover C l from C l-1 and D l-1 . This process can be written under the form of the following equation:

C l = P l .C l-1 + Q l .D l-1 ,
where P l and Q l are the synthesis filters determined by the previously described process.

This scheme corresponds to the lazy wavelet transform. It just separates the low and high frequency terms during the compression and reassemble them during the decompression.

Lifting step

To improve the R-D performance of our coder, after a decimation conquest the position of the mesh remaining vertices (C l-1 ) are moved in function of the positions of the removed vertices (see Figure 12 (b)). Each vertex position p v j is locally modified as follows:

p v j = p v j + γ 1 |F v j | f i ∈F v j r f i ,
where F v j is the set of the neighbor faces of v j . If a face f i does not have a removed vertex, then r f i = 0. In our experiments, we set γ = 0.5 because it provided the best results with our datasets. This local process can be formulated as a lifting step in our global wavelet formulation. So, (3) and (4) become:

C l-1 = A l .C l + γL l .B l .C l , D l-1 = B l .C l ,
where L l is the matrix that computes the average residual of the neighbor faces (second term of the formula 5.1.1). During the decompression it is possible to rebuild the matrix L l and then to restore the vertex positions with the decoded residuals using the following formula:

C l = P l .(C l-1 -γL l .D l-1 ) + Q l .D l-1 .

Curvature prediction and residual projection

As explained in Section 4.2, the entropy coder does not directly encode the residuals r f i but instead encodes the g f i values projected in the local Frenet frames. Given [START_REF] Isenburg | Face fixer: compressing polygon meshes with properties[END_REF], the global matrix of the symbols values S l-1 is determined by the following equation:

S l-1 = U l-1 .(D l-1 -αG l-1 .C l-1 ),
where G l-1 is the matrix used to compute the l f i values and U l-1 is the matrix used to perform the local projections. During the decompression, the wavelet coefficients matrix can be restored with

D l-1 = V l-1 .(S l-1 + αG l-1 .C l-1 ),
where V l-1 is the matrix that performs the reverse local projections, before applying the rest of the lifting scheme. The whole process is summarized by Figure 13.

Figure 13: One level wavelet analysis and synthesis lifting scheme.

Adaptive quantization

To get the best rate-distortion performance when compressing a mesh, two variables can be played with: the number of vertices V and the number of geometry quantization bits B. For the single rate compression of triangle mesh, King and Rossignac proposed in [START_REF] King | Optimal bit allocation in compressed 3d models[END_REF] methods to optimize the choice of V and B to minimize either the approximation error or the file size. For the progressive compression of triangle meshes, Lee et al. [START_REF] Lee | Rate-distortion optimization for progressive compression of 3d mesh with color attributes[END_REF] showed that the R-D performance of the original AD coder [START_REF] Alliez | Progressive compression for lossless transmission of triangle meshes[END_REF] can be significantly improved by interleaving decimation conquests with vertex global quantization operations (see Figure 14). The main rationale is that a precise level of the initial quantization is not needed for low LODs which have very few vertices. The decimation contests are encoded with the AD coder and the quantization contests are encoded with the Peng and Kuo geometry coder [START_REF] Peng | Geometry-guided progressive lossless 3d mesh coding with octree (ot) decomposition[END_REF].

The authors propose two methods to choose at each iteration whether to decimate the mesh or to quantize its vertex positions. In the first method, denoted by optimal the decimated mesh and the quantized mesh are generated with their compressed data. The two R-D ratio are then compared. The chosen operation is the one that yields the lowest ratio. This method provides the best R-D performance for all levels of details but is computationally intensive. It requires to generate both meshes, to determine the size of the encoded data and to measure the two distortions with the initial mesh. Therefore, the authors recorded the choices made by their optimal coder on a mesh corpus to learn the parameters µ and β of a function q G (K G ) that provides the best number of quantization bits according to the level of decimation:

q G (K G ) = round(µ * log(K G ) -β)
where round() corresponds to the nearest integer rounding function. K g is defined as:

K G = volume o f bounding box area × number o f vertices . µ = -1.
248 and β = -0.954 are reported as the best parameter values for the selected triangle mesh corpus. By this way, a second method, denoted by quasi-optimal, was proposed to choose at each iteration whether to decimate or to quantize. If the current q G (K G ) value is lower than the current number of quantization bits, then the mesh is quantized. Else, it is decimated. The authors experimentally demonstrated that the quasi-optimal and optimal methods yield similar results. We implemented Lee's quasi-optimal method to improve the R-D performances of our coder at low rates.

Experimental results

The experimental results shown are obtained with an implementation based upon the halfedge data structure of the CGAL library [36]. We observe in our experiments that the computation times are approximately linear in the number of vertices. A 10M face mesh is compressed in 1m48s and decompressed in 1m22s on a desktop computer equipped with an Intel Core i7 CPU clocked at 2.80GHz and 8GB of RAM. We measure the distortion through the METRO software tool [START_REF] Cignoni | Metro: Measuring error on simplified surfaces[END_REF], after triangulating each polygon through inserting a vertex at the barycenter of its vertices and connecting it with all the polygon vertices. All results are obtained with the predictions described in Section 4.

Progressive compression of polygon meshes

This section presents polygon mesh compression results. Figure 15, shows that the lifting scheme presented in Section 5.1.1 clearly improves the rate distortion curve at low rates without increasing the final compression rate. The R-D optimization algorithm provides even better distortion at low rates but increases the overall compression rate. Table 1 lists compression rates on polygon models depicted in Figure 19. In order to compare the effectiveness of our method against simple triangulation prior to compression (see Section 2.6), we triangulate polygonal models by choosing for each polygon an arbitrary vertex as pivot and adding edges between this vertex and all the others. The triangulated models are then compressed with the state of the art progressive compression method specialized to triangle meshes [START_REF] Lee | Rate-distortion optimization for progressive compression of 3d mesh with color attributes[END_REF]. We add to the obtained compression rates the cost of the edge flag encoding required to restore the mesh initial connectivity after decompression (see Section 2.6). This cost is computed for each mesh with the Shannon entropy of the Boolean symbol sequence. The obtained values clearly improve over the trivial one. For our polygon mesh corpus, the trivial method costs on average 4 more bpv than our method. Our codec results are also compared with the results of the single-rate coder from [START_REF] Isenburg | Compressing polygon mesh geometry with parallelogram prediction[END_REF] to evidence the cost of the progressiveness. This cost is high for the simple models (Shark, Teapot, Triceratops, Beethoven, Fandisk, Elephant) as their regularity gets rapidly worse during the decimation. Complex, regular models (Neptune, Chinese lion, Gargoyle, Rabbit) are efficiently compressed. Our algorithm performs well with irregular models (Horse, Fertility, Ramesses, Dinosaur) com-Table 1: Compression rates in bits per vertex, without any R-D optimizations. The first part of the table contains meshes with arbitrary face degrees. The second part contains only triangle meshes. C. stands for connectivity. G. stands for geometry. We first triangulate the polygon meshes before compressing them with [START_REF] Lee | Rate-distortion optimization for progressive compression of 3d mesh with color attributes[END_REF]. We add to the obtained compression rates the cost of the edge flags required to restore the original connectivity.

For our algorithm, we provide a curve with the lifting scheme and a curve with the adaptive quantization algorithm. We notice that several coders achieve better results in term of compression ratio and R-D curves than our algorithm. In particular, [START_REF] Mamou | Shape approximation for efficient progressive mesh compression[END_REF] performs well in terms of compression ratio and R-D distortion, at the price of high compression and decompression times (3m for a mesh with 20,000 vertices). The approach of Valette et al. [START_REF] Valette | Progressive lossless mesh compression via incremental parametric refinement[END_REF] also provides very good results: the compression ratio are high and the R-D curve is excellent with the rabbit model, at the price of not guaranteeing the restoration of the initial connectivity. The octree coder [START_REF] Peng | Geometry-guided progressive lossless 3d mesh coding with octree (ot) decomposition[END_REF] gives good compression rates but the distortion is high at low rates. The R-D optimized coder from [START_REF] Ahn | R-d optimized progressive compression of 3d meshes using prioritized gate selection and curvature prediction[END_REF] also performs well in terms of compression ratio but is weaker than our coder in terms of distortion at low rates. For triangle meshes, our algorithm provides similar distortion at low rates than the approach described in [START_REF] Lee | Rate-distortion optimization for progressive compression of 3d mesh with color attributes[END_REF]. It also yields better compression rates as shown in the second part of Table 1. For the two presented triangle meshes the lifting scheme yields slightly better final compression rates than the adaptive quantization method. For the irregular horse model however, the R-D performance is worse at some point.

While being more general than progressive coders specialized to triangle meshes, our coder achieves competitive results for the compression of triangle meshes. It works with any 2-manifold mesh and exhibits good R-D performances at low rate and average compression rates.

Conclusion

We introduced a new progressive mesh compression algorithm. One distinctive property of our method is that it can handle surface meshes with arbitrary face degree unlike previous approaches which only implement triangle mesh compression. Starting from a simple al-gorithm based on decimation traversals, we propose solutions to improve the compression of both geometry and connectivity. We then incorporate two methods designed to optimize the R-D performance: one based on wavelet lifting scheme, and the other based on adaptive global quantization. Experimental results show the effectiveness of our technical choices. Beside being more general than previous approaches, our method is also competitive for the compression of surface triangle meshes.

As future work, we wish to further improve the compression rates by optimizing the selection of the patch decimation operations to maximize the number of vertices removed per decimation step. We will also investigate the best parameters for the adaptive quantization method specialized to non-triangle meshes, and will handle polygon meshes with boundaries.

Figure 1 :

 1 Figure 1: Levels of details of the quadrangle elephant model generated by our compression algorithm.

Figure 2 :

 2 Figure 2: (a) One example decimation step. 1.A seed vertex (in red) is chosen to form the dark gray patch. 2. The non-triangular faces of the patch are split by inserting edges (in blue) between their two vertices connected to the center vertex. 3, 4, 5. The patch center vertex is removed. The face generated by the vertex removal is marked as having a center vertex removed (in green). Its vertices are marked as visited (in blue). All its unvisited adjacent vertices are added to a FIFO queue. A new patch center vertex (red) is popped out of the FIFO. It may not be possible to build a valid patch around this vertex (green) because, for example, its decimation would not preserve the manifold property. In this case, its unvisited adjacent vertices are added to the FIFO and an other vertex is popped. 6. No more patches can be decimated. The current decimation step is finished. The results are: a set of faces with a center vertex removed (green), a set of faces without a removed center vertex (pink), a set of inserted edges (blue) and a set of original edges (black).

Figure 3 :

 3 Figure 3: Decimation of a patch with non-triangle faces. (a) The active patch is in blue. (b) The faces are split by the re-edging operation with the inserted blue edges in order to create a new polygon around v j . (c) v j is removed. f i is marked as having a center vertex removed.

Figure 4 :

 4 Figure 4: Decimation of a patch with triangle faces. (a) The active patch is in blue. (b) v j is removed. f i is marked as having a center vertex removed.

Figure 5 :

 5 Figure 5: Decimation of an intermediate level of detail of the bunny model. Left: the initial level of detail. Right: the new level of detail after decimation. The inserted edges are depicted in blue. Faces with a removed vertex are depicted in green.

Figure 6 :

 6 Figure 6: Examples of two successive decimations of regular connectivities: a regular grid (a), a regular triangle mesh (b) and a regular hexagonal mesh (c).

Figure 7 :

 7 Figure 7: Encoding of the geometry residual r f i in the local Frenet frame (t 1 , t 2 , n). t 1 takes the direction of the gate edge (in blue). n is the patch normal. And t 2 = -t 1 ∧ n.

Figure 9 :

 9 Figure 9: Structure of the compressed data generated by our coder.

Figure 8 :

 8 Figure 8: (a) Example of a patch encoding conquest. Faces with, resp. without, a removed center vertex are depicted in green, resp. pink. Conquered faces are depicted in gray. The arrows represent the gates. (b) Example of an edge encoding conquest. The inserted edges are depicted in blue. The edges in the FIFO queue are depicted in bold. The conquered edges are depicted in gray. During decompression the same conquest orders are followed to determine the faces with a removed center vertex and the inserted edges.

Figure 10 :

 10 Figure 10: Result of a decimation conquest. The average surface of the face with a removed center vertex (green) is larger than the average surface of the other faces (red). The average length of the inserted edges (blue) is also superior to the average length of the others (black).

Figure 11 :

 11 Figure 11: Example of connectivity symbol distributions after one decimation step shown on the right mesh. The first row shows the distribution of S f and S e . The second row shows the distributions of the prediction symbols.

Figure 12 :

 12 Figure 12: Transversal views of a mesh during the compression with the lifting scheme enabled. (a) The mesh before decimation is depicted by gray dotted lines. The mesh after decimation is depicted by black continuous lines. The detail coefficients r f i are the vectors between the barycenter of the face vertices (b f 1 , b f 2 and b f 3 ) and the removed vertices (v 3 , v 4 and v 5 ). (b) The mesh before the lifting step is in gray dotted lines. The mesh after is in black continuous lines. The lifting step moves the position of the remaining vertices v 1 , v 2 and v 3 according to the neighbor face r f i values to improve the R-D distortion performance.

Figure 14 :

 14 Figure 14: Progressive mesh traditional simplification vs. Lee's method for adaptive quantization [7]. Figure inspired from [7].

Figure 16 :

 16 Figure 16: Decompression of the Bimba model (15770 quads) with the lifting scheme. The final compression rate is 13.8bpv with 12 bits quantization.

Figure 18 :

 18 Figure 18: R-D curves for the compression of the Rabbit model with 12 bits quantization.

Figure 19 :

 19 Figure 19: Input meshes fromTable1and their tessellations:(a) Beethoven (b) Bunny (c) Triceratops (d) Teapot (e) Elephant (f) Shark (g) Neptune (h) Gargoyle (i) Lucy VSA (j) Chinese lion (k) Hippo (l) Fandisk (m) Rabbit (n) Dinosaur (o) Venusbody (p) Fertility (q) Horse (r) Ramesses.

Acknowledgements

The authors would like to thank Ho Lee for providing some experimental compression results first published in [7] and Jae-Kyun Ahn for providing the experimental compression results of its encoder [20]. Some models are courtesy of the AIM@SHAPE Repository. This work has been founded by French National Research Agency (ANR) through COSINUS program 12 (project COLLAVIZ n • ANR-08-COSI-003) and the European Research Council (ERC Starting Grant "Robust Geometry Processing", Grant agreement 257474).

pared to the single-rate approach. It even improves on very irregular models such as the VSA-remeshed [START_REF] Cohen-Steiner | Variational shape approximation[END_REF] Lucy or the Hippo models. Figure 16 depicts the decompression of the Bimba quadrangle surface mesh. We now compare our coder with state-of-the-art compression methods specialized to triangle meshes.

Progressive compression of triangle meshes


Paul Escande
email: paul.escande@gmail.com

Pierre Weiss
email: pierre.armand.weiss@gmail.com

Franc ¸ois Malgouyres
email: francois.malgouyres@math.univ-toulouse.fr

SPATIALLY VARYING BLUR RECOVERY Diagonal Approximations in the Wavelet Domain

Keywords: Image deblurring, spatially varying blur, operator approximation, wavelet packet transform, bi-harmonic spline interpolation, convex optimization

. We provide various theoretical and numerical results in order to justify the proposed methods. We also show preliminary deblurring results illustrating the relevance of our formalism.

INTRODUCTION

Image restoration in the presence of spatially varying blur is a problem of increasing importance. It was first studied in the context of satellite imaging with Hubble space telescope [START_REF] Nagy | Restoring images degraded by spatially variant blur[END_REF]. It is now becoming increasingly important with the emergence of new fluorescence microscopes, producing highly deteriorated images, since light interacts with the biological tissues. In microscopy, it is often possible to incorporate micro-beads in the medium surrounding the sample or even in the sample itself, giving access to the point spread function (PSF) of the system at some known locations (see e.g. [START_REF] Preibisch | Software for bead-based registration of selective plane illumination microscopy data[END_REF][START_REF] Temerinac-Ott | Multiview deblurring for 3-d images from light sheet based fluorescence microscopy[END_REF]). This information allows to interpolate the PSF in the whole space and thus to get approximations of the degradation operator for further processing.

In the case of spatially invariant blur, fast deconvolution algorithms can be devised since the convolution is diagonal in the Fourier domain. This allows using O(dn d log(n)) algorithms (where d denotes the space dimension and n d denotes the number of pixels) based on the fast Fourier transform. These approaches are unsuitable in the case of spatially varying blurs and it appeals for the development of new fast numerical algorithms. Our aim in this paper is to propose fast O(n d ) algorithms based on the wavelet or wavelet packet transforms.

We consider a blurring operator H in R d and defined for any u ∈ L 2 (Ω) as the following integral operator : ∀x ∈ Ω, (Hu)(x) := y∈Ω K(x, y)u(y)dy,

where Ω ⊆ R d is the image domain. The function K(x, •) is a spatially varying kernel defining the PSF at each location x. In all the following, we assume that H is a bounded linear operator from L 2 (Ω) to L 2 (Ω).

The most naive approach to compute Hu numerically consists in discretizing (1) by : ∀x ∈ X, Hu(x) = ∑ y∈X K(x, y)u(y),

where X ⊂ Ω denotes the set of pixels locations. This approach is simple to implement, but costs O(n 2d ) arithmetic operations. This is unsuitable for large 2D images or medium sized 3D images. Two alternative approaches are commonly used :

-The first one consists in approximating K(x, •) by a tensor product of kind :

K(x, (y 1 , • • • , y d )) = d ∏ k=1 K k (x, y k ).
FIGURE 1 -An orthogonal view of a Variable Refractive Gibson and Lanni PSF obtained with the PSF Generator plugin for ImageJ [START_REF] Kirshner | 3D PSF models for fluorescence microscopy in ImageJ[END_REF] This reduces the computational cost to O(dn d+1 ) operations which is usually tractable even in large scale scenarii. Moreover this model is exact for Gaussian PSF which sometimes accurately describe perfect microscopy systems [START_REF] Zhang | Gaussian approximations of fluorescence microscope point-spread function models[END_REF]. Unfortunately, it is too rough to describe more complex patterns commonly encountered in optical diffraction or sample induced degradations. Figure 1 shows a typical PSF in three dimensions, which cannot be approximated by separable functions. -The second one consists in using piecewise constant blurs using local FFT [START_REF] Nagy | Restoring images degraded by spatially variant blur[END_REF][START_REF] Hansen | Deblurring images : matrices, spectra, and filtering[END_REF]. Its complexity is roughly the same as that of spatially invariant blurs in O(dn d log(n)). Moreover, this approach also allows to use linear interpolations of the PSF. This might be interesting in our case since the PSF is known only at a few positions and linear interpolations allow an operator reconstruction on the whole image domain. Unfortunately, piecewise constant blurs or linear interpolations are too rough to describe some practical settings. This is illustrated on Figure 9(b), where we can observe that the linear interpolation gives a cross like PSF in the middle, that would be undesirable in practical cases. In this work, we propose to approximate H by operators diagonal in wavelet or wavelet packet transforms. More precisely, we show that H ΨΣΨ * , where Ψ denotes the wavelet transform and Σ is a diagonal matrix. The computation of Hu is thus reduced to O(n d ) operations.

The structure of the paper is as follows. In section 3, we justify the use of such a structure by theoretical and numerical results. In section 4, we propose an algorithm to reconstruct the diagonal operator Σ when the impulse response of H is given at some known locations. Finally, we present a deconvolution algorithm and some results in section 5.

NOTATION

In order to simplify the notation, we consider wavelet transforms and not wavelet packet transforms. We present the theoretical results in 1D for the sake of simplicity and clarity and the experimental results in 2D. The proposed approaches can be extended to any dimensions and would be particularly suited to large 3D problems.

We consider an orthogonal wavelet basis of L 2 (R) :

{φ l 0 ,n } n∈Z ∪ {ψ j,n } j≤l 0 ,n∈Z , where ψ j,n (t) = √ 2 -j ψ(2 -j t -n), (2)
and ψ is the mother wavelet. The function φ l 0 ,n is defined by

φ l 0 ,n (t) = √ 2 -l 0 φ(2 -l 0 t -n),
where φ is the scaling function. In all the paper, Ψ * denotes the forward wavelet transform and Ψ denotes its inverse (in the discrete and in the continuous setting). F * denotes the Fourier transform (discrete or continuous) and F denotes its inverse. The convolution between u and h is denoted h u. The Fourier transform of u is denoted û or F * u

The indicator function of a convex and closed set C ⊆ R n is denoted χ and defined as :

χ C (x) = 0 if x ∈ C +∞ otherwise.
The proximal operator or resolvent of a convex, closed function F : R n → R ∪ {+∞} is defined for all x 0 ∈ R n by :

(I + ∂F * ) -1 (x 0 ) = arg min x∈R n F(x) + 1 2 x -x 0 2 2 .
For a discrete image in R d , we define the discrete partial derivative in direction i by :

∂ i u(•, k, •) = u(•, k + 1, •) -u(•, k, •) if 1 ≤ k < n 0 if k = n.
where the indice k is that the i-th position in the array.

The discrete gradient operator in R d is defined by :

∇ =     ∂ 1 ∂ 2 . . . ∂ d     Let q ∈ (R n d ) d represent a discrete vector field. We set q =     q 1 q 2 . . . q d     .
The isotropic l 1 -norm in (R n d ) d is defined by :

q 1,2 = n d ∑ i=1 d ∑ j=1 q j (i) 2 .
Finally, the discrete total variation of u ∈ R n d is defined by : TV (u) = ∇u 1,2 .

DIAGONALIZATION OF THE VARIABLE BLUR OPERATOR IN A WAVELET BASIS

The main ingredient allowing the design of efficient deconvolution algorithms is the fact that a convolution is diagonalized in the Fourier domain. For any kernel h, h u = F ΣF * u where Σ can be considered as a diagonal operator that multiplies F * u by F * h.

The main idea of this paper is to mimic this property for spatially varying blur operators. We propose to approximate H by an operator H diagonal in the wavelet domain :

Hu Hu := ΨΣΨ * u = ∑ n∈Z φ l 0 ,n , u φ l 0 ,n + ∑ j≤l 0 ,n∈Z σ j,n ψ j,n , u ψ j,n
where (σ j,n ) j,n is a sequence of weights that will be described later. This particular structure allows to approximate Hu in O(n d ) arithmetic operations, which is doable even for very large scale problems. Such operators have been deeply analyzed from a theoretical point of view in various articles or monographs (see e.g. [START_REF] Beylkin | Fast wavelet transform and numerical algorithm[END_REF][START_REF] Coifman | Wavelets, calderónzygmund and multilinear operators[END_REF]). However, we found very few image processing applications in the literature. To our knowledge, the closest practical application is dedicated to the fast computation of image foveation [START_REF] Chang | Wavelet foveation[END_REF]. However, this work is only adapted to very particular kind of kernels K met in foveation that do not correspond to our practical problems.

Since H is a linear operator in a Hilbert space, it can be written as :

H = ΨΘΨ * ,
where Θ : l 2 → l 2 is characterized by the coefficients,

(θ j,m,k,n ) j,m,k,n := ( Hψ j,m , ψ k,n ) j,m,k,n .
In order to justify the proposed approach, we first recall some theoretical results presented in [START_REF] Beylkin | Fast wavelet transform and numerical algorithm[END_REF] that assess the decrease of θ j,m,k,n away from the diagonal (i.e. when |m -n| > 0 and | j -k| > 0). Then we provide an interpretation of the coefficients σ j,n in terms of amplitudes of the Fourier coefficients of the local PSF.

Decay of Θ Away from the Diagonal

In [START_REF] Beylkin | Fast wavelet transform and numerical algorithm[END_REF], it has been proved that, for compactly supported wavelets possessing M vanishing moments and smoothly varying kernels, the values of Θ are small away from the diagonal in the one and two-dimensional cases. Typical results are as follow : Theorem 1 ( [START_REF] Beylkin | Fast wavelet transform and numerical algorithm[END_REF]). Suppose that |K(x, y)| ≤ 1 |x-y| and that K(x, y) is of class C M+1 with , ) , where M denotes the number of vanishing moments of ψ. Then θ j,m,k,n satisfies the following inequality :

|∂ M x K(x, y) + ∂ M y K(x, y)| ≤ C M |x -y| (1+M
|θ j,m,k,n | ≤ O 1 1 + | j -k| M+1 .
Moreover, for compactly supported kernels K :

|θ j,m,k,n | = 0, for sufficiently large |m -n|.
The authors also show that the operator norm H -Ψ ΘΨ * can be made arbitrarily small if Θ is obtained by thresholding Θ in such a way that only O(n d ) coefficients are kept. It roughly means that if K is a smooth kernel, computing Hu can be performed in O(n d ) operations, rather than O(n 2d ), by making use of the wavelet transform. In this work, rather than considering sparse matrices Θ, we use simpler diagonal matrices.

We illustrate these results experimentally in the discrete setting on Figure 3. We consider an operator H whose kernel is a two-dimensional Gaussian with variances linearly increasing in the vertical direction, see Figure 2(c). This operator applied to the mandrill image results in the image Figure 2(b). The matrix Θ is shown on Figure 3. It is seen that Θ is dominated by its diagonal entries and that the coefficients away from the diagonal decrease extremely fast (actually much faster than the result in Theorem 1).

Interpretation of the Diagonal Values

In this paragraph, we show that the values σ j,n can be interpreted as local frequency responses of H. We assume that ψ is a compactly supported wavelet on the interval [-β, β].

Let us analyze the impulse response of H at point x :

Hδ x = ΨΣΨ * δ x = ∑ n∈Z φ l 0 ,n (x)φ l 0 ,n + ∑ j≤l 0 ,n∈Z σ j,n ψ j,n (x)ψ j,n = ∑ n∈Z φ l 0 ,n (x)φ l 0 ,n + ∑ j≤l 0 , n∈k(x, j) σ j,n ψ j,n (x)ψ j,n ,
where k(x, j) := n ∈ Z such that 2 -j xn < β . The sets k(x, j) are represented in Figure 4 in the two-dimensional case. They contain at most 2β elements. Now, if we assume that σ j,n varies little in k(x, j) and satisfies σ j,n σ j,x we obtain :

ΨΣΨ T δ x ∑ n∈k(x,l 0 ) φ l 0 ,n (x)φ l 0 ,n + ∑ j≤l 0 σ j,x ∑ n∈k(x, j) ψ j,n (x)ψ j,n .
The local frequency response of H is thus where α j,x is a complex coefficient that depends on the choice of ψ and ψ j (x) = √ 2 -j ψ(2 -j x). Since ψ j is well localized in the frequency domain, the coefficient σ j,x α j,x can be interpreted as a local frequency attenuation in a certain frequency band that depends solely on the scale j. This principle is illustrated in Figure 5.

ΨΣΨ * δ x ∑ n∈Z φ l 0 ,n (x) φ l 0 ,n + ∑ j≤l 0 σ j,x ∑ n∈k(x, j) ψ j,n (x) ψ j,n = ∑ n∈Z φ l 0 ,n (x) φ l 0 ,n + ∑ j≤l 0 σ j,x α j,x ψ j ,

Spatial Regularity of the Eigenvalues

A simple way to find a matrix Σ such that H H consists in setting Σ = Diag(Θ). If the kernel K varies sufficiently smoothly in space, the discrete values (σ j,n ) n∈Z also vary smoothly, meaning that σ j,n σ j,n+1 . This can be verified experimentally : Figure 6 represents the diagonal of Θ for an operator H displayed in Figure 2(c) in the usual wavelet domain.

The eigenvalues vary smoothly in each sub-band. This remark is central to understand the interpolation algorithm proposed in the next section. Also notice that the coefficients σ j,n decrease from the top to the bottom of the image at each scale. It means that the high-frequencies are attenuated on the image bottom. This clearly corresponds to the operator H shown in Figure 2(c).

Wavelet and Basis Choice

An important parameter in the proposed algorithm is the mother wavelet ψ. According to Theorem 1, ψ should have many vanishing moments for H to approximate H correctly. This condition can be verified experimentally. In the following, we will always consider Daubechies 8 orthogonal wavelet bases which appeared to produce good practical results. It was shown in [START_REF] Malgouyres | A framework for image deblurring using wavelet packet bases[END_REF][START_REF] Kalifa | Deconvolution by thresholdong in mirror wavelet bases[END_REF], that wavelet packet bases are more adapted to deconvolution problems than standard wavelet bases. The reason is that they provide more precise fre-quency tilings, notably in the high frequencies. This is illustrated in Figure 7. The frequency tiling of the standard wavelet decomposition is much coarser than that of a wavelet packet transform.

Finally, the non translation invariance of discrete wavelets tend to produce ringing artifacts. They can be reduced using redundant bases such as the translation invariant wavelets. In this work we use cycle spinning [START_REF] Coifman | Translation-invariant (a) Original Image (b) Degraded Image SNR = 28.14 (c) Deblurred Image[END_REF], which consists in averaging wavelet transforms translated in the set {0, 1} d . This transform is a tight frame with a redundancy of 2 d .

OPERATOR RECONSTRUCTION FROM LOCALLY KNOWN PSF

In this section, we propose a method to recover the matrix Σ from the knowledge of local impulse responses. This setting corresponds to various practical applications. In astronomy, stars may sometimes be considered as Diracs. Their observation thus provides the impulse response of the system K(x, •), where x denotes the star location. In microscopy, micro-beads may be inserted in the sample and provide the impulse responses at locations spread in the whole image domain.

The problem tackled in this section is the reconstruction of K everywhere, from the knowledge of K(x i , •) at a few locations (x i ) i∈{1,••• ,m} . We assume that two images are available :

-An image

u = m ∑ i=1 δ x i ,
that describes the Dirac locations. -An image u o = Hu which provides the impulse responses at locations x i . Figure 8 illustrates two images u and u o . The Diracs could be randomly located on the image rather than on a uniform grid. We considered this simple setting for experimental reasons. The number of known impulse responses can also be considerably reduced as will be shown later. The knowledge of u o allows to reconstruct the eigenvalues σ j,n of H only close to the known locations x i in each sub-band. These eigenvalues should thus be interpolated in order to recover K everywhere. Note that this problem is not standard since is consists in interpolating an operator eigenvalues and not an image.

Since the eigenvalues vary smoothly in space, we propose to use bi-harmonic splines which are well adapted to scattered data interpolation [START_REF] Wahba | Spline models for observational data[END_REF].

The approximation problem we propose formulates as the following variational problem :

Find Σ ∈ arg min ΨΣΨ * u -u o 2 2 + λR (Σ) , (3)
where λ > 0 is a regularization parameter. We also set :

R(Σ) = ∑ j≤l 0 ∆σ j 2 2 ,
where ∆ denotes the discrete Laplacian and σ j denotes the set of eigenvalues at scale j. This energy provides the approximation of minimal curvature. It is equivalent to using bi-harmonic splines [START_REF] Wahba | Spline models for observational data[END_REF].

The quadratic structure of problem (3) allows the use of conjugate gradient like methods for the minimization. We are currently investigating the use of preconditionners in the wavelet domain for accelerating the convergence.

We present approximation results in Figures 9 and 10. Figure 10 displays a interpolated matrix Σ. This result can be evaluated by comparing it with the true diagonal of Θ presented in Figure 6. Overall, the reconstruction leads to near perfect results. Figure 9 compares the interpolation provided by Fourier based methods such as [START_REF] Nagy | Restoring images degraded by spatially variant blur[END_REF] with the proposed approach. Our method produces some artifacts, however, the proposed interpolation is rather close to the reality in the image center. Note that this result is obtained using knowing the PSF at only 4 locations in the plane. Deblurring an image with kernel 9(b) would be disastrous, since horizontal and vertical frequencies would be enhanced, leading to strong ringing artifacts.

DEBLURRING ALGORITHM AND RESULTS

Deblurring Method

In this section, we assume that the diagonal Σ has been reconstructed using the method proposed in section 4. We propose a total variation (TV) based algorithm to tackle the deblurring problem. We suppose that a degraded image v o is obtained according to the following discrete model :

v o = Hv + η,
where α > 0 is a user fixed parameter and TV (v) is the isotropic total variation of v defined in the notation.

In settings where H is perfectly known, users should set α = σ 2 n. The proposed approach slightly differs since total variation serves as a regularizer for both the noise and the errors in the operator approximation. In practice we found that setting α = (1 + ε)σ 2 n where ε > 0 is a small parameter provides good experimental results.

Problem (4) can be rewritten as :

min v∈R n d F(Av) + G(v), (5)
where

A : R n d → R dn d × R n d v → ∇v Hv , F : R dn d × R n d → R + (y 1 , y 2 ) → y 1 1,2 +X {y, y 2 -v o 2 ≤α} (y 2 ) and G = 0.
This reformulation allows to use the primal-dual algorithm proposed in 1 detailed in algorithm 1. k . This rate is somehow optimal in the class of first order methods [START_REF] Nemirovskii | Complexity of problems and efficiency of optimization methods[END_REF]. The proximal operators of F * and G are easy to compute analytically and we refer the reader to [START_REF] Combettes | Signal recovery by proximal forward-backward splitting[END_REF][START_REF] Chambolle | A first-order primaldual algorithm for convex problems with applications to imaging[END_REF] for more details. The steps sizes τ and σ should satisfy στ A 2 < 1 for convergence, where A = max j,n |σ j,n |. In practice, this algorithm requires around 100 low cost iterations to provide satisfactory results for the visual system. Note that H and H must be applied to a vector at each iteration. This is performed in O(n d ) arithmetic operations due to the special structure of H.

Results for Different Noise Variances

Now we will present some results of the deblurring algorithm. We used the Mandrill Figure 11(a) rescaled in [0, 1] and blurred with an operator having a twodimensional Gaussian PSF with variance increasing linearly in the vertical direction. It impulse responses are displayed in 11(b). In Figures 12 and13 we respectively added a noise of variance σ η = 0 and σ η = 3.10 -2 .

In the case σ η = 0, Figure 12 shows that the algorithm is able to recover thin details of the image even in the coat and the beard of the Mandrill in the bottom of the image. This highlights the fact that the approximation of H by H is sufficiently good for the sake of deblurring.

In the case of a larger noise, σ η = 3.10 -2 , Figure 13 shows that the image quality is improved but suffers from the standard defects of total-variation based regularizations : stair-case appears and thin details are not recovered. Overall, these results confirm that the proposed approximation is capable of producing nice reconstruction results with low computational costs.

CONCLUSION

This paper contains various contributions. First, we showed that spatially varying blur operators can be well approximated by operators diagonal in the wavelet domain. Second, we proposed an original approach that allows to reconstruct blur operators when their PSF is known only at a few discrete locations. Compared to previously proposed approaches, this method allows a much wider class of interpolation algorithms. In particular, we showed that bi-harmonic splines interpolations produce near perfect results on a few examples and outperform standard approaches. Finally we presented preliminary deblurring results. These results outline that the proposed approximations provide results similar to what would be obtained with a perfect knowledge of the blur operator.

The proposed algorithm will be validated on real data coming from emergent microscopes such as the selective plane illumination microscope.

  PSF at various locations FIGURE 2 -Image blurred using the operator H. The kernel K of the operator is a Gaussian which grows linearly in the vertical direction.

  (a) In a linear scale (b) In a log 10 scale FIGURE 3 -Matrix Θ for the variable operator illustrated in Figure 2. This matrix is obtained using Daubechies 8 wavelets and a decomposition level J = 2.

FIGURE 4 -

 4 FIGURE 4 -The sets k(x, j) are indicated in orange at each scale.

  FIGURE 5 -Local Fourier attenuation are determined by the coefficients σ j,x α j,x .

  (a) In a wavelet wasis (b) In a wavelet packet basis FIGURE 7 -Frequency tilings for Daubechies-8 wavelets in the wavelet and wavelet packet basis. Each color corresponds to the Fourier transform of a wavelet at a given scale.

  (a) The dirac map u (b) u o = Hu FIGURE 8 -Dirac map and the associated impulse responses. This information is used to reconstruct an approximation of the blurring operator H.

  where H : R n d → R n d is the spatially varying blur operator, v ∈ R n d is the unaltered image and η ∈ R n d is a white Gaussian noise, η ∼ N (0, σ η Id n d ). Our aim is to recover v knowing v 0 . Since H and H are usually compact, the inverse problem of recovering v should be regularized. We propose to use a standard total variation based reconstruction approach. It reads :

  Operator reconstruction using different methods. The operator is reconstructed using the information available in the red rectangles.

FIGURE 10 -

 10 FIGURE 10 -The matrix Σ reconstructed using bi-harmonic splines. It should be compared to the real diagonal presented on Figure 6.

Algorithm 1 :

 1 Primal-Dual Algorithm (Chambolle and Pock, 2011) Input: ε : the desired precision ; (x 0 , y 0 ) : a starting primal-dual pair ; Output: x : an approximate solution to the problem. begin Choose, τ, σ > 0, θ ∈ [0, 1] ; k = 0; v0 = v 0 ; while Convergence Criterion > ε do y k+1 = (I + σ∂F * ) -1 (y k + σA xk ) ; x k+1 = (I + τ∂G) -1 (x k + τA * y k+1 ) ; xk+1 = x k+1 + θ(x k+1x k ) ; k = k + 1; end end This algorithm is easy to implement and it can be shown to converge in O 1

  The original Image and the impulse response of the blurring operator.

  Degraded Image SNR = 22.51 (c) Deblurred Image, SNR = 29.02 FIGURE 12 -Restoration results for σ η = 0.

  FIGURE 13 -Restoration results for σ η = 3.10 -2 .

  The original image and the and the impulse response of the blurring operator.

ACKNOWLEDGEMENTS

This work was partially funded by the Mission pour l'interdisciplinarité from CNRS, Région Midi Pyrénées and the ANR SPHIM3D.

The authors wish to thank J. Bigot and J. Fehrenbach for interesting discussions. They also thank all the ITAV staff, Charlotte Emery and Emmanuel Soubies for their tireless support during this project. (a) Original Image (b) Degraded Image SNR = 30.30 (c) Deblurred Image, SNR = 36.13 FIGURE 15 -Restoration for σ η = 0.