Matthieu Dreher
email: mdreher@anl.gov

Tom Peterka
email: tpeterka@mcs.anl.gov

Bredala: Semantic Data Redistribution for In Situ Applications

In situ processing is a promising solution to the problem of imbalance between computational capabilities and I/O bandwidth in current and future supercomputers. Initially designed for staging I/O, in situ middleware now can support a wide range of domains such as visualization, machine learning, filtering, and feature tracking. Doing so requires in situ middleware to manage complex heterogeneous codes using different data structures. Data need to be transformed and reorganized along the data path to fit the analysis needs. However, redistributing complex data structures is difficult. In many cases, arbitrarily splitting the arrays of a data structure destroys the semantic integrity of the data. We present Bredala, a lightweight library to annotate a data model with enough information to preserve the semantic integrity of the data during a redistribution. Bredala allows developers to describe how to split and merge a data model safely, operations usually done by in situ middleware. We evaluate the cost and performance of our library in a molecular dynamics application. We show that our data model can simplify the workflow graph of largescale applications, improve the reusability of tasks, and offer an efficient alternative to redistribute the data.

I. INTRODUCTION

In each new generation of top 500 supercomputers, the gap between compute capabilities and I/O bandwidth grows. Extrapolating this trend to exascale, researchers estimate that less than 1% of data produced by simulation will be stored to disk [START_REF] Moreland | #! Exascale! The Effect of Emerging Architectures on Scientific Discovery[END_REF]. At this point, data analysis done in postprocessing mode will not be practical anymore because of a significant loss of data and poor I/O bandwidth. Moreover, processing this amount of data may require a similar amount of computational power as was used to generate it.

In situ data analysis is a promising solution to this problem. Treating data while they are still in memory avoids the I/O bottleneck, uses the supercomputer to perform analysis, and shortens the overall time to discovery.

Initially designed for staging I/O, in situ middleware manages workflows that merge simulations and analysis. For example cosmology simulations produce a large set of particles [START_REF] Peterka | High-Performance Computation of Distributed-Memory Parallel 3D Voronoi and Delaunay Tessellation[END_REF] that can be tessellated in order to compute a density estimation of the initial particles. At each of these steps, the original input data are transformed and possibly redistributed across different resources. That is the case when a producer hosted on N processes sends data to a consumer hosted on M processes. Scientists expect that these data transformations will become more common as the workflows become more complex.

Several in situ middleware infrastructures have been proposed. Some come from the I/O community [START_REF] Liu | Hello ADIOS: the Challenges and Lessons of Developing Leadership Class I/O Frameworks[END_REF], others from the visualization community [START_REF] Ahrens | 36 paraview: An end-user tool for large-data visualization[END_REF], [START_REF] Ahern | Visit: Experiences with sustainable software[END_REF]. Simulation data generally are composed of large arrays representing a set of points or cells with their associated properties such as speed or temperature. I/O interfaces such as NetCDF [START_REF]Unidata netcdf home page[END_REF] or HDF5 [START_REF] The | Hierarchical Data Format, version 5[END_REF] are usually well suited to represent these data. The redistribution of such data can be done by simply splitting the arrays in equal parts.

However, the data structures used in analysis can be more diverse than those of simulations. Graphs, meshes, trees, or structures with multiple levels of indirection are frequent outside the simulation world. Such structures raise two problems for in situ middleware: (1) It can be difficult to express these structures with current data interfaces inherited from traditional I/O interfaces, and (2) Redistributing such data models without breaking the semantics of the data is specific to each data structure.

Some systems allow transforming the data along the in situ pipeline [START_REF] Boyuka | Transparent I Situ Data Transformations in ADIOS[END_REF], [START_REF] Abbasi | Just in time: Adding value to the io pipelines of high performance applications with jitstaging[END_REF] as part of their communication functionalities. For most in situ middleware, the redistribution is part of the system itself. Therefore, redistributing specific data structures that require complex data manipulation can necessitate the developer to modify the in situ middleware itself for each specific data structure. In the end, it becomes difficult for in situ middleware to redistribute all the possible data models.

We present Bredala 1 , a lightweight library designed to build a data model 2 with enough information to preserve the semantics of the data during a redistribution. The key design decision of our library is to include split and merge operations inside the data model. For each field of a data structure, the developer can specify how to split a field into chunks or merge two or more chunks of the same data field. Setting the split/merge rules at the data model level instead of by the underlying in situ middleware provides more flexibility for the user to define complex data models. It also allows the same model to be used in different in situ middleware infrastructures. Bredala also provides redistribution components for common strategies such as round-robin or contiguous strategies that can be used on top of existing in situ middleware or in the middleware itself.

We evaluate Bredala inside FlowVR, middleware designed to create asynchronous in situ workflows. We improved a previous analysis algorithm called Quicksurf, a distributed isosurface algorithm designed for molecular surface visualization. In particular, we show that Bredala allows simplifying the workflow graph of the FlowVR especially at large scale, automates the redistribution of the data, gains performance during the redistribution, and improves the reusability of the workflow tasks.

The reminder of this article is organized as follows. Section II reviews related work. Section III showcases some data structures that are difficult to redistribute with current in situ systems. Section IV describes our data model, and Section V describes our redistribution components. Section VI describes Quicksurf and discusses our implementations. Section VII evaluates the performance and usability of our library. Section VIII summerizes our work and briefly describes future work.

II. RELATED WORK

We briefly survey common interfaces to describe data and transform them in common in situ middleware and I/O systems.

One of the first examples of in situ middleware is visualization. VisIt [START_REF] Ahern | Visit: Experiences with sustainable software[END_REF] and Paraview [START_REF] Ahrens | 36 paraview: An end-user tool for large-data visualization[END_REF] are widely adopted tools to visualize and analyze datasets. Both contain an in situ library, respectively libsim [START_REF] Whitlock | Parallel In Situ Coupling of Simulation with a Fully Featured Visualization System[END_REF] and Catalyst [START_REF] Fabian | The Paraview Coprocessing Library: a Scalable, General Purpose In Situ Visualization Library[END_REF], to connect simulations to their visualization pipeline. In both cases, the library serves as an adapter to convert the data format of the simulation to the VTK format [START_REF] Schroeder | The Visualization Toolkit-An Object-Oriented Approach To 3D Graphics[END_REF] at the basis of both tools.

ADIOS [START_REF] Liu | Hello ADIOS: the Challenges and Lessons of Developing Leadership Class I/O Frameworks[END_REF] is a common interface for several in situ infrastructures [START_REF] Zheng | GoldRush: Resource Efficient in Situ Scientific Data Analytics Using Fine-grained Interference Aware Execution[END_REF], [START_REF] Zheng | FlexIO: I/O Middleware for Location-Flexible Scientific Data Analytics[END_REF], [START_REF] Abbasi | Just in time: Adding value to the io pipelines of high performance applications with jitstaging[END_REF], [START_REF] Dayal | Flexpath: Type-based publish/subscribe system for large-scale science analytics[END_REF]. Data are described in a separate XML file with a name, a type, and a layout. Data can be organized by groups with a transport method for each group. The transport method can be an I/O method such as MPI-IO or an in situ system. Some of these systems enable communications from N to M processes. Several systems [START_REF] Zheng | PreDatA -Preparatory Data Analytics on Peta-Scale Machines[END_REF], [START_REF] Abbasi | Just in time: Adding value to the io pipelines of high performance applications with jitstaging[END_REF], [START_REF] Boyuka | Transparent I Situ Data Transformations in ADIOS[END_REF], [START_REF] Lofstead | Xchange: high performance data morphing in distributed applications[END_REF] also allow transforming data along the I/O path.

DataSpaces [START_REF] Docan | DataSpaces: an Interaction and Coordination Framework for Coupled Simulation Workflows[END_REF] provides a distributed shared-memory space to act as a staging area between applications. Producers push data into the shared space. The data are then indexed, and the system determines where to store them. Data can be organized following a space-filling curve, but other policies are also possible. Clients make a request for a subpart of the global data. Tong and al. [START_REF] Jin | Exploring data staging across deep memory hierarchies for coupled data intensive simulation workflows[END_REF] build on top of DataSpaces a framework staging data in memory and on disk.

Legion [START_REF] Bauer | Legion: Expressing Locality and Independence with Logical Regions[END_REF] is a data-centric programming model based on logical regions, typed collections that can be further partitioned recursively. The developer also can map where the data are placed in memory. Bauer et al. [START_REF] Bauer | Structure slicing: Extending logical regions with fields[END_REF] introduce the notion of a field in the logical region data model that allows slicing fields from the data structure.

Damaris [START_REF] Dorier | Damaris: How to Efficiently Leverage Multicore Parallelism to Achieve Scalable, Jitter-free I/O[END_REF] allocates dedicated cores or nodes to perform I/O or analytics. Data are described similar to ADIOS with a separate XML file. Data are gathered on the dedicated resources; however, the user merges the data. Damaris/Viz [START_REF] Dorier | /Viz: a Nonintrusive, Adaptable and User-Friendly In Situ Visualization Framework[END_REF] extends the XML to describe meshes and automatically transfers the meshes to Visit and Paraview for visualization.

FlowVR [START_REF] Dreher | A Flexible Framework for Asynchronous In Situ and In Transit Analytics for Scientific Simulations[END_REF] creates a graph between parallel applications. FlowVR does not provide any support for data typing; the developer writes data into untyped buffers. Special components called filters can be used to modify and route the data between tasks.

FFS [START_REF] Eisenhauer | A type system for high performance communication and computation[END_REF] is a serialization library for HPC applications. An FFS data type is closely linked to a C structure. Similar to an MPI Datatype, each field of the data structure is described in a map with a name, a size, a type, and a memory offset in the data structure. Additionally, the data types can be imported/exported to XML to allow more flexibility. EVPath [START_REF] Eisenhauer | Event-based systems: Opportunities and challenges at exascale[END_REF] is designed to be an event transport middleware layer. It creates a network of modules to compute, transform, and route the data. EVPath uses FFS to manage the data and enables data transformation at runtime. HDF5 [START_REF] The | Hierarchical Data Format, version 5[END_REF] and NetCDF [START_REF]Unidata netcdf home page[END_REF], as well as their parallel interfaces [START_REF] The | Parallel Hierarchical Data Format, version 5[END_REF], [START_REF]Parallel netcdf home page[END_REF], are two common data formats for scientific applications. Both formats are based on multidimensional arrays. Data are described in a hierarchical manner with the declaration of groups and subgroups. The user can read and write data with a path representing the placement of the data set in the hierarchy. The interfaces can be specialized for specific application domains. For example, H5hut [START_REF] Howison | H5hut: A High-Performance I/O Library for Particle-Based Simulations[END_REF] is designed to facilitate the manipulation of particle data sets.

III. EXAMPLES OF PROBLEMATIC DATA STRUCTURES

We describe in this section two different data structures from simulation and analysis. We show in particular why traditional array-splitting methods do not preserve the semantics of the data during a redistribution operation. In both cases, the goal is to redistribute the data from N to M processes following a block redistribution. We assume that N processes have data within a global spatial domain. The redistribution splits the global domain into M subdomains (or blocks) and sends data to their corresponding subdomain. After the redistribution, each of the M processes has a subset (or chunk) of the initial particles.

A. Particles structure in a molecular dynamics simulation

Figure 1 shows a subset of the global data structure from Gromacs, a well-known parallel molecular dynamics simulation package [START_REF] Pronk | Gromacs 4.5: a high-throughput and highly parallel open source molecular simulation toolkit[END_REF].

In this structure, a particle is composed of three floatingpoint values for the particle position, three floating-point values for the velocity, and one integer for the particle ID. A valid redistribution preserving the semantics requires the Fig. 2: Subset of the data structure of Tess. The structure contains the number of particles, the positions, the number of tetrahedras generated by the tessellation, and the array of tetrahedra following guarantees: (1) the position or velocity array can be split only every three floating-point values; [START_REF]The boost collection[END_REF] the association between individual positions, velocities, and particle IDs must be preserved; and (3) nbParticles must represent the current number of particles in the structure. Four steps are required to split this structure properly:

• Identify the field pos as the field representing the particle position in space. • Split the atom position so that particles of the same block are in the same chunk. • Apply the same partitioning as the atom position on the fields vel and globalId. • For each chunk update the field nbParticles with the number of particles present in the chunk.

These four operations cannot be done by simply spliting each array independently at the granularity of its elements. Splitting requires semantic information and coordinated manipulations on the arrays. By following the above procedure, however, splitting the array in N chunks does actually preserve the semantics because it does not split within a position or velocity and it does not separate atom positions from their respective velocities and global ID.

B. Tetrahedral mesh in a Delaunay tesselation

Another example is a data structure extracted from a tessellation analysis code [START_REF] Peterka | High-Performance Computation of Distributed-Memory Parallel 3D Voronoi and Delaunay Tessellation[END_REF] (Figure 2).

In this example, tets contains offsets to the array particles and neighboring tetrahedra. Therefore, simply splitting all arrays into M chunks breaks the semantics of the data structure because the offsets of the tetrehedra become invalid.

Depending on the application, we can redistribute by tetrahedra or by particles. The splitting key is the field that is first redistributed, and the rest of the fields follow the split of the key. Considering the particles as the splitting key, a proper redistribution requires three steps:

• Split the particles in M chunks according to their position. • For each tetrahedron, assign it to a subdomain only if its four vertices are assigned to the same subdomain, and adjust the offsets of the vertices. • In a second pass over the tetrahedra, adjust the offsets of the neighboring tetrahedra.

No generic in situ middleware can properly redistribute this data structure. The developer needs an interface to describe to the middleware how to split and merge this structure properly.

IV. BREDALA DATA MODEL DESCRIPTION

We describe our library, Bredala, for data redistribution. It allows the programmer to describe complex data structures with enough information to guarantee the semantic integrity of the data when it is redistributed.

We chose to integrate splitting and merging functions inside the data model instead of an underlying in situ infrastructure. The library provides common splitting and merging strategies and can be extended by the user to fit specific needs.

A. Definition of a semantic item

Bredala's goal is to protect data model semantics during a split or merge. Data structures are often a form of container for smaller pieces of data called semantic items. We define a semantic item as the smallest subset of data that contains all the fields of the original data structure and still preserves its semantics. For instance the structure in Figure 1 is a collection of particles. A semantic item in this structure is a particle. It is the combination of the number of particles (equal to one), three floating-point values for the position, three for the velocity, and one integer for the particle ID. Bredala provides a safe way to access, extract, and merge semantic items within a data model. Note that in some cases a semantic item in a data structure can be interpreted differently depending on the application. In Figure 2, both particles and tetrahedra are valid semantic items.

B. Building a data model

We define a data model as a collection of fields, where each field has a name, a data type, a layout, and split/merge policy. A field is the association of a storage container and an elementary data type. Currently, a storage container can be a singleton (only one object allowed), an array (1D, 2D, or 3D), or an STL vector. An elementary data type is a type that will not be divided during a split. It can be a basic type such as an integer or a complete data structure. We create each field with enough information to identify the subset of data for each semantic item. For instance, in Figure 1 a position is composed of three floats that form one indivisible element for this field. The split and merge policies will be described in the next section.

We build a data model for a particular data structure in two steps. First, similar to MPI data types or FFS [START_REF] Eisenhauer | A type system for high performance communication and computation[END_REF], we describe each field of the data structure (see Figure 3). However, we use a higher-level interface that simplifies the description of the data structure compared with MPI data types. For instance, we do not need to indicate the offset address of a field in a C data structure.

Bredala provides classes for simple types (float, double, int, etc.); structures; and 1D, 2D, and 3D arrays. Future versions of Bredala will also provide support for recursive data types.

Once all the fields are declared, the developer inserts them into a handler ConstructData that represents the data model. When inserting a field, the user provides (among others) a name for the field, a semantic flag, and a split and merge policy. The handler can be seen as a map of fields with the field name as key and with additional information per field. The semantic flag allows the user to indicate that a field has a particular meaning in the data model that can be used during a split/merge. We will show an example in Section V-C. We describe the split and merge policies in the next section.

Since our data model is not tightly linked to a C/C++ structure as with MPI datatypes, one can construct a model incrementally at runtime. The developer can, for instance, merge several C structures into one data model. Moreover, the map-based data model enables a form of introspection. When reading a data model, the developer does not have to know its exact structure, only the names of the needed fields and their data types. The map-based representation also enables slicing: the user can simply remove a field by its name. The user receives a pointer to the removed data field, which can be inserted into another data model or destroyed. The user can insert the same field into multiple data models. In that case, the data models have access to the same data. Bredala manages the fields with shared pointers. A field will be destroyed only if there are no more references to the field. This approach allows the memory footprint to remain as low as possible.

C. Splitting and merging the data model

The handler provides an interface to split and merge the data model safely. The data model can be split by range of items, list of individual items, or spatially. Splitting the data model results in N subdata models. A subdata model has the same fields as the original data model, but each field has only a subset of data from its original field. Merging N identical data models, that is, models with the same field names, aggregates the data from each field.

Splitting and merging are performed field by field by using the selected policies. We implemented several common split and merge policies for each collection of elementary data types. For instance, in Figure 3 the policy SPLIT DEFAULT, used for the field pos, will divide the array of positions while guaranteeing that three consecutive floating-point values will never be separated. The policy MERGE APPEND VALUES will concatenate the two arrays of position. Figure 4 gives an example of split and merge of a set of particles.

D. Extending the data model

Developers might want to extend our data model to (a) add or modify a split or merge policy to an existing elementary data type or (b) to support a new elementary data type. Adding a split or merge policy to an existing elementary data type simply requires adding the new constant name and the split or merge code into the appropriate elementary data type class.

Modifying an existing split or merge policy is done by extending the elementary data type class and overloading the split or merge function. This can be desired when creating specific data types with particular behaviors. For example, Bredala has a bounding box data type built on top of the array data type. Both the split and merge functions are overloaded to modify the behavior of the default array split and merge policies.

The developer can add the support for new data types in two ways. The first is to extend an existing elementary data type as we did to create the bounding box type; this requires little effort. The second is to extend the abstract base class used for every elementary data type. We deliberately kept its interface simple, with only five functions to overload, plus the constructor.

E. Serialization

We built Bredala as a bridge between parallel codes involving network communication. The handler and all the elementary data types implement a serialize() and deserialize() function. Our current implementation uses Boost [START_REF]The boost collection[END_REF] to serialize the data model. This choice allows Bredala to not rely on any communication library such as MPI for its serialization and facilitates its integration into a wide range of infrastructures.

V. DATA REDISTRIBUTION COMPONENTS

The core of the problem addressed by Bredala is to safely distribute a data model from N to M processes without breaking the semantics of the data model. Bredala provides redistribution components for three common strategies: roundrobin, contiguous, and block (Figure 5). The round-robin component distributes the data item by item in a cyclic manner. The contiguous component redistributes all the items while preserving their order. The block component divides an n-dimensional global domain into subdomains that are contiguous in each dimension and assigns each item to a particular subdomain. A. General structure

The redistribution components have three main functions: assigning a destination process for each item, transporting the split data models from N to M processes, and merging the received data models on the M processes. Each component is composed of sources and receivers of data. Sources and receivers can be the same processes or can be separated.

We adopted an algorithm similar to the one used by FlexIO [START_REF] Zheng | FlexIO: I/O Middleware for Location-Flexible Scientific Data Analytics[END_REF] to redistribute the data. Each component implements four functions: computeGlobal(), split(), redistribute(), and merge(). The computeGlobal() function computes the global parameters necessary to determine how to split the data. It is executed on the source side only. The split() function assigns a destination for each item, splits the local data model into subdata models, and serializes them. The redistribute() function transports the serialized data models. The transport is done in five steps:

1) Each source stores in an array the ranks of all the destinations it has to send a message. 2) The arrays of all the sources are gathered at rank 0 of the sources.

3) The aggregated array is sent to the root of the destinations. 4) The array is scattered among all the destination. At this point, all the destinations know how many messages they will receive. 5) The data messages are exchanged point to point. The components are implemented with MPI, but they could be implemented with any communication library. The only collective communications occur within the sources (MPI Gather, step 2) and within the destination (MPI Scatter, step 4) but not between the sources and destinations. We also have a point-to-point implementation available that requires no collectives. By default, steps 1 to 4 are performed at each iteration. The developer can skip these steps on demand if the data distribution has not changed between two consecutive iterations.

The merge() function aggregates all the received subdata models into one data model. It is called every time a destination receives a subdata model. It has two implementations. The first one directly deserializes and merges the received data model with the accumulated data model. The second implementation deserializes the received data model and stores it. When all the subdata models are received, they are merged all at once. The first implementation has a lower memory footprint, but it can require multiple memory reallocations when a new data model is received. On the other hand, merging all the sub data models at once requires one memory reallocation. Therefore, the second implementation runs faster at the cost of a higher memory footprint. The user must select the desired implementation.

B. Round-robin and contiguous redistribution components

The round-robin component does not require any global computations. Each semantic item is distributed according to its rank in the local data model.

The contiguous component requires two global pieces of information to perform the redistribution: the total number of semantic items and the global rank of the first local semantic item, that is, the rank of the semantic item when considering all the semantic items in all the sources. The component then computes the destination for each semantic item.

C. Block redistribution component

The block component redistributes items into subdomains with optional ghost regions. The component assumes that all the data on the sources belong to a global spatial domain. The component then assigns to each receiver a subdomain of the global domain and sends each item to their corresponding subdomains.

The basic structure to describe a domain in Bredala is a block. A block can be seen as a hexahedral region in a 3D grid. A block contains three pieces of information.

• Global grid, that is, the size of the grid of the complete space of the data • Local grid, that is, the size of the grid managed by the local data model including ghost regions, and • Own grid, that is, the size of the grid managed by the local data model excluding ghost regions.

The size of the regular grid cell is stored as well.

The computeGlobal() function reads the field domain_block added by the user from the data model and computes the subblocks for the destinations; that is, it generates new blocks and fills their local and own grid information. The split() function determines for each semantic item in which subblock(s) it belongs based on its position (semantic flag ZCURVE). The redistribute() and merge() functions follow the same pattern as do the other redistribution components.

VI. APPLICATION TO THE QUICKSURF ALGORITHM

We tested our data model with an in situ pipeline developed in previous work [START_REF] Dreher | A Flexible Framework for Asynchronous In Situ and In Transit Analytics for Scientific Simulations[END_REF], based on the Quicksurf algorithm [START_REF] Krone | Fast Visualization of Gaussian Density Surfaces for Molecular Dynamics and Particle System Trajectories[END_REF]. Quicksurf is an isosurface extraction algorithm. We integrated Bredala with FlowVR and adapted the Quicksurf pipeline to use our data model and redistribution methods. Our goal is to evaluate how our data model can be integrated into an in situ middleware infrastructure to improve its usability and performance.

A. FlowVR

The FlowVR middleware is designed to build in situ workflows. A FlowVR workflow is a graph where nodes are parallel programs (called modules) performing computations and edges are communication channels. Communication channels are created between input and output ports; each module can have several input and output ports. The developer describes each communication channel in a Python script that allows creating virtually any communication patterns. However, FlowVR provides no support for data types. Instead, the developer is given a buffer to data that must be serialized manually.

The same problem arises with data redistribution. The developer has to provide the necessary redistribution components that are usually hard coded for certain types or structures. A common practice with FlowVR to mitigate this problem is to separate different data structures into individual communication channels. For example, a structure composed of several arrays is sent through a different communication channel for each array, resulting in a more complex graph.

B. Quicksurf algorithm

The Quicksurf algorithm is a common rendering method in biology to visualize the surface of a molecule. It is a variation of the isosurface extraction method [START_REF] Lorensen | Marching cubes: A high resolution 3d surface construction algorithm[END_REF]. Quicksurf is composed of three steps. a All the atoms are inserted into a regular 3D grid according to the position of their centers. b For each grid cell, we compute a density value based on the atoms present in the cell and those in its neighbors. c We perform a marching cube algorithm on the density grid.

Gromacs uses an irregular domain decomposition to distribute atom positions across multiple processes. However, the Quicksurf algorithm cannot use directly the same redistribution because Quicksurf works on a regular domain decomposition and the isosurface requires ghost regions to compute the density values at the borders of the blocks. Consequently, atom positions are exchanged between steps (a) and (b), or density values are exchanged between steps (b) and (c) to transition between the distribution from Gromacs to the distribution necessary for the isosurface.

C. Pipeline improvement with Bredala

In [START_REF] Dreher | A Flexible Framework for Asynchronous In Situ and In Transit Analytics for Scientific Simulations[END_REF], the Quicksurf pipeline was implemented with FlowVR. Each step was a separate module. The first module computed a Morton code [START_REF]A computer oriented geodetic data base and a new technique in file sequencing[END_REF] 3 for each atom and sorted the atom positions by their Morton codes. The second module computed the grid density, and the third module performed the marching cubes. The full pipeline includes three data streams: positions, Morton codes, and grid density. A specialized redistribution module was also designed to redistribute the atom positions or grid density. The module is composed of two submodules: routers and mergers. The routers are the equivalent of the sources of the redistribution components, and the mergers are the receivers. After the redistribution, each merger has all the atom positions or density values from a subdomain. The redistribution communications are performed through FlowVR using an NxM pattern. This pattern is reproduced for each data stream.

We integrated Bredala into this pipeline. First, each data structure is reassembled into a single data model. Consequently, only one data stream is necessary between the modules. In terms of code complexity, both implementations require the same number of lines of code to create a data model: two instructions for each data field.

Besides reducing the number of communication channels, another improvement is the management of the metadata. Along the pipeline, we maintain metadata about the domain decomposition such as the global domain of the simulation, the local domain, and the own domain. With FlowVR, such metadata are managed through a different API from the Message API, adding complexity. With Bredala, the metadata are simply another field pushed into the data model.

We also reimplemented the redistribution module using the block redistribution component from Bredala. The FlowVR implementation required approximately 300 lines of code, two separate modules, and an NxM communication pattern for each data stream that is hard coded for the particular data structures used in the pipeline. The Bredala implementation requires 25 lines of code (see Figure 6) with only one module that hides all the costly communications from FlowVR. The communications are directly performed through MPI inside the component instead of the FlowVR communication channels. The module is fully generic and can be used with any data model that can be expressed with Bredala. We evaluate the performance of these implementation in Section VII.

D. Extensibility

A possible extension to Quicksurf is to color the surface based on the speed of the atoms. This requires adding the atom speed into the pipeline and modifying the grid density to incorporate a color component in addition to the density value.

Most of the pipeline is hard coded with particular data types in the FlowVR-only version. To add the atom speeds, we have to modify the graph of the workflow to add communication channels for the speed and modify the code in all the modules to forward the atom speeds along the pipeline. Extra modifications are also required in the density grid module to compute the color of a cell. The FlowVR with Bredala version does not need to modify the graph workflow. The speed is simply pushed inside the data model and forwarded across all the modules automatically. The same applies for the redistribution component. We only need to modify the module computing the density grid to compute the color of a cell.

Thus, Bredala simplified the graph of the workflow, improved the reusability of the pipeline, and simplified code of the modules.

VII. PERFORMANCE EVALUATION

In this section we evaluate the performance of Bredala when building and redistributing a data model, and we compare the results the performance of with FlowVR. We used Gromacs as our application driver.

A. Experimental context

The experiments ran on Froggy, a 190 compute node cluster from the Ciment infrastructure. 4 Each compute node has 2 eight-core processors, Sandy Bridge-EP E5-2670 at 2.6 GHz, with 64 GB of memory. Nodes are interconnected through a FDR InfiniBand network. FlowVR 2.1 and Gromacs 4.6.7 are compiled with OpenMPI 1.8.3 (no OpenMP). For all experiments using Gromacs, we run a Martini simulation with a patch of 54,000 lipids representing about 2,100,000 particles in coarse grain [1].

B. Data model cost

We modified the code of Gromacs to output data from the simulation using three different methods. In all cases, we extracted the atom positions and the atom IDs from the simulation. The first configuration uses only Bredala to measure the cost of building a data model. We built a data model composed of one array of the form [x 1 , y 1 , z 1 , id 1 , ..., x N , y N , z N , id N ] that is then serialized. The second configuration uses only FlowVR to isolate the cost of our in situ middleware. We created the same array as the Bredala configuration and copied it in a FlowVR message. The third setup uses both Bredala and FlowVR. This evaluates the cost of both building the data 4 https://ciment.ujf-grenoble.fr model and the cost of the in situ middleware and compares these costs with the cost of the two previous configurations. We built the same data model as in the Bredala setup, serialized it, and copied the serial buffer in a FlowVR message. We measured the performance of these three configurations and compared the results with the performance of the unmodified version of Gromacs (Figure 7). We used 16 cores per node for the simulation. Data were extracted every 100 iterations of the simulation.

Molecular dynamic simulations have a very high frequency, commonly reaching 1000 Hz. At this rate, the slightest perturbation on the simulation can have a noticeable impact on the simulation performance. In this context, extracting data from the simulation while impacting as little as possible is challenging. Our goal with these experiments is to evaluate the additional cost of Bredala when combined with an in situ middleware.

At all scales, Bredala is the less costly setup, with a maximum cost of 2.9% on the simulation performance at 32 nodes while FlowVR has an impact of 5.5%. This difference of cost is explained by the synchronizations that each MPI process of the simulation has to make with the FlowVR daemon every output iteration. In all the cases, Bredala adds a limited extra overhead when combined with FlowVR. In the worst case, the FlowVR and Bredala configuration adds an extra 0.4% of performance drop compared with the FlowVRonly setup on 32 nodes. These results show that Bredala can be used in conjunction with another in situ middleware framework for a slight additional cost compared with the cost of the in situ middleware itself. Note also that adding Bredala with FlowVR does not introduce more variability in performance compare to FlowVR alone.

C. Redistribution component performance

Next we benchmarked our block redistribution component and compared it with our hard-coded version using FlowVR developed for the Quicksurf pipeline. We reproduced the same data model as FlowVR with Bredala to test our redistribution component. The data model is composed of two arrays. The first one has the same structure as in Section VII-B (atom positions with atom IDs). The second array is composed of a Morton code for each atom. Upon initialization, we generated

x atoms, computed their Morton codes, and sorted the atoms by their Morton code. Atom positions are generated uniformly in a global domain space.

We first benchmarked our component in an in situ scenario. N processes are emitting data (sources) and the same processes receive the data (receivers). The sources emit atoms from a global domain, but the receivers only receive atoms from a subdomain of the global domain. With our uniform distribution of atoms, each source sends atoms to each receiver. To measure the performance of our components, we run an infinite loop redistributing the atoms. Atoms are generated only once at initialization. The components can start the next iteration only when all the receivers have received and merged their data. The component's performance is then given by the number of redistribution they can perform per second.

We measured our component's speed in both weak (Figure 8(a)) and strong scaling (Figure 8(b)). In weak scaling, each source emits 200,000 atoms. In strong scaling, 2,100,000 atoms are distributed among the sources. For each scenario, we placed one source and receiver per node. Each component computes all the necessary information for the redistribution at each iteration.

In the weak-scaling scenario, FlowVR displays the best speed at low scale with 68 Hz on eight nodes compared with 53 Hz with Bredala. At larger scale, however, the gap between both implementations shrinks, and their performance drops to 13 Hz with Bredala and 16 Hz with FlowVR at 64 nodes. With Bredala, the main reason for the drop of performance when scaling is the cost of the split() function (Figure 8(c)), which is the most time-consuming operation (58% of the total execution time at 64 nodes). When increasing the number of nodes, we increase the number of subdata models to generate, one subdata model for each receiver. This increases the amount of operations necessary to divide and concatenate a data model. FlowVR also suffers from the increased complexity to split the initial data (Figure 8(e)). However, the time spent doing communication increases significantly as well. At 32 nodes, the split and the communications represent 46% and 40%of the total time. With FlowVR, each communication is described in the workflow graph. Adding a node generates 2xN new communication channels with N the number of nodes (two data structures to transport). This increases the amount of work done by the FlowVR daemon and slows the communications managed by the daemon.

In the strong-scaling scenario, FlowVR also displays a better speed at low scale, with 56 Hz on eight nodes against 43 Hz for Bredala. However, the performance of FlowVR starts to drop beyond 16 nodes to reach 29 Hz on 64 nodes. We observe that the time spent in communication increases significantly from 5 ms at 8 nodes to 25 ms at 64 nodes (Figure 8(f)). As with the weak-scaling case, the increase in communication time is explained by the increased complexity of the workflow graph. Even though the time to split the data is decreasing when scaling up, it is not enough to compensate for the increased communications. Bredala is able to scale up to 90 Hz at 32 nodes but drops to 35 Hz at 64 nodes. We observed that both the times of splitting the data and communications are decreasing up to 32 nodes (Figure 8(d)) but then increase at 64 nodes. Future work will improve the manipulation of array to avoid unnecessary copies during split operations.

Overall, we see that except at low scale, the Bredala redistribution component scales better than FlowVR while providing a fully generic redistribution pattern.

D. Quicksurf pipeline

We also compared the performance of the hard coded version of the Quicksurf algorithm with our implementation using Bredala as described in Section VI-C. Gromacs runs on 15 cores per node, and we use the last core to compute Quicksurf. We triggered the in situ pipeline every 100 iterations of Gromacs. We configured the pipeline so that the redistribution component exchanges atom positions.

Our first experiment runs a full in situ scenario where the complete analysis pipeline is hosted on the simulation nodes and the redistribution is performed between the simulation nodes. Figure 9 compares the performance of both implementations with the performance of Gromacs without modifications. Except at 16 nodes, both implementations have a similar cost for the simulation. At 64 nodes, both implementation costs 30% of the simulation performance.

Our second experiment uses staging nodes to compute the analysis after the redistribution. After the computation of the Morton codes, the data are routed to the staging nodes to compute the grid density and the marching cubes. Results are summarized in Figure 10. We allocate one node as a staging node for every 32 simulation nodes. As previously, the performance of both implementations is close, with the Bredala implementation slightly less costly in all the cases. At 64 nodes (62 simulations, 2 staging nodes), the FlowVR implementation costs 23% of the simulation performance while the Bredala implementation costs 16%.

The performance improvement brought by Bredala is overall quite small: 7% performance improvement in the best case and less than 2% otherwise. This is surprising considering the performance difference between the two redistribution component alone (Figure 8). Two main facts explain these results. First, most of the impact on the simulation performance comes from the data extraction from the simulation and the extra in situ computations that impact the cache. Second the simulation extracts data at only six iteration per second on 64 nodes. That data rate is slower than the peak data rate manageable by the components.

We showed that converting the pipeline to use Bredala does not add a significant cost to the simulation performance compared with what FlowVR already costs. Moreover, Bredala enables more generic and maintainable code that is easier to extend. In some cases, Bredala can even improve the impact on the simulation performance with more efficient data redistribution components.

VIII. CONCLUSION

We introduced Bredala, a lightweight library designed to protect the semantics of data during parallel redistribution in in situ applications. The user describes a data model field by field and adds extra information to allow the library to isolate individual semantic items and manipulate them during parallel redistribution. Bredala comes with several common redistribution patterns such as spatial redistribution. We compared our library with FlowVR, and we studied the benefits of combining FlowVR with Bredala and applied them to the Quicksurf algorithm. We showed that combining FlowVR with Bredala does not add a significant overhead compared with the overhead of FlowVR itself while improving greatly the reusability of the FlowVR modules.

So far, Bredala uses only the main memory for certain operations such as serialization. Several in situ middlewares such as FlowVR or Damaris expose specialized allocators to transfer data to the middleware. Future version of Bredala will provide an allocator interface to use directly the allocator provided by in situ middleware and avoid unnecessary copies. We will also extend the current set of basic types for data field to n-dimensional arrays.

Fig. 1 :

 1 Fig.1: Subset of the data structure of Gromacs. The structure contains the number of particles, the position, velocity, and the global ID of each particle.

Fig. 3 :

 3 Fig. 3: Code snapshot to build the data model corresponding to the data structure of Figure 1. First, an object is created for each field. Then the objects are pushed in the container representing the data model. Each field is associated with a name, optionally a semantic flag (ZCURVEKEY for instance), and a splitting/merging policy.

Fig. 4 :

 4 Fig.4: Initial data model D0 containing three particles. We split D0 in two sub data models D1 and D2. D1 contains the particle 1 and D2 contains the particles 0 and 2. D0 remains unmodified. We then merge in place D2 into D1, but the destination could be another data model as well.

Fig. 5 :

 5 Fig. 5: Redistribution components available in Bredala. (a) The round-robin component distributes the data item by item in a cyclic manner. (b) The contiguous component redistributes all the items while preserving their order. (c) The block component divides a global domain into subdomains and attributes each item to a particular subdomain.

Fig. 6 :

 6 Fig. 6: Code snapshot to build the redistribution module with Bredala.

Fig. 7 :

 7 Fig. 7: Performance measurements of Gromacs with several instrumentation methods.

Fig. 9 :

 9 Fig. 9: Performance measurements of Gromacs with the Quicksurf pipeline implemented with FlowVR and Bredala with an in situ strategy. Gromacs performance stand-alone performance is given as reference.

Fig. 8 :

 8 Fig. 8: Redistribution components performance and time decomposition with Bredala and FlowVR.

Fig. 10 :

 10 Fig. 10: Performance measurements of Gromacs with the Quicksurf pipeline implemented with FlowVR and Bredala with an in transit strategy. Gromacs performance stand-alone performance is given as reference.

The traditional cookie of the lead author's home and a favorite of both authors

A Morton code is a hash of a cell coordinate.

ACKNOWLEDGMENT

This material was based upon work supported by Advanced Scientific Computing Research, Office of Science, U.S. Department of Energy, under Contract DE-AC02-06CH11357, program manager Lucy Nowell. All the computations presented in this paper were performed using the Froggy platform of the CIMENT infrastructure (https://ciment.ujf-grenoble.fr), supported by the Rhône-Alpes region (GRANT CPER07 13 CIRA) and the Equip@Meso project (reference ANR-10-EQPX-29-01) of the programme Investissements d'Avenir supervised by the ANR.


Ph.D Huazhong Z Shu

J Zhou

G N Han

L M Luo

J L Coatrieux

Image reconstruction from limited range projections using orthogonal moments

Keywords: Image reconstruction, Radon transform, projection moments, image moments, orthonormal polynomials

A set of orthonormal polynomials is proposed for image reconstruction from projection data. The relationship between the projection moments and image moments is discussed in detail, and some interesting properties are demonstrated. Simulation results are provided to validate the method and to compare its performance with previous works.

Introduction

In its classical formulation, computerized tomography (CT) deals with the reconstruction of an object from measurements which are line integrals of that object at some known orientations. This formulation has found many applications in the fields of medical imaging, synthetic aperture radar, electron-microscopy based tomography, etc [START_REF] Deans | The Radon transform and some of its applications[END_REF]. The Radon transformation, due to its explicit geometric meaning, has played an important role.

In the past decades, a number of studies have been carried out on this subject. Lewitt has summarized a series of projection theorems [START_REF] Lewitt | Reconstruction algorithm: transform methods[END_REF], a detailed description about the properties of Radon transform and its relationship to other transforms have been given by Deans [START_REF] Deans | The Radon transform and some of its applications[END_REF]. Major contributions have been reported in the biomedical engineering literature, one of the most active fields for the reconstruction problem, in order to study new X-ray source-detector trajectories (typically a spiral for CT) and to deal with truncated [START_REF] Defrise | Image reconstruction from truncated, two-dimensional, parallel projections[END_REF] and cone-beam projections [START_REF] Grangeat | Indirect cone-beam three dimensional image reconstruction[END_REF].

Recently, the moment-based approaches to tomographic reconstruction have attracted considerable attention of several research groups. Salzman [START_REF] Salzman | A method of general moments for orienting 2D projections of unknown 3D objects[END_REF] and Goncharev [START_REF] Goncharev | Methods of integral geometry and recovering a function with compact support from its projections in unknown directions[END_REF] respectively proposed the methods based on the moments to find the view angle from the projection data. Milanfar et al. [START_REF] Milanfar | A moment-based variational approach to tomographic reconstruction[END_REF] described a variational framework for the tomographic reconstruction of an image from the maximum likelihood estimates of its orthogonal moments. Basu and Bresler [START_REF] Basu | Uniqueness of tomography with unknown view angles[END_REF][START_REF] Basu | Feasibility of tomography with unknown view angles[END_REF] discussed the problem of recovering the view angles using moments of the projections. Wang and Sze [START_REF] Wang | The image moment method for the limited range CT image reconstruction and pattern recognition[END_REF] proposed an approach based on the relationship between the projection moments and the image moments to reconstruct the CT images from limited range projections. In the Wang's algorithm, the geometric moments were used and interesting results have been obtained. However, the use of the geometric moments has the following disadvantages: [START_REF] Deans | The Radon transform and some of its applications[END_REF] the geometric moments of an image are integrals of the field shape over space, and the image can be uniquely determined by the geometric moments of all orders. They are sensitive to digitization error and minor shape deformations [START_REF] Shen | Discriminative wavelet shape descriptors for recognition of 2D pattern[END_REF]; [START_REF] Lewitt | Reconstruction algorithm: transform methods[END_REF] the geometric moments are basically projections of the image function onto the monomials x n y m . Unfortunately, the basis set {x n y m } is not orthogonal. These moments are therefore not optimal with regard to the information redundancy and other useful properties that may result from using orthogonal basis functions.

To overcome these inconvenients, we propose in this paper, by extending Wang's algorithm, a new moment-based approach using the orthogonal basis set to reconstruct the image from limited range projection.

The paper is organized as follows. A brief review of Radon transform and the definition of projection moments and image moments are given in Section 2. In Section 2, we also establish the relationship between projection moments and image moments and discuss how to estimate the projection moments at any specific view from image moments. Simulation results are provided in Section 3. Section 4 concludes the paper and provides some additional perspectives.

Method

We first sketch the basics of Radon transform. The orthogonal projection moments, defined in terms of normalized polynomials, are then introduced in subsection 2.2. Some theorems relating projection and image moments are reported and demonstrated in subsection 2.3. In the last subsection, we show how to estimate the unknown projections from the calculated image moments.

Radon transform

Let f(x, y) ! L 2 (D) be a square-integrable function with support inside the unit circle D in the plane and g(s, !) be the Radon transform of f(x, y) defined as follows

"" ! + = D dxdy s y x y x f s g ) sin cos ( ) , ( ) , ( # # $ # (1)
where "(") denotes the Dirac delta function, s is the distance from the origin to the ray, and ! is the angle between the x-axis and the ray.

For a given view !, the two-dimensional (2D) function g(s, !) becomes a one-variable of s, denoted by g ! (s). Since g ! (s) represents a collection of integrals along a set of parallel rays, it is also called parallel projection of g(s, !) at view ! [START_REF] Wang | The image moment method for the limited range CT image reconstruction and pattern recognition[END_REF]. The Radon transform given by Eq. ( 1) can be rewritten as

"" ! + = D dxdy s y x y x f s g ) sin cos ( ) , ( ) ( # # $ # (2)

Orthogonal projection moments and image moments

The moments of g ! (s) are called projection moments in the Radon domain [START_REF] Lewitt | Reconstruction algorithm: transform methods[END_REF]. In this paper, we use a set of orthonormal polynomials instead of the set of monomials {s p } to define the projection moments. Let {P p (s)}, p = 0, 1, 2, … #, be a set of orthonormal polynomials defined on the interval [-1, 1], the pth order orthonormal projection moment of g ! (s) is defined as

ds s g s P L p p ) ( ) ( ) ( 1 1 ! ! " # = (3)
Let # nm be the (n + m)th order orthogonal moment of f(x, y) defined as [

! + = ! + = "" " " "" ! ! # # $ # # $ # (5)
Using the property of delta function, we have

!! + = D p p dxdy y x f y x P L ) , ( ) sin cos ( ) ( " " " (6)
This last equation will allow us to establish a relationship between the orthogonal projection moments defined by Eq. ( 3) and the orthogonal moments of f(x, y) defined by Eq. ( 4). This is the objective of the next subsection.

Relationship between projection moments and image moments

Let us first introduce some basic definitions. Let the pth order normalized polynomials P p (t) be

! = = p r r pr p t c t P 0 ) ( (7)
and let

=
where the subscript T indicates the transposition, then we have

) ( ) ( M V p p p = ( 8
)
where C p = (c kr ), with 0 $ r $ k $ p, is a (p + 1) % (p + 1) lower triangular matrix.

Since all the diagonal elements of C p , c kk , are not zero, the matrix C p is non-singular, thus

) ( ) ( ) ( 1 V D V C t M p p p p p = = ! (9)
where D p = (d kr ), with 0 $ r $ k $ p, is the inverse matrix of C p .

Eq. ( 9) can be rewritten as

! = = k r r kr k t P d t 0 ) ( , for 0 $ k $ p (10)
We then have the following theorem about the projection moments and image moments:

Theorem 1. The orthogonal projection moment of order p at given view !, L p (!), can be expressed as a linear combination of image moments of same order and lower, i.e.,

!!

= " = = p n n p m nm nm p p L 0 0 ) , ( ) ( # $ µ $ (11)
where

( ) ( ) r q m r n m n p q q r m r q m n r n m n q p nm r n q m n d d c p " + + + " = = " + + + + " " # # $ % & & ' ( + + + = ) ) ) µ sin cos ) , ( ) ( 0 0 , , , (12)
The proof of Theorem 1 is deferred to Appendix A.

Theorem 1 can also be expressed in matrix form. To do this, let us introduce the notations:

[ ] T M M L L L ) ( ),..., ( ), ( ) ( 1 0 ! ! ! ! = " , [ ] T k k k k k 0 1 , 1 1 , 1 0 ) ( , ,..., , ! ! ! ! ! " " = , and
[ ] M M ) ( ) 1 ( ) 0 ( ,..., , ! ! ! = "
where M denotes the maximum order of moments we want to utilize. Then, we have

M M M ! ! = " ) ( ) ( # # (13)
where T M (!) is a matrix of size (M+1) % (M+1)(M+2)/2 whose definition is given by

! ! ! ! ! ! " # $ $ $ $ $ $ % & = ! ! ! ) , ( ) , ( ) , ( ) , ( ) , ( ) , ( ) , 2 ( ) , 2 ( ) , 2 ( ) , 2 ( ) , 2 ( ) , 2 ( ) , 1 ( ) , 1 ( ) , 1 ( ) , 0 ( )
! i in [0, $), 1 $ i $ m, one can uniquely determine the first m moment vectors # (k) of f(x, y), 0 $ k < m.
This can be done using only the first m orthogonal moments L p (! i ), 0 $ p < m of the projections.

For the proof of Proposition 1, we refer to Milanfar et al. [START_REF] Milanfar | A moment-based variational approach to tomographic reconstruction[END_REF].

Note that Theorem 1 is valid for any type of polynomials and therefore the reconstruction based on geometric moments is a particular case of the problem considered here. In fact, for the geometric moments, the coefficients c kr and d kr in Eqs. [START_REF] Milanfar | A moment-based variational approach to tomographic reconstruction[END_REF] and [START_REF] Wang | The image moment method for the limited range CT image reconstruction and pattern recognition[END_REF] are given by

kr kr kr d c ! = = (15)
where " kr is the Kronecker symbol.

In this case, Eq. ( 12) becomes

( ) ( ) ( ) ( ) n p n m n m n p nm n p c n m n p ! + " " # $ % % & ' = " " # $ % % & ' + = ( ( ( ( ( µ sin cos sin cos ) , ( , (16)
In the remaining part of the paper, we focus on the use of orthogonal polynomials. In particular, we are interested in the normalized Legendre polynomials. The pth order normalized Legendre polynomial P p (x) over [-1, 1] is defined by

" = = ! + = p r r pr p p p p p x c x dx d p p x P 0 2 ) 1 ( ! 2 1 2 1 2 ) ( (17)
with

[ ][ ] ! " ! # $ % % + % + + % = % odd for 0, even for , ! ! 2 / ) ( ! 2 / ) ( )! ( 2 1 2 1 2 ) 1 ( 2 / ) ( r p r p r r p r p r p p c p r p pr (18)
An essential step when applying a special polynomial to the reconstruction problem is to find the inverse matrix D p in Eq. ( 9). For the normalized Legendre polynomials, we have Proposition 2. For the lower triangular matrix C p whose elements c kr are defined by Eq. ( 18), the elements of the inverse matrix D p are given as follows

[ ] ! ! " ! ! # $ % % + + % + = & % = % odd for 0, even for , )! 2 ( ! ! ) 1 2 2 ( ! 2 / ) ( 2 1 2 2 2 / ) ( 1 2 1 2 3 r k r k r r k j r r k r d r k j k r kr , 0 $ r $ k $ p (19)
The proof of Proposition 2 is deferred to Appendix A.

Based on the above proposition, we are now ready to give an explicit expression of µ nm (p, !) in Eq. ( 12) for orthonormal Legendre polynomials. Let Eq. ( 12) can thus be rewritten as

! ! " # $ $ % & + + + = ' + + + + r n q m n d d c m n p m r q m n r n m n q p qr , , ,
( ) ( ) ! ! + " = = " + + = ) ( 0 0 sin cos ) , , ( ) , ( m n p q q r r q m r n qr nm m n p p # # $ # µ (22)
where the notation "

+ ! = ) ( 0 m n p q
stands for the summation with respect to the even value of q varying from 0 to p -(n+m).

Substitution of Eqs. (18) and (19) into (20) yields

[ ] [ ][ ][ ]! 2 / ) ( ! 2 / ) ( ! 2 / ) ( )! 2 / ( )! ( )! 2 ( )! 2 ( ! ! ) 1 2 2 ( ) 1 2 2 ( 1 2 ) 1 ( ) 1 2 )( 1 2 ( ) 1 2 ( 2 ) , , ( 2 / 1 ) ( 1 2 / 2 / ) ( q m n p q m n p r q r q m n p m n m n j m i n m n p m n p r i r q j q m n p q m n p qr + + + ! ! ! ! + + + " + + + + ! + + + = # # = ! = + ! ! + + ! " (23)
The apparent complexity of Eq. ( 23) can be omitted by the following recurrence relations

) , , ( ) 3 2 )( 2 ( ) )( 1 ( ) , , ( , 2 m n p r q m r q q m n p q m n p m n p qr r q " " + ! + + ! ! ! ! + + + + ! = + (24) ) , , ( ) 3 2 )( 2 ( ) 1 2 )( ( ) , , ( 2 ,
m n p r n r r q m r q m n p qr r q " "

+ + + + ! + ! ! = + (25) ) , , ( 3 2
+ + + + ! ! + + ! = + (26) ) , , ( ) 3 2 )( 1 2 ( ) 1 2 )( 1 2 ( ) 1 , 1 , ( 00 00 m n p n n m m m n p " " + + + ! = ! + (27) [ ][ ]! 2 / ) ( ! 2 / ) ( )! ( )! 2 ( )! 2 ( ! ! 2 ) 1 ( ) 1 2 )( 1 2 ( ) 1 2 ( 2 ) , , ( 2 /
+ + ! ! + + ! + + + = ! ! ! ! " (28)

Estimation of the projection moments at any view

In real situations, projections from certain views are known while projections from other views are unknown. To distinguish them, we adopt here the same terminology used in [START_REF] Wang | The image moment method for the limited range CT image reconstruction and pattern recognition[END_REF]. The projections are called given projections if the views are given views; the unknown projections correspond to the projections whose views are unknown. In the previous subsection, we have discussed the relationship between the given projection moments and the image moments. The following shows how to compute the unknown projections from the image moments.

When the orthonormal polynomials P p (s) are used in Eq. ( 3), the orthogonality property leads to the following inverse transform

! " = = 0 ) ( ) ( ) ( p p p s P L s g # # (29)
Eq. (29) calculates g ! (s), or g(s, !), from its projection moments. If only the projection moments of order up to M are used, Eq. ( 29) is approximated by

! = = M p p p s P L s g 0 ) ( ) ( ) ( ~" " (30)
Substituting Eq. ( 11) into (30), we have

!! ! = " = = = p n n p m nm nm M p p p s P s g 0 0 0 ) , ( ) ( ) ( ~# $ µ $ ( 31
)
where µ nm (p, !) are defined by Eq. [START_REF] Teague | Image analysis via the general theory of moments[END_REF]. Note that when the orthonormal Legendre polynomials are used, Eqs. ( 22) and (23) allow to compute µ nm (p, !).

The above formula establishes the connection between projections from any specific view and image moments. To summarize, the projections of unknown views can be estimated from given projections as follows 1) Compute the image moments of order up to M from given M+1 views;

2) Calculate the unknown projections from the image moments of order up to M by using Eq. (31).

Results

Simulations

We first apply the proposed solution to a simulated phantom f(x, y), shown in Fig. 1, in order to compare the two methods for computing the image moments. The first method is directly based on Eq. ( 4) since the image is known, the orthonormal Legendre moments, # nm , of order up to 15 are calculated using Eq. ( 4). The other one is achieved through the computation of projection moments. A projection simulation program is used to generate all the projections from 15 given views and the image moments can then be computed using Proposition 1. The moment values using these two methods are shown in Table 1, where the values without parentheses are computed from the image itself, and the values with parentheses are obtained from projection moments. It can be seen from this table that the errors between the real image moments and calculated image moments are small.

When the image moments of order up to M are calculated, Eq. ( 31) is used to estimate the projection moments, ) ( ~s g ! , at any specified view !. In order to test the robustness of the proposed method, we first consider the case where all the projections are known in the interval [0, $). The estimation of ) ( ~s g ! from the image moments with maximum order M = 15 at views of 0 0 , 15 0 , 60 0 , 90 0 , 120 0 and 165 0 is shown in Fig. 2. For comparison, the real projections (computed by a projection simulation program) and the estimation results from these views using the geometric moment-based method proposed by Wang and Sze [START_REF] Wang | The image moment method for the limited range CT image reconstruction and pattern recognition[END_REF] are also shown in this figure. Fig. 2 shows that the results obtained by the proposed method and Wang's algorithm are almost the same. We then assume that projections are available in the range of 20 0 $ ! $ 160 0 , but unknown in the ranges 0 0 $ ! < 20 0 and 160 0 < ! $ 179 0 . The results using both the proposed method and Wang's method are illustrated in Fig. 3. It can be observed that our method performs better than Wang's method, especially at views 0 0 , 15 0 and 165 0 .

In the second example, we consider the image reconstruction problem from limited range projections. The filtered backprojection (FBP) algorithm is used in the reconstruction process. Since this algorithm requires the projection, g ! (s) or g(s, !), for all s and all ! in the whole interval [0, $), we use the image moment method to estimate the unknown projections. The original grey level image of size 128 % 128 is depicted in Fig. 4(a) (it is the same as that used in the previous example), the projections being assumed known from 25 0 to 155 0 . A projection simulation program is applied to compute the given projection in which the following configuration is used: (a) the total view varies from 0 0 to 179 0 ; (b) the angular sampling rate is 1 0 per view; (c) the ray sampling rate is 128 rays per view. Fig. 4(b) shows the reconstructed image from incomplete projections (25 0 -155 0 ) using the FBP method.

We then use the method proposed in the previous section as well as the Wang's method for calculating the image moments with maximum order M equaling to 5, 10, 15, and 20, respectively. These moments are then used to estimate the unknown projections.

Based on the given projections and the estimated projections (0 0 -25 0 and 155 0 -179 0 for this example), the FBP method is applied to reconstruct the image. The reconstructed images based on both the proposed method and Wang's method are illustrated in Fig.

4(c)-(i).

It can be visually seen that our method leads to better results for the same value of M. Note that the quality of the reconstructed image for Wang's method with M = 20 is poor due to the numerical instability, so that we have not shown this result in Fig. 4. We also use the mean square error (MSE) to qualitatively measure the difference between the original image and the reconstructed image. The MSE is defined as follows

% 100 (%) 2 2 * ! " = f f f MSE ( 32
)
where f is the original image, and f * is the reconstructed image.

Table 2 shows the MSEs for the proposed method and Wang's method with different values of M. The values of MSE again demonstrate the better performance of the proposed method compared to Wang's method.

Evaluation on "pseudo-real" data

The projection data are not made available by the imaging manufacturers, thus the experimentation has been conducted using a reconstructed Computed Tomography image depicted Fig 5(a). The image intensity values (e.g. density levels) have been inverted in order to enhance the internal details of the brain slice and the tumor features. We have applied a two steps procedure: (i) the projection set has been computed by using a parallel geometry (thus the so-called "pseudo-real" data) over the limited angular range; (ii) the above reconstruction algorithms have then been used.

In this example, the projections are assumed to be known on the interval [%, 180 0 -%]

where % is an adjustable parameter. We apply the proposed method to compute the image moments of order up to M = 25 for % equaling to 5 0 , 10 0 , 15 0 , 20 0 , 25 0 and 30 0 , respectively. These moments are then used to estimate the unknown projections according to (31). The reconstructed images based on both the proposed method and FBP method are depicted in Fig. 5(b)-(m). Table 3 shows the corresponding MSEs for both methods. Note that in the images reconstructed using FBP method (Fig. 5(b)-(g)), the image intensity value of the unknown projections is assigned to 0. It can be observed from Fig. 5(h)-(m) that the proposed method leads to good results even for large % .

In order to compare the performance of the proposed method with Wang's method, we calculate the image moments of order up to M = 15 for different values of % (% = 5 0 , 10 0 , 15 0 , 20 0 , 25 0 and 30 0 ); these moments are then used to estimate the unknown projections. The MSEs of the reconstructed imaged using FBP method for both methods are illustrated in Fig. 6. Fig. 6 points out the better behavior of our method over the geometric moment-based method presented in [START_REF] Wang | The image moment method for the limited range CT image reconstruction and pattern recognition[END_REF] in terms of the reconstruction capability from limited range projections.

Discussion and Conclusion

In this paper, by extending Wang's method, we have proposed a new approach based on the orthogonal moments to solve the tomographic reconstruction problem from limited range projections. By demonstrating some properties of Legendre polynomials, we have established the relationship between the orthogonal projection moments and orthogonal image moments. This relationship permits us to estimate the unknown projections from the computed image moments so that the reconstruction problem from limited range projections can be well solved.

It has been shown that orthogonal moments have some merits in comparison with the geometric moments. First, the geometric moments, especially at high order, are sensitive to noise and digitization error. Second, the orthogonal moments have simple inverse transform, thus, the image can be more easily reconstructed from the orthogonal moments.

This new approach was evaluated through a fully known phantom and a sample of "pseudo-real" data set. Experimental results have shown that our method is efficient, and provides better reconstruction results compared to Wang's method.

Reconstruction from limited range projections is of relevance for a number of medical imaging modalities, in particular for X-ray devices. Systems such as Electron Beam Computed Tomography (EBCT) and Rotational-X (Rot-X) make use of a reduced set of projections and allow decreasing the radiation dose for the patient, an important health care constraint. However, the use of sparse data leads to more visible reconstruction artifacts. A compromise is then required in order to preserve the relevant information. We have not addressed in this paper the issue regarding the minimal number of range projections from which a reasonable reconstruction could be obtained using the above methods. This is a difficult problem because it can not be addressed by MSE computation alone; subjective analysis by physicians in order to evaluate the medical consequence of a loss of informational contents is also of major importance. It also depends on the targeted use: when diagnosis requires high reconstruction performance, interventions may accept lower precision in reconstruction as far as faster algorithms are available. The works in progress are therefore devoted to such issues and also the comparison with other reconstruction algorithms (in complexity for instance) and the performance that can be expected for vascular reconstruction when dealing with a modality like Rot X. They also focus on the extension of the method to other projection modes such as cone-beam projections.

Appendix A

Proof of Theorem 1. Using Eq. ( 7), Eq. ( 6) can be rewritten as ( ) ( )

dxdy y x r q c y x f dxdy y x c y x f L D p q q r r q r r q r pq D p q q pq p !! "" !! " = = # # = $ $ % & ' ' ( ) = + = 0 0 0 sin cos ) , ( ) sin cos ( ) , ( ) ( * * * * * (A1) where )! ( ! ! r q r q r q ! = " " # $ % % & '
is the combination number.

Using Eq. ( 10), we have

( ) ( ) ( ) ( ) ( ) ( ) !! "" " " !! "" "" "" "" !! = # = + = # = # # = = = # = # # = # = # # = = $ $ % & ' ' ( ) = $ $ % & ' ' ( ) = $ $ % & ' ' ( ) =
q r m r q n r m m r q n n r m n q p p n r m n q d d c L ! " " " ## # # = $ = + $ = = + $ + + $ + + + % % & ' ( ( ) * + + + = 0 0 ) ( 0 0 , , , sin cos ) ( (A3)
The pr oof is now achieved.

Proof of

)! 2 ( )! ( )! ( )! 2 2 ( 2 ) 1 ( 2 1 4 2 2 , 2 v v u v u v u u c c u v u v u kr + ! + ! + = = ! (A5) )! 2 2 2 ( )! ( )! 1 ( )! 2 ( 2 ) 1 4 ( 2 )! 4 ( )! 2 ( )! 2 ( ) 1 2 4 ( )! ( 2 1 4 2 1 2 ) ( 1 3 2 , 2 + + ! + + + = + + ! + = = + ! = ! " w v w v w v v w w w v j w w v w d d w w v j w v w v rl (A6) thus " " " = + ! = = + + !

Proposition 1 .

 1 Given line integral projections of f(x, y) at m different angles

!

  from Eq. (19), it is obvious that for q or r being odd nu mber (21)

Proposition 2 .

 2 To prove the proposition, we need to demonstrate the following relation k, l $ p (A4) From Eqs. (18) and (19), it is obvious that Eq. (A4) is true when the two integers k and l have different parity. Thus, we need only to demonstrate it for the following two cases: (1)k, r and l are all even numbers; (2) k, r and l are all odd numbers. We consider the first case. Let k = 2u, r = 2v and l = 2w, we have

Figure 1 :

 1 Figure 1: Computer simulated phantom.

Figure 2 :Figure 3 :Figure 4 .Figure 5 .Figure 6 .

 23456 Figure2: Results when all the projections are known. The projections are compared at different views: 0 0 , 15 0 , 60 0 , 90 0 , 120 0 , 165 0 . Image moments of order up to 15 are used. (solid line: original projections; dash line: projections estimated from Legendre moments; cross: projections estimated from geometric moments)

Table 1 :

 1 Comparison of image moment values obtained from Eq. (4) (without parentheses) and image moment values estimated from projection moments (with parentheses)

	n	m	0	2	5	7
		0	0.80200 (0.80200)	-0.59407 (-0.59396)	-0.079439 (-0.079385)	0.037834 (0.037785)
		2	-0.67940 (-0.67929)	0.47723 (0.47706)	-0.086657 (-0.086584)	-0.052171 (0.052094)
		4	0.33773 (0.33757)	-0.19324 (-0.19310)	-0.082248 (-0.082153)	0.049279 (0.049188)
		8	0.059616 (0.059581)	-0.02436 (-0.02436)	-0.065381 (-0.065233)	0.038325 (0.038214)

Table 2 :

 2 Comparison of reconstruction MSE (%) using geometric moment method and orthonormal moment method

	Moment order	5	10	15	20
	MSE (Geometric)	11.7358	9.8997	27.4587	--
	MSE (Legendre)	11.7344	9.8863	6.8392	6.5655

Table 3 :

 3 Comparison of reconstruction MSE (%) using FBP method and orthonormal moment method with maximum order M = 25

	Value of %	5	10	15	20	25	30
	MSE (FBP)	2.3462	7.0099	13.6023	22.1117	32.7264	45.7178
	MSE (Legendre)	0.8022	1.4430	2.1014	2.7697	3.4405	4.0936

Acknowledgements: This work was supported by National Basic Research Program of China under grant, No. 2003CB716102, the National Natural Science Foundation of China under grant No. 60272045 and Program for New Century Excellent Talents in University under grant No. NCET-04-0477.

Here µ nm (p, !), with 0 $ n, m, p $ M, are defined by Eq. [START_REF] Teague | Image analysis via the general theory of moments[END_REF].

Theorem 1 not only provides the solution for finding the projection moments from image moments, but also the solution for finding image moments from projection moments. Concerning this latter point, we have the following proposition:

where

it can be verified

It can be easily deduced from Eqs. (A9) and (A11) that

The proof of the proposition when k, r and l are even numbers is now complete. For k, r and l being odd numbers, the proposition can be demonstrated in a similar way and is not given here.

Note that the proof of Proposition 1 was inspired by a technique proposed by Zeilberger [START_REF] Petkovsek | The book is available on line at the University of Pennsylvania[END_REF]. His experience is related to 3D images, signal processing, pattern recognition, computational modeling and complex systems with applications in integrative biomedicine. He published more than 300 papers in journals and conferences and edited many books in these areas. He has served as the Editor-in-Chief of the IEEE Transactions on Biomedical Engineering (1996)(1997)(1998)(1999)(2000) and is in the Boards of several journals. Dr. COATRIEUX is a fellow member of IEEE. He has received several awards from IEEE (among which the EMBS Service Award, 1999 and the Third Millennium Award, 2000) and he is Doctor Honoris Causa from the Southeast University, China.

About the