Yi Peng
email: 15pengyi@gmail.com

Li Chen

Quan Li

Zhu Zhu

Jie Dong

Jun-Hai Yong

A Multi-Perspective Visual Retrieval Toolkit for Multi-dimensional Data

Keywords: perspective visual retrieval toolkit CR Categories: I.3.3 [Computer Graphics]: Three-Dimensional Graphics and Realism-Display Algorithms information retrieval, multi-perspective, visual navigation mode, focus+context

Rapid retrieval of specified items in a large amount of multidimensional data has become an urgent and challenging topic for improving the effectiveness of data searching and analysis approaches. In order to help users quickly re-find the desired items, we present a novel multi-perspective visual retrieval method which is totally different from traditional visualizations for searching results, such as list-based techniques. In our prototype, customized dimensions are generated adaptively which enables users to fast approach to the desired targets from multiple perspectives. We also design several display modes to layout the items by introducing degree-of-interest(DOI) ranking, which is implicitly defined by the

INTRODUCTION

Data Navigation for a large collection of information items is a difficult problem. It the information is well-defined, the process becomes a simple retrieval task. However, we may know vaguely about what we are searching for, as is the usual case. If the rapid retrieval of specified items involves complex logical thinking and memory recall, improving the effectiveness of data searching and analysis becomes urgent and challenging.

Existing traditional retrieval methods, such as list-based techniques which aim to rearrange items based on only one or two dimensions are not adequate for complex retrieval task. Moreover, the process of exploration of desired items is vital to analysis, which greatly helps users to understand the data in different applications.

Several existing techniques are proposed to enhance information retrieval, such as list-based item view, thumbnails of items to provide an overview, or re-sorting based on a certain attribute, etc. Meanwhile, researchers have found other ways to make full use of visualization to aid data retrieval, such as tree-based view, contentdependent bookmark [START_REF] Nakajima | Context-dependent web bookmarks and their usage as queries[END_REF], timeline-based toolkit [START_REF] Hoeber | Browseline: 2d timeline visualization of web browsing histories[END_REF], semantic map [ Šimko et al. 2010], tree map [Clarkson et al. 2009] and so on. However, they often focus on the hierarchical structure of data while ignoring other characteristics, such as multi-dimensional, entry date, etc., which may lead to ineffectiveness and impractical manipulation.

In order to solve these problems, we present a novel Multi-Perspective Visual Retrieval method (MPVR) which is different from traditional methods and enables users to re-find the desired items in high efficiency. Compared with other techniques, our method has four primary advantages.

• Rapid Retrieval

Customized dimensions are adaptively generated in a multiperspective view, making re-finding easy and efficient.

• Progressive Searching

Collaborative Filtering is provided to help user fast approach to the target progressively and self-feedback technology enables users to see details in fine scale.

• Visual Navigation Modes

In view of different favorites, we provide four visual navigation modes with adjustable parameters to analyze the data and understand their underlying relationships.

• Coordinated Highlighting

For better analysis, we integrate several visualization techniques, such as Focus + Context, Visual Cues , hierarchical view and animations, with our display modes to recommend results according to user interests.

Primitive evaluations in three scenarios: text message retrospection, desktop search and web browsing history retrieval are carried out, showing that our prototype is flexible, extensible, and efficient for re-finding and data analysis.

The remainder of this paper is organized as follows. We discuss the background and motivation behind MPVR before illustrating the functional capabilities of MPVR in Section 2. In Section 3 we make a comprehensive description of our work. The results of evaluation in several scenarios and different applications are given in Section 4. Section 5 summarizes the characteristics of our method and present the planned future work.

RELATED WORK

Research on visualization of information items has been conducted for several years, and related techniques that influence our design can be divided into several groups as follows.

Layout Approaches

Researchers have designed a variety of intelligent techniques to improve the effectiveness of information retrieval tasks. Crnovrsanin et al. propose a network visualization and exploration system which recommends relevant nodes to users based on different metrics. [Crnovrsanin et al. 2011] We extend this idea by using different metaphor and automatically the satisfied results will appear within an elaborate designed layout. Yi et al. present a novel multivariate visualization techniques that are easier to understand and more accessible by using a magnet metaphor and various interactive techniques, which enables the public to deal with multivariate data effectively [START_REF] Yi | Dust & magnet: multivariate information visualization using a magnet metaphor[END_REF].

In regard of retrieving previously revisited information, several techniques are proposed. Simko, Tvarozek and Bielikova present a novel approach to exploit web search and navigation history by collecting streams of user actions, identifying user goals and constructing visual trees to represent visiting history [ Šimko et al. 2010]. The result is based on semantic trees. However, this work requires a lot of efforts on semantic interpretation, which is not the usual case. Some works are based on tree layouts, such as Gandhi's work that implements a Domain Tree Browser to enable users to find their way back to the previously visited pages by building a tree structured visual navigation using domain names and thumbnails. [START_REF] Gandhi | Domain name based visualization of web histories in a zoomable user interface[END_REF]] Also, Mase and Yamada propose a structurebased hierarchical clustering method to show the topics and the relationships between them by extracting topic maps from histories. [START_REF] Mase | Extracting topic maps from web histories by clustering with web structure and contents[END_REF] To break the limit of two-dimensional screen, Frecon and Smith gives a virtual reality (VR) based application to be used alongside traditional web browsers, which provides users with a flexibly tailorable real-time visualization of their history. [START_REF] Vartiainen | Graphical history list with multi-window support on a mobile web browser[END_REF]

Focus + Context

In PaperVis Chou and Yang modify the existing Radial Space Filling (RSF) and Bullseye View techniques to arrange involved papers that better depict the relationships [Chou and Yang 2011]. Leung and Apperley give a survey on distortion-oriented visualization techniques including Polyfocal Display, Bifocal Display, Fisheye View and Perspective Wall [START_REF] Leung | A review and taxonomy of distortion-oriented presentation techniques[END_REF]. We use the same idea to highlight the focused targets while maintaining context information.

Treemap

Jarke J. Van and Wijk Huub Van de Wetering conceive the concept of Cushion Treemaps for hierarchical information visualization. They use a simple shading model to render the subdivision ridges providing insight in the hierarchical structures [START_REF] Van Wijk | Cushion treemaps: visualization of hierarchical information[END_REF]. ClarkSon, Desai and James apply Result Maps based on treemap techniques on online search interfaces to enhance digital library search engine result pages, which leads to ancillary benefits of data understanding [Van Wijk and Van de Wetering 1999]. However, the hierarchy property doesn't necessary lies in data itself.

Multi-Resolution

Hoeber and Gorner employ a two-dimensional timeline based visual interface to visualize browsing histories [START_REF] Hoeber | Browseline: 2d timeline visualization of web browsing histories[END_REF], which allows users to identify temporal patterns visually and re-find previously viewed pages effectively. However, their method is mainly based on time dimension and is not efficient enough to deal with large amount of data. Hightower et al. have implemented a browser companion which builds a graphical history-map of visited web pages dynamically using minimal screen space to enhance the usability of web browsers [START_REF] Hightower | links, objects, time and space-structure in hypermedia systems: links, objects, time and space-structure in hypermedia systems[END_REF]].

Other Related Visual Techniques

Several other related techniques also provide important references for our design.

Roberts surveys the techniques on coordinated and multiple views [START_REF] Roberts | State of the art: Coordinated multiple views in exploratory visualization[END_REF]], which we also integrate into our prototype. In order to help users analyze the histories from different aspects, we design visual navigation modes which are highly configurable and easy to operate.

Liang and Huang summarize the existing highlighting techniques that can be applied in information visualization and propose a new three-layer model of highlighting [START_REF] Liang | Highlighting in information visualization: A survey[END_REF]. In our method, we make a comprehensive use of highlighting techniques. The distortions of opacity, color, size and shape are used for view control, tooltips including lighting effects, directed edges and thumbnails are used as interaction control, and an extended fisheye lens is used to expose the most relevant data in the view center for visual recommendation.

THE PROPOSED VISUAL WEB HISTORY RETRIEVAL TOOLKIT

Our toolkit is an integrated system with three basic components (like abstract class): a data collection engine, a data viewer and an interactive visual interface. In different scenarios, the three components are application-specific derived from the basic one (like derived class).

Framework of MPVR

• Data Collection Engine

We design an abstract data collection engine by connecting to different data-sources. For example, in the application of web browsing history analysis, we develop a plugin using the kernel of Internet Explorer which provides users basic functionalities for surfing the Internet and connect it to our data collection engine (Figure2(a)). In another application of text message retrospection, this engine is connected to an export interface of smart phones.

• Data Viewer

After data collection, users can get an overview of the raw data and its statistics shown by four interdependent relational tables. The statistical data(Figure2(b) and 2(c)) is carried out automatically by a derived member function for data preprocessing in different scenarios. We can also sort them based on one attribute.

• Interactive Visual Interface

With the statistics, the general visual interface can directly show their relationships(Figure2(d)) with four features: rapid retrieval, progressive searching, visual navigation modes and coordinated highlighting techniques. The interface consists of three parts as seen in Figure2(d). The left area is for primary visualization. The top right region contains the configuration options for data dimensions and navigation modes, while the bottom right is used for adjustment of layout parameters.

Multi-Perspective View

The primary idea of our design is to make information item retrieval easy and efficient. So we design a multi-perspective view based on two important observed facts during the searching process.

Fact 1 Target items are multi-dimensional and can be represented by some major dimensions.

Fact 2 For users the general impression of an information item is always a description of several keywords.

With customized perspectives, retrieval and analysis of information items becomes easy due to the following three reasons.

Reason 1 With tags in different dimensions, data filtering can be much faster than before.

Reason 2 Searching process becomes progressive when different dimensions cooperate with each other as seen in Figure4.

Reason 3 By introducing related links, a hierarchical structure on each dimension is clearly presented.

Figure1 shows the customization and inner-relationship display features of multi-perspective view under different conditions. Fig- ure1(a), 1(b) and 1(c) show that the items are projected into five dimensions and a hierarchical structure is presented by interaction. Figure1(d), 1(e) and 1(f) show three other conditions restricted by three dimensions, two dimensions and only one dimension, respectively.

In web browsing history analysis application we elaborately choose five primary dimensions listed in Table 1 for history retrieval, which means each of the histories is projected into five dimensions with a value indicating its importance. We set up a local dictionary to filter out and extract those keywords related to semantic aspect, such as theme and the name of website, rather than manually labeling those history items.

Progressive Searching

Progressive searching provides users searching benefits by adding searching conditions implicitly through adding incremental con-

Collaborative Filtering

By selecting more dimensions and tags, more items will be filtered out since there exists more constraints, which means by using these filterings methods alternately, the number of result items may decrease significantly, therefore, users can approach to the targets rapidly.

(a) (b)

Figure 3: Collaborative Filtering for progressive searching.

Figure3 shows an example of Collaborative Filtering for progressive searching. For instance, the retrieval of the web page "yahoo mail" consists of only two steps. By selecting the "email" element around the "theme", the number of clusters greatly reduced as seen in Figure3(a). Further progress is made by specifying a "U-Z" tag to filter the results as shown in Figure3(b). Now the clusters are even fewer than before and the web pages "yahoo.mail.best.web" are already at the view center.

Self-feedback

Although Collaborative Filtering contributes a lot to the searching process, the number of filtered items can still be very large, which leads to more overlapping and chaotic layout. Here we use selffeedback technology to solve this problem.

This idea utilizes focus + context and takes hierarchical clustering into consideration. By interacting with overlapping items, we recalculate the position of each item in fine scale, keeping them separate from each other. Figure4 shows an searching example using self-feedback feature.

The original clusters are shown in Figure4(a). By right-clicking on the cluster "us.fires.wave.of", the overlapping items diffuse into different positions along a spiral path, as seen in Figure4(b). The location of each scattered item is calculated by its importance using fine scale statistics.

Visual Navigation Modes

Visual navigation modes are designed to meet the needs of different scenarios for searching. First, we will clarify some basic concepts.

Three display modes are designed to meet various needs of different stages in searching process. Before the introduction, some basic concepts need to be clarified.

• Item, Attribute and Statistic

An item indicates a record in the raw dataset. Each dimension is named an attribute. A statistic s is a set of items with the same attributes. For instance, a history item e may have at least two attributes: access frequency A f and residence (how long the user stays on this page).

These three concepts have corresponding visual elements in our prototype: Item Visual Element (IVE), Attribute Visual Element(AVE) and Statistic Visual Element(SVE).

• Statistic Importance

For each item, we use several attributes to define its importance. For example, giving a weight vector w, the importance of a history item I(e) is defined by its frequency A f and residence Ar. Then the statistic importance I(s) can be derived based on the item importance as shown in Formula 1. (1)

• Element Location

The location of a visual element is defined by (ρ(e),θ(e)) in polar coordinates system. The reason why we choose polar coordinates system is that the more central place an element locates, the more important it is. In this way, we can directly use the radius ρ(e) to measure its importance.

• AVE Layout Algorithm

All the AVEs are distributed evenly along the circumference of an ellipse, which means the included angles of two adjacent AVEs are the same. The location of an AVE is calculated by Formula 2.

θi =( i+0.5 na Figure5 shows the layout of AVEs which are evenly distributed on the boundary of an ellipse. Since the included angles of adjacent AVEs are the same, the areas of all the sectors are also the same(according to the geometrical characteristic of ellipse). The layout can be changed if the attribute importance is defined. Which means, the area of a sector will be weighted by the importance.

) • 2π ρi = a 2 cos 2 θi + b 2 sin 2 θi i =0, 1, ••• ,na -1 (2)

• SVE Layout Algorithm

We place a SVE around the AVE it belongs to. The relative polar coordinates to the corresponding AVE of a SVE is defined by Formula 3, where cs is a constant value, ns indicates the number of visible SVEs. If ns < 3 , the field angle defined by all the SVEs is set to π, otherwise it is equal to the inner angle of a regular polygon with ns vertices as shown in Figure6.

⎧ ⎨ ⎩ θi = ( i+0.5 ns )πn a < 3 ( i+0.5 ns )π • (1 -2 na ) na ≥ 3 ρi = cs i =0, 1, ••• ,ns-1 (3)

• IVE layout algorithm

The same principle applies when we layout the IVEs. In order to place IVEs with greater importance closer to the view center, we compute the radius ρ(e) carefully. The ρ(e) should take the statistic importance into consideration. For each SVE an IVE belongs to, we first compute a candidate location (ρc,θc), which is calculated by Formula 4:

⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ θc = θs + π 2 ρc = (xv -xs) • cos θs sin θs • cos θs sin θs T -(xv -xs) (4)
where xv denotes the location of the view center, xs denotes the location of a SVE, θs denotes the angle of the SVE. In specific, the candidate location is the projection of the view center on the line passing through the SVE in its direction.

wc =1-1 -a Is 2 1+a Is (5)

Formula 5 defines the importance of the SVE correspond to a IVE, which means SVE of less importance will have a weightvalue close to 1, while greater important SVE may have a weight close to 0. In our paper we choose a =0.05.

Cluster Mode

The Cluster Mode is designed to provide users an overview of the data, where IVEs belong to the same SVEs are clustered. In this mode, the distribution of information items is clearer. The position of an IVE is defined by Formula 6.

ρ(e)= {s|e∈s} wc(s) • ρc(s) θ(e)= {s|e∈s} wc(s) • αc(s) (6) Formula 6 indicates that the distance ρ(e) between IVE and view center is mainly determined by the least important SVE, which means if and only if all the corresponding SVEs are relatively important, the location of IVE will be central as shown in Figure7. Figure7 shows that the IVEs belonging to the same SVEs form a cluster.

Radial Mode

The radial mode diffuses the IVEs which belong to the same cluster into separate locations, which help users see detailed information.

ρ(e)= {s|e∈s} wc(s)

• ρc(s) θ(e)=( i+0.5 n i ) • 2π i =0, 1, ••• ,ni -1 (7)
This mode follows the computation of ρ(e) in the Cluster Mode, while calculating θ(e) to make the IVEs distributed evenly along the boundary of a sector as seen in Figure8.

We sort the IVEs based on their importance to determine their positions on the sector (Formula 7). Figure8 shows that the IVEs in a same cluster are diffused into different positions along the boundary of an sector.

Uniform Mode

The Uniform Mode is based on the Archimedean Spiral, which uses the same criteria with the Cluster Mode to compute ρ(e) and sort them in an ascending order. Finally, we use formula 8 to compute the index of an IVE in the ordered list.

t(e)= i+0.5 n i ,i =0, 1, ••• ,ni -1 ρ(e)=c2 • t e 2 + c3 θ(e)=c1 • t e 1 (8)
In this mode we set e1 = e2 =0 .5, c1 = r r i , c2 = rri and c3 = ri, where ri denotes the radius of IVEs and r is the size of the window.

Figure9 shows the Uniform Mode with the same data shown in Fig- ure7. In this mode the histories with the name "quan.li" are diffused along a spiral path, which means these histories lay on an ellipse boundary approximately.

Custom Mode

The custom mode is an extension of Uniform Mode. It allows advanced users to adjust layout parameters, such as c1, c2, c3, e1 and e2.

Coordinated Highlighting

For better user experience, we integrate several visualization techniques, such as Focus + Context, dynamic tooltips to present detailed information while retaining the context, therefore, users can concentrate on the specific items he/she is interested in.

For all the IVEs, the colors are defined by a user specified dimension. And for each visual element, its size and opacity are determined by its importance and when focused both of them will increase to show detailed information.

We also design dynamic tooltips to give a thumbnail and other important details. Moreover, when a user moves the mouse over a visual element, several linked light beams will be generated automatically to show the relationship between the focused one and other related elements. A color bar is also provided at the top-left corner, helping users understand the color mapping schema. Figure10 shows that we use a tooltip to show its detailed information for each visual element. Especially for information item, a thumbnail is presented to help users recapture memories on it. Meanwhile, a color bar is provided to show the correspondence of items and their colors.

EVALUATION METHODOLOGY

In this section we will demonstrate the effectiveness of MPVR through some potential scenarios in three familiar applications: we-b browsing history analysis, mobile text message retrospect, and desktop searching. Our case studies are based on the following three datasets: Web browsing histories can be used to help users retrieve important web pages which are used frequently. It also plays a significant role in analyzing user behaviors. It has been studied that over 80% of web pages have been previously visited [COCKBURN and MCKENZIE 2001]. Several techniques are integrated into web browsers to enhance the history utilization, such as list-based history view, thumbnails of most visited sites, URL auto-complete based on browsing histories, management of favorite links and so on [Ayers and Stasko 1995[START_REF] Hightower | links, objects, time and space-structure in hypermedia systems: links, objects, time and space-structure in hypermedia systems[END_REF]] [Cockburn et al. 1999[START_REF] Gandhi | Domain name based visualization of web histories in a zoomable user interface[END_REF][START_REF] Mase | Extracting topic maps from web histories by clustering with web structure and contents[END_REF]. However, they often focus on the hierarchical structure of data while ignoring their multi-dimensional properties. Therefore, in this scenario, we concern about how to effectively retrieve those web pages from a large amount of data by introducing multi-dimensional information. We show two cases that can be easily encountered in browsing experience to evaluate our method.

• Case Study 1: Finding webpages with high access frequency, these webpages are visited a certain time ago.

When serving the internet, users tend to re-visit a web page on a specific topic, such as shopping, online learning and so no. In this case, we carry out an experiment on dataset 1(Fig- ure13(a)). Suppose the goal is to re-find the shopping webpages of amazon, which are visited three month ago.

Since the target is clear defined, we first click the "shopping" element around "theme". Unfortunately, there are still no IVEs corresponding to "shopping" (as seen in Figure13(b)) due to the current websites having no relations with shopping.

We then select the "A-F" tag and click the "amazon" element, and we can see that more relevant items appear (see and 13(d)). We then choose the tag "¡3m" and click the "3 months" element to see how many webpages have been visited three month ago(see Figure13(e)).

Although the majority of the relevant items are presented, they overlap heavily. So we choose the Radial Mode and click the "¡50" element to see the most frequently visited pages(Figure13(f)). Then, we focus on the pages with the longest time of residence by choosing the "¡600" tag and clicking the "¡600" element. After that, there is only one page which meets our conditions (Figure11g). We can move the mouse over it to see the details(Figure11h). So we get the final page which is about "guitar", and maybe the user is a guitar fun. Actually what we guess is right.

Scenario 2: Mobile text message retrospect

SMS (Short Message Service) has already been a popular communication technology for mobile users. However, the visualization technology for text message retrieval and retrospect is not highly developed. There are many works that explore the social and temporal interactions in instant messaging which may involve different people, such as CrystalChat. This tool visualizes personal chat history and reveals the communication patterns [START_REF] Tat | Crystalchat: Visualizing personal chat history[END_REF]. In this scenario, we focus on how to aid users to retrospect their messages from the visualization aspect.

• Case Study 2: Finding important text messages contained key information from a certain person.

The scenarios is that a user collects recent mobile text messages and hope to check some key information, however he can't remember clearly. Suppose we know the name and location of the contact(a person). In general, the fastest way to find these messages is to find out all candidate messages from this person and check them one by one. Here, we use our toolkit to have a try. • Case Study 3: Searching for a document modified long time ago with specific extension.

Usually users takes a long time to do documents management, especially for the searching and maintenance of important documents such as resumes, schedules etc. In this case, we use dataset 3(Figure12(a)) to carry out an experiment on important document searching. Assume that the target document with "docx" file extension contains characters of "p"-"s" and is modified a long time ago but accessed recently.

We notice the particular file extension. So we choose the "office" tag and click the "document" element(as seen in Fig- ure12(b)). we then specify its creation time, write time and access time by clicking the "1 year" element around "creation time" and "write time", "1 month" element around "access time", and choose the "O-T" tag. After that many old documents appear as seen in Figure12(c). Since there are too many documents, we choose the Uniform Mode to diffuse them as seen in Figure12(d). Because the desired document is created a long time ago, we can make a further assumption that it may be placed in a deep directory. So we click the "4" element around "file path". Fortunately, the overlapping is eliminated(see Figure12(e)). Then we can check these 4 documents with the extension of "docx" one by one to see whether they are what we want(Figure12(f)).

Evaluation Results

The evaluation shows that our toolkit based on multi-perspective visual retrieval outweighs traditional retrieval methods. Meanwhile, we apply our method to three domains, web browsing history analysis, text message retrospect and desktop searching. And the three case studies turn out that our toolkit performs quite well in high dimensional data retrieval and the participants involved in the tasks think that our approach is attractive, efficient and quite interactive, while some of the participants maintain that they still need training before carrying out retrieval tasks.

CONCLUSION

In this paper, we present a novel approach named Multi-Perspective Visual Retrieval(MPVR) to help users quickly re-find the desired items of their interests. By providing four key features: Multi-Perspective View, Progressive searching, Visual Navigation Modes and Coordinated Highlighting, users can fast approach to the targets. Although the preliminary investigations give us positive feedback, there still remain a number of unsolved problems such as overlapping in Uniform Mode, unsatisfied layout with much more dimensions, less functionalities of coordinated views and so on.

Other future work includes improving user interface, applying our method to other domains and analyzing user's behaviors.

Figure 2 :

 2 Figure 2: The framework of MPVR.

Figure 4 :

 4 Figure 4: Self-feedback for progressive searching.

  A f (e)+wr •Ar(e),w f +wr =1

Figure 5 :

 5 Figure 5: The layout of customized dimensions.

Figure 6 :

 6 Figure 6: Different layouts of SVE.

Figure6

  Figure6 shows different layouts of SVEs. If the number of SVEs is more than a default value(we choose 5 in our toolkit), users need to specify a sliding window, which means only 5 SVEs are visible. The visible SVEs are evenly distributed on a circle around the AVE they belong to.

Figure 7 :

 7 Figure 7: The Cluster Mode for IVE layout.

Figure 8 :

 8 Figure 8: The Radial Mode for IVE layout.

Figure 9 :

 9 Figure 9: The Uniform Mode for IVE layout.

Figure 10 :

 10 Figure 10: The coordinated highlighting techniques.

Figure 11 :

 11 Figure 11: Finding important text messages contained key information from a certain person.

Figure 13 :

 13 Figure 13: Finding specific webpages of high access frequency visited a certain time ago.

Table 1 :

 1 Dimensions of Browsing History

	Dimension Symbol	Description
	Website	w	URL of a web page
	Theme	m	Topic to which a web page belongs
	Time	t	Visit time of a web page
	Hits	n	Number of visits on this page
	Residence	r	Time spent on a web page
	Thumbnail		Detail information of a web page

straints, thus reducing the number of result items and helps users quickly find what they want. This mechanism is composed of two steps below.

Table 2 :

 2 Three datasets for method evaluation

	ID	Name	Items Dimensions
	1	Web browsing history 1218	12
	2	Mobile text message	1886	14
	3	Local Computer Files 1857	9
	4.1 Scenario 1: Browsing history retrieval

Acknowledgements

This research is supported by National Science Foundation of China(61272225, 5126112037691315302), Chinese 863 Program(2012AA041606, 2012AA040902), Chinese 973 Program(2010CB328001) and the Important National Science & Technology Specific Projects (2011ZX02403)

AYERS, E. Z., AND STASKO, J. T. 1995. Using graphic history in browsing the world wide web.


K.-L Shi

J.-H Yong

L T Ang

J.-G Sun

J.-C Paul

Polar NURBS Surface with Curvature Continuity

Keywords: (according to ACM CCS): I.3.5 [Computer Graphics]: Computational Geometry and Object Modeling-Curve, surface, solid, and object representations; Splines '2, FJI

Polar NURBS surface is a kind of periodic NURBS surface, one boundary of which shrinks to a degenerate polar point. The specific topology of its control-point mesh offers the ability to represent a cap-like surface, which is common in geometric modeling. However, there is a critical and challenging problem that hinders its application: curvature continuity at the extraordinary singular pole. We first propose a sufficient and necessary condition of curvature continuity at the pole. Then, we present constructive methods for the two key problems respectively: how to construct a polar NURBS surface with curvature continuity and how to reform an ordinary polar NURBS surface to curvature continuous. The algorithms only depend on the symbolic representation and operations of NURBS, and they introduce no restrictions on the degree or the knot vectors. Examples and comparisons demonstrate the applications of the curvature-continuous polar NURBS surface in hole-filling and free-shape modeling.

Introduction

Most CAD systems and many geometric modeling systems in computer graphics represent free-form surfaces by NURBS (non-uniform rational B-splines) and its special case NUPBS (non-uniform polynomial B-splines). They have become de facto industrial standards for model representation, processing and data exchange [START_REF] Farin | Curves and surfaces for computer-aided geometric design: a pratical guide[END_REF]PT97a]. However, the rectangular parametric domain also has been pointed out not to be well-suited for capping disc-shaped circular regions [START_REF] Peters | C 2 splines covering polar configurations[END_REF][START_REF] Shi K.-L | G n filling orbicular N-sided holes using periodic B-spline surfaces[END_REF], which are common in industrial geometric modeling, for example, the head of an airplane, the end of an aerofoil or a screw propeller (see Figure 1). A natural way to adapt the quadrilateral untrimmed NURBS surface to a cap-like shape is as follows. First, roll up the surface and glue two opposite boundary edges together to form a periodic NURBS surface. Then, shrink a loop edge to a single polar point as shown in Figure 2. The † The research was partly supported by Chinese 973 Program (2010CB328001) and Chinese 863 Program (2012AA040902). The second author was supported by the NSFC (61035002, 61272235). The fifth author was supported by the NSFC (61063029, 61173077). However, at the pole of the polar NURBS surface, there is still a challenging issue about curvature continuity (G 2 continuity), which is essential and sometimes the minimum standard in geometric modeling such as class-A car-body design. At an interior point of a polar NURBS surface, the infimum order of parametric continuity, which is determined by the multiplicity of knots, has nothing different from ordinary periodic NURBS surface. However at the pole, which is actually a degenerate boundary edge, continuity also depends on the surrounding control points. The polar curvature of a non-well-formed polar NURBS surface is usually multivalued and that violates curvature continuity at the pole (see Figure 7). This problem prevents the polar NURBS surface from achieving higher-order continuity in practice. In order to solve the problem, we analyze the local differential feature of the singular pole and then propose methods to construct polar NURBS surfaces with curvature continuity. To be specific, the contributions of this paper are as follows. For the polar NURBS surface, we -propose a sufficient and necessary condition of the polar curvature continuity; -present a method to construct the control points that surround the pole to satisfy curvature continuity; -present a method to reform a non-curvature-continuous polar NURBS surface to curvature-continuous, preserving the overall shape of the surface.

Not like other subdivision or spline-based methods, our work focuses on the relation between curvature continuity and radial cross-boundary derivatives in the NURBS form.

All the procedures are based on the symbolic representation and operations of NURBS [START_REF] Piegl | Symbolic operators for NURBS[END_REF], which introduce no restrictions on the input, for example, the degree and the knot vectors can be arbitrary.

Overview. After reviewing the related work in Section 2 and preliminaries in Section 3, we start from a sufficient and necessary condition of curvature continuity at the pole in Section 4. Based on the condition, Sections 5 and 6 then give a construction method and a reformation algorithm, respectively. The applications of n-sided hole filling and free-shape design in Section 7 illustrate the practicability of the curvature-continuous polar NURBS surfaces in CG and CAD. Section 8 concludes the paper.

Related work

Some methods of shrinking one edge of a NURBS surface to a single point have been proposed since the introduction of NURBS. For example, in Ref.

[PT97a] §8.7, a method was presented to construct a three-sided spherical surface with one degenerate boundary. Ref. [START_REF] Shi K.-L | G n blending multiple surfaces in polar coordinates[END_REF] introduced a polar parametric domain in blending multiple surfaces around an n-sided region. In Ref. [START_REF] Shi K.-L | G n filling orbicular N-sided holes using periodic B-spline surfaces[END_REF], the periodic B-spline surface was used to fill n-sided holes. They assumed that the polar curvature is zero to avoid the problem of curvature continuity at the pole. Besides the NURBS form, studies on G 2 continuity at the extraordinary vertices have also yielded substantial results by subdivision [KMP06,KP07,Lev06,MP09]. However as discussed in Ref. [START_REF] Peters | C 2 splines covering polar configurations[END_REF], those subdivision methods do not readily fit into existing CAD frameworks, and the construction process, which may consist of an infinite sequence of surfaces, is unfriendly in further geometric processing such as intersection and trimming. There are also many methods to generate curvature-continuous splines from coarse quad meshes, for example, Hartmut and Prautzsch's free-form splines [Pra97], Karčiauskas and Peters' guided spline surfaces [KP09] and Reif's TURBS [START_REF] Reif | TURBS -topologically unrestricted rational Bsplines[END_REF]. Recently, Myles and Peters [START_REF] Peters | C 2 splines covering polar configurations[END_REF] gave a spline solution to the problem of curvature continuity for polar configurations in quad-dominant meshes. Their C 2 polar spline surface consists of one tensor-product spline of degree (3, 6), which is bi-3 C 2 compatible. It is more CADfriendly than polar subdivision. However, in real industrial CAD/CAM applications, we normally cannot limit the knot vectors, the degrees and the weights to this single pattern, and moreover, especially for mechanical modeling, the inputs are usually trimmed NURBS surface patches instead of two-manifold quad-dominant meshes. The degrees of the input patches may be different, and therefore the control points of two adjacent surfaces may not coincide even though they share the same edge geometrically. From a different aspect, this paper presents a more adaptive solution especially for CG/CAD/CAM modeling. The methods are derived from the symbolic representation of the cross-boundary derivative curves, which will not be limited by the inconsistent control meshes, degrees or knot vectors. In implementation, the symbolic operators can be applicable to both NURBS and NUPBS without changing the formulae.

Definitions

The periodic NURBS surface discussed in this paper is a NURBS surface that has a periodic knot sequence in only one parametric coordinate. It can be formulated as follows. Parameter (u, v) is defined in a horizontal band in Euclidean space (see Figure 3). P w [i, j] are the control points. In this paper, the superscript w always denotes four-dimensional vectors in homogeneous space, and the corresponding notations without w represent their three-dimensional Euclidean projections. A pair of square brackets denotes the index of the control points of NURBS. The B-spline basis functions Ñi,p (u) and N j,q (v) have the same recursive form with the standard NURBS in Ref.

S w (u, v)= ∑ i, j Ñi,p (u)N j,q (v)P w [i, j] (1)
[PT97a]. The only difference is that in u-direction, or the so-called periodic direction, the knot sequence Ũ is infinite and periodic. We denote a complete nondecreasing subsequence of it in one period by U = {u 0 , u 1 ,...,um}. In the other direction, the knot sequence V = {v 0 , v 1 ,...,v l } is a standard knot vector. Without loss of generality, we assume u 0 = 0, v 0 = v 1 = ... = vq = 0, v l-q = ... = v l-1 = v l = 1 and the period of Ũ is 1.

The topology of a non-degenerate periodic NURBS surface, as well as its control quads, is usually isomorphic with a cylindrical surface. It has two opposite periodic boundary edges as shown in Figure 2(a). We can construct a capshaped polar NURBS surface by shrinking a boundary edge of it to a single degenerate point in the center. The control mesh of it is in the polar configuration [START_REF] Peters | C 2 splines covering polar configurations[END_REF], which has triangular fans surrounding the the central pole. Formally, we denote the pole by P 0 = P[i, 0] for all valid i. As shown in Figure 2(b), a polar NURBS surface are isomorphic with a cap or a disk, and therefore, it can be used to represent some cap-shaped features, which are usually difficult to be represented by untrimmed standard NURBS surfaces with four boundary edges. Figure 1 shows an example.

Further analysis indicates that a p-degree B-spline basis function is C p-k continuous at a knot with multiplicity k, which denotes the number of overlapping knots of the same value. However the polar NURBS surface has an exception at the degenerate pole: the polar curvature continuity also depends on the two rings of the surrounding control points besides the multiplicities of the knots. In order to focus on this

1 'YX < ; 'Y 'X 7 &W 3 6X 3>@ 3>@ 3>@ 3>@ 3>@ 3>@ 3>@ 3>@ 3>@ 3>@ 3>@ 3>@ IJ Figure 4: Polar continuity. P[i, j] (i = 0, 1,...,5;j= 1, 2)
are the control points. The surrounding black arrows show the radial cross-boundary derivatives. X and Y span the tangent plane τ of the pole P 0 . T is an arbitrary direction on the local plane spanned by Du and Dv at parametric point (u 0 , 0), and the radial curve C(t) is in that direction.

special case and study the relation between polar curvature continuity and the zeroth-, the first-and the second-order derivatives determined by these control points, we assume that the multiplicities of the knots satisfy C 2 condition of the corresponding basis functions.

Condition of curvature continuity at the pole

First, we introduce a simplified notation of the derivatives.

Du ...u r v ...v s (u)= ∂ r+s S (l, v) ∂l r ∂v s l=u v=0
We omit the parameter "(u)" in the following derivation. A polar NURBS surface S(u, v) has a degenerate isoparametric curve S(u, 0)=P 0 for all valid u. At point S(u, 0), normal (G 1 ) continuity requires that all the derivatives T = αDu + βDv (∀α, β ∈ R) should be on the polar tangent plane τ (see Figure 4). So there exists a unit normal vector N that satisfies T • N = 0. Since P[ * , 0] (" * " denotes all valid indices) shrink to P 0 ,w eh a v eDu ≡ Duu ≡ 0 according to the derivative formula on the boundary

[PT97a]. It follows that N • Dv = 0.
The NURBS derivative formula also presents that the radial derivative Dv is a linear combination of the pole and the first ring of the control points P[ * , 1]. So we have the following sufficient and necessary condition of G 1 polar continuity.

Condition 1 The polar NURBS surface is G 1 continuous at the pole if and only if the control points P[i, 0] and P[i, 1] are coplanar for all valid i.

Remark. Strictly, this paper only considers ordinary manifolds that have neither overlap, fold or normal flip. We can simply guarantee this by assuming that

((P[i, 1] -P 0 ) × (P[(i + 1) mod n, 1] -P 0 )) • N > 0,
where n is the number of control points in the periodic direction. For curvature continuity, we first assume P[ * , 1] = P 0 to avoid the degenerate case: Dv = 0. The G 1 continuity condition guarantees the existence of two vectors X and Y satisfying Dv(u)=α(u)X + β(u)Y, where α(u) and β(u) are scalar functions (see Figure 4). As a manifold without overlaps, any direction on the polar tangent plane has a unique corresponding parameter u satisfying that equation. So on the parametric point (u, 0) with the corresponding spatial direction Dv(u), the normal curvature K of the curve C(t)= S(γ •t + u, δ •t) can be formulated as follows according to the definition and the expanding rule of the derivative of composite functions.

K(u, γ, δ)= ∂ 2 C(t)/∂t 2 • N ∂C(t)/∂t • ∂C(t)/∂t t=0 = γ 2 Duu • N + 2γδDuv • N + δ 2 Dvv • N γ 2 Du • Du + 2γδDu • Dv + δ 2 Dv • Dv .
Note that we already have

Du = Duu = 0. It follows that Duu • N = 0, Duv • N = 0, Du • Du = 0 and Du • Dv = 0. So we have K (u, γ, δ)=( Dvv • N) (Dv • Dv).
The directional normal curvature K is a constant at S(u, 0). No matter how to choose γ and δ, the polar curvature K only depends on the parameter u, as well as the corresponding derivatives Dv(u) and Dvv(u). That is,

K (u)=(Dvv(u) • N) (Dv(u) • Dv(u)). (2)
According to the definition of geometric continuity, a surface is G k continuous at a specified point if and only if there exists a reparameterization that can reform the surface to C k continuous locally. So the sufficient and necessary condition of curvature continuity at the pole is that there exists a C 2 continuous parametric surface Ŝ(x, y) that osculates the pole in order two. We suppose that the local frame of the osculation surface is (P 0 , X, Y). That is, ∂ Ŝ/∂x = X and ∂ Ŝ/∂y = Y at P 0 . Then, in the spatial direction Dv(u)=α(u)X + β(u)Y, the directional normal curvature is as follows.

K (u)= α 2 (u) L • N + 2α (u) β (u) M • N + β 2 (u) R • N α 2 (u) X • X + 2α (u) β (u) X • Y + β 2 (u) Y • Y , (3)
where L = ∂ 2 Ŝ/∂x 2 , M = ∂ 2 Ŝ/∂x∂y and R = ∂ 2 Ŝ/∂y 2 at P 0 . The osculation surface is C 2 continuous if and only if the three second-order partial derivatives L, M and R are unique, respectively. Note that the denominator part equals Dv(u) • Dv(u), and therefore, combining (2) and the nondegenerate assumption of Dv(u),wehave

α 2 (u) L • N + 2α (u) β (u) M • N + β 2 (u) R • N = Dvv (u) • N.
We introduce two scalar functions ξ(u) and φ(u) to remove "•N" in each items. That is,

Dvv (u)=α 2 (u) L + 2α (u) β (u) M + β 2 (u) R + ξ (u) X + φ (u) Y, (4)
or the equivalent form with three scalars L, M and R

Dvv (u)=( α 2 (u) L + 2α (u) β (u) M + β 2 (u) R) • N + ξ (u) X + φ (u) Y. (5)
So we have the following sufficient and necessary condition of curvature continuity at the pole.

Condition 2 The polar NURBS surface S(u, v) is curvature continuous at the pole if and only if 1. the surface satisfies the G 1 polar continuity condition; 2. with independent vectors X and Y (typically the two principle directions) satisfying Dv(u)=α(u)X + β(u)Y, there exist three vectors satisfying (4), or three scalars satisfying (5).

Construction of a polar NURBS surface with curvature continuity

A (p, q)-degree polar NURBS surface is G 2 continuous if it is a) C 2 continuous and b) G 2 continuous at the pole. For the first condition, we only need to ensure that the multiplicity of each knot of U is not greater than p -2 and of V , not greater than q -2(p, q ≥ 3). Besides, as what geometric modeling systems usually do, we should also ensure that the control-point mesh has no superposition or flip. Then, the key problem is to construct the control points to satisfy the condition of polar curvature continuity.

According to the derivative properties of NURBS on the boundary, once P w [ * , 0] are fixed, the radial derivative curve Dv(u) and the first ring of control points P w [ * , 1] are interconvertible exactly (see Ref.

[PT97a] §4.5, §5.3 and §5.5). For the same reason, Dvv(u) and P w [ * , 2] are interconvertible once P w [ * , 0] and P w [ * , 1] are fixed. So instead of the control points, we discuss the first-and the second-order radial derivatives Dv(u) and Dvv(u) as symbols, since they are equivalent to the corresponding control points. This symbolic method helps to unify the formulae of NURBS and NUPBS.

In practice, there are two cases when the pole is prescribed: the first ring of the control points (or Dv(u)) are specified, or not. For the case that Dv(u) is not given, we first discuss how to construct Dv(u) to satisfy G 1 continuity. We freely specify N, the unit normal vector of the pole. Then, we choose two linearly independent vectors X and Y to satisfy X • N = 0 and Y • N = 0, which span the tangent plane τ at the pole (see Figure 4). Condition 1 implies that

Dv(u)=α(u)X + β(u)Y, (6)
where α(u) and β(u) are scalar periodic NURBS functions, and the constant vectors X and Y also can be treated as onedegree one-control-point periodic NURBS curves. With the symbolic operators of addition and scalar multiplication of NURBS [START_REF] Piegl | Symbolic operators for NURBS[END_REF], we directly compute the NURBS form of Dv(u). Its degree is the maximum of the degrees of α(u) and β(u). The introduction of the symbolic operators here makes the formulae simple and clear, and it prevents the readers from being lost in the detailed implementation. α(u) and β(u) can be arbitrary. An ellipse-shaped radial derivative curve like the one in Figure 4 can be formulated by Dv(u)=cos(u)X + sin(u)Y. Compared with (6), we choose α(u) ≈ cos(u) and β(u) ≈ sin(u). The quality of approximation does not interfere with continuity however only has slight consequences for the shape. So for simplicity, we usually use a uniform knot vector and the same weight of control points to construct α(u) and β(u). For example, the control points of α(u) and β(u) can be {0, 1, 0, -1} and {1, 0, -1, 0}, respectively. More strictly, the trigonometric functions can be exactly represented by NURBS after the following reparameterization: sin(u)=2t/(1 + t 2 ), cos(u)=( 1t 2 )/(1 + t 2 ) where t = tan(u/2) (see Ref.

[PT97a] §7.5). Of course, Dv(u) also can be designed manually satisfying Condition 1. Backwardly, α(u) and β(u) are retrieved from Dv(u) (or indirectly from the first ring of control points). Once Dv(u), X and Y are determined, we have

α[i]X + β[i]Y = Dv[i] for each control point with index i.It follows that [α[i], β[i], 0] t =[XYN] -1 Dv[i].
Both α(u) and β(u) share the same degree, weights and knots with Dv(u).

Then we discuss the construction of Dvv(u) assuming that Dv(u) (or P w [ * , 1]) has been determined. Condition 2 gives the construction formulae (4) and ( 5), where ξ(u) and φ(u) are arbitrary scalar periodic NURBS curves. They involve the symbolic operations of addition, inner product and scalar multiplicity of NURBS [START_REF] Piegl | Symbolic operators for NURBS[END_REF], and the final result of Dvv(u) still preserves the NURBS form. Its degree is the maximum of the twice of the degree of Dv(u), the degree of ξ(u) and the degree of φ(u). Note if {X, Y, N} form an orthonormal basis, (3) can be reduced to the following equation by assuming that M = 0.

K(u)= α 2 (u)L + β 2 (u)R α 2 (u)+β 2 (u) .
At a prescribed parametric point u, the curvature K(u) is a positive combination of L and R. In non-degenerate cases, the continuous curves α(u) and β(u) have different zero points. So the assumption of M ensures that L and R are the maximum and the minimum of the curvature function K(u), which are the so-called principle curvatures, respectively. Thus, we can use that formula to construct Dvv(u) by principle curvatures. For ξ(u) and φ(u), we usually first ignore them in the construction of Dvv(u), and then move the control points of Dvv(u) along the tangent plane spanned by X and Y. This method assumes that ∆(u)=ξ (u) X + φ (u) Y shares the same degree, weights and knot vector with Dvv(u). So the addition of ∆(u) to Dvv (u) can be reduced to pointwise additions of the control points.

Once P 0 , Dv(u) and Dvv(u) are determined, the three central rings of control points can be constructed according to the derivative formulae of NURBS surface on the boundary [PT97a]. In practice, the radial derivatives are usually polynomial (for example, n-sided hole filling), and in these cases, the conversion algorithm from radial derivative curves to control points can be as simple as follows.

,I WKH SULQFLSOH FXUYDWXUH VSHFLILHG XVH RUWKRQRUPDO EDVLV ; < 1 DQG OHW ) -Form a one-degree periodic NURBS curve Do(u) with only one control point P 0 . -Elevate the degrees of the curves Do(u), Dv(u) and Dvv(u) to the maximum degree of them, denoted by p. -Compute the union of the knot vectors of the three radial derivative curves, denoted by Ū, and then ensure they share the same knot vector Ū by knot refinement (so the numbers of control points are also the same). -For the kth ring of the control points of the surface, we use the same weight w k . For each column i, the three control points P[i, 0], P[i, 1] and P[i, 2], as well as the corresponding weights w 0 , w 1 and w 2 , form the starting part of a NURBS curve with degree q and knot vector V that come from the radial direction of the polar NURBS surface. Then, the three control points can be solved by Hermite interpolation using Do[i], Dv[i] and Dvv[i] as boundary conditions at v = 0.

,QSXW 3 3 Z >@ VSHFLILHG" &RQVWUXFW 3 Z >@ DFFRUGLQJ WR &RQGLWLRQ &RQVWUXFW 'YXZLWK ĮX DQG ȕX DFFRUGLQJ WR &RQGLWLRQ &DVH 5HWULHYH 'YX ĮX DQG ȕX IURP 3 Z > @ 1 &RQVWUXFW 'YYXXVLQJ )RUPXOD RU < &RQVWUXFW 3> @ IURP 3 'YX DQG 'YYX 2XWSXW 3 Z > @ DQG 3 Z > @ VHH >37D@ V\PEROLF RSHUDWRUV RI 185%6 VHH >37E@ VXUIDFH FRQVWUXFWLRQ ZLWK ERXQGDU\ FRQGLWLRQV VHH >37D@ &DVH
The entire construction process is summarized in Figure 5. Our method employs the Hermite interpolation algorithm of NURBS surface, which aims to compute the control points of the result surface column by column, with the control points of the specified cross-boundary derivatives. If Dv(u) and Dvv(u) are rational, the formulae in the last step will be more complicated in implementation. We do not discuss the details of it here since it is a fundamental operation of construction of NURBS surface [PT97a].

Reforming an ordinary polar NURBS surface to curvature continuous

The reformation algorithm aims to achieve curvature continuity at the pole by moving the control points locally. First we assume that the pole P 0 and the first ring of the control points are coplanar so that G 1 continuity is guaranteed according to Condition 1. If not, we only need to project the control points onto an estimated polar tangent plane. So the local frame {X, Y, N} at the pole can be constructed by choosing two independent vectors X and Y that are perpendicular to the polar unit normal vector N. Then, Dv(u) and Dvv(u), as well as α(u) and β(u), can be retrieved from the three central rings of the control points. We decompose Dvv(u) into two orthogonal parts: the scalar normal part Dvv(u) and the tangential part ∆(u). That is,

Dvv(u)=Dvv(u) • N + ∆(u), ∆(u) • N ≡ 0. ( 7
)
The NURBS curves Dvv(u) and ∆(u) share the same degree, weights and knot vector with Dvv(u), and pointwise,

Dvv[i]=Dvv[i] • N and ∆[i]=Dvv[i] -Dvv[i]
• N for all valid i. Compared with (5), the surface is curvature continuous if and only if there exist three scalar functions satisfying

Dvv(u) ≡ α 2 (u)L + 2α(u)β(u)M + β 2 (u)R. (8)
(8) is an alternative form of (5) in Condition 2. A randomly constructed polar NURBS surface usually fails to achieve curvature continuity because this equation can be rarely satisfied without careful construction of the second-order radial derivative curve Dvv(u).

The satisfaction of (8) depends on all valid u. Fortunately, the NURBS form of the scalar functions allows this problem to be converted to finite times of verification of the following equation about the control points.

Dvv[i]=(α 2 )[i] • L + 2(αβ)[i] • M +(β 2 )[i] • R, ∀i, (9)
where the NURBS functions Dvv, (α 2 ), (αβ) and (β 2 ) share the same degree, weights and knot vector after degree elevation, knot refinement and uniformization of denominators. Dv(u) of an ordinary surface has at least three independent control points (or, the surface will be ill conditioned). So the number of rows of the linear system (9) is no less than the number of the unknowns (L, M and R). Thus, the surface is curvature continuous at the pole if and only if the linear system has a unique solution (where all equations are satisfied).

Otherwise, the overdetermined system only has a leastsquare solution { L, M, R} by Moore-Penrose pseudoinverse [Moo20], and that gives an approximation to the discontinuous pole. Substitute them to (8) and then to (7), we can construct Dvv(u) to satisfy polar curvature continuity. The final control points can be computed as the previous section. A shortcoming of this method is the introduction of degree elevation and knot insertion, which may cause the increment of the number of control points. However this is necessary if we need to exactly preserve the control points of the first ring, as well as the pole. Otherwise, we should simplify α(u) and β(u) before the computation of (8). Furthermore, if we want to preserve the number of the control points of the original surface, as well as their weights and knot vectors, it is also practical to slightly move the control points of the second ring to fit Dvv(u) to the ideal Dvv(u) that satisfies polar curvature continuity. Both these fitting and simplification methods may introduce an error. The choice of strategy depends on the application and the balance between precision and simplicity. We do not expand the very details here since they are surely beyond the limit of this paper.

Examples and applications

We first present examples of the results by our method. The constructive method proposed in Section 5 is used to generate free-shape surfaces. Figure 6 demonstrates three cases, which have an ellipsoidal (two different principle curvatures with the same sign), a saddle (principle curvatures with different signs) and a spherical (same principle curvatures) pole, respectively. The curvature charts illustrate the smoothness and continuity of the shape. In order to ensure C 2 continuity in the periodic direction, the first-order radial derivative curves are cubic. So the degrees of the final surfaces are all (6, 3), which are acceptable in CG and CAD.

In Figure 7, we reform a non-G 2 polar NURBS surface to G 2 continuous. The algorithm uses the pole and the first ring of the control points to derive the zeroth-and the first-order radial derivative curves. The second-order radial derivative curve is then constructed according to (7) and (8). After ensuring the consistence of the degree and the knots of it with the periodic direction of the surface, the new positions of the control points of the second ring are then solved. In this case, the degree of the first-order radial derivative curve can be reduced to 3. Formula (8) doubles it to 6, so the final surface has degree 6 in the periodic direction. Our G 2 cap-like surface is also proved to be a powerful tool in filling holes, which aims to provide one or more surface patches to interpolate the hole's boundary and crossboundary derivatives with continuity constraints. For example, a recent work [START_REF] Shi K.-L | G n filling orbicular N-sided holes using periodic B-spline surfaces[END_REF] used a polar NURBS surface to cap a hole that has a loop-shaped boundary. This method avoids the trimmed domain and is more CADfriendly than the other methods based on blending or subdivision [GLZ90, NSY09, Rei98, SYSP10]. However it also suffers from the problem of polar curvature continuity. As a practical application, the following example demonstrates how our method solves this continuity problem. In Figure 8, we combine the methods in [START_REF] Shi K.-L | G n filling orbicular N-sided holes using periodic B-spline surfaces[END_REF] and in this paper to generate a G 2 polar NURBS surface that contacts the specified five boundary surfaces with curvature continuity. We achieves the curvature continuous pole with non-zero Gaussian curvature compared to the artificial flat pole by their method. Figure 9 shows a mechanical model in which we fill the hole of the rounded tip of each propeller blade with a single G 2 polar NURBS surface. The boundary of the hole, which is reconstructed from sectional views, is a quartic NURBS curve that has non-uniform knots to represent the sharp feature of the blade. ply as we can, we form the first-order radial derivative curve Dv(u) in periodic NURBS using the same degree and knot vector with the boundary Bo(u). After that, we estimate the control points of Dv(u) to achieve a better shape. For each control point of Dv(u), we project the corresponding control point of Bo(u) onto the polar tangent plane, and let it be the vector from O to the projection point. That is,

Dv[i]=(Bo[i] -O) -((Bo[i] -O) • N) • N, (10)
for all valid i. After choosing an arbitrary orthonormal basis {X, Y, N}, we can retrieve the scalar functions α(u) and β(u), which have the same degree, weights and knot vector with Dv(u). The three parameters L, M and R can be specified or estimated by finding a least-square fitting to the normal curvatures of each parabola through

Bo[i], B τ o [i] and O. That is, K[i]=2((Bo[i] -O) • N)/(Dv[i] • Dv[i]).
According to the formula of directional curvature (3), we simply have the following linear system of L, M and R.

K[i](α[i]+β[i]) 2 = α[i] 2 • L + 2α[i]β[i] • M + β[i] 2 • R,
for all valid i. Use Moore-Penrose pseudoinverse to get the least-square solution { L, M, R} and then construct Dvv(u) as the previous section. After degree elevation and knot refinement to ensure the consistence of the curves of boundary conditions D h (u) and B h (u) (h = o, v, vv), we can generate the final surface by the method in Ref. [START_REF] Shi K.-L | G n filling orbicular N-sided holes using periodic B-spline surfaces[END_REF].

The algorithm in Ref. [START_REF] Shi K.-L | G n filling orbicular N-sided holes using periodic B-spline surfaces[END_REF] assumes the zero polar curvature to guarantee polar curvature continuity. It leads to a flat region at the pole and a high-curvature area around the pole (see Figure 8(d)). The comparison of the distribution of Gaussian curvature illustrates that the non-flat pole improves the quality (smoothness) of the surface. Note the control meshes of the boundary surfaces may not coincide with the filling surface although they share the same boundary geometrically. This is common in NURBS surface modeling for two possible reasons: 1) different degrees of NURBS, and 2) different parameterization. The final surface is friendlier to CAD since periodic NURBS surfaces are interconvertible with standard NURBS exactly (the control points in all figures are of the corresponding standard NURBS).

Conclusion

We proposes a sufficient and necessary condition of curvature continuity at the pole of a polar NURBS surface. Based on that criterion, we present two constructive algorithms to generate a polar NURBS surface with curvature continuity and to reform a polar NURBS surface to be curvature continuous. These methods focus on more general cases compared to [MP09, MP11] since for our methods, there are no limitations on the degree, the knot vector, as well as the weights, of the input data. Formulation in this paper only involves the symbolic representation and operations of NURBS rather than control points or subdivision formulae. Besides, a polar NURBS surface can be directly and precisely converted to standard NURBS without elevating the degree or introducing a tolerance. The methods are compatible with existing CAD/CAM systems.

Figure 1 :Figure 2 :

 12 Figure 1: A transitional structure of a mechanical part designed by SolidWorks 2007 in (a). The 14 surface patches in (a) can be replaced by a single polar NURBS surface in (b).

cFigure 3 :

 3 Figure 3: The periodic knot sequence in u-direction of the domain (blue band) and the corresponding basis functions.

  c 2013 The Author(s) Computer Graphics Forum c 2013 The Eurographics Association and John Wiley & Sons Ltd

c

  2013 The Author(s) Computer Graphics Forum c 2013 The Eurographics Association and John Wiley & Sons Ltd

Figure 5 :

 5 Figure 5: The flow chart of the construction method.

c

  2013 The Author(s) Computer Graphics Forum c 2013 The Eurographics Association and John Wiley & Sons Ltd

Figure 6 :

 6 Figure 6: Curvature continuous polar NURBS surfaces. They have (a) ellipsoidal, (b) saddle and (c)(d) spherical poles, respectively. The red curves are the first order radial directive curves Dv(u). The Gaussian curvature charts in (a), (b) and (c) show their continuity.

Figure 7 :

 7 Figure 7: Reforming a non-G 2 polar NURBS surface to G 2 . (a) shows the control points of the original discontinuous surface, while (b) presents its Gaussian curvature chart. Note that around the pole, the Gaussian curvatures converge to different values in different directions. (c) shows the control points of the reformed surface. The blue dots are the original control points before moving, and the red curve represents the first-order radial derivative. The Gaussian curvature chart in (d) illustrates its G 2 continuity.

Figure 8 :

 8 Figure 8: Capping a loop-shaped five-sided hole with a G 2 continuous periodic NURBS surface. (a) shows the control points of the boundary surfaces (chromatic) and the filling surface (white), respectively; the Gaussian curvature chart in (b) presents G 2 continuity; the following comparison is between our algorithm (c) and Ref. [SYSP11] (d).

Figure 9 :

 9 Figure 9: The tip of each blade of the propeller in (a) is filled with a single G 2 polar NURBS surface (see the continuous reflection stripes in (c)) which has non-uniform knots in the periodic direction (denser near the sharp part), (40 × 6) control points (c) and degrees (8, 5).


Hai-Chuan Song

Xin Xu

Kan-Le Shi

Jun-Hai Yong

Projecting points onto planar parametric curves by local biarc approximation

Keywords: point project, parametric curves, biarc interpolation, local approximation et al. [15]

This paper proposes a geometric iteration algorithm for computing point projection and inversion on planar parametric curves based on local biarc approximation. The iteration begins with initial estimation of the projection of the prescribed test point. For each iteration, we construct a biarc that locally approximates a segment on the original curve starting from the current projective point. Then we compute the projective point for the next iteration, as well as the parameter corresponding to it, by projecting the test point onto this biarc. The iterative process terminates when the projective point satisfies the required precision. Examples demonstrate that our algorithm converges faster and is less dependent on the choice of the initial value compared to the traditional geometric iteration algorithms based on single-point approximation.

Introduction

Projection of a test point on a curve or surface aims to find the closest point, as well as the corresponding parameter, on the curve or surface. Specially, when the test point lies on the curve or surface, the problem of point projection becomes the problem of point inversion. This operation has been extensively used in geometric processing algorithms such as surface intersection [1], interactive object selection and shape registration [2,3,4]. Moreover, it is a fundamental component of the algorithms of curve and surface projection as well [5,6]. In this paper, we address the problem of point projection and point inversion on planar parametric curves. We provide a geometric iteration algorithm, which approximates a segment on the original curve by a biarc. Compared with traditional single-point approximation algorithms [7,8,4], our algorithm converges faster and is less dependent on the choice of the initial value.

Related work

The problem of point projection and inversion can be translated to solving the minimum distance equation (Q-P)×n = 0, where P is the prescribed test point, Q is the point closest to P on the original curve or surface and n is the normal vector of the original curve or surface at 25 Q. In most of the early work, Newton-Raphson method, 26 which involves the first and second order derivatives, 27 was used to solve this minimum distance equation and 28 get the projective point [1,8]. Piegl and Tiller [9]gavea 29 detailed description on this method for point projection 30 and inversion.

31

In order to achieve a good initial value, which is im-32 portant for Newton-Raphson method to converge reli-33 ably, subdivision methods were introduced [10, 11, 12, 34 [START_REF] Johnson | Distance extrema for spline models 587 using tangent cones[END_REF][START_REF] Oh | Efficient point 590 projection to freeform curves and surfaces[END_REF][START_REF] Chen | Computing 593 the minimum distance between a point and a NURBS curve[END_REF][START_REF] Selimovic | Improved algorithms for the projection of points 596 on NURBS curves and surfaces[END_REF]. The key point of this kind of algorithms 35 is to eliminate the curve segments or the surface patches 36 which do not contain the nearest points. Ma and He-37 witt [START_REF] Ma | Point inversion and projection for NURBS 584 curve and surface: Control polygon approach[END_REF] divided the NURBS surface into several Bézier 38 patches and checked the relationship between the test 39 point and the control point nets of these Bézier patches.

40

However their elimination criterion may fail in some 41 cases [START_REF] Chen | A coun-599 terexample on point inversion and projection for NURBS curve[END_REF]. Johnson and Cohen utilized the tangent cone 42 to search for the portions of the surface contain the pro-43 jective points [START_REF] Johnson | Distance extrema for spline models 587 using tangent cones[END_REF]. A more practical exclusion criterion 44 based on Voronoi cell test was proposed in [START_REF] Selimovic | Improved algorithms for the projection of points 596 on NURBS curves and surfaces[END_REF]. Chen

Circle P Q 0 Q 4 Q 3 Q 2 Q 1 Biarc T angent Q 5 Q 6 (a) Circle Q 0 Q 4 Q 3 Q 2 Q 1 Biarc T angent (b)
Figure 1: A comparison of the approximation precision of the first order algorithm [7,8], the second order algorithm [4] and our algorithm: P is the test point, the black curve is the original curve, Q 0 is the initial point, orange point Q 6 , yellow point Q 5 , blue point Q 2 and red point Q 4 are projective points obtained by Newton-Raphson method [9], the first order algorithm [7,8], the second order algorithm [4] and our algorithm after the first iteration, respectively. Q 3 is the exact closest point. (a) the whole view of the projection; (b) the zoom view of (a). query points onto planar spiral curves [18]. Seong et al.

[19] dealt with this problem in another way. By elevating the dimension of the problem, they transformed the point projection onto planar parametric curve into the intersection of an implicit surface and a straight line.

Besides algebraic methods (Newton-Raphson method), geometric methods were also proposed, which only involve the geometric information that is common to all possible parameterizations. Hoschek and Lasser [7], Hartmann [8] introduced the first order geometric iteration method. Hu and Wallner [4] proposed a second order geometric iteration method, in which they generated an osculating circle (a circle possessing the same curvature with the original curve at the osculating point) and projected the test point on it instead of the original curve or surface. Liu et al.

[20] improved their method by replacing the circle of normal curvature with a second order osculating torus patch to the surface.

Note that the algebraic method and geometric methods generally converge (if can converge) to the local minimum projective point nearest to the initial value [START_REF] Johnson | Distance extrema for spline models 587 using tangent cones[END_REF]. All the subdivision methods we introduced above [10,11,[START_REF] Ma | Point inversion and projection for NURBS 584 curve and surface: Control polygon approach[END_REF][START_REF] Johnson | Distance extrema for spline models 587 using tangent cones[END_REF][START_REF] Oh | Efficient point 590 projection to freeform curves and surfaces[END_REF][START_REF] Chen | Computing 593 the minimum distance between a point and a NURBS curve[END_REF][START_REF] Selimovic | Improved algorithms for the projection of points 596 on NURBS curves and surfaces[END_REF] can generate the initial value near to the global minimum projective point. Therefore, generally, a complete point projection include two parts:

1. An algorithm to generate initial values; 2. An iteration algorithm to compute the precise projective point.

In this paper, we mainly focus on the iteration algorithm in the second part. It means that we assume the initial value has been provided by some kind of initial value 81 generating algorithm. Newton-Raphson method [9], first order algorithm 163 [7,8], or second order algorithm [4]). (1)

2. Newton-Raphson method [9], and

∆t = C ′ (t 0 ) • (C(t 0 ) -P) C ′′ (t 0 ) • (C(t 0 ) -P) + |C ′ (t 0 )| 2 .
(2)

3. First order algorithm [7,8], and

∆t = C ′ (t 0 ) • (Q -C(t 0 )) C ′ (t 0 ) • C ′ (t 0 ) , (3)
where Q is the projective point of P on the tangent 4. Second order algorithm [4], and

∆t = (Q -C(t 0 )) × C ′′ (t 0 ) κ C ′ (t 0 ) 3 , ( 4
)
where Q is the projective point of P on the osculat-201 ing circle at C(t 0 ), κ is the curvature of C(t 0 ), and

202 we have κ = (C ′ (t 0 ) × C ′′ (t 0 ))/ C ′ (t 0 )

246

The locus of all possible joint J is a circle (C in Fig- the "equal chord" biarc and the "parallel tangent" biarc.

251

The former one is constructed so that the two segments 252 P 0 J and JP 1 have equal arc lengths on C, while the latter one ensures that the tangent at point J is parallel to segment P 0 P 1 . In this paper we choose the "equal chord" biarc as the interpolation tools, owning to its simplicity in implementation and better approximation precision according to our experiments.

Given boundary data P 0 = C(t 0 ), D 0 = C ′ (t 0 ) and

P 1 = C(t 0 +∆t), D 1 = C ′ (t 0 +∆t),
we generate the biarc BA(s) in following steps:

1. Find the center of the joint circle C, by intersecting the bisectors of P 0 P 1 and of (P 0 + U 0 )(P 1 + U 1 ),

where U 0 and U 1 are the unit vectors of D 0 and D 1 , respectively.

2. Find J by intersecting the bisector of the P 0 P 1 and the minor arc on C bounded by P 0 and P 1 .

3. Construct the unique arcs A 0 , A 1 satisfying properties 1, 2 of Definition 1.

4. Generate a piecewise parametric biarc BA(s),

where s ∈ [0, 1]. It follows that BA(s) = A 0 when s ∈ [0, s J ], and
BA(s) = A 1 when s ∈ [s J , 1],
where

s J = ArcLength(A 0 )/(ArcLength(A 0 ) + ArcLength(A 1 )
).

An example is shown in Figure 1, where we use the second order algorithm to generate the approximation interval. In this figure, the biarc is much more closer to the original curve C(t) compared to the approximation geometries of the first and the second order algorithms.

It means that its approximation precision is higher than the first and the second order algorithms.

The reasons why we choose biarc as our approximation curve are as follows:

1. The construction of the biarc is simple, and can be done in constant time. There are four projective points on the two circles.

300

For the projective points we obtained, we only choose 301 one valid projective point by the following method: We apply the convergence criteria provided by Piegl and Tiller [9], which are

|(t i+1 -t 1 )C ′ (t i )|≤ε 1 .
(5)

|C(t i ) -P|≤ε 1 . (6) |C ′ (t i ) • (C(t i ) -P)| |C ′ (t i )||C(t i ) -P| ≤ ε 2 . (7)

Examples and comparisons

We present five examples for point projection, and make comparisons with NRA (Newton-Raphson method [9]), FOA (first order algorithm [7,8]), SOA (second order algorithm [4]), NRIBA (Newton-Raphson-integrated biarc algorithm: our algorithm whose ∆t is generated by Newton-Raphson algorithm [9] as introduced in subsection 2.1), FOIBA (first-orderintegrated biarc algorithm: our algorithm whose ∆t is generated by first order algorithm [7,8] as introduced in subsection 2.1), SOIBA (second-order-integrated biarc algorithm: our algorithm whose ∆t is generated by second order algorithm [4] as introduced in subsection 2.1).

All the experiments are implemented with Intel Core i5 CPU 3.0 GHz, 8G Memory. In all of our experiments ε 1 ,ε 2 are both the convergence tolerances introduced in subsection 3.2.

There are three main criteria to evaluate point projection iteration methods.

409

SOIBA converges faster than SOA, and the gap is get- ting larger as we decrease the tolerance.

410 Q P (a) Q P (b) Q P (c)

411

Example 2. We project point P = (381, 252) onto a 412 cubic Bézier curve, and we set t 0 = 0.53 (see Figure 6).

413

In this example, we set ε 1 = ε 2 = 10 -6 . Table 4 shows 414 the iteration steps of the six algorithms. (see the first steps of SOA in Table 4). This case always occurs at the special point whose curvature is relatively much bigger than its neighboring region, leading to a very small osculating circle, the approximation region of which is small. With the help of our local biarc approximation in SOIBA, the approximation region is enlarged, and the iteration can "jump" away from the special point (see the first two steps of SOIBA in Table 4), and converges to the correct projective point even with the "bad" initial value. In this example, all the single-point algorithms fail to some extent, while our integrated algorithms all converge to the correct projective point.

Example 3. We project 134 points on a spurious off- and the CPU time of our integrated algorithms are also less than the corresponding single-point algorithms.

Example 4. We project 100 points on a bounding circle (loose bounding) onto a hexagonal star (designed by hand, not an exact hexagonal star), which has sharp features (corresponding to the special points mentioned in Example 2) on its six angles (see Figure 8), using the six algorithms, respectively. In this example, we set ε 1 = ε 2 = 10 -10 , and the average initial value error is 5.48 × 10 -03 . The statistic data of the projection is shown in Table 6. Note that the average iterations only record the correct projections, and the CPU time records all the projections (including both correct and incorrect projections).

Experimental results show that, our SOIBA and NRIBA find all correct solutions, FOA fails in all projection, while the successful ration of FOIBA is 79.0%.

The successful ratios of other algorithms are NRA: 62.0%, SOA: 67.0%, respectively. The average number of iterations and the CPU time of our integrated algorithms are also less than the corresponding single-point algorithms.

Example 5. We project 266 points on the bounding and the CPU time of our integrated algorithms are also less than the corresponding single-point algorithms.

Example 4. We project 100 points on a bounding circle (loose bounding) onto a hexagonal star (designed by hand, not an exact hexagonal star), which has sharp features (corresponding to the special points mentioned in Example 2) on its six angles (see Figure 8), using the six algorithms, respectively. In this example, we set ε 1 = ε 2 = 10 -10 , and the average initial value error is 5.48 × 10 -03 . The statistic data of the projection is shown in Table 6. Note that the average iterations only record the correct projections, and the CPU time records all the projections (including both correct and incorrect projections).

Experimental results show that, our SOIBA and NRIBA find all correct solutions, FOA fails in all projection, while the successful ration of FOIBA is 79.0%.

The successful ratios of other algorithms are NRA: 62.0%, SOA: 67.0%, respectively. The average number of iterations and the CPU time of our integrated algorithms are also less than the corresponding single-point algorithms.

Example 5. We project 266 points on the bounding

82 1 . 2 .

 12 Our contributions 83 The main idea of geometric methods is to locally ap-84 proximate the original curve by a special curve (first or-85 der algorithm uses tangent line, and second order algo-86 rithm uses osculating circle). The next projective point 87 is estimated by projecting the test point onto the ap-88 proximation curve instead of the original curve. It is 89 shown in the evolution of the geometric methods that 90 higher approximation precision generally means higher 91 convergence speed and better stability. However, during 92 each iteration step, the traditional geometric algorithms 93 approximate the original curve only with curve deriva-94 tives computed at a single point. So the approximation 95 region is limited around this point, the approximation 96 precision will reduce when moving away from this point 97 on the original curve. 98 In order to improve the convergence speed and sta-99 bility of the point projection and inversion, we provide 100 a geometric iteration algorithm based on local biarc ap-101 proximation. Our method approximates the correspond-102 ing segment on the original curve by a biarc rather than 103 only one point during each iteration. Our local biarc ap-104 proximation has larger approximation region and higher 105 approximation precision compared to traditional single-106 point approximation. According to the experimental re-107 sults in Section 4, our algorithm converges faster and is 108 more robust than the traditional geometric algorithms.

  109

Figure 2 : 1 .

 21 Figure 2: Approximation deviation comparison in Figure 1. Abscissa is the parameter of the moving point on the original curve. Ordinate is the distance from the moving point to approximation curve.

146 2 . 3 .

 23 We present a framework that adapts to any single-147 point approximation algorithm. If a single-point 148 algorithm is integrated in our framework, the con-149 vergence speed and independence of the initial 150 value of this single-point algorithm will be im-Outline of our algorithm 153 Given a test point P, a planar parametric curve C(t) 154 and the parameter value t 0 of the roughly estimated pro-155 jective point Q 0 , as illustrated in Figure 1, we need to 156 compute the parameter of the exact projective point. 157 Our algorithm framework can be described in summary 158 as follows: 159 1. According to the initial projection parameter t 0 , 160 compute the interval width ∆t using any step length 161 strategy (we can use constant parameter increment, 162

164 2 .

 2 Compute the tangent vectors C ′ (t 0 ) and C ′ (t 0 +∆t), 165 respectively. Interpolate the boundary conditions 166 C(t 0 ), C ′ (t 0 ) and C(t 0 +∆t), C ′ (t 0 +∆t) with a biarc 167 BA(s) (the red curve in Figure1), which is used to 168 approximate the curve segment from C(t 0 )toC(t 0 + 169 ∆t) on the original curve C(t).

170 3 .

 3 Project the test point P onto the biarc BA(s) and 171 compute a new estimated parameter of the projec-172 tive point (Q 5 in Figure 1) on the original curve 173 C(t).

174 4 .

 4 Use the new parameter as the initial value t 0 and 175 repeat steps 1-3 until the corresponding projective 176 point satisfies the precision requirement. 177 The rest paper is organized as follows. Section 2 178 presents the method for local curve approximation by 179 biarc interpolation. In Section 3, the methods of point 180 projection onto the biarc and parameter inversion to the 181 original curve are described. The experimental results 182 including the evaluation of performance data are given 183 in Section 4. Finally, Section 5 concludes the paper.

184 2 . 1 .

 21 Local biarc approximation of a curve segment 185 We approximate a curve segment on the original 186 curve C(t) by a biarc according to the initial projection 187 parameter t 0 , which consists of the following steps: 188 Compute the interval width ∆t.189 2. Compute the interpolation boundary conditions 190 C(t 0 ), C ′ (t 0 ) and C(t 0 +∆t), C ′ (t 0 +∆t).

191 3 3 .

 3 Interpolate the boundary conditions with a biarc. 192 We present the three steps in details in the following 193 subsections, respectively.

194 2 . 1 .

 21 Compute the interval width 195 The interval width ∆t determines which curve seg-196 ment on C(t) is to be approximated. We select any step 197 length strategy from following ones: 198 1. User-defined constant parametric increment, and ∆t = const t.

199

  line at C(t 0 ).

  200

247 ure 3 )

 3 passing through P 0 and P 1 [21]. Various biarc 248 interpolation schemes are distinguished by the choice 249 of the joint J. Among the most important ones are 250

2 . 3 . 3 . 3 . 1 . 1 Figure 4 :

 2333114 Figure 4: Point projection on the biarc.

302 1 .

 1 If there is only one projective point in the paramet-303 ric domain [0, 1], this point is chosen as the valid 304 projective point.

305 2 . 3 . 3 . 2 .

 2332 If there are more than one projective point in the 306 parametric domain [0, 1], we choose the projective 307 point, which is closest to the test point P, from the 308 projective points in the parametric domain [0, 1] as 309 the valid projective point. 310 If there is no projective point in the parametric do-311 main [0, 1], we choose the projective point, whose 312 parameter is closest to the boundary of the para-313 metric domain [0, 1], as the valid projective point. 314 The parameter of the valid projective point is recorded 315 in variable param biarc. 316 Estimate the next projective point 317 After we derived the projective parameter 318 param biarc on the biarc, we refine ∆t with it. 319 Recall that the curve segment BA(0) ∼ BA(1) on 320 the biarc is corresponding with the curve segment 321 C(t 0 ) ∼ C(t 0 +∆t) on the original curve. So ∆t can 322 be refined by ∆t = param biarc × ∆t (note that, for 323 extreme high-order Bézier curves, this linear parameter 324 interpolation may be unstable). Then the parameter of 325 the next projective point is estimated by t 0 = t 0 +∆t, 326 which will be used for the next iteration.

  327

Figure 6 :

 6 Figure 6: Illustration of Example 2. The yellow point is the initial projection. The control points of the Bézier curve are {(0, 0), (110, 1000), (90, 1000), (200, 0)}: (a) projection result Q of NRA and its tangent; (b) projection result Q of SOA and its tangent; (c) projection result Q of SOIBA.

415Figure 6

 6 Figure 6 (c)). This is because (Q -C(t 0 )) is nearly par-441

Figure 7 :

 7 Figure 7: Illustration of Example 3: point projections on a smooth curve.

516Figure 8 :

 8 Figure 8: Illustration of Example 4: point projections on a star with sharp features.

Figure 9 :

 9 Figure 9: Illustration of Example 5: point projections on font characters.

Table 1 :

 1 Convergence comparisons for Example 1.

	P 1 = (381, 252), t 0 = 0.75
	Step	1234	56 C P U time (ms)
	∆t NRA	2.01 × 10 -02	-4.38 × 10 -04	-2.81 × 10 -07	0.00	t=0.769514		1.15 × 10 -02
	∆t FOA	3.24 × 10 -02	-2.38 × 10 -02	1.88 × 10 -02	-1.43 × 10 -02	1.13 × 10 -02	-8.75 × 10 -02	1.54 × 10 -02
	∆t SOA	2.07 × 10 -02	-1.21 × 10 -03	-4.13 × 10 -06	0.00	t=0.769514		2.46 × 10 -03
	∆t NRIBA	1.95 × 10 -02	3.02 × 10 -06	0.00	t=0.769514			2.11 × 10 -03
	∆t FOIBA	1.97 × 10 -02	-1.89 × 10 -04	0.00	t=0.769514			2.21 × 10 -03
	∆t SOIBA	1.95 × 10 -02	6.59 × 10 -06	0.00	t=0.769514			2.96 × 10 -03
	P 2 = (332, 200), t 0 = 0.5
	Step	1234	56 C P U time (ms)
	∆t NRA	1.53 × 10 -01	-2.87 × 10 -02	-2.29 × 10 -03	-2.05 × 10 -05	0.00	t=0.6223419	1.65 × 10 -02
	∆t FOA	8.53 × 10 -02	2.86 × 10 -02	6.88 × 10 -03	1.29 × 10 -03	2.25 × 10 -04	3.86 × 10 -05	3.46 × 10 -02
	∆t SOA	1.21 × 10 -01	1.24 × 10 -03	-3.76 × 10 -06	0.00 t=0.6223419		2.98 × 10 -03
	∆t NRIBA	1.19 × 10 -01	3.84 × 10 -03	3.08 × 10 -07	0.00 t=0.6223419		2.41 × 10 -03
	∆t FOIBA	1.27 × 10 -01	-4.33 × 10 -03	-4.81 × 10 -06	0.00 t=0.6223419		2.89 × 10 -03
	∆t SOIBA	1.23 × 10 -01	-2.02 × 10 -04	0.00	t=0.6223419			2.96 × 10 -03

t i is the parameter obtained at the ith iteration, and ε 1 ,ε 2 are two zero tolerances of Euclidean distance and cosine. The iteration is halted if any of the three conditions above is satisfied.

Table 2 :

 2 Comparisons of the number of iterations for Example 1, where P 1 = (381, 252).

		376
		377	P 2 = (332, 200) are projected onto a cubic B-spline
		378	curve C(t) with initial parameter 0.75 and 0.5, respec-
		379	tively. In this example, we set ε 1 = ε 2 = 10 -6 , which is
		380	the same with [4].
	1. Correctness. If the distance between the computed	381	In Table 1, we compare the convergence of the six
	projective point and the exact closest point satisfies	382	algorithms. If any method can converge within 10 iter-
	a given precision, it is treated as a correct solution.		ations, we will provide the convergence parameter in
		385

2. Speed of convergence. We measure the convergence speed by two kinds of experimental data: the number of iterations and the CPU time. We record 356 the average and the worst numbers of iterations in 357 each computation. 358 3. Independence on the initial value. The initial value 359 has a significant impact on the correctness of the 360 Newton-like iteration algorithms. Although there 361 were lots of methods finding initial value (as de-362 scribed in subsection 1.1), no method can claim 363 that its initial value is good enough to make the it-364 eration always converge. So if the iteration method 365 is less dependent on the initial value, the method 366 will be more robust. 367 In the following examples, all the initial values are 368 set by hand (except Example 1 which uses the same 369 initial value in [4] to compare with its second or-370 der algorithm) in order to test the performance of 371 different iteration methods. In practice, we recom-372 mend to estimate the initial value with the method 373 introduced in [15]. 374 Example 1. We first test Example 2 of [4] (see 375 Figure 5), where two test points P 1 = (381, 252) and 383 this

table. As shown in this table, our integrated al-384 gorithms (NRIBA, FOIBA and SOIBA) converge with 386 rithms (NRA, FOA and SOA). The processing times of 387 SOA, NRIBA, FOIBA and SOIBA are comparable, and 388 are less than those of NRA and FOA. 389 In Table 2, we compare the robustness of the six algo-390 rithms with respect to the choice of the initial parameter 391 t 0 . We use P 1 and C(t) in Example 1, and choose 101 392 different t 0 from 0, 0.01, 0.02, ..., 1. Table 2 shows the 393 number of iterations to converge correctly for the six 394 algorithms. As shown in this table, our integrated algo-395 rithms averagely converge with less iterations than the 396 corresponding single-point algorithms. 397 In

Table 3 ,

 3 we compare the convergence of the six al-

	398
	399	gorithms under different tolerances. As shown in this ta-
	400	ble, when we decrease the tolerance, the numbers of it-
	401	erations of FOA, SOA, FOIBA increase; however NRA,
	402	NRIBA and SOIBA remain. NRA performs better than
	403	FOA and SOA, because the initial value is not too far
	404	from the exact projective point. However, NRIBA con-
	405	verges much faster than NRA. FOIBA converges much
	406	faster than FOA. In the third tolerance case, the number
	407	of iteration of FOA jumped significantly from 66 to 91.
		Influenced by FOA, the number of iteration of FOIBA

408

jumped from 5 to 18, which is still relatively acceptable.

Table 3 :

 3 Comparisons of number of iterations / CPU time (ms) under different tolerances for Example 1, where P 1 = (381, 252), t 0 = 0.75.

	Tolerance (ε 1 = ε 2 )1 0 -6	10 -7	10 -8	10 -9	10 -10
	∆t NRA	5 / 1.15 × 10 -02	5 / 1.17 × 10 -02	5 / 1.11 × 10 -02	5 / 1.21 × 10 -02	5 / 1.22 × 10 -02
	∆t FOA	57 / 1.54 × 10 -02	66 / 1.95 × 10 -02	91 / 2.96 × 10 -02	101 / 3.43 × 10 -02	112 / 3.64 × 10 -02
	∆t SOA	5 / 2.46 × 10 -03	5 / 2.97 × 10 -03	11 / 6.48 × 10 -03	11 / 6.51 × 10 -03	11 / 6.50 × 10 -03
	∆t NRIBA	4 / 2.11 × 10 -03	4 / 2.14 × 10 -03	4 / 2.12 × 10 -03	4 / 2.21 × 10 -03	4 / 2.23 × 10 -03
	∆t FOIBA	4 / 2.21 × 10 -03	5 / 2.35 × 10 -03	18 / 1.50 × 10 -02	18 / 1.51 × 10 -02	18 / 1.50 × 10 -02
	∆t SOIBA	4 / 2.96 × 10 -03	4 / 3.06 × 10 -03	4 / 4.06 × 10 -03	4 / 4.07 × 10 -03	4 / 4.06 × 10 -03

Table 4 :

 4 Convergence comparisons for Example 2, where P = (381, 252).

	S t e p12345	6	7
	∆t NRA	-4.32 × 10 -02	3.98 × 10 -04	-4.98 × 10 -08	0.00	t=0.4872014
	∆t FOA	2.49	-7.64 × 10 -02	-7.05 × 10 -03	-6.12 × 10 -05	0.00 t=0.9164463
	∆t SOA	-3.18 × 10 -02	1.30 × 10 -02	1.36 × 10 -03	4.71 × 10 -05	-1.08 × 10 -07	0.00	t=0.5126524
	∆t NRIBA	8.78 × 10 -02	2.14 × 10 -01	8.44 × 10 -02	-9.45 × 10 -05	0.00 t=0.9164463
	∆t FOIBA	3.11 × 10 -01	7.53 × 10 -02	8.64 × 10 -04	-8.98 × 10 -07	0.00 t=0.9164463
	∆t SOIBA	1.16 × 10 -01	1.85 × 10 -01	8.64 × 10 -02	-3.90 × 10 -04	0.00 t=0.9164463

Table 5 :

 5 Statistic data for Example 3. = ε 2 = 10 -10 , and the average initial value er-

		Methods	Correct	Worst	Average	CPU time (ms)
			solutions iterations iterations
		NRA	118	35	4.38	1.47
		FOA	91	179	18.10	2.27
		SOA	131	7	4.25	1.29
		NRIBA	131	11	4.03	0.44
		FOIBA	127	37	5.27	1.38
		SOIBA	134	5	3.23	0.37
	458	set of a smooth curve onto the curve itself (see Figure 7),
	459	using the six algorithms, respectively. In this example,
	we set ε 1 460
	461	ror is 2.27 × 10 -02 . The statistic data of the projection is
		shown in

Table 5 .

 5 Note that the average iterations only

	462
	463	record the correct projections, and the CPU time records
	464	all the projections (including both correct and incorrect
	465	projections).
	466	Experimental results show that, our SOIBA finds all
	467	correct solutions, while the successful ratio of SOA is
	468	97.8%, even if there is no special points mentioned in
	469	Example 2 (because all initial values are not too far
	470	from the exact projective points and there is no sharp
	471	features). The successful ratios of other algorithms are
	472	NRA: 88.1%, FOA: 67.9%, NRIBA: 97.8%, FOIBA:
	473	94.8%, respectively. The average number of iterations

Table 6 :

 6 Statistic data for Example 4.

	Methods	Correct	Worst	Average	CPU time (ms)
		solutions iterations iterations
	NRA	62	8	6.77	5.35
	F O A0--	-
	SOA	67	176	19.11	7.66
	NRIBA	100	5	4.7	0.24
	FOIBA	79	43	5.98	4.85
	SOIBA	100	5	4.04	0.43

Table 7 :

 7 Statistic data for Example 5.

		Methods	Correct	Worst	Average	CPU time (ms)
			solutions iterations iterations
		NRA	237	8	4.77	3.03
		FOA	75	167	11.23	5.38
		SOA	256	9	4.72	0.782
		NRIBA	251	6	4.38	1.41
		FOIBA	226	30	5.98	0.79
		SOIBA	266	6	3.63	0.72
		Figure 9: Illustration of Example 5: point projections on font charac-
		ters.
	497	box (loose bounding) of several characters (represented
		by B-spline curves) onto these characters (see Figure 9),
	506
	507	correct solutions, while the successful ratio of SOA is
	508	96.2%, even if there are no special points mentioned
	509	in Example 2 (because all initial values are not too far
	510	from the exact projective points and there is no sharp
	511	features). The successful ratios of other algorithms are
	512	NRA: 89.1%, FOA: 28.2%, NRIBA: 94.7%, FOIBA:
	513	85.0%, respectively. The average number of iterations
	514	and the CPU time of our integrated algorithms are also
		less than the corresponding single-point algorithms.

498

using the six algorithms, respectively. In this example, 499 we set ε 1 = ε 2 = 10 -10 , and the average initial value 500 error is 7.01 × 10 -02 . The statistic data of the projection 501 is shown in Table 7. Note that the average iterations 502 only record the correct projections, and the CPU time 503 records all the projections (including both correct and 504 incorrect projections). 505 Experimental results show that, our SOIBA finds all

Table 6 :

 6 Statistic data for Example 4.

	Methods	Correct	Worst	Average	CPU time (ms)
		solutions iterations iterations
	NRA	62	8	6.77	5.35
	F O A0--	-
	SOA	67	176	19.11	7.66
	NRIBA	100	5	4.7	0.24
	FOIBA	79	43	5.98	4.85
	SOIBA	100	5	4.04	0.43

Table 7 :

 7 Statistic data for Example 5.

	Methods	Correct	Worst	Average	CPU time (ms)
		solutions iterations iterations
	NRA	237	8	4.77	3.03
	FOA	75	167	11.23	5.38
	SOA	256	9	4.72	0.782
	NRIBA	251	6	4.38	1.41
	FOIBA	226	30	5.98	0.79
	SOIBA	266	6	3.63	0.72

Table 7 .

 7 Note that the average iterations

	506
	507	correct solutions, while the successful ratio of SOA is
	508	96.2%, even if there are no special points mentioned
	509	in Example 2 (because all initial values are not too far
	510	from the exact projective points and there is no sharp
	511	features). The successful ratios of other algorithms are
	512	NRA: 89.1%, FOA: 28.2%, NRIBA: 94.7%, FOIBA:
	513	85.0%, respectively. The average number of iterations
	514	and the CPU time of our integrated algorithms are also
	515	less than the corresponding single-point algorithms.
		5. Conclusion
	518

502 only record the correct projections, and the CPU time 503 records all the projections (including both correct and 504 incorrect projections). 505 Experimental results show that, our SOIBA finds all 516

We present a geometric iteration algorithm to com-517 pute the projection and inversion of a point onto pla-

Acknowledgements 543 Acknowledgements 543

approximate the curve locally, and this local approxi-520 mation achieves both higher precision and larger fitting 521 region as compared to the single-point approximation 522 using tangents or osculating circles [7,8,4]. Given the