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

Abstract

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. ¹ 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 ² 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. ³ 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. ⁴ 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. ⁵ 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 + (r− 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 ⁶ 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. ² 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 ⁷ and Sadeghi and Samavati ⁸ constructed biorthogonal non-uniform B-spline 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. ⁹ 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 ¹⁰ 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 ¹¹ and Fortes and Rodríguez ¹² 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 non-uniform 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 ¹³ 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 ¹⁴ constructed a PR filter by setting a rational sampling factor, which realized the non-uniform splitting. On this basis, Blu ¹⁵ presented a design algorithm of orthonormal rational filter banks which allow rational wavelet to transform. Munoz et al. ¹⁶ 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 ¹⁷ 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. ¹⁸ and Li and Tian ¹⁹ proposed a new curve multi-resolution representation method without the limitation of wavelet multi-resolution 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. ²⁰ presented a multi-resolution 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 ²¹ 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 φ(t) = N_r(t) = N₄(t). The knot vector corresponding to the scaling factor m can be expressed as [0, 0, 0, 0, 1/m, 2/m, …, 1 − (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, φ(mt−k) = N₄(mt−k), where k is translation. Set ${\phi }_k^m=\phi (\mathrm{mt}-k)$ , it follows that a linear space V _m can be spanned by ${\mathit{\Phi }}_m=({\phi }_0^m,{\phi }_1^m,\dots ,{\phi }_{m+2}^m)$ , expressed as ${\mathit{V}}_m=\mathrm{span}\{{\phi }_k^m\}$ . It is found that there are m+3 B-splines in space V _m , so m+3 control vertices ${\mathit{C}}_m=(c_0^m,c_1^m,\dots ,c_{m+2}^m{)}^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 ${\phi }_l^n=\phi (\mathrm{nt}-l)$ , where l is translation. It follows that a linear space V _n can be also spanned by ${\mathit{\Phi }}_n=({\phi }_0^n,{\phi }_1^n,\dots ,{\phi }_{n+2}^n)$ , expressed as ${\mathit{V}}_n=\mathrm{span}\{{\phi }_l^n\}$ . It is found that there are n+3 B-splines in space V _n , so n+3 control vertices ${\mathit{C}}_n=(c_0^n,c_1^n,\dots ,c_{n+2}^n{)}^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 low-frequency smooth curve f_n by the constructed low pass filter banks. That is

$$f_m=f_n+g_n$$

(1)

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.

Figure 1.

Diagram of multi-resolution wavelet fairing.

For ${\mathit{D}}_n=(d_0^n,d_1^n,\dots ,d_{m-n-1}^n{)}^T$ , its dimension is m−n and it is used to determine the detail curve g_n. As a matter of fact, the corresponding basis functions ${\psi }_j^n$ of g_n are very the wavelets which would be constructed in this article. These basis functions can be expressed as ${\mathit{\Psi }}_n=({\psi }_0^n,{\psi }_1^n,\dots ,{\psi }_{m-n-1}^n)$ , and can be also spanned into a liner space W _n , expressed as ${\mathit{W}}_n=\mathrm{span}\{{\psi }_j^n\}$ . 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

$${\mathit{V}}_m={\mathit{V}}_n\oplus {\mathit{W}}_n$$

(2)

In equation (2), ⊕ means direct sum. Because any two translational B-spline basis functions at the same scale are not orthogonal, only the orthogonality between Φ _n and Ψ _n at the same scale can be guaranteed in the process of wavelet construction, that is

$$〈{\phi }_i^n,{\psi }_j^n〉=0,i\in [0,n+2];j\in [0,m-n-1]$$

(3)

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

$$\sum _{k=0}^{m+2}{c}_k^m{\varphi }_k^m=\sum _{l=0}^{n+2}{c}_l^n{\varphi }_l^n+\sum _{l=0}^{m-n-1}{d}_l^n{\psi }_l^n$$

(4)

According to the author’s literature ²²

$${\mathit{A}}_{\mathrm{mn}}={[〈{\mathit{\Phi }}_n|{\mathit{\Phi }}_n〉]}^{-1}[〈{\mathit{\Phi }}_n|{\mathit{\Phi }}_m〉]$$

(5)

$${\mathit{B}}_{\mathrm{mn}}={[〈{\mathit{\Psi }}_n|{\mathit{\Psi }}_n〉]}^{-1}[〈{\mathit{\Psi }}_n|{\mathit{\Phi }}_m〉]$$

(6)

where A _mn is a low-frequency decomposition matrix whose dimension is (n+3) × (m+3), $[〈{\mathit{\Phi }}_n|{\mathit{\Phi }}_m〉]$ is a (n+3) × (m+3) matrix which is calculated by inner product $〈{\phi }_i^n,{\phi }_j^m〉$ , and $[〈{\mathit{\Phi }}_n|{\mathit{\Phi }}_n〉]$ is a (n+3) × (n+3) matrix which is calculated by inner product $〈{\phi }_i^n,{\phi }_j^n〉$ .

B _mn is a high-frequency decomposition matrix whose dimension is (m−n) × (m+3), $[〈{\mathit{\Psi }}_n|{\mathit{\Phi }}_m〉]$ is a (m−n) × (n+3) matrix which is calculated by inner product $〈{\psi }_i^n,{\phi }_j^m〉$ , and $[〈{\mathit{\Psi }}_n|{\mathit{\Psi }}_n〉]$ is a (m−n) × (m−n) matrix which is calculated by inner product $〈{\psi }_i^n,{\psi }_j^n〉$ .

In equations (5) and (6), Φ _m and Φ _n are already known for given scale m and n. If Ψ _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

$${\psi }_l^n=\sum _{k=0}^{m+2}{q}_{k,l}^m{\varphi }_k^m$$

(7)

Taking Q _mn as a (m+3) × (m−n) matrix determined by elements q_k,l, then equation (7) could be expressed as equation (8) in matrix form

$${\mathit{\Psi }}_n={\mathit{\Phi }}_m{\mathit{Q}}_{\mathrm{mn}}$$

(8)

If a reasonable matrix Q _mn could be figured out, wavelets Ψ _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 Φ _n on both sides to construct an inner product matrix $[〈{\mathit{\Phi }}_n|{\mathit{\Psi }}_n〉]$ . Because of their orthogonality between Φ _n and Ψ _n , equation (8) could be re-expressed as equation (9)

$$[〈{\mathit{\Phi }}_n|{\mathit{\Psi }}_n〉]=[〈{\mathit{\Phi }}_n|{\mathit{\Phi }}_m〉]{\mathit{Q}}_{\mathrm{mn}}=0$$

(9)

Furthermore, because B-spline functions Φ _m and Φ _n are already known, if scale factors m and n are determined, then a (n+3) × (m+3) inner product matrix $[〈{\mathit{\Phi }}_n|{\mathit{\Phi }}_m〉]$ could be calculated. Assume that $\mathit{E}=[〈{\mathit{\Phi }}_n|{\mathit{\Phi }}_m〉]$ 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

$$\mathit{EX}=\mathit{0}$$

(10)

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 − 3 = 10, n = 8 − 3 = 5. So, a perfect corresponding coefficient matrix E _8×13 could be achieved in Table 1.

Table 1.

Coefficient matrix E_8×13.

0.018661	0.019821	0.009658	0.001845	1.49E-05	0	0	0	0	0	0	0	0
0.005781	0.023438	0.036533	0.02756	0.006224	0.000461	3.72E-06	0	0	0	0	0	0
0.000543	0.006295	0.024799	0.050685	0.043761	0.019449	0.004158	0.000308	2.48E-06	0	0	0	0
1.49E-05	0.000446	0.004008	0.019603	0.045838	0.060179	0.045838	0.019603	0.004159	0.000308	2.48E-06	0	0
0	0	2.48E-06	0.000308	0.004159	0.019603	0.045838	0.060179	0.045838	0.019603	0.004008	0.000446	1.49E-05
0	0	0	0	2.48E-06	0.000308	0.004158	0.019449	0.043761	0.050685	0.024799	0.006295	0.000543
0	0	0	0	0	0	3.72E-06	0.000461	0.006224	0.02756	0.036533	0.023438	0.005781
0	0	0	0	0	0	0	0	1.49E-05	0.001845	0.009658	0.019821	0.018661

The elements in Table 1 show that

1. Because B-spline basis functions are compactly supported, and Φ _m and Φ _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) −i, (m+ 4) −j) is equal to E (i, j). In other words, the inversion of ith column in Table 1 is just the ((n+ 1) −i)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.

Table 2.

One of the solutions of null space X.

0.190971	−0.142955	−0.109745	0.307595	−0.427767	0.312066	0.079211	−0.455296	0.501944	−0.280361	0.044982	0.063030	−0.056984
−0.042320	0.154388	−0.286404	0.336401	−0.384424	0.275330	−0.013186	−0.014713	−0.353398	0.568249	−0.304127	−0.020774	0.098023
−0.186505	0.380583	−0.457885	0.310808	−0.069431	−0.072300	−0.130414	0.301190	−0.019149	−0.365775	0.418777	−0.256564	0.107181
−0.165787	0.306760	−0.320821	0.145864	0.154940	−0.398121	0.312519	−0.120691	0.108290	−0.023766	−0.278813	0.457552	−0.362345
0.121919	−0.074798	−0.074155	0.096828	0.084995	−0.388223	0.523423	−0.315416	−0.047995	0.140437	0.123988	−0.374187	0.438245

For scale factor m = 10 and n = 5, m−n = 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.

Figure 2.

Five unreasonable wavelets constructed based on Table 2: (a) 1st wavelet, (b) 2nd wavelet, (c) 3rd wavelet, (d) 4th wavelet and (e) 5th wavelet.

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) −i)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.

3. For an equation set E _{(n+3) × (m+3)} X _{(m+3) × (m−n)} = 0_{(n+3) × (m−n)}, there are n+3 equations and m+3 unknown variables; therefore, there are at least m−n (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)

$$\begin{array}{c}\left[\begin{array}{cccc}{e}_{1,1}\hfill & e_{1,2}\hfill & \cdots \hfill & e_{1,m+3}\hfill \\ \multicolumn{1}{c}{e_{2,1}}& e_{2,2}\hfill & \cdots \hfill & e_{2,m+3}\hfill \\ \multicolumn{1}{c}{⋮}& ⋮\hfill & \ddots \hfill & ⋮\hfill \\ \multicolumn{1}{c}{e_{n+3,1}}& e_{n+3,2}\hfill & \cdots \hfill & e_{n+3,m+3}\hfill \end{array}\right]\hfill \\ \left[\begin{array}{cccc}{q}_{1,1}\hfill & q_{1,2}\hfill & \cdots \hfill & q_{1,m-n}\hfill \\ \multicolumn{1}{c}{q_{2,1}}& q_{2,2}\hfill & \cdots \hfill & q_{2,m-n}\hfill \\ \multicolumn{1}{c}{⋮}& ⋮\hfill & \ddots \hfill & ⋮\hfill \\ \multicolumn{1}{c}{q_{m+3,1}}& q_{m+3,2}\hfill & \cdots \hfill & q_{m+3,m-n}\hfill \end{array}\right]\hfill \\ =\left[\begin{array}{cccc}{0}_{1,1}\hfill & 0_{1,2}\hfill & \cdots \hfill & 0_{1,m-n}\hfill \\ \multicolumn{1}{c}{0_{2,1}}& 0_{2,2}\hfill & \cdots \hfill & 0_{2,m-n}\hfill \\ \multicolumn{1}{c}{⋮}& ⋮\hfill & \ddots \hfill & ⋮\hfill \\ \multicolumn{1}{c}{0_{n+3,1}}& 0_{n+3,2}\hfill & \cdots \hfill & 0_{n+3,m-n}\hfill \end{array}\right]\hfill \end{array}$$

(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 (i− 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 m−n equation sets E _i X _i = 0, as shown in equation (12)

$$\begin{array}{c}\left[\begin{array}{cccc}{e}_{1,i}\hfill & e_{1,i+1}\hfill & \cdots \hfill & e_{1,i+n+3}\hfill \\ \multicolumn{1}{c}{e_{2,i}}& e_{2,i+1}\hfill & \cdots \hfill & e_{2,i+n+3}\hfill \\ \multicolumn{1}{c}{⋮}& ⋮\hfill & \ddots \hfill & ⋮\hfill \\ \multicolumn{1}{c}{e_{n+3,i}}& e_{n+3,i+1}\hfill & \cdots \hfill & e_{n+3,i+n+3}\hfill \end{array}\right]\hfill \\ \left[\begin{array}{c}{q}_{i,i}\hfill \\ \multicolumn{1}{c}{q_{i+1,i}}\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{q_{i+n+3,i}}\end{array}\right]=0\hfill \end{array}$$

(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 ) = Column− 1, the non-zero solution of X _i is unique. If Rank( E _i ) < Column− 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 ) < Column− 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

$${\mathit{E}}_i=\left[\begin{array}{c}0\hfill \\ \multicolumn{1}{c}{\begin{array}{cccc}{e}_{j,i}\hfill & e_{j,i+1}\hfill & \cdots \hfill & e_{j,i+n+3}\hfill \\ \multicolumn{1}{c}{e_{j+1,i}}& e_{j+1,i+1}\hfill & \cdots \hfill & e_{j+1,i+n+3}\hfill \\ \multicolumn{1}{c}{e_{j+2,i}}& e_{j+2,i+1}\hfill & \ddots \hfill & e_{j+2,i+n+3}\hfill \\ \multicolumn{1}{c}{⋮}& ⋮\hfill & ⋮\hfill & ⋮\hfill \end{array}}\\ \multicolumn{1}{c}{0}\end{array}\right]=\left[\begin{array}{c}0\hfill \\ \multicolumn{1}{c}{E_{\mathrm{ij}}}\\ \multicolumn{1}{c}{0}\end{array}\right]$$

(13)

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− 1, the non-zero solution of X _j is unique. If Rank( E _j ) < Column− 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.

Figure 3.

Frame diagram for extracting coefficient matrix.

Figure 4.

Solution flowchart of null space.

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 m−n = 5 equation sets could be formulated and 5 compactly supported wavelets could be constructed.

The specific steps are as follows:

First, 1st–9th, 2nd–10th, 3rd–11th, 4th–12th, and 5th–13th columns of coefficient matrix are taken, respectively, to structure five equation sets.

Second, these five equation sets are, respectively, solved and five untrivial solutions q₁, q₂, q₃, q₄, q₅ are obtained.

Finally, one of the perfect matrix X = [q₁, q₂, q₃, q₄, q₅] has been constructed, as shown in Table 3.

Table 3.

Solution X of null space.

−0.4553583	0	0	0	0
0.6629917	0.14765904	0	0	0
−0.5291775	−0.4040594	0.0101278	0	0
0.25356687	0.53330248	−0.0547781	−2.1616E-5	0
−0.0917187	−0.5760568	0.21925884	0.002677	0
0.01875206	0.40984569	−0.4929074	−0.0362088	−1.1165E-5
−0.0013845	−0.1714297	0.64166143	0.17142972	0.0013845
1.1165E-5	0.03620881	−0.4929074	−0.4098457	−0.0187521
0	−0.002677	0.21925884	0.57605676	0.09171869
0	2.1616E-5	−0.0547781	−0.5333025	−0.2535669
0	0	0.0101278	0.4040594	0.52917749
0	0	0	−0.14765904	−0.6629917
0	0	0	0	0.45535825

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 .

Table 4.

The final reconstruction matrix Q_mn.

−0.68682345	0	0	0	0
1	−0.25632724	0	0	0
−0.79816608	0.701422905	0.015783715	0	0
0.382458587	−0.92578114	−0.08536914	−3.75246E-5	0
−0.13834062	1	0.341704872	0.004647102	0
0.028284005	−0.71146756	−0.76817357	−0.06285633	1.68396E-05
−0.00208825	0.297591712	1	0.297591712	−0.00208825
1.68396E-05	−0.06285633	−0.76817357	−0.71146756	0.028284005
0	0.004647102	0.341704872	1	−0.13834062
0	−3.75246E-5	−0.08536914	−0.92578114	0.382458587
0	0	0.015783715	0.701422905	−0.79816608
0	0	0	−0.25632724	1
0	0	0	0	−0.68682345

According to equation (8), the final reasonable wavelets was constructed, as shown in Figure 5.

Figure 5.

B-splines and constructed wavelets at given scale: (a) B-spline functions at scale factor m = 10, (b) B-spline functions at scale factor n = 5 and (c) wavelet functions at scale factor n = 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 m−n = 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.

Example of multi-resolution fairing for general curve: (a) original curve at m = 117, (b) faired curve at n = 27, (c) detail curve at n = 27, (d) faired curve at n = 17, (e) detail curve at n = 17, (f) faired curve at n = 7, and (g) detail curve at n = 7.

Figure 6(a) is the original curve determined by 120 control vertices, so the related scale factor is m = 120 − 3 = 117. 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 − 3 = 27 and m−n = 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 − 3 = 17 and m−n = 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 −3=7 and m−n=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.

Figure 7.

Analysis of computational efficiency of algorithm.

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.

Figure 8.

(a) Supercharger impeller and (b) a surface of supercharger impeller.

The corresponding original surface about Figure 8(b), as shown in Figure 9(a), has 4000 (50 × 80) control vertices, where there are 50 control vertices in the u direction and 80 control vertices in the v direction, that is, m₁ = 47, m₂ = 77.

Figure 9.

Original surface and faired surface: (a) original surface S_50×80 and (b) faired surface S_18×26.

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.

Figure 10.

Surface fairing framework.

Set target scale n₁ = 15, n₂ = 23, after ARWFA, a faired low-resolution surface determined by 468 (18 × 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 ${\mathit{D}}_{15\times 23}^1$ , mentioned in Figure 10, is one of the detail surface in u direction, as shown in Figure 11(a). ${\mathit{D}}_{15\times 23}^2$ is one of the detail surface in v direction, as shown in Figure 11(b). ${\mathit{D}}_{15\times 23}^3$ is one of the detail surface in uv direction, as shown in Figure 11(c).

Figure 11.

Three high-resolution surfaces: (a) details in u direction ( D ¹); (b) details in v direction ( D ²); and (c) details in uv direction ( D ³).

What needs to be emphasized is that original surface C _m1×m2 could be accurately reconstructed by faired surface C _n1×n2 and three detail surfaces ${\mathit{D}}_{n1\times n2}^1$ , ${\mathit{D}}_{n1\times n2}^2$ and ${\mathit{D}}_{n1\times n2}^3$ .

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 high-quality reconstruction matrix Q _mn could be solved directly without modification of knot vector, different from literature.^18–21 The advantage is that there is no approximate calculation in the process of MRF. Certainly, the related compactly supported wavelets Ψ _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.