Clément Jamin
email: clement.jamin@liris.cnrs.fr

Pierre-Marie Gandoin
email: pierre-marie.gandoin@liris.cnrs.fr

Samir Akkouche
email: samir.akkouche@liris.cnrs.fr

CHuMI Viewer: Compressive Huge Mesh Interactive Viewer

Keywords: Lossless compression, Interactive visualization, Large meshes, Out-of-core

come   L'archive ouverte pluridisciplinaire

Introduction

Mesh compression and mesh visualization are two fields of computer graphics that are particularly active today, whose constraints and goals are usually incompatible, even contradictory. The reduction of redundancy often goes through signal complexification, because of prediction mechanisms whose efficiency is directly related to the depth of the analysis. This additional logical layer inevitably slows down the data access and creates conflicts with the speed requirements of real-time visualization. Conversely, for efficient navigation through a mesh integrating the user's interactions dynamically, the signal must be carefully prepared and hierarchically structured. This generally introduces a high level of redundancy and sometimes comes with data loss, if the original vertices and polyhedra are approximated by simpler geometrical primitives. In addition to the on-disk storage and network transmission issues, the data growth implied by this kind of preprocessing could become detrimental to the visualization itself since the gap between the processing times and the access times from external memory is increasing.

In this article, we attempt to reconcile compression and interactive visualization by proposing a method that combines good performance in terms of both the compression ratio and the visualization frame rates. The interaction requirements of a viewer and a virtual reality system are quite similar in terms of the computational challenge but the two applications use different criteria to decide how a model is viewed and updated. Our goal here is to take advantage of excellent view independent compression and introduce distance-dependent updates of the model during the visualization. As a starting point, the in-core progressive and lossless compression algorithm introduced by Gandoin and Devillers [START_REF] Gandoin | Progressive lossless compression of arbitrary simplicial complexes[END_REF] has been chosen. On top of its competitive compression ratios, out-of-core and LOD capabilities have been added to handle meshes with no size limitations and allow local refinements on demand by loading the necessary and sufficient data for an accurate real-time rendering of any subset of the mesh. To meet these goals, the basic idea consists in subdividing the original object into a tree of independent meshes. This partitioning is undertaken by introducing a primary hierarchical structure (an nSP-tree) in which the original data structures (a kd-tree coupled to a simplicial complex) are embedded, in a way that optimizes the bit distribution between geometry and connectivity and removes the undesirable block effects of kd-tree approaches.

After a study of related works (Sec. 2) focusing on the method chosen as the starting point (Sec. 3), our contribution is introduced by an example (Sec. 4), then detailed in two complementary sections. First, the algorithms and data structures of the out-of-core construction and compression are introduced (Sec. 5 and 6), then the visualization point of view is adopted to complete the description (Sec. 7). Finally, experimental results are presented and the method is compared to prior art (Sec. 8).

Previous Works

Compression

Mesh compression is a domain situated between computational geometry and standard data compression. It consists in efficiently coupling geometry encoding (the vertex positions) and connectivity encoding (the relations between vertices), and often addresses manifold triangular surface models. We distinguish single-rate algorithms, which require full decoding to visualize the object, and multiresolution methods that allow one to see the model progressively refined while decoding. Although historically geometric compression began with singlerate algorithms, we have chosen not to detail these methods here. The reader can refer to surveys [START_REF] Gotsman | Simplification and compression of 3d meshes[END_REF][START_REF] Alliez | Recent advances in compression of 3d meshes[END_REF] for more information. For comparison purposes, since progressive and outof-core compression methods are very rare, Sec. 8 will refer to Isenburg and Gumhold's well-known single-rate out-of-core algorithm [START_REF] Isenburg | Out-of-core compression for gigantic polygon meshes[END_REF]. Similarly, lossy compression, whose general principle consists in a frequential analysis of the mesh, is excluded from this study.

Progressive compression is based on the notion of refinement. At any time of the decoding process, it is possible to obtain a global approximation of the original model, which can be useful for large meshes or for network transmission. This research field has been very productive for approximately 10 years, and rather than being exhaustive, here the choice has been to adopt a historical point of view. Early techniques of progressive visualization, based on mesh simplification, were not compression-oriented and often induced a significant increase in the file size, due to the additional storing cost of a hierarchical structure [START_REF] Hoppe | Mesh optimization[END_REF][START_REF] Popović | Progressive simplicial complexes[END_REF]. Afterward, several single-resolution methods were extended to progressive compression. For example, Taubin et al. [START_REF] Taubin | Progressive forest split compression[END_REF] proposed a progressive encoding based on Taubin and Rossignac's algorithm [START_REF] Taubin | Geometric compression through topological surgery[END_REF]. Cohen-Or et al. [START_REF] Cohen-Or | Progressive compression of arbitrary triangular meshes[END_REF] used techniques of sequential simplification by vertex suppression for connectivity, combined with position prediction for geometry. Alliez and Desbrun [START_REF] Alliez | Progressive compression for lossless transmission of triangle meshes[END_REF] proposed an algorithm based on progressive removal of independent vertices, with a retriangulation step under the constraint of maintaining the vertex degrees around 6. Contrary to the majority of compression methods, Gandoin and Devillers [START_REF] Gandoin | Progressive lossless compression of arbitrary simplicial complexes[END_REF] gave the priority to geometry encoding. Their algorithm, detailed in Sec. 3, gives competitive compression rates and can handle simplicial complexes in any dimension, from manifold regular meshes to triangle soups. Peng and Kuo [START_REF] Peng | Geometry-guided progressive lossless 3d mesh coding with octree decomposition[END_REF] took this paper as a basis to improve the compression ratios using efficient prediction schemes (sizes decrease by roughly 15%), still limiting the scope to triangular models.

Large meshes processing

Working with huge datasets implies to develop out-of-core solutions for managing, compressing or editing meshes. Lindstrom and Silva [START_REF] Lindstrom | A memory insensitive technique for large model simplification[END_REF] proposed a simplication algorithm based on [START_REF] Lindstrom | Out-of-core simplification of large polygonal models[END_REF] whose memory complexity is neither dependent on the size of the input nor on the size of the output. Cignoni et al. [START_REF] Cignoni | External memory management and simplification of huge meshes[END_REF] introduced a data structure called Octree-based External Memory Mesh (OEMM), which enable external memory management and simplification of very large triangle meshes, by loading only selected sections in main memory. Isenburg et al. [START_REF] Isenburg | Streaming meshes[END_REF] developed a streaming file format for polygon meshes designed to work with large data sets. Isenburg et al. [START_REF] Isenburg | Streaming compression of triangle meshes[END_REF] proposed a single-rate streaming compression scheme that encodes and decodes arbitrary large meshes using only minimal memory resources. Cai et al. [START_REF] Cai | Progressive out-of-core compression based on multi-level adaptive octree[END_REF] introduced the first progressive compression method adapted to very large meshes, offering a way to add out-of-core capability to most of the existing progressive algorithms based on octrees. The next section pays particular attention to out-of-core visualization techniques.

Visualization

Fast and interactive visualization of large meshes is a very active research field. Generally, a tree or a graph is built to handle a hierarchical structure which makes it possible to design out-of-core algorithms: at any moment, only necessary and sufficient data are loaded into memory to render the mesh. Level of detail, view frustum and occlusion culling are widely used to adapt displayed data to the viewport. Rusinkiewicz and Levoy [START_REF] Rusinkiewicz | Qsplat: a multiresolution point rendering system for large meshes[END_REF] introduced QSplat, the first out-of-core point-based rendering system, where the points are spread in a hierarchical structure of bounding spheres. This structure can easily handle levels of detail and is well-suited to visibility and occlusion tests. Several million points per second can thus be displayed using adaptive rendering. El-Sana and Chiang [START_REF] El-Sana | External memory view-dependent simplification[END_REF] proposed a technique to segment triangle meshes into view-dependence trees, allowing external memory simplification of models containing a few millions polygons, while preserving an optimal edge collapse order to enhance the image quality. Later, Lindstrom [START_REF] Lindstrom | Out-of-core construction and visualization of multiresolution surfaces[END_REF] developed a method for interactive visualization of huge meshes. An octree is used to dispatch the triangles into clusters and to build a multiresolution hierarchy. A quadric error metric is used to choose the representative point positions for each level of detail, and the refinement is guided by visibility and screen space error. Yoon et al. [START_REF] Yoon | Quick-vdr: Interactive view-dependent rendering of massive models[END_REF] proposed a similar algorithm with a bounded memory footprint: a cluster hierarchy is built, each cluster containing a progressive submesh to smooth the transition between the levels of detail. Cignoni et al. [START_REF] Cignoni | Adaptive tetrapuzzles: Efficient out-of-core construction and visualization of gigantic multiresolution polygonal models[END_REF] used a hierarchy based on the recursive subdivision of tetrahedra in order to partition space and guarantee varying borders between clusters during refinement. The initial construction phase is parallelizable, and GPU is efficiently used to improve frame rates. Gobbetti and Marton [START_REF] Gobbetti | Far voxels: a multiresolution framework for interactive rendering of huge complex 3d models on commodity graphics platforms[END_REF] introduced the far voxels, capable of rendering regular meshes as well as triangles soups. The principle is to transform volumetric subparts of the model into compact direction-dependent approximations of their appearance when viewed from a distance. A BSP tree is built, and nodes are discretized into cubic voxels containing these approximations. Again, the GPU is widely used to lighten the CPU load and improve performance. Cignoni et al. [START_REF] Cignoni | Batched multi triangulation[END_REF] proposed a general formalized framework that encompasses all these visualization methods based on batched rendering. Recently, Hu et al. [START_REF] Hu | Parallel view-dependent refinement of progressive meshes[END_REF] introduced the first highly parallel scheme for view-dependent visualization that is implemented entirely on programmable GPU. Although it uses a static data structure requiring 57% more memory than an index triangle list, it allows real-time exploration of huge meshes.

Combined Compression and Visualization

Progressive compression methods are now mature (the rates obtained are close to theoretical bounds) and interactive visualization of huge meshes has been possible for several years. However, even if the combination of compression and visualization is often mentioned as a perspective, very few papers deal with this problem, and the files created by visualization algorithms are often much larger than the original ones. In fact, compression favors a small file size to the detriment of fast data access, whereas visualization methods focus on rendering speed: the two goals oppose and compete with each other. It must be noted that there are a few counter-examples such as [START_REF] Oliveira | Pnorms: platonic derived normals for error bound compression[END_REF] where there is a reduction of redundancy without any runtime penalty, but it remains unusual. Among the few works that introduce compression into visualization methods, Namane et al. [START_REF] Namane | Qsplat compression[END_REF] proposed a QSplat compressed version by using Huffman and differential coding for spheres (position, radius) and normals. The compression ratio is approximately 50% compared to the original QSplat files, but the scope is limited to point-based rendering. More recently, Yoon and Lindstrom [START_REF] Yoon | Random-accessible compressed triangle meshes[END_REF] introduced a triangle mesh compression algorithm that supports random access to the underlying compressed mesh. However, their file format is not progressive and therefore inappropriate for global views of a model since it would require loading the full mesh.

The Starting Point: kd-tree Compression of Vertex Position and Hierarchical Connectivity Coding

In this section, we briefly present the method originally proposed by Gandoin and Devillers [START_REF] Gandoin | Progressive lossless compression of arbitrary simplicial complexes[END_REF]. The algorithm, valid in any dimension, is based on a top-down kd-tree construction by cell subdivision. The root cell is the bounding box of the point set and the first output is the total number of points. Then, for each cell subdivision, only the number of points in the first half-cell has to be coded, the content of the second half-cell being easily deduced. As the algorithm progresses, the cell size decreases and the data transmitted provide a more accurate localization of the points. The subdivision process iterates until there is no longer a non-empty cell with a side greater than 1 spatial unit, such that every point is localized with the full mesh precision.

As soon as the points are separated, the connectivity of the model is embedded (one vertex per cell) and the splitting process is run backwards: the cells are merged bottom-up and their connectivity is updated. The connectivity changes between two successive versions of the model and is encoded by symbols inserted between the geometry codes. The vertices are merged by the two following decimating operators:

• Edge collapse, originally defined by Hoppe et al. [START_REF] Hoppe | Mesh optimization[END_REF] and widely used in surface simplification, merges two adjacent cells under certain hypotheses. The endpoints of the edge are merged, which leads to the deletion of the two adjacent triangles (only one if the edge belongs to a mesh boundary) (see Fig. 1a).

• Vertex unification, as defined by Popović and Hoppe [START_REF] Popović | Progressive simplicial complexes[END_REF], is a more general operation that merges any two cells even if they are not adjacent in the current connectivity. The result is non-manifold in general (see Fig. 1b) and the corresponding coding sequence is approximately twice as big as an edge collapse coding sequence.

The way the connectivity evolves during these decimating operations is encoded by a sequence that enables a lossless reconstruction using the reverse operators (edge expansion and vertex split). If this lossless compression method reaches competitive ratios and can handle arbitrary simplicial complexes, it is not appropriate for interactive navigation through large meshes. Not only does its memory footprint make it strictly impracticable for meshes over one million vertices, but, even more constraining, the intrinsic design of the data structure imposes a hierarchical coding of the neighborhood relations that prevents any local refinement: indeed, given 2 connected vertices v and w in any intermediate level of the kd-tree, it is possible to find a descendant v i of v that is connected to a descendant w j of w. Therefore, in terms of connectivity, the cell containing v cannot be refined without refining the cell containing w. Consequently, random access and selective loading of mesh subsets is impossible: to visualize a single vertex and its neighborhood at a given LOD, the previous LOD must be entirely decoded. In the following, algorithms and data structures are presented that are based on the same kd-tree approach but remove these limitations.

Overview

The compression process of CHuMI begins by partitioning the space, creating a hierarchical structure called an nSP-tree in which the simplices are dispatched. In order to avoid SP-cells containing only a few triangles, an SP-cell must contain at least N min (user defined) triangles to be split. Moreover, to allow the independent decoding of any vertex with its neighborhood, the simplices whose vertices belong to several SP-cells are duplicated into each cell. We then traverse the SP-tree in postorder (a cell is processed after all its children have been processed).

In each SP-cell, incident kd-cells are progressively merged as described in Sec. 3, and the corresponding coding sequence is constructed, where geometry and connectivity codes are interleaved. Merging proceeds until the minimal precision of the SP-cell is reached: in the example given in Fig. 2, this corresponds to one kd-cell merging in each dimension for a regular SP-cell, and 3 mergings in each dimension for the root SP-cell. When all the children of an SP-cell s have been treated, the submeshes they contain are merged together (the duplicated simplices collapse) to form the submesh associated with s, which is then ready to be treated. This process runs until a single vertex remains in the root SP-cell.

The refinement process runs backward. Figure 2 illustrates its main stages: steps 1 to 6 are similar to [START_REF] Gandoin | Progressive lossless compression of arbitrary simplicial complexes[END_REF]: we start with a null precision and a single centered point. Then we sequentially read the geometry and connectivity codes from the compressed

0-bit precision 1-bit precision 2-bit precision 3-bit precision 5-bit precision 4-bit precision 1-bit precision on Y axis 2 bits on X axis 0-bit precision on Y axis 1 bit on X axis

Figure 2: Overview of the method file, which progressively refine the entire mesh. In this example, when precision reaches 3 bits (in practical cases, we generally wait until precision reaches 7 or 8 bits), the space is split into 4 subspaces, and the mesh subdivided into these subspaces. Remember that border triangles and their associated kd-cells have been duplicated during the compression stage. The 4 submeshes obtained become independent, each one corresponding to a subspace called an SP-cell whose content is directly accessible from an index pointing in the compressed file. Thus we can select one or several SP-cells and read their associated refinement codes in the compressed file in order to continue the refinement process locally, depending on the camera moves. When one of them reaches its maximum precision (4 bits in our example), it is split again. This process is repeated recursively until each SP-cell in the view frustum reaches a level of precision determined by its distance to the camera. Finally, we connect the SP-cells together, handling the possible differences in resolution to prevent border triangles from overlapping, and the content of these cells is sent to the graphic pipeline.

Out-of-Core Construction

Before describing the construction process itself, we begin with a few definitions:

An nSP-tree is a space-partitioning tree with n subdivisions per axis. Each node delimits a cubic subspace called the nSPcell. Since our meshes are embedded in 3D, an internal nSPcell has n 3 children (wherever it is in the tree), and in our case, the nSP-cells are evenly split, i.e. their n 3 children are the same size.

Each nSP-cell c contains vertices whose precision at runtime can vary according to the distance from c to the camera. The boundaries of the vertex precision in c are called minimal and maximal precision of c.

A top-simplex is a simplex (vertex, edge, triangle or tetrahedron) that does not own any parent, i.e. that does not compose any other simplex. In a triangle-based mesh, the top-simplices are mainly triangles, but in a general simplicial complex, they can be of any type.

Quantizing point coordinates

First, the points contained in the input file are parsed. A cubic bounding box is extracted, then the points are quantized to obtain integer coordinates relative to the box, according to the maximum level of precision p (number of bits per coordinate) requested by the user. It should be noted that this quantization step does not contradict the lossless characteristic of the method since p can be chosen as high as needed to match the object's initial precision. These points are written to a raw file (RWP) using a classical binary encoding with constant size codes, such that the i th point can be directly accessed by file seeking.

nSP-tree construction

An nSP-tree is built so that each leaf contains the set of the top-simplices lying in the corresponding subspace. The height of the nSP-tree depends on the mesh precision p (the number of bits per coordinate): the user can choose the maximal precision p r of the root, then the precision p l gained by each next level is computed from n (p l = log 2 n). For example, for p = 12 bits, p r = 6 bits and n = 4, the nSP-tree height is 4. Another condition is applied: an nSP-cell can be split only if the number of contained top-simplices is at least equal to a threshold N min . Consequently, the minimal precision of the leaves can take any value, whereas their maximal precision is always p. The content of the leaves is stored in a raw index file (RWI) where the d + 1 vertex indices that compose a d-simplex are written using a raw binary format to minimize the disk footprint while maintaining a fast and simple reading. Fig. 4 presents a 2D example of such an nSP-tree.

Since an nSP-cell is only subdivided when it contains more than N min triangles, the nSP-tree takes into account the local density of the mesh: sparse regions are hardly divided, whereas high density areas are heavily splitted. This simple scheme has been chosen because it combines adaptability and low encoding cost.

RWP and RWI file management

Even if the points and triangles are not ordered in the input file, it is likely that the different accesses to a given point of the RWP file during the nSP-tree construction stage will be close in time. Consequently, the algorithm uses a cache in the main memory, based on a hash table, to minimize the disk accesses to vertex positions. The hash keys are modulos of the point indices and guarantee a constant time access to cached points.

Since the final number of the triangles contained in an nSPcell is not known in advance, memory cannot be allocated statically for each cell in the RWI file. Consequently, the algorithm has to use block list management: the RWI file is split into constant size blocks (containing 100 triangles for instance), and the content of an nSP-cell is stored in a chained list of such blocks. As shown in Fig. 3, each nSP-cell keeps 3 numbers in main memory: the offset of its head block in the RWI file, the offset of the current position (the location where the next triangle will be written in the RWI file), and the number of triangles contained in the current block (to know when a new block must be allocated). Each block contains the offset of its successor (or -1 for the tail block of a list) and a triangle list composed of vertex indices pointing to the RWP file. Rather than writing in the RWI file for each triangle read in the input file, each nSP-cell manages a buffer containing the last added triangles, and the disk is accessed only when the cumulative size of all the buffers exceeds some fixed value b max . Then the buffers are flushed in the RWI file, from the largest to the smallest, to favor sequential writing on the disk, until their cumulative size is brought back to b max /2. When the content of an nSP-cell c reaches N min triangles, c is split and the triangles (in the buffer and in the block list of the RWI file) are spread out in the newly created nSP-cells. The blocks of c are released and added to a list of available blocks used for future allocations.

Choosing n

In order to obtain an integer for p l , n must be chosen as a power of 2. This guarantees that any descendant of a kd-cell contained in an nSP-cell c remains inside c. Furthermore, it prevents the kd-tree cells from overlapping two nSP-cells.

Simplex duplication

The first way to determine whether a simplex belongs to an nSP-cell c, consists in checking whether at least one of its incident vertices lies inside c. This test is fast and suitable for most cases, but it is imperfect since a cell can be overlapped by a simplex without containing any of its vertices. This occurs very rarely and can be perceived visually even more infrequently, but this could be an issue for certain specific applications, in which case a more accurate simplex to axis-aligned box intersection test could be performed.

Independently of the chosen test, a simplex of the original model can simultaneously belong to several nSP-cells; consequently, some simplices and their associated kd-cells must be duplicated during the distribution in the leaves to prevent holes during visualization. Although this technique degrades slightly compression ratios, it has been chosen for efficiency purposes. Drawings 7-8 and 10-11 in Fig. 2 show the space partitioning of an nSP-cell.

Out-of-Core Compression

In this section, we detail the successive steps of the compression process, then discuss various aspects of the implementation that contribute to reaching good performance.

To achieve competitive compression ratios, the algorithm requires a two-pass coding that makes use of two files: a semicompressed temporary file (SCT) and the final entropy coded file (FEC).

Structures Example

The RWI file contains a sequence of : Writing the header of the FEC file. Regardless of the statistical data that are gathered for the final compression, some parameters characterizing the original model and the nSP-tree can be written to the FEC file at this stage. The file begins with a header that contains all general data required by the decoder: the coordinate system (origin and grid resolution, to reconstitute the exact original vertex positions), the dimension d of the mesh, the number p of bits per coordinate, the number n of subdivisions per axis for the space partitioning, as well as the root and level precision, p r and p l .

Treating each cell of the nSP-tree

Then the nSP-tree is traversed in postorder (in red in Fig. 4) and the following treatment is applied to each nSP-cell c.

(a) Building the kd-tree and the simplicial complex. If c is a leaf: the top-simplices belonging to c are read from the RWI file and the associated vertices from the RWP file; then a kd-tree is built that separates the vertices (top-down process) and finally the simplicial complex based on the kd-tree leaves is generated.

(b) Computing the geometry and connectivity codes. The leaves of the kd-tree are sequentially merged following a bottom-up process. For each fusion of two leaves, a code sequence is computed that combines geometry and connectivity information. It should be noted that for rendering needs, this sequence is enriched by specifying the orientation of the triangles that collapse with merging. It is not a complete normal coding, but only a binary description of each triangle orientation. The vertex normals are computed in real-time during the rendering stage by averaging the triangle normals. The process is stopped when the minimal precision of c is reached. In large part, the ideas of [START_REF] Gandoin | Progressive lossless compression of arbitrary simplicial complexes[END_REF] are used to compute geometry and connectivity codes, with some adaptations or simplifications to guarantee high frame rates. As for the geometry, each cell subdivision is encoded to describe one of the 3 following events: (0) the first half-cell is non-empty, (1) the second half-cell is nonempty, or (2) both cells are non-empty. The theoretical cost of log 2 3 = 1.58 bits per subdivision is reduced to an average of 1.2 bits using entropy coding with per-level statistical data (high levels are composed essentially of 2 codes, whereas low levels contain mainly 0 and 1 codes). As for the geometry, we separate edge expansion from vertex split. Encoding both operations is the same as in [START_REF] Gandoin | Progressive lossless compression of arbitrary simplicial complexes[END_REF] (see Sec. 3), with the exception of an additional code describing the orientation of each newly created triangle. Using a simple prediction based on the triangle orientation in the neighborhood, the overcost is made negligible (around 0.01 bit per vertex). It is worth noting that none of the prediction schemes of [START_REF] Gandoin | Progressive lossless compression of arbitrary simplicial complexes[END_REF] were used, thus promoting decoding speed and visualization efficiency. However, for a compressiondriven application (or in a context where meshes are visualized on high-end computers), the original prediction schemes can definitely be included.

(c) Writing the codes in a temporary file. The codes obtained are written in a semi-compressed temporary file (SCT) using an arithmetic coder [START_REF] Amir | Introduction to arithmetic coding theory and practice[END_REF]. At this stage, we cannot write in the final entropy coded file (FEC) since statistical data are not yet available.

(d) Merging into the parent nSP-cell. We are now left with a kd-tree and a simplicial complex whose vertices have the minimal precision of c. To continue the bottom-up postorder traversal of the nSP-tree, the content of c must be moved to its parent. If c is the first child, its content is simply transferred to the parent. Otherwise, the kd-tree and simplicial complex of c are combined with the parent's current content, which implies detecting and merging the duplicated simplices (see Sec. 7.5).

Final Output

Once all nSP-cells have been processed, statistical data over the whole stream are available and efficient entropy coding can be applied. The probability tables are added to the header of the FEC file, then the nSP-cells code sequences of the SCT file are sequentially passed through an arithmetic coder. As geometry code frequencies vary with respect to the level of detail, the probabilities are computed independently for each level. In our case, dynamic probabilities would not improve this static per-level probability scheme since the independence of nSPcells would require the statistical data to be reinitialized too frequently. At the end of the file, a description of the nSP-tree structure is added that indicates the position of each nSP-cell in the file. This table allows the decompression algorithm to rebuild an empty nSP-tree structure that provides direct access to the code sequence of any nSP-cell.

Memory and Efficiency Issues

Memory management. To minimize memory occupation, each loaded nSP-cell only stores the leaves of its kd-tree. The internal nodes are easily rebuilt by recursively merging the sibling leaves. Furthermore, the successive refinements and coarsenings imply a great deal of creation and destruction of objects such as kd-cells, simplices and lists. A naive memory management would involve numerous memory allocations and deallocations, causing fragmentation and substantial slowdowns. To avoid this, memory pools are widely used to preallocate thousands of objects in one operation, and preference is given to object overwriting rather than deletion and recreation.

Common implementations of the simplicial complex build a structure where each d-simplex (d > 1) is composed of d+1 (d-1)-simplices. This structure is costly to maintain, and since we do not particularly need to keep this parent-to-child structure, we have chosen to store the following minimal data structure: simplicial complexes are represented as a list of top-simplices, each simplex containing references to the corresponding kdcells that give a direct link to its incident vertices. Conversely, each kd-cell stores a list of its incident top-simplices.

Finally, the external memory accesses are optimized to avoid frequent displacements of the hard disk heads. As seen in section 5.3, pagination and buffering techniques are used to write in the RWI file. During visualization, each time an nSP-cell is needed, its complete code is read from the FEC file and put into a buffer to anticipate upcoming accesses.

Multi-core parallelization. The nSP-tree structure is intrinsically favorable to parallelization. The compression has been multithreaded to benefit from the emerging multi-core processors. While file I/O remain single-threaded to avoid costly random accesses to the hard disk, vertex splits and unifications, edge expansions and collapses, nSP-cell splits and mergings can be executed in parallel on different nSP-cells. To favor effective sequential access to the hard drive and to avoid frequent locks of the other computing threads (file I/O are exclusive), the processing of an SP-cell starts by loading the list of contained triangles from RWI and RWP files into a buffer. Then the following steps (see Sec. 6.1) only require CPU and RAM operations and can be performed by one thread running parallel to the other threads. Experimentally, a global increase in performance is observed between 1.5 and 2 with a dual-core CPU, and between 2 and 3 with a quad-core CPU. This parallelization similarly benefits the decompression and visualization stage, although the performance increase is more difficult to estimate.

Decompression and View-Dependent Visualization

Decompression and visualization are strongly interlinked, since the data must be efficiently decoded and selectively loaded in accordance with the visualization context. In this section, each component of the rendering process is detailed.

Initialization

The first step consists in reading general informations about the mesh from the FEC file header, then using the offset table situated at the end of the file to build the nSP-tree structure. Once done, each nSP-cell only contains its position in the FEC file. To complete the initialization process, the kd-tree of the root nSP-cell is created at its simplest form, i.e. one root, and likewise the associated simplicial complex is initialized to one vertex.

Real-Time Adaptive Refinement

At each time step of the rendering, the mesh is updated to match the visualization frame. The nSP-cells in the view frustum are refined or coarsened so that the maximum error on the coordinates of any vertex guarantees that its projection respects the imposed display precision. The other cells are coarsened to minimize memory occupation, and are not sent to the graphic pipeline to reduce GPU computing time. According to this viewpoint criterion, a list of nSP-cells to be rendered with their associated level of precision is maintained. A cell that reaches its maximal precision is split, and the content of each child is efficiently accessed thanks to the offset contained in the nSPtree structure. Conversely, if an nSP-cell needs to be coarsened at a level below its maximal precision, its possible children are merged together.

Multiresolution Transitions

Once the nSP-cells that need to be rendered are known, two main issues arise, both concerning the borders between nSPcells. First, the simplices being duplicated in different nSP-cells may overlap. Second, visual artifacts can be caused by adjacent nSP-cells with different levels of precision. Only border topsimplices can cause these problems; the others can be straight rendered using their representative points. Border top-simplices are composed of at least one vertex belonging to another nSPcell. The following algorithm is applied to each border topsimplex s:

1. Let c be the nSP-cell to which s belongs, p c the current precision of c, and l the list of nSP-cells to render. Let N be the number of kd-cells (i.e. the number of vertices) composing s, c i (for i in 1, . . . , N) the nSP-cell in which lies the i th kd-cell k i composing s, and p i the current precision of c i (if c i l then p i = p c ). 2. If there is an i in 1, . . . , N such that p i > p c or (p i = p c and c i index > c index), s is removed and will not be rendered.

3. Else, for each kd-cell k i such that c i ∈ l and c i c, the kdcell k i of c i containing k i is searched. The representative point of k i is used to render s.

Thus, a smooth and well-formed transition is obtained between nSP-cells that have different levels of precision (see Fig. 5).

Encoding Geometry Before Connectivity

The main drawback of this algorithm is its tendency to produce many small triangles in the intermediate levels of precision. These triangles, which usually have a size near the requested screen precision, slow down the rendering without resolving the well-known block effects inherent to the kd-tree or octree approaches (see Fig. 9, left). This issue is addressed by encoding geometry earlier than connectivity in the FEC file: the decoder will hear of the point splits m bits before the associated connectivity changes. Consequently, a kd-cell c will be composed of a connectivity code that describes the creation of the representative vertex of c, followed by multiple geometry codes that describe the positions of the child vertices of c. The level of the descendants coded in a kd-cell depends on m. For a 3D mesh and a geometry advance of m bits, 3 m levels of children are contained in each kd-cell. For each vertex displayed, its future children are thus known and their center of mass can be taken as the new refined position of the representative vertex (see Fig. 6).

(c) (b) (a)

Figure 6: Example of rendering with the geometry advance: comparison between (a) 2 bit precision for both geometry and connectivity, (b) 2 bits for connectivity and 3 bits for geometry (2 + 1 early bit), (c) 3 bits for both Figure 7 shows an example of an nSP-tree with a geometry advance of 1 bit. During the split of an nSP-cell, some kd-cells are duplicated in several children, such as kd-cell 11 in this example. The duplicates of the same kd-cell may evolve differently, since their neighborhood depends on the nSP-cell they belong to. Thus, their codes may differ, such as G11 and G11b. The parent kd-tree can store the early geometry codes of only one nSP-cell child. We choose the child containing the original kd-cell. For the other nSP-cells, the missing geometry codes are added at the beginning of the sequence. For instance, in our example, nSP-cell #0 stores G11, G23 and G24, corresponding to kd-cells 11, 23 and 24 of nSP-cell #2; missing G11b and G23b are stored at the beginning of nSP-cell #1 (bold red box). Similarly, since the root nSP-cell (#0) has no parents, a few geometry codes must be added at its beginning. Figures 8 and9 illustrate the benefits of this technique in terms of rendering. The rate-distortion curve, constructed using the L2 distance computed by the METRO tool [START_REF] Cignoni | Metro: Measuring error on simplified surfaces[END_REF], shows a significant improvement in the intermediate levels of detail. Ci / Gj+k i Kd-cell n°i containing connectivity codes of kd-cell i and geometry codes of kd-cells j and k Gm+n+... At the beginning of an nSP cell , the geometry codes that were not in the parent nSP -cell and that are not in the following kd -tree are stored

C0 / G3+5 nSP-cell #0 C1 / G8 C2 / G11 C3 / G17 C5 / G23+24 G0+1+2 nSP-cell #1 nSP-cell #2 0 1 2 3

Efficient Adaptive Rendering

To attain good performance in terms of smooth real-time rendering, particular attention was paid to a number of aspects of the implementation. During an nSP-cell split, kd-cells and top-simplices are transferred and sometimes duplicated in the children. A costly part of this split consists in determining the child to which each kd-cell must be transferred and the possible children where it must be copied (when it is incident to a border top-simplex). Rather than computing this information for numerous kd-cells just before the split, we compute it as soon as the kd-cell precision suffices to determine its future containing nSP-cell, store it and transfer it directly to the descendant kd-cells. Likewise, from this information, the progression of a top-simplex during the next nSP-cell split can be deduced. To smooth the computation load over time, each top-simplex is tested either after its creation (due to a kd-cell split) or after its duplication (due to an nSP-cell split).

Conversely, merging two nSP-cells implies deleting all the duplicated objects. To optimize this step and quickly determine these objects, each kd-cell and top-simplex stores the kd-tree level it was created in and whether it was created after an nSPcell split. Another issue is that top-simplices moved into the parent nSP-cell may refer to a duplicate kd-cell that must be deleted. Since we cannot afford to look for the corresponding original kd-cell among all the existing kd-cells, each duplicate stores a pointer to its original cell. To keep this working properly, the coarsening and refinement processes must be reversible: by refining an nSP-cell and then coarsening it back to its minimal precision, the same kd-cells must be obtained.

Computing the connectivity changes associated with kd-cell splitting or merging are costly operations on which we have focused to ensure a fast rendering. First, an efficient order relation is needed over the kd-cells: rather than the geometric position, the kd-index order has been chosen, which costs a single comparison and gives a natural order relation over the top-simplices. In addition, it has been shown that each topsimplex stores the kd-tree level it was created in. This value is also used to speed up edge collapse and vertex unification during coarsening: since it points out the top-simplices that must be removed after merging 2 kd-cells, we only need to replace the removed kd-cell with the remaining kd-cell in the remaining incident top-simplices, then check whether any triangles degenerate into edges.

In order to improve the visualization speed, we implemented prefetching, which can be used either to preload the whole file into memory and thus avoid any subsequent hard drive access (if the file fits into main memory), or to prefetch upcoming nSP-cells from disk by trying to predict the camera position at the next time step. However, this hardly improves the overall performance: for example, the computing time to obtain a view of a close-up on the St. Matthew face is 5.8s without any prefetching, whereas it is 5.6s with full prefetching (whole file preloaded). This experiment shows that, with a cumulative cost of 0.2s over 5.8s (3.4%), the loading process is almost negligi-ble. This can be explained by two factors: the first one is that we carefully optimized external data access (buffered and sequential), the second one is that compression drastically reduces the amount of data we need to read. This last point argues in favor of compression, and we can assume that it will become more and more important, in view of the increasing gap between processing units and data access performances [START_REF] Hennessy | Computer Architecture -A Quantitative Approach[END_REF].

Experimental Results

We cannot fairly compare our method to previous work, since, to our knowledge, no other method combines the features of compression and interactive rendering. However, for reference, we provide comparisons to state-of-the-art compression methods (single-rate or progressive, in-core or out-of-core) in Sec. 8.1, then to state-of-the-art visualization methods in Sec. 8.2. The results presented were obtained from a C++/OpenGL implementation of the method running on a PC with an Intel Q6600 QuadCore 2.4Ghz CPU, 4GB DDR2 RAM, 2 RAID0 74GB 10000 RPM hard disks, and an NVIDIA GeForce 8800 GT 512MB video card. Regarding the parameter settings, n must be big enough to allow the selection of small parts of the mesh in a few nSP-cell splits, but small enough to avoid displaying too much nSP-cells simultaneously (in which case borders computation, for instance, could be excessively expensive). A small N min increases the number of duplicated simplices and degrades the compression performance, whereas a big N min leads to a smaller nSP-tree, to the detriment of the multiresolution capability. n = 4 and 5000 ≤ N min ≤ 8000 appear to be good trade-offs.

Compression

Table 1 presents results from the out-of-core compression stage. p r was set at 7 bits so that the number of triangles in the root nSP-cell remains small enough to allow this coarse version of the mesh to be loaded in main memory. For each model, we first indicate the number of triangles, the PLY file size, and the size of the raw binary coding, which is the most compact version of the naive coding: if v is the number of vertices, t the number of triangles and p the precision in bits for each vertex coordinate, this size in bits is given by 3pv for the geometry +3 t log 2 v for the connectivity. Then the compression rate is given (as the ratio between the raw size and the FEC size) with the total number of nSP-cells. Compression times are also provided, the column headers referring to the sections that describe the corresponding steps of the algorithm. As expected, the computation of code sequences (see section 6.1) is the most expensive part. However, this step, which is directly proportional to the number of vertices and triangles, benefits from our multithreaded implementation. Finally, the table details the memory usage (including temporary disk usage), attesting that the method is actually out-of-core. It is worth noting that the spatial coherence of the input model is important and can explain, for instance, why the step 6.1 is shorter for Lucy than for the T8 model, where Lucy contains more triangles. However, we found that the scanned models we used showed in general a good spatial coherence. Table 2 shows the results of our algorithm in terms of bits per vertex compared to those of [START_REF] Touma | Triangle mesh compression[END_REF] (a reference in-core singlerate method), the original in-core multiresolution algorithm [START_REF] Gandoin | Progressive lossless compression of arbitrary simplicial complexes[END_REF], the single-rate out-of-core method [START_REF] Isenburg | Out-of-core compression for gigantic polygon meshes[END_REF] and the single-rate outof-core random access method [START_REF] Yoon | Random-accessible compressed triangle meshes[END_REF]. The first three are pure compression methods, while the fourth provides a trade-off between compression and random accessibility, very useful for applications that need mesh traversal but not well-suited to visualization. The two non-compressive visualization methods [START_REF] Cignoni | Adaptive tetrapuzzles: Efficient out-of-core construction and visualization of gigantic multiresolution polygonal models[END_REF] and [START_REF] Gobbetti | Far voxels: a multiresolution framework for interactive rendering of huge complex 3d models on commodity graphics platforms[END_REF] have been added for comparison purposes. Unfortunately, we were not able to compare with [START_REF] Cai | Progressive out-of-core compression based on multi-level adaptive octree[END_REF] (to our knowledge, the sole out-of-core and progressive compression algorithm) since the authors only give connectivity rates. Compared to compression methods that do not allow interactive visualization (first three columns), an average extra cost ranging from 15% to 80% is observed, mainly due to the redundancy induced by space partitioning (5-7% of the original triangles are duplicated) and the fact that complex prediction is bypassed to guarantee high frame rates during rendering.

To show that the method behaves well with point clouds (possibly unbalanced) and does not rely on the mesh connectivity, two examples of unstructured models have been added to our experimental pool: Bali and Wallis.

Decompression and Visualization

Figure 10 presents a few examples of different LODs of the St. Matthew model with the times observed to obtain a given view distance by refinement or coarsening of the previous one (or from the beginning of the visualization for view (a)). They may seem high but it must be noted that during navigation, the [START_REF] Jamin | CHuMI videos[END_REF] transition from one view to the next one is progressive and the user's moves are usually slow enough to allow smooth refinements. This is why the best way to appreciate the algorithm's behavior is to manipulate the viewer or watch the accompanying videos [START_REF] Jamin | CHuMI videos[END_REF]. All the results have been produced with a maximal screen space error tolerance of 1 pixel. To appreciate the data's decoding speed, Fig. 11 presents the number of triangles rendered per second during the St. Matthew video. Our decoder can render up to 200 million triangles per second (tps), with an average of about 173 million when it runs at full speed (between t = 30 to t = 65sec. in Fig. 11). For comparison purposes, [START_REF] Lindstrom | Out-of-core construction and visualization of multiresolution surfaces[END_REF] renders up to 3 Mtps, [START_REF] Cignoni | Adaptive tetrapuzzles: Efficient out-of-core construction and visualization of gigantic multiresolution polygonal models[END_REF] an average of 70 Mtps, and [START_REF] Gobbetti | Far voxels: a multiresolution framework for interactive rendering of huge complex 3d models on commodity graphics platforms[END_REF] an average of 45 Mtps. Of course, this comparison would not be fair without taking into account each method's year of publication and the progress in hardware performance. In addition, [START_REF] Cignoni | Adaptive tetrapuzzles: Efficient out-of-core construction and visualization of gigantic multiresolution polygonal models[END_REF] and [START_REF] Gobbetti | Far voxels: a multiresolution framework for interactive rendering of huge complex 3d models on commodity graphics platforms[END_REF] update the model according to viewing distance and direction, whereas CHuMI uses a distance criterion only. On the other hand, it must be remembered that none of these visualization methods achieves compression: for instance, [START_REF] Cignoni | Adaptive tetrapuzzles: Efficient out-of-core construction and visualization of gigantic multiresolution polygonal models[END_REF] and [START_REF] Gobbetti | Far voxels: a multiresolution framework for interactive rendering of huge complex 3d models on commodity graphics platforms[END_REF] enlarge the original raw files, respectively, by 10% and 80% on average, as shown in Tab. 2. Fig. 12 shows the number of triangles which are created and deleted per second; zooming in leads to the creation of up to 310,000 triangles per second, while up to 565,000 triangles per second can be destroyed during zooming out. It must be added that the small polyhedral holes that can be observed on the captures and especially on the videos are not artifacts of the viewer but come from the original meshes. It is worth noting that the transitions between the k and k + 1 levels of precision of an nSP-cell, easily perceptible on the video, occur quite rarely (more precisely, when the zoom factor on the SP-cell is doubled, since one additional bit on point coordinates doubles the rendering accuracy). Consequently, in most cases, if such a transition does not occur during a zoom, this does not mean that the refinement process is stuck but that the threshold is not reached.

Conclusion and Future Work

As seen in Sec. 8, the data structures and algorithms presented in this article achieve good results compared to the state of the art, at the same time as an out-of-core lossless compression method, and as an interactive visualization method for arbitrary untextured meshes (in terms of type and size), thus offering a valuable combination for numerous applications. However, we are currently improving several points, which should remove the method's main limitations. Our first perspective aims at reducing latency in the navigation stage. To achieve better performance in this respect, at least three paths can be explored: cache-oblivious optimizations, more GPU-friendly data structures such as patch-based LOD, and the intensive use of the GPU via shader programming. Furthermore, we believe that it is still possible to improve the compression ratios by introducing new predictive schemes. Secondly, although theoretically, the method can address meshes in any dimension, some important optimizations must be done to handle volume meshes in good conditions. For instance, the edge collapse and vertex unification operators have to be efficiently extended to tetrahedra. In addition, occlusion culling becomes a major issue in the volume case. Finally, CHuMI is based so far on a pure distance criterion to set the precision of displayed triangles. We think that it could be favorably extended to use a hierarchy of surface patch normals and thus enable updates according to view direction as well as distance.

supported by the French Centre National de la Recherche Scientifique (CNRS).

Figure 1 :

 1 Figure 1: (a) Edge collapse, (b) Vertex unification

  Refinement of the root nSP -cell... Split of the root nSP -cell... the nSP-cell #4 only intersects the view frustum ... If the nSP-cell #18 only intersects the view frustum ...

Figure 4 :

 4 Figure 4: 2D nSP-tree example with p = 8 bits, p r = 6 bits, n = 2 and N min = 1000

Figure 3 :

 3 Figure 3: RWI file structure

Figure 5 :

 5 Figure 5: Rendering of the border top-simplices between two nSP-cells with different LODs

Figure 7 :

 7 Figure 7: Example of a kd-tree with geometry before connectivity: 2D case, 1 bit earlier (i.e. 2 kd-tree levels)

Figure 8 :Figure 9 :

 89 Figure 8: Comparison of rate-distortion curves

Figure 11 :

 11 Figure 11: Triangles per second (TPS) and frames per second (FPS) over the time during the accompanying St. Matthew video[START_REF] Jamin | CHuMI videos[END_REF]

Figure 10 :

 10 Figure 10: St. Matthew refinement times: 0.7 s to display (a), 1.2 s from (a) to (b), 2.0 s from (b) to (c), 1.7 s from (c) to (d); coarsening times: 0.8 s from (d) to (c), 0.6 s from (c) to (b), 0.4 s from (b) to (a).

Figure 12 :

 12 Figure 12: Stacked column chart showing the number of triangles decoded/deleted per second over the time during the St. Matthew video[START_REF] Jamin | CHuMI videos[END_REF]

Table 2 :

 2 Comparison of compression ratios in bpv

	Models	Vertices	p	[32]	[4]	[1]	[28]	[22]	[23]	CHuMI
	Bunny	35 947	12	-	-17.8	-	-	-	27.0
	Horse	19 851	12	19.3	-20.3	-	-	-	29.3
	Dino	14 050	12	19.8	-	-	31.8	-	-	28.7
	Igea	67 173	12	17.2	-	-	25.0	-	-	25.9
	David 2mm	4 128 028	16	-14.0	-	-	306.2	-	22.3
	UNC Powerplant	10 890 300	16	-14.1	-	-	-	-	16.3
	Lucy	13 797 912	16	-16.5	-	-	-	-	25.9
	David 1mm	27 405 599	16	-13.1	-	-	282.3	-	22.2
	St. Matthew	166 933 776	16	-10.7	-	22.9	282.1	508.0	19.4

Acknowledgments

The models are provided courtesy of The Digital Michelangelo Project, Stanford University Computer Graphics Laboratory and UNC Chapel Hill. The EDF Tx are parts of the Dancers column project: Electricité de France sponsored the Dancers column project for the benefit of Ecole Française d'Athènes by offering its technical support and financing the partnerships. The Greek authorities of the Museum of Delphi gave access to the column. The Insight team was selected for the survey and the registration of the data, with a contribution of Paolo Cignoni from Visual Computing Lab of the Inst. of Information Science and Technologies of the Italian National Research Council. The Department of Computer Science at University of Bristol kindly provided a scanner. This research was part of the ACI project Eros3D and the ANR project Triangles,


Jianghao Chang
email: jianghao.chang@umontreal.ca

Benoît Alain
email: benoit.alain@umontreal.ca

Victor Ostromoukhov

Structure-Aware Error Diffusion

: Halftone images produced with different techniques. Visual quality of our images is comparable with that produced with structureaware halftoning. At the same time, our method is two to three orders of magnitude faster than structure-aware halftoning (less than one second vs. hundreds of seconds). This difference is even bigger for larger images. Important notice: all images in this article have been prepared to be visualized on a graphic display, at 300 dots per inch. To visualize them correctly using Acrobat reader, please set the menu Edit/Preferences/Page Display/Resolution to "Custom resolution 100 pixels/inch", and visualize using 300% zoom.

Introduction

Digital halftoning is a technique for rendering continuous-tone images on devices having a limited number of available tone values. A typical application for digital halftoning algorithms is printing. In fact, millions of printers rely on digital halftoning in order to reproduce on paper a wide range of images: photos, business graphics, synthetic images, etc. Image quality and speed are two crucial issues to which any digital halftoning must be very sensitive. The algorithm that we introduce in this paper addresses both issues: it combines the quality of one of the best known techniques with the speed of conventional error-diffusion methods.

According to [START_REF] Pappas | Model-Based Digital Halftoning[END_REF][START_REF] Lee | The hybrid screen: improving the breed[END_REF], there are three classes of digital halftoning algorithms:

• Class I algorithms employ point processes, also called dithering [START_REF] Bayer | An optimum method for two-level rendition of continuous-tone pictures[END_REF][START_REF] Ulichney | Digital Halftoning[END_REF]]. • Class II algorithms use neighborhood processes, such as error diffusion [START_REF] Floyd | An adaptive algorithm for spatial grey scale[END_REF]]. • Class III algorithms make use of iterative algorithms such as global search or direct binary search (DBS) [START_REF] Analoui | Model based halftoning using direct binary search[END_REF][START_REF] Baqai | Halftoning via direct binary search using analytical and stochastic printer models[END_REF]].

Class I algorithms are the fastest, but their quality is insufficient for many practical applications, such as driving popular ink-jet printers or visualization on displays. Class III algorithms are the best from the quality viewpoint. Nevertheless, their complexity is too large, and they are impractical for most real-world applications.

Class II algorithms are the most popular ones: they combine quality and speed. Nevertheless, the quality of images produced with Class II algorithms is far from perfect; it is clearly inferior to that of Class III.

Recently, impressive structure-aware halftoning has been pro-posed [START_REF] Pang | Structure-aware halftoning[END_REF]. This Class III algorithm performs iterative optimization according to a perceptual texture-sensitive metric.

The method produces excellent results, but at very high cost: images of practical interest require minutes of processing time, which is obviously unacceptable for most practical applications.

The goal of the present article is to show that a Class II structureaware algorithm can be built. The name of our system-structureaware error diffusion-reflects its essential nature: on one hand, it inherits the essential features of structure-aware halftoning. On the other hand, it is basically a Class II (error-diffusion) algorithm. Consequently, our algorithm benefits from the best features of both: it has quality comparable to that of Class III algorithms, and the speed of Class II ones.

The rest of the article is organized as follows: In Section 2 we give a brief review of the main relevant prior art. In Section 3 we give an overview of the system. In Section 4 the results are presented and discussed. Finally, in Sections 5 and 6 we discuss limitations and draw conclusions.

Previous work

The error-diffusion algorithm, first introduced by Floyd and Steinberg [1976], has been extensively studied in the literature [START_REF] Ulichney | Digital Halftoning[END_REF][START_REF] Kang | Digital Color Halftoning[END_REF][START_REF] Sharma | Digital Color Imaging Handbook[END_REF]. We consider that this research area is well-established and well-documented. For this reason, we explore in this section only the bibliography directly related to our work. We refer interested readers to excellent and very complete textbooks mentioned hereabout.

The main advantage of the initial version of the algorithm is its simplicity combined with fairly good overall visual quality of the produced binary images. The basic algorithm has a number of inherent drawbacks: first, the algorithm produces identifiable visually harmful artifacts in highlights and in dark areas (worm artifacts). Moreover, at certain intensity levels, close to 1/2, 1/3, 2/3, etc., patches of regular structure may appear. These problems have been addressed in [START_REF] Ostromoukhov | A simple and efficient error-diffusion algorithm[END_REF]] by introducing intensity-dependent variable diffusion coefficients. This technique has been further improved in [START_REF] Zhou | Improving mid-tone quality of variable-coefficient error diffusion using threshold modulation[END_REF], by introducing intensity-dependent noise, in order to break regularities near critical intensity levels such as 1/2, 1/3, 2/3, etc. Both papers use very similar optimization techniques for finding near-optimal intensitydependent diffusion coefficients. This approach became a popular standard base (see [START_REF] Pang | Structure-aware halftoning[END_REF]). For this reason, we shall designate by the term "standard error diffusion" a generic implementation of the method described in [START_REF] Ostromoukhov | A simple and efficient error-diffusion algorithm[END_REF]] and improved in [START_REF] Zhou | Improving mid-tone quality of variable-coefficient error diffusion using threshold modulation[END_REF]. In fact, we build our system on top of this "standard error diffusion": we effectively use intensity-dependent diffusion coefficients. In addition to this "standard" component, we use a set of modified diffusion filters, as explained in Section 3.2. Modified diffusion filters depend on local spectral content of the image (local frequency/orientation/contrast); the proportion between "standard" and modified diffusion filters is adjusted during the calibration phase, as explained in Section 3.3. [START_REF] Eschbach | Error-diffusion algorithm with edge enhancement[END_REF] improved the basic error diffusion mechanism by modulating the threshold process with the edge response, in order to improve the standard error diffusion in the textured area. They used the inverse of the input image as the threshold. Their approach (also known as edge-enhancement) has been further improved and generalized by [START_REF] Hwang | Improved edge enhanced error diffusion based on first-order gradient shaping filter[END_REF][START_REF] Kwak | Edge-enhanced error diffusion halftoning using human visual properties[END_REF]]. The latter considers both local luminance average and variation. However, even these improved versions of edgeenhancement contain a major drawback: the strength of the enhancement is controlled by a unique global coefficient. Consequently, the method is sensitive to only one particular sub-range of frequencies and contrasts. In complex images with mixed frequency/contrast content this method can fail. This particular drawback has been explained in [START_REF] Pang | Structure-aware halftoning[END_REF]]. In accordance with [START_REF] Pang | Structure-aware halftoning[END_REF], we give more evidence of fundamental insufficiency of the edge-enhancement technique, in the most advanced version described in [START_REF] Kwak | Edge-enhanced error diffusion halftoning using human visual properties[END_REF]]. Basic and generalized edge-enhancement techniques show one additional drawback when repetitive structures are present in the image. In fact, harmful undesirable Moiré-like patterns may appear . These harmful artifacts are clearly visible in Figure 2, no matter whether strong or weak enhancement coefficients were used. Once again, our method does not suffer from this harmful effect, as shown in Figure 2.

In the present paper, we further improve the ideas of [START_REF] Kwak | Edge-enhanced error diffusion halftoning using human visual properties[END_REF]]. We do use image-dependent threshold modulation, but we replace the unsharp filter by a frequency-and orientation-dependent Gabor filter as explained in Section 3.1. Moreover, we modulate the strength of the threshold modulation by an amount determined during the calibration phase, as explained in Section 3.3. Our method provides a consistent improvement over [START_REF] Kwak | Edge-enhanced error diffusion halftoning using human visual properties[END_REF]], even though visual difference between the results produced with both methods may appear rather subtle in some cases.

More recently, [START_REF] Pang | Structure-aware halftoning[END_REF] introduced a new, more sophisticated structure-aware halftoning method. This method preserves both the global tone and the texture information, by using, within a complex optimization process, a combined metric which takes into account tonal and structural fidelity. This method generates results of excellent quality, one of the best known to date. However, as we mentioned before, this method belongs to Class III, and is therefore very slow. We use [START_REF] Pang | Structure-aware halftoning[END_REF] as the target reference for quality for our Class II structure-aware error diffusion.

Local frequency analysis

Classical Fourier analysis [START_REF] Bracewell | Fourier Transform and Its Applications[END_REF]] represents an image as a sum of sine waves of different orientations, periods, amplitudes, and phases. Usually, a Discrete Fourier Transform (DFT) [START_REF] Bracewell | Fourier Transform and Its Applications[END_REF]] is performed on images represented by discrete sets of intensities. Frequencies describe the images' essential nature: high frequencies correspond to fine details whereas low frequencies represent large-scale variations. Natural images may have heterogeneous frequency content.

While DFT provides a global analysis, many image processing applications need local analysis. Some existing techniques have been designed to that end [START_REF] Turner | Texture discrimination by gabor functions[END_REF][START_REF] Simoncelli | The steerable pyramid: A flexible architecture for multi-scale derivative computation[END_REF].

Our error-diffusion method uses local frequency content (dominant frequency/orientation/contrast within a small window) in order to get better sensitivity to structures present in the image. In presence of a complex frequency content, our method takes into account the most prominent frequency component only.

Local spectral analysis is not part of our contribution. Any analysis tool can be used as long as it detects the main frequency content.

Local frequency content analysis described in Chapter 13 of [START_REF] Jahne | Practical Handbook on Image Processing for Scientific and Technical Applications[END_REF]] is well-documented and satisfies our needs. Let us briefly summarize their technique. A gradient vector ∆g is computed for each pixel of the image. Finding the dominant orientation within a window is done by maximizing the term ΣN w(∆g • k0) 2 over the neighborhood, where w is a windowing function and k0 is the local orientation vector. Then, the DFT of the neighborhood window is computed and the dominant frequency and contrast are extracted, using known dominant orientation. For more details, refer to [START_REF] Jahne | Practical Handbook on Image Processing for Scientific and Technical Applications[END_REF]].

Original

Std. error diffusion Structure-Aware Halftoning [START_REF] Zhou | Improving mid-tone quality of variable-coefficient error diffusion using threshold modulation[END_REF] [ [START_REF] Pang | Structure-aware halftoning[END_REF] Edge-enhancement (strong coeff.) Edge-enhancement (weak coeff.)

Our method [START_REF] Kwak | Edge-enhanced error diffusion halftoning using human visual properties[END_REF] [ [START_REF] Kwak | Edge-enhanced error diffusion halftoning using human visual properties[END_REF]] 3 Structure-aware error diffusion

In this section, we provide an overview of our structure-aware errordiffusion algorithm.

The conventional error-diffusion algorithm processes each pixel of an input image in a particular sequential order (typically scanline or serpentine). The intensity of the pixel is compared with a thresh-old (constant or position-dependent). The result of that comparison determines the intensity of the output pixel, 0 or 1; the quantization error is distributed to the unprocessed neighboring pixels according to a particular diffusion filter.

In our method, uniform thresholding is adaptively replaced by thresholding from a parametrized anisotropic Gabor filter, and standard diffusion filter is adaptively replaced by a parametrized Gaussian filter. The adaptive transition from standard to structure-driven error diffusion is controlled by a weight parameter w.

During the offline calibration step, we select a perceptually optimal set of parameters for each of a large number of combinations of frequency/orientation/contrast (typically, 6 3 combinations). In the run-time part of our method, analyzed frequency content of each pixel then serves as index to the calibration table, and error diffusion is performed in accordance with the retrieved parameters. Figure 3 gives a schematic representation of our run-time architecture; pseudo-code is given in Appendix A.

Threshold modulation

Threshold modulation is a powerful tool for edge enhancement in error diffusion [START_REF] Knox | Edge enhancement in error diffusion[END_REF][START_REF] Eschbach | Error-diffusion algorithm with edge enhancement[END_REF][START_REF] Hwang | Improved edge enhanced error diffusion based on first-order gradient shaping filter[END_REF][START_REF] Kwak | Edge-enhanced error diffusion halftoning using human visual properties[END_REF]]. In this paper, we bring two important modifications to previous threshold modulation schemes. First, our modulation depends on the local frequency content of the image (local dominant frequency/orientation/contrast). Second, we calibrate the amount of threshold modulation in such a way that for each specific combination of dominant frequency/orientation/contrast, our error-diffusion system produces halftone output which is perceptually close to the source image.

[ [START_REF] Kwak | Edge-enhanced error diffusion halftoning using human visual properties[END_REF]] applied an isotropic unsharp mask filter on the input image for calculating the threshold matrix. This improves the original edge enhancement [START_REF] Knox | Edge enhancement in error diffusion[END_REF][START_REF] Eschbach | Error-diffusion algorithm with edge enhancement[END_REF], but, as [START_REF] Pang | Structure-aware halftoning[END_REF]] pointed out, the method suffers from blurriness because important and unimportant frequencies are equally amplified.

More recently, [START_REF] Kang | Flow-based image abstraction[END_REF]] remarked that feature enhancement using isotropic filters such as difference-of-Gaussians does not preserve the sense of "directedness", making aggregates of light edge pixels look less like lines.

Our method relies on the hypothesis that anisotropic filters should perform better than their isotropic counterpart in presence of anisotropic image structure. We chose the Gabor filter [START_REF] Gabor | Theory of communication[END_REF]] to accomplish this task for its very simple spectral properties. Since anisotropy is a property that varies throughout the image, the parameters of the Gabor filter change locally. We checked our hypothesis, and found that our method produced very satisfactory results.

Let us describe the threshold modulation part of our structure-aware error diffusion. Standard error diffusion uses constant threshold.

We use threshold modulation determined as a linear interpolation between two terms: (1) frequency-content-dependent filtered version FG of input image Img, and (2) noise modulation m, as in [START_REF] Zhou | Improving mid-tone quality of variable-coefficient error diffusion using threshold modulation[END_REF].

The term FG is the result of a local convolution with a suitable Gabor kernel and input image Img modulated by a variable strength parameter β:

FG(x, y) = β • i j Img(x -i, y -j) • Gabor θ,f (i, j)
where the function Gabor θ,f is defined from local frequency f and orientation θ given by image analysis:

Gabor θ,f (x, y) = exp - x 2 + y 2 2σ 2 G cos f x + c, x = x cos θ + y sin θ, y = -x sin θ + y cos θ,
and c is a constant that nulls the sum of the filter coefficients. In our current implementation, a fixed value σG = 1.6 was used.

The amount of threshold modulation β constitutes the main balancing factor for sensitivity to different natures of structural content. It is controlled carefully at the calibration step for the whole range of frequencies and contrasts. Since the image analysis described in Section 2.1 provides a finite number of combinations of frequency/orientation, all Gabor kernels were pre-calculated for fast run-time access (see Appendix A).

Modified diffusion filter

The method described in [START_REF] Ostromoukhov | A simple and efficient error-diffusion algorithm[END_REF]] provides a collection of 3-coefficient error-diffusion filters that produce best results on regions of constant tone. However, this method was conceived for low-frequency images and may blur high-frequency details. Notice that the midtone-improved version of the standard error diffusion [START_REF] Zhou | Improving mid-tone quality of variable-coefficient error diffusion using threshold modulation[END_REF] suffers from the same drawback, as shown in Figures 1 and7. Ulichney mentioned in his book [Ulichney 1987] that larger filter size can allow sharper feature preservation. Also, [START_REF] Marcu | Halftoning by back error compensation[END_REF] explored the use of analytical coefficient distribution, in particular that of a Gaussian kernel.

According to our hypothesis, diffusion coefficients following a distribution based on frequency, orientation, and contrast should preserve structure sensitivity better than standard coefficients in presence of strong frequency content. We checked this hypothesis, and found that best results were obtained from a 12-coefficient Gaussian distribution with adjustable variance σ 2 and anisotropy a.

Both σ and a are controlled by calibration. In the absence of local frequency content, standard coefficients H d are used. Transition from standard to structure-driven diffusion coefficients H b is controlled by linear interpolation with the calibration parameter w. Gaussian kernels are pre-calculated for fast run-time access (see Appendix A).

Calibration

Perhaps the most crucial step of our method is its calibration, the process by which we ensure that any local frequency information is associated with proper threshold modulation and error diffusion coefficients.

Our method relies on one fundamental hypothesis: by obtaining optimal halftoning results on an artificial image containing a single characteristic frequency, we should obtain near-optimal halftoning results on a region of a natural image that shows the same local frequency.

Automatic calibration was used by [START_REF] Ostromoukhov | A simple and efficient error-diffusion algorithm[END_REF]] and [START_REF] Zhou | Improving mid-tone quality of variable-coefficient error diffusion using threshold modulation[END_REF] for patches of uniform pixel luminance. Even though such an approach could be applied to our case, we chose to rely instead on the well-established method of side-byside comparisons used, for example, in the classical color matching experiment [START_REF] Wyszecki | Color Science: Concepts and Methods, Quantitative Data and Formulae[END_REF]. In our method, a uniform reference patch is shown side-by-side with a halftoned calibration patch whose parameters are adjusted by the observer until the best possible match between the two patches is reached. The experiment is repeated for several values of frequency/orientation/contrast of the reference patch.

Figure 4 illustrates the calibration process. To calibrate a given frequency, orientation, and contrast, we create a reference square patch with that frequency and orientation displaying two close but different contrasts: the central part, delimited by the red circle, shows the contrast to calibrate; the peripheral area shows a contrast value which has already been successfully calibrated.

First, the observer adjusts viewing distance until sine waves in the peripheral area are just imperceptible. At that distance, sine waves in the continuous-tone side of the central area should be at the threshold of visibility. Then, sliders are used to adjust the parameters for uncalibrated contrast until sine waves in the halftone side of the central area are also barely perceptible. The patches are dynamically updated, permitting the observer to act interactively.

Four parameters are required to control the behavior of our halftoning process: β is the strength of the applied threshold modulation, which is the amplitude of the Gabor convolution filter; σ and a are the parameters of the Gaussian error diffusion filter; w is a linear blending factor between no-contrast standard coefficients and noise-modulated threshold, and parametrized Gaussian coefficients and Gabor-filtered threshold. Once the calibration is done, a 3D lookup table is filled by interpolation of the calibration parameters for all frequency/orientation/contrast. They are used at run-time for a lookup table search of an arbitrary required frequency/orientation/contrast parameter.

All images shown in this paper use the same calibration. In our implementation, we use a total of around 1.3MB of lookup data split among five tables: interpolated calibration parameters for all frequency/orientation/contrast; Gabor filters for a variety of frequencies and orientations; Gaussian diffusion filters for a variety of σ and a; no-contrast standard diffusion filters for all tones; mid-tone improvement noise weights for all tones [START_REF] Zhou | Improving mid-tone quality of variable-coefficient error diffusion using threshold modulation[END_REF]. The pseudo-code in Appendix A indicates how each table is used.

Kernel size analysis

In Section 2.1, local frequency content is extracted from a pixel's neighborhood. The quality and computation time of the analysis depend on the size of the neighborhood window. To obtain the best quality/computation-time tradeoff, a systematic study was done on the choice of kernel size.

Table 1 shows how kernel size relates to computation time. Since local frequency analysis is parallel by nature, we implemented it both on GPU and CPU.

Figure 5 illustrates how results produced by our method can be degraded due to a poor frequency analysis. Some important structure details were lost with the 8 × 8 kernel, even more with the 4 × 4 one. Best results were observed for the 32×32 kernel, but the quality difference between the 16 × 16 and the 32 × 32 kernels is very small. The best tradeoff between speed and quality was reached by the 16 × 16 kernel on our hardware.

Results

We evaluate the results produced with our method using both subjective (visual) comparisons and objective measurements.

Visual comparison

For visual evaluation, we use the set of test images already shown in previous articles, namely in [START_REF] Pang | Structure-aware halftoning[END_REF]. This set comprises a wide range of important features: small details of various contrasts, smooth tone gradation, visually identifiable structures (leaves, hair, tissue or brush structure), etc.

Images in Figures 1 and7 show a clear and undeniable advantage of our method over standard error diffusion: small details are better preserved; our images appear sharper and globally nicer. Visual comparison of our method with structure-aware halftoning [START_REF] Pang | Structure-aware halftoning[END_REF] shows that both methods produce very comparable quality.

We added comparisons using a few artificial images: a grayscale ramp in Figure 6 and patches of sine waves of different frequencies and contrasts in Figure 2. The visual quality of the ramps is not a surprise: in absence of strong frequency content, our error diffusion becomes a mid-tone improved standard error diffusion. The patches in Figure 2 show a clear advantage of our method over all others, except the structure-aware halftoning. Compared to standard error diffusion, our method better reproduces the range of contrasts and frequencies. Our method performs much better than edge enhancement, too. Notice excellent visual coherence of the results produced with our method: the progression of relative contrast is visually close to that of the original.

Objective quality evaluation

Objective metrics were applied to all test images to confirm our method's capacity to preserve tones and structures.

The peak-to-signal noise ratio (PSNR) is a standard measure of reconstruction quality for lossy image compression. [START_REF] Pang | Structure-aware halftoning[END_REF]] used this metric on blurred versions of grayscale and halftone images to evaluate tone preservation. PSNR for Gaussian blurred (σ = 2) grayscale and halftone image pairs are given in Table 2. Results show that PSNR for our method remains consistently close to that of standard error diffusion and structure-aware halftoning.

The mean structure similarity measure (MSSIM) is a visual quality metric first introduced by [START_REF] Wang | Image quality assessment: from error visibility to structural similarity[END_REF]] and used by [START_REF] Pang | Structure-aware halftoning[END_REF]] as a central criterion for their convergence algorithm.

From the results shown in Table 3, we find that our method is close to structure-aware halftoning, which confirms visual observation.

Computation time

Table 4 shows computation time for our structure-aware error diffusion compared to previous methods. Although the quality of our results is close to that of [START_REF] Pang | Structure-aware halftoning[END_REF]], our method is dramatically faster. In fact, execution time under one second makes our method suitable for real-world applications. Note that our code is not optimal, especially its analysis part. With a relatively modest effort, our execution time could be further improved.

Limitations and future work

Although our method produces very satisfactory results in most cases, there are several issues which could be improved. First, our method works very well in presence of one clearly dominant local frequency. In rare cases when two or more equally important frequencies are present, our method is presently misguided by the analysis. We plan to carry out research which would incorporate treatment of multiple equally important frequencies in the framework described in this article. A better, specific analysis tool should be developed for this purpose.

In this paper, we limited our consideration to one-channel images. Extension to multi-channel (3-channel RGB or 4-channel CMYK images) may require additional work. A naïve straightforward extension such as mutually independent per-channel halftoning works rather well. But, in order to get even better quality for color images, a more sophisticated collaborative multi-channel error-diffusion scheme may be required. We plan to carry out research in this direction too.

The images shown in this paper have been prepared using a calibration with respect to a "dummy" LCD display. It would be quite worthwhile to integrate the calibration process into the framework defined by ICC [START_REF] Sharma | Digital Color Imaging Handbook[END_REF], to be used with a concrete printing device having strong hardware constraints (e.g. strong non-linear dot gain, depending on the nature of the printer: electrostatic, inkjet, etc.).

Conclusion

The error-diffusion method presented in this article combines the advantages of several halftoning methods known to date. On one hand, the overall visual quality of the images produced with our method is comparable with that of the best known error-diffusion algorithms, such as [START_REF] Ostromoukhov | A simple and efficient error-diffusion algorithm[END_REF][START_REF] Zhou | Improving mid-tone quality of variable-coefficient error diffusion using threshold modulation[END_REF]. On the other hand, our algorithm is sensitive to structures present in the source images, in a wide range of frequencies and contrasts. From this point of view, our algorithm approaches the quality of [START_REF] Pang | Structure-aware halftoning[END_REF], without being prohibitively slow. We have tested our algorithm with a wide variety of artificial and natural images, producing very satisfactory results. Our algorithm is conceptually simple and easy to implement.

Thus, we have shown that a structure-aware digital halftoning algorithm of Class II mentioned in Section 1 can be efficiently built.

We think that thanks to its unique combination of visual quality, sensitivity to structures, and speed, our algorithm may be used in a wide range of graphical applications in which digital halftoning is required. Typical applications may be printing, visualization, or various sampling techniques.

Figure 2 :

 2 Figure 2: Patches of sine waves. In each set of 6x6 patches, frequency changes along the x axis; contrast -along the y axis. Standard error diffusion spoils low-contrast patches, especially for high frequencies. Edge-enhancement is unable to reproduce faithfully all frequencies and contrasts. Please notice the excellent visual coherence of the results produced with our method: the progression of relative contrast is visually close to that of the original.

Figure 3 :

 3 Figure 3: Schematic representation of the architecture of our structure-aware error diffusion. Please refer to Sections 3.1 and 3.2 for the meaning of all terms.

Figure 4 :

 4 Figure 4: The calibration patches. Left side: continuous-tone sine waves of fixed frequency and angle. Right side: halftone images produced using adjustable set of parameters. The peripheral area delimited by the red circle serves to determine the appropriate observation distance. For the purpose of illustration, we exaggerated the contrast difference between inner and outer areas.

Figure 5 :Figure 6 :

 56 Figure 5: Structure-aware error diffusion using analysis tool of different kernel sizes.StandardMid-tone improved standard Structure-aware halftoning Our Method[START_REF] Ostromoukhov | A simple and efficient error-diffusion algorithm[END_REF][START_REF] Zhou | Improving mid-tone quality of variable-coefficient error diffusion using threshold modulation[END_REF][START_REF] Pang | Structure-aware halftoning[END_REF]

Table 2 :

 2 PSNR comparison.

Table 3 :

 3 MSSIM comparison.

	Our method	Our method	Structure-	Edge-	Std.
		+ Analysis	aware	enhanc.	error-
			halftoning		diff.
	0.22	0.97	135.41	0.10	0.04

Table 4 :

 4 Computation time (in seconds) for a typical 512 × 512 image

Acknowledgements

We would like to thank Franc ¸ois Duranleau, Neil Stewart, Jean Vaucher, and Pierre Poulin for their precious help. We are grateful to the authors of [Pang et al. 2008] for insightful discussions. The funding of this work was partially provided by the Natural Sciences and Engineering Research Council of Canada.

APPENDIX A: Pseudo-code of the Structure-Aware Error Diffusion STRUCTUREAWAREERRORDIFFUSION(Img, FreqContent) 1 Err ← NIL £ Array of quantization error values 2 for each processed pixel (i, j) 3 £ Retrieve threshold and diffusion filter from lookup tables 4 β, σ, a, w ← GETCALIBRATIONPARAMS(FreqContenti,j ) 5

Gabor ← GETGABORFILTER(FreqContenti,j )

Outi,j ← QUANTIZE(Imgi,j , modulatedThreshold ) 17

Erri,j ← Imgi,j -Outi,j 18 DISTRIBUTEERROR(i, j , Err , modifiedDiffusionFilter ) 19 return Out Lerp(v1, v2, t) is a linear interpolation function, which calculates (1 -t)v1 + tv2 [START_REF] Pharr | Physically Based Rendering: From Theory to Implementation[END_REF].