Kai-Mo Hu

Bin Wang

Jun-Hai Yong

Jean-Claude Paul

Relaxed lightweight assembly retrieval using vector space model

Keywords:

Assembly searching technologies are important for the improvement of design reusability. However, existing methods require that assemblies possess high-level information, and thus cannot be applied in lightweight assemblies. In this paper, we propose a novel relaxed lightweight assembly retrieval approach based on a vector space model (VSM). By decomposing the assemblies represented in a watertight polygon mesh into bags of parts, and considering the queries as a vague specification of a set of parts, the resilient ranking strategy in VSM is successfully applied in the assembly retrieval. Furthermore, we take the scale-sensitive similarities between parts into the evaluation of matching values, and extend the original VSM to a relaxed matching framework. This framework allows users to input any fuzzy queries, is capable of measuring the results quantitatively, and performs well in retrieving assemblies with specified characteristics. To accelerate the online matching procedure, a typical parts based matching process, as well as a greedy strategy based matching algorithm is presented and integrated in the framework, which makes our system achieve interactive performance. We demonstrate the efficiency and effectiveness of our approach through various experiments on the prototype system.

Introduction

Over the last twenty years, Computer Aided Design (CAD) systems have become popular in product development organizations, and this results in the creation of large databases of assemblies. For such organizations, the efficient management and effective reuse of these assemblies become a key factor to stay ahead in the completion.

Searching technologies of mechanical components are the key to getting the required assemblies. However, to the best of our knowledge, existing search methods mainly focus on part retrieval [START_REF] Hong | Similarity comparison of mechanical parts to reuse existing designs[END_REF][START_REF] Bespalov | Local feature extraction and matching partial objects[END_REF][START_REF] Jayanti | Shape-based clustering for 3D cad objects: a comparative study of effectiveness[END_REF][START_REF] Bai | Design reuse oriented partial retrieval of CAD models[END_REF][START_REF] Iyer | Three-dimensional shape searching: state-of-the-art review and future trends[END_REF]. Though part retrieval methods can be applied to find such assemblies whose overall shapes are similar, they all fail to take into consideration relationships, compositions and structures that exist in the assemblies. Recently, some assembly retrieval methods have been proposed [START_REF] Deshmukh | Content-based assembly search: a step towards assembly reuse[END_REF][START_REF] Chen | A flexible assembly retrieval approach for model reuse[END_REF]. They assume that the assemblies possess full information, such as the topological structures, the orientation relationships and the joint information, and thus provide some ingenious mechanisms to retrieve assemblies with characterized topologies or features. However, if the product data in PDM or PLM systems is in lightweight file format (e.g. watertight polygon mesh), the existing methods cannot be applied.

To provide a practical way for the retrieval tasks of assemblies represented in lightweight file format, we propose a new framework, in which only the geometrical information of parts and composition information in assemblies is used. The VSM, which has been used successfully in the information retrieval (IR) field, is applied to get the partial matching and quantitative ranking of assemblies. To overcome the limitation of VSM that fuzzy queries are not supported, we extend the traditional VSM to a relaxed matching framework by integrating the similarities between parts into the evaluation of matching values, and thus provide more flexible means for query specifications and result rankings.

Since the efficiency is essential for online retrieval, we design a typical parts based matching process as well as a greedy strategy based matching algorithm, and apply them in our relaxed retrieval framework. Experiment results demonstrate the efficiency and effectiveness of our work in a variety of manners, and verify its utility in retrieving the lightweight assemblies with desired characteristics.

Additionally, we propose a normal compatibility based assembly decomposition algorithm and a scale-sensitive similarity based parts merging algorithm. Using these two algorithms, the assemblies represented in lightweight formats can be segmented and merged as a set of representative single parts accurately, such that the extended VSM can be applied directly.

In summary, we make the following contributions:

• Propose a normal compatibility based assembly decomposition algorithm, as well as a scale-sensitive similarity based parts merging algorithm, so as to disassemble the assemblies into sets of representative parts.

• Apply the classical VSM into the assembly retrieval, and extend the traditional VSM to a relaxed matching framework to support fuzzy retrieval.

• Present a typical parts based matching process and a greedy strategy based matching algorithm for the system acceleration. By integrating them into the matching algorithm, our system could achieve interactive performance.

Compared to the state-of-the-art work, we provide a well defined mechanism for the retrieval of assemblies represented in polygon mesh format. By introducing and extending the classical VSM in IR into assembly retrieval, our proposed framework allows users to retrieve assemblies through inputting parts sets, and is able to measure the retrieval results quantitatively.

The remainder of this paper is organized as follows. Section 2 reviews the relevant research on assembly retrieval, followed by an overview of our framework in Section 3. In Sections 4 and 5, we describe the preprocessing algorithms and online retrieval algorithms respectively, and then the experiment results and some discussions are presented in Section 6. Finally, we conclude our work and give some future work in Section 7.

Pre knowledge and related work

Related concepts in information retrieval

Since we try to apply the VSM in the assembly retrieval task, some related concepts in IR are presented. Generally, the classical models consider that each document is described by a set of representative keywords called index terms. An index term is simply a word whose semantics helps in remembering the documents' main themes. Since the distinct index terms have varying relevance, this effect is captured through the assignment of numerical weights to each index term. Based on different promises, three classic models, namely the Boolean model [START_REF] Wartick | Boolean operations[END_REF], the probabilistic model [START_REF] Robertson | Relevance weighting of search items[END_REF] and the VSM [START_REF] Raghavan | A critical analysis of vector space model for information[END_REF], are regarded as the basic models in IR systems.

The hypothesis of VSM is that two documents are similar if they contain some of the same terms. By assigning non-binary weights to index terms in queries and in documents, the degree of similarity between each document is stored in the system and the user query can be computed to support partial matching. Compared with the Boolean model and the probabilistic model, the main benefits of the VSM are: (1) its term-weighting scheme improves retrieval performance; (2) its partial matching strategy allows retrieval of documents that approximate the query conditions; and (3) its cosine ranking formula sorts the documents according to their degrees of similarity to the query. Theoretically, the VSM has the disadvantage that index terms are assumed to be mutually independent, which is not true in practice. However, modeling the dependencies between index terms is a very difficult task. Since most of the term dependencies are local, and only appear in a few special cases, if we apply them in applications indiscriminately, the overall performance does not improve significantly [START_REF] Ricardo | Modern information retrieval[END_REF].

In assembly retrieval, if we decompose the assemblies into single parts, and treat an assembly as a set of representative parts, then each assembly can be regarded as a document, and the parts in the assembly can be considered as index terms. By this transformation, the classic VSM can be applied to perform the retrieval task of assemblies. In addition, the geometrical information of single parts can be easily embedded into the matching evaluation. This framework is promising to provide an easy way to retrieve the required assemblies, and thus to access the associated data for the support of downstream manufacturing and service operations.

Retrieval methods for models and assemblies

Due to its wide applications in a variety of areas, 3D model retrieval has drawn much attention in recent years. Some early research was aimed at proposing powerful descriptors, such as Shape Distribution [START_REF] Osada | Shape distributions[END_REF], Shape Histogram [START_REF] Ankerst | Nearest neighbor classification in 3D protein databases[END_REF], Extended Gaussian Image [START_REF] Horn | Extended Gaussian images[END_REF], Spherical Harmonic Descriptor [START_REF] Kazhdan | Rotation invariant spherical harmonic representation of 3D shape descriptors[END_REF] and Light Field Descriptor (LFD) [START_REF] Chen | On visual similarity based 3D model retrieval[END_REF]. Shilane et al. [START_REF] Shilane | The princeton shape benchmark[END_REF] tested some of the existing descriptors on the Princeton Shape Benchmark, and found that LFD provides the best retrieval precision. With the development of mesh segmentation technologies, some part-in-whole shape retrieval algorithms have been proposed [START_REF] Gal | Salient geometric features for partial shape matching and similarity[END_REF][START_REF] Funkhouser | Partial matching of 3D shapes with priority-driven search[END_REF]. These methods first extract the descriptors of local regions, and then match them using some special techniques. Recently, some deformationinvariant shape descriptors have been introduced [START_REF] Zhang | Deformation-driven shape correspondence[END_REF][START_REF] Rustamov | Laplace-Beltrami eigenfunctions for deformation invariant shape representation[END_REF][START_REF] Gal | Pose-oblivious shape signature[END_REF]. Since the pair-to-pair matching of local features is time consuming, the bag-of-words method proposed in IR has been employed in recent 3D model retrieval tasks. For instance, Bronstein et al. [START_REF] Bronstein | Shape google: geometric words and expressions for invariant shape retrieval[END_REF] used multi-scale diffusion heat kernels as ''geometric words'', and constructed shape descriptors by means of ''bag of features''. This method is efficient and can deal with the deformable objects.

In CAD systems, the retrieval tasks are performed to improve the reusability of existing designs, and the work in literature can be classified into two categories. The first category tries to discover and extract common design patterns from databases. For example, Ma et al. [START_REF] Ma | Extracting common design patterns from a set of solid models[END_REF][START_REF] Ma | Automatic discovery of common design structures in CAD models[END_REF] developed some approaches to extract common local structures as design patterns from boundary representation (B-rep) models. They first translated the models into some graphs, and then extracted the frequent subgraphs as common design structures. The second category targets partial retrieval of CAD models, such that new parts do not have to be designed from scratch. Hong et al. [START_REF] Hong | Similarity comparison of mechanical parts to reuse existing designs[END_REF] presented a two-step similarity comparison method for B-rep, the overall appearances are compared first, and then the detailed features are matched. To solve the partial matching problem of CAD models, Bespalov et al. [START_REF] Bespalov | Local feature extraction and matching partial objects[END_REF] proposed a scale-space feature extraction technique based on the recursive decomposition of polyhedral surfaces into surface patches. Their idea of decomposition and retrieval is in some way similar to our framework. Another partial matching work was reported by Bai et al. [START_REF] Bai | Design reuse oriented partial retrieval of CAD models[END_REF]. They first defined the criteria for determining whether a subpart of CAD models is reusable, and then extracted and stored the reusable subparts in the library for later retrieval. El-Mehalawi and Miller [START_REF] El-Mehalawi | A database system of mechanical components based on geometric and topological similarity. Part I: representation[END_REF][START_REF] El-Mehalawi | A database system of mechanical components based on geometric and topological similarity. Part II: indexing, retrieval, matching, and similarity assessment[END_REF] provided the representation scheme of CAD models, including representation, indexing, retrieval, matching and similarity assessment.

In recent years, with the wide usage of PDM, more and more assemblies have been accumulated. To reuse the design knowledge, some researchers propose to retrieve similar cases when designing new assemblies. For instance, Kim et al. [START_REF] Kim | Case-based design for assembly[END_REF] described the case-based reasoning (CBR) method. By searching through the case library, the previous designs whose identified features are similar to the current case can be found. However, the knowledge in the case library only contains product features and process features. Based on CBR, a recent work on producing interesting redesigns through general design knowledge was reported in [START_REF] Chen | A framework for an automotive body assembly process design system[END_REF]. Chao et al. [START_REF] Chao | Case retrieval in body-in-white parts based on similarities of welding and assembly process[END_REF] also proposed a case retrieval method, in which the welding and assembly processing information are employed. Provided with assembly drawings, Liu et al. [START_REF] Liu | Component parts extraction from assembly drawings for content based retrieval[END_REF][START_REF] Liu | Device parts retrieval from assembly drawings with svm based active relevance feedback[END_REF] proposed to extract component parts, and then perform the assembly drawing Fig. 1. The overall workflow and sub-procedures of our relaxed retrieval framework. The first row shows the preprocessing phase, while the 2-4 rows depict the online retrieval phase. For each sub-procedure, we present the action and its corresponding description section in our paper around the arrow.

retrieval by matching the extracted component parts. The support vector machine based relevance feedback was adopted to improve the performance. The idea of extracting component parts and performing partial retrieval resembles ours in some way, but they can only deal with 2D assembly drawings, and fuzzy retrieval is not supported.

The approaches that are directly related to assembly retrieval are introduced by Deshmukh et al. [START_REF] Deshmukh | Content-based assembly search: a step towards assembly reuse[END_REF] and Chen et al. [START_REF] Chen | A flexible assembly retrieval approach for model reuse[END_REF]. In [START_REF] Deshmukh | Content-based assembly search: a step towards assembly reuse[END_REF], the search criteria, search algorithm and search strategies for content-based assembly search are stated in detail, and a simple system is developed to support content-based searches. Though the five assembly search examples have illustrated the promising power of assembly retrieval, their current implementation mainly focuses on the global information of the assemblies, such as the names of the constituent parts, the mating conditions between different parts, the material and mass properties of the parts and so on, and the quantitative rankings are not provided. To overcome the limitations of existing methods, Chen et al. [START_REF] Chen | A flexible assembly retrieval approach for model reuse[END_REF] proposed a flexible and effective approach, in which the multilevel descriptor, the hierarchical assembly structure and the semantic assembly interface are presented to preserve the implicit highlevel design semantics. Additionally, they provided an efficient indexing mechanism to accelerate the retrieval process.

While these two methods provide goodish solutions for assembly retrieval, they cannot be applied directly to the assemblies represented in lightweight file format, since in such format, only the low-level geometrical information is preserved, and the high-level feature information is lost. Unfortunately, to the best of our knowledge, most of the product data in PDM and PLM systems are stored in lightweight file format, so as to save storage space.

Overview

To solve the assembly retrieval problems existing in literature, we propose a relaxed assembly retrieval framework, targeting at solving the retrieval tasks of assemblies represented in a watertight polygon mesh. In this framework, we treat the assemblies as ''bags of parts'', and apply the well-designed ranking strategies of VSM in our retrieval algorithm. To support fuzzy queries, we take the similarities between parts into the evaluation of matching values, and thus extend the original VSM to a relaxed retrieval framework. A typical parts based matching process and a greedy strategy based matching algorithm are also integrated, so as to make our system achieve interactive performance. The overall workflow is illustrated in Fig. 1.

In the following sections, we first introduce the preprocessing algorithms of assemblies in Section 4, including the assembly decomposition algorithm and the parts merging algorithm, and then describe our relaxed matching framework in detail in Section 5, including the application of VSM, the typical parts based matching process and the greedy matching algorithm.

The preprocessing of assemblies

Normal compatibility based assembly decomposition

Suppose the parts are exported from CAD modeling systems in regular form, in which the dangling edges and dangling faces are not allowed, and one edge can only be adjacent to exactly two faces. Furthermore, we require that the mesh models must be watertight, and the normals of adjacent faces are compatible.

Based on these two constraints, we have designed an algorithm that can accurately decompose an assembly into a set of partial meshes, each of which corresponds to a single part. In this algorithm, a breadth-first search (BFS) strategy is applied to construct the decomposed mesh progressively. Generally, the faces in the assembly A can be regarded as the nodes in a graph, and the adjacent relationships between faces can be regarded as edges. In the decomposing process, we first construct the set of adjacent faces to u, and then for each adjacent face v, we check whether it is compatible to u. If so, it is placed into the mesh with u. After a complete traversal of BFS started from u, all the faces that belong to the same single part are marked as ''visited''. Hence, the newly constructed mesh M is added to the set S, and we continue to search the ''not visited'' faces in A until all the faces have been marked as ''visited''. The implementation details are shown in Algorithm 1.

The most time consuming step in Algorithm 1 is the construction of adjacent relationships between faces in A, which requires

Algorithm 1: Assembly decomposition algorithm

Input: A: the assembly represented by mesh; Output: S: the set that stores the decomposed meshes.

1 foreach face f in A do 2
Mark f as ''not visited''; ⊲Q is a queue. O(V 2 ) time complexity, and V is the number of faces in A. However, this decomposition algorithm can be performed offline, so it does not hurt the online retrieval performance at all.

7 ENQUEUE(Q , f ); ⊲ Place f at the tail of Q . 8 while Q = ∅ do 9 u ← DEQUEUE(Q ); ⊲ Extract the head of Q . 10 foreach v that is adjacent to u in A do 11 if v is marked as ''not

Scale-sensitive similarity based parts merging

Fig. 2 gives an example showing how the assemblies are decomposed into single parts using Algorithm 1. However, a single part may appear several times in the same assembly, such as the flanges appearing in Fig. 2. To recognize and merge the identical parts that have been separately decomposed by Algorithm 1, we propose to use robust 3D shape descriptors.

Since LFD [START_REF] Chen | On visual similarity based 3D model retrieval[END_REF] performs best according to Shilane's test [START_REF] Shilane | The princeton shape benchmark[END_REF], we employ their program to characterize the parts' global geometries. However, LFD returns small values if two models are similar. To make it fit for our algorithm, we normalize the similarity s LFD returned by Chen's program as follows:

s geometry (u,v) = ⎧ ⎨ ⎩ 1 - s LFD (u,v)-s min s max -s min if s max > s min s LFD if s max = s min (1)
where s min and s max are the minimum and maximum values returned from all the models by LFD. Using s geometry alone is not enough, since the relative scales between single parts are also important characteristics in assemblies. However, LFD is a scale-invariant descriptor. To solve this problem, we introduce another invariable s scale to measure the scale differences between two models:

s scale (u,v) = 1 - |r u -r v | r max -r min if r max > r min 1i f r max = r min (2)
where r u and r v are the radii of models u and v, while r min and r max are the minimum and maximum values of radii in the part database, respectively. Using the definitions of Eqs. ( 1) and ( 2), we define the scalesensitive similarity measurement as follows:

s(u,v) = s geometry (u,v) * s scale (u,v), (3)
which characterizes both the geometry and scale properties.

Based on this similarity measurement, a parts merging algorithm is designed in Algorithm 2. In this algorithm, we first initialize each part as an individual set, and then traverse every pair (u,v) in the set of decomposed parts. If the similarity between u and v measured by Eq. ( 3) is larger than a threshold, they are regarded as identical, and the two sets are merged in the case that they are previously in different sets. In the LFD program, the geometry(u,v)always returns the same value if two parts are the same. In our implementation, we consider two parts have the same scale if the difference of their scale(u,v) is less than 1%. Hence, the threshold ǫ was set to 0.99 in our experiment. The merging experiments also verified that this configuration could get good results.

Algorithm 2: Parts merging algorithm

Input: ǫ: threshold that identifies the sameness of parts; S: the set that stores the decomposed parts. Output: C: the set that stores the merged parts. The time complexity of Algorithm 2 is analyzed as follows.

Suppose the element number in S is V , the double loops in lines 4-15 take O(V 2 ) times. In our prototype system, the FIND-SET and UNION operations have been implemented on the disjoint-set forest, which take O [START_REF] Hong | Similarity comparison of mechanical parts to reuse existing designs[END_REF] and O(mα(n)) respectively. Here, α(n) is a very slowly growing function, and can be regarded as α(n) ≤ 4in any conceivable applications [START_REF] Cormen | Introduction to algorithms[END_REF]. So the upper bound of running time is O(V 3 ).

Since this algorithm can be executed offline, the high time complexity does not affect the online retrieval performance. In industrial practice, we first run Algorithm 2 online, and then store the merged results for later utilization. When a new part is added, we do not have to run Algorithm 2 again, but only need to scan in the merged parts. If there exists a part whose similarity with the new part is larger than the threshold ǫ, we merge the new part with that part; otherwise, the new part is added into the merged parts set directly. Fig. 3. Two assemblies that are represented in 3-dimensional part vector space (w 1 ,w 2 ,w 3 ). The similarity is measured by the angle between their vectors in the same part vector space.

Relaxed assembly retrieval framework

By considering assemblies as''bags of parts'', and regarding the query parts as a special type of assembly, the classic IR model, such as VSM can be applied directly to the assembly retrieval task. In addition, the geometrical information of single parts can be easily embedded into this retrieval framework. In this section, we will show how to apply the VSM in assembly retrieval, and how the application can be extended to a relaxed matching framework.

Vector space model in assembly retrieval

Similar to the document retrieval, we treat the assemblies as a collection C of parts and regard the user query as a vague specification of a set A of parts. In this scenario, the assembly retrieval problem can be reduced to the problem of determining which assemblies are in the set A and which ones are not. This can be partially viewed as a clustering problem, and two main issues have to be resolved. First, we need to determine what are the features that better describe the assemblies in the set A. Second, we need to determine what are the features that better distinguish the objects in the set A from the remaining parts in the collection C. The first set of features provides for quantification of intracluster similarity, while the second set of features provides for quantification of inter-cluster dissimilarity. We try to balance these two effects.

Basically, we have the hypothesis that two assemblies are similar if they contain some of the same parts. The possible measures of similarity that we might take into consideration are: (a) the number of parts in the assemblies; (b) the number of parts in common; (c) whether the parts are common or unusual and (d) how many times each part appears. To encode this idea, we transform the assemblies into an n-dimensional space, where n is the number of different parts used to index the set of assemblies.

An assembly j is represented by a vector d j = (w 1j ,w 2j ,...,w nj ), and its magnitude in dimension i is defined as w ij , where

w ij > 0, if part i occurs in assembly j; w ij = 0, otherwise.
Here, w ij denotes the weight of part i in assembly j.

By transforming the assemblies into part vector space, the similarity between two assemblies is defined as a function of the angle between their vectors in the part vector space. We take two assemblies represented in 3-dimensional part vector space as an example, as shown in Fig. 3. Generally, the correlation between vectors can be quantified by the cosine of the angle between these two vectors. That is,

sim(x 1 , x 2 ) = cos(θ ) = x 1 • x 2 |x 1 ||x 2 | (4)
where x 1 • x 2 is the inner product of vectors x 1 and x 2 .

Using Eq. ( 4), the similarity between two assemblies can be calculated easily. However, different assignments of part weights may result in different rankings. We will show how to achieve the best performance by balancing the effects of intra-cluster similarity and inter-cluster dissimilarity.

The intra-clustering similarity is quantified by measuring the raw frequency of a term k i inside an assembly d j , since a part that appears many times within an assembly is likely to be more important than a part that appears only once. Suppose the part i appears f ij times in assembly j, a simple method is to use f ij as the weight of part i. However, parts are likely to appear more often in complex assemblies. Therefore, f ij should be formalized by some variable related to the complexity of the assembly. A standard method in IR is to formalize f ij so that it is relative to the frequencies of other parts in the assembly. Let m j = max(f ij ) in assembly j, then the Term Frequency (tf) is formalized as

tf ij = f ij /m j , when f ij > 0. ( 5
)
Since m j presents a simple way to reflect the complexity of the assembly, and thus normalizes the tf ij approximately, Eq. ( 5) works well in our experiments.

The inter-cluster dissimilarity is quantified by measuring the inverse of the frequency of a part i among all the assemblies in the collection. This factor is usually referred to as the inverse document frequency (idf) in IR. In our intuition, a part that occurs in a few assemblies is likely to be a better discriminator than a part that appears in most or all of the assemblies. To model this concept, we suppose there are n assemblies in database and the number of assemblies in which part i occurs is n i , then a possible method is to use n/n i as idf . However, this over-emphasizes small differences. Therefore, the logarithm is applied, and the idf can be revised as

idf i = log n n i + 1, when n i > 0. (6)
We add 1 to the logarithm such that the idf i is always positive.

Generally, the intra-clustering similarity and the inter-cluster dissimilarity depict the importance of a specified part from two different aspects. For an assembly, the higher value a part's tf and idf values are, the more important the part is. Thus, we can model the weighting scheme of a part by the product of tf and idf . Practical experience in IR has also demonstrated that this weighting scheme performs well in a wide variety of circumstances [START_REF] Ricardo | Modern information retrieval[END_REF]. By applying this observation, and using Eqs. ( 5) and ( 6), the weight of the i-th part in the j-th assembly can be formulized as

w ij = tf ij * idf i = f ij m j * log n n i + 1 , when n i > 0. (7)
In our VSM framework, a query is treated as a special type of assembly. By substituting Eq. ( 7) to Eq. ( 4), the similarity between query q and assembly d j is given by sim

(q, d j ) = cos(d q , d j ) = n k=1 w kq w kj |d q ||d j | , (8)
where d q and d j are the corresponding weighted part vectors. For the query part weight w kq , we use the suggestion of Salton and Buckley [START_REF] Salton | Term-weighting approaches in automatic retrieval[END_REF]:

w kq = 0.5 + 0.5 * f kq m q * log n n k + 1 , when f kq > 0, (9)
and for the assembly part, the weight w kj is given by Eq. ( 7).

Typical parts based matching process

Suppose there are m parts in the query and n parts in the assembly to be matched. For each online matching process, we need to extract the descriptors for all the m query parts, and then get their pair-to-pair similarities to the n assembly parts, which requires at least m * n similarity calculations. This is insupportable in practice, especially when the scale of the assembly database becomes large. To accelerate the online matching process, we propose a typical parts based matching process, as shown in the following two steps.

In the offline step, we select a number of typical parts, all of which are representative in CAD model databases, both in geometries and scales. Then, all the unique parts that participate in the assemblies are decomposed, and their pair-to-pair scalesensitive similarities to the typical parts are pre-computed and stored in the similarity table.

In the online step, given a part q provided by the user, we first search the most similar part t from the list of typical parts, and then get the similarity between t and the assembly part p by simply looking up in the similarity table. The final similarity between q and p can be approximately given by: s(q, p) = s(q, t) * s(t, p), [START_REF] Raghavan | A critical analysis of vector space model for information[END_REF] where s(q, t) is the scale-sensitive similarity between q and t, and s(t, p) is the scale-sensitive similarity between t and p, which has been calculated in the offline step. Since the count of query parts and typical parts are limited, s(q, t) can be calculated fast in the online step.

The relaxed matching algorithm

In most cases, we do not know what exact parts the assemblies contain, and still would like to find the assemblies with the desired characteristics by providing some fuzzy query parts. In this situation, simply applying the VSM into the assembly retrieval may result in poor performance, since for most of the fuzzy query parts, we may fail to find their corresponding parts in the assembly database, such that the participation relationships between the query parts and the assemblies can be constructed.

To provide a more universal framework for assembly retrieval, we propose to use the relaxed retrieval mode, which was introduced in [START_REF] Bai | Design reuse oriented partial retrieval of CAD models[END_REF]. In this mode, the query parts do not necessarily match the assembly parts exactly. Instead of using the Boolean values to indicate whether the query parts match the assembly parts exactly, we take the similarity measurement between query parts and assembly parts into consideration, and present a bipartite graph matching based framework to solve the relaxed retrieval problem.

We consider the m query parts and n assembly parts as the nodes in graph G, and the pair-to-pair similarities between query parts and assembly parts as edges. By constructing such a graph, the similarity measurement between the query and the assembly can be rendered as a bipartite graph matching problem. In this bipartite graph G, the weights of parts are assigned as the associated values of the nodes, and the scale-sensitive similarities between query parts and assembly parts are assigned as the weights of edges, as illustrated in Fig. 4. Suppose that the i-th part in the query is matched to the j-th part in the assembly, then the matching value is defined as:

S ij = w iq * s(i, j) * w jp , (11)
where s(i, j) is the approximate scale-sensitive similarity defined in Eq. ( 10), and w iq and w jp are the weights of the query part and assembly part, respectively. The final similarity can be computed according to a match M,in which the sum of matching values is maximized:

M = arg max (i,j) ΣS ij , 1 ≤ i ≤ m, 1 ≤ j ≤ n. ( 12
)
Given M, the similarity between query q and assembly d j can be defined similar to Eq. ( 8):

sim(q, d j ) = cos(d q , d j ) = S |d q ||d j | = ΣS ij |d q ||d j | , ( 13
)
where S is the summation of matching values according to M.

The bipartite graph matching problem can be solved using the Kuhn-Munkras (KM) algorithm. However, the time complexity of its naive implementation is O(n 4 ), and even the improved implementation still requires O(n 3 ). This is impractical for large assembly databases. To make the application achieve interactive performance for users, we propose and implement an approximate matching algorithm using greedy strategy. Suppose the number of query parts is smaller than that of the assembly parts, then the matching procedure can be designed as Algorithm 3.

Algorithm 3: Relaxed matching algorithm

Input: Q : the list of query parts; P: the list of assembly parts in assembly A. Output: S: the similarity between Q and A.

1 Calculate the weights of parts in Q using Eq. ( 9); 2 Sort the parts in Q by weights in descending order; 3 Mark all the parts in P as ''not matched''; In this algorithm, we first sort the parts in Q by their weights in descending order, so as to ensure the parts with larger weights are matched first. Then for each part in the sorted Q , we find such a part in the assembly that its matching value is maximized, as

4 S ← 0; 5 for i ← 1 to size(Q ) do 6 maxMatchValue ← 0, maxMatchPart ← 0;

Table 1

The snapshot of our benchmark. The ''Total parts in assemblies'' are all the parts that participate in the assemblies (with duplication), and the ''Shared parts'' are such parts that occur in at least two different assemblies. shown in lines 5-17. In our implementation, the matching value S ij is evaluated by Eq. ( 11).

If the number of query parts is larger than that of the assembly parts, we first sort the parts in P, and then seek the best match for each part in the sorted P, which is similar to the procedure presented in Algorithm 3. Suppose the number of parts in Q and P are m and n, respectively, then the time complexity of Algorithm 3 is O(mn). In practice, the number of query parts is usually small and can be regarded as a constant in most cases, so the matching process can be done in approximate O(n) time complexity, which is linear to the complexity of the assembly to be matched.

Experiment results and discussions

We implemented a prototype system, and the snapshot of the user interface is shown in Fig. 5. To make users easy to refine the fuzzy queries, the prototype system provides a function of retrieving similar parts. By providing an initial query part, users can get its scale-sensitive similar parts from a typical parts list by clicking the button ''Query Similar Parts''. If the required part is returned, users can add it to the query list, modify its frequency, and perform fuzzy retrieval by clicking ''Query Fuzzy Assemblies''. In the following subsections, we will first introduce the benchmark constructed in our experiment, and then present and discuss the relaxed retrieval results in several manners.

Experiment benchmark

There is no publicly available assembly database currently in the literature. To test the effectiveness of our approach, we gathered a large number of assemblies that had been modeled in industry, including household appliances, transportation, buildings, mechanical products and so on. 614 assemblies containing 5100 parts were finally selected as the test benchmark. By removing the duplications from the 5100 parts, we obtained 2814 unique parts for descriptor extraction. Then, the typical parts were selected from those 2814 unique parts in the following manner:

Initially, the scale-sensitive similarity defined in Section 3 was employed as similarity measurement, so that the pair-to-pair similarities among those unique parts can be computed. Then we applied the classical k-means clustering algorithm, and selected the mean of points in each cluster as typical parts. To get the parameter k and the initial partition before executing the k-means algorithm, we utilized Algorithm 2 to cluster similar parts, in which the threshold ǫ was set to 0.6. Table 1 gives the summary of our final benchmark after preprocessing.

Relaxed retrieval results

Comparison with the original VSM

To illustrate the advantages of our relaxed retrieval framework, a comparison with the performance of the original VSM in assembly retrieval is provided, and we call it exact retrieval for convenience. In this example, we take a flange that participates in seven assemblies in our benchmark as input, and then perform the exact retrieval and relaxed retrieval respectively. In the exact retrieval, the system only returns seven assemblies, in which the provided flange participates exactly, as shown in Fig. 6(a). On the contrary, the relaxed retrieval algorithm returns more relevant assemblies when provided with the same query, and Fig. 6(b) gives the corresponding retrieval result. We see that the top assemblies in Fig. 6(b) contain more flanges than that of Fig. 6(a), since the high frequencies of the flanges make their weights high, and some assemblies containing similar flanges in high frequencies are also ranked in the front. Some other comparisons have been conducted on our prototype system as well. Compared to the exact retrieval, our relaxed retrieval framework allows users to input any fuzzy queries, takes both the part similarities and assembly compositions into the evaluations of matching values, and thus provides more powerful means for assembly retrieval.

Retrieval performance on query variations

Similar to that of document retrieval, the variations of queries in the relaxed retrieval also lead to the changes of retrieval results. For simplicity, we again take the flange in Fig. 6 as an example, and furthermore, a pipe fitting is added into the query list, as shown in the left side of Fig. 7. Different from the result shown in Fig. 6(b), in which the assemblies are ranked by the number of similar flanges, most of the assemblies that contain two flanges and one pipe fitting are ranked in the front in Fig. 7.

We have tested a number of queries by modifying the number and kind of input parts. The conclusion is, the more kinds of parts we provide, the larger possibility that the required assembly is found, since more kinds of parts tend to characterize more features in the required assembly. Moreover, the relative frequencies between query parts also influence the retrieval results. The higher frequency a query part is specified, the more assemblies containing it are ranked in the front.

Ranking characteristics on different assembly levels

Due to the complexity of assemblies, designers often tend to produce simple sub-assemblies first, and then assemble these subproductions into a more complicated one. All the relationships between assemblies, sub-assemblies and single parts form a treelike model, as shown in Fig. 9. In this situation, the parts participate in the assemblies on different assembly levels. Using our relaxed retrieval framework, the sub-assemblies on different assembly levels can be retrieved and sorted well.

We take the three sub-assemblies in a sedan as an example, as shown in Fig. 9. In this assembly tree, the tyre assembly is composed of parts a and b in Fig. 8; the wheel assembly contains parts c to g in Fig. 8, and one type assembly. As depicted in the top of Fig. 9, the front axle assembly constitutes two wheel assemblies and some other subsidiary parts.

By taking different parts in Fig. 8 as input, the similarities calculated by our algorithm vary accordingly, as listed in Table 2.In the first step, we only input model a, and the tyre, the wheel and the front axle assemblies are returned in descending similarities. Though all three assemblies contain a, the tyre contains the least parts, so it matches the query best. In the second step, we add b as input. Since the tyre is composed of a and b exactly, it matches the queries perfectly, and thus has similarity 1.00. The similarities of wheel and front axle are also higher than that in the first retrieval, since both of them contain a and b. In the third step, we take a, b and c as input. For the reason that only the wheel and the front axle contain c, their similarities rise, and the similarity of the tyre decreases. In the fourth step, the models a to d are employed in the query. Different from the ranking in the first three rounds, the wheel and the front axle are ranked top, because they match the queries better than that of the tyre. The difference between Step five and Step four is that we replace d with g in the query list. This modification does not affect the matchings of the wheel and the front axle, while making the matching of the tyre worse, so its similarity decreases significantly. Furthermore, we have tested our method on silencer, mounting bracket and other assemblies, all of which are composed of sub- 8. The parts that participate in the front axle assembly and its sub-assemblies. We denote them in different colors. The assemblies that these parts participate in are shown in Fig. 9. Fig. 9. The three-level assembly tree that is composed of the parts shown in Fig. 8. In this assembly tree, the assemblies from top to bottom are: (a) the front axle assembly, (b) the wheel assembly and (c) the tyre assembly, respectively. All the assemblies and sub-assemblies are in blue boxes. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) assemblies and form the multi-level assembly trees. By adjusting the queries, the similarities and rankings of the sub-assemblies change accordingly, and the variation tendencies are similar to Table 2.

Sub-assemblies are higher granularity reused units compared to single parts. Using our relaxed retrieval method, the higher reusable sub-assemblies can be easily retrieved by providing only a few of their participating parts. Our method provides an easy way

Table 2

The rankings of the three assemblies by specifying different combinations of parts in Fig. 8. We use indices a to g to represent the queries. The three associated assemblies are listed in columns 2-4. In each row, the assembly names are followed by their corresponding similarity measurements defined in Eq. ( 9). to organize and reuse those assemblies in PDM databases. This is of great help for assembly analysis and modifications.

Queries

Influence of fuzzy queries on retrieval results

To illustrate the variations of retrieval results caused by fuzzy queries intuitively, we take a tyre as an experimental case, as shown in Fig. 10. The three most similar parts measured by Eq. ( 3) are taken as queries in sequence, and their top 5 assemblies ranked by our relaxed retrieval algorithm are presented in rows 2-4, respectively. The simple result indicates that, with a limited modification of similar queries, the retrieval results ranked on the top do not change distinctly.

A statistical experiment was designed, targeting at testing the influence of fuzzy queries. For each of the 1755 typical parts, we queried the top k assemblies, and set them as benchmarks. Then we took the top m similar parts as inputs in sequence, and queried their top k assemblies, respectively. Finally, the assemblies returned by the m similar parts were compared with the benchmark, and the percentages of coincident assemblies were calculated. Fig. 11 illustrates the average percentages of coincident assemblies with the variations of k and m.

We have the conclusion that by setting the most similar parts as input, more than 95% of the assemblies coincide with the benchmark. With the increase of m, the percentages of coincidence assemblies decrease accordingly. However, even if the fifth similar part is set as the query, there are still more than 85% coincident assemblies in the top 20 assemblies. Another observation is that the percentages of coincident assemblies increase quickly when k is smaller than 5, while they grow slowly when k is larger than 6. This indicates that we should pay more attention to all the top 5 assemblies, rather than being confined to the first one when seeking our required results.

Since our relaxed retrieval strategy is not sensitive to similar parts, we can retrieve fuzzy assemblies by inputting fuzzy queries. This is useful if users do not know the exact parts that an assembly contains, or want to retrieve fuzzy assemblies for conceptual design. For instance, if we input a set of parts that participate in the automobiles, the automobiles whose parts are in small variations can also be retrieved.

Performance statistics of the relaxed retrieval

Another statistical experiment was designed to test the effectiveness of our relaxed retrieval algorithm. For each assembly, we selected q parts randomly, and then performed the relaxed retrieval by setting them as input. The percentages that the original assemblies were ranked in the top k were evaluated. In our experiment, k varied from 1 to 20, and q varied from 1 to 10. To reduce the errors introduced by random selection, the experiment was performed 10 times, and the average percentages were recorded. In Fig. 12(a), the percentages are presented when q is set to 1, 3, 5, 7 and 9.

The results show that, with the increase of q, the performance gets better and better, since the more parts we provide in the query, the better the characteristics of the required assemblies can be carried. When 9 (which is close to the average count of parts in an assembly, according to Table 1) randomly selected query parts are provided, about 74% of the required assemblies can be found in the top 20 retrieval results. The statistical results also show that it is not wise to input only one part as query, since the percentage is only about 15% by statistics.

To illustrate the power of important parts in assemblies, another contrast experiment was conducted, and the result is depicted in Fig. 12(b). The only difference is that we selected the top q parts with highest weights in the assembly, rather than using the random selection tactic. The results show that the performance improves about 20% on average compared to the randomly selected queries. This demonstrates that high weight parts play more important roles in distinguishing their associated assemblies from others; and in addition, it verifies that the standard tf * idf weighting scheme characterizes the importance of parts well in assembly retrieval.

Conclusion and future work

In this paper, we have presented a novel framework for lightweight assembly retrieval. The most significant new idea of our approach is that the VSM has been employed in assembly retrieval successfully, and the original VSM in assembly retrieval has been extended into a relaxed matching framework. Since the VSM provides a graceful manner for partial matching, and the relaxed retrieval framework takes both the part similarities and assembly compositions into consideration, the system works pretty well in assembly retrieval, which is demonstrated by a mass of experiments.

To make the assembly retrieval problem suitable for the application of classical IR models, the face normal compatibility based assembly decomposition algorithm and the scale-sensitive similarity measurement based parts merging algorithm have been presented for the preprocessing of assemblies.

Furthermore, the typical parts based matching process, as well as the greedy strategy based matching algorithm for the solution of bipartite graph matching problems is proposed and integrated in the framework. Though both of them are approximate approaches, they preserve the retrieval performance well, while make our system achieve interactive performance.

Compared to Deshmukh's [START_REF] Deshmukh | Content-based assembly search: a step towards assembly reuse[END_REF] and Chen's [START_REF] Chen | A flexible assembly retrieval approach for model reuse[END_REF] methods, our algorithm only utilizes the geometrical information of parts and the composition information of assemblies, and thus can be applied Experiment results showing the effectiveness of the relaxed retrieval algorithm. For each assembly, the relaxed retrieval task is performed by setting its selected q parts as inputs. The percentage that the original assembly is ranked in the top k results is calculated, where k varies from 1 to 20. In the randomly selecting group, the experiment was performed for 10 times, and the results were averaged to reduce the errors introduced by random selection.

to the assemblies represented in lightweight file format in PDM and PLM frameworks. In addition, the extension of VSM enables our algorithm to support fuzzy queries and partial matching. This is especially helpful for the designers if they want to query assemblies by simply providing parts or sub-assemblies. However, there are still some limitations in our current work. For example, while the LFD descriptor is good at capturing the global geometry, it fails to match local details. On this occasion, two similar parts having different functionalities and assembling requirements may be matched false positively, such as the lock nuts and washers. We will seek more powerful 3D CAD part descriptors in our future work, with which the local details of parts, especially the key features for assembling can be characterized. Moreover, by understanding the functionalities of parts, some semantic based retrieval approaches may also be supposed to improve the assembly retrieval performance.

Another limitation is the lack of considering relative positions between query parts. To some extent, relative positions between important parts also provide useful information for the specification of the required assemblies. We shall provide more powerful user interfaces, such that the relative positions between parts can be specified by users and characterized by our system. It is promising that by integrating the position measurement into the matching evaluation, the performance of the relaxed retrieval would be further improved.

3 end 4 while 5 Initialize an new mesh M; 6 Q

 3456 Exists a ''not visited'' face f in A do ←∅;

visited'' and is compatible with u then 12 ENQUEUE 15 Add u to M; 16 Mark 18 AddFig. 2 .

 121516182 Fig. 2. The illustration of the decomposed parts from an assembly. The pipe fitting is composed of 2 flanges and 1 adapter. The black line segments indicate the normals of their corresponding faces. Note that the normals in each watertight mesh are compatible.

1 foreach part u ∈ S do 2 C

 2 ← MAKE-SET(u); 3 end 4 for i ← 1 to size(S) -1 do 5 u ← the ith part in S; 6 for j ← i + 1 to size(S) do 7 v ← the jth part in S; 8 Calculate s(u,v)using Eq. (3); 9 if s(u,v) ≥ ǫ then 10 if FIND-SET(u) = FIND-SET(v) then 11 UNION(u,v); ⊲ Union the two sets.

Fig. 4 .

 4 Fig. 4. The illustration of the relaxed matching framework. w 1q ,w 2q ,...,w mq represent the weights of the m query parts, while w 1p ,w 2p ,...,w np represent the weights of the n assembly parts. The values s 11 , s 12 ,...,s mn represent the scalesensitive similarities between query parts and assembly parts.

7 for j ← 1 9 if S ij > maxMatchValue then 10 maxMatchValue ← S ij ; 11 maxMatchPart 15 S

 719101115 to size(P) do 8 if P[j] is marked as ''not matched'' then ← S + maxMatchValue; 16 Mark the maxMatchPart-th part as matched; 17 end

Fig. 5 .

 5 Fig.5. The user interface of our prototype system. The left side shows the query list, and the right side shows the similar parts of the specified part or the returned assemblies. In this snapshot, the first part is a pipe fitting, and the second one is a flange. The right side shows the list of parts that are similar to the pipe fitting.

Fig. 6 .

 6 Fig.6. The comparison between exact retrieval and relaxed retrieval by providing the same flange. The left side shows the query parts, and the integers below are their frequencies; the retrieved results are depicted in the right side, with similarities listed below the corresponding pictures.

Fig. 7 .

 7 Fig. 7. The relaxed retrieval result of 2 flanges and 1 pipe fitting. The left side is the two inputs, while the right side is the top 10 returned assemblies.

Fig.

  Fig.8. The parts that participate in the front axle assembly and its sub-assemblies. We denote them in different colors. The assemblies that these parts participate in are shown in Fig.9.

Fig. 10 .

 10 Fig.10. The test of the relaxed retrieval with similar parts. On the left side, the first row is the original tyre, and 2-4 are its top 3 similar parts. The top 5 assemblies returned by our relaxed retrieval algorithm are presented on the right side of their corresponding rows.

Fig. 11 .

 11 Fig.11. The percentage of coincident assemblies, according to the similar input variations. The set of assemblies queried by the original part is set as the benchmark. m in the legend stands for the rank of the according input measured by the original part.

  (a) Retrieval results with randomly selected queries.(b) Retrieval results with high weight queries.

Fig. 12 .

 12 Fig. 12.Experiment results showing the effectiveness of the relaxed retrieval algorithm. For each assembly, the relaxed retrieval task is performed by setting its selected q parts as inputs. The percentage that the original assembly is ranked in the top k results is calculated, where k varies from 1 to 20. In the randomly selecting group, the experiment was performed for 10 times, and the results were averaged to reduce the errors introduced by random selection.

Acknowledgements

The research was supported by National Basic Research Program of China (2010CB328001), National Science Foundation of China (61003096), and National High-tech R&D Program (2012AA040902).


Hai-Chuan Song

Jun-Hai Yong

Algorithm for directional projection of parametric curves onto B-spline surfaces

Keywords: Directional projection, B-spline, Curves on surfaces, Approximation 1

This paper proposes an algorithm for calculating the directional projection of parametric curves onto B-spline surfaces. It consists of a second order tracing method with which we construct a polyline to approximate the pre-image curve of the directional projection curve in the parametric domain of the base surface. The final 3D approximate curve is obtained by mapping the approximate polyline onto the base surface. The Hausdorff distance between the exact directional projection curve and the approximate curve is controlled less than the user-specified distance tolerance. And the continuity of the approximate curve is 1 T G F , where T F is the user-specified angle tolerance. Experiments demonstrate that our algorithm is faster than the existing first order algorithm.

INTRODUCTION

Directional projection of curves onto surfaces is widely implemented in CAD (computer aided design) systems [START_REF]CATIA[END_REF], [START_REF]SolidWorks[END_REF]. It plays fundamental roles in many operations, for example surface cutting, model designing, etc. It can also be utilized as a tool for constructing curves on surfaces. Given parametric curve () Pt and surface (,) Suv in 3D space, the directional projection of () Pt onto S can be defined by the directional projection of the points of () Pt onto S . Letting p denote an arbitrary point on () Pt , which is called the test point, the directional projection of p onto S is a set of points q on S such that the vector () qp is parallel and in the same direction with a given vector Dir , which can be expressed as follows:

() 0 , qpD i r q and () 0 qp D i r ¸ holds. As we move the test point p along () Pt , the motion of q results in a set of points on S , and that is the directional projection curve (see Fig. 1.). Denote the directional projection curve by () Qt . To get () Qt , first order algorithm were proposed in [START_REF] Wang | Projecting curves onto free-form surfaces[END_REF]. They derived the first order differential equation systems. By solving the system with numerical methods, they got a sequence of points along () Qt . Then the approximate projection curve could be obtained using the interpolation method on the point sequence.

There are some drawbacks in the first order algorithm above. First, the step length of the points they got along () Qt is not well controlled (just by a constant parametric increment), which results in the uneven distribution of the points, and directly impact on the approximation result. Second, the approximate curve does not lie completely on S , which is not acceptable for many CAD applications such as surface blending and surface-surface intersection [START_REF] Yang | Approximate computation of curves on B-spline surfaces[END_REF]. The approximation precision of the approximate curve cannot be controlled. Moreover, they cannot deal with the "jumping" projection case (this will be introduced in Subsection 3.4).

In order to overcome the drawbacks of the first order algorithm, we provide an algorithm dealing with the directional projection of parametric curves onto B-spline surfaces, which approximates the exact directional projection curve ()

Qt with a piecewise curve on S . And this is similar to our published method [START_REF] Song | Algorithm for orthogonal projection of parametric curves onto B-spline surfaces[END_REF], in which we deal with the orthogonal projection of parametric curves onto Bspline surfaces.

According to the definition of the directional projection, the projection of a curve onto a surface can be regarded as the projections of every point of the curve onto the surface. And this is actually related to the intersections of lines and surfaces, which is also one of the key techniques that we will apply in our algorithm, and of which we will make brief survey below.

Line Surface Intersection

The intersection of line and surface is the fundamental problem in curve surface computation, and is of great use in geometric modeling, robotic techniques, collision detections, etc [START_REF] Srijuntongsiri | A condition number analysis of a line-surface intersection algorithm[END_REF]. The subdivision method proposed by Whitted [START_REF] Rubin | 3-Dimensional Representation For Fast Rendering Of Complex Scenes[END_REF], [START_REF] Whitted | Improved Illumination Model For Shaded Display[END_REF] was firstly utilized to solve this problem, which is very easy and of high robustness. However, when elevating the precision requirement, the efficiency of this method will become low. As a transformation of the subdivision method, Nishita et al. [START_REF] Nishita | Ray tracing trimmed rational surface patches[END_REF] proposed the Bézier clipping method, which recursively divide the surface into regions that may contain the solution, and regions that do not contain the solution.

Generally speaking, the intersection of line and surface can be reduced into the solving of a nonlinear equation system. As a classic method, Newton-Raphson method was widely used in calculating the intersection of line and surface [START_REF] Alain | Chebyshev polynomials for boxing and intersections of parametric curves and surfaces[END_REF], [START_REF] Barth | Efficient ray tracing for Bézier and B-spline surfaces[END_REF], [START_REF] Daniel | On Ray Tracing Parametric Surfaces[END_REF], [START_REF] Michael | Ray Tracing Free-Form B-spline Surfaces[END_REF], [START_REF] Qin | New method for speeding up ray tracing NURBS surfaces[END_REF], [START_REF] Srijuntongsiri | A condition number analysis of a line-surface intersection algorithm[END_REF], [START_REF] Wang | An efficient and stable ray tracing algorithm for parametric surfaces[END_REF], [START_REF] Yang | On Speeding Up Ray Tracing Of B-spline Surfaces[END_REF]. The choice of the initial is very important for Newton-Raphson method to converge reliably. As a result, on one hand, some researches on whether the Newton-Raphson will converge with a given initial value were proposed. Toth [START_REF] Daniel | On Ray Tracing Parametric Surfaces[END_REF] determined the converging interval of the Newton-Raphson method by using interval analyzing techniques. Lischinski and Gonczarowski [START_REF] Daniel | Improved techniques for ray tracing parametric surfaces[END_REF] made use of the consistency of the surface and improved Toth's method [START_REF] Daniel | On Ray Tracing Parametric Surfaces[END_REF]. With the help of Kantorovich detection, Srijuntongsiri and Vavasis [START_REF] Srijuntongsiri | A condition number analysis of a line-surface intersection algorithm[END_REF] can determine whether the convergence order of Newton-Raphson method is second order or not. On the other hand, some researchers considered how to provide a good initial value. Generally, they subdivide the surface recursively until every surface patches are flat enough, and the intersection points of the line and the bounding boxes of the surface patches are chosen as the initial value [START_REF] Alain | Chebyshev polynomials for boxing and intersections of parametric curves and surfaces[END_REF], [START_REF] Barth | Efficient ray tracing for Bézier and B-spline surfaces[END_REF], [START_REF] Michael | Ray Tracing Free-Form B-spline Surfaces[END_REF], [START_REF] Qin | New method for speeding up ray tracing NURBS surfaces[END_REF], [START_REF] Yang | On Speeding Up Ray Tracing Of B-spline Surfaces[END_REF]. Wang et al. [START_REF] Wang | An efficient and stable ray tracing algorithm for parametric surfaces[END_REF] combined the Bézier clipping method and Newton-Raphson method together.

OUTLINE OF OUR ALGORITHM

First, some symbols which will be used all through the paper will be introduced. In 3D

p i
Nu and ()

q j
Nv are the pth-degree and qth-degree B-spline basis functions, respectively (see Fig. 2.).

Fig. 2: Surface of the stomach of Venus: the original B-spline surface, the test curve and the projection direction.

We denote the directional projection curve of () Pt onto S by () Qt . Since () Qt lies on S , there exists a pre-image curve () ( () , () ) T qt ut vt of () Qt in the parametric domain of S and () ( () , () ) Qt Sut vt holds. We suppose that () Qt is at least 1 G continuous and () qt is in Bézier representation.

We provide an algorithm which generates an approximate curve of () Qt on S . First, we construct a polyline (a continuous curve composed of one or more line segments) to approximate () qt in the parametric domain of S , taking the Hausdorff distance and angle tolerances into account. Then the final 3D piecewise approximate curve is obtained by mapping the approximate polyline onto S . The Hausdorff distance between () Qt and the final approximate curve is less than D F . And the continuity of the final approximate curve is 1 T G F , which means that for any pair of adjacent curves in a piecewise curve, the angle between the end derivatives at the mutual point of them is less than T F .

Here D F and T F are the user-specified tolerances. The main algorithm flow is described as follows:

1. Compute the directional projection of the starting point of () Pt onto S with the line surface intersection method [START_REF] Wang | An efficient and stable ray tracing algorithm for parametric surfaces[END_REF]. 2. Based on the projection point, with a second order tracing method, an array of points along () qt is derived. Set the array of points as the vertices, and the initial approximate polyline is obtained in the parametric domain of S . 3. According to the user-specified Hausdorff distance tolerance D F , discrete points are sampled at specified positions on () qt and added into the approximate polyline so as to guarantee that the Hausdorff distance between the final approximate curve and () Qt is less than D F .

4. According to the user-specified angle tolerance T F , split () qt at specified position to make sure that the final approximate curve is at least

1 T G F continuous.
5. After the steps above, the final approximate polyline in the parametric domain of S is obtained. Map the polyline onto S using blossoming techniques [START_REF] Derose | Composing Bézier simplexes[END_REF], [START_REF] Derose | Functional composition algorithm via blossoming[END_REF] and the 3D approximate curve of () Qt is obtained.

The rest of the paper is organized as follows. Section 3 introduces the second order tracing method, and constructs the initial approximate polyline in the parametric domain of S . In Section 4, the approximation precision and the continuity of the final approximate projection curve are controlled, using our published method [START_REF] Song | Algorithm for orthogonal projection of parametric curves onto B-spline surfaces[END_REF]. Section 5 shows the experimental results and Section 6 concludes the paper.

CONSTRUCTING THE INITIAL APPROXIMATE POLYLINE OF THE DIRECTIONAL PROJECTION IN THE PARAMETRIC DOMAIN OF S

First, the initial approximate polyline is to be constructed in the parametric domain of S with our second order tracing method, which can be summarized in the following steps: 1. Compute the first directional projection of the points on () Pt onto S , and take the projection point as the current tracing point. 2. Compute the first and second order derivatives of () qt at the current tracing point in the parametric domain of S and locally approximate () qt with the second order Taylor polynomial. 3. Get the parametric increment using the constant arc length increment strategy. Then the new test point on () Pt and the estimated position of the next tracing point are obtained. 4. Refine the accuracy of the estimation with the Newton-Raphson method. Connect the result point of the refinement and the current tracing point with a line segment, and add it at the end of the approximate polyline. Take this refinement result point as the new tracing point. 5. Repeat steps 2 ~ 4 until the test point we get reaches the end of () Pt and the initial approximate polyline is obtained.

Preparation for Tracing

Suppose 0 P to be the starting point of () Pt . To start the tracing we need to compute the directional projection points of 0 P onto (,) Suv . For these points lies on a line L , which passes 0 P and is parallel with Dir , so we derive these points utilizing line surface intersection method. We first divide the B-spline surface (, ) Suv into several Bézier surface patches (see Fig. 3.). Then we can just take advantage of the division result to compute the directional projection of 0 P . And it can be figured out in the following three steps:

1. First we take the surface patches, whose control polygons intersect with L , as the candidate surface patches. According to the convex hull ability of Bézier surfaces, other surfaces can be pruned. So we can reduce the range of the computation. 2. Then for every candidate surface patch, use the line surface intersection method [START_REF] Wang | An efficient and stable ray tracing algorithm for parametric surfaces[END_REF] to compute the intersection of L with it. And each intersection point is taken as one of the candidate intersection points of L with the base surface. So we can get the intersection points of L and the base surface. Assume that one of the projection points of 0 P onto S is 0 Q , meanwhile in the parametric domain of S , we can get the corresponding pre-image point 0 q on () qt .

To control the tracing process, we preserve two points After the directional projection of 0 P , the point 0 q is taken as the current tracing point _ cur Q and 0 P is taken as the current test point _ cur P . To keep on tracing, we need to get the next tracing point according to _ cur Q . Note that there may exist more than one intersection point of L and the base surface, and for each intersection point, we take it as a seed point, and trace the corresponding projection curve, respectively. Taking one projection curve for example, now we will introduce how to construct the approximate curve of the projection curve.

3.2

The First and Second Order Deriv ativ es of the Directional Projection Curv e

We estimate the next tracing point with the second order Taylor expansion of () qt at _ cur Q , so we need to compute the first and second order derivatives of () qt at _ cur Q . Now we will deduce the first and second order differential equation systems satisfied by the directional projection. Assume that S is at least 2 C continuous locally. First, recall the definition of the directional projection:

() 0 , qpD i r q (3.1)
where p is the point on the test curve () Pt , q is the corresponding point on the projection curve () Qt , Dir is the projection direction vector, (,) Suv is the base surface, and () ( () , () ) Qt Sut vt holds.

Take the derivative of Eqn. (3.1) on both sides with respect to t , and it follows that:

0. uv du dv dp SS D i r dt dt dt ¬ q ®
According to the distributive law of cross product, we have:

() () , uv du dv dp S Dir
S Dir Dir dt dt dt q q q multiply u S and v S on both sides of the equation above, respectively, and it follows that:

• . • vu v uv u du dp SS D i r S D i r dt dt dv dp SSD i r S D i r dt dt £ ¬ ¦ ¦ q q ¦ ¦ ¦ ® ¦ ¤ ¬ ¦ ¦ q q ¦ ¦ ¦ ® ¥ ¦ ¸
Solve the system above, we have:

• . • v vu u uv dp SD i r dt du dt SS D i r dp SD i r dt dv dt SSD i r £ ¬ ¦ ¦ q ¦ ¦ ® ¦ ¦ ¦ ¦ q ¦ ¦ ¤ ¬ ¦ ¦ q ¦ ¦ ® ¦ ¦ q ¥ ¦ ¦ ¦ ¦ ¸ (3.2) Eqn. (3.
2) is just the first order differential equation system of the directional projection, and this was first derived by [START_REF] Wang | Projecting curves onto free-form surfaces[END_REF]. Moreover, take the derivative of Eqn. (3.2) on both sides with respect to t , we have:

2 2 2 2 2 2 2 2 • vv u v v vu u u v v u v v
¬ ¬ ¬ q q ® ® ® ¬ ¬ ¬ q q ® ® ® ¬ q ¸ ® . • uu uv u v uv v u
£ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¤ ¦ ¬ ¬ ¦ ¦ q ¦ ¦ ® ® ¦ ¦ ¦ ¬ ¬ ¬ ¦ ¦ q q ¦ ¦ ® ¸® ¦ ® ¦ ¥ (3.3)
And Eqn. (3.3) is just the second order differential equation system of the direction projection.

Searching For the Next Tracing Point

We estimate the position of the next tracing point with the second order Taylor expansion of () qt at _ cur Q , whose expression is as follows: % , which directly impacts on the distribution of the sampling points and the approximation result.

2 2 0 00 0 () ( ) () () ( ) 2 ,
In order to obtain an evenly distributed sampling during the generation process of the initial approximate polyline, we use the step length controlling method based on constant arc length increment used in [START_REF] Song | Algorithm for orthogonal projection of parametric curves onto B-spline surfaces[END_REF] as follows.

Let ds and dt denote the arc differential and parameter differential of () Qt , respectively. It follows that:

() . ds Q t dt a
Further, according to the chain rule, we can get the expression of () Qt a :

) This step length controlling method is incomplete in [START_REF] Song | Algorithm for orthogonal projection of parametric curves onto B-spline surfaces[END_REF], because it is hard to estimate the total arc length of the projection curve before we can generate it. And improper arc length increments may leads to the wrong results. In the case of directional projection, the arc length of the projection curve () Qt can be estimated like this. We can first sample on the test curve () Pt , link the sample points with a polyline, with which we approximate () Pt . Then we directionally project this polyline onto the control polygons, which are only composed of planar triangles, of the Bézier surface patches we derived in Subsection 3.1. The directional projections of polyline onto triangles can be computed in linear time and are much easier to realize, whose projection result curve is also a polyline. Then the length of the projection result polyline is taken as the approximate arc length of the projection curve () Qt , and we can choose a proper arc length increment according to the estimated total arc length.

Usually, the estimated point deviates from the precise point we desire, so an error exists. According to the definition of the directional projection in Eqn. ( (,) ( ((,) ) ) (,) (,) 0

uu vv fuv Suv p D i r S uv r S uv g u v S u v p D i r Su v rSu v £ ¦ q ¦ ¦ ¤ ¦ q ¦ ¦ ¥ (3.9)
Solve the equation system above with Newton-Raphson method, where we set:

1 1 , ii i ii uu u vv v E ¯ % ¡°¡ ° ¡°¡ °% ¡°¡ °¢±¢ ± ) • , ( () u v uu uv u v i uv u v v u v v ff Sr Sr SD i r S J gg SrSD i r S S r ¯ q ¡°¡ ° ¡°¡ °q ¡°¸¡°¢ ± ¢ ± ¸ (,) . (,) ii i ii fu v k gu v ¡ ° ¡°¡ °¢±
Then the parameter increment of the iteration can be obtained by solving:

. ii i Jk E
Then the convergence criteria of the iteration are as follows:

1. where (,) ii uv is the parameter of the ith iteration, 1 F and 2 F two zero tolerances of Euclidean distance and sine. The iteration is halted if any of the three conditions above is satisfied. In the iterations above, we set the initial parameter value equals to (, ) uv , and set _ pc u rP in the above Newton-Raphson method. The convergence of the iteration depends on the errors of the estimated initial values. As shown in the experimental results in Section 5, the average estimation error of our algorithm is comparable with that of the first order algorithms [START_REF] Wang | Projecting curves onto free-form surfaces[END_REF] implemented with Runge-Kutta method, which have been widely used. The tracing will not stop until the test point we obtain reaches the end of () Pt . Then we can get the initial approximate polyline of () qt in the parametric domain of S .

11 1 () ( , ) () ( , ) , ii u i i ii v i i uu S u vvv S u vF b 2. 1 (,) _ , ii Su v c u r P F b 3.

Projection Jumping Checking

In the generating process of the approximate projection curve, the directional projection points of the test point on () Pt onto S may "jump", with the increment of the parameter (see Fig. L is bigger than s % .

The "jumping" may happen if any of the four conditions above is satisfied. Under this condition, we need to compute the intersection points of 2 L with the base surface (, ) Suv using line surface intersection method [START_REF] Wang | An efficient and stable ray tracing algorithm for parametric surfaces[END_REF] as we did in Subsection 3.1. For each intersection point derived, match it with the estimated point (we got in Subsection 3.3) of current tracing points _ i cur Q (each estimated point can only match to one nearest intersection point within the distance s % ). Each matched intersection point is taken as the new tracing point of the corresponding projection curve. And each unmatched intersection point is taken as a new seed point which is the start point of a new segment of projection curve.

Another special case of "jumping" is the projection curve crosses the boundary of the surface, where the number of intersection points of 2 L and (, ) Suv is less than that of 1 L . This can be detected like this.

1. The estimated point (we got in Subsection 3.3) of current tracing points _ i cur Q exists the boundary of the surface.

The number of intersection points of 2

L and the control polygons of the surface patches is less than that of 1 L .

The "jumping" may happen if any of the two conditions above is satisfied. Under this condition, we need to compute the intersection points of 2 L with the base surface (, ) Suv using line surface intersection method [START_REF] Wang | An efficient and stable ray tracing algorithm for parametric surfaces[END_REF] as we did in Subsection 3.1. If 2 L does not intersect the surface at any projection curve, the tracing of the corresponding projection curve should be stopped. Fig. 5. shows our processing result of this "jumping" case. Recall that the first order algorithm [START_REF] Wang | Projecting curves onto free-form surfaces[END_REF] only solve the first differential equation systems with numerical method, so it cannot deal with this "jumping" case, which frequently appears in practical.

CONTROLLING THE APPROXIMATION PRECISION AND THE CONTINUITY OF THE FINAL APPROXIMATE PROJECTION CURVE

After the construction of the initial approximate polyline in the parametric domain of (,) Suv , we control the approximation precision and the continuity of the final approximate curve with the same method we used in [START_REF] Song | Algorithm for orthogonal projection of parametric curves onto B-spline surfaces[END_REF], which adds more sampling points of () Qt into the approximate polyline and guarantees that the Hausdorff distance between the final 3D approximate curve and the exact projection curve () Qt is less than the user-specified Hausdorff distance tolerance D F , and the continuity of the final 3D approximate curve is

1 T G F , where T
F is the user-specified angle tolerance.

More details can be found in [START_REF] Song | Algorithm for orthogonal projection of parametric curves onto B-spline surfaces[END_REF].

With the final approximate polyline in the parametric domain of S , we map the polyline to the Bézier surfaces we get in Subsection 3.1 using blossoming techniques [START_REF] Derose | Composing Bézier simplexes[END_REF], [START_REF] Derose | Functional composition algorithm via blossoming[END_REF], and the corresponding 3D approximate curve of the exact directional projection curve () Qt on S is obtained, which is in Bézier form and lies completely on S . The degree of the 3D approximate curve is pq , where p and q are the u-direction degree and v-direction degree of the B-spline surface S , respectively.

EXPERIMENTAL RESULTS AND COMPARISONS

We We implement the first order algorithms in [START_REF] Wang | Projecting curves onto free-form surfaces[END_REF] with fourth order Runge-Kutta method (R-K4), which is also the recommended method in [START_REF] Wang | Projecting curves onto free-form surfaces[END_REF]. Comparisons between the first order algorithms [START_REF] Wang | Projecting curves onto free-form surfaces[END_REF] and our second order tracing method on efficiency and accuracy are performed during the process of generating the initial approximate polyline in Example 1, where a cubic B-spline curve () Pt is directionally projected onto a 51 q order Bézier surface, as illustrated in Fig. 6.. Specifically, the 6 control points of () Pt are: And the knot vector is: (0:00; 0:00; 0:00; 0:00; 0:26; 0:68; 1:00; 1:00; 1:00; 1:00): The 61 q control points of the surface are:

£ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¤ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¥
For the purpose of equity, constant parametric increment strategy, as Runge-Kutta method does, are utilized in both of the two kinds of algorithms. Further, according to this paper, both of the two kinds of algorithms generate the initial approximate polyline by estimating the position of the tracing point. In order to measure the accuracy of the two kinds of algorithms, we record the estimation errors, defined in Eqn. (3.8), for each estimated point of the two kinds of algorithms, and compute the average estimation errors for each parametric increment of the two kinds of algorithms. And the refinement based on the Newton-Raphson method introduced in Subsection 3. As shown in Tab. 1., as we decrease the parametric increment, the average errors decrease, the time costs and the numbers of discrete points increase. We first consider the errors. For each parametric increment, the average errors of our algorithm are comparable with those of the first order Computer-Aided Design & Applications, 10(a), 2013, bbb-ccc © 2013 CAD Solutions, LLC, http://www.cadanda.com 13 algorithm with R-K4, and are even smaller than the first order algorithm when the parametric increment decreases to 5 11 0 q . As for the time costs, the difference of the two kinds of algorithms can be clearly illustrated in Fig. 7. Our algorithm is about 2 times faster than the first order algorithms implemented in R-K4 method. And the gaps are getting larger as the parametric increment decreases. Besides the comparisons with the first order algorithm during the process of generating the initial approximate polyline, we also test the performance of our algorithm for the whole projection process, which uses the constant arc length increment strategy, performs "jumping" checking and controls the approximation precision and the continuity of the final approximate projection curve. As we can see from the results, the efficiency of the algorithm meets the real-time requirement. Moreover, the approximation precision and the continuity of the final approximate projection curve can be controlled.

Then we present two examples for "jumping" case projection, which we introduced in Subsection 3.4. Fig. 9(a). and 9(b). plot the "jumping" case projections, where two cubic B-spline curves are directionally projected onto a tour (transformed into NURBS form) and a free form B-spline surface, respectively. The solid blue curves are the test curve, while the curve with green break points is the 3D approximate projection curve. As we can see from the results, the "jumping" is correctly detected, and new segment of projection curves are added into the final projection result, and the projection curve stops at the boundary of the surface.

CONCLUSION

This paper presents an approximation algorithm for directional projection of parametric curves onto B-spline surfaces. The second order tracing method is applied to construct the initial approximate polyline in the parametric domain of the surface. With the initial approximate polyline, further sampling on the exact directional projection curve is performed to control the approximation precision and the continuity of the final 3D approximate curve.

The main contribution of this paper is the derivation of the second order differential equation system of the directional projection, with which we can compute the farthest points and the Hausdorff distances between the curve segments and their corresponding line segments in the parametric domain using iteration methods, which enables us to control the approximation precision and the continuity of the approximate curve [START_REF] Song | Algorithm for orthogonal projection of parametric curves onto B-spline surfaces[END_REF].

Experimental results indicate that the accuracy of our algorithm is comparable with that of the first order algorithm [START_REF] Wang | Projecting curves onto free-form surfaces[END_REF], and at the same time our algorithm is faster than it. Moreover, our algorithm can deal with the "jumping" projection case, and control the approximation precision and the continuity of the approximate curve, while the first order algorithm [START_REF] Wang | Projecting curves onto free-form surfaces[END_REF] cannot do this.

Fig. 1 :

 1 Fig. 1: The definition of the directional projection.

3 .Fig. 3 :

 33 Fig. 3: Surface of the stomach of Venus: (a) the original B-spline surface; (b) the division result and the candidate surface patch.

3 _ cur P R and 2 _

 32 cur Q R to store the values of the current test point on () Pt and the corresponding current tracing point on () qt . _ cur P and _ cur Q have the same parameter value on () Pt and () qt , respectively.

Let 2 _

 2 next Q R denote the result point of the refinement in the parametric domain of S . Connect _ cur Q and _ next Q with a line segment L , and add L into the approximate polyline. Take _ next Q as the new tracing point and we continue seeking for the next tracing point.

2 Q

 2 Pt only intersects with surface in the back part (the point 1 Q in Fig. 4.). However the projection line 2 L of the next tracing point in front and back parts of the surface (the points 2 Q and ' in Fig. 4.). And this is called the "jumping" of the intersection points. If we still only trace the projection curve in the back part of the surface, the should-exist projection curve in the front part of the surface will be missed. We deal with this case as follows.

Fig. 4 :

 4 Fig. 4: The jumping of the intersection points.

, but the distance from any of intersection points of 2 L

 2 and the control polygon of _ Inter Patch to the corresponding intersection point of 1

Fig. 5 :

 5 Fig. 5: The processing result of jumping.

10 FF , where 1 F , 2 F

 1012 present several examples for the directional projection. All the experiments are implemented with Microsoft Visual Studio 2005 with windows XP on the same PC with Intel Core2 Duo CPU 2.53 GHz, 1GB Memory. In all of our experiments below 10 12 are two convergence tolerances for Newton-Raphson method introduced in Subsection 3.3.

Fig. 6 : 1 ,

 61 Fig. 6: Example 1, Directional projection of a cubic B-spline curve onto the bicubic B-spline free form surface where the projection direction is opposite to z axis.

  Design & Applications, 10(a), 2013, bbb-ccc © 2013 CAD Solutions, LLC, http://www.cadanda.com

Fig. 7 :

 7 Fig. 7: Efficiency comparison between our second order tracing method and the first order algorithms for Example 1.

Fig. 8 :

 8 Fig. 8: Directional projection of a cubic B-spline curve onto the bicubic B-spline free form surface where the projection direction is opposite to z axis, 3 11 0 D F q , 10 T F

  Fig. 8. shows the directional projection of cubic B-spline curves onto the bicubic B-spline free form surfaces. The results of our algorithm on Example 2 and Example 3 are shown in Tab. 2..

Fig. 9 :

 9 Fig.9: "jumping" case projection.

Fig. 10 :

 10 Fig. 10: The design of the front cover of a car.

Fig. 11 :

 11 Fig. 11: Construction of curves on surface.

  , (

						Qt a	S	uv du S dt	dv dt	(3.6)
	where the values of u S , v S ,	du dt	,	dv dt	can all be calculated with the given value of 0 t , respectively.
	Hence, with the user-specified arc length increment s % on () Qt , the corresponding approximate
	parametric increment t % is determined by:
						t %	() s Qt a %	.	(3.7)
	Substitute t % into Eqn. (3.5), the estimated position of the next tracing point on () qt can be
	obtained as	0 () ( , ) qt t uv %	.

  Suv is just the projection point of it.As a result, to get the precise point, a refinement is needed. Take the point

	denominator equals to 0. In this situation, we let the error	0 F , because	_ cur P lies on S now and
	(, )
							0 () Pt t % as the new
	test point	_ cur P . We introduce the following Newton-Raphson method based refinement method.
	Multiply u S and v S on both side of Eqn. (3.1), respectively, and set	(( ( , ) rS u v	) p D i r q	)	. We have
	the following equation system:
						(,) ( ((,)	)	)	(,)	(,) 0	,
							3.1), we define the error of the
	estimation as follows:
						F	(, ) (, ) Dir S u v Suv c u r P _ _ cur P q	,	(3.8)
	where	(, ) uv	( qt	0	) % is the estimated point generated by our tracing method in the parametric t
	domain, (, ) Suv is the 3D image point of (, ) uv on S , Dir is the normalized projection direction. The
	geometric significance of F is the absolute value of the sine of the angle between Dir and
	((, ) Suv	_ ) c u r P	. Specially, when (, ) Suv coincides with	_ cur P , Eqn. (3.8) becomes invalid for the

  3 is performed in both of the two kinds of algorithms. The comparison results of the two kinds of algorithms on Example 1 are recorded in Tab. 1..

	Increment	Time (ms)	Errors	Number of discrete points
	t %		R-K4 [16]	ours	R-K4 [16]	ours
	11 0 q	2	1.927	1.486	2.14749e-010	2.1971e-010	101
	11 0 q	3	14.443	9.372	1.02919e-012	2.27651e-010	1001
	11 0 q	4	151.814	96.347	1.71291e-014	9.19197e-014	10001
	11 0 q	5	1698.425	1076.167	1.9516e-014	1.67022e-014	100001
	11 0 q	6	17359.394	7194.017	1.49422e-014	1.42424e-014	1000001

Tab. 1: Time cost and error comparison for Example 1.

© 2013 CAD Solutions, LLC, http://www.cadanda.com

Computer-Aided Design & Applications, 10(a), 2013, bbb-ccc © 2013 CAD Solutions, LLC, http://www.cadanda.com

ACKNOWLEDGEMENTS

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