Charl Botha

Gordon Kindlmann

Wiro Niessen

Bernhard Preim

Hans-Jörg Schulz

Mathias John

Andrea Unger

Heidrun Schumann

Visual Analysis of Bipartite Biological Networks

Keywords: Descriptors (according to ACM CCS): I, 3, 8 [Computing Methodologies]: Computer Graph-icsApplications; H, 5, 2 [Information Systems]: Information Interfaces and PresentationUser Interfaces

In life sciences, the importance of complex network visualization is ever increasing. Yet, existing approaches for the visualization of networks are general purpose techniques that are often not suited to support the specific needs of researchers in the life sciences, or to handle the large network sizes and specific network characteristics that are prevalent in the field. Examples for such networks are biomedical ontologies and biochemical reaction networks, which are bipartite networks -a particular graph class which is rarely addressed in visualization. Our table-based approach allows to visualize large bipartite networks alongside with a multitude of attributes and hyperlinks to biological databases. To explore complex network motifs and perform intricate selections within the visualized network data, we introduce a new script-based brushing mechanism that integrates naturally with the interlinked, tabular representation. A prototype for exploring bipartite graphs, which uses the proposed visualization and interaction techniques, is also presented and used on real data sets from the application domain.

Introduction

In the fields of biology and medicine, large databases have been set up over the last years to allow access to the ever increasing flood of publicly available biological data (http://biodatabase.org). These databases are partially interconnected through cross references, and provide a wealth of information to the research community -as long as they can find it. While many of the databases can easily be queried through special interfaces and some of them even provide static depictions of their data (e.g. drawings of pathways or heatmaps for gene expression data), interactive means of visualizing and querying biomedical data bases are still a current research topic [START_REF] Klukas | Dynamic exploration and editing of KEGG pathway diagrams[END_REF][START_REF] Streit | Interactive visualization of metabolic pathways[END_REF].

In this paper, we present an interactive visualization approach especially designed for bipartite networks [START_REF] Asratian | Bipartite Graphs and Their Applications[END_REF]. These are graphs that allow their set of vertices to be partitioned into two independent sets, so that within each set, there are no adjacent nodes. Bipartite networks occur frequently in biology and other application domains: e.g. as ontologies with concepts and instances, in modeling and simulation as petri nets, and as biochemical reaction networks.

Only recently, their importance has been noticed and first visualization techniques appeared: a well adapted node-link-visualization called Anchored Maps [START_REF] Misue | Drawing bipartite graphs as anchored maps[END_REF] as well as some solutions that have been specifically designed for movieactor-networks as in the InfoVis 2007 contest [START_REF] Kosara R | Information visualization contest 2007[END_REF]. Yet none of these works are concerned with biomedical data and users from this field.

Our visualization approach has been developed within our research training school dIEM oSiRiS (http://diemosiris.de), where biologists, medical researchers, and computer scientists work together to develop modeling and simulation methods as an experimental methodology in biology to achieve new insights into the functioning of biological cell systems. This allows integration of many suggestions and helpful comments from real world users to influence the design process from the very early stages on. During the discussions with the domain experts, two main themes emerged as the most important ones: scalability to large data sets and familiarity to the end users. Because of these two design constraints, we adapted the idea of a tabular display for the node sets. Without cluttering the display, it scales well up to a couple of thousand nodes and even more if fo-cus+context techniques like the Table Lens [START_REF] Rao R | The table lens: Merging graphical and symbolic representations in an interactive focus+context visualization for tabular information[END_REF] are used. Additionally, every user familiar with spreadsheet applications is able to interact with a table naturally and can focus on learning to handle the additional controls for the exploc The Eurographics Association 2008. ration of the graph structure. Furthermore, table-based approaches have recently been used in other domains like visual analytics [START_REF] Stasko | Jigsaw: Supporting investigative analysis through interactive visualization[END_REF] or for the exploration of transition graphs [START_REF] Pretorius | Visual inspection of multivariate graphs[END_REF]. Yet, these approaches are targeting different application domains which pose different challenges and requirements than the life sciences.

We enhance traditional spreadsheet-based visualizations by showing two tables side-by-side, one for each node set in the bipartite graph, which requires additional screen real estate as well as the handling and intuitive communication of two different focus regions within each table. Also, both tables are connected by the lines of the edge set, which inevitably leads to edge crossings and visual ambiguities. When handling the large node and edge sets of real-world bipartite graphs, our technique must be able to minimize the occurring edge clutter and improve the orientation and navigation of the user within such huge networks. This is achieved by combining numerous interaction techniques (clickable edges, clickable selection markers on the scroll bars, etc.) with a new selection mechanism. This allows the user to trigger and to parameterize automated selection scripts, which are able to encapsulate predefined, complex selection logic.

In Section 2, we present our table-based visualization technique for bipartite graphs, as it forms the basis for the complex structural selection techniques presented in Section 3. To underline the usefulness of our visualization technique in combination with the proposed scriptable selections, Section 4 discusses our software implementation in the context of two large datasets: a reconstruction of the human metabolic network [DBJ * 07] from the BiGG Database [SBRG] and a snapshot of the Gene Ontology Database [START_REF] Ontology Consortium | Gene ontology: tool for the unification of biology[END_REF] together with the Homo Sapiens annotation from the Gene Ontology Annotation Database [START_REF] Camon | The gene ontology annotation (GOA) database[END_REF]. The last Section 5 concludes our paper and covers future work.

The Visualization Technique

A bipartite graph G(V 1 ,V 2 , E) consists of two independent node sets V 1 and V 2 and an edge set E connecting both node sets, containing only edges (x, y) with X ∈ V i and y ∈ V j =i . At the very basis, our visualization consists of two tables that represent the two independent node sets. Both tables have columns for the attribute set A i of the respective node set V i that they are displaying. Each row represents one node and displays the node attributes. The connecting edges between the two node sets are represented as lines, which run from a node's row within one table to a row of the other table. Edge weights are mapped onto the width of the individual edges. Additionally, the two 1-mode projections P V 1 and P V 2 are computed by adding all edges (x, z) with x, z ∈ V i , where ∃y ∈ V j =i so that {(x, y), (y, z)} ⊆ E. The number of different intermediate nodes y that connect x and z is used as a weight on the resulting projected edge (x, z). These projections are basically additional edge sets that connect only nodes of either V 1 or V 2 and can be seen as shortcuts that directly show dependencies between nodes of the same set. They are depicted as arcs at the side of the two tables. Apparently, this way of presenting the graph structure is specifically fitted for bipartite graphs, whose edges always connect two nodes from the two different node sets, but never from the same. The overall layout of the described tables, edges, and arcs is shown in Figure 1.

The basic setup of our table-based visualization technique is able to represent all features of a bipartite graph, including derived ones like the projections, in a compact way. Yet, especially large and dense bipartite networks demand some additional features to minimize visual cluttering and to enhance the accessibility of our representation. There are basically three problems to target here:

-Large node sets result in very long tables, which are tedious to browse and navigate. This is especially true for our two tables which are linked by edges, as the user who is investigating the graph's structure is not browsing the tables sequentially, but is rather following edges from one side to another. This results in a lot of back and forth browsing within both tables, which is more time consuming as tables get larger. -Large edge sets result in a lot of lines running in between the two tables, which in turn produces a lot of edge crossings. This makes it harder to follow the course of an individual edge. For this reason, edge crossings are usually to be avoided by any standards of graph drawing aesthetics. -Large attribute sets with dozens of attributes per node, as they occur in real world data, would result in very wide tables with equally dozens of columns. Since we do not want to introduce another scrolling axis by allowing horizontal scrolling, the visualization is limited in its width by the current screen width, whereas it can be arbitrarily extended downwards. The situation gets even worse if node attributes consist of textual descriptions or even images and figures, as it is not uncommon for real world data sets.

Therefore, the following three sections give an overview of additional features which have been included to reduce visual clutter and to maintain the orientation of the user.

Additional Features for Large Node Sets

Focus+Context (A): Using a focus+context technique like the table lens [START_REF] Rao R | The table lens: Merging graphical and symbolic representations in an interactive focus+context visualization for tabular information[END_REF] reduces the height of the table by minimizing all rows that are not part of the current region of interest. In our case, the region of interest within a table is defined by the position of the mouse cursor. The row under the cursor, as well as its neighboring rows, will be zoomed and can then be read and investigated, as demonstrated in the accompanying video. The reduction of the row height in the context area now allows for faster scrolling even in large tables.

Figure 1: This screenshot shows the two node sets as tables, the connecting edges in between, and both 1-mode projections at the sides. The markers point to the special visualization features we added to the basic concept. (A) -Focus+Context in table, (B) -Fisheye scrollbars with selection markers, (C) -Hide unselected rows, (D) -Minimization of edge crossings, (E) -Clickable edges, (F) -URL-references, (G) -Columns for the two different selections, (H) -Maximum level of script, (I) -Highlighting of traversed edges and 1-mode projections.

Fisheye scrollbar with selection markers (B): As rows of interest can also be selected (see Section 3), such selections can span over both of the tables and be scattered all over them. To easily find regions with selected rows in large tables, we have integrated additional selection markers into the scrollbars at the sides. They indicate where selected rows are located in a table. The user can either use the scrollbar to scroll up/down to a selection or directly click on a selection marker to jump instantly to the respective row. Because the selection markers can be placed quite densely and are hard to pinpoint for clicking, we have also added a fisheye lens to the scrollbar. This lens follows the mouse cursor and spreads out the focus area so that, even in crowded regions, individual selection markers can be clicked. A tooltip displays information about the row to which a selection marker belongs. This feature is also shown in detail in Figure 2.

Hide unselected rows (C):

In very large data sets, even with the help of the selection markers, the exploration of a scattered selection can be tiresome. Therefore, we allow the user to reduce the view of the tables to show only selected rows. In this condensed view, unselected rows will be substituted blockwise by a single row that gives information about how many rows have been hidden at that point. An example is given in Figure 2, where only the selected rows of a table are shown.

Additional Features for Large Edge Sets

Minimization of edge crossings (D): As the node sets in the tables can be freely ordered, the edges running in between the tables just follow the ordering of the rows. Thus, for a minimization of edge crossings, at least one of the tables needs to be reordered. A barycentric crossing minimization heuristic [JM97] can be called with just one mouse click to rearrange one of the tables. This immensely reduces the visual clutter in most real-world cases.

Clickable edges (E):

The visual tracking of edges is hampered by their crossings. This makes it hard to discern the c The Eurographics Association 2008. edges and increases the cognitive load for the user. Thus, it is possible to move the mouse over an edge of interest that in turn gets highlighted and becomes clickable if one of its end nodes lies outside of the screen. Additional information on the off-screen endpoint is displayed in a tooltip and a simple mouse click carries the user to the off-screen node.

Additional Features for Large Attribute Sets

URL-references (F): If the set of node attributes is large, it makes sense to divide it into primary, relevant attributes and secondary, supplementary attributes. Primary attributes are shown in the table where the columns can be resized to provide more space in case of lengthy attribute data and tooltips provide quick access to even more details on demand. Secondary attributes on the other hand are outsourced to HTML pages, from where they can be retrieved on demand. A mouse click on the hyperlinked name or label of a node opens up an internet browser window that displays the HTML page (see Figure 3). The choice of the secondary attributes is constrained by the following two requirements:

-Even with the reduced set of primary attributes, it should be possible for the user to identify individual nodes. So, any name or label that is part of the attribute set should be a part of the primary attribute set. -The primary attribute set should still allow to perform the exploration in the desired way: e.g., if a comparison of nodes with respect to some attribute is to be performed, this attribute should be part of the attributes shown in the table.

Structure-based Selection Mechanisms for Bipartite Networks

Especially in large bipartite networks, the ability to derive a subset of the entire data set is important for further analysis or in-depth visualization. Mechanisms for deriving such a subset can be categorized as:

Addition vs. deletion: These are the two dual modes of constructing a subset. One can either start with an empty subset and add the desired data values, or with the complete data set and consecutively delete undesired data values. For interactive selection, often both modes are combined to allow for the refinement of a previously constructed subset. The Toggle-functionality can be seen as the simplest example for that, but also more complex mechanisms have been shown to c The Eurographics Association 2008.

be useful [START_REF] Wills | Selection: 524,288 ways to say "this is interesting[END_REF]. Our structure-based selection mechanism makes use of selection by addition. This is mostly motivated by the fact that it is usually easier to define what kind of data is of interest, instead of defining its inverse.

Automatic vs. interactive:

The subset can either be constructed automatically or be specified interactively by the user. In visualization, selection is generally understood as the interactive counterpart to the more automatic, preprocessing-oriented sampling (constructed by addition) and filtering (constructed by deletion). In this section, we present a 2-step selection mechanism for bipartite graphs. Its two steps divide the selection process into an interactive and an automatic part, which allows to combine the best parts of both worlds: the possibility for the user to influence the selection by interactively triggering and parameterizing a complex automated selection logic.

Binary vs. discrete vs. continuous: In the binary case, a data value either belongs to a subset or not. Whereas in the discrete and continuous cases, each data value is mapped to a Degree of Interest (DoI) between 0 and 1 [START_REF] Doleisch | Smooth brushing for focus+context visualization of simulation data in 3d[END_REF], which determines "how much" a data value belongs to the subset. Discrete and continuous selection methods come naturally with focus+context techniques. Here, the focus region belongs to the selection and is assigned a DoI of 1. Whatever lies outside of this focus region is mapped to a DoI between 0 and 1. When defining such a mapping, there is usually some notion of relatedness to the focus region involved. This relatedness can simply be spatial proximity in data or image space [START_REF] Furnas | Generalized fisheye views[END_REF], or more complex semantic concepts like similarity [START_REF] Novotný | Similarity brushing for exploring multidimensional relations[END_REF] or relevance [SCH * 06]. The selection mechanism presented in this section uses a discrete DoI based on the user's interest in specific parts of the subset.

Combining Interactive and Automatic Selection

Accounting for structural relatedness, selection mechanisms have already been presented for trees and hierarchies [START_REF] Ward | Navigating hierarchies with structure-based brushes[END_REF]. There, they are implemented as a highly interactive brushing mechanism. A number of different parameterizations allows to refine the scope of the selection within the tree. As we now take structural selection further to the more complex case of bipartite graphs, the number of possible configurations in the selection scope becomes even larger. This is because there is more than one way of being related in a bipartite graph: being adjacent through the edge set E, being adjacent through one of the 1-mode-projections P V1 or P V2 , or even not being directly adjacent, but still being connected through a path of two other nodes. Our structural selection method adapts the view of spatial proximity in the graph structure as shown by these examples. This means that we are measuring the distance of a node to a focus node within the graph structure itself and not within their screen representation. Thus, two nodes might have been laid out far apart on the screen, but if they are connected by an edge, their structural proximity is still considered high.

In order to allow the full power of these configurations without reducing the ease of use of the selection, we decided to split the selection process into an interactive and an automatic part. For this, the user selects a set of focus nodes in an interactive first step. This is passed to the automated part, in a second step, which gathers successively all nodes that can be reached by its predefined, set-based selection logic.

1st Step: The interactive selection. In our visualization, we chose a toggling mechanism to mark focus nodes, as it can be seen in Figure 1. Additionally, the user can just click and drag along the selection columns to select entire regions or intervals of rows. Together with the ability to sort the table with respect to any attribute, this provides a brushing mechanism that allows to quickly select items within certain attribute ranges. To enable the comparison of results from two different selections, two columns are added in both tables (marked by (G) in Figure 1), where focus nodes can be picked independently. A unique color is assigned to each of the two individual selections and everything related to them (icons, menu entries, highlighted rows and edges, etc.). By default, red and green are used, as biologists are familiar with this color scheme. If it happens that a node is affected by both selections, the color yellow is used for its respective row in the table, which is inspired by traffic lights, where yellow lies between red and green. However, the colors can be adjusted according to the user's preferences.

2nd

Step: The automatic selection. The result of the interactive selection in the first step is passed to the loaded selection script. Starting with a set that contains only the focus nodes, the scripted selection logic proceeds to add new node sets script line by script line in a breadth first manner by traversing along edges or projections. The maximum number of selection levels is determined by the number of lines within the selection script. It can be interactively lowered from within the visualization for both available scripts, as shown in Figure 1 at marker (H). Our selection logic allows to define a DoI-value in every script line, which is then assigned to the node set added to the selection by the script line. That way, the user's interest in the result of individual script lines can be specified more precisely compared to an automatically derived DoI-value, for example by proximity to the focus nodes. The values lie between 0 and 100, where 0 means no accentuation or amplification of a row, and 100 means maximum accentuation. The value ranges have been chosen because we found that the "percentage-thinking" is intuitively understood throughout all application domains. In our table-based visualization, the DoI-values define how the nodes of the set should appear in the tables: it is used to define the color saturation and the height of the respective row. An example of such a selection script is given in the context of the examples discussed in Section 4.

To make our concept more flexible, we extended it to allow the filtering of node sets by the nodes' attributes values. These can be attributes from the dataset itself or some c The Eurographics Association 2008.

of the commonly derived structural attributes, like indegree, outdegree or the clustering coefficient, which is used in its adapted version for bipartite graphs [START_REF] Robins | Small worlds among interlocking directors: Network structure and distance in bipartite graphs[END_REF]. The filtering of node sets can be used to narrow down the traversed paths to only those paths of interest. The paths along which the script traverses the graph are highlighted in the visual representation of the bipartite graph (marker (I) in Figure 1), the edges along the paths are colored with the same saturation as the nodes they emanate from. This is very helpful for understanding the inner workings of selection scripts, as one is not only given the script result (nodes with DoI-values attached), but also some information on how this result came about. It ties in very nicely with the possibility to interactively define the maximum level up to which a script is computed.

Discussion

This section will give an overview of the usage of the visualization and selection technique by discussing two real world data sets. The first data set, a biomedical ontology, is foremost used to describe the visualization technique, and the second data set, a biochemical reaction network, is then used to illustrate the selection mechanism.

A Biomedical Ontology

Ontologies are a collection of interrelated concepts or terms that aim at providing a basis for knowledge management in a certain area. Data sets can use the provided terms, to declare their data instances of these concepts, annotating them and thus providing semantics alongside the pure numbers. These annotations form a bipartite graph, with the terms as one node set and the instances as the other. This bipartite model is well known from research on Semantic Web technologies like folksonomies or lightweight ontologies.

One of the largest biomedical ontologies freely available today is the Gene Ontology Database [START_REF] Ontology Consortium | Gene ontology: tool for the unification of biology[END_REF] (OBO snapshot from 12-JUN-2008), which we used together with one of the largest annotated data set of gene products, the Homo Sapiens annotation from the Gene Ontology Annotation Database [START_REF] Camon | The gene ontology annotation (GOA) database[END_REF] (GOA snapshot from 07-JUN-2008). This data set consists of 26389 terms from the ontology and 35043 gene products annotated with these terms. These annotations are encoded in 377132 undirected links in between these two sets. By projecting these annotations onto the set of terms, 122379 additional edges are produced. The projection onto the gene products ran into a combinatorial explosion and produced literally billions of additional edges. Since this number is well beyond the scope of our tool, we filtered this huge set of projected edges down to the 291586 edges, with weight 16 or higher. The weight of 16 is basically an arbitrary choice, as it presents the lowest filtering threshold for this data set that yields an edge set with less than a million edges, which is more or less the limit of our implementation. The data sets also have a variety of attributes, as short natural language definitions for the terms and bibliographical references for the gene products. To make all of this information available, the terms have been hyperlinked to the Gene Ontology Database and the gene products to the five databases they originally stem from: UniProtKB, RefSeq, ENSEMBL, H-invDB, and VEGA. So, this is a good example of how our tool ties together individual biological databases.

When exploring this rather large data set with our tool, it soon becomes apparent that scrolling through the tables is a very time consuming and slow method to access the network. It actually takes more than 700 turns on the mouse wheel to scroll from the top of the terms table to its bottom. So, this is exactly one of the cases, where all the other means of rapid navigation come into play: clickable edges, fisheye scrollbar and folding of unselected nodes to compress the view. As shown in Figure 2, the folded view of the terms table compresses several thousand unselected rows into onerow-placeholders. Navigating in the compressed view and only switching back to the full view if a new selection is to be made, speeds up the exploration process immensely.

The ontology data set exhibits an interesting way of using the projections: as indicators of similarity. That is because the number of shared terms is often taken as a measure of similarity for two instances. Also, the more instances two terms have in common, the more alike they probably are. This notion is captured exactly by the projections and their weights: the higher the weight of a projected edge, the more disjoint paths exist between its two incident nodes. And the more disjoint paths there are, the more joint intermediary nodes must lie within the other node set. Hence, clicking and thus following the projections during the exploration makes it very easy to find similar nodes within this large data set -something that would not have been possible without the projections on the side.

A Biochemical Reaction Network

Reaction networks are usually described as hyper graphs, in which each directed hyper edge encodes a biochemical reaction connecting its reactants with its products. We use a common transformation of hyper graphs into the so-called König's representation [START_REF] Temkin | Chemical Reaction Networks: A Graph-Theoretical Approach[END_REF] to map the reaction network to a bipartite graph. Hereby, the hyper edges (reactions) are transformed into nodes themselves and new directed edges are introduced from each reactant to its reaction and from each reaction to its products. This transformation yields two node sets: the chemical substances and the reactions.

As one of the largest biochemical reaction networks that has been reconstructed so far, we used the human metabolic network [DBJ * 07] from the BiGG Database [SBRG] (SBML snapshot from 20-DEC-2007). It comprises 2764 chemical compounds, 3311 chemical reactions and 17519 directed links with stoichiometric factors encoded as weights between compounds and reactions. Additionally, the projecc The Eurographics Association 2008. tion on the compounds yields 121118 edges and the projection on the reactions 295922 edges. The data set includes an abundance of attributes, from charges for the compounds to bibliographical references for most of the reactions. As not all of them could be displayed, the names of the compounds and reactions are hyperlinked to the BiGG database with all the additional information, as seen in Figure 3 To illustrate the script-based selection mechanism, we derived a selection logic for a dependency analysis, as this kind of analysis is often performed by biologists on reaction networks. This means, for a focussed chemical compound the selection script should automatically determine on which other compounds it directly depends. As there are no edges between the nodes of one node set, direct dependence means that there exists at least one directed path of length 2 between them. So, such a script must follow down all reactions that produce this very compound and select their reactants. This selection logic needs to be transformed into a set-based notation, which could look like the following, with

V = V 1 ∪V 2 : Level 0 (DoI = 100%) = focus nodes Level 1 (DoI = 0%) = {u ∈ V : ∃(u, v) ∈ E : v ∈ Level 0 } Level 2 (DoI = 75%) = {u ∈ V : ∃(u, v) ∈ E : v ∈ Level 1 }
Here, with each step, the logic is traversing backward edges, starting from the focussed compound and traversing backwards twice: once to reach all reactions that produce the compound and another time to reach their reactants. Thereby, it is guaranteed to return in Level 2 only nodes that are also chemical compounds. The example also shows that the DoI is not necessarily decreasing with each level. Instead, each level can be assigned its individual DoI according to the user's goal of analysis. In our example, we are only interested in the compounds, but not in the intermediate reactions and assign them therefore a DoI of 0. Given the intended logic in a set-based notation, it can now be written as a selection script, which is shown in Figure 4. Once defined, this script can be loaded and used with a mouse click on whichever compound the biologist likes to run the dependency analysis on. This is the interactive part of the selection, which toggles nodes to build up a set of focus nodes to be passed to the script for dependency analysis.

It would also be possible to shorten the shown selection script by traversing the backward 1-projection on the compounds in just one step, instead of going two steps along the edges themselves as before. This illustrates nicely what we meant when we introduced the 1-mode-projections as shortcuts within the bipartite graph. It is, of course, possible to create much more complex selection scripts, as we have done for the example in Figure 1, which shows a part of the discussed metabolic network.

Conclusion and Future Work

We have presented a compact visualization technique for the interactive exploration of large bipartite biomedical networks that adapts the intuitive and familiar tabular display to the characteristics of bipartite graphs. Displaying the two node sets in separate tables offers a natural way to lay out the bipartite graph and to combine structural information (edges and 1-mode projections) and multiple attributes in one display. Furthermore, the tabular display is scalable to very large data sets and provides features to easily rearrange the node sets based on their attribute values. These are major advantages of our concept in comparison to standard nodelink graph layouts with their cumbersome node placement and label arrangement. A rich set of interaction techniques allows for a convenient exploration of the bipartite graph, e.g. lens functions to provide focus+context and functions to jump directly to interesting parts of the data outside of the display. Another important aspect in this regard is the combination of interactive and automated selection of nodes: an interactively defined selection is extended by automatic, script-based processing. While the technique in its whole is tailored to the concrete application domain, the underlying concepts, e.g. the multi-tabular display, the edge-based navigation, or the selection scripts can be applied to the general area of visualizing multipartite graphs.

The applicability of our technique to real world data was ensured by including the feedback of domain experts from biology and medical sciences into the design process. Their overall impression of the resulting visualization and interaction techniques, as we have presented them in this paper, was mostly positive. They considered the tabular display an intuitive way to analyze the data and found the 2-step structural selection useful to find features of interest after some time of getting used to its concept.

In addition to the features already present, they would have wished for an even tighter integration with the databases in the sense that navigation through the database interfaces links back and adapts the visualization likewise. This would give them all the power of the visualization while still being able to use the textual query interfaces of the databases as they are used to do. Such a "hybrid" approach will probably also further the acceptance of our tool, as it would result in a more moderate learning curve.

As for enhancements to the visualization itself, we plan to improve the navigation along edges as this has turned out to be the preferred way of the biologists to traverse the network. Here, the problem is, that edges are often collinear and c The Eurographics Association 2008.

cluttered, which makes it hard to pick and click a specific edge for traversal. This effect can be reduced by adding edge lenses to locally declutter the view [START_REF] Wong | Edgelens: An interactive method for managing edge congestion in graphs[END_REF][START_REF] Tominski | Fisheye treeviews and lenses for graph visualization[END_REF] and allow a precise selection of individual edges.

Another hurdle to overcome is the initial skepticism towards the selection scripts by new users, as their concept seems very abstract at first. Yet, we experienced that after a few minutes of training, this attitude changes. Therefore, we imagine to construct a whole library of different selection scripts and to investigate mechanisms that provide a suitable subset of this library for different exploration goals and data domains. These scripts would provide a good starting point for new users, who could readily use them without having to worry about their inner workings. As scripts can be shared, it is also possible for more advanced users to encapsulate their selection logic and pass their self-made scripts on to others.

Figure 2 :

 2 Figure2: A table in which unselected rows are collapsed into placeholder rows. The example also illustrates the possibility to adapt the color scheme to the user's personal preferences.

Figure 3 :

 3 Figure 3: Linked HTML page -the hyperlinked name of a focussed node was clicked such that a browser starts showing an HTML page with secondary attributes.

Figure 4 :

 4 Figure 4: An example of a selection script.

c The Eurographics Association 2008.

Acknowledgements

The authors wish to thank Steffen Hadlak for his work on the implementation and the video presentation. This work is supported by the DFG graduate school dIEM oSiRiS.


Yang Chen
email: chenwf@fimmu.com

Bayesian Image Restoration Using A Large-scale Total Patch Variation Prior

Keywords: Bayesian restoration, Large-scale Total Patch Variation (LS-TPV) Prior, Parallelization, Convergence

. However, such SA prior restoration suffers from high computational cost and the unguaranteed convergence. Concerning these issues, this paper proposes a Large-scale Total Patch Variation (LS-TPV) Prior model for Bayesian image restorations. In this model, the prior for each pixel is defined as a singleton conditional probability, which is in a mixture prior form of one patch similarity prior and one weight entropy prior. A joint MAP estimation is thus built to ensure the iteration monotonicity. The intensive calculation of patch distances is greatly alleviated by the parallelization of Compute Unified Device Architecture(CUDA). Experiments with both simulated and real data validate the good performance of the proposed restoration.

Introduction

Image restorations have wide applications in remote sensing, radar imaging, tomographic imaging, microscopic imaging, astronomic imaging, digital photography, etc [START_REF] Srivastava | On advances in statistical modeling of natural images[END_REF][START_REF] Tan | A unified iterative denoising algorithm based on natural image statistical models: derivation and examples[END_REF][START_REF] Molina | On the hierarchical Bayesian approach to image restoration. Applications to astronomical images[END_REF][START_REF] Archer | On some Bayesian/regularization methods for image restoration[END_REF]. For linear and shift invariant imaging systems, the transformation from f to g is well described by the following additive linear degradation model [START_REF] Molina | On the hierarchical Bayesian approach to image restoration. Applications to astronomical images[END_REF][START_REF] Archer | On some Bayesian/regularization methods for image restoration[END_REF]:

g = A * f +  (1)
where g, f and  represent the degraded observed image, the original true image and the corrupting white Gaussian noise with variance  2 , respectively. A is the PSF (Point Spread function) of the imaging system and * is the linear convolution operator. Throughout this paper, we assume that the PSF A is known or could be numerically determined or estimated. Bayesian or Maximum a posteriori (MAP) approach, within Markov Random Fields (MRFs) framework, can provide a stable solution of the ill-posed inverse image restorations through incorporating a priori information about the geometrical properties of an image [START_REF] Molina | On the hierarchical Bayesian approach to image restoration. Applications to astronomical images[END_REF][START_REF] Archer | On some Bayesian/regularization methods for image restoration[END_REF][START_REF] Geman | Stochastic relaxation, Gibbs distributions and the Bayesian restoration of images[END_REF][START_REF] Galatsanos | Methods for choosing the regularization parameter and estimating the noise variance in image restoration and their relation[END_REF][START_REF] Bioucas-Dias | A New TwIST: Two-Step Iterative Shrinkage/ Thresholding Algorithms for Image Restoration[END_REF][START_REF] Wan | A Bayesian multiscale approach to joint image restoration and edge detection[END_REF]. Via modelling the unknown parameters in the prior probability density functions, prior information can also be interpreted as the regularization metrics, which measures the extent to which the contextual constraint assumption is violated. With the incorporated prior information, Bayesian approach is able to distinguish solutions from less desirable ones by transforming the original ill-posed problem into a well-posed one. Generally, we can build the following posterior probability P ( f |g) for image restorations.

P ( f |g)  P (g| f ) P ( f ) (2)
P (g| f ) = exp ( L(g, f )) = exp - 1 2 g -A * f 2 (3)
P ( f ) = Z -1 × exp(-U ( f )) = Z -1 × exp -  j U j ( f ) (4)
Here, P (g| f ) and P( f ) denote the likelihood distribution and the prior distribution, respectively. L(g, f ) and U( f ) denote the corresponding likelihood energy and the prior energy function. U j ( f ) is the notation of the energy function U evaluated at pixel index j. The partition function Z is a normalizing constant.  is the global parameter controlling the balance between the likelihood energy and the prior energy. We can build the posterior energy function in the logarithm form:

  ( f ) = log P ( f |g) = L(g, f ) -U( f ) = - 1 2 g -A * f 2 -  j U j ( f ) (5)
The restored or deconvolved image f can be obtained through maximization of function   ( f ).

The widely used quadratic membrane (QM) prior or the Tikhonov L2 regularization tend to smooth both noise and edge details, and often lead to unfavorable oversmoothing in restorations [START_REF] Geman | Stochastic relaxation, Gibbs distributions and the Bayesian restoration of images[END_REF]. On the other side, some edge-preserving Bayesian restoration methods have been proposed, which can be classified into three main categories: wavelet relevant regularization restorations, Bayesian restorations with nonquadratic prior energies and Bayesian restorations with total variation(TV) prior. Wavelet relevant regularization restorations, with multi-scale stochastic prior models, have been proposed in the research of edge-preserving restorations [START_REF] Wan | A Bayesian multiscale approach to joint image restoration and edge detection[END_REF][START_REF] Figueiredo | An EM algorithm for wavelet based image restoration[END_REF][START_REF] Bioucas-Dias | Bayesian wavelet-based image deconvolution: A GEM algorithm exploiting a class of heavytailed priors[END_REF][START_REF] Neelamani | Forward: Fourier-wavelet regularized deconvolution for ill-conditioned systems[END_REF]. In [START_REF] Bioucas-Dias | Bayesian wavelet-based image deconvolution: A GEM algorithm exploiting a class of heavytailed priors[END_REF], priors based on wavelet decompositions and heavy-tailed pdfs were incorporated into the EM restoration algorithm. In [START_REF] Neelamani | Forward: Fourier-wavelet regularized deconvolution for ill-conditioned systems[END_REF], R. Neelamani et al. proposed a Fourier-Wavelet hybrid restoration, which exploits the potentials of the sparse representation of noise for Fourier transform and the sparse representation of the coefficients in wavelet transform. All these wavelet regularization methods are based on the distribution modeling and thresholding of the decomposed wavelet coefficients [START_REF] Neelamani | Forward: Fourier-wavelet regularized deconvolution for ill-conditioned systems[END_REF]. Edge-preserving priors with nonquadratic energies were also proposed to preserve the edge details by using nonlinear inverse-proportional functions between edge existences and intensity differences [START_REF] Geman | Stochastic relaxation, Gibbs distributions and the Bayesian restoration of images[END_REF][START_REF] Bouman | A generalized Gaussian image model for edge-preserving MAP estimation[END_REF][START_REF] Yu | Edge-Preserving Tomographic Reconstruction with Nonlocal Regularization[END_REF]. The weighting matrices in nonquadratic prior Bayesian restorations preserve the edges by turning off or suppressing smoothing at appropriate locations. Another impressive research direction in this area is the total variation (TV) based image restorations [START_REF] Bioucas-Dias | A New TwIST: Two-Step Iterative Shrinkage/ Thresholding Algorithms for Image Restoration[END_REF][START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF][START_REF] Chambolle | An algorithm for total variation minimization and applications[END_REF]. This kind of approach can also be viewed as TV priors with the half-quadratic regularization prior energies. With no enforced global image continuity, TV prior restorations often demonstrate good property of edge-preserving.

The nonquadratic priors preserve the edges by determining the pixel-wise regularizations based on the intensity-difference information within local fixed neighborhoods. The regularizations on TV priors also come from the information of the local intensity gradients. For all these nonquadratic and TV priors, the local information of intensity difference can by no means provide effective regularizations to discriminate noise from image textures in different scales. So these edge-preserving approaches often fail to cope well with some complicated restorations with high noise contaminations, and tend to produce block artifacts around the continuous edges in those situations [START_REF] Mignotte | A segmentation-based regularization term for image deconvolution[END_REF][START_REF] Chandas | Bayesian Restoration Using a New Non-Stationary Edge-Preserving Image Prior[END_REF][START_REF] Chandas | Bayesian Image Restoration Based on Variational Inference and a Product of Student-t Priors[END_REF][START_REF] Chen | Nonlocal prior Bayesian tomographic reconstruction[END_REF][START_REF] Buades | A nonlocal algorithm for image denoising[END_REF][START_REF] Chen | Effective Image Restorations Using a Novel Spatial Adaptive Prior[END_REF].

Objects in the images are reasonably assumed to be a composition of edges and backgrounds with intensities changing coherently over different scales. Thus, for the traditional prior model, the intensity differences between individual pixels within local neighborhoods are often insufficient to characterize objects [START_REF] Chen | Nonlocal prior Bayesian tomographic reconstruction[END_REF][START_REF] Buades | A nonlocal algorithm for image denoising[END_REF]. This limit elicited our previous work in [START_REF] Chen | Effective Image Restorations Using a Novel Spatial Adaptive Prior[END_REF] which introduced a spatial adaptive (SA) prior restoration approach. This SA prior works by adaptively including the relevant neighboring pixels and excluding the negative irrelevant ones within a large neighborhood. Though showing expressive results, the iterative restoration in [START_REF] Chen | Effective Image Restorations Using a Novel Spatial Adaptive Prior[END_REF] updates the prior weight using the latest image estimate in a One Step Late (OSL) way, which leads to inconsistent posterior energy function and non-convergent iterations. In addition, the approach in [START_REF] Chen | Effective Image Restorations Using a Novel Spatial Adaptive Prior[END_REF] is greatly limited by the high computational cost in calculating all the patch distances for all the pixel pairs in each prior neighborhood. In this article, to overcome the first problem, we applied a Large-scale Total Patch Variation (LS-TPV) Prior model, which is built through the specification of the conditional probability for each pixel [START_REF] Lyu | An implicit Markov random field model for the multi-scale oriented representations of natural images[END_REF][START_REF] Winkler | Image Analysis, Random Fields And Markov Chain Monte Carlo Methods[END_REF][START_REF] Peyr | Non-local regularization of inverse problems[END_REF][START_REF] Zhu | FRAME: Filters, Random field And Maximum Entropy:Towards a Unified Theory for Texture Modeling[END_REF]. The joint maximization estimation of this constrained entropy function will lead to convergent image/weight restorations. To overcome the second problem, we introduce the parallel computation structure of CUDA to calculate the patch distances [START_REF]NVIDIA CUDA Programming Guide[END_REF][START_REF] Fatica | Accelerating matlab with cuda[END_REF][START_REF] Chen | Improving Low-dose X-ray CT Images by Weighted Intensity Averaging over Large-scale Neighborhood[END_REF]. Furthermore, another benefit from the CUDAaccelerated parallelization is that it allows further performance enhancement for the LS-TPV prior restoration by using larger search neighborhoods.

In section 2, a review of some previous prior models is illustrated, and after that we introduce the proposed LS-TPV prior model. In section 3, we give a joint estimation algorithm for the proposed restoration. In section 4, we perform comparative experiments with both simulated and real data. Relevant visual/quantitative results are also presented. Conclusion and discussion are given in section 5.

Prior Model

Conventionally, the value of U j ( f ) is commonly computed through a weighted sum of potential functions v of the differences between the pixels in the neighborhood N j :

U j ( f ) =  b∈N j w b j v ( f b -f j ) (6)
Generally, different choices of the potential function v lead to different priors. For the simple space-invariant QM prior, the potential function has the form v(t) = t 2 /2. Some Edgepreserving nonquadratic priors could be chosen by adopting a nonquadratic potential function v, such as the Huber potential function:

v(t) = t 2 2, |t| ≤   |t| - 2 2, |t| >  ( 7
)
where  is the threshold parameter [START_REF] Geman | Stochastic relaxation, Gibbs distributions and the Bayesian restoration of images[END_REF][START_REF] Yu | Edge-Preserving Tomographic Reconstruction with Nonlocal Regularization[END_REF]. Such edge-preserving nonquadratic priors preserve structure information by choosing nonquadratic potential functions that increase less as the differences between the adjacent pixels become bigger. Weight w b j is a positive value that denotes the interaction degree between pixel b and pixel j. In traditional prior models, it is simply considered to be inversely proportional to the distance between pixels b and j. So on a square lattice of image f , in which the 2D positions of pixels b and j (b = j) are respectively (b x , b y ) and ( j x , j y ), w b j is usually calculated by the geometric distance 1

(b x -j x ) 2 + (b y -j y ) 2 .
Restorations using Total Variation (TV) prior take the prior energy as the following expression:

U TV ( f ) =  j  h j f 2 +  v j f 2 (8)
where  h j and  v j are linear operators corresponding to the first-order horizontal and vertical intensity differences at pixel j, respectively [START_REF] Bioucas-Dias | A New TwIST: Two-Step Iterative Shrinkage/ Thresholding Algorithms for Image Restoration[END_REF][START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF][START_REF] Chambolle | An algorithm for total variation minimization and applications[END_REF]. Since both the smooth and sharp edges tend to have similar prior energies, this TV prior energy given by ( 8) does not penalize weak discontinuities, and often favors images with bounded variations.

Traditional quadratic priors, nonquadratic priors and TV prior only provide local intensity information for Bayesian restorations. And the edge-preserving nonquadratic and TV priors, though being able to preserve edge information, tend to produce unfavorable block artifacts or false edges in high noise situations.

To improve restoration quality, in [START_REF] Chen | Effective Image Restorations Using a Novel Spatial Adaptive Prior[END_REF] we proposed a Spatial Adaptive (SA) prior model, in which a spatial adaptive prior neighborhood is implemented by setting the weight w b j to an 1-0 binary function to classify the pixels in neighborhood N into the pixels relevant to the center pixel j and those not. The building of the SA prior can be formalized as follows:

U SA ( f ) =  j U j ( f ) =  j  b∈N j w b j ( f b -f j ) 2 Nr j ( 9
)
w b j = 1 |n b ( f ) -n j ( f )| 2 <=  0 |n b ( f ) -n j ( f )| 2 >  ( 10
)
|n b ( f ) -n j ( f )| 2 =  l f l∈n b -f l∈n j 2 (11)
Here, U SA is the energy function for the SA prior. w b j , which represents the classification of the neighboring pixels in the search neighborhood N j , and can be computed via equations ( 10)- [START_REF] Neelamani | Forward: Fourier-wavelet regularized deconvolution for ill-conditioned systems[END_REF]. Nr j , the number of the neighboring pixels with nonzero w b j in the neighborhood N j , is the normalization factor for the different sizes of non-1 neighboring pixels in each N j . Parameter  in [START_REF] Bioucas-Dias | Bayesian wavelet-based image deconvolution: A GEM algorithm exploiting a class of heavytailed priors[END_REF] is the threshold parameter. The distance |n b ( f )n j ( f )| 2 is determined by the distance measure [START_REF] Neelamani | Forward: Fourier-wavelet regularized deconvolution for ill-conditioned systems[END_REF] between the two translated comparing patches n b and n j , whose centers are respectively located at pixel b and pixel j. Though proved to be effective in suppressing both noise and block artifacts, all the weights w in SA prior need to be heuristically inferred using current image estimate in the way of one step late (OSL), which often leads to inconsistent posterior energy function with no guaranteed convergence [START_REF] Wan | A Bayesian multiscale approach to joint image restoration and edge detection[END_REF]. To overcome this problem, in this paper we proposed a LS-TPV prior model, in which regularization takes effect through penalizing the weighted distances/variations between the patches surrounding each neighboring pixel pair. The weights w in each N are considered variables with a nonlinear dependence on image f . For the central pixel j and its neighboring pixel b in each prior neighborhood N j , w b j is considered the similarity probability between the two patches n b and n j surrounding pixels b and j, respectively [START_REF] Peyr | Non-local regularization of inverse problems[END_REF]. A entropy prior of all the weights w b j for each N j is introduced into our model in the form of - j  b∈N j w b j lnw b j . Here as the "Frame" model in [START_REF] Zhu | FRAME: Filters, Random field And Maximum Entropy:Towards a Unified Theory for Texture Modeling[END_REF], the constraint is  b∈N j w b j = 1. Then, based on joint MAP theorem [START_REF] Bouman | A generalized Gaussian image model for edge-preserving MAP estimation[END_REF], we can obtain

P( f , w, |g)  P(g| f )P( f |w, )P(w, ) (12)
Here,  is the introduced pixel-wise Lagrange-multiplier parameter. For each pixel j, with f denoting the current image estimate, the prior energy U LS-T PV ( f j , w, | f ) of the LS-TPV prior is in the combined form of a weighted patch similarity energy, a entropy energy -

  ( f , w, |g) = log P(g| f ) + logP( f |w, ) + log P(w , ) = L(g, f ) - (  j U LS-T PV ( f j , w, | f )) = - 1 2 ||g -A * f || 2 -  j ((  b∈N j w b j d b j + h  b∈N j w b j ln w b j ) +  j (  b∈N j w b j -1)) (13)
Here, as (4),  is the global parameter controlling the balance between the likelihood energy and the prior energy. The patch distance d b j is calculated as 2 , in which G  denotes the Gaussian convolution kernel with standard deviation .  j is the value of Lagrangemultiplier parameter at pixel j. The parameter h controls the maximum entropy constraint on w, which is routinely modulated with respect to the noise levels (large h is routinely set to suppress more noise for high noise situations).

||G  * n b -G  * n j ||

Restoration Algorithm

To obtain tractable maximization of the posterior probability using the proposed LS-TPV prior energy U LS-T PV in (13), we use the algorithm of Iterative Coordinate Descent (ICD), in which each element is sequentially updated by maximizing the log-style local conditional probability function [START_REF] Bouman | A unified approach statistical tomography using coordinate descent optimization[END_REF]. In this way the objective image f and the weight w are jointly estimated through alternative maximization of the   ( f , w, |g) with respect to f and w: 1) arg max w -U LS-T PV (w, | f ): f denotes the fixed current image estimation, and each weight w b j for each pixel j is updated by solving the following singleton optimization problem: arg max

w b j (- b∈N j w b j ||G  * n b -G  * n j || 2 -h  b∈N j w b j ln w b j - j (  b∈N j w b j -1)) (14)
With Lagrange multiplier  j introduced, by joint solving [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF] with respect to each w b j and  j , we can then obtain w b j :

w b j = exp(-G  * n b -G  * n j 2 /h)  b ∈N j exp(-G  * n b -G  * n j 2 /h) ( 15
)
where, G  * n b -G  * n j 2 can be easily calculated using currently available f . Each w b j is determined not only by the similarity metrics between patches n b and n j , but also by the sum of all the similarity metrics between patches with the central pixels located in N j . So w b j is different from w jb .

2) arg max f   ( f , ŵ, |g): Here ŵ denotes the fixed current weight estimate. Omitting the terms with no f , the   ( f , ŵ|g) becomes

  ( f , ŵ|g) = L(g, f ) -  j (  b∈N j ŵbj ||G  * n b -G  * n j || 2 )) (16)
To make a tractable maximization of ( 13), based on the theories of pseudo-likelihood and Iterated Conditional Mode (ICM) in [START_REF] Bouman | A unified approach statistical tomography using coordinate descent optimization[END_REF] and [32], we can restore each f j by solving the following factorized problem:

argmax f j   ( f j , ŵbj |g, f ) ⇒ arg max f j L ( f , g) - (  b∈N j ŵbj ||G  * n b -G  * n j || 2 )) (17)
Considering the convexity for ||g - 2 , we conclude that, given the fixed weight ŵ, the Hessian matrix of each   ( f j , ŵbj |g, f ) in ( 17) is a negative definite. In this step of image updating, we applied the restoration algorithm in [START_REF] Archer | On some Bayesian/regularization methods for image restoration[END_REF] to solve [START_REF] Chandas | Bayesian Restoration Using a New Non-Stationary Edge-Preserving Image Prior[END_REF] to obtain the image estimate of each pixel f j . The   ( f , w, |g) in ( 13) is separately convex with respect to each f and w. Denoting f n , ŵn and f n+1 , ŵn+1 as the n th and n + 1 th iterated estimates, we can obtain following relation:   ( f n+1 , ŵn+1 |g) >=   ( f n , ŵn+1 |g) >=   ( f n , ŵn |g), which implies convergence to one local maximum can be obtained for this joint updating strategy [START_REF] Chen | Adaptively regularized constrained total least-squares Image Restoration[END_REF].

A * f || 2 and ||G  * n b -G  * n j ||
In each iteration in the proposed restoration, it is very expensive to perform the pixel-wise calculations of all the d b j and w b j for each pair of translated patches n in each N over the whole image region. In the experiments, to save computation cost, we do not extend the sizes of N over the whole image region, and set them to some appropriate sizes. We achieve significant acceleration by using parallel Compute Unified Device Architecture(CUDA) [START_REF]NVIDIA CUDA Programming Guide[END_REF][START_REF] Fatica | Accelerating matlab with cuda[END_REF][START_REF] Chen | Improving Low-dose X-ray CT Images by Weighted Intensity Averaging over Large-scale Neighborhood[END_REF]. For each iteration in the restoration using the proposed LS-TPV prior model, all threads in this block-grid structure execute simultaneously to perform all the involved pixel-wise operations. Furthermore, considering the symmetry property that distance d b j equals to distance d jb , we can further save one half computation cost by only calculating one of the two distances d b j and d jb .

EXPERIMENTS

Both simulated and real data are used in the experiments. To provide comparative results, we also performed Wiener deconvolution, ForWaRD restoration, Bayesian restoration using TV prior, Huber prior and the SA prior in [START_REF] Chen | Effective Image Restorations Using a Novel Spatial Adaptive Prior[END_REF]. The Wiener deconvolution was performed by using the degraded image to estimate the power spectrum in Fourier domain. The TWIST algorithm in [START_REF] Bioucas-Dias | A New TwIST: Two-Step Iterative Shrinkage/ Thresholding Algorithms for Image Restoration[END_REF] and Newton-Raphson algorithm in [START_REF] Bouman | A generalized Gaussian image model for edge-preserving MAP estimation[END_REF] are used in the restorations using TV prior and Huber prior, respectively. Restored images from above Wiener filter method are taken as the initial images for all the Bayesian restorations. In experiment, the parameters needed to be preset include the  for all Bayesian restorations, the threshold  for Huber prior, the wavelet threshold and decomposition level in the ForWaRD method, the  for the SA prior, the h and the n and N settings for the proposed LS-TPV prior. In the simulation experiment, all the parameters are set based on SNR maximization criteria with SNR calculated as:

SNR = 10 log 10 ⎛ ⎜ ⎜ ⎜ ⎝ M  j (F j -F) 2 M  j F j -f j 2 ⎞ ⎟ ⎟ ⎟ ⎠ (18)
Here M, f , F, and F denote the total pixel number in image, the restored image, the original true image and the mean of the original true image, respectively. As to the real data experiment, which has no available true image for SNR calculations, all the parameters are set based on the tradeoff between edge preservation and noise suppression. N and n of sizes (11×11 N and 7×7 n) is used in our work. A full study of the effects of the sizes of N and n on the resulting restorations is performed in the below Section 4.5. Generally, the solution to the restoration is considered optimized when the iteration goes stable. So for the Bayesian restorations using TV prior, Huber prior and SA prior in which the energy functions are maximized with respect only to image, we stop the iterations when the current image estimate satisfy the condition:

  ( f n+1 |g)-  ( f n |g) < 0.995(  ( f n |g)-   ( f n-1 |g)).
For the restorations using the proposed LS-TPV prior in which the energy functions are maximized with respect to both image and weights, we stop the iterations when the current image and weight estimates satisfy the condition:

  ( f n+1 , ŵn+1 |g)-  ( f n , ŵn |g) < 0.995(  ( f n , ŵn |g)-  ( f n-1 , ŵn-1 |g)).
In practical experiments, to reach above iteration conditions, it is found that 411, 433, 475 and 486 iterations are required for the Bayesian restorations using TV prior, Huber prior, SA prior and the proposed LS-TPV prior, respectively.

Simulation Experiment

In this section, simulated experiments with 256 × 256 "Lena" image (intensity range [-98 117]) are performed. In the simulation, a uniform 9 × 9 PSF is assumed, and Gaussian noise with variance = 0.1663 is added. The PSF used in this experiment is normalized to 1. Fig. 1(b) is the degraded image simulated by (1) using the PSF with the white Gaussian noise. The parameters for different restorations are listed in Table 1. All the restored images are illustrated in Fig. 1(c)-(h), respectively. The restored images from TV prior restoration and Huber prior restoration are respectively shown in Fig. 1(d) and Fig. 1(e). Fig. 1(f) shows the restored image using the ForWaRD method. The result of Bayesian restoration using the SA prior is illustrated in Fig. 1(g). For the Bayesian restoration using the proposed LS-TPV prior, the algorithm in above Section 3 is used, and the corresponding restored image is shown in Fig. 1(h).

Through the results, we find the ForWaRD method and all the Bayesian restorations (Fig. 1(d)-(h)) can overcome the ring effects in Wiener deconvolution (Fig. 1(c)). The block artifacts (the left arrows), which are observed in the ForWaRD approaches and the restorations using Huber prior and TV prior, are effectively suppressed using the restorations using the SA prior and the proposed LS-TPV prior. The LS-TPV prior (Fig. 1(h)) presents a further visual enhancement over the SA prior (Fig. 1(g)) in block-suppressing (the left arrows) and edgepreserving (the right arrows). Fig. 2 plots the the profiles along one specified track (the red line in the hair region in "Lena" image in Fig. 2) for different restorations, and we can see that the profile from the proposed restoration (the blue profile) has the closest match with the profile of the true "Lena" image (the red profile) than those from other restorations. Table 2 shows that, compared to other methods, the proposed LS-TPV prior restoration can lead to restored images with a higher SNR. (g)SA prior restoration [START_REF] Chen | Effective Image Restorations Using a Novel Spatial Adaptive Prior[END_REF] (h) LS-TPV prior restoration

Real Data Experiment

In this section, a real blurred astronomical 256 × 256 "moon" image with intensity range [-77 170] is obtained from 'http://www.iceinspace.com.au/ '. Here we use a 5 × 5 normalized Gaussian PSF, which was tested to be effective. The parameters for different restoration approaches are listed in Table 3. 3(d). Fig. 3(e) and Fig. 3(f) show the restored images from the ForWaRD method and the Bayesian restoration using SA prior. As to the Bayesian restoration using the LS-TPV prior, the restored image is shown in Fig. 3(g). Also, we can note that the ForWaRD method and all the Bayesian restorations (Fig. 3(c)-(g))are free of the ring effects in Wiener deconvolution (Fig. 3(b)). Also, the block artifacts observed in the restorations using Huber prior and TV prior are effectively suppressed in the restorations using SA prior and the LS-TPV prior(the lower red arrows in Fig. 3(f) and (g)). The oversmoothing in the ForWaRD restoration can not be observed in the restorations using the SA prior and the proposed LS-TPV prior (the upper arrows in Fig. 3(e), Fig. 3(f) and (g)). Compared to the restoration using the SA prior, the restoration using the LS-TPV prior also shows a better visual performance in preserving edges (the upper arrows in Fig. 3(f) and (g)).

Computation Cost Comparisons

Table 4 lists the recorded CPU times (in seconds) of different restorations for the experiments with simulated data and real data, respectively. Please note that the notations simu and real in Table 4 correspond to experiments with simulated data and real data, respectively. We can see that the restorations using SA prior and the LS-TPV prior are much computationally intensive than other methods. We define LS-TPV prior res1, LS-TPV prior res2 and LS-TPV prior res3 (f)SA prior restoration [START_REF] Chen | Effective Image Restorations Using a Novel Spatial Adaptive Prior[END_REF] (g) LS-TPV prior restoration as the original serial restorations, the restorations after optimization, and the restorations after both optimization and parallelization. Table 4 show that, with respect to the the original serial version, the proposed restorations after optimization are about 3 times faster, and the restorations after both optimization and parallelization are about 50 times faster. But we can not claim that the LS-TPV prior restoration needs less computation costs than other methods since that we did not design the parallelized versions for them.

Algorithm Monotonicity

Here, to analyze the monotonicity of the joint image/weight algorithm proposed in Section 3, we calculated the posterior function energy   ( f , w, |g) and SNR with respect to the iteration numbers for the proposed restorations. The total iteration number is set to 500. As illustrated in Fig. 4, we can see that both the calculated posterior function energy and the SNR increase monotonically during the whole iteration.

Study on Neighborhood N and n of Different Sizes

To study the roles of the sizes of N and n in restorations, we perform restorations using the proposed prior with N and n of different sizes in above simulation experiment with "Lena" image. To characterize local image structures , we do not set n larger than 7×7. Considering the fact that the total pixel numbers are different for different comparing patches n (49 for 7×7 n, 9 for 3×3 n and 1 for 1×1 n), parameter h should also change with the restorations using different patches n. The hyperparameter  and parameter h for the proposed prior are chosen by hand to produce the best stable images in terms of SNR maximization. Combinations of N and n with different sizes, and the corresponding computed SNR and the recorded CPU time cost (with CUDA parallelization) are all listed in Table 5.

We can see that higher SNR can be obtained when enlarging N and n from 1×1 to 11×11 and 7×7. We also note in Table 5 that more CPU time is needed for the restorations using the proposed prior with larger N and n. We should set the sizes of N and n based on the trade-off between performance and computation cost. We find in Table 5 that no significantly quantitative SNR differences are made between 11×11 N and larger 13×13 N. In fact,as the iteration proceeds, more global information beyond the 11×11 N will be incorporated into the regularization of the pixels.

Conclusions And Future Work Plan

In this article on image restoration, we proposed a LS-TPV prior which penalizes the total distances between neighboring patches through a constrained mixture prior mode. A conver-gent joint image/weight updating algorithm, which estimates image and weights sequentially, is proposed to overcome the heuristical OSL weight determination for the SA prior restoration in [START_REF] Chen | Effective Image Restorations Using a Novel Spatial Adaptive Prior[END_REF]. We can see that, in addition to providing effective regularization, the proposed approach can lead to stable iteration with convergence guaranteed. The application of the proposed prior model needs a pixel-wise computation of the distances between neighboring patches over a large region, which implies high computation cost of restorations. In this paper, remarkable advance in shortening computation time is achieved by optimizing patch distance computation and replacing the original pixel-wise calculation of patch distances by CUDA framework.

Further work includes applying the proposed approach in blind image deconvolutions, further justifying the proposed approach by using more image-quality measures based on human visual system, and further accelerating the computation by parallelize the PSF convolution in restorations.

  lnw b j and a normalization constraint  b∈N j w b j = 1. We thus can thus obtain the following posterior energy function:

  , = 0.5, h = 150,7×7n, 11×11N

  (a) Original image (b) Degraded image (c) Wiener restoration (d) TV prior restoration (e) Huber prior restoration (f) ForWaRD restoration [11]

Fig. 1 .

 1 Fig. 1. (a), original "Lena" image. (b), degraded "Lena" image (SNR=6.65) with 9 × 9 uniform blur and additive white Gaussian noise (variance = 0.1663). (c), Wiener deconvolution. (d), TV prior Bayesian restoration. (e), nonquadratic Huber prior Bayesian restoration. (f), the ForWaRD restoration in [11].(g), the SA prior Bayesian restoration proposed in [21].(h), the proposed LS-TPV prior Bayesian restoration.(the left arrows: blocksuppressing illustration, the right arrows: edge-preserving illustration)

Table 3 .Fig. 3 (

 33 Fig. 3(b) is the restored result from Wiener deconvolution. The restored images from TV

Fig. 2 .

 2 Fig. 2. Profile illustrations of different Restorations

  (a) Degraded image (b) Wiener restoration (c) TV prior restoration (d) Huber prior restoration (e) ForWaRD restoration [11]

Fig. 3 .

 3 Fig. 3. (a), original degraded "moon" image. (b), Wiener deconvolution. (c), TV prior Bayesian restoration. (d), nonquadratic Huber prior Bayesian restoration. (e), the ForWaRD restoration in [11].(f), the SA prior Bayesian restoration proposed in [21].(g), the proposed LS-TPV prior Bayesian restoration.(the lower arrows: block-suppressing illustration, the upper arrows: edge-preserving illustration)

Fig. 4 .

 4 Fig.[START_REF] Archer | On some Bayesian/regularization methods for image restoration[END_REF]. For the proposed restoration approach, the calculated SNR (left) and the posterior function energy (right) with respect to iteration number (1-500 iterations)

Table 1 .

 1 Parameter settings in simulated data experiment

Table 2 .

 2 Signal to noise ratio (SNR) of the observed degraded images and the restored images in the simulated experiment with "Lena" image

	Degraded Wiener TV prior Huber Prior ForWaRD SA prior LS-TPV prior
	6.65	15.06	15.54	15.59	15.50	15.84	16.02

Table 4 .

 4 CPU times (seconds) needed for different restorations for real data experiment.

Wiener filtering: simu : 58.

[START_REF] Mignotte | A segmentation-based regularization term for image deconvolution[END_REF]

, real : 64.48 TV prior res: simu : 214.96, real : 228.92 Huber prior res: simu : 270.52, real : 292.36 ForWaRD res: simu : 2128.27, real : 2411.15 SA prior res: simu : 2128.27, real : 2382.43 LS-TPV prior res1:simu : 2885.61, real : 3058.07 LS-TPV prior res2: simu : 910.48, real : 954.48 LS-TPV prior res3: simu : 58.61, real : 64.45

Table 5 .

 5 SNR and CPU time cost(seconds) for the simulation restorations using the proposed method.

	N and n	SNR	CPU times	N and n	SNR	CPU times	N and n	SNR	CPU times
	(a),N 3×3 n 1×1	15.48	17.98 sec	(b),N 3×3 n 3×3	15.65	28.59 sec	(c),N 3×3 n 7×7	15.79	53.69 sec
	(d),N 7×7 n 1×1	15.54	21.16 sec	(e),N 7×7 n 3×3	15.70	31.37 sec	( f ),N 7×7 n 7×7	15.88	50.91 sec
	(g),N 11×11 n 1×1	15.55	24.59 sec	(h),N 11×11 n 3×3	15.72	38.42 sec	(i),N 11×11 n 7×7	16.02	58.61 sec
	( j),N 13×13 n 1×1	15.55	29.54 sec	(k),N 13×13 n 3×3	15.75	77.67 sec	(l),N 13×13 n 7×7	16.07	95.31 sec

Acknowledgement

This research was supported by National Basic Research Program of China under grant (2010CB732503),National Natural Science Foundation under grant (81000636), and the Project supported by Natural Science Foundation of Jiangsu Province (BK2009012).


Olivier Baudon
email: baudon@labri.fr

Florent Foucaud
email: foucaud@labri.fr

Monika Pilśniak
email: monika.pilsniak@agh.edu.pl

A note on online and recursively arbitrarily vertex-partitionable balloons

In this note, we study arbitrarily vertex-partitionable balloons. In particular, we show that an online arbitrarily vertex-partitionable balloon has maximum degree at most 5. Moreover, we exhibit some infinite families of online and recursively vertex-partitionable balloons with four and five branches.

Preliminaries

In this note, the specific class of balloons is studied:

) is called a k-balloon with roots r 1 and r 2 and branches P 1 , . . . , P k .

We will also need the following definition. Definition 2 We denote by P B(b 1 , . . . , b i , . . . , b j , . . . , b k ) the partial k-balloon with branches of sizes b 1 , . . . , b k and roots r 1 , r 2 where the branches denoted b i are not adjacent to r 1 and the ones denoted b j are not adjacent ro r 2 .

Let n, τ 1 , . . . , τ k be positive integers such that 1≤i≤k τ i = n. Then τ = (τ 1 , . . . , τ k ) is called a decomposition of n.

We study the problem of partitioning the vertex set of graphs according to certain conditions. The following definitions have been introduced in [1], [4] and [3], respectively.

Definition 3 Let G be a connected graph on n vertices.

• G is said to be arbitrarily vertex-partitionable (AP for short) if for every decomposition τ = (τ 1 , . . . , τ k ) of n, there exists a partition V 1 . . . , V k of V (G) such that each part V i has order τ i and induces a connected subgraph of G.

• G is said to be online arbitrarily vertex-partitionable (OLAP for short) if for every decomposition τ = (τ 1 , . . . , τ k ) of n, there exists a partition V 1 . . . , V k of V (G) such that each part V i has order τ i , induces a connected subgraph of G, and V (G) \ V i induces an OLAP subgraph of G.

• G is said to be recursively arbitrarily vertex-partitionable (RAP for short) if for every decomposition τ = (τ 1 , . . . , τ k ) of n, there exists a partition V 1 . . . , V k of V (G) such that each part V i has order τ i , and both

V i and V (G) \ V i induce a RAP subgraph of G.
It is worth noting that every graph with a hamiltonian path is RAP, every RAP graph is OLAP, every OLAP graph is AP, and every AP graph contains a (quasi)-perfect matching.

OLAP balloons have maximum degree at most 5

It is shown in [START_REF] Baudon | Recursively arbitrarily vertex-decomposable graphs[END_REF] that every RAP balloon has maximum degree at most 5. We show that this holds as well for OLAP balloons.

Theorem 6 Let G be an OLAP k-balloon. Then k ≤ 5.

Proof The proof is by contradiction. Let G be a minimum counterexample to the theorem: G is an OLAP k-balloon (k ≥ 6) with least possible order.

For each λ < n, it is possible to find a subset of vertices

V λ with |V λ | = λ, G[V λ ] is connected, and G[V \ V λ ] is OLAP. Since by the minimality of G, G ′ = G[V \ V λ ]
cannot be a partial k-balloon (otherwise adding an edge to G ′ would provide an OLAP k-balloon with order strictly less than G), for each λ, one of the following conditions must be fulfilled:

1. G[V \ V λ ] is a k ′ -balloon with k ′ = k -1 (i.e. V λ is a branch of G), or 2. G[V \ V λ ] is a path (i.e. V λ consists of a root of G together with k -2 branches), or 3. G[V \ V λ ] is an OLAP caterpillar (i.e. V λ consists of a root of G together with k -2 branches but one vertex), or 4. G[V \ V λ ] is the tripode T 3 (2, 4, 6) (i.e. V λ consists of a root of G together with k -2

branches but two vertices of the same branch)

We claim that G has branches of sizes respectively 1, 2, 3, 4, 5 and 6. We show this using contradiction. Let λ ∈ {1, 2, 3, 4, 5, 6} and suppose there is no branch of size λ in G.

• λ ∈ {1, 2, 3}: for any choice of V λ , G[V \V λ ] is either a partial k-balloon or a tree with a vertex of degree 4 (since at most two branches can be completely included in V λ ), a contradiction.

• λ = 4: similarly as in the previous case, for any choice of V λ , G[V \ V λ ] is either a partial k-balloon or a tree with a vertex of degree 3. By Observation 4, G[V \ V λ ] cannot be a caterpillar, thus it is T 3 (2, 4, 6) and G = B 6 (1, 1, 1, 2, 4, 6), a contradiction since G has a branch of size 4.

• λ = 5: by the previous cases we know that G has branches of sizes 1, 2, 3 and 4. G[V \ V λ ] cannot be a path since by Observation 4, five vertices cannot cover four branches and a root.

If G[V \ V λ ] is an OLAP caterpillar, then necessarily k = 6, the caterpillar is Cat(4, 5) and G contains three branches of size 1:

G = B 6 (1, 1, 1, 2, 3, 4). But then following Observation 4, G is not OLAP. Finally, if G[V \ V λ ] is T 3 (2, 4, 6
), V λ must either totally cover the branch of size 3, or the branch of size 2 plus a vertex of the one of size 3. But then it is only possible to totally cover two branches, and the remaining part has a vertex of degree at least 4 and is not T 3 (2, 4, 6).

• λ = 6: by the previous cases we know that G has branches of sizes 1, 2, 3, 4 and 5. Now suppose there is one additional branch b i (i ≤ 6) such that b i ≤ 5. By Lemma 5 there can be at most one such additional branch, and then b i = 1. Thus G = B k (1, 1, 2, 3, 4, 5, ...). But now the only potential choice for V λ would be the three smallest branches, a root and one additional vertex.

If k > 6, G[V \ V λ ]
is a tree with a vertex of degree at least 4, which is not OLAP. Otherwise it is T 3 (2, 4, 5), which is not OLAP either. Thus, G = B k (1, 2, 3, 4, 5, b 6 , ...) with b 6 ≥ 7. But again for λ = 6 it is not possible to decompose it. Thus G is not OLAP, a contradiction.

So, G = B k (1, 2, 3, 4, 5, 6, ...). But using Lemma 5, G is not AP, a contradiction.

Families of OLAP balloons

In this section, we show that there are infinitely many OLAP k-balloons, for each k ≤ 5. The case where k ≤ 3 is trivial since then any k-balloon has a hamiltonian path and is therefore RAP.

OLAP 4-balloons

Proposition 7 B(1, 1, n, m) is OLAP iff at most one of n and m is odd.

Proof Easily partitionable for λ = 1, 3. For λ = 2, one can prove by induction that B(1, 1, 2k, m) and B(1, 1, n, 2k) is OLAP.

Proposition 8 Let a ≤ b ≤ c be such that Cat(b + 1, a + c + 2) or Cat(c + 1, a + b + 2) or Cat(a + 1, b + c + 2) is OLAP. Then B(1, a, b, c) is OLAP.
Proof These caterpillars are spanning subgraphs of B(1, a, b, c).

Proposition 9 B(1, 2, 2, k), k ≥ 3, is OLAP iff k ≡ 0, 1 mod 3.
Proof Easily partitionable for λ = 3. If k ≡ 0, 1 mod 3, also OK for λ = 3 because of Cat(3, k + 1) which is OLAP. Now suppose k ≡ 2 mod 3. For λ = 3, for the balloon to be partitionable we should have either:

• P B(1, 2, 2, k -3) OLAP. • P B(1, 2, 2, a, b) (a + b = k -3) OLAP
But in both cases these partial balloons are not (3, . . . , 3)-partitionable.

Proposition 10 B(1, 2, 3, k), k ≥ 3, is OLAP.
Proof Observe that it is true for k ≤ 3 since Cat(3, 5) and Cat(4, 5) are spanning subgraphs of P B(1, 1, 2, 1), P B(1, 1, 2, 2) and P B(1, 1, 2, 3), respectively. Now, suppose the proposition is true for all i ≤ k -1 and consider P B(1, 1, 2, k). For λ ≤ 4, P B(1, 1, 2, k) can be partitioned into P λ and P n-λ . For λ = 5, it can be partitioned into Cat(2, 3) and P k+1 . For λ = 6, it can be partitioned into B(1, 1, 2) and P k . Finally, for λ ≥ 7, it can be partitioned into P B(1, 1, 2, λ -6) and P k-λ+6 , the first being RAP by induction hypothesis. Proof For any λ ≤ n 2 , P B(1, 2, 3, k) must be partitionable into two vertex-disjoint RAP subgraphs of orders λ and nλ.

• λ = 1: P 1 + P n-1

• λ = 2: P 2 + P n-2

• λ = 3: P 3 + P n-3

• λ = 4: P 4 + P k+4

• λ = 5: P 5 + P k+3

• λ = 6: P 6 + P k+2

• λ = 7: Cat(3, 4) + P k+1

• λ = 8: B(1, 2, 3) + P k Note that by the previous cases, the proposition is true for n ≤ 16 (i.e. k ≤ 8). Now, let λ ≥ 9, suppose it is true for any i ≤ k -1, and consider P B(1, 1, 2, 3, k). Now one can partition the graph into P n-λ and P B(1, 2, 3, λ -8), which is RAP by induction hypothesis. Note that by the previous cases, the proposition is true for n ≤ 14 (i.e. k ≤ 2). Now, let λ ≥ 8, suppose it is true for any i ≤ k -1, and consider P B(1, 1, 2, 3, 2k).

RAP 5-balloons

• λ even:

if λ ≤ 2k, one can partition the graph into P λ and P B(1, 1, 2, 3, 2kλ) by placing P λ at the end of the last branch, and use the induction hypothesis.

if λ > 2k, then nλ ≤ 7 (since both 2k and λ are even), thus one can use the cases where λ ≤ 7.

• λ odd. Then partition the graph into P B(1, 1, 2, 3, λ -9) (note that it is RAP by induction hypothesis since λ -9 is even and positive) and P 2k-λ+9 .

This proves the existence of an infinite family of RAP 5-balloons:

Corollary 19 B(1, 1, 2, 3, 2k) is RAP for any k ≥ 1.

In [START_REF] Baudon | Recursively arbitrarily vertex-decomposable graphs[END_REF], it is shown that B(1, 1, 1, 2, 4) and B(1, 1, 2, 2, 3) are RAP.

Proposition 20 Other RAP 5-balloons are B(1, 2, 2, 3, 4) and B(2, 2, 3, 4, 6).

Observation 4 (

 4 [3]) Let B = B(b 1 , ..., b k ) be an AP balloon with b 1 ≤ ... ≤ b k . If n is odd, B contains at most three branches of odd size. If n is even, B contains at most two branches of odd size. Lemma 5 ([2]) Let B(b 1 , ..., b k ) be an AP balloon with b 1 ≤ ... ≤ b k . Then for each i ≤ k we have 2b i ≥ j<i b j .

  Proposition 18 P B(1, 1, 2, 3, 2k) is RAP for any k ≥ 0.Proof For any λ ≤ n 2 , P B(1, 1, 2, 3, 2k) must be partitionable into two vertex-disjoint RAP subgraphs of orders λ and nλ.• λ = 1: P 1 + B(1, 2,3, 2k) (RAP by Proposition 17) • λ = 2: P 2 + B(1, 1, 3, 2k) (RAP by Proposition 14) • λ = 3: P 3 + B(1, 1, 2, 2k) (RAP by Proposition 15) • λ = 4: P 4 + Cat(4, 2k + 1) • λ = 5: Cat(2, 3) + P 2k+4 • λ = 6: P 6 + Cat(2, 2k + 1) • λ = 7: Cat(3, 4) + P 2k+2

  Proposition 16 P B(1, 2, 2, k) and B(1, 2, 2, k) are RAP iff k ≡ 1, 2 mod 3.Proof On the one hand, observe that Cat(3, k + 4) is a spanning subgraph of P B(1, 2, 2, k) and B(1, 2, 2, k), which is RAP when k ≡ 1, 2 mod 3. On the other hand, if k ≡ 0 mod 3, they are not (3, . . . , 3)-partitionable.

	Proposition 17 P B(1, 2, 3, k) is RAP for any k ≥ 3.

OLAP 5-balloons

The proof of the following result is omitted.

Proposition 13 Let G be a graph of Table 1. Then G is OLAP.

Families of RAP balloons

In this section, we show that there are infinitely many RAP k-balloons, for each k ≤ 5.

RAP 4-balloons