Rapid wavelet construction of multi-resolution fairing for curves and surfaces with any number of control vertices

In the field of curves and surfaces fairing, arbitrary resolution wavelet fairing algorithm made wavelet fairing technology widely extended to general curves and surfaces, which are determined by any number of control vertices. Unfortunately, accurate wavelet construction algorithm for general curves and surfaces still has not been perfect now. In this article, a concrete algorithm for reconstruction matrix and wavelet construction was emphatically studied, which would be used in the multi-resolution fairing process for curves and surfaces with any number of control vertices. The essence of this algorithm is to generalize wavelet construction into the solution of null space, which could be solved gradually and rapidly by the procedures of decomposition and simplification of coefficient matrix. Certainly, the related compactly supported wavelets could be constructed efficiently and accurately, too. In the last of the article, a complex curve and a complex surface case were provided to verify the stability, high performance, and robustness of this algorithm.


Introduction
In the field of reverse engineering (RE), some kinds of complex parts have high requirements for the fairness of their computer-aided design (CAD) models, such as turbine engine blades, supercharger impellers, compressor rotors, aircrafts, high-speed railway, and so on which must meet the aerodynamics requirements, and household appliances, motor vehicles, and so on which must meet the artistic characteristic requirements. The two primary error sources that affect the fairness of CAD models are measuring accuracy, mainly caused by measuring instruments and measuring methods, and the manufacturing accuracy of the parts themselves.
Digitization of these complex parts mentioned above is an important part of RE. Associated measuring accuracy will greatly determine the final accuracy and fairness of reconstructed curves and surfaces. At the same time, the manufacturing accuracy of these parts will also have a great impact to final measuring accuracy. In order to improve the fairness of reconstructed curves and surfaces, fairness operation has always been a good solution and a research focus in the field of RE.
After years of research, there are many kinds of different ways to fair curves and surfaces, such as energy method, the least square method, the circular rate method, the grid-based spline method, and the rebound method. With the development of fairing algorithm and Wavelets theory, a brand-new fairness method, which is usually called multi-resolution fairing (MRF), appeared and provided a new idea to solve this problem.
The multi-resolution analysis (MRA) is a new mathematical tool that cuts up data, functions, or operators into different frequency components, and then studies each component with a resolution matched to its scale. MRA is usually applied to the field of signal and image processing. 1 If taken curves and surfaces as raw signals, they could be faired certainly by MRA theory. Because of this, a technology called MRF appears. The first paper on this topic is published by Quak and Weyrich 2 in 1994. The core idea is to decompose curves into two parts by reconstruction matrixes P j and Q j . The first one is low-frequency component, and the second one is high-frequency component. Here, low-frequency component is just the faired curves and high-frequency component is the noise exactly that is filtered out. In order to reduce the difficulty of derivation and calculation of the reconstruction matrixes, the curves are specified as cubic B-spline ones and their number of control vertices must be 2 j + 3. Based on the levels of detail (LOD) representation method, Ceruti et al. 3 represented a curve by the least square coefficient and the detail coefficient. So, the detail part was extracted by constructed LOD filter that could determine the geometric defects of the curve, resulting in a faired Class-A curve. Li et al. 4 proposed an improved curve fitting method for the resolution of overlapping peaks. By continuous wavelet transform (WT) of sharp peaks, the reasonable initial values of each peak parameter were obtained to guarantee the curve accuracy on improved fitting conditions. By studying on traditional WT, Hussein et al. 5 proposed a second-generation wavelet fairing technology based on lifting scheme, which mainly composed of splitting, prediction, and updating. Taking quasi-uniform cubic B-spline curve as example, the second-generation wavelet representation method on B-spline curves was deduced successfully. Then, he generalized the lifting wavelet scheme to the fairing operation for irregular complex free-form surfaces.
Because the scaling translation series of scaling functions used in the above literatures are dyadic, the number of control vertices of curves and surfaces must be strictly limited to 2 j + (r2 1) (r is order). Based on this, reconstruction matrixes P j and Q j can be figured out in advance by dyadic wavelet. This kind of algorithm is called dyadic wavelet fairing algorithm (DWFA). But in practice, the number of control vertices of curves and surfaces is arbitrary usually, DWFA mentioned above is not applicable any more. In short, the application of DWFA is theoretically limited in many cases.
In order to extend the applied range of wavelet fairing algorithm to more common curves and surfaces with any number of control vertices, a lot of scholars paid close attention to it. Wang and Zhang 6 constructed the orthogonal non-uniform B-spline wavelets by removing nodes from knot vector according to the given weighting formula. The idea of this algorithm is to change the knot vector to meet the requirements of literature. 2 So as to realize the wavelet fairing and simplification for non-uniform rational B-splines (NURBS) curves and surfaces by DWFA. Based on the discrete norm, Pan and Yao 7 and Sadeghi and Samavati 8 constructed biorthogonal non-uniform Bspline wavelets using least-squares fitting algorithm according to the given knot vector and then achieved the MRF on general free-form curves. Because these algorithms avoid inner product operation, the calculation efficiency has been improved greatly. By constructing the length-preserving projection of curves, Wang et al. 9 deduced the double-scale subdivision equations of multi-scale functions and multi-scale wavelets, and constructed the multi-wavelets on some smooth plane curves. In order to control the better deformation of the curves and obtain the more accurate modeling, Wang and Zhao 10 initially introduced a new concept of reference index of modeling (RIOM) based on wavelet theory. Compared to traditional algorithms, this one takes RIOM as objective function and objective shape, which can quantitatively evaluate objective function value determined by RIOM and ensure the overall shape of the curve at the same time. Fortes and Moncayo 11 and Fortes and Rodrı´guez 12 overcame the limitation that the control vertices of the curves must be uniformly distributed, and proposed a new MRA method for non-uniform surfaces based on the nonuniform triangulation.
While all these algorithms mentioned above contain approximate calculation, faired curves and surfaces could not be reconstructed accurately. This is not consistent with the essence of MRA. The reason for this is all these algorithms put emphasis on the curves and surfaces, not wavelet construction itself.
Arbitrary scale factors MRA provides a new solution from the view of MRA itself. This is another study field of MRA. Lazar and Bruton 13 discovered the relationship between orthogonal wavelet decomposition and lossless perfect reconstruction (PR). Based on the M-band PR filter, a regular and orthogonal wavelet method for generating arbitrary integer scale factor M was proposed. On the premise that the sum of all the branch sampling factors equals to one, Kovacevic and Vetterli 14 constructed a PR filter by setting a rational sampling factor, which realized the non-uniform splitting. On this basis, Blu 15 presented a design algorithm of orthonormal rational filter banks which allow rational wavelet to transform. Munoz et al. 16 proposed a method-based B-spline to allow continuous WT at any scale, and the complexity of the method is not limited by dyadic or integer scales.
In fact, arbitrary resolution wavelet decomposition and reconstruction of general curves and surfaces is a typical MRA at rational scaling factor, and many scholars try to apply the research results of MRA to fair the general curves and surfaces. Kazinnik and Elber 17 proposed a multi-resolution orthogonal decomposition method for non-uniform B-spline (non-uniform rational (NUR)) spaces. The decomposition of discontinuous surfaces can preserve the global shape of the surface as well as the local characteristics of the surface. An interactive tool was developed that allows editing of NUR surfaces, which can reflect the effects of different strength of the multi-resolution decomposition. Based on the available non-uniform semi-orthogonal B-spline wavelet algorithm, Li et al. 18 and Li and Tian 19 proposed a new curve multi-resolution representation method without the limitation of wavelet multiresolution representation. This algorithm can be applied to the MRA for any uniform or non-uniform B-spline curves, and reduces the computational cost of construction of non-uniform B-spline wavelets in the process of MRF. Hamm et al. 20 presented a multiresolution representation method to spline curves with arbitrary knots for non-equally distributed geometric data. By modifying knot vector of NURBS surface through newly proposed T-spline control meshes, Wang and Zhao 21 realized the multi-resolution representation of NURBS surface. Generally, the fairing algorithm that has no special requirement on the number of vertices of curves and surfaces is usually called arbitrary resolution wavelet fairing algorithm (ARWFA).
To ARWFA, the author has achieved the concrete MRF algorithm with any rational scaling factor by direct mathematical derivation, according to the definition of wavelet in literature. 22,23 However, the efficiency, quality, stability of wavelet construction are still unsatisfactory and need a further discussion. In this article, we mainly discuss the solution of wavelet construction matrix based on null space.

Basis of MRF
In this article, quasi-uniform cubic B-spline curves were taken as the research object, so the order of curves is r = 4. On the basis of Cox-de Boor recurrence formula, the simplified expression of cubic B-spline basis function can be represented as u(t) = N r (t) = N 4 (t). The knot vector corresponding to the scaling factor m can be expressed as [0, 0, 0, 0, 1/m, 2/m, ..., 1 2 (1/m), 1, 1, 1, 1]. A set of B-spline basis functions defined by knot vector can be made up of their scaling translation series; that is, u(mt2k) = N 4 (mt2k), where k is translation. Set u m k = u(mt À k), it follows that a linear space . . , c m m + 2 ) T can be used to uniquely determine one curve.
Similarly, for the objective scale factor n, which is small than m (n \ m), there are also a set of B-spline basis functions u n l = u(nt À l), where l is translation. It follows that a linear space V n can be also spanned by F n = (u n 0 , u n 1 , . . . , u n n + 2 ), expressed as V n = spanfu n l g. It is found that there are n + 3 B-splines in space V n , so n + 3 control vertices C n = (c n 0 , c n 1 , . . . , c n n + 2 ) T can be used to uniquely determine one curve. Assume that f m is a curve determined by C m in space V m , and f n is a curve determined by C n in space V n .
The essence of MRF based on MRA is to filter out the high-frequency detail curves g n from f m to get a lowfrequency smooth curve f n by the constructed low pass filter banks. That is The tower structure diagram of decomposition and reconstruction of MRF is shown in Figure 1. Here, C m is the control vertices of original curve, C n is the control vertices of faired curve, D n is the control vertices of detail curve, A mn and B mn are the decomposition matrixes, P mn and Q mn are reconstruction matrixes.
For D n = (d n 0 , d n 1 , . . . , d n mÀnÀ1 ) T , its dimension is m2n and it is used to determine the detail curve g n . As a matter of fact, the corresponding basis functions c n j of g n are very the wavelets which would be constructed in this article. These basis functions can be expressed as C n = (c n 0 , c n 1 , . . . , c n mÀnÀ1 ), and can be also spanned into a liner space W n , expressed as W n = spanfc n j g. Usually, W n is commonly referred to as the wavelet subspace of V m about V n . In order to make wavelets to attenuate to zero as quickly as possible, V n 'W n is prerequisite when wavelets are constructed. Then, we have In equation (2), 4 means direct sum. Because any two translational B-spline basis functions at the same scale are not orthogonal, only the orthogonality between F n and C n at the same scale can be guaranteed in the process of wavelet construction, that is

Construction of multi-resolution wavelet
According to equation (1) and the description in the second section, the MRF process of B-spline curve can be expressed as According to the author's literature 22 where A mn is a low-frequency decomposition matrix whose dimension is (n + 3) 3 (m + 3), ½hF n jF m i is a (n + 3) 3 (m + 3) matrix which is calculated by inner product hu n i , u m j i, and ½hF n jF n i is a (n + 3) 3 (n + 3) matrix which is calculated by inner product hu n i , u n j i. B mn is a high-frequency decomposition matrix whose dimension is (m2n) 3 (m + 3), ½hC n jF m i is a (m2n) 3 (n + 3) matrix which is calculated by inner product hc n i , u m j i, and ½hC n jC n i is a (m2n) 3 (m2n) matrix which is calculated by inner product hc n i , c n j i. In equations (5) and (6), F m and F n are already known for given scale m and n. If C n is also known, then wavelet decomposition and reconstruction at any scale could be worked out as a matter of course.

Description of wavelet construction
On account of W n & V m , each wavelet in wavelet space W n is a linear representation of B-spline basis functions in linear space V n , so we have Taking Q mn as a (m + 3) 3 (m2n) matrix determined by elements q k,l , then equation (7) could be expressed as equation (8) in matrix form If a reasonable matrix Q mn could be figured out, wavelets C n then could be constructed by equation (8). As a result, wavelet construction is transformed into the solution of reconstruction matrix Q mn .
According to equation (8), it can be pre-multiplied by F n on both sides to construct an inner product matrix ½hF n jC n i. Because of their orthogonality between F n and C n , equation (8) could be re-expressed as equation (9) hF n jC n i ½ = hF n jF m i ½ Furthermore, because B-spline functions F m and F n are already known, if scale factors m and n are determined, then a (n + 3) 3 (m + 3) inner product matrix ½hF n jF m i could be calculated. Assume that E = ½hF n jF m i and X = Q mn , then equation (9) is transformed into the solution of null space, as expressed in equation (10). Here, E is the coefficient matrix of null space Solution of null space At the beginning of solution, the characteristics of coefficient matrix E need to be analyzed detailedly.
For example, if a quasi-uniform cubic B-spline curve f m with 13 control vertices want to be faired to f n with 8 control vertices, then we have m = 13 2 3 = 10, n = 8 2 3 = 5. So, a perfect corresponding coefficient matrix E 8 3 13 could be achieved in Table 1.
The elements in Table 1 show that 1. Because B-spline basis functions are compactly supported, and F m and F n are scaling translation series of themselves, there are massive zero elements in the bottom left corner and top right corner of coefficient matrix E in the process of inner product calculation. At the same time, non-zero elements focus near the diagonal of the matrix E. 2. According to the wavelet fairing theory, n should be less than m and all non-zero elements focus near the diagonal, so the rank of matrix E is Rank(E) = n + 3.  3. Because n \ m, the coefficient matrix E is flat rectangular. 4. Coefficient matrix E is strictly diagonally dominant. 5. Coefficient matrix E is symmetrical. If the element in ith row and jth column of matrix E could be expressed as E(i, j), then E((n + 4) 2i, (m + 4) 2j) is equal to E(i, j). In other words, the inversion of ith column in Table 1 is just the ((n + 1) 2i)th column.
Because Rank(E) = (n + 3) \ (m + 3), the common solution X of null space (Equation 10) is non-unique and is not strictly diagonally dominant according to linear algebra. Table 2 is one set of common solutions X of above null space.
For scale factor m = 10 and n = 5, m2n = 5 wavelets could be constructed by equation (8), as shown in Figure 2. Obviously, all these wavelets are not compactly supported and do not conform to the characteristics of wavelets.
In order to construct more reasonable wavelets, the solution X of null space must meet some specific requirements.
Considering that wavelet functions are compactly supported, the solution X of null space should also be a sparse and diagonally dominant matrix. That is there should be as more zero elements as possible both in the bottom left corner and top right corner of X. In order to take the solution quality under control, each column of null space should be worked out, respectively. The concrete calculation idea is as follows: 1. Because the solution X of null space based on general numerical computation method is usually orthogonal, there are almost no zero elements in X generally, as shown in Table 2. So, the solution X is not suitable for constructing the compactly supported wavelets. 2. Because of the characteristic of coefficient matrix E, that is, the last column is same as the inversion of the first one, the last column but one is same as the inversion of the second one, and so on, so the ((n + 1) 2i)th column of the solution X is same as the inversion of the ith one. That is, if a series of solutions of EX = 0 are figured out, then the inversion of them are also the solutions of EX = 0.

For an equation set E (n + 3) 3 (m + 3)
X (m + 3) 3 (m2n) = 0 (n + 3) 3 (m2n) , there are n + 3 equations and m + 3 unknown variables; Table 2. One of the solutions of null space X.  therefore, there are at least m2n (n \ m) free unknown variables. In order to make the solution vector contain as more zero elements as possible and be able to construct a band matrix, the good idea is to search the fundamental system of solution of coefficient matrix E. The first step is to take the first n + 1 columns of E to constitute an equation set and obtain a fundamental system of solution. Then, in the second step, take columns from the second column to the (n + 2)th column to constitute another equation set and obtain another fundamental system of solution, and so on.
Expand the equation EX = 0, then it could be described as equation (11) . . .
So, for q i , each ith column vector of X, an equation set could be correspondingly extracted from equation (11). At the same time, initialize the specific elements of column vector q i , which is from the first one to the (i2 1)th one and from the (i + n + 4)th one to the (m + 3)th one, to be zero.
So, the solution procedure of the equation EX = 0 is resolved into the solutions of m2n equation sets E i X i = 0, as shown in equation (12 Set the rank of matrix E i is Rank(E i ) and the column number of E i is Column. If Rank(E i ) = Column2 1, the non-zero solution of X i is unique. If Rank(E i ) \ Column2 1, the non-zero solution of X i is not unique, then the matrix E i must be further decomposed and simplified.
For the situation of Rank(E i ) \ Column2 1, research shows that there are several zero rows (all elements of corresponding rows are zero) at the top and bottom parts of coefficient matrix E i sometimes, as shown in equation (13). So, we can eliminate all zero rows and extract a new coefficient matrix E ij Extract Rank(E ij ) + 1 columns from the left part of E ij to constitute a new coefficient matrix E j . Renewedly, set the rank of matrix E j is Rank(E j ) and the column number of E j is Column. If Rank(E j ) = Column 2 1, the non-zero solution of X j is unique. If Rank(E j ) \ Column2 1, the non-zero solution of X j is not unique, then E j must be further decomposed and simplified in the same way.
In conclusion, the extraction process of coefficient matrix can be represented by the frame diagram as shown in Figure 3 and the solution procedure of null space could be designed and shown in Figure 4.

Verification of algorithm
Now taking the coefficient matrix in Table 1 as an example, on the basis of the above algorithm, the null space and the corresponding construction matrix Q mn are calculated. According to Table 1, there are m = 10 and n = 5, so m2n = 5 equation sets could be formulated and 5 compactly supported wavelets could be constructed.
Second, these five equation sets are, respectively, solved and five untrivial solutions q 1 , q 2 , q 3 , q 4 , q 5 are obtained.
Finally, one of the perfect matrix X = [q 1 , q 2 , q 3 , q 4 , q 5 ] has been constructed, as shown in Table 3.
Obviously, the resulting matrix X is a sparse and diagonally dominant matrix and will meet the requirements of wavelet construction faultlessly.
Standardize Table 3 further to make the maximal elements of every column to be 1, as shown in Table 4. Table 4 is just the finally expected construction matrix Q mn .
According to equation (8), the final reasonable wavelets was constructed, as shown in Figure 5. Figure 5(a) presents m + 3 = 13 B-splines corresponding to scale factor m = 10, Figure 5(b) presents n + 3 = 8 B-splines corresponding to scale factor n = 5, and Figure 5(c) presents m2n = 5 B-spline wavelets, which are constructed in the process of fairing from scale factor m to n. Figure 5 indicates that 1. The constructed wavelets are compactly supported. They can attenuate rapidly to zero after finite oscillations. The reason is that the elements in the bottom left corner and top right corner of construction matrix Q mn acquired by algorithm in this article are all zero. 2. The constructed wavelets are symmetrical, except for the finite wavelets at both ends. The reason of asymmetry of wavelets at the both ends is that the each end point of knot vector of quasi-uniform cubic B-spline curve is quadruple. 3. The constructed wavelets in Figure 5(c) indicate that the construction matrix Q mn obtained by the above algorithm is correct, reasonable, and high quality.
Arbitrary resolution wavelet fairing for general curves and surfaces

ARWFA for curve
In order to verify the correctness, accuracy and stability of the algorithm proposed in this article, a complex curve with 120 control vertices, which is no longer conform to the DWFA, was taken to be faired repeatedly and its detail curves were extracted synchronously. The faired curves and corresponding detail curves are shown in Figure 6.       Figure 6(b) is a directly faired curve determined by only 30 control vertices from 120 control vertices in a single fairing process in 0.636 s. Now, the objective scale factor is n = 30 2 3 = 27 and m2n = 90 wavelets were constructed. Figure 6(c) is the detail curve, which is decomposed from the original curve in the process of wavelet fairing. In order to show the detail curve clearly, Figure 6(c) is magnified 16 times. Figure 6(d) is a faired curve determined by only 20 control vertices, which is faired directly from the original curve in Figure 6(a) in a single process in 0.769 s. Now, the objective scale factor is n = 20 2 3 = 17 and m2n = 100 wavelets were constructed. Figure 6(e) is the detail curve, which is decomposed from the original curve in the process of wavelet fairing. In order to show the detail curve clearly, Figure 6(e) is magnified six times. Figure 6(f) is a faired curve determined by only 10 control vertices, which is faired directly from the    original curve in Figure 6(a) in a single process in 0.943 s. Now, the objective scale factor is n=10 23=7 and m2n=110 wavelets were constructed. Figure 6(g) is the detail curve, which is decomposed from the original curve in the process of wavelet fairing. In order to show the detail curve clearly, Figure 6(g) is magnified two times.
The computational efficiency related to different objective scale in this case was researched in Figure 7. The smaller the objective scale is, the more wavelets should be constructed, and the more time it takes. Our surveys indicate that if the fairing process must meet the requirement of man-machine interaction on the computer, fairing time should be limited in 1 s (control line in Figure 7). Fortunately, all operation times shown in Figure 7 are limited to 1 s usually, so the computational efficiency is acceptable and satisfactory.

ARWFA for surface
Now, another surface fairing case will be provided in this section. Figure 8(a) is a supercharger impeller, and Figure 8(b) is one of the blade surface of Figure 8(a) that will be faired.
The corresponding original surface about Figure  8(b), as shown in Figure 9(a), has 4000 (50 3 80) control vertices, where there are 50 control vertices in the u direction and 80 control vertices in the v direction, that is, m 1 = 47, m 2 = 77.
According to surface fairing framework based on tensor product method in Figure 10, one low-resolution surface and three high-resolution surfaces could be figured out in the process of ARWFA.
Set target scale n 1 = 15, n 2 = 23, after ARWFA, a faired low-resolution surface determined by 468 (18 3 26) control vertices could be figured out by one step, as shown in Figure 9(b), where there are 18 control vertices in the u direction and 26 control vertices in the v direction.
According to surface faired framework in Figure 10, three high-resolution surfaces could be also figured out at the same time, where D 1 15 3 23 , mentioned in Figure  10, is one of the detail surface in u direction, as shown in Figure 11(a). D 2 15 3 23 is one of the detail surface in v direction, as shown in Figure 11(b). D 3 15 3 23 is one of the detail surface in uv direction, as shown in Figure  11(c).
What needs to be emphasized is that original surface C m1 3 m2 could be accurately reconstructed by faired surface C n1 3 n2 and three detail surfaces D 1 n1 3 n2 , D 2 n1 3 n2 and D 3 n1 3 n2 .

Conclusion
For the faultiness of current reconstruction algorithm of arbitrary resolution wavelet fairing for curves and surfaces, a specific algorithm for matrix construction in the process of MRF to curves and surfaces with any number of control vertices was focused on and achieved successfully in this article. The article has carried on the research mainly from the following several aspects: 1. The essence of this algorithm is to generalize wavelet construction into the solution of null space. By means of decomposition and simplification of coefficient matrix E, an expected highquality reconstruction matrix Q mn could be solved directly without modification of knot vector, different from literature. [18][19][20][21] The advantage is that there is no approximate calculation in the process of MRF. Certainly, the related compactly supported wavelets C n could be constructed efficiently and accurately, too. 2. The algorithm could be directly applied to the MRA and wavelet fairing for curves and surfaces with any number of control vertices, and truly realized the PR by the ARWFA, which reflects the essence of accurate reconstruction of MRA. 3. A professional calculation program based on this algorithm has been programmed and tested.
The experiment results show that this algorithm has a high executing efficiency, computational stability, and good practical value.

Declaration of conflicting interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.