Luca Castelli Aleardi

Olivier Devillers

Jarek Rossignac

ESQ: Editable SQuad representation for triangle meshes

Keywords: triangle meshes, compact representations, dynamic data structures, maillages triangulaires, représentations compactes, structures de données dynamiques

We consider the problem of designing space efficient solutions for representing the connectivity information of manifold triangle meshes. Most mesh data structures are quite redundant, storing a large amount of information in order to efficiently support mesh traversal operators. Several compact data structures have been proposed to reduce storage cost while supporting constant-time mesh traversal. Some recent solutions are based on a global re-ordering approach, which allows to implicitly encode a map between vertices and faces. Unfortunately, these compact representations do not support efficient updates, because local connectivity changes (such as edge-contractions, edge-flips or vertex insertions) require re-ordering the entire mesh. Our main contribution is to propose a new way of designing compact data structures which can be dynamically maintained. In our solution, we push further the limits of the re-ordering approaches: the main novelty is to allow to re-order vertex data (such as vertex coordinates), and to exploit this vertex permutation to easily maintain the connectivity under local changes. We describe a new class of data structures, called Editable SQuad (ESQ), offering the same navigational and storage performance as previous works, while supporting local editing in amortized constant time. As far as we know, our solution provides the most compact dynamic data structure for triangle meshes. We propose a linear-time and linear-space construction algorithm, and provide worst-case bounds for storage and time cost.

Introduction

Many graphics, animation and modeling applications use triangle meshes. Different data structures have been proposed for representing them [START_REF] Boissonnat | Triangulations in CGAL[END_REF][START_REF] Baumgart | Winged-edge polyhedron representation[END_REF][START_REF] Bruce | A polyhedron representation for computer vision[END_REF][START_REF] Campagna | Direct edges -a scalable representation for triangle meshes[END_REF][START_REF] Guibas | Primitives for the manipulation of general subdivisions and computation of voronoi diagrams[END_REF]. In general one associates consecutive positive integer IDs with the different vertices, store the vertex data (locations, normals, textures) in an array, and use arrays to store incident relations between faces and vertices. In some cases the traversability (the ability of accessing neighboring cells at a constant time cost) is not required: for example, for applications such as rasterization and integral properties evaluation, the information directly available from the indexed face set suffices. However, many mesh processing functions (such as identifying silhouette edges, rendering with on-the-fly adaptive subdivision, or walking along the intersection with a plane or another mesh) require access to adjacent elements. Furthermore, some applications require a set of operators in order to locally update the mesh structure: this occurs, for example, in the incremental computation of Delaunay triangulation [START_REF] Lawson | Software for C 1 surface interpolation[END_REF][START_REF] Devillers | The Delaunay hierarchy[END_REF]. Others require changing the connectivity dynamically (such as adaptive mesh refinement [START_REF] Botsch | Polygon Mesh Processing[END_REF]). Usual mesh representations contain redundant information in order to achieve the prescribed requirements. Finally, the compactness has been considered more recently: the data structure should take as little storage as possible, so as to reduce page faults and cache misses. Although some applications store a large amount of information associated to mesh elements (such as vertex or faces), most applications only store vertex coordinates (often quantized to 16 bits or less) and possibly normal and texture coordinates (which can also be quantized). In general the cost of connectivity is really expensive and always dominates storage. For example, in the case of a triangle mesh, there are typically twice as many triangles as vertices: therefore the storage cost of connectivity in common data structures (such as Corner Table or Half-edge), ranges from ((6 × 2) + 1) × 32 = 416 to ((3 × 6) + 1) × 32 = 608 bits per vertex, while the storage of the quantized vertex locations is usually between 3 × 16 = 48 and 3 × 32 = 96 bits per vertex.

Traversable, modifiable and compact meshes

Two natural requirements for a mesh data structure are the accessability and navigability: traversable meshes should provide a basic set of operators allowing efficient (possibly constanttime) access and navigation between mesh elements. For example, the implementation of mesh traversals requires to perform the walk on the mesh around a face or to retrieve the neighbors of a vertex. Data structures should be also indexable. The application should be able to associate consecutive integer IDs to all triangles, and to access the IDs of their vertices in constant time, so as to support various marking and book-keeping functionalities. The reverse indexability may also be important to applications where one wishes to obtain the IDs of all triangles that are incident upon a given vertex. Modifiable data structures have the ability to dynamically maintain the mesh under local updates. Most applications need to perform local connectivity changes (attach/remove a triangle, insert or split a vertex, collapse or flip an edge) in constant time. A further requirement is the simplicity: a data structure should be simple to implement, so as to facilitate the development and testing of operators for mesh traversal and modification. The operators for changing the connectivity and for indexing or traversing the mesh should be fast in practice, so that their cost is negligible when compared to the geometric or photometric processing cost performed by typical applications. Because of the increasing complexity of surface meshes used in applications, a number of works propose compact data structures: the compactness of a mesh representation is a quality of measure, which can be defined as the average number of references stored per vertex. For dealing with meshes having huge size one has to eliminate, as much as possible, the redundancies of data structures: this is a crucial point, to reduce overall storage and hence memory thrashing.

(d) O(d) O(d) O(d) O(1) O(1) O(d) O(d) O(d) triangle split O(1) O(1) O(1) - - - - O(1) O(d 2 ) O(1) O(1) O(1) edge flip O(1) O(1) O(1) - - - - O(1) O(d 2 ) O(1) O(1) O(1) vertex delete O(1) O(1) O(1) - - - - O(1) O(d 2 ) O(d) O(1) O(d)

Prior art: mesh encodings and data structures

Mesh compression As already observed, representing a triangulation in a basic way induces a lot of redundancies. There exists several schemes allowing to efficiently compress the connectivity of a surface mesh into a few bits per vertex: this is the case, for example, of the Edgebreaker scheme [START_REF] Rossignac | Edgebreaker: Connectivity compression for triangle meshes[END_REF][START_REF] Lopes | Efficient edgebreaker for surfaces of arbitrary topology[END_REF], which encodes a planar triangulation of size n with at most 3.67n bits. A recent scheme [START_REF] Poulalhon | Optimal coding and sampling of triangulations[END_REF] is even more compact, allowing to compress into the optimal number of bits (3.24n), matching asymptotically Tutte's entropy bound [START_REF] Tutte | A census of planar triangulations[END_REF]. Although compressed formats reduce connectivity storage, they are not useful for mesh processing, since to support traversal operators they require full (or at least partial) decompression.

Classical data structures Most popular geometric data structures store a large number of references in order to describe incidence relations. Their classical implementations, in most programming environments are explicit representations: references allow to navigate in the data structure through address indirection. For example, in edge-based representations such as Halfedge [START_REF] Baumgart | Winged-edge polyhedron representation[END_REF], for each basic element (the half-edge) one stores 2 references to incident half-edges (the next one in ccw order in the same face, and the opposite one, lying in the neighboring face), plus a reference to a bounding vertex. Moreover, each vertex stores an incident half-edge, which leads to 3e + n references for a general manifold mesh having e edges and n vertices. Thus, in the case of triangle meshes, the Half-edge data structure uses 19 rpv (references per vertex ): other edge-based representations, such as Quad-edge [START_REF] Guibas | Primitives for the manipulation of general subdivisions and computation of voronoi diagrams[END_REF] or Winged-edge [START_REF] Bruce | A polyhedron representation for computer vision[END_REF], have similar storage requirements. Face-based representations (Corner Table, or triangle DS [START_REF] Boissonnat | Triangulations in CGAL[END_REF]) have slightly smaller requirements. They store, for each triangle, only 3 references to the neighboring faces (or opposite corners) and 3 references to the bounding vertices, plus a reference to an incident face for each vertex: this leads to a storage cost of (2 × (3 + 3)) + 1 = 13 rpv. Classical data structures are modifiable as they support updates in constant time: implementing a local modification (such as a triangle split or an edge flip) in the mesh involves a constant number of memory access and reference updates.

Compact data structures

The main goal of compact (practical) data structures ( [START_REF] Alumbaugh | Compact array-based mesh data structures[END_REF][START_REF] Castelli Aleardi | Catalog-based representation of 2D triangulations[END_REF][START_REF] Castelli Aleardi | Explicit array-based compact data structures for triangulations[END_REF][START_REF] Gurung | SOT: compact representation for tetrahedral meshes[END_REF][START_REF] Gurung | SQuad: Compact representation for triangle meshes[END_REF][START_REF] Gurung | LR: compact connectivity representation for triangle meshes[END_REF][START_REF] Lage | CHF: A scalable topological data structure for tetrahedral meshes[END_REF][START_REF] Snoeyink | Tripod: a minimalist data structure for embedded triangulations[END_REF]) is to reduce the number of references, while still guaranteeing an efficient implementation of navigational operations as in usual representations. This goal can be achieved,

Inria

for example, by grouping adjacent mesh elements (as neighboring triangles in [START_REF] Castelli Aleardi | Catalog-based representation of 2D triangulations[END_REF]): this allows to save some internal references, while still maintaining the representation under local modifications.

A more recent approach consists in re-ordering mesh elements (typically faces, or edges) according to a given permutation. Adopting this strategy, the SOT data structures [START_REF] Gurung | SOT: compact representation for tetrahedral meshes[END_REF] allows to efficiently navigate in the mesh, using at most 6 rpv. Combining this idea with a triangle pairing and a matching between triangles and vertices, it is possible to obtain a more compact version (SQUAD data structure [START_REF] Gurung | SQuad: Compact representation for triangle meshes[END_REF]) which uses about 4 rpv in practice (the LR data structure is even more compact, requiring about 2.16 rpv for tested meshes). Performing a re-ordering of input points, a recent edge-based representation [START_REF] Castelli Aleardi | Explicit array-based compact data structures for triangulations[END_REF] (referred to as sorted TRIPOD) guarantees the same navigation performances, providing a space bound of 4 rpv in the worst case (as in the TRIPOD data structure [START_REF] Snoeyink | Tripod: a minimalist data structure for embedded triangulations[END_REF], the main ingredient are Schnyder woods [START_REF] Schnyder | Embedding planar graphs on the grid[END_REF]).

Theoretical solutions: succinct representations For completeness, we also mention that succinct representations [START_REF] Castelli Aleardi | Succinct representations of planar maps[END_REF][START_REF] Castelli Aleardi | Succinct representation of triangulations with a boundary[END_REF][START_REF] Castelli Aleardi | Dynamic updates of succinct triangulations[END_REF] are successful in matching optimal asymptotic bounds for many class of meshes, running under the well known word-Ram model. In this case a memory word can store an arbitrary (small) number of references (or service bits) of tiny size: typically, one may stores up to O( lg n lg lg n ) sub-words of length l = O(lg lg n) each. Succinct data structures are of theoretical interest, mainly because the amount of auxiliary bits required, even if asymptotically negligible, remains important for triangulations of practical size. [START_REF] Kallmann | Star-vertices: a compact representation for planar meshes with adjacency information[END_REF] uses 7 rpv to the price of navigation in O(d) time, where d is the degree of the vertex involved. This structure is dynamic, as the mesh is editable: update operations can be performed in O(d 2 ) time, where d is the degree of involved vertices. The catalog-based representation described in [START_REF] Castelli Aleardi | Catalog-based representation of 2D triangulations[END_REF] has a storage cost slightly larger (the most compact version requires 7.67 rpv), while being fully dynamic: standard local modifications can be performed in O(1) time (as for [START_REF] Boissonnat | Triangulations in CGAL[END_REF], the number of reference updates is still constant). This catalog-based approach is somehow a practical version of [START_REF] Castelli Aleardi | Dynamic updates of succinct triangulations[END_REF], a succinct dynamic representation allowing to encode a triangulation having f faces with at most 2.17f +o(f ) bits, while supporting local modifications in O(log 2 f ) amortized time.

Modifiable compact solutions Star vertices

Our contributions

Most recent data structures for the connectivity of manifold triangle meshes are either compact, modifiable, or traversable, but not all three. Our contribution is to introduce new data structures achieving all the requirements above. The main idea is based on a new sorting strategy, making use of a matching from vertices to triangles: our novel approach for storing connectivity could lead others to invent new interesting dynamic mesh encoding schemes. We are able to describe a new class of mesh data structures supporting standard local edits in constant (amortized) time. As in prior works ( [START_REF] Castelli Aleardi | Explicit array-based compact data structures for triangulations[END_REF][START_REF] Gurung | SQuad: Compact representation for triangle meshes[END_REF][START_REF] Gurung | SOT: compact representation for tetrahedral meshes[END_REF][START_REF] Gurung | LR: compact connectivity representation for triangle meshes[END_REF]), local navigation operators are supported in worst case O(1) time, while the access to vertex data can require O(d) time (for a degree d vertex). Our data structures are simple to implement and provided with an analysis of worst case storage bounds. The simplest solution uses 6 rpv, while supporting updates operations in O(1) amortized time: this is obtained with a re-ordering of input data which allows to encode the map from triangles to vertices. Our most compact data structure makes use of a grouping strategy between adjacent triangles, and uses only 4.8 rpv, while still supporting efficient navigation and update operations. We provide experimental results concerning storage requirements and time-cost performances, obtained with implementations, which confirm the practical interest of our approach (experimental comparisons are reported in Table 2, while Table 1 shows the space requirements and worst-case performances for traversal and update operators of existing mesh representations).

v △ i i+1 i+2 g1 = neighbor(△, ccw(i)) g 0 g 1 g 2 g2 = neighbor(△, cw(i)) g0 = neighbor(△, i) △ = face(v) z z = vertex(g2, faceIndex(g2, △)) v = vertex(△, i) int degree(int v) { int d = 1; int f = face(v); int g = neighbor(f, cw(vertexIndex(v, f))); while (g ! = f) { i = vertexIndex(v, △) int cw(int i) {return (i + 2)%3; } int ccw(int i) {return (i + 1)%3; } int next = neighbor(g, cw(faceIndex(f, g))); int i = faceIndex(g, next); g = next; d + +; } return d; }

Preliminaries

Meshes and triangulations

Here we consider manifold triangle meshes, whose combinatorics correspond to simple triangulations embedded on a surface of arbitrary genus (without loops and multiple edges). Given a triangulation T , we denote by n (resp. by f ) the number of its vertices (resp. faces). Let us denote by V = {v 0 , v 1 , . . . , v n-1 } the set of its vertices, and by F = {△ 0 , △ 1 , . . . , △ f -1 } the set of its triangles.

Dynamic data structures

As in previous works [START_REF] Brodnik | Resizable arrays in optimal time and space[END_REF][START_REF] Castelli Aleardi | Dynamic updates of succinct triangulations[END_REF][START_REF] Raman | Succint dynamic dictionaries and trees[END_REF][START_REF] Castelli Aleardi | Catalog-based representation of 2D triangulations[END_REF] we assume that the programming environment provides a system memory manager which allows to allocate and free memory. The problem of designing and implementing efficient dynamic data structures is of both theoretical and practical interest. From the design and analysis point of view resizable arrays [START_REF] Brodnik | Resizable arrays in optimal time and space[END_REF] provide worst case O(1) time random access to data (reading and writing operations), while supporting updates (shrinking and growing operators) in amortized O(1) time, where additional storage cost is only O( √ n) (where n is the number of elements stored). Similar time-cost performances are also matched by practical implementations provided in classic programming languages (such as Java or C++), the wasted storage due to deallocation being of order Θ(n). In the rest of this work we assume to have an implementation of an abstract data type (ADT) for dynamic arrays, supporting the following operators with the performances above:

read(T, i): return the value stored at T [i], write(T, i, e): store element e at position i in T , grow(T ): increase the size of T , shrink(T ): decrease the size of T .

As usual, we assume elements have indices between 0 and n -1 and are all of equal length, coinciding with the size of memory words (32 or 64 bits in practice).

Inria

Abstract data type for dynamic meshes

Our representations are face-based and implement the same interface as in [START_REF] Castelli Aleardi | Catalog-based representation of 2D triangulations[END_REF] and [START_REF] Boissonnat | Triangulations in CGAL[END_REF] (the Corner Table data structure provides a slightly different interface). More precisely, we have the following operators: neighbor(△, i): retrieve the i-th neighbor of △, vertex(△, i): retrieve the vertex with index i of △, face(v): access one triangle incident to v.

With a combination of these operators it is also possible to define other functions, such as:

-faceIndex(△ 1 , △ 2 ): return the index of △ 1 among the neighbors of △ 2 . -vertexIndex(v, △): return the index of vertex v in face △.
The combination of the operators listed above allows to turn around a face in cw (or ccw) order, and thus to compute the degree of a vertex or to list its neighbors (as illustrated in Fig. 1). The mesh can be dynamically modified under the following update operators: defining a complete set of operators for editing a mesh, as in previous works: split(△, p): subdivide a given triangle △, into three new sub-triangles by inserting a new point p, flip(e): perform the flip of a given edge e, delete(v): remove a degree 3 vertex v, together with its three incident triangles.

3 Dynamic compact meshes: general scheme

Overview of our approach

As in previous compact representations, we exploit a re-ordering approach to reduce the number of references representing the map from vertices to faces and the map from faces to vertices. Here, we further explore the power of such re-ordering strategies. Instead of re-ordering mesh elements (typically faces [START_REF] Gurung | SQuad: Compact representation for triangle meshes[END_REF] or edges [START_REF] Castelli Aleardi | Explicit array-based compact data structures for triangulations[END_REF]), we rather exploit a permutation of geometric data associated with vertices. We avoid the indirection between triangles and vertices, by storing the vertex coordinates together with the data associated with one of its incident triangles: in this way we save 3 references for each triangle. The main novelty is to show that a coherent numbering of triangles and vertices can be maintained efficiently after local changes. Furthermore, as suggested in SQUAD [START_REF] Gurung | SQuad: Compact representation for triangle meshes[END_REF] and [START_REF] Castelli Aleardi | Catalog-based representation of 2D triangulations[END_REF], we also perform a grouping of few adjacent triangles in patches: this further reduces the number of references between neighboring elements.

Scheme description

Given a (triangular) mesh T having f faces, we consider a partition of its faces into patches. A patch is a face-connected triangulation with some marked vertices, namely a set of adjacent triangles whose dual graph is connected: in practice we will deal with patches which are simply connected and with a simple boundary. Given a face partition of T into patches, we say that such a partition is valid if the two following conditions are satisfied. Any face belongs to one and only one patch: patches can only share edges and vertices. A vertex v ∈ T may belong to several distinct patches, but is matched only once, and is marked in the corresponding patch.Two patches have the same shape if there is a bijection between their vertices that preserve triangles and marks. Two patches are of the same type (s, b, m), if they have the same number s of triangles, the same number m of marked vertices and boundaries of the same size b (b is the number of boundary edges). A catalog is the set of distinct (types of) patches defining the partition of the mesh: the size |C| of the catalog is the number of (different) patches it contains. ) and m j columns (one for each vertex marked in a patch stored in T j , if m j = 0 the table may be empty). The entries of P j are in correspondence with the patches stored in the associated table T j . More precisely, the a-th row P j [a] stores the coordinates corresponding to vertices w 1 , . . . w mj which are matched by the a-th patch in table T j . When required, additional informations associated to the faces can also be stored in T j .

Catalog-based representation

References encoding Recall that, according to our decomposition, each triangle belongs to exactly a patch, and each vertex is associated to a patch. Thus, in order to be able to refer to a given triangle △, a reference is composed of the following data: -a patch type: a code for indicating the type of patch p containing △, -a patch shape: a code for indicating the shape of patch p containing △, -a triangle index : a code for indicating the face △ among the triangles contained in p, -a patch index : the index of p in table T j (T j being the table containing the patches of same type as p).

A reference is encoded as an integer value r (typically on 32 bits): first bits (also referred to as service bits) are reserved for storing the three first codes, while remaining bits encode the patch index. Since in our representations the number of arrays, as well as the number of triangles per patch, are fixed, the first three codes require a constant number of bits (between 1 and 4 bits suffice in practice). The main part of a reference (representing the index of patch p, stored in T j ), requires ⌈log 2 |T j |⌉ bits. Retrieving service bits and patch indices can be performed with a combination of bit shift and bit mask operations (see Fig. 2). Observe that we do not need to store references to vertex data, which actually can be retrieved by the correspondence between vertices and matched patches.

Limitations As one can observe, our approach apply to meshes for which we can design an efficient matching correspondence (between vertices and patches) that can be easily maintained under dynamic updates. In the next sections we describe a few examples of data structures instantiating the general scheme described above, for the case of triangle meshes without boundaries (the case of boundaries is more involved and discussed in the last section). The choice of different catalogs of patches has a deep impact on the performances of mesh representations: for example, a small catalog (such as C 1 consisting of few types of patches) is simpler to implement and faster in practice compared to a bigger catalog (such as C 2 ) with better time-cost guarantees.

Inria T S u 0 u 1 u 2 u 3 u f -n-1
int neighbor(int r, int i) { if (patchType(r) == S) } return tableS[patchIndex(r) * 3 + i]; else return tableU[patchIndex(r) * 3 + i]; int vertex(int r, int i) { if (patchType(r) == S && i == 0) } return patchIndex(r); int f = neighbor(r, cw(i)); int j = faceIndex(r, f); while (f ! = r) { if (j == 1 && patchType(f) == S) return patchIndex(f); int next = neighbor(f, ccw(j)); j = faceIndex(f, next); f = next; } v f r next i j cw(i)

Instantiating the general scheme

The simplest catalog

The smallest catalog, denoted by C 1 , contains two types of patches, having only one triangle with zero or one matched vertex: we call U the patch with one triangle and no matched vertex, and S the triangle associated with one vertex (refer to Fig. 3).

Storage requirements

In such a case we have exactly two arrays of size n and fn respectively, one for patches of type U and one for patches of type S: with some abuse of notation, we denote them by T U and T S . Since the triangles in T U are unmatched, we need only one table P S to store vertex coordinates: P S has one column and n rows. P S [i] stores the coordinates of v i , while T S [i] stores the 3 neighbors of the triangle (denoted by s i ) matched with v i . As each patch stores 3 references to adjacent triangles, the cost for connectivity is 6 rpv (recall that there about 2n faces in a triangulation). Finally, observe that for this simple catalog we need only one service bit, for distinguishing between patches of type U and S: no need of more service bits, since all patches consist of exactly one triangle.

Matching triangles and vertices

The construction of our data structure relies on a matching between faces and vertices. This matching is dynamically maintained by the update operation described later, but if a triangulation is already given, then this matching can be computed as in [START_REF] Gurung | SOT: compact representation for tetrahedral meshes[END_REF] visiting the faces, according to an Edgebreaker traversal of the dual graph. The algorithm visits triangles (and vertices) in a depth-first manner (visiting right neighbor first), starting from a seed face (green face in Fig. 3): each time an unvisited vertex v is reached, it is matched to the current incident triangle. As for the Edgebreaker coder, the traversal ends when all faces and vertices have been visited, producing a valid partition (for more details see [START_REF] Gurung | SOT: compact representation for tetrahedral meshes[END_REF]).

Implementing navigational operators For the sake of completeness, we briefly sketch how to implement the navigational operators. Let us denote by r (resp. v) a reference to a given triangle (resp. vertex): recall that v ranges between 0 and n -1, while r is an integer whose first bit (service) describes the type of patch (the remaining bits encode the triangle index, which ranges between 0 and at most n). As one could expect, the implementation of neighbor operator is straightforward and fast, as it requires only one access in tables T U or T S . The implementation of vertex operator is more involved (as in [START_REF] Castelli Aleardi | Explicit array-based compact data structures for triangulations[END_REF][START_REF] Gurung | SOT: compact representation for tetrahedral meshes[END_REF][START_REF] Gurung | SQuad: Compact representation for triangle meshes[END_REF][START_REF] Gurung | LR: compact connectivity representation for triangle meshes[END_REF]): as we did not explicitly stored references . . . to vertices, we have to iteratively turn around a given vertex until the corresponding matched triangle is reached (the pseudocode is shown in Fig. 4).

(
Update operators and performance analysis Update operations for catalog C 1 are described in Fig. 5. Most update operations are quite easy to implement, involving a constant number of references updates in arrays T U and T S : the only exception concerns the case of vertex deletion.

The only difficult operation is the deletion of a degree 3 vertex, when all the three incident triangles have a marked vertex. In one triangle the marked vertex is the deleted one v, but after the deletion two vertices a and b are unmarked and we have a single triangle abc to assign it. The first workaround is found if one of the neighboring triangle has no marked vertex (Fig 5-bottomleft), then this triangle can be assigned one of the two orphan vertices. If the three neighboring triangles already have a marked vertex, only one of the triangles incident to edges bc and ca can have c as a marked vertex, assume without loss of generality that c is not marked in the triangle incident to bc; then a is marked in abc and among all the triangles incident to b at least one has to be free of mark (since c is marked in a triangle not incident to b) thus b can be assigned to that triangle (Fig 5 -bottom-right).

Although this last operation cannot be guaranteed to be in constant time and may require a time linear in the degree of b, it can be observed that this search for an unmarked triangle incident to b can be stopped as soon as such a triangle is found. Since half of the triangles are unmarked this exploration in average will not need to explore all the neighbors of b even if the degree of b is large.

A catalog with guaranteed constant time vertex removal

A slight modification of catalog C 1 allows to obtain a data structure that ensures constant time vertex removal, to the price to a less efficient edge flip (but still in constant time) (refer to sibgrapi/Figure 6). We add a new type patch to the catalog, called D, consisting of one triangle Performance analysis Since each patch contains one triangle, the storage analysis is the same as for catalog C 1 , leading to a cost of 6 rpv. Triangle split can be easily performed as before. The main differences concern the support of edge flip and vertex deletion.

Catalog C 3 U S Q 0 Q 1a Q 1b Q 2a Q 2b Q 2c Q 2d
Edge flips As for catalog C 1 , the correctness relies on an exhaustive case analysis. In addition to the cases already considered in Fig. 5, we have to deal with few more configurations, involving patches of type D. As depicted in Fig. 6, edge flips can always be supported with a constant number of memory changes. In most cases we need only some few read and write operations, which can be performed in worst case O(1) time on dynamic arrays. Unfortunately, there remains one case (refer to Fig. 6), where the edge flip implies the creation of two new patches of type S, and two removals (of patches D and U respectively): it still leads amortized O(1) time (the cost for updating tables T S , T D and T U ) but is clearly less efficient compared to catalog C 1 where the patches are just modified in place in the tables T U and T S .

Vertex removals The advantage of catalog C 2 is in the way we can perform vertex deletions more efficiently. As depicted in Fig. 6, we can use an adjacent triangle (say triangle (c, b, d)), in order to distribute the remaining vertices. For example, vertices a and b can be associated to triangle (a, b, c) of type D, while the remaining vertex c is matched by triangle (c, b, d): observe that such triangle is originally of type U or S, and becomes of type S or D after the removal of v. The removal/addition of a constant number of patches in tables T S , T U and T D can be supported in amortized O(1) time.

Triangles and quads: more compact catalogs

Here we consider catalogs containing patches which are simple triangles or quads (groups of paired adjacent triangles): as before they may have zero or more marked vertices.

Catalog C 3 : guaranteed upper bound The catalog C 3 (see Fig. 7) contains triangles with zero or one marked vertex and quads (groups of two triangles with zero or one marked vertex each). We need five catalogs T U , T S , T Q0 , T Q ' 1 , and T Q2 . For the type Q 1 and Q 2 there are several possible shapes depending on the position of the marked vertices with respect to the quad's diagonal and the triangle matched by the marked vertices. The gray pictures on Fig. 7 can be avoided yielding 16 different triangles on that pictures allowing 4 service bits to distinguish reports the experimental results obtained by TDS and ESQ representations on the tested 3D models: we compare both navigational and update operators. Our tests confirm the practical interest of ESQ approach: time-cost performances of ESQ are close to the ones of TDS, while the storage gain is significant (TDS is not compact, using 13rpv).

Navigational operators As one could expect, the ESQ implementation of the vertex operator is slower than TDS (about 12 times slower): recall that to retrieve a bounding vertex we have to turn around a vertex until the matched triangle is found (while TDS needs a memory access). Quite surprisingly, ESQ is (slightly) faster than TDS when comparing degree operator. The main reason is that ESQ can take advantage of the matching correspondence between triangles and faces. Given a vertex v we know that its index in s v (the matched triangle) is 0, so the clock-wise neighboring face around v is given by neighbor(s v , 2): which allows to save some computations compared to the TDS implementation.

Update operators When evaluating the split operators we only consider the combinatorial cost: we do not perform geometric calculations, and we split triangles just performing connectivity changes. As shown by our experiments, ESQ has (slightly) better performances than TDS: observe that both data structures perform the same reference updates to neighboring faces, while TDS has to update also vertex references (which do not exist in ESQ). Concerning the flip operator, ESQ is slower than TDS (about 1.7 times slower, in our tests). This comes from the number of cases to consider (as illustrated in Fig. 5): in order to distinguish all different cases ESQ has to perform a number of tests and bit-wise operations (to retrieve the patch type of the two neighboring faces).

Construction The preprocessing time for our ESQ data structure is also shown: it includes the matching phase, as well as the construction step (memory allocation and references setting). In all our tests, vertices and faces are accessed consecutively. We run our experiments on a Dell XT2, equipped with a Core 2 Duo 1.6GHz, Java 1.6 (the JVM using 1GB heap memory), under Windows 7 (32 bits).

Concluding remarks and extensions

We have proposed compact representations for triangle meshes supporting local editing operations. Our data structures are provided with guaranteed bounds (between 4.8 rpv and 6 rpv) and experimental evaluation. Further improvements and extensions are possible, according to the remarks below.

Inria

Meshes with boundaries Our results hold for triangle meshes which are manifold and without boundaries (and with arbitrary fixed genus). Dealing with boundaries is more involved and requires to design new matching rules (for assigning vertices to faces after local changes). Furthermore, our bounds for the triangle-quad catalogs should be updated to take into account boundaries of arbitrary size (relying on a slightly different counting argument).

Total vs. partial pairing In the static setting, one way to obtain better bounds (for the triangle-quad catalog) is to compute a perfect matching between pairs of adjacent triangles. Unfortunately, it is not clear how to maintain such a pairing when local changes are performed (without re-processing the entire mesh). A better trade-off between time-cost performances and lower memory occupancy could be achieved, in practice, with a smaller triangle-quad catalog C 4 (which would be simple to implement, involving only three types of patches, S, U single triangles and Q for quads). According to experimental tests (see [START_REF] Gurung | SQuad: Compact representation for triangle meshes[END_REF]), the matching and pairing strategy of the SQuad representation allows to regroup most of triangles into quads, leading to a cost of about 4 rpv.

Figure 1 :

 1 Figure1: Abstract data type for triangle meshes: our representations support the same operators as in[START_REF] Boissonnat | Triangulations in CGAL[END_REF][START_REF] Castelli Aleardi | Catalog-based representation of 2D triangulations[END_REF].

Figure 2 :

 2 Figure 2: References encoding.

Figure 4 :

 4 Figure 4: Implementation of navigational operators.

Figure 5 :

 5 Figure 5: Update operations for the catalog C 1 .

Inria

3 u 5 u 1 d 0 d 1 s 3 s 7 .

 17

Figure 6 :

 6 Figure 6: All updates operations can be perform in amortized O(1) time, with the slightly richer catalog C 2 .

Figure 7 :

 7 Figure 7: The triangle/quad catalog C 3 requires at most 4.8rpv.

Figure 9 :

 9 Figure 9: Left: triangle-quad catalog C 4 . Right: our tested meshes are manifold with no boundary (genus 0, 1 and 2).

Table 1 :

 1 Comparison between data structures for triangle meshes. Memory requirements are expressed in terms of references per vertex (rpv ) and hold in the worst case (at the exception of SQUAD and LR, whose performances come from experimental tests). Navigation and vertex access time hold in the worst case. The update cost is given in terms of the number of reference updates. We denote by d the degree of involved vertices.

  As in previous array-based representations[START_REF] Castelli Aleardi | Explicit array-based compact data structures for triangulations[END_REF][START_REF] Gurung | SOT: compact representation for tetrahedral meshes[END_REF][START_REF] Gurung | SQuad: Compact representation for triangle meshes[END_REF][START_REF] Gurung | LR: compact connectivity representation for triangle meshes[END_REF][START_REF] Alumbaugh | Compact array-based mesh data structures[END_REF], our data structure represents the mesh connectivity by storing a given number of references describing incident relations between mesh elements. Let us consider a partition of T based on a catalog C with patches of k different types. We represent the connectivity of T storing the references in k tables T 0 , T 1 , . . . T k-1 : table T i stores the connectivity of a patch with s i triangles and m i marked vertices. More precisely T i stores, for each patch of type (s i , b i , m

i ), a block of b i integers in consecutive memory words: they represent the b i references to faces in neighboring patches. Geometric information (and other information associated to vertices or triangles), are stored in a collection {P j } 0<j≤k of tables containing vertex coordinates, each associated to a table T j . P j has |T j | rows (one for each row in T j

  Construction of the ESQ data structure (catalog C 1 ). Triangles and vertices are matched according to a Edgebreaker traversal, starting from the seed face (green triangle). Patches of type S (blue faces) are stored in table T S : they are ordered according to the index of the matched vertex (table P S stores vertex coordinates). Table T U stores the unmatched patches (gray faces).

	Catalog C 1									s 0 s 1	s 5 s 12	s 2 u 2 s 4 . . .	(x 0 , y 0 , z 0 ) (x 1 , y 1 , z 1 )	s 3 s 14 s 11 s 6 s 8 . . .
	U	v S	11 6 s 6 s14	3 u 0 u 1 s11	8 s 3 s 8 s13 4	7 1 s 7 s 5 s 4	5 s 0 s 2 s 1 2	0 s10 u 2 s12	10 9 u 3	s 2 s 3 s 8 s 7 s 4 s 5 s 6	s 0 s 7 s 3 s 4 s 1 s 8 s 14	s 7 s 10 s 8 u 0 s 5 u 1 s 2 s 3 s 7 s 13 s 0 . . . u 1 . . .	(x 2 , y 2 , z 2 ) . . . . . . . . . . . . . . . . . .	. . . . . . s 0 s 10 . . . . . .
		P S				T U				s 9	. . .	u 2	. . .	. . .
										s n-1				. . .
	Figure 3:

RR n°8066

route des Lucioles -BP 93 06902 Sophia Antipolis Cedex Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399

This work is partially supported by ERC under the agreement ERC StG 208471 -ExploreMap.

all cases. Update operations are sketched in Fig. 8. The way the vertices are marked is not described and is identical to the one in Section 4.1.

Performance analysis If we enforce that two neighboring triangles are always merged in a quad, then we can lower bound the number of quads by 3 5 n and upper bound the number of isolated triangles by 4 5 n and get that we need only 4.8 rpv, as shown in [START_REF] Castelli Aleardi | Catalog-based representation of 2D triangulations[END_REF]. Table 2: Experimental results. We compare the time cost performances of our ESQ data structure (Catalog C 1 , using 6 rpv), with the triangle based data structure [START_REF] Boissonnat | Triangulations in CGAL[END_REF] (TDS), which is not compact and requires 13 rpv. Timings are expressed in terms of nanoseconds(ns) per operation. The preprocessing time (in seconds) for the construction of ESQ is also reported.

In order to evaluate practical performances, we have written Java array-based implementations of the triangle data structure [START_REF] Boissonnat | Triangulations in CGAL[END_REF] (TDS) and of our ESQ data structure (catalog C 1 ). Table 2


Sébastien Ménigot
email: sebastien.menigot@univ-tours.fr

Member, IEEE Jean-Marc Girault
email: jean-marc.girault@univ-tours.fr

Iulian Voicu

Member, IEEE Anthony Novell

Optimization of Contrast to Tissue Ratio by Frequency Adaptation in Pulse Inversion Imaging

Keywords: Adaptive system, closed loop system, contrast enhancement, microbubbles, optimal control, optimization, pulse inversion technique, signal processing, transmitted pulse, ultrasound imaging

Contrast imaging has significantly improved clinical diagnosis by increasing the contrast-to-tissue ratio after microbubble injection. Pulse inversion imaging is the most commonly used contrast imaging technique, as it greatly increases the contrast-to-tissue ratio by extracting microbubble nonlinearities. The main purpose of our study was to propose an automatic technique providing the best contrast-to-tissue ratio throughout the experiment. For reasons of simplicity, we proposed to maximize the contrast-to-tissue ratio with an appropriate choice of the transmit frequency.

The contrast-to-tissue ratio was maximized by a closed loop system including the pulse inversion technique. An algorithm based on the gradient provided iterative determination of the optimal transmit frequency. The optimization method converged quickly after six iterations. This optimal control method is easy to implement and it optimizes the contrast-to-tissue ratio by selecting the transmit frequency adaptively.

I. INTRODUCTION

O VER the past twenty years, improvements in the sensitivity of medical ultrasound imaging systems have provided more accurate medical diagnoses through intravenous injection of ultrasound contrast agents containing microbubbles. The perfusion imaging thus obtained, as in the myocardium or in tumors for example, has provided physiological and pathological information [START_REF] Frinking | Ultrasound Contrast Imaging: Current and New Potential Methods[END_REF]. Initially, the linear interactions between the microbubbles and the ultrasound waves were only operated in B-mode, to increase the sensitivity between the tissue and the microbubbles. However, the use of ultrasound contrast imaging was revolutionized in clinical practice when the nonlinear interaction was taken into account.

The nonlinearity of contrast agent responses has become a major focus of research to obtain the best contrast-to-tissue ratio (CTR). However, obtaining an ideal method has been limited by two factors. First, good separation of the harmonic components requires a limited pulse bandwidth [START_REF] Averkiou | Tissue Harmonic Imaging[END_REF], which reduces the axial resolution as in second harmonic imaging [START_REF] Burns | Instrumentation for Contrast Echocardiography[END_REF], and secondly the effects of the ultrasound wave propagation limit the CTR because of the presence of nonlinear components generated in tissue [START_REF] Frinking | Ultrasound Contrast Imaging: Current and New Potential Methods[END_REF].

Several imaging methods have been proposed to improve CTR and/or resolution. Some of the best known techniques such as second harmonic imaging [START_REF] Burns | Instrumentation for Contrast Echocardiography[END_REF], subharmonic imaging [START_REF] Forsberg | Subharmonic Imaging of Contrast Agents[END_REF], super harmonic imaging [START_REF] Bouakaz | Super Harmonic Imaging: A New Imaging Technique for Improved Contrast Detection[END_REF], imaging using second order Volterra filter [START_REF] Phukpattaranont | Post-Beamforming Second-Order Volterra Filter for Pulse-Echo Ultrasonic Imaging[END_REF] and attenuation correction [START_REF] Tang | Attenuation Correction in Ultrasound Contrast Agent Imaging: Elementary Theory and Preliminary Experimental Evaluation[END_REF] have been mainly based on post-processing. Some techniques have been based on post-processing with discrete or continuous encoding of the amplitude, the phase or the frequency of the transmitted ultrasound wave. They use differences in nonlinear acoustic signatures of microbubbles and tissue, such as pulse inversion [START_REF] Simpson | Pulse Inversion Doppler: A New Method for Detecting Nonlinear Echoes from Microbubble Contrast Agents[END_REF], power modulation [START_REF] Brock-Fisher | Means for Increasing Sensitivity in Non-linear Ultrasound Imaging Systems[END_REF], contrast pulse sequencing [START_REF] Phillips | Contrast-Agent Detection and Quantification[END_REF], pulse subtraction [START_REF] Borsboom | Pulse Subtraction Time Delay Imaging Method for Ultrasound Contrast Agent Detection[END_REF] and harmonic chirp imaging [START_REF] Borsboom | Harmonic Chirp Imaging Method for Ultrasound Contrast Agent[END_REF].

The setting parameters must be correctly adjusted for optimal use of these methods. However, there is usually no optimization process, the pulse settings to date being manual and pre-selected. The aim of the study presented here was to find the setting parameters of the technique to be used to provide the best CTR, the best resolution or the best compromise between CTR and resolution. This step is crucial.

Unfortunately, solution of the problem often requires a priori knowledge of the medium and transducer that is inaccessible. Moreover, the existing methods cannot adjust to different variations throughout a qualitative medical examination, as may be the case for a variation in microbubble concentration [START_REF] Becher | Handbook of Contrast Echocardiography : Left Ventrical Function and Myocardial Perfusion[END_REF] or microbubble size [START_REF] Soetanto | Study on the Lifetime and Attenuation Properties of Microbubbles Coated with Carboxylic Acid Salts[END_REF]. No method has been able to solve this problem to date.

In this study, we aimed to solve this problem through the concept of the optimal command: the system parameters θ were set by optimizing a cost function J (e.g. the CTR):

θ ⋆ = arg max θ (J (θ)) , (1)
where θ ⋆ are the optimal system parameters. We therefore replaced the current system with a closed loop system whose transmitted pulse was modified by feedback. Implementation of optimization required two steps: specification of the cost function J and specification of the setting parameters θ.

In the first step, the cost function J was chosen to take into account the user's needs and the medical application. Here, in contrast harmonic imaging, the cost function was the CTR. This choice was also justified by the microbubble sensitivity of the CTR. Moreover, to complement our approach, microbubble detection was performed by pulse inversion imaging [START_REF] Simpson | Pulse Inversion Doppler: A New Method for Detecting Nonlinear Echoes from Microbubble Contrast Agents[END_REF], since it is one of the most common methods used to increase CTR while ensuring good axial resolution.

In the second step, the CTR optimization must be achieved without a priori knowledge of the imaging system [START_REF] Ménigot | Improvement of the Power Response in Contrast Imaging with Transmit Frequency Optimization[END_REF] (i.e. explored medium, the transducer and pulse parameters), since such information is inaccessible in practice. As we wanted to optimize the CTR by modifying the signal transmitted, we could propose adjusting parameters characterizing the transmitted signal such as frequency, amplitude, phase or pulse duration. It should be noted that both the amplitude and the frequency are relevant pulse parameters which contribute to a significant increase in microbubble response, as in the case of power modulation or chirp imaging, respectively. To simplify our approach, we therefore chose to optimize only the transmit frequency f , for a pulse inversion imaging system, thereby producing a new device that can adapt to changing conditions.

Consequently, the optimization problem can be written as follows:

f ⋆ = arg max f (CTR(f )) , (2)
where f ⋆ is the optimal transmit frequency which provides the best CTR. We proposed an iterative approach to find the optimal transmit frequency f ⋆ and we applied it in simulations and in experiment.

II. CLOSED-LOOP IMAGING SYSTEM

The principle of pulse inversion imaging including feedback is described in Fig.

x k,p (t) = A • w k,p (t). (3)
The sinus modulated by a Gaussian function [START_REF] Borsboom | Harmonic Chirp Imaging Method for Ultrasound Contrast Agent[END_REF] w k,p (t) was constructed as follows:

w k,p (t) = exp - (t -t 0 ) 2 Nc 2f k sin(2πf k t + φ p ), (4)
where t is the time, t 0 the time for which the Gaussian function is maximum, N c the cycle number and φ p the phase equal to 0°if p = 1 and π if p = 2.

The amplitude of the driving pressure A was then adjusted so that the power of the pulse x k,p (t) was constant for all iterations:

A = A 2 0 • P x ref P w , (5)
where A 0 is the driving pressure amplitude of the reference signal x ref . This signal x ref was calculated at the central frequency of the transducer. Its power P x ref constituted the reference power, while P w was the power of the signal w k,p . The power of the transmitted wave thus remained constant by adjusting the amplitude signal A.

B. Cost Function

In the receiver, CTR k was computed from a line z k (t) of pulse inversion image:

z k (t) = y k,1 (t) + y k,2 (t), (6)
where y k,p (t) is the echo of the transmitted pulse x k,p (t). It is defined as the ratio of the power P b,k backscattered by the area of the perfused medium to the power P t,k backscattered by the area of the non-perfused medium [START_REF] Phukpattaranont | Post-Beamforming Second-Order Volterra Filter for Pulse-Echo Ultrasonic Imaging[END_REF] as follows:

CTR k = 10 • log 10 P b,k P t,k , (7)
These powers were measured from the lines z k (t) of the pulse inversion image at iteration k. The areas were determined manually before the optimization process, but a segmentation step could be implemented.

The contrast gain G dB was also defined between the optimized system and the non-optimized system.

The CTR obtained with the non-optimized system was determined at the central frequency of the transducer [START_REF] Ma | Improvement of Tissue Harmonic Imaging using the Pulse-Inversion Technique[END_REF]. The contrast gain G dB is obtained by the next equation:

G dB = CTR(f ⋆ ) CTR(f 0 ) , (8)
where CTR(f 0 ) is the CTR obtained at the central frequency of the transducer.

C. Iterative Optimization Algorithm

The algorithm was based on the principle of the gradient descent [START_REF] Widrow | Adaptive Signal Processing[END_REF]. It determined a new transmit frequency f k+1 for the next pulse sequence to optimize the CTR k+1 by the following recurrence relation:

f k+1 = f k + µ k • d k , (9)
The first coefficient µ k set the speed of convergence as follows:

µ k =                  0 if k 3; ∆f if k = 4; µ k-1 if sgn(∇CTR(f k )) = sgn(∇CTR(f k-1 )); µ k-1 2 if sgn(∇CTR(f k )) = sgn(∇CTR(f k-1 )). (10)
where ∆f fixed at 100 kHz provided the best compromise between speed of convergence and robustness, the sign function sgn(t) is equal to 1 if t > 0, 0 if t = 0 and -1 if t < 0, and the CTR gradient defined by:

∇CTR(f k ) = CTR k -CTR k-1 f k -f k-1 . (11)
The second coefficient d k set the direction as follows:

d k =          1 if k 3; 1 if sgn(∇CTR(f k )) = sgn(∇CTR(f k-1 )); -1 if sgn(∇CTR(f k )) = sgn(∇CTR(f k-1 )). (12)
In order to compute µ k and d k , the system operated in open loop for the first three iterations (k = {1, 2, 3}). The first three frequencies f 1 , f 2 and f 3 were chosen initially. The appropriate choice could increase the speed of convergence, but it was not essential to reach the optimal CTR, when the cost function was concave.

III. SIMULATION EVALUATION

The optimization principle was initially applied in simulation. Several simulations were performed to demonstrate the feasibility of our novel method.

A. Simulation Model

The simulation model followed the same process as described in the experimental setup (Fig. 1). It was composed of different phases: transmission, 2D nonlinear forward propagation, nonlinear oscillations of microbubbles, 2D nonlinear backward propagation and reception. A pulse wave was nonlinearly propagated into an attenuating intermediate medium without microbubbles. This wave, composed of harmonic components, excited a microbubble in the vascular system. The nonlinear oscillations of this microbubble were backscattered and measured by the receiver.

1) Nonlinear Propagation in Tissue:

The core of the model was based on the free simulation program developed by Anderson [START_REF] Anderson | A 2D Nonlinear Wave Propagation Solver Written in Open-Source MATLAB Code[END_REF], and e.g., used in [START_REF] Calle | Temporal Analysis of Tissue Displacement Induced by a Transient Ultrasound Radiation Force[END_REF]. It consisted of digitally solving the 2D nonlinear wave propagation into an attenuating medium by using a pseudo-spectral derivative and a time-domain integration algorithm.

A pulse signal was generated digitally at iteration k and filtered by the transfer function of the ultrasound probe used in the experiment, centred at 4 MHz with a fractional bandwidth of 53% at -3 dB. Here, to reduce the computation time, only 8 elements of the ultrasound probe composed of 128 elements were used. Note that linear beamforming was used in transmission and in reception [START_REF] Szabo | Diagnostic Ultrasound Imaging: Inside Out[END_REF] so that the beam was focused at 15 mm from the transducer. The nonlinear wave propagation into the intermediate medium was obtained by solving Anderson's model where the physical parameters were the following [START_REF] Szabo | Diagnostic Ultrasound Imaging: Inside Out[END_REF]:

• density of 928 kg•m -3 ;

• speed of sound of 1578 m•s -1 ;

• B/A nonlinearity parameter of 6.7;

• attenuation of 0.45

dB•MHz -1.05 •cm -1 .
This driving pressure at the focus point was included into the microbubble model described below. Finally, the wave backscattered by the microbubble was nonlinearly backpropagated up to the receiver. This signal was filtered by the transfer function of the same ultrasound probe used in experiment. To take into account imperfections in our simulation, a white noise ε(t) was added to x k,p (t), and a signal-to-noise ratio (SNR) of 50 dB was chosen.

2) Microbubble: The simulated ultrasound contrast agent had the properties of encapsulated microbubbles of SonoVue (Bracco Research SpA, Geneva, Switzerland), with a phospholipid monolayer imprisoning sulfur hexafluoride gas (SF 6 ) [START_REF] Greis | Technology Overview: SonoVue (Bracco, Milan)[END_REF] where the polytropic gas exponent κ was 1.095. The microbubbles had the following properties:

• mean diameter: 2.5 µm [START_REF] Greis | Technology Overview: SonoVue (Bracco, Milan)[END_REF];

• resonance frequency: 2.6 MHz [START_REF] Van Der Meer | The Resonance Frequency of SonoVue(TM) as Observed by High-SpeedOptical Imaging[END_REF].

The acoustic response was computed for one microbubble from the model of Marmottant [START_REF] Marmottant | A Model for Large Amplitude Oscillations of Coated Bubbles Accounting for Buckling and Rupture[END_REF] based on the Rayleigh-Plesset equation [START_REF] Plesset | The Dynamics of Cavitation Bubbles[END_REF] and the polytropic transformation. The solver used was the fourth-order Runge-Kutta method. Finally the surface pressure of the microbubble was transmitted to the nonlinear propagation model in order to deduce the backscattered signal. In order to simulate the mean behavior of a microbubble cloud, we hypothesized that the response of a cloud of N microbubbles was N times the response of a single microbubble with the mean properties. To be more realistic, the attenuation effects due to the high concentration of microbubbles were taken into account [START_REF] Gorce | Influence of bubble size distribution on the echogenicity of ultrasound contrast agents: A study of sonovue[END_REF] for a dilution of 1/2000.

B. Simulation Results

1) Empirical Optimization: Fig. 2 represents the results of the first simulation, i.e. the CTR as a function of the transmit frequency for different pressure levels A 0 (from 240 to 400 kPa). The number of cycles N c was set at 2.3 cycles so that the relative bandwidth 1 was 100% at the central frequency of the transducer.

The first results indicated that the CTR had a global maximum whatever the pressure level A 0 , and showed that the CTR can be improved by choosing the appropriate transmit frequency, thus confirming the validity of our study. This property could be also interesting, because an automatic search could more easily be achieved by a gradient algorithm. The second result demonstrated that the higher the pressure, the higher the global maximum of the CTR, because the power backscattered by the microbubble increased more quickly than the power backscattered by the tissue particle when the pressure level increased. The maximum values of the CTR ranged between 21.2 dB and 22.2 dB for pressure levels A 0 ranging from 240 to 400 kPa, and the contrast gain G dB was from 14.3 to 9.1 dB, respectively. These results revealed that the best transmit frequency was not the central frequency of the transducer, again confirming the need to optimize the imaging process.

To summarize, the results shown in Fig. 2 confirmed the need to optimize the imaging system by seeking the best transmit frequency which maximized the CTR function. Indeed, Fig. 2 shows that the CTR depended on transmit frequency f and could reach a single global maximum of approximately 22

dB. The single maximum of the CTR could be detected automatically by a simple technique such as the gradient. Note that the optimal transmit frequency was lower than the transmit frequency which produced the best axial resolution.

[Fig. 2 about here.]

2) Automatic Optimization: The maximum CTR was reached automatically using the gradient algorithm. Fig. 3b shows the CTR measured at each iteration k, and Fig. 3a shows variations in the transmit frequency during iterations.

The transmit frequency converged to a stable value after six iterations, whatever the pressure level A 0 . As an illustration, the black solid line in Fig. 2 shows the first twenty iterations which confirmed the convergence after the first six iterations. Moreover, the CTR reached its maximum when the transmit frequency converged. Note that the CTR and the contrast gain G dB obtained automatically were equal to those obtained empirically in the first simulation.

As depicted in Fig. 2, the optimal transmit frequency was lower than the central frequency of the transducer. Moreover, this optimal transmit frequency was equal to neither the central frequency of the 1 The relative bandwidth was defined as the percentage of the bandwidth of the signal in the fractional bandwidth of the transducer.

transducer nor to the microbubble resonance frequency. However, this optimal transmit frequency obtained after convergence of the algorithm made it possible to receive the harmonic components necessary to optimize the CTR.

In summary, the results in Fig. 3 show that it was possible to find the transmit frequency which maximized the CTR automatically. No a priori knowledge was required, except for the choice of the first three transmit frequencies which impacted on the speed of convergence.

[Fig. 3 about here.]

3) Adaptive Optimization: Fig. 4 depicts the CTR and the transmit frequency as a function of iterations when:

1) the microbubble radius was fixed at 2.5 µm throughout the simulation;

2) the microbubble radius was fixed at 1.25 µm throughout the simulation;

3) the microbubble radius decreased with time, increasing its linear resonance frequency. In the simulation, the radius changed from 2.5 to 1.25 at iteration 10.

The results showed that the first optimization (i.e. with a microbubble of 2.5 µm) converged at a transmit frequency of 1.9 MHz for a CTR of 22.1 dB. The second optimization converged at a transmit frequency of 2.1 MHz for a CTR of 16.6 dB. For the third optimization, the values of the CTR converged to the previously obtained values. Our system provided a quasi-instantaneous adaptation of the CTR during simulation. In all three simulations, convergence was reached after ten iterations. Note that the optimal transmit frequency was different for different sizes of microbubbles. Furthermore, when the oscillations increased as the size of the microbubble increased and thus the backscattered power increased. Hence, the CTR was dependent on the microbubble size.

To summarize, our optimization method adjusted the transmit frequency throughout the simulation by taking into account the changes in the microbubble properties. When the contrast agent did not change, the CTR did not change. This adaptability was also demonstrated when the tissue or the transducer changed with time. Note that when the optimization was not quick or robust enough, another algorithm could be used. However, for the gradient algorithm, the speed of convergence could enable us to reset the optimization parameter µ k if the environment changed a great deal.

[Fig. 4 about here.]

4) Note on Resolution: In the previous results, the axial resolution varied with the transmit frequency since the number of cycles N c was constant. An alternative solution can ensure a constant axial resolution if the variations in transmit frequency are compensated for by changing the number of cycles N c (k) for each iteration k. The length of the pulse x k,p (t) thus remained constant whatever the transmit frequency as follows:

2σ t (k) 2 = N c (k) 2f k = constant, ( 13
)
where 2σ t (k) 2 is the pulse duration at iteration k. In this case, when f k decreased, N c (k) decreased.

When the axial resolution was 100% of the relative bandwidth (i.e. 0.45 mm), the maximum values of CTR ranged between 17.9 dB and 19.8 dB for pressure levels A 0 increasing from 240 to 400 kPa, and the contrast gains G dB ranged between 10.9 and 6.6 dB. Note that the CTR and the gain were smaller with this setting compared to the previous results when was N c constant. Indeed, the small number of cycles reduced the good separation of harmonic components. Consequently, increasing the axial resolution reduced the CTR. Finally, the system converged at the maximum CTR with the same speed of convergence as indicated by Fig. 3.

IV. EXPERIMENTAL VALIDATION

The aim of this section was to confirm experimentally the results obtained in the simulation.

A. Experimental Setup

The experimental setup is presented in Fig. 1. The signals transmitted (x k,1 (t) and x k,2 (t)) were first generated digitally by a personal computer using equation 3. They were sent from an ultrasound scanner to the medium via an ultrasound probe. The waves insonified the medium which was composed of tissues and microbubbles. The reception system collected the echoes y k,1 (t) and y k,2 (t) and computed a pulse inversion image line.

1) Ultrasound Scanner and Transducers:

The transmitted signal x k,p (t) was sent to an "open" ultrasound scanner (MultiX WM, M2M, Les Ulis, France). This ultrasound scanner automatically duplicated the signal x k,p (t) for each element of the ultrasound probe. It applied the delays necessary to obtain phased-array beamforming [START_REF] Szabo | Diagnostic Ultrasound Imaging: Inside Out[END_REF]. Then the signals were transmitted to a linear array of 128 elements (Vermon SA, Tours, France), centred at 4 MHz with a fractional bandwidth of 53% at -3 dB.

The beam was focused at 28 mm from the surface.

The pulse transfer time for a single focusing was high, because it took approximately two seconds per radio frequency line. Consequently, we proposed optimization carried out only on an ultrasound image of five lines sweeping an angle of 0°25 ′ . The duration of the experiment was reduced to limit microbubble destruction.

2) Medium Explored: The wave propagated through a tissue-mimicking phantom (model 524, peripheral vascular Doppler flow phantom, ATS Laboratories Inc, Bridgeport, CT, USA), crossed by a 4 mm-diameter tube in which a 1/2000 diluted solution of SonoVue circulated.

The pulse was chosen with a cycle number corresponding to 55% of the relative bandwidth at the central frequency of the transducer (i.e. N c = 4) and with a pressure level A 0 of 400 kPa at the focal point.

B. Experimental Results

The experimental results presented in Fig. 5 show the transmit frequency and the CTR during the iterations. Fig. 5b shows that the CTR converged to its optimal value after six iterations. The mean CTR achieved after convergence was around 12.22 dB ±2.4 dB, i.e. a mean contrast gain of 7.77 dB. Fig. 5a shows that the transmit frequency converged to 2.9 MHz ±0.1 MHz. The temporal fluctuation of the transmit frequency showed that the system adjusted itself to changes. These variations may have been caused by two effects: (i) statistical fluctuation due to microbubble movements and (ii) fluctuation related to changes in the size distribution of microbubbles and number of microbubbles. However, these results did not permit definitive conclusion regarding the nature of the fluctuation and we could only conclude that our system automatically adjusted itself.

[Fig. 5 about here.] Furthermore, the experimental results were in accordance with our simulation results: we observed that the optimal transmit frequency was lower than the central frequency of the transducer. This optimal transmit frequency obtained after the algorithm convergence also made it possible to receive the harmonic components required in the CTR optimization. Note that the difference between the gain value in our simulation and that obtained in our experiment might be due to the more simplistic hypothesis of the model used.

Finally, these experimental results confirmed the feasibility of our method.

V. DISCUSSION AND CONCLUSIONS CTR optimization in pulse inversion imaging was performed automatically, without taking into account a priori knowledge of the medium or the transducer, except for the first three values of the transmit frequency chosen only for their impact on the speed of convergence. The optimal transmit frequency was reached automatically by feedback within only a few iterations. The algorithm itself adjusted the transmit frequency throughout the experiment by each time taking into account the changes in microbubble properties. This closed loop system gave the best compromise between the transducer bandwidth and the frequency responses of microbubbles and tissue. This trade-off has usually been calculated empirically to date. To make such automatic trade-off possible, the proposed algorithm itself adjusted the transmit frequency to maximize the power backscattered by microbubbles while minimizing the power backscattered by the tissue within the transducer bandwidth.

Our method was feasible for two reasons. Optimization was iteratively fulfilled by using (i) an easily implemented algorithm and (ii) a single setting parameter (here the transmit frequency). One major advantage of our approach was that it adjusted itself to any medium explored by taking into account the effects of attenuation and nonlinear propagation. This advantage was due to the fact that the cost function, on which the optimization was based, was exclusively the result of the input and the output measurements of our system. One interesting consequence is that our method can be applied to any imaging system and to any medium to be explored, since our algorithm should converge to the maximum of the cost-function. Note that for a robust optimization, this maximum must be exclusive.

We identified three major points for discussion, for future integration in an imaging system:

• the first concerned the gradient algorithm for which six basic operations were necessary to compute the transmit frequency for the next iteration. This low number of operations should not significantly change the frame rate;

• the second concerned the CTR computation from regions of interest (L × L size) in the image of size M ×M. For efficient optimization, it is important to determine the correct positions of the perfused and non-perfused areas. Indeed, poor assessment of the CTR could lead to poor optimization. Moreover, the CTR computation required 2(2L + 1) 2 + 1 operations. It should not considerably change the frame rate;

• the third concerned the transfer time between the computer and the programmable analogue transmitter. Note that for our "open" ultrasound scanner, this transfer time was not negligible compared to the frame rate. Indeed, this had led us to limit the image size. However we do not think that is is a problem for future integration of our method in an imaging system, since the current development of new imaging methods based on chirp or time reversal also requires such instrumentation.

Finally, our method only focused on qualitative imaging. Although our technique could offer an optimal frequency for each image line, it was preferable to perform optimization on the whole image, because the image had a single resolution. For instance, our optimization without resolution constraints could offer the best tradeoff between the CTR and axial resolution. We therefore believe that the method would be particularly appropriate for contrast echocardiography where the tradeoff must favor contrast [START_REF] Burns | Instrumentation for Contrast Echocardiography[END_REF].

Nevertheless, our method was suitable for quantitative contrast ultrasound imaging, since it adjusted itself to microbubble variations, although a possible solution may be to cease optimization during the quantification step.

To conclude, the method ensured optimal CTR throughout the experiments by adaptively selecting the transmit frequency. With our new approach, manufacturers and clinicians would not themselves need to tune the transmit frequency. The method should automatically adapt the transmit frequency to the medical examination conditions and maintain optimalCTR. Finally, our closed-loop method should be adapted using a larger number of contrast imaging techniques.

Fig. 1 :

 1 Fig. 1: Block diagram of adaptive imaging by pulse inversion technique.

Fig. 2 :

 2 Fig.2: Simulation of CTR for a transmit frequency of 1 to 4 MHz, a pressure level A 0 of 240 to 400 kPa and a constant number of cycles N c of 2.3. The simulation included the transducer. Note that for transmit frequencies above the central frequency of the transducer, the axial resolution could not be reduced because of the transducer impulse response.

Fig. 3 : 2 . 3 .

 323 Fig. 3: Simulation of automatic optimization of the CTR by iterative searching for the optimal transmit frequency for different pressure levels A 0 from 240 to 400 kPa and a constant cycle number N c of 2.3. The simulation included the transducer.

Fig. 4 :

 4 Fig.4: Simulation of CTR optimization for a pressure level A 0 of 400 kPa and a constant cycle number N c of 2.3 in the presence of changing microbubble properties. The first optimization was performed for a microbubble with a constant radius of 2.5 µm; the second for a microbubble with a constant radius of 1.25 µm. The third optimization began with a microbubble of 2.5 µm radius until iteration 10, the microbubble then changed to 1.25 µm radius.

Fig. 5 :

 5 Fig. 5: Experiment of CTR optimization for a pressure level A 0 of 400 kPa at the focal length and a constant cycle number N c of 4.

  (t) of the two respective echoes y k,1 (t) and y k,2 (t) formed a radiofrequency line of the image I k . From the CTR k measured on this image I k , a new transmit frequency f k+1 was computed by the algorithm to optimize the CTR k+1 on the next image I k+1 .

	[Fig. 1 about here.]
	A. Transmitted Signal
	The pulse signal x k,p (t) at transmit frequency f k was computed digitally with Matlab (Mathworks,
	Natick, MA, USA):

1. At the iteration k, two pulses x k,1 (t) and x k,2 (t) with opposite phases and with a frequency f k were transmitted. The sum z k

ACKNOWLEDGEMENTS

The authors would like to thank A. Zaylaa and A. Bouakaz for helpful discussions and for language editing of the manuscript, and the Clinical Investigation Centre for Innovative Technology of Tours (CIC-IT 806 CHRU de Tours, Tours, France) for the ultrasound contrast agents.

The authors would also like to thank the anonymous reviewers for their valuable comments.

The authors thank the Agence Nationale de la Recherche (Projects ANR-07-TECSAN-015 MONITHER and ANR-07-TECSAN-023 SURFOETUS) for financial support.

The authors are with Université Franc ¸ois-Rabelais de Tours, UMR-S930, Tours, France. They are also with Inserm,

Generator and Beamforming

Transducer Tx


Sébastien Ménigot
email: sebastien.menigot@univ-tours.fr

Jean-Marc Girault
email: jean-marc.girault@univ-tours.fr

Transmit Frequency Adaptation for Contrast to Tissue Ratio Optimization in Second Harmonic Imaging

Keywords: Adaptive system, closed loop system, contrast enhancement, microbubbles, optimal control, optimization, second harmonic, signal processing, transmitted pulse, ultrasound imaging

Conventional ultrasound contrast imaging systems use sequences of sine waves at a fixed transmit frequency. The excitation of the system is thus a sine wave whose frequency is often fixed around two-thirds of the center frequency of the transducer in the second harmonic imaging. However, this choice requires a knowledge of the transducer and does not take into account the specificities of microbubbles and medium. The main purpose of our study was to propose an automatic technique providing the best contrast-to-tissue ratio by selecting the optimal transmit frequency. The simulations and the experiment showed that the optimal transmit frequency value did not correspond to the usual value. By using the optimized value, the contrast was improved by around 6 dB. Finally, by providing a closed loop system, the system automatically proposed the optimal control without any a priori knowledge of the system or the medium explored.

I. INTRODUCTION

Over the past twenty years, research in medical ultrasound imaging systems have enabled to improve the image quality in terms of sensitivity, contrast and resolution. To provide additional qualitative information about physiology and pathology [START_REF] Frinking | Ultrasound Contrast Imaging: Current and New Potential Methods[END_REF], clinicians have used intravenous injection of ultrasound contrast agents containing micro-bubbles. These microbubbles can generate harmonic components of the transmitted ultrasound wave explained from the equations of nonlinear dynamics of the microbubble [START_REF] Leighton | The Acoustic Bubble[END_REF]. Moreover, these harmonic components generated by the microbubbles are more intense than those generated by soft tissues when the pressure level of the transmitted wave is small [START_REF] Borsboom | Nonlinear Coded Excitation Method for Ultrasound Contrast Imaging[END_REF]. The contrast is thus enhanced by increasing the contrast-to-tissue ratio (CTR) and particularly the harmonic power P b backscattered by the microbubbles compared to the harmonic power P t backscattered by the tissue [START_REF] Phukpattaranont | Post-Beamforming Second-Order Volterra Filter for Pulse-Echo Ultrasonic Imaging[END_REF]:

CTR (2f ) = P b (2f ) P t (2f ) (1)
where f is the transmit central frequency. However, obtaining an ideal imaging method has been limited by two factors. Firstly, a good separation of the harmonic components requires a limited pulse bandwidth [START_REF] Averkiou | Tissue Harmonic Imaging[END_REF], which reduces the axial resolution as in the second harmonic imaging [START_REF] Burns | Instrumentation for Contrast Echocardiography[END_REF]. Secondly, the effects of the ultrasound wave propagation limit the CTR. Several imaging methods have been proposed to improve the contrast. Some techniques are based on post-processings such as second harmonic imaging [START_REF] Burns | Instrumentation for Contrast Echocardiography[END_REF], subharmonic imaging [START_REF] Forsberg | Subharmonic Imaging of Contrast Agents[END_REF], super harmonic imaging [START_REF] Bouakaz | Super Harmonic Imaging: A New Imaging Technique for Improved Contrast Detection[END_REF], imaging using second order Volterra filter [START_REF] Phukpattaranont | Post-Beamforming Second-Order Volterra Filter for Pulse-Echo Ultrasonic Imaging[END_REF] and attenuation correction [START_REF] Tang | Attenuation Correction in Ultrasound Contrast Agent Imaging: Elementary Theory and Preliminary Experimental Evaluation[END_REF]. Other techniques are based on discrete or continuous encoding of amplitude, phase or frequency of the ultrasound wave transmitted such as pulse inversion [START_REF] Simpson | Pulse Inversion Doppler: A New Method for Detecting Nonlinear Echoes from Microbubble Contrast Agents[END_REF], power modulation [START_REF] Brock-Fisher | Means for Increasing Sensitivity in Non-linear Ultrasound Imaging Systems[END_REF], contrast pulse sequencing [START_REF] Phillips | Contrast-Agent Detection and Quantification[END_REF], pulse subtraction [START_REF] Borsboom | Pulse Subtraction Time Delay Imaging Method for Ultrasound Contrast Agent Detection[END_REF] and har-monic chirp imaging [START_REF] Borsboom | Nonlinear Coded Excitation Method for Ultrasound Contrast Imaging[END_REF], [START_REF] Borsboom | Harmonic Chirp Imaging Method for Ultrasound Contrast Agent[END_REF]. Since the easiest and the most commonly used technique is the second harmonic imaging, we focused our study on this technique. The optimal use of such technique is only ensured if the setting parameters, e.g. the transmit frequency, is adjusted to its optimal value. However, this optimal setting is not assured in commercial devices to date, since the solution of our problem often requires inaccessible a priori knowledge of the medium and the transducer. Consequently, the transmit frequency was only set to the two-thirds of the transducer central frequency [START_REF] Hossack | A High Bandwidth Transducer Optimized for Harmonic Imaging[END_REF] from empirical inference. Thus the question that we can naturally ask is: is there a choice of a transmit frequency which maximizes the criterion of contrast as the CTR?

In this study, we aimed to solve this problem through the concept of the optimal command: the transmit frequency f was set by optimizing CTR in second harmonic imaging:

f ⋆ = arg max f (CTR(f )) , (2)
where f ⋆ is the optimal transmit frequency which provides the best CTR. Therefore we replaced the current system with a closed loop system whose transmitted pulse was modified by feedback. We proposed an iterative approach to find the optimal transmit frequency f ⋆ and we tested the originality of our solution both by simulation and experiment.

II. OPEN LOOP IMAGING SYSTEM As previously mentioned, the studied imaging system was a second harmonic imaging system [START_REF] Burns | Instrumentation for Contrast Echocardiography[END_REF]. The principle of this system was described in Fig. 1. It consists of transmitting an incident ultrasound wave of frequency f to the medium, and then receiving the harmonic component generated by the medium at the frequency 2f . The contrast is thus enhanced by improving the microbubble detection compared to the tissue. Note that the bandwidth of the transducer was shared by the two filters (emission and reception filters).

The pulse signal x k (t) at transmit frequency f k was computed digitally with Matlab (Mathworks, Natick, MA, USA):

x k (t) = A • w k (t). (3)
Note that in the case of the open loop system, the iteration was always the same with k = 1. The Gaussian-modulated sinusoidal pulse [START_REF] Borsboom | Harmonic Chirp Imaging Method for Ultrasound Contrast Agent[END_REF] w k (t) was constructed as follows:

w k (t) = exp - (t -t 0 ) 2 Nc 2f k sin(2πf k ), ( 4
)
where t is the time, t 0 the time for which the Gaussian function is maximum and N c the cycle number. The value of N c must be significant enough in order to limit the transmitted bandwidth. In this case, direct transmissions at the harmonic frequencies could be avoided and thus measured the harmonic components provided by a nonlinear behavior [START_REF] Averkiou | Tissue Harmonic Imaging[END_REF].

Finally, the power P w of the transmitted signal was equal to the power P x ref of the reference signal defined to the two-thirds of the central frequency 2/3f c of the transducer and with a driving pressure A 0 . The amplitude of the driving pressure A was thus adjusted as follows:

A = A 2 0 • P x ref P w . (5)
At the receiver, the backscattered signal y k (t) was filtered by the transducer and by the reception filter at the frequency 2f . The final signal z k (t) represented a line of the second harmonic image.

CTR k was then computed from the line z k (t). However, in order to compute the powers P b and P t in eq. 1, the area of the microbubbles and the area of the tissue were determined manually before the optimization process.

III. CLOSED LOOP IMAGING SYSTEM

The closed loop system was composed of the open loop system and a feedback. This feedback ensures the optimization of the cost function, which is here the CTR in the second harmonic imaging (eq. 2). This cost function was mainly based on the upcoming hypotheses. It must:

• depend on the transmit frequency f ;

• be, in absolute, independent of the simulation model or the experiment;

• be concave in the range of the transmit frequency f to increase the robustness of the algorithm.

A. Iterative Optimization Algorithm

Based on the last hypotheses, a simple algorithm can be based on the principle of the gradient descent [START_REF] Widrow | Adaptive Signal Processing. Englewood Cliffs[END_REF]. It determined a new transmit frequency f k+1 at the iteration k+1 for the next pulse sequence to optimize the CTR k+1 according to the following recurrence relation:

f k+1 = f k + µ k • d k , (6)
The first coefficient µ k set the speed of convergence as follows:

µ k =          0 if k 3; ∆f if k = 4; µ k-1 if sgn(∇CTR(f k )) = sgn(∇CTR(f k-1 )); µ k-1 2 if sgn(∇CTR(f k )) = sgn(∇CTR(f k-1 )). (7)
where ∆f was fixed at 100 kHz providing the best compromise between the speed of the convergence and the robustness, the sign function sgn(t) is equal to 1 if t > 0, 0 if t = 0 and -1 if t < 0, and the CTR gradient defined by:

∇CTR(f k ) = CTR k -CTR k-1 f k -f k-1 . ( 8
)
where f ⋆ is the optimal transmit frequency. The second coefficient d k set the direction as follows:

d k =    1 if k 3; 1 if sgn(∇CTR(f k )) = sgn(∇CTR(f k-1 )); -1 if sgn(∇CTR(f k )) = sgn(∇CTR(f k-1 )). (9)
In order to compute µ k and d k , the system operated in open loop for the first three iterations (k = {1, 2, 3}). The first three frequencies f 1 , f 2 and f 3 were chosen initially. The appropriate choice could increase the speed of convergence, but it was not essential to reach the optimal CTR, when the cost function was concave.

B. Gain

Finally, the contrast gain G dB was defined between the optimized system and the non-optimized system. The CTR obtained with the non-optimized system was determined at the two-thirds of the central frequency 2/3f c of the transducer [START_REF] Hossack | A High Bandwidth Transducer Optimized for Harmonic Imaging[END_REF]. The contrast gain G dB was obtained by the following equation:

G dB = CTR(f ⋆ ) CTR(2/3f c ) . (10)
IV. SIMULATION EVALUATION

The optimization principle was initially applied in simulation. Several simulations were performed to demonstrate the feasibility of our new method.

A. Simulation Model

The simulation model followed the same process as in the experimental setup (Fig. 1). A pulse signal was digitally generated at the iteration k and filtered by the emission filter at 2.25 MHz with a fractional bandwidth of 74% at -3 dB. It was then transmitted simultaneously in microbubble and tissue models to deduce the surface pressures. In reception, the signal measured by the transducer was filtered around 3.5 MHz with a fractional bandwidth of 63% at -3 dB. Note that the transducer was centred at f c = 3 MHz and with a fraction bandwidth of 160% at -3 dB.

1) Microbubbles:

The simulated ultrasound contrast agent had the properties of encapsulated microbubbles of SonoVue (Bracco Research SpA, Geneva, Switzerland). The microbubble gas was composed of a phospholipid monolayer imprisoning sulfur hexafluoride gas (SF 6 ) [START_REF] Greis | Technology Overview: SonoVue (Bracco, Milan)[END_REF]. The mechanical shell properties were the shear modulus G s 46 MPa [START_REF] Vos | Method for Microbubble Characterization Using Primary Radiation Force[END_REF] and the resonance frequency of the microbubble was 2.1 MHz [START_REF] Ménigot | Improvement of the Power Response in Contrast Imaging with Transmit Frequency Optimization[END_REF].

To carry out the simulations, the free simulation program BUBBLESIM by Hoff [START_REF] Hoff | Acoustic Characterization of Contrast Agents for Medical Ultrasound Imaging[END_REF] was used to calculate the oscillations and the scattered echoes for a specified contrast agent and excitation pulse. A modified version of the Rayleigh-Plesset model was chosen. The model was based on the theoretical description of microbubbles as air-filled particles with surface layers of elastic solids presented by Church [START_REF] Church | The Effects of an Elastic Solid Surface Layer on the Radial Pulsations of Gas Bubbles[END_REF] and later modified by Hoff [START_REF] Hoff | Acoustic Characterization of Contrast Agents for Medical Ultrasound Imaging[END_REF]. The properties of the surrounding medium were those of blood, since the microbubble is supposed to be moving in the vascular system. In order to simulate the mean behavior of a microbubble cloud, we hypothesized that the response of a cloud of N microbubbles was N times the response of a single microbubble with the mean properties.

2) Tissue:

The tissue response was simulated by fat globules with a density of 928 kg/m 3 [START_REF] Szabo | Diagnostic Ultrasound Imaging: Inside Out[END_REF]. The computation of their response was based on Rayleigh backscattering [START_REF] Rayleigh | The Theory of Sound[END_REF] for a small fat ball of 10 µm, the size that was chosen to approximate the small size of fat cells. We hypothesized that the response of N particles was N times the response of a single particle.

B. Simulation Results

1) Empirical Optimization: Fig. 2 represents the results of the first simulation, i.e. the CTR as a function of the transmit frequency for different pressure levels A 0 (240 and 370 kPa) and for 4 cycles.

The first result indicated that the CTR had a global maximum whatever the pressure level A 0 . These maxima were obtained neither with the twothirds of the central frequency of the transducer, nor with the microbubble resonance frequency. This result showed that the CTR can be improved by choosing the appropriate transmit frequency, thus confirming the validity of our study. Secondly, the global maximum was the single maximum in the studied range. This property was interesting, because an automatic search could be achieved more easily by a gradient algorithm. Thirdly, the maximum values of the CTR ranged between 17.4 dB and 20.5 dB for pressure levels A 0 240 and 370 kPa, and the contrast gain G dB ranged from 5 to 5.9 dB.

To summarize, the results shown in Fig. 2 confirmed the need to optimize the imaging system by seeking the best transmit frequency which maximized the CTR function. Indeed, Fig. 2 shows that the CTR depended on transmit frequency f and could reach a single global maximum of approximately 20 dB. The single maximum of the CTR could be detected automatically by a simple technique such as the gradient.

2) Automatic Optimization:

The maximum CTR was reached automatically using the gradient algorithm. Fig. 3b shows the CTR measured at each iteration k, and Fig. 3a shows variations in the transmit frequency during iterations. The transmit frequency converged to a stable value after six iterations, whatever the pressure level A 0 . As an illustration, the black solid line in Fig. 2 shows the first twenty iterations which confirmed the convergence after the first six iterations. Moreover, the CTR reached its maximum when the transmit frequency converged. Note that the CTR and the contrast gain G dB obtained automatically were equal to those obtained empirically in the first simulation.

In summary, the results in Fig. 3 showed that it was possible to find the transmit frequency which maximized the CTR automatically. No a priori knowledge was required, except for the choice of the first three transmit frequencies which impacted on the speed of convergence.

V. EXPERIMENTAL VALIDATION

The aim of this section was to confirm experimentally the results obtained in simulation.

A. Experimental Setup

The experimental setup was presented in Fig. 1. The transmitted signal was firstly generated digitally by a personal computer using equation 3. It was sent to the medium via a transducer. The wave insonified the medium which was composed of tissues and microbubbles. The receiver collected the echo y k (t) and computed a second harmonic image line.

1) Measurements and Transducers:

The pulse was chosen with a cycle number N c of 4 and with a pressure level A 0 of 370 kPa at the focal point. This transmitted signal x k (t) was transmitted through a GPIB port (National Instruments, Austin, TX) to an arbitrary function generator (33220A, Agilent, Palo Alto, CA). The signal was then amplified using a power amplifier (Amplifier Research 150A100B, Souderton, PA) and transmitted to a 2.25 MHz PZT single element transducer (Sofranel, Sartrouville, France) focused at 55 mm and with a fractional bandwidth of 74%.

The backscattered signals were measured by a 3.5 MHz PZT single element transducer with a fractional bandwidth of 63%, placed perpendicularly to the transmitter, also focused at 55 mm and used in the receiver mode. Note that this setup configuration can look upon as a single transducer centred at f c = 3 MHz and with a bandwidth of 160% at -3 dB.

The echoes measured were amplified by 30 dB (Pulser-Receivers 5072R, Olympus Panametrics, Waltham, MA, USA) and then visualized on a digital oscilloscope (Tektronix, Beaverton, OR). Signals were transferred to a personal computer through a GPIB port for further analysis.

2) Medium Explored:

The wave propagated through a blood mimicking fluid ("EU-DFS-BMF-ver.1.2", Dansk Fantom Ser-vice, Jyllinge, Denmark). This latter was crossed by a 1cm-side tank in which there was a 1/2000 diluted solution of Sonovue™ microbubbles.

B. Experimental Results

The experimental results presented in Fig. 4 showed the transmit frequency and the CTR during the iterations. The CTR converged to its optimal value after six iterations for a transmit frequency of 1.2 MHz. The mean CTR achieved after convergence was around 11.7 dB, i.e. a mean contrast gain of 6.8 dB. Our simulation was in accordance with the experimental results: we observed that the optimal transmit frequency was lower than the central frequency of the transducer. Furthermore, this optimal transmit frequency obtained after the algorithm convergence also made it possible to receive the second harmonic component required in the CTR optimization.

Although the CTR was time-varying due to the non-stationarity of the contrast agents, these experimental results confirmed the feasibility of our method. Note that there was a difference between the CTR value in our simulation and the one obtained in our experiment. This difference might be due to the simplified hypotheses of the model.

VI. DISCUSSION AND CONCLUSIONS

CTR optimization in second harmonic imaging was performed automatically by estimating the transmit frequency f . This optimization was obtained through the imaging system without taking into account a priori knowledge of the medium or the transducer, except for the first three values of the transmit frequency chosen only for their impact on the speed of convergence. The optimal transmit frequency was reached automatically by a feedback within only a few iterations.

CTR had a single maximum in simulation. An easily implemented algorithm based on gradient was thus robust. Although this property was demonstrated in experiment, the proposed algorithm itself adjusted the transmit frequency to maximize the power backscattered by microbubbles while minimizing the power backscattered by the tissue within the transducer bandwidth.

One major advantage of our approach was its independence of the medium explored since the cost function was exclusively based on the input and the output measurements of our system. One interesting consequence is that our method can be applied to any imaging system.

Note that a real-time implementation was possible, since the computation time was insignificant. However, the method requires a programmable analogue transmitter. Moreover, for an efficient optimization, it is important to determine the correct position of the perfused and non-perfused areas.

To conclude, the method described ensured optimal CTR by selecting the transmit frequency. With our new approach, manufacturers and clinicians do not need themselves to tune the transmit frequency. The method automatically adapts the transmit frequency to help the diagnosis improvement.

Our closed-loop method can be adapted using a larger number of contrast imaging techniques. The only remaining difficulty is in the instrumentation. However, the current development of new imaging methods based on chirp or time reversal also requires this instrumentation.

Fig. 2 .

 2 Fig. 2. Simulation of CTR for a transmit frequency of 1 to 2.5 MHz, a pressure level A0 of 240 and 370 kPa and a constant number of cycles Nc of 4. The simulation included the transducer

Fig. 3 .

 3 Fig. 3. Simulation of automatic optimization of the CTR by iterative searching for the optimal transmit frequency for different pressure levels A0 of 240 and 370 kPa and a constant cycle number Nc of 4. The simulation included the transducer

Fig. 4 .

 4 Fig. 4. Experiment of CTR optimization for a pressure level A0 of 370 kPa at the focal length and a constant cycle number Nc of 4

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20

			a
	Frequency (MHz)	1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2 2.1 2.2	Experiment Measures at 2 MHz Average after convergence
		Iterations

ACKNOWLEDGEMENTS

The authors thank the Agence Nationale de la Recherche (Projects ANR-07-TECSAN-015 MONITHER and ANR-07-TECSAN-023 SUR-FOETUS) for financial support.

The authors would like to thank A. Zaylaa and A. Novell for helpful discussions and for language editing of the manuscript, and the Clinical Investigation Centre for Innovative Technology of Tours (CIC-IT 806 CHRU de Tours, Tours, France) for the ultrasound contrast agents.