Quan Li

Li Chen

Junhai Yong
email: yongjh@tsinghua.edu.cn

Huamin Qu
email: huamin@cse.ust.hk

Robert Wang

Detan Si

Visual Analysis of Retweeting Propagation Network in a Microblogging Platform

Keywords: H.4 [Information Systems Applications]: General, H.5.2 [Information Interfaces and Presentation]: User Interfaces-Graphical user interfaces (GUI), E.1[Data]: Graphs and Networks retweeting propagation, dynamic graph visualization, interaction, visual mining, micro-blogging

As a novel type of real-time social networking service, microblogging has already become ubiquitous and an irreplaceable tool. Tracking in the pulse of retweeting propagation is important and meaningful. In this paper, we investigate how information propagation in a specific microblogging platform evolves to identify relevant patterns and understand dynamic attributes of information propagation and the underlying sociological motivations. More specifically, based on the node-link diagram, we propose three efficient strategies to map the multiple attributes of information propagation graph to appropriate visual elements. For revealing the dynamic attributes, we propose two models: the depthvarying and the time-varying parallel data model to illustrate the temporal evolution efficiently. We also present a novel method by combining the traditional scatter plot with Hough transformation to represent the distribution of propagation instances and trace the propagation speeds. We integrate our methods to a visual mining tool and develop several interactive features. We demonstrate how our approaches improve the understanding of the propagation graph from a visual perspective by employing propagation datasets collected from Sina Weibo, the largest microblogging service provider in mainland China. Meanwhile, this visual mining tool has been evaluated by data analysts and successfully used in Sina Corporation as a helpful assistant to them.

INTRODUCTION

As a type of social networking service, microblogging, such as Twitter, Sina Weibo has already become ubiquitous in Web 2.0 era. Tracking in the pulse of retweeting propagation is important. For corporations, it enables them to get feedback of user coverage and gain insight into how to improve and market better. For individuals, the abundance of information and opinions from diverse sources in microblogging platforms helps them tap into the wisdom of crowds, understand the retweeting behavior and aid in making more informed decisions [START_REF] Leskovec | Social media analytics: tracking, modeling and predicting the flow of information through networks[END_REF]. Information propagation in social networks takes on different structures and scopes. However, little work has been done to reveal the coupled dynamics of both the structure and the information propagation. Most researchers have focused on the qualitative study or mathematical modeling of information propagations in social networks [START_REF] Grabowicz | Social features of online networks: the strength of weak ties in online social media[END_REF]. By contrast, in this paper, our objective is to track information propagation and mine the patterns from a visual perspective.

How to represent the graph of social network has been heavily researched [START_REF] Görg | Dynamic graph drawing of sequences of orthogonal and hierarchical graphs[END_REF]. The major challenge is the computing cost of graph drawing criteria. As the information propagation network involves multiple attributes and dynamic evolution, how to model the multi-attributes and map them efficiently to visual elements is a big challenge. A prevailing approach to deal with the dynamic attribute is to use animation [START_REF] D I E H La N Dc | Graphs, they are changing[END_REF] [START_REF] Frishman | Online dynamic graph drawing[END_REF]. However, animation mapping leads to high cognitive efforts due to human's limited short term memory. It is difficult to observe trends over longer periods or compare two non-subsequent time steps. Abrupt changes of the graph might furthermore destroy the viewer's mental image of the graph [START_REF] Burch | Parallel edge splatting for scalable dynamic graph visualization[END_REF].

In this paper, we investigate how information propagation in a specific microblogging platform evolves and carry out visual mining tasks on the social network graph. Our visualization tool provides several interactive features to explore dynamic graph data. We demonstrate how the above techniques and methods increase the understanding of the information propagation from a visual perspective in several application scenarios. Meanwhile, this system has been evalu-atedandusedintheW eibodepartmentofSinaCorporation. The major contributions of this paper are: (1): A comprehensive visual analysis system for mining retweeting propagation patterns in a microblogging platform; (2): Three novel visualization schemes, hierarchy-aware layouts, an accumulated filtering scheme and micro-aggregation, for multiattribute propagation graphs; (3): Novel depth-varying and time-varying data models and visualization for scalable graph visualizations; (4): Case studies with real data to demonstrate the effectiveness of the system.

The remainder of this paper is organized as follows. We introduce the related backgrounds and existing work and then we describe our data and the multi-attribute graph. After that, we propose three novel visualization schemes to represent multiple social attributes regarding to specific business requirements. Then, we deal with the dynamic attributes. Following are several interactive features. Finally, two case studies are introduced to demonstrate the effectiveness of our system.

RELATED WORK

The related work falls into two categories: information propagation on social network and how to visualize the dynamic propagation patterns.

Information propagation on social network

Information propagation in social networks has been heavily studied in the past few years, especially with the rapid emergence of microblogging services, and has drawn considerable research attention in the area of data mining. In 2009, Sun et al. [START_REF] Sun | Gesundheit! modeling contagion through facebook news feed[END_REF] investigated the propagation of Facebook's News Feed and created a social network based on users and fans of Facebook for analysis. Kwak et al. [START_REF] Kwak | What is twitter, a social network or a news media[END_REF] used the relationships existed in the followers and fans in Twitter to construct a social network for analysis. Sakaki et al. [START_REF] Sakaki | Earthquake shakes twitter users: real-time event detection by social sensors[END_REF] proposed a real-time event detection system using Twitter data. Diakopoulous and Shamma demonstrated the use of timeline analytics to explore the 2008 Presidential debates through Twitter sentiment [START_REF] Diakopoulos | Characterizing debate performance via aggregated twitter sentiment[END_REF]. Shamma et al. later carried out a similar analysis at the President Obama's Inauguration [START_REF] Shamma | Conversational shadows: Describing live media events using short messages[END_REF]. Gaffney performed a retrospective analysis of the recent Iran election, in which Twitter played a pivotal role in news reporting [START_REF] Gaffney | # iranelection: quantifying online activism[END_REF]. Jansen et al. found that tweet sentiment analysis could provide useful insights into product opinion, an idea leveraged for the sentiment visualization in TwitInfo [START_REF] Jansen | T w i t t e rp o w e r : T w e e t sa se l e c t r o n i cw o r do fm o u t h[END_REF]. Ho et al. focused on ways to measure how information is propagated and considered how to quantify a person's capability to disseminate ideas via a micro-blog and measure the extent of propagation of a concept in a micro-blog [START_REF] Ho | Modeling and visualizing information propagation in a micro-blogging platform[END_REF]. The above works use microblogging data for mathematical modeling or statistical analysis, without considering the dynamic and multiple attributes. Meanwhile, to the best of our knowledge, works on visualization and analysis for retweeting propagation are scarce, many of which use statistical representations, timeline, or social graphs to visualize the propagation graph.

Dynamic Propagation Visualization

The traditional node-link diagrams are effective to reveal propagation paths from a visual perspective. However, this graph layout cannot deal with massive graph data efficiently. The Force-Directed visualization layout model has been a widely adopted model, such as the classic Force-Directed algorithm [START_REF] Eades | A heuristic for graph drawing[END_REF], the KK algorithm [START_REF] Kamada | An algorithm for drawing general undirected graphs[END_REF], and the FR algorithm [START_REF] Fruchterman | Graph drawing by force-directed placement[END_REF]. Researchers have done a lot of simulation to reduce the algorithmic complexity [START_REF] Barnes | A hierarchical 0 (n log iv) force-calculation algorithm[END_REF]. In our improved algorithm, we employ the traditional node-link diagram and the spring model. However, we optimize the initial position of each vertex before iterations.

As for the time dimension in dynamic graphs, traditional approaches use animation [START_REF] Frishman | Online dynamic graph drawing[END_REF] [START_REF] Von Landesberger | Visual analysis of large graphs[END_REF]. However, they require a heavy user's cognitive load. Burch, et al. proposed a novel method, mapping time to space and showing the subsequent graphs in a static diagram [START_REF] Burch | Parallel edge splatting for scalable dynamic graph visualization[END_REF]. Our approach is inspired by their work. However, the graph model and mapping strategy in our work are quite different from theirs.

MULTI-ATTRIBUTE GRAPH MODEL

Difference Between Weibo and Twitter

Even though there existed several difference between the weibo and the western Twitter, we focus on the difference of replies, comments, and retweeting. For Twitter, the replies and comments appear independently in the feed, while on a weibo, they are listed under the entry, like a traditional blog. Replies are response to what someone else has written on Twitter whereas ReTweets are when you want to share what someone else has written with others so you copy and paste it as a Tweet update. The action of "retweeting", which is central in this paper, is a little bit different from that of the famous Twitter network. If you want to retweet one's tweet, you can directly retweet it or add your comments in addition to the original tweet and retweet them, provided that it does not exceed the maximum length of a tweet. Weibo supports the function of Reply and ReTweet simultaneously and you can add some comments in addition to the tweet. In this paper, we record the retweeting action, including retweeting without additional comments and with other comments.

Multi-Attribute Graphs

The propagation record consists of eight attributes. For all the records, the rmid (root message id) and ruid (root user id) are the same since all the retweeting behaviors aiming at a specific message posted by a specific user. In the propagation network, the message id and user id are unique. The pmid (parent message id) refers to the message id of the source vertex of the current retweeting occurrence and the puid (parent user id) represents the corresponding user id. The cmid (current message id) is the current message id of the retweeting occurrence and the cuid (current user id) is the corresponding user id. To map the multi-attribute information propagation graph into meaningful visual elements, we design a multi-attribute graph model. In this paper, we represent an information retweeting instance as a single node and an occurrence of retweeting behavior as a single edge. A user may retweet more than once, given three, for a single micro-blog, resulting in three vertices, representing three retweeting instances. Therefore, we consider our model as a Tree T with one vertex r as root node: T =( V, E, r). The root node is the very beginning of propagation. Re-lationships between vertices are expressed by T :V e r t i c e s sharing a common parent in T belong to the same "group". The time attribute represents the occurrence time of this retweeting behavior and the optional gid (group id) means the group id associated with this record. The retweeting propagation graph may become very large and is dynamic in nature, which pose special challenges for visualizations. In t h en e x tt w os e c t i o n s ,w ew i l li n t r o d u c es o m en e wm e t h o d s we develop to address these two issues.

MULTI-ATTRIBUTE GRAPH

We present our methods to visualize large multi-attribute propagation graphs. First, we introduce a hierarchy-aware layout which can better reveal the retweeting relationship between nodes. Then we present a simplified layout to reduce the visual clutter for large data. Finally, a micro-aggregation scheme is proposed to visualize the groups in information propagation networks.

Hierarchy Aware Force-Directed Layout

In this section, we propose a Hierarchy Aware Force-Directed layout. We adopt the spring model and the Nbody algorithm and use a Quad-Tree data structure [START_REF] Barnes | A hierarchical 0 (n log iv) force-calculation algorithm[END_REF] to realize the Force-Directed layout. Each particle is stored in the tree structure, and corresponds to a leaf node in the tree. Each non-leaf node records the center coordinates of its children nodes and the sum of quality of all child nodes. When calculating the force of particles, we first traverse from the root node of the hierarchical tree. If the distance between the node and the particle exceeds a preset threshold value, we use the quality and the location of the node to calculate the force effect for particles. Otherwise, we continually traverse along the tree. More importantly, we change the initial position of each node according to its parent node's position in the hierarchy before each iteration, thus greatly decreasing the time complexity. The detailed framework is shown in the following steps.

(1) Instantiate the layout algorithm and add two sublayout algorithms: Nbody algorithm and spring model. We use a dataset to compare the visualization results by adopting different algorithms. Firstly, we, respectively use the original Force-Directed layout algorithm and our improved method to show the graph with group clustering representation, the first time-step plot and the last time-step plot. We can clearly see our Hierarchy Aware Force-Directed layout is much clearer and has less visual clutter, especially when it reaches to the final time-step. Meanwhile, the group structures in the network are better preserved (Fig. 1).

Simplified Layout

Visualizing graphs in node-link diagrams may lead to visual clutter in the final display for large graphs even with an Visualization with the Hierarchy Aware Force-Directed layout algorithm (Bottom, From left to right: graph with group colors, first time-step, final time-step). Colors represents different retweeting groups appropriate layout algorithm. In a microblogging retweeting propagation network, the most concerned part is the key nodes for retweeting. In our model, an edge represents a single retweeting behavior; therefore, nodes with higher outdegrees should be considered more significant according to the sociology and graph theory. However, it is not adequate to just consider the direct retweeting number (out-degree) that one single vertex brings in. The accumulated retweeting number that propagation paths starting from this vertex should be taken into account.

Data Evaluation

We evaluate the simplified layout method based on one microblogging retweeting propagation dataset. It, which was firstly published by China Newsweek, covers 594 retweeting propagation instances. The longest path of propagation is 3 and the average path is 2.75. It starts at 2011-11-14 10:33:00 and ends at 2011-11-14 23:56:24, which lasts about 803 minutes. The original propagation graph and the simplified graph are shown in Fig. 2. When we visualize the top one team for each depth, we get the corresponding information of the most related users, as shown in Fig. 2. The number of fans and followers of "China Newsweek" are 1572, 079 and 1, 793, respectively. "Direct" means the number of nodes that directly retweet the message published by their superior parent in hierarchy. "Accumulation" refers to the number of all the nodes that directly or indirectly retweet the message. Of all the 8 nodes, the node of China Newsweek appears 4 times. One is functioned as the root node of this microblogging retweeting propagation network. The other three lay in the first depth of this network and act as retweeting nodes. We establish that the retweeting behavior of the root node itself plays a vital role in retweeting information.

Propagation Analysis

It is significant to trace the information flow and diffusion along the propagation path, especially those retweeting behavior and retweeting messages published by them. These key retweeting behaviors and the additional messages attached on the original message are essential for tracking the flow of sentiment through networks and emergence of polarization. Take this propagation dataset for example, at 11:21:35, China Newsweek retweeted its message published at 10:33:00. At 13:49:01, it retweeted again and the same went to the time of 15:20:19. We set the threshold accumulation value higher and visualizing the top three at each depth, and we get these messages: We conclude that at this circumstance, each key retweeting behavior is another extension to the original topic/message. In other instances, the emergence of polarization may exists, since different polarizations arguing on a certain topic are likely to trigger large portion of web users' discussion and bring in new retweeting propagations. In this propagation network, "China Newsweek" is a public news media verified by Sina Weibo, therefore, it more likely retweets a summary of the related or similar topics, or the follow-up reports, which are likely to be proposed by single individuals.

Micro-aggregation

An important phenomenon for microblogging users is that the emergence of groups. The groups could be all the users having the same hobbies, or the same labels and all the corresponding topics gathered inside. When visualizing and mining these propagation records with "group" information, it is not intuitive to just plot them in the traditional node-link layout without giving more visual elements to help analysts better understand how the information flows within and between different groups. In this section, we propose a novel visual representation to visualize and analyze the microblogging propagation with aggregation information and evaluate with a real propagation dataset. This dataset covers 2000 instances. The number of micro-group is 38. We first collect all unique group ids and assign different color to each group and represent different groups as aggregation circles, then iterate all vertices of Set(V )i nGraphV ertex, assign the vertex to the corresponding aggregation circle based on group id. The color of each vertex is based on its current user id (cui). We initially partition the screen according to the size of aggregation circles. For each aggregation group, we apply the Force-Directed layout algorithm to the vertices inside the group. The movement of the vertices is limited to the boundary of the aggregation group. For the connections between aggregation groups, we just link the vertices and cancel the force between these vertices. By observing Fig. 3, we find almost all the intra-connections happen from one specific aggregation group to another. By observing the color of vertices, we find that although there are 2000 vertices, they only correspond to 7 user ids, which means all the propagations are due to these 7 users. Among the 7 users, we find that one particular user is excessively active and occupies most of the propagations. By analyzing the information of this user, we find that this user has joined many different groups and has retweeted this microblogging message to all the groups. By observing the largest group, represented by yellow color, we find that this user had continually retweeted this message within this group, and no other users had published any messages, which means this user is continually refreshing the screen by retweeting the message inside the aggregation group. From this example, we can clearly see that micro-aggregation can significantly reduce visual clutter and reveal the groups in the propagation network.

SCALABLE DYNAMIC GRAPH VISUAL-IZATION

Another important and special feature of the propagation network is the temporal evolution. In this section, we propose a visualization method to use a single-image to explore dynamic and hierarchical propagation graphs. Two data models, the depth-varying data model and the time-varying data model have been proposed. In the depth-varying data model, we map each depth of the graph to a single axis on screen. The vertices belonged to each depth are arranged along the vertical axis according to their retweeting time. A line that connects the vertices starts from the left most side and ends at the right most side, as long as this propagation path has reached that depth. One link represents one propagation record. In the time-varying data model, we m a pe a c ht i m es t e po ft h eg r a p ht oas i n g l ea x i so ns c r e e n . The vertices belonged to each time step lies at certain posi-tion according to their hierarchy in the graph. To give the distribution of clusters in graphs, we apply histogram view along the vertical axes.

Depth-varying Parallel Data Visualization

Based on the idea of parallel coordinates [START_REF] Inselberg | Parallel coordinates: a tool for visualizing multi-dimensional geometry[END_REF] and the depth-varying attribute in the propagation graph, we plot the depth-varying parallel data visualization. The result is a ss h o w ni nF i g .4 .

Time-varying Parallel Data Visualization

Different from the depth-varying parallel data model, we extract all the time steps and plot them as vertical axes. Since in most cases, the number of time steps is much higher, and thus more axes would be generated. In this model, parallel axes represent time steps ranging from the beginning time to the ending time. However, not all the propagation records continually cover all the time steps. Most of the propagation paths just jump from one time step to another. Therefore, we need to deal with this "jumping" and the saltation elegantly. For those axes on which there are no propagation occurrences, we simply ignore these axes and connect those neighboring propagation instances together. The visualization plot of the same dataset is shown in Fig. 5. In Fig. 5, it is easy to detect how the depths of propagation records distribute on a single axis from bottom to up, and from left to right, and thus we can get an overview of how all the propagation records change their paths as the time goes by. Meanwhile, we can add interactions to this plot and facilitate analysts to conveniently understand the retweeting propagation dataset better.

Hough-scatter Diagram

In this section we consider the speed of propagation and the patterns in the retweeting propagation of the depthvarying parallel data visualization plot. The traditional scatter plot only serves to represent dataset according to the distribution of values on dimensions, rather than present the underlying patterns in the dataset. By Hough transformation [START_REF] Tatu | Combining automated analysis and visualization techniques for effective exploration of high-dimensional data[END_REF], the trend of lines between two neighboring axes can be easily detected. In our case, the transformation from the depth-varying parallel data model to the Hough plot makes the propagation speed apparent. In addition, in order to detect the "propagation clusters" in the Hough transformation plot, we add a density-based pixel plot, indicating the propagation summits from one depth to the successive one. Fig. 6 demonstrates the density-based pixel plot.

In Fig. 6, the upper part is the traditional scatter plot with density and the lower part is the Hough transformation plot with density. The green line indicates the "missing" values and we fill up them with -1. By observing the neighboring space between two adjacent axes, the "diagonal phenomenon" is common for scatter plot and the "horizon phenomenon" happens frequently. The reason for explaining this phenomenon is that a large number of the retweeting propagations happen instantly. In the upper panels, vertical clusters of red points can be witnessed as they are propagation summits.

INTERACTIVE FEATURES

We combine several necessary interactive features to help users carry out those tasks and understand the properties of retweeting propagation networks.

Depth Filtering: Due to the large scale of certain retweeting propagation networks, we add a depth filtering strategy to control the data loading process and make the system perform smoothly and easily (Fig. 7). Retweeting Filtering/Time Controlling: We provide users with freedom to adjust the retweeting number slider to filter out those nodes whose retweeting number is less than the threshold (see Fig. 8). As the network evolves with time, it is also quite intuitive to add a time slider to control the evolution of retweeting propagation graph, as shown in Fig. 9.

Lens: The lens lifts the focused items up while keeping other items staying still in the view plane. Meanwhile,

CASE STUDY 7.1 Viral Online Marketing

In the first case, we collect the retweeting propagation records of a message published by a famous hotpot restaurant. This retweeting propagation dataset covers a large number of propagation records and lasts for a fairly long time, about ten days. But the number of its followers in the follower-fan network is limited. The success of micromarketing in this case has attracted the attention of the relevant departments. We first plot the statistics of this Next, we employ the Hierarchy Aware Force-Directed layout to plot all the records. We slightly adjust the time slider to watch the animation from the beginning to the current time and get the following four layouts after different lasting time (Fig. 13). In the final visualization, we then increase the value of depth slider to add edge restraint. The fifth and sixth plots of Fig. 13 show the initial step and the final result by adding edge restraint. Depth information of each vertex i se n c od e dw i t ht h ec o l o ro fe a c hv e r t e x . T h e r e f o r e ,v e r t i c e s with the same color belong to the same level. We turn on the display of group information and level displaying to see the boundary outline of the retweeting propagation structure, as shown in the seventh and eighth plots of Fig. 13. Most propagation paths are slim since most propagation instances just directly retweeted the original message from the root user.

By "Link + View" technique, we transform the fetched information from the node-link diagram to the depth-varying parallel data model and emphasize the highly focused items with blue lines (the second plot of Fig. 14). We extract these key instances and find out that there are few intraconnections between them, telling us that those key nodes may not be strongly linked and may be connected by those unimportant nodes.

In the depth-varying parallel data model, we detect some interesting phenomenon: in propagation section 1 ranging from time 0 to time 3230 minutes (Fig. 15), the longest path has taken place. The longer the lasting time lasts, the shorter the propagation path is. One propagation summit happens in about 200 minutes, and it brings in three major propagation paths from depth 1 to depth 2, even to depth 8. The reason to explain for the long lasting time of propagation is that users on depth 1 continually retweeted the message after several breaks.

To take a deeper investigation about the relevant web users and their retweets, we select several interesting patterns, as illustrated in Fig. 16. All the selected patterns are propagated to distant depth of the retweeting propagation network and cover a wide range of web users. We are curious about why the whole propagation covers so many users while the number of its followers in the follower-fan network is limited. After crawling the relevant content of their retweets, we find that the effects of Circle of Friends have dramatic consequences. The web users who had retweeted this micro-blog were continually recommending this message to and interacting with their friends. By observing their mutual interactions, we conclude that it's quite important and necessary to promote discussion and interaction among web users. In Fig. 17, the longest path is found in 647 minutes after the beginning of publishing the message. From depth 0 to depth 1, two distinguished parts are formed after a long time break. Therefore, we are interested in the following questions: (1) How to explain the two distinguished parts when propagation from depth 0 to depth 1? (2) Why does the longest path appear after a long time? In most cases, the first half part of the longest path always appears at the bottom of the visualization plot. (3) Why does not the propagation from depth 0 to depth 1 at an early time last to a deeper level?

To filter out those unimportant nodes and keep the main structure of the graph, we highlight the center of each cluster in a node-link diagram, map them to the depth-varying model and then employ our simplified layout algorithm to keep the top two accumulation degree of each depth. We fetch the related information of users and tweets (Fig. 18). Then, we find out the reasons: The original message is about article recommendation published by root user A. This arti-clewaswrittenbyuserB.However,theURLlinkintheoriginal message was incorrect until in about 647 minutes when the article's author-user B retweeted this message and corrected it with the right URL link, which led to a great rise of propagation from depth 0 to depth 1 again, since this corrected message appeared on the retweeting webpage of root user A. At the early time from depth 0 to depth 1 the propagation did not reach to a deeper level, although there was a key figure C having retweeted this message and brought in some propagations. However, we cannot find that that key figure's retweeting record by API connection. Therefore, we conclude that user C must have deleted this original message that contains the wrong URL link, resulting in no more propagations. And the propagations that passed the article author-user B increased significantly and formed the longest path in the overall retweeting propagation network.

Then, we are interested in how the propagation paths split, emerge and vanish with the spread of the scope of the constantly propagation and deepening. In the scatter plot and Hough plot, we can clearly see the changes of propagation paths. We select the interested pattern with a blue circle that presents a long propagation path, and this will simultaneously reflect in other relevant panels and depth-varying data view. The selected pattern happened from propagation depth 1 to depth 2 splits into three parts. Two of them are still spreading and the other part has vanished, as show in Fig. 19 (top). As we have demonstrated previously, the "diagonal phenomenon" in the scatter plot and "horizon phenomenon" in the Hough plot are apparent, which indicate the instantly retweeting behaviors from parent nodes to child nodes. We select the "horizon line" in the Hough plot and check out the corresponding patterns. In this case, two thirds of the propagation instances from depth 2 to depth 3 are due to the instantly retweeting behaviors from depth 1 After a long time break, the blue line in "Level 1-Level 2" panel becomes the major propagation in the next panel and then breaks into three major branches (Fig. 20 (left)). In the time-varying parallel data visualization plot, the distribution of depth and time scale is apparent. The structure of this graph presents a major trend: In about 647 minutes after the beginning of propagation, retweeting from the root user still occupies a lot, much denser than propagations from depth 1 to depth 2, or depth 2 to higher depths (Fig. 20 (right)).

Interview with Data Analysts and Technical Support Engineer

We conducted extended one-on-one interviews with twenty data analysts and one technical support engineer with very strong domain expertise. Most of the data analysts were majoring in statistics. The technical support engineer is responsible for providing technical assistance between the service platform and our visualization system.

All the data analysts were intrigued by the interactive features that our system provides for examining propagation patterns and highlighting those key propagation instances. Before adopting our system, they used a simple tool which generates a static image by directly plotting all the propagation records. They complained that the previous tool can only deal with relatively small data and does not support any user interactions. For our system, the data analysts were attracted to the interactivity and felt that the multilayouts provided significant value. They liked the mechanism of the system that the displaying of propagation nodes is triggered by interactions or time-steps, instead of displaying all of them. They thought that the improved Hierarchy Aware Force-Directed layout is more in line with their intuitions.

Referring to the detailed designs of our system, one data analyst said "it gives me an intuitive way to examine propagation instances just by taking a glance at the color encoding distribution". And the animation provided in the system makes the data analysts convenient to control their interested time-steps. Another data analyst responsible for What they liked most, is the two proposed models: the depth-varying and the time-varying parallel data model which intuitively present the dynamic attributes. "I can get a lot of information", said one of the data analysts, "these two models represent the distribution of time, depth, propagation intervals and even propagation clusters clearly". When asked whether this system helped them understand the underlying information of propagation data and improve their work efficiency, they responded "of course". We also design

CONCLUSION AND FUTURE WORK

In this paper we have investigated how information propagation in a specific microblogging platform evolves from a visual perspective. Scenario-oriented layouts have been proposed to map the multi-attribute propagation network into appropriate visual elements. We have also introduced two parallel data models to depict the dynamic retweeting propagation graph in a single image, which has clear advantages over the traditional presentation for dynamic retweeting propagation graph. The proposed methods and several interactive techniques have been integrated into our visual mining system, which has been evaluated and used in the Weibo department of Sina Corporation. In the future, we will exploit some NLP techniques, such as opinion mining and clustering to analyze the content of messages and provide a deeper understanding of propagation.

( 2 )

 2 Calculate filtering condition depth: current time and depth value. (3) Filter all vertices to invisible and set those satisfying the filtering conditions visible and put them into a set. (4) Initial the vertex coordinates and positions in Set(V). (5) Assign v's father node's position to v's current position. (6) Iterate on the set and get all the visible edges. (7) Calculate the force between the two vertex of each edge and initial the length of a spring to represent the property of the edge. (8) Update nodes' position.

Figure 1 :

 1 Figure 1: Visualization with the original Force-Directed layout algorithm (Top, From left to right: graph with group colors, first time-step, final timestep). Visualization with the Hierarchy Aware Force-Directed layout algorithm (Bottom, From left to right: graph with group colors, first time-step, final time-step). Colors represents different retweeting groups

Figure 2 :

 2 Figure 2: The original propagation graph and the simplified graphs (From top left to bottom right: the original propagation graph, the simplified graph with top retweeting number for each depth, the simplified graph with user's information, and the simplified graph with top three retweeting number for each depth)

  [10:33:00: kinship fall: who is "robbing" your donations?], [11:21:35 Dalian City: Parents hijack 80,000 donation; a two-year-old scalded child struggling in hospital], [13:49:01 Jiangsu Province: Teacher privately takes possession of donation; Students suffered and helpless], [14:20:58 Charity Collapse], [16:25:36 Every donation should be a hot potato on the hand], and [17:20:16 Charitable organization should not be bypassed].

Figure 3 :

 3 Figure 3: Visualization of micro-aggregation algorithm for the dataset with micro-aggregation information (From left to right: initial positions for each aggregation, inter/intra-connection between aggregations)

Figure 4 :

 4 Figure 4: Depth-varying parallel data visualization (top) and interactions (bottom): dragging the depth axis up and down and highlight the propagation path By a deep visual analysis with Fig. 4 (top), we can easily detect the longest propagation path, as the green line indicating. From the tooltip we can get more information along the longest propagation path, such as when the microblogging message arrives and the current depth. By showing the histograms of each axis, the distribution of vertices is quite intuitive. Meanwhile, we facilitate analysts by combining several interactive techniques and "Link + View" techniques. Fig. 4 (bottom) shows the related supporting interactive techniques: dragging the depth axis up and down to highlight the propagation path.

Figure 5 :

 5 Figure 5: Time-varying parallel data visualization (After dealing with "jumps") and interactive exploration

Figure 6 :

 6 Figure 6: Density-based pixel plot for scatter plot and Hough plot. The green line indicates the "missing" values and the red points are the places that propagations occur

Figure 7 :

 7 Figure 7: Depth controlling incremental visualization (depthV alue =1 ,2 ,5 ,1 2 ) ( O n l yd i s p l a ye d g e so f nodes of which the depth is 1, 2, 5, and 12)

Figure 8 :Figure 9 :

 89 Figure 8: Retweeting filtering visualization (From left to right:retweetingN umber ≥ 0, retweetingN umber ≥ 1 and retweetingN umber ≥ 5 a tt h efi n a lt i m es t e p )

Figure 10 :

 10 Figure 10: Lens: Lift the focused items up and keep others in the original view plane. (From left to right: threshold retweeting number is 5, 15, and 50, respectively)

Figure 11 :

 11 Figure 11: Visualization results by clustering retweeting groups and show different propagation coverage with different depths

Figure 12 :

 12 Figure 12: Statistics of the first retweeting propagation dataset

Figure 13 :

 13 Figure 13: Hierarchy aware node-link diagram and some interactive techniques

Figure 14 :

 14 Figure 14: Link + View: Transform the fetched information from the node-link diagram to the depthvarying parallel data model and emphasize the highly focused items with blue lines

Figure 15 :

 15 Figure 15: Propagation section 1 (Retweeting time displaying and related web users and their retweeting message)

Figure 16 :

 16 Figure 16: Relevant web users and their retweets

17

 17

  .

Figure 17 :

 17 Figure 17: From left to right: node-link diagram and depth-varying parallel data visualization

Figure 18 :

 18 Figure 18: Simplified layout with displaying the identification of each vertex and the relevant patterns of them

Figure 19 :

 19 Figure 19: Propagation splitting and instantly retweeting behaviors presentation

Figure 20 :

 20 Figure 20: Hough scatter plot and time-varying parallel data visualization

Figure 21 :

 21 Figure 21: Quantifiable indicators to compare the performance between the previously used analysis tool and our visualization system

ACKNOWLEDGMENTS

This work was supported by Tsinghua-Sina Project, National Science Foundation of China(61272225,51261120376) and Chinese 863 Program(2012AA041606,2012AA040902).


Hai-Chuan Song

Jun-Hai Yong

H.-C Song

J.-H Yong

H.-C Song

H.-C Song

J.-H Yong

H.-C Song

J.-H Yong

Calculating Jacobian coefficients of primitive constraints with respect to Euler parameters Yo n g L i u •

Keywords: Assembly constraint, Primitive constraint, Jacobian coefficient, Virtual rotation, Variational method

It is a fundamental problem to calculate Jacobian coefficients of constraint equations in assembly constraint solving because most approaches to solving an assembly constraint system will finally resort to a numerical iterative method that requires the first-order derivatives of the constraint equations. The most-used method of deriving the Jacobian coefficients is to use virtual rotation which is originally presented to derive the equations of motion of constrained mechanical systems. However, when Euler parameters are adopted as the state variables to represent the transformation matrix, using the virtual rotation will yield erroneous formulae of Jacobian coefficients. The reason is that Euler parameters are incompatible with virtual rotation. In this paper, correct formulae of Jacobian coefficients of geometric constraints with respect to Euler parameters are presented in both Cartesian coordinates and relative generalized coordinates. Experimental results show that our proposed formulae make Newton-Raphson iterative method converge faster and more stable.

Introduction

Assembly modeling plays an important role in product design activities because most products are composed of a number of parts. Instead of designating the position and orientation of each part, a designer usually specifies geometric constraints or kinematic constraints between parts, which involves solving those constraints to obtain the postures of parts [START_REF] Kim | An assembly modelling system for dynamic and kinematic analysis[END_REF][START_REF] Anantha | Assembly modeling by geometric constraint satisfaction[END_REF][3]. Numerous approaches to solving 2D or 3D geometric constraint problems have been proposed in recent decades. The main idea of them is to split a large geometric constraint problem into smaller and easier ones using graph-based decomposing algorithm, rule-based strategy, algebraic theory, or numerical methods. Except for a small set of subproblems that can be solved analytically [START_REF] Kramer | A geometric constraint engine[END_REF][START_REF] Kim | Solving 3D geometric constraints for assembly modelling[END_REF], most of them will resort to numerical methods, e.g., a Newton-type iterative method or mathematical optimization method [START_REF] Zou | Computer-aided design using the method of cut-joint kinematic constraints[END_REF][START_REF] Li | A constructive approach to solving 3-D geometric constraint systems using dependence analysis[END_REF][START_REF] Kim | Solving 3D geometric constraints for closed-loop assemblies[END_REF][START_REF] Peng | A geometric constraint solver for 3-D assembly modeling[END_REF][START_REF] Shi | An angular constraints solving approach for assembly modeling based on spherical geometry[END_REF][START_REF] Xia | 3D geometric constraint solving using the method of kinematic analysis[END_REF].

To use a numerical method, the geometric constraints are converted into a set of nonlinear equations first. Those nonlinear equations are called primitive constraints or basic constraints. The selections of state variables can be classified into two categories: Cartesian coordinates and generalized coordinates. The Cartesian coordinate formulation yields a maximal set of highly sparse equations, but the form of the equations is independent on the topological structure of assembly systems. On the other hand, recursive generalized coordinate formulation can dramatically reduce the number of constraint equations and variables, but the form of the equations varies with the topological structure of assembly systems. In both cases, to compute the Jacobian coefficients of constraint equations, a variationalvector calculus approach is widely adopted in literature, which is originally presented for deriving the equations of motions in dynamics [START_REF] Haug | A variational-vector calculus approach to machine dynamics[END_REF][START_REF] Bae | A recursive formulation for constrained mechanical system dynamics: Part I. Open loop systems[END_REF][START_REF] Bae | A recursive formulation for constrained mechanical system dynamics: Part II. Closed loop systems[END_REF][START_REF] Haug | Computer aided kinematics and dynamics of mechanical systems: basic method[END_REF][16]. The approach introduces a concept called virtual rotation, with which, it is convenient to derive the Jacobian coefficients of primitive constraints with respect to the state variables. However, this concept should not be applied when Euler parameters are taken as the state variables, because Euler parameters break its prerequisites. Otherwise, taking the concept of the virtual rotation will yield wrong formulae of Jacobian coefficients, as presented in literature [START_REF] Zou | Computer-aided design using the method of cut-joint kinematic constraints[END_REF][START_REF] Peng | A geometric constraint solver for 3-D assembly modeling[END_REF]. In this paper, the prerequisites of virtual rotation are discussed. The correct Jacobian coefficients in Cartesian coordinates are derived using variational-vector calculus method. Moreover, a recursive formulation of Jacobian coefficients of constraint equations in relative generalized coordinates is also derived, which can be applied to any type of relative generalized coordinates, compared with the restricted formulation presented in [START_REF] Zou | Computer-aided design using the method of cut-joint kinematic constraints[END_REF].

The rest of this paper is organized as follows. In Section 2, we first generalize primitive constraints in literature to a formal representation. The reason of the wrong formulae presented in literature is discussed. Then the correct formulae of Jacobian coefficients with respect to Euler parameters in both Cartesian coordinates and recursive generalized coordinates are derived. In Section 3, four examples are provided to illustrate the correctness and efficiency of presented formulae. Finally, conclusions are made in Section 4.

Problem statement

In an assembly system, the posture of a rigid body i, could be identified by r i and A i ,wherer i is the origin and A i is the transformation matrix from the local reference frame to the global reference frame. The values in the global reference frame, of a vector v i and point P i fixed on body i can be obtained from the following transformations

v i = A i v ′ i P i = r i + v P i = r i + A i v P ′ i ,
where the symbols with a superscript ′ denote the values of the vector or point in the local reference frame. There exists an orthogonal condition

A i A T i = I (1)
and then the variation δA i satisfies

δA i A T i =-A i δA T i (2)
A tilde (∼) operator on a vector v = v x v y v z T ,forms

a skew symmetric matrix ṽ = ⎡ ⎣ 0 -v z v y v z 0 -v x -v y v x 0 ⎤ ⎦ .F r o m
Eq. 2,thevirtual rotation δπ i is defined by δ π i = δA i A T i and δπ ′ i = A T i δπ i . First, we generalize primitive constraints introduced in [START_REF] Haug | Computer aided kinematics and dynamics of mechanical systems: basic method[END_REF] to the following formal definition.

Definition 1 A primitive constraint is a scalar function f = f(r i 1 , r i 2 , ••• , r i M , v j 1 , v j 2 , ••• , v j N )
,w h e r er i k denotes the origin of body i k ,andv j k be any vector fixed on body j k .

In this paper, we use this formal definition to illustrate the derivations of Jacobian coefficients of primitive constraints, instead of a concrete constraint. The derivations of Jacobian coefficients with respect to both Cartesian coordinates and relative generalized coordinates are presented.

Cartesian coordinates and Euler parameters

A spatial mechanism, or assembly, constructed using Cartesian coordinates and Euler parameters, is defined in terms of the state variables q = r 1 p 1

p i = 1. Let e i = e 1i e 2i e 3i
T , E i = -e i , ẽi + e 0i I ,a n dG i = -e i , -ẽ i + e 0i I . A i can be expressed as

A i = E i G T i .I n [15], δπ ′ i is derived as δπ ′ i = 2G i δp i (3)
Using the virtual rotation, Peng et al. [START_REF] Peng | A geometric constraint solver for 3-D assembly modeling[END_REF] derived the variation of a vector v i as

δv i = δA i v ′ i = A i δ π ′ i v ′ i =-A i ṽ′ i δπ ′ i =-2A i ṽ′ i G i δp i . (4)
Then the variation of the function f is derived as

δf = k ∂f ∂r i k δr i k + k ∂f ∂v j k δv j k = k ∂f ∂r i k δr i k -2 k ∂f ∂v j k A j k ṽ′ j k G j k δp j k (5)
from which the Jacobian coefficients of f are obtained as

f r i k = ∂f ∂r i k ,k = 1, ••• ,M (6a)
f p j k =-2 ∂f ∂v j k A j k ṽ′ j k G j k ,k = 1, ••• ,N, (6b)
where the partial derivatives ∂f ∂r i k and ∂f ∂v j k are only determined by the form of f itself and easy to compute.

As is pointed out by Yen et al. [START_REF] Yen | Automatic generation and numerical integration of differential-algebraic equations of multibody dynamics[END_REF], however, Eq. 6b is not the true derivative of f with respect to Euler parameters p j k , because the definition of the virtual rotation implies the establishment of Eqs. 1 and 2, which requires p T i p i = 1and p T i δp i = 0. It results in a contradiction. Yen et al. presented a formula of the variation of a vector v i to compute the correct Jacobian coefficients of primitive constraints, but his formula is complicated. In this paper, using the variationalvector calculus method, we derive this variation in a more concise form as follows.

δv i = δ A i v ′ i = δ E i G T i v ′ i = δ e i e T i + ẽi ẽi + 2e 0i ẽi + e 2 0i I v ′ i = δ e i e T i v ′ i + δ ẽi ẽi v ′ i + δ 2e 0i ẽi v ′ i + δ e 2 0i v ′ i = δe i e T i v ′ i + e i δe T i v ′ i + δẽ i ẽi v ′ i + ẽi δẽ i v ′ i + 2δe 0i ẽi v ′ i + 2e 0i δẽ i v ′ i + 2e 0i v ′ i δe 0i = e T i v ′ i δe i + e i v ′ T i δe i -ẽi v ′ i δe i -ẽi ṽ′ i δe i + 2ẽ i v ′ i δe 0i -2e 0i ṽ′ i δe i + 2e 0i v ′ i δe 0i = 2 ẽi + e 0i I v ′ i δe 0i + e T i v ′ i I + e i v ′ T i -ẽi v ′ i -ẽi ṽ′ i -2e 0i ṽ′ i δe i = K v ′ i , p i δp i , (7)
where

K(u, p) = 2(ẽ + e 0 I)ue T uI + eu T -ẽ ũ -ẽu -2e 0 ũ .
Substituting Eq. 7 into Eq. 5 yields

f p j k = ∂f ∂v j k K v ′ j k , p j k ,k = 1, ••• ,N. (8)

Relative generalized coordinates

Using relative generalized coordinates, r i and A i can be represented recursively instead. The relation of a pair of coupled bodies, bodies i and j , is depicted in Fig. 1. Vectors that locate joint attachment points in bodies i and j are denoted by s ij and s ji , respectively. Orthogonal matrices C ij , C ji ,andA ′′ ij are transformations from the joint definition frames ( ′′ ) to the body frames( ′ ) on bodies i and j and from the joint definition frame on body j to the joint definition frame on body i, respectively. From Fig. 1,wehavethe following relations:

r j = r i + s ij + d ij -s ji (9a) A j = A i C ij A ′′ ij C T ji (9b)
and

d ij = A i C ij d ′′ ij (10a)
Fig. 1 The representation of relative motions of two bodies

s ij = A i s ′ ij ( 10b
)
s ji = A j s ′ ji (10c)
where s ′ ij and s ′ ji are fixed vectors on each body frame separately. In the equations above, A ′′ ij and d ′′ ij are only the functions of relative generalized coordinates q ij between bodies i and j .

A mechanical system can be represented by a graph G = G(V , E), in which each node of V represents a rigid body, and each edge e ∈ E represents the joint between two bodies. Both an open-loop and closed-loop mechanical systems can be abstracted as a topological tree structure with additional geometric constraints. Consider any body i. There exists a unique path 0 → ••• → i -1 → i from base body 0 to body i. Using Eqs. 9 and 10, r i and A i can be represented recursively by the relative generalized coordinates q 0,1 , q 1,2 , •••, q i-1,i [START_REF] Bae | A recursive formulation for constrained mechanical system dynamics: Part I. Open loop systems[END_REF].

Let q denote the collection of all relative generalized coordinates. To obtain the Jacobian coefficients of f ,t h e variations δr i and δv i should be derived first. Let δr i and δv i have the following forms

δr i = J r i δq (11a)
δv i = J v i δq (11b)
Therecursi veformulaeofJ r i and J v i were derived using the virtual rotation in [START_REF] Zou | Computer-aided design using the method of cut-joint kinematic constraints[END_REF][START_REF] Xia | 3D geometric constraint solving using the method of kinematic analysis[END_REF]. However, similarly, when a free rotational joint is involved and Euler parameters are chosen as the relative generalized coordinates, the virtual rotation should not be applied. Otherwise, using the virtual rotation will yield erroneous formulae. In this paper, we derive the recursive formulae directly instead. First, two auxiliary coefficient matrices are introduced.

Definition 2 Let u be any 3D vector. The auxiliary coefficient matrices J i-1,i (u) and J i (u) satisfy J i-1,i (u)δq i-1,i = δA ′′ i-1,i u and J i (u)δq = δA i u separately.

From the definition, J i-1,i is only determined by relative generalized coordinates q i-1,i and can be deduced from the joint definition. For example, if q i-1,i is the Euler parameters, we can obtain its formula from Eq. 7.

J i-1,i (u)δq i-1,i = δA ′′ i-1,i u = K(u, q i-1,i )δq i-1,i (12)
TherecursiveformulaofJ i (u) is derived as follows

J i (u)δq = δA i u = δ A i-1 C i-1,i A ′′ i-1,i C T i,i-1 u = δA i-1 C i-1,i A ′′ i-1,i C T i,i-1 u + A i-1 C i-1,i δA ′′ i-1,i C T i,i-1 u = J i-1 C i-1,i A ′′ i-1,i C T i,i-1 u δq + A i-1 C i-1,i J i-1,i C T i,i-1 u δq i-1,i = J i-1 C i-1,i A ′′ i-1,i C T i,i-1 u δq + 0A i-1 C i-1,i J i-1,i C T i,i-1 u 0 δq from which we obtain J i (u) = J i-1 C i-1,i A ′′ i-1,i C T i,i-1 u + 0A i-1 i-1,i J i-1,i C T i,i-1 u 0 . ( 13
)
The variation δv i can be obtained as

δv i = δA i v ′ i = J i v ′ i δq (14)
from which we obtain

J v i = J i v ′ i . (15)
The variation r i is derived as follows

δr i = δ r i-1 + s i-1,i + d i-1,i -s i,i-1 = δr i-1 + δA i-1 s ′ i-1,i + δ A i-1 C i-1,i d ′′ i-1,i -δA i s ′ i,i-1 = δr i-1 + J i s ′ i-1,i δq + J i-1 C i-1,i d ′′ i-1,i δq + A i-1 C i-1,i δd ′′ i-1,i -J i s ′ i,i-1 δq = δr i-1 + J i s ′ i-1,i δq + J i-1 C i-1,i d ′′ i-1,i δq + A i-1 C i-1,i D i-1,i δq i-1,i -J i s ′ i,i-1 δq (16)
from which we obtain

J r i = J r i-1 + J i s ′ i-1,i + J i-1 C i-1,i d ′′ i-1,i + 0A i-1 C i-1,i D i-1,i 0 -J i s ′ i,i-1 ,
where D i-1,i is the Jacobian coefficient of d ′′ i-1,i with respect to local relative generalized coordinates q i-1,i . Substituting Eqs. 14 and 16 into Eq. 5 yields

δf = k ∂f ∂r i k δr i k + k ∂f ∂v j k δv j k = k ∂f ∂r i k J r i k δq + k ∂f ∂v j k J i v ′ j k δq
from which, the Jacobian matrix of f is obtained as

f q = k ∂f ∂r i k J r i k + k ∂f ∂v j k J i v ′ j k . (17)

Jacobian coefficients of primitive constraints

Based on Haug's work [START_REF] Haug | Computer aided kinematics and dynamics of mechanical systems: basic method[END_REF], some primitive geometric constraints were presented to build the geometric constraint library [START_REF] Xia | 3D geometric constraint solving using the method of kinematic analysis[END_REF]. Most important primitive constraints are dot-1 constraint, dot-2 constraint, angle constraint, and distant constraint.

d1 (v i , v j ) = v T i v j = 0 (18a) d2 (v i , d ij ) = v T i d ij = v T i d ij = 0 (18b) ang (v i , v j ) = v T i v j -cos α 0 = 0 (18c) dist (d ij ) = d T ij d ij -d 2 0 = 0, ( 18d
)
where

d ij = Q j -P i = r j + A j s Q ′ j -r i -A i s P ′ i .
The Jacobian coefficients of these constraints are derived as follows. For Cartesian coordinates, those are

d1 r i = d1 r j = 0 (19a) d1 p i = v T j K v ′ i , p i (19b) d1 p j = v T i K v ′ j , p j (19c) d2 r i =-v T i (19d) d2 r j = v T i (19e) d2 p i = d T ij K v ′ i , p i -v T i K s p ′ i , p i (19f) d2 p j = v T i K s p ′ j , p j (19g) dist r i =-2d T ij (19h) dist r j = 2d T ij (19i) dist p i =-2d T ij K s p ′ i , p i (19j) dist p j = 2d T ij K s p ′ j , p j , (19k)
and for relative generalized coordinates,

d1 q = v T j J i v ′ i + v T i J j v ′ j (20a) d2 q = d T ij J i v ′ i + v T i J r j + J j s p ′ j -J r i -J i s p ′ i (20b) dist q = 2d T ij r j + J j s p ′ j -J r i -J i s p ′ i (20c)
In the next section, experimental results will be presented to show that the formulae presented in this paper are more feasible than those proposed in [START_REF] Peng | A geometric constraint solver for 3-D assembly modeling[END_REF].

Experimental results

In this section, we use the Newton-Raphson iteration method listed in Algorithm 1 to evaluate the effectiveness of the presented formulae. Four examples are used to illustrate the advantages of our proposed formulae over those presented by Peng et al. [START_REF] Peng | A geometric constraint solver for 3-D assembly modeling[END_REF]. Without loss of generality, all the examples are two-body constraint systems with fixed origins so that each of them only has three rotational degrees of freedom left. First, we use Example 1 to demonstrate that the iteration using Eq. 6b breaks the assumption of orthogonality of the transformation matrix, and then use these four examples to illustrate that our formulae are more applicable.

Example 1 Let origins of bodies 1 and 2 be coincident. There are three dot-2 constraints ( d2 (v i , P i , Q i ), i = 1 •••3) between the bodies. This constraint system is well constrained.

With the parameters given in Table 1, this constraint system has real solutions. The unknowns of the constraint system is the q = p 2 = e 0 e 1 e 2 e 3 T . The equation set of the constraint system is represented by

(q) = ⎧ ⎨ ⎩ q T q -1 v T i (Q i -P i ) ,i = 1 •••3
Algorithm 1: Newton iteration procedure input : The initial estimate q 0 the set of equations output: The solved status q ← q 0 n ← 0 while n<Nand (q) <εdo J ← q (q) if J is invertible then q ← q -J -1 (q) end else / * Use Penrose-Moore inverse instead * / q ← q -J + (q) end n ← n + 1 end if n == N then return Failed end return Success of which the Jacobian matrix q is derived from Eqs. 19f and 19g as

q (q) = ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ 2q T (Q 1 -P 1 ) T K v ′ 1 , q -v T 1 K P ′ 1 , q (Q 2 -P 2 ) T K v ′ 2 , q -v T 2 K P ′ 2 , q (Q 3 -P 3 ) T K v ′ 3 , q -v T 3 K P ′ 3 , q ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ (21)
while Eq. 6b gives

ˆ q (q) = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 2q T 2 v ′ T 1 P′ 1 -(Q 1 -P 1 ) T Aṽ ′ 1 G 2 v ′ T 2 P′ 2 -(Q 2 -P 2 ) T Aṽ ′ 2 G 2 v ′ T 3 P′ 3 -(Q 3 -P 3 ) T Aṽ ′ 3 G ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ . (22)
Given an initial estimate q 0 = 0.5 -0.5 -0.50 .5 T , which satisfies the normalization constraint, we can get

A = ⎡ ⎣ 00-1 10 0 0 -10 ⎤ ⎦ and AA T = I,
Substituting A into Eq. 22 yields

ˆ q 0 = ⎡ ⎢ ⎢ ⎣

Tabl e 1 The parameters used in examples

Example Parameters

1 v ′ 1 = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 1 √ 5 1 √ 5 1 √ 5 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , v ′ 2 = ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ 0 2 √ 5 1 √ 5 ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ , v ′ 3 = ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ 1 √ 2 1 √ 2 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ , P ′ 1 = ⎡ ⎢ ⎣ 4 4 -2 ⎤ ⎥ ⎦, P ′ 2 = ⎡ ⎢ ⎣ 10 0 7 ⎤ ⎥ ⎦, P ′ 3 = ⎡ ⎢ ⎣ 0 3 10 ⎤ ⎥ ⎦, Q ′ 1 = ⎡ ⎢ ⎣ 0 5 1 ⎤ ⎥ ⎦, Q ′ 2 = ⎡ ⎢ ⎣ 3 1 -1 ⎤ ⎥ ⎦, Q ′ 3 = ⎡ ⎢ ⎣ 1 2 2 ⎤ ⎥ ⎦ 2 P ′ 1 = ⎡ ⎢ ⎣ 0 0 0 ⎤ ⎥ ⎦, P ′ 2 = ⎡ ⎢ ⎣ 1 0 0 ⎤ ⎥ ⎦, P ′ 3 = ⎡ ⎢ ⎣ 0 1 0 ⎤ ⎥ ⎦, Q ′ 1 = ⎡ ⎢ ⎣ 0 0 1 ⎤ ⎥ ⎦, Q ′ 2 = ⎡ ⎢ ⎣ 2 1 0 ⎤ ⎥ ⎦, Q ′ 3 = ⎡ ⎢ ⎣ 1 0 2 ⎤ ⎥ ⎦, d 1 = 1, d 2 = √ 6, d 3 = √ 5 3 u ′ 1 = ⎡ ⎢ ⎣ 1 0 0 ⎤ ⎥ ⎦, u ′ 2 = ⎡ ⎢ ⎣ 0 1 0 ⎤ ⎥ ⎦, v ′ 1 = ⎡ ⎢ ⎣ 1 0 0 ⎤ ⎥ ⎦, v ′ 2 = ⎡ ⎢ ⎣ 0 1 0 ⎤ ⎥ ⎦, u ′ 3 = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 1 √ 3 1 √ 3 1 √ 3 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , P ′ 1 = ⎡ ⎢ ⎣ 0 1 1 ⎤ ⎥ ⎦, Q ′ 1 = ⎡ ⎢ ⎣ 1 0 2 ⎤ ⎥ ⎦ 4 u ′ 1 = ⎡ ⎢ ⎣ 1 0 0 ⎤ ⎥ ⎦, u ′ 2 = ⎡ ⎢ ⎣ 0 1 0 ⎤ ⎥ ⎦, v ′ 1 = ⎡ ⎢ ⎣ 1 0 0 ⎤ ⎥ ⎦, v ′ 2 = ⎡ ⎢ ⎣ 0 1 0 ⎤ ⎥ ⎦, P ′ 1 = ⎡ ⎢ ⎣ 1 -1 0 ⎤ ⎥ ⎦, Q ′ 1 = ⎡ ⎢ ⎣ 0 1 2 ⎤ ⎥ ⎦, α 1 = π 6 , α 2 = π 3 , d 1 = 2
From Algorithm 1, we can get the next iteration point q 1 = 0.3750 -0.6250 -0.8750 0.1250

T , and the trans- The difference between these two regions is that the latter is the set where the initial estimate satisfies the normalization condition.

formation matrix A = ⎡ ⎣ -0.
Other used examples are listed as follows:

Example 2 The constraint system consists of three distant constraints ( dist (P i , Q i ,d i ), i = 1 •••3).

Example 3 The constraint system consists of two dot-1 constraints ( d1 (u i , v i ), i = 1, 2) and one dot-2 constraint ( d2 (u 3 , P 1 , Q 1 )).

Example 4 The constraint system consists of two angle constraints ( ang (u i , v i ,α i ), i = 1, 2) and one distant constraint ( dist (P 1 , Q 1 ,d 1 )).

The parameters used in the examples are given in Table 1. The simulation results are listed in Table 2. The results show that the iteration methods adopting our presented formulae as Jacobian coefficients are more stable. As depicted in Table 2, there are more convergent cases in general, using proposed formulae. Moreover, our proposed formulae will make the iteration method converge faster in fewer iterations. The experimental results demonstrate the feasibility of the method presented in this paper.

Conclusion

In this paper, the erroneous derivation of Jacobian coefficients of primitive constraint equations with respect to Euler parameters presented in literature is discussed, and correct formulae are derived using variational-vector calculus method. A recursive formulation of Jacobian coefficients of primitive constraints in relative generalized coordinates is also proposed. Compared with the formulae presented in literature, this formulation can be applied to all types of joint variables. Experimental results are presented to illustrate the correctness and computational efficiency of the formulae presented in this paper.

  Continuing this procedure, we will find that the iteration fails to converge. On the other hand, if we adopt Eq. 21, the iteration converges in about 10 times given a tolerance 10 -14 .To study the convergency generally, many initial estimates are generated to test whether Algorithm 1 converges, and if it does, in how many iterations. Our tests are For each example, we generate 10, 000 initial estimates uniformly at random from an estimate region to test whether Algorithm 1 converges, using our and compared formulae separately. We count the convergent cases and calculate the average iterations needed for each successful case. The two estimate regions to concern are 1. {(e 0 ,e 1 ,e 2 ,e 3 )| |e i | ≤ 1} (denoted as hyper box) 2. {q 0 |q T 0 q 0 = 1} (denoted as hyper sphere)

	performed as follows.
	which breaks Eq. 1. Tabl e 2 The simulation results	2500 1.0000 -0.8125 1.1875 0.5000 0.2500 ⎤ ⎦ 0.5000 -0.6875 -1.0000
		Experiment	Estimate region	Jacobian formulae	Convergent cases	Average iterations
		Exam. 1,	hyper box	Proposed	10,000	9.43
		body 2 movable		Peng's	9,919	23.38
			hyper sphere	Proposed	10,000	10.65
				Peng's	9,824	27.19
		Exam. 1,	hyper box	Proposed	3,439	10.98
		body 1 movable		Peng's	30	12.77
			hyper sphere	Proposed	2,666	10.07
				Peng's	238	10.83
		Exam. 2	hyper box	Proposed	8,120	8.92
				Peng's	651	58.74
			hyper sphere	Proposed	6,446	7.42
				Peng's	1,252	62.76
		Exam. 3	hyper box	Proposed	1,592	8.88
				Peng's	1,681	13.65
			hyper sphere	Proposed	1,311	8.49
				Peng's	1,531	13.63
		Exam. 4	hyper box	Proposed	8,250	8.42
				Peng's	7,748	10.94
			hyper sphere	Proposed	8,547	8.43
				Peng's	8,082	10.99

Acknowledgments

The research was supported by Chinese 973 Program (2010CB328001) and Chinese 863 Program (2012AA040902). The first author was supported by the NSFC(61035002, 61173077). The second author was supported by the NSFC (61063029). The last author was supported by the NSFC (61272235).


Yi Peng
email: 15pengyi@gmail.com

Li Chen

Jun-Hai Yong
email: yongjh@tsinghua.edu.cn

Importance-Driven Isosurface Decimation for Visualization of Large Simulation Data Based on OpenCL

Keywords: Isosurface Visualization, Importance Driven, Transfer Function, Marching Cubes, Mesh Simplification, OpenCL ✦

For large simulation data, Parallel Marching Cubes algorithm is efficient and commonly used to extract isosurfaces in 3D scalar field. However, the isosurface meshes are sometimes too dense and it is difficult for scientists to specify the areas they are interested in. In this paper, we provide them a new way to define mesh importance for decimation using transfer functions and visualize large simulation data in case the normal visualization methods cannot handle due to memory limit. We also introduce a parallel isosurface simplification framework which uses pyramid peeling to extract the decimated meshes progressively without generating the original surface. Since the implementation uses OpenCL which is oriented to heterogeneous computing, our method can be applied to different parallel systems and scientists can see the visualization results while doing simulations. Finally, we evaluate the performances of our algorithm and use different scientific datasets to show the efficiency of our method.

there are many mesh simplification methods such as edge collapse, vertex removal and polygon merging [START_REF] Luebke | A developer's survey of polygonal simplification algorithms[END_REF] which can preserve geometric features. It is still difficult for users to define mesh importance according to their interests. In order to solve this problem, we propose an importance-driven mesh decimation strategy, which provides users an interactive way to specify the importance of meshes according to features using transfer functions and allow them to directly see the feature distributions and evolutions on the decimated surfaces. The most important contributions we have made are listed below:

1) Using transfer functions to define the importance of meshes for isosurface generation interactively; 2) Simplifying isosurface before it is generated; 3) The entire algorithm is parallelized based on OpenCL which includes feature extraction, mesh generation and model rendering;

Our method consists of three parts. Data Pre-processing classifies the cubes and extracts user-specified features. Pyramid Construction builds the histograms and feature pyramids which contain decimation information. Pyramid Peeling generates and renders the decimated meshes in a top-down way. These meshes are generated while being patched and refined. We also test the performance of our algorithm and use several datasets including medical datadsets and scientific datasets to evaluate our method. The results show that our method is efficient and interactive.

The remainder of this paper is organized as follows. Section 2 introduces the related work. Section 3 defines the mesh importance which is used as the decimation criteria. Section 4 gives a detail introduction of our algorithm which is highly parallelized. Section 5 shows the results of our importance-driven decimation strategy based on several datasets. Section 6 summarizes the advantages of our method and outlines the future work.

2R ELATED WORK

The research related to our method can be divided into three groups. First, our idea comes from importance-driven visualization method and this work is based on some latest research in Parallel Marching Cubes. We also benefits from traditional mesh simplification methods.

Importance-Driven Visualization Method

In scientific visualization, there are two ways to define data importance. The feature-based methods need to extract features from the original data while the visibility-based methods pay attention to user perceptions.

For feature-based methods, researchers focus on texture, size, importance curves for featuretemporal analysis [START_REF] Wang | Importance-driven timevarying data visualization[END_REF], shape [START_REF] Peng | Transfer function design using shape features in medical visualization[END_REF] and so on. They also do mesh deformation to reduce volumetric data based on voxel importance, which provides users a focus + context visualization [START_REF] Wang | Featurepreserving volume data reduction and focus+context visualization[END_REF].

From the point of users, visibility is introduced by Bordoloi et al. [START_REF] Bordoloi | View selection for volume rendering[END_REF]. Further research computes the visibility histograms [START_REF] Correa | Visibility histograms and visibility-driven transfer functions[END_REF] and informational divergence [START_REF] Ruiz | Automatic transfer functions based on informational divergence[END_REF] which is used to maximize the user-defined importance.

In our method, we use geometric features which are extracted from the original data and allow users to specify mesh importance using transfer functions. So users can interactively define the mesh importance and directly see the decimated meshes and the feature distributions on that isosurface.

Marching Cubes on GPU

Marching Cubes is introduced by Lorensen et al. [START_REF] Lorensen | Marching cubes: A high resolution 3d surface construction algorithm[END_REF] and there are two different ways to parallelize it. On one hand, Goetz et al. [START_REF] Goetz | Real-time marching cubes on the vertex shader[END_REF] extends Marching Cubes using vertex shader which is compatible to OpenGL fix function. After that, the work is improved by introducing preclassification and histogram pyramid. On the other hand, Smistad [START_REF] Smistad | Fast surface extraction and visualization of medical images using opencl and gpus[END_REF] implements the original algorithm in OpenCL, so that the algorithm runs entirely on the GPU, which means their work is oriented to general computing.

Our work is based on OpenCL since we want to unify the computing process and rendering process. So our visualization techniques can be easily integrated into different research fields.

Mesh Simplification

For traditional mesh simplification methods, incremental decimation methods such as edge collapse, vertex collapse, vertex removal and polygon merging are difficult for parallelism due to data sorting [START_REF] Luebke | A developer's survey of polygonal simplification algorithms[END_REF]. Only vertex clustering can be easily parallelized [START_REF] Vetter | Scalability issues of in-situ visualization in parallel simulation of unsteady flows[END_REF].

For marching cubes surfaces, early work reads adjacent slices of volumetric data and merges the mesh vertices or uses octree to simplify the main continuous surface [START_REF] Shekhar | Octreebased decimation of marching cubes surfaces[END_REF]. Attali et al. [START_REF] Attali | Extraction and simplification of iso-surfaces in tandem[END_REF] introduce tandem algorithm which extends the Marching Cubes by sweeping the dataset and extracting meshes on the accumulated partial surface. Dupuy et al. [START_REF] Dupuy | [END_REF] extend this algorithm to a parallel system. However, these methods are hard to be implemented on GPU since they are streamed.

In this paper, we introduce a parallel decimation method based on mesh patching and refining, which benefits from crack patching [START_REF] Shekhar | Octreebased decimation of marching cubes surfaces[END_REF] and vertex clustering.

3I MPORTANCE DEFINITION

Before mesh simplification, we use both geometric properties and user-specified features to define mesh importance. The geometric properties are mainly based on gradient which is also used to control the quality of meshes. The userspecified features can interactively be changed through transfer function editor and be applied to the decimation process, so that the generated meshes are important to users.

Geometric Properties

For isosurface, three most important geometric properties are continuity, smoothness and topology. In this paper, we consider continuity and smoothness because the topology is changed during simplification.

Since Marching Cubes algorithm guarantees mesh quality for a same resolution, the main challenge is to maintain the continuity between two different resolutions. We follow the idea of [START_REF] Shekhar | Octreebased decimation of marching cubes surfaces[END_REF] which performs crack patching in a top-down way. During this process, we use a histogram pyramid(see [START_REF] Smistad | Fast surface extraction and visualization of medical images using opencl and gpus[END_REF]) like an octree to record the patching information, so the finer meshes generated by a child cube are consistent to the coarser ones correspond to its parent cube. The parent-children relationship of cubes is defined by their relationship in the octree, which means a parent cube has 8 child cubes at the next resolution(finer resolution).

On the other hand, we choose the smoothness of gradient as one of the decimation criteria. For each cube c, we first compute its gradient using central difference and then calculate the smoothness feature s c according to the gradients of its 8 child cubes as follow:

s (l) c =1- ⎛ ⎜ ⎜ ⎜ ⎜ ⎝ 1+ N i=1 n c i g (l-1) c i 2 N i=1 n c i - N i=1 n c i g (l-1) c i N i=1 n c i 2 ⎞ ⎟ ⎟ ⎟ ⎟ ⎠ -1
(1) where for volumetric data, N =8which is the total number of child cubes of a parent cube, l is the resolution level, g c denotes normalized gradient vector and and n c i denotes the number of meshes generated by cube c.I n implementation, for each cube we only store two gradient moments and one mesh number. So the cubes with lower smoothness are more likely to be decimated.

User-Specified Features

For user-defined features such as texture, pressure, velocity and so on, we first construct a feature pyramid in a same way as smoothness feature:

f (l) c = N i=1 n c i f (l-1) c i N i=1 n c i (2)
where N =8see Equation. 1, f c denotes the features which are domain-specific derivatives.

After that the features are mapped to importance values by transfer functions. So users can define mesh importance according to the features they are interested in interactively. In this paper, we provide users a one-dimensional transfer function editor to design the mapping functions as shown in Fig. 1. However the mapping can easily be extended to multidimensional, if different domains are considered. In this transfer function editor, users can edit the control points to generate the mapping functions. For each control point, the horizontal coordinate indicates the corresponding feature which is obtained by blending geometric feature with user-specified features. The vertical coordinate and user-specified color denotes the importance and color of that feature respectively.

The advantages of using transfer functions to define mesh importance are obvious. First, this kind of interactions are widely used in scientific visualization. Also, it is easy to learn. For one-dimensional transfer function, users only need to add(delete) control points or just modify their positions and colors. Then the system will generate the mapping functions automatically. More importantly, the background histograms of features help users design the mapping functions, since they connect feature with color and mesh importance.

Mesh Importance

The final importance of a mesh combines the geometric properties and user-defined features. For each cube c, the importance I c is defined as follow:

I (l) c = t (1 -α) s c (l) + α •|f c (l) | ( 3
)
where α is a blend factor and t is the transfer function specified by user. And in the next stage, the importance is used to decimate the meshes.

4M ESH SIMPLIFICATION

Based on [START_REF] Smistad | Fast surface extraction and visualization of medical images using opencl and gpus[END_REF]'s work, our mesh simplification method is highly parallelized. Also we have made several improvements on it, such as lower memory cost and progressive rendering support.

Simplification Framework

The framework of our parallel simplification algorithm includes 3 stages as shown in Fig. 2. All these stages are carried out in GPU.

The first stage computes both the geometric and user-defined features. After that, the construction stage build histogram and feature pyramids while recording the decimation information. The final stage performs a topdown mesh peeling. In the last stage, the meshes are rendered shortly after they have been assembled which means they are generated progressively. So that users can choose which resolution should be rendered. The advantages of decimating meshes when generating are less memory cost since we only need to store decimation information and final meshes instead of the full resolution meshes, faster rendering due to less meshes and progressive exploration.

Decimation Criteria

With mesh importance the decimation criteria can be defined as follow:

d (l) c = 1 I (l) c ≤ β 0 I (l) c >β (4)
where β is a threshold used to control the degree of mesh decimation. When a parent cube is decimated, the same operation should be done for all the descendants. However, since for each resolution we need to compute the number of meshes, which means we will do a bottom-up scan, so we only need to record the decimation information while scanning the pyramid.

Data Pre-processing

For cube classification, we use the infrastructures provided by [START_REF] Smistad | Fast surface extraction and visualization of medical images using opencl and gpus[END_REF]. Although it leads to topological ambiguity, we find that the performances are better to be measured when using the same look up table. Meanwhile, we also compute the gradient feature and extract userdefined features during the same time.

Pyramid Construction

In this stage, histogram and feature pyramids are constructed in a bottom-up way, so that cubes with higher resolution are processed earlier. For each resolution, we use smoothness and user-define features to compute the imof cubes. After that, the importance is used to decide whether this cube should be decimated. If it is true, the meshes belong to this cube will be updated, which means this cube will be chosen as a representative of all its descendants.

Pyramid Peeling

The last stage generates the meshes progressively, so the coarsest meshes are rendered first. For each resolution there is a corresponding peeling process which consists of three steps. The gathering step traverses the histogram pyramid to get the total number of meshes. The assembling step computes mesh positions, normals and color attributes. The last step just renders all the meshes. In the second step, there are three key techniques which are used to improve mesh quality.

Mesh Refinement

For meshes in different resolutions, cracks may be generated due to discontinuous meshes as shown in Fig. 3.

Since coarse meshes are prominent, we first improve their qualities by refining them according to detail information. That is to say, we need to re-compute mesh positions, normal vectors using the information of their descendants. In implementation, we use binary search to calculate the finest value in the full resolution.

Fig. 4 shows the results of mesh refinement. Comparing Fig. 4(a) with the original results shown in Fig. 3(a), the results are not that different.

Crack Patching

In order to patch the cracks, we re-compute the finer meshes according their neighboring ancestors, which means for all the vertices on a finer mesh, we compute their new properties. For a vertex, we first find its 8 neighboring cubes. Then all of their ancestors will be checked in a top-down way, so coarser cracks are patched first. If there are several ancestors which may cause cracks, we choose one at the coarsest resolution. After that, we project this vertex to the edges generated by that ancestor. If there are several edges which are coplanar with this vertex, we choose one with the shortest distance. Fig. 5 shows the results of crack patching.

From Fig. 5 we can see that, the crack patching technique gets remarkable results, since the boundary consistency of meshes in different

Progressive Rendering

We also render the meshes progressively. So users can choose which resolution to be rendered as shown in Fig. 6.

5R ESULTS

In this section, we analyze the performance and memory cost of our algorithm and use several datasets to test our method.

Algorithm Performance

The performance of our method is shown in TABLE1. From this table, we can see that the time cost is strongly related to the number of decimated meshes. Larger decimation threshold leads to better performance. These tests are done on a standard PC: CPU Intel(R) Core i7-2670QM @ 2.20GHz, RAM 2 * 4.00GB @ 1333MHz, GPU Nvidia GeForce GT 555M @ 2GB and OS Windos 7 Ultimate SP1. TABLE2 shows the memory costs of our method and Smistad's method [START_REF] Smistad | Fast surface extraction and visualization of medical images using opencl and gpus[END_REF]. Compared with their method, we use buffer objects instead of 3D textures and support progressive and multi-pass rendering which can deal with large data. Since their method relies on the maximum size of the data, they can not handle inhomogeneous data. Even though we need addition space to record decimation information and features, since our pyramid only relies on the total size of the data, we can still handle large datasets. If we only consider geometric property, the feature size will be four times the data size(three for gradient and one for mag- nitude), see "Head" dataset. If another scalar feature is used(see "Hurricane" dataset), this size can be five times the data size.

Mesh Decimation Results

We use both medical datasets and scientific datasets to evaluate our method. The evaluation includes two parts. The first part focuses on the geometric features, while the second part shows the decimation results combing geometric properties and user-specified features. Fig. 7 and Fig. 8 show the results of mesh simplification using only geometric features. That is to say, the blending parameter α is set to zero. In Fig. 7, the meshes around tooth have higher smoothness values, so they are preserved during simplification while the meshes on forehead are decimated. In Fig. 8, meshes which are far from the eyes are decimated since they are smoother than the central meshes. So that the eye of the hurricane is preserved during mesh decimation.

Fig. 9 and Fig. 10 show the results of mesh simplification using user defined features. In From these figures, we can see that our method preserves both geometric feature and user-defined features well while performing mesh simplification.

6C ONCLUSION

In this paper, we introduce an importancedriven parallel mesh simplification method which has three main advantages:

1) mesh importance definition based on transfer functions allows users to interactively see the distributions of features with their interests; 2) mesh decimation before generation costs less memory and leads to higher rendering performance, which also support progressive mesh generation; 3) oriented to general computing ensures that our method can be easily integrated into other parallel systems; We also evaluate the performance of our algorithm and use several datasets to show the efficiency of our method.

Other future work includes improving current method to eliminate topological ambiguity, enhancing the performance of pyramid peeling stage since it is not so satisfactory, integrating more features and designing better interface for users to design transfer functions.

Fig. 1 .

 1 Fig. 1. One-dimensional transfer function editor which is used to specify mesh importance. The height of a control point defines the feature importance while its color corresponds to the color mapping functions used in surface rendering. The background curves show the distributions of data and user-specified feature.

Fig. 2 .

 2 Fig. 2. Parallel Marching Cubes surfaces decimation framework. The framework includes three parts: pyramid initialization, construction and peeling. In the first stage, the cubes are classified and features including geometric properties and user-specified features are extracted. In the next construction stage, the isosurface meshes are generated while decimated, and in the last stage, meshes with different resolutions are rendered top-down progressively.

Fig. 3 .

 3 Fig. 3. Cracks generated by Marching Cubes Algorithm. (a) shows the original isosurface. (b) shows the meshes in different resolutions.

Fig. 4 .

 4 Fig. 4. Mesh refining results. (a) shows the results of mesh refinement technique. (b) draws the wireframes of the same model.

Fig. 5 .

 5 Fig. 5. Crack patching results. (a) shows the results of crack patching technique. (b) draws the wireframes of the same model.

Fig. 6 .

 6 Fig. 6. Progressive rendering results. (a) shows the meshes of all resolutions, in which higher lightness values are set to coarser meshes. (c) is their corresponding transfer functions. (b) shows the meshes at coarse resolutions and (d) is their transfer functions.

Fig. 9 ,Fig. 7 .

 97 Fig. 7. Results of Head dataset with geometric features.

Fig. 8 .

 8 Fig. 8. Results of Hurricane dataset with geometric features.

Fig. 9 .

 9 Fig. 9. Results of Hurricane dataset with user defined features.

Fig. 10 .

 10 Fig. 10. Results of Fluid dataset with user defined features.

TABLE 1

 1

		Algorithm Performance
	Dataset	Parameters	Mesh Number		Data Pre-	Pyramid Pyramid Total
	Name	Iso αβ Original Decimated % processing Construction Peeling Time
		0.25 0.00 0.0035 703196 175122 24.9	40ms	41ms	115ms 196ms
	Head	0.25 0.00 0.0130 703196	58868	8.37	39ms	41ms	63ms 143ms
	256x256x225	0.38 0.00 0.0040 827884 180752 21.83 40ms	41ms	97ms 178ms
		0.38 0.00 0.0120 827884	63606	7.68	39ms	40ms	45ms 124ms
		0.25 0.00 0.0015 612936 131598 21.47 35ms	36ms	67ms 138ms
	Hurricane(18 th ) 0.25 0.00 0.0055 612936	44860	7.32	35ms	35ms	44ms 114ms
	661x401x42	0.25 0.35 0.0225 612936 118344 19.31 36ms	35ms	59ms 130ms
		0.25 0.35 0.0240 612936	53498	8.73	36ms	35ms	38ms 109ms
		0.50 0.00 0.0150 460062 107732 23.42 103ms	83ms	136ms 322ms
	Fluid Data(100 th ) 0.50 0.00 0.0350 460062	61096 13.28 103ms	82ms	80ms 265ms
	2048x56x239	0.50 0.50 0.0400 460062 140120 30.46 104ms	83ms	93ms 280ms
		0.50 0.50 0.0500 460062	65320	14.2 104ms	82ms	63ms 249ms

TABLE 2

 2 Memory Cost

	Dataset	Parameters	SEL11 [10]	Our Method
	Discription	Iso αβMesh Pyramid(Texture) Mesh Pyramid(Buffer) Feature
	Head	0.25 0.00 0.0035 72.4MB	18 MB	18.0MB	63.3MB	224MB
	256x256x225x8b 0.25 0.00 0.0130 72.4MB	18 MB	6MB	63.3MB	224MB
	14.0MB	0.38 0.00 0.0040 85.3MB	18 MB	18.6MB	63.3MB	224MB
	-	0.38 0.00 0.0120 85.3MB	18 MB	6.6MB	63.3MB	224MB
	Hurricane	0.25 0.00 0.0015 63.1MB	1188 MB	13.6MB	49.3MB	212MB
	661x401x42x32b 0.25 0.00 0.0055 63.1MB	1188 MB	4.6MB	49.3MB	212MB
	42.4MB	0.25 0.35 0.0225 63.1MB	1188 MB	12.2MB	49.3MB	212MB
	18 th	0.25 0.35 0.0240 63.1MB	1188 MB	5.5MB	49.3MB	212MB
	Fluid Data	0.50 0.00 0.0150 47.4MB	9509MB	11.1MB	121.4MB	520MB
	2048x56x239x32b 0.50 0.00 0.0350 47.4MB	9509MB	6.3MB	121.4MB	520MB
	104MB	0.50 0.50 0.0400 47.4MB	9509MB	14.4MB	121.4MB	520MB
	100 th	0.50 0.50 0.0500 47.4MB	9509MB	6.7MB	121.4MB	520MB

This article has been accepted for publication in Computing in Science and Engineering but has not yet been fully edited.Some content may change prior to final publication.

ACKNOWLEDGMENTS

We would like to thank Hongsen Liao, Min Wu, Zhu Zhu, Kan Wu for their enthusiastic help and helpful suggestions. This research is supported by National Science Foundation of China(61272225, 51261120376), Chinese 863 Program(2012AA041606, 2012AA040902), Chinese 973 Program(2010CB328001) and the Important National Science & Technology Specific Projects (2011ZX02403).