#### Abstract

Geometric (or shape) distortion may occur in the data acquisition phase in information systems, and it can be characterized by geometric transformation model. Once the distorted image is approximated by a certain geometric transformation model, we can apply its inverse transformation to remove the distortion for the geometric restoration. Consequently, finding a mathematical form to approximate the distorted image plays a key role in the restoration. A harmonic transformation cannot be described by any fixed functions in mathematics. In fact, it is represented by partial differential equation (PDE) with boundary conditions. Therefore, to develop an efficient method to solve such a PDE is extremely significant in the geometric restoration. In this paper, a novel wavelet-based method is presented, which consists of three phases. In phase 1, the partial differential equation is converted into boundary integral equation and representation by an indirect method. In phase 2, the boundary integral equation and representation are changed to plane integral equation and representation by boundary measure formula. In phase 3, the plane integral equation and representation are then solved by a method we call wavelet collocation. The performance of our method is evaluated by numerical experiments.

#### 1. Introduction

Geometric (or shape) distortion may be produced in the data acquisition phase in pattern recognition, computer vision, and robot vision systems, and it can be characterized by geometric transformation model [1, 2]. An obvious example of such a distortion can be found in a photograph taken by a camera in a computer vision system. In the acquisition phase shown in Figure 1, the trade mark “Coke” is printed on a Coca-Cola bottle, due to the cylindrical shape of the bottle, and the square shape of the trade mark has been changed. This kind of distortion can be characterized by a geometric transformation model, specifically, the biquadratic geometric transformation model in this example [2].

An image is displayed by a set of coordinate points; hence, a geometric transformation can be viewed as the procedure for calculating new coordinate positions of these points under a certain model. The geometric transformation is defined by such that Through this transformation, an image in Cartesian coordinates is transformed into a new image in Cartesian coordinates as shown in Figure 2. The properties of the transformation are determined by functions and ; that is different functions can produce different kinds of geometric transformations.

There are many models of the geometric transformations, which have been widely used in many disciplines [14]. In our previous work, they are categorized into two main classes [1, 2]. If the function is fixed, we call it a fixed transformation model; otherwise we call it an elastic transformation model. The former consists of linear and nonlinear models including bilinear, quadratic, biquadratic, cubic, and bicubic models, while the latter comprises Coons model and harmonic model. These transformation models can be summarized as follows [2]:(A)fixed models:(1)linear models(a)translation,(b)rotation,(c)scaling,(2)nonlinear models:(a)bilinear model,(b)quadratic model,(c)biquadratic model,(d)cubic model,(e)bicubic model,(B)elastic models(1)coons model,(2)harmonic model.

Generally speaking, by the fixed models, a standard image can be transformed into some special shapes. An example of the fixed transformation can be found in Figure 1. More precisely, a fixed transformation has a certain mathematical form to approximate it, which can be described by a fixed class of mathematical functions. In our previous work [1, 2], some significant forms (models) and algorithms have been developed to handle these geometric transformations. For example, the geometric transformation in Figure 1 can be approximated by the biquadratic model, and its mathematical function can be written as where and    denote the new coordinates of the vertices in the quadrangle, which is a distorted image of the trade mark “Coke” in Figure 1.

Once the distorted image is approximated by a certain geometric transformation model, its inverse transformation can be applied to remove the distortion for the geometric restoration. Look at Figure 3, the original trade mark “Coke” is shown in Figure 3(a), and its distorted images are displayed in Figure 3(b). When the inverse biquadratic transformation is applied to these images, the normalized images can be produced and illustrated in Figure 3(c), which are much more easy to be recognized by a recognition system.

Consequently, finding a mathematical form to approximate the distortion of a distorted image plays a key role for the restoration.

Unfortunately, the harmonic geometric transformation model, which converts an image into arbitrary shape, does not have any fixed mathematical forms. An example can be found in Figure 4, where the image of Canadian flag is distorted as shown in Figure 4(b). Note that the shape of the flag is changed, which is so complex that it cannot be described by any fixed model; that is, it cannot be represented by any fixed functions in mathematics. In fact, this model is characterized by other kinds of mathematical formula, that is, partial differential equation. Unlike solving a fixed mathematical formula, solving a partial differential equation is difficult.

The harmonic transformation model is the most important and most complicated one in the geometric transformation models. Actually, all of the other models can be in it. The harmonic model is represented by the partial differential equation (PDE) with boundary conditions.

Let be the region of the elastic plane, where the image is located, and let be its boundary. Suppose that the functions of the transformation of boundary are and . The harmonic transformation satisfies the partial differential equation: where is Laplace's operator; thus the above partial differential equation is called Laplace's equation or harmonic equation.

Accordingly, the task of the restoration is solving the above harmonic equation. We return to the previous example as shown in Figure 4, and the distorted image in Figure 4(b) can be approximated by harmonic transformation. As the corresponding harmonic equation is solved and its inverse transformation is utilized, the restored image can be obtained in Figure 4(c). Therefore, solving the harmonic equation (5) plays a key role in the geometric restoration.

This paper proposes a novel approach based on wavelet analysis to handle the harmonic transformation.

The paper is organized as follows. The existing methods are reviewed in Section 2. The core of the proposed wavelet-based approach is presented in Section 3. As the first step of this method, in Section 3.1, the conversion of the partial differential equation into boundary integral equation and representation is discussed. Two integral methods can be applied; here the indirect method is chosen and the boundary integral equation of the first kind is produced. In Section 3.2, the boundary measure formula, which can change the forms of integral equation and integral representation from boundary to the plane, is presented. Based on the new forms, the wavelet collocation method is used to solve the equation, which is the main task in Section 3.3. A couple of algorithms of IEWC are provided in Section 4. The experiments are illustrated in Section 5. Finally, the conclusions are given in Section 6.

#### 2. Review of the Existing Methods

Two successful approaches have been used to solve this kind of partial differential equation, namely, finite element method [5] and finite difference method [6]. In our previous work [2], the finite element method has been employed solving the equation to handle the harmonic transformation, which gives the following algorithm.

Algorithm 1 (finite element method). One has the following.

Step 1 (discrete region ). Divide region by many small triangular elements such that . For example, region in Figure 5 is divided into twenty two triangular elements based on twelve dots, which produce a lattice .

Step 2 (discrete solution ). Use piecewise linear function to approximate the original solution . Because is linear in each triangular element, therefore is fully determined by the values at lattice . To determine these values, replacing with in the variational form of the equation produces a set of linear equations where the elements in vector are the values of at lattice , is a known matrix, and is a right-hand side vector.

Step 3 (consider the boundary conditions). Find the boundary dots on lattice , and thereafter, the values of on these dots are assigned to satisfy boundary condition , which can be done by modifying the matrix and the right-hand side in the equations .

Step 4 (solve equations). Solve the linear algebraic equations and finally obtain , where indicates the value of at lattice .

In our example, shown in Figure 5, twelve values at lattice are to be determined, seven dots are on the boundary, and the correct values of at these seven boundary dots (dots 1 to 7) are given due to the boundary condition . The values at remainder five inner dots (dots 8 to 12) will be obtained by solving the linear equations.

The finite difference method is another common way to solve partial differential equation numerically. It can change the partial differential equation into a set of corresponding algebraic linear equations. The details can be found in [6].

Both the finite element method and finite difference method have some defects when they will be used in our cases. For example, two issues of the weakness will occur when the finite element method is applied to the harmonic transformation.(i)It depends on the lattice . In finite element method, the values of all points at lattice are solved, even though some points do not lie on the image. For example, in Figure 5, the values of all points at lattice, say points 8–12, are solved. However, only three points (8–10) are the pixels of the pattern, English letter “A.” The values of these points are required to be transformed, while the values of the remainder two points (11 and 12) are not required to be transformed.On the other hand, the values of some pixels on the pattern but not at the lattice are still unknown after solving . For example, in Figure 5, most of the pixels of letter “A,” which are required to be transformed, are not at the lattice. Thus, the values of these pixels are unknown in solution of . Consequently, the interpolation will be used to approximate these values. In this way, the extra cost and error will be brought in.(ii)More attention must be paid to deal with the boundary conditions. Specifically, for a given lattice, we should know which dots are on the boundary and assign them the values satisfying the boundary conditions, as mentioned in Step  3 of Algorithm 1. Therefore, the program code of the lattice and the linear equations is required to be modified, when the domain or lattice is changed.

Thus, more efficient method must be developed to handle the geometric transformation. In this paper, a novel approach called Integral Equation-Wavelet Collocation (IEWC) is presented.

#### 3. Integral Equation-Wavelet Collocation (IEWC) Approach

This is the core section of the paper; a novel approach based on the integral equation and wavelets, called Integral Equation-Wavelet Collocation (IEWC), is presented in this section.

The diagram of the new method is shown in Figure 6. The basic idea of the method is briefly described in Figure 6.

Phase 1. First, the partial differential equation (PDE) (Laplace's equation) is changed into the form of integral equation and integral representation on boundary , which are called boundary integral equation (BIE) and boundary integral representation (BIR), respectively. There are two ways to do so, namely, direct method and indirect method [7, 8]. In this paper, the indirect method is utilized. Mathematically, the process for solving can be written as follows: where is a unknown function, which can be solved in BIE and will be used in BIE. indicates the curvilinear integrate with respect to variable . The first shortcoming, which arises in the finite element method, can be overcome by this way.

Similarly, can also be obtained in the same way, which will be omitted to save the space. In the remainder of this paper, we only discuss .

Phase 2. In order to solve the boundary integral equation (BIE) efficiently and to get rid of the second defect in the finite element method, the boundary measure formula (BMF) is used. It changes the boundary integral equation and boundary integral representation into an integral equation and an integral representation on the whole plane rather than the special boundary . They are called plane integral equation (PIE) and plane integral representation (PIR), respectively. In this way, when is changed, the program code will not to be modified at all. In mathematics, this process can be presented in the following: where stands for the boundary measure, which will be discussed in Section 3.2.

Phase 3. Then, wavelet collocation technique is used to solve the plane integral equation. In the integral equation, the integrand has a discontinuity across boundary . Hence, there is a kind of singularity in it. Fortunately, wavelets have a good property to approximate this kind of singularity. More specifically, suppose that we have a function , which has a discontinuity across curve ; otherwise it is smooth. When a standard Fourier representation is applied to approximate with the best nonzero Fourier terms, , we have When a wavelet representation is used to approximate with the best nonzero wavelet terms, , we can obtain Note that if we compare the right-hand sides of the above two equations, that is, and , it is clear that the rate of the approximation is very slow when Fourier approach is utilized, and the rate of wavelet approximation is better than that of Fourier approximation. Moreover, until now, the wavelet approach is the best result for a fixed nonadaptive method [9].

The last step is the choice of the points to be transformed in the domain and calculate the new coordinates of them by integral representation. The mathematical representation of these operations can be illustrated as follows:

The principal advantages of our method are as follows.(i)The algorithm is divided into two parts, integral equation and integral representation. After solving the plane integral equation, we can choose the pixels to be transformed in the domain arbitrarily and use plane integral representation to evaluate their new coordinates. Therefore, only the pixels on the pattern are transformed to the new coordinate space. In Figure 2, for example, we will choose all pixels on characters “Coke,” and find the new coordinates of them by the integral representation.(ii)The program code is independent of the domain considered; that is, the program code will require no change for the different kind of boundaries. It benefits from the boundary measure formula. In fact, we do not need the function of boundary at all. What we really need are the original coordinates of the pixels at the boundary and the coordinates of the new ones.Finally, the summary of this novel approach can be illustrated in the following mathematical expression:

##### 3.1. Boundary Integral Method

According to the boundary integral method, the first phase of the IEWC converts the partial differential equation (PDE) (5) which is recalled below, into a boundary integral equation (BIE) and a boundary integral representation (BIR).

There are two kinds of boundary integral equations, namely, (1) integral equation of the first kind and (2) integral equation of the second kind.

If the equation takes the form of it is called integral equation of the first kind, where and are known functions and indicates the integrate with variable . Otherwise, the equation with the form is called an integral equation of the second kind, where is a known constant.

Most of the researchers work on the integral equations of the second kind in both of the theories and applications, because some integral equations of the first kind are quite ill-conditioned; that is, the speed of the convergent will be slow if an iterative algorithm is used. However, the integral equations of the first kind have been an increasingly popular approach to solve various boundary value problems [7]. In this paper, the boundary integral equations of the first kind are applied to facilitate the use of the boundary measure formula, as well as the further application of the wavelet collocation method.

Furthermore, there are two ways to convert the partial differential equation (PDE), into the boundary integral equation of the first kind, namely, (1) direct method and (2) indirect method [7, 8]. They will briefly be presented in the following.

(1) Direct Method. In this method, based on Green's formula, where , , , and is the unit outer normal vector; the partial differential equation (5) can be changed to a system with a boundary integral equation and a boundary integral representation. Thereafter, the following main steps are done. Step 1. Solve the integral equation to find on such that Step 2. , calculate by formula

(2) Indirect Method. This method is based on the potential theory, and the procedure of solving (5) is presented in the following.Step 1. Solve the integral equation to find on such that Step 2. Calculate by formula

It is obvious that either the pair of (18), (19) or (20), (21) is equivalent to (5). The latter is utilized in this paper, because the right-hand sides of these equations are very simple and can be solved easily. Furthermore, to ensure the uniqueness of the solution, we tacitly assume that the interior of domain has the property of As to the existence of the solution, [10] has proved that (20) is equivalent to and for arbitrary function and constant , (23) has unique solutions and , which ensures that (20) is of viability (i.e., unique solvable).

The above discussion gives us a hint that we can first obtain the data from (20) and then use (21) to calculate at any pixel in the domain . Note that the pixel is chosen arbitrarily in the integral representation, which makes us free from the lattice built in either the finite element method or finite difference one.

If the function of the boundary can be determined, that is, is parameterized, (20) can be changed into one-dimensional integral equation on a closed interval. Thereafter, it can be solved by the classical methods or periodic wavelet method [1116]. In this way, the program code is dependent on the representation of curve . Unfortunately, in most of the cases, the function of the boundary is unknown as we mentioned in Section 1. In order to develop a new method, which can avoid knowing the function of boundary , the boundary measure formula is utilized in our study to change the forms of (20) and (21) into other suitable forms.

##### 3.2. Boundary Measure Formula

Based on the boundary measure formula, the second phase of the IEWC converts the boundary integral equation (BIE) and boundary integral representation (BIR) into the plane integral equation (PIE) and plane integral representation (PIR).

Assume that is a bounded domain in , whose boundary can be presented by Lipschitz function , and denotes its characteristic function, which has value 1 if the point is in domain ; otherwise is 0. The unit normal vector along can be written as where is the gradient of and stands for its norm. They can generalize the vector-valued measure and the Radon measure, respectively, if is only of bounded variation over domain [17]. Hence, the boundary measure formula can be described below.

Theorem 2. For any integrable function defined on , after extending from to , we have where is called boundary measure.

Proof. In fact, we can derive This establishes (25).

Reference [18] has proved that the gradient of the characteristic functions and can be approximated by Daubechies scale or wavelet function in the Sobolev space or space , respectively. That ensures that can be approximated by Daubechies scale or wavelet function if is integrable. We will use this formula in our new approach. It should be noted that has singularities along boundary , therefore the wavelet representation is more effective to handle this problem than the Fourier representation due to its sparse representation of singularities. For example, to represent an edge, a type of singularity, with error , roughly speaking, requires sinusoids in Fourier form but needs only about wavelet items in wavelet representation. As we discussed in Section 1, the wavelets outperform the classical method. That is the main reason we use wavelet here.

Using the boundary measure formula, the boundary integral equation (20) and the boundary integral representation (21) become the following plane integral equation and plane integral representation:

Theorem 3. Plane integral representation (29) is the solution of partial differential equation (5), where is the solution of plane integral equation (28).

Proof. The proof of this theorem is omitted in this paper to avoid the complicated mathematics and to save space.

##### 3.3. Wavelet Collocation Technique

Wavelet analysis has been widely applied in image processing [1922]. In this paper, the wavelet theory is used in IEWC method, in which, after phase 2, the plane integral equation (PIE) and plane integral representation (PIR) have been obtained; thereafter, in phase 3, the wavelet collocation technique is employed to arrive the solution.

Let be the Daubechies scale function and Based on the boundary measure formula (25), the task now is to solve plane integral equation: Recall that the support of the gradient of the characteristic function is a compact domain in ; in fact, it is a tubular neighborhood of . Therefore, we need only finite number of scale functions to represent . Let be an index set, and let be the cardinal of . The key point here is solving the product instead of solving the unknown function . In fact, when we know , then from integral representation (29), we can obtain for any . That is why we use integral equation of the first kind in this paper. As can be approximated by substituting in (28) with produces the error representation : Choose collocation points , and let , and then we have That is, That is, Equation (36) is a set of linear algebraic equations, which can be rewritten in form of from which, is obtained, where matrix with entries of is and the right hand side vector is

Remark 4. Note that matrix does not depend on , except points that should be chosen on . Meantime, the right-hand side in (37) does not need the representation of function , except values of the discrete pixels on the boundary . These points are called tiepoints. It indicates that, in the implementation, the program code is independent of the boundary.

Remark 5. The condition number of matrix may be large, because it is from the integral equation of the first kind; thus, usually the Tikhonov's regularize method is used to solve it [23]. In our new approach, linear algebraic equations (37) are generated by the wavelet collocation technique. Therefore, the diagonal preconditioning method [14] can be employed to reduce the condition number of matrix . It leads us to use a very easy way to solve (37) instead of using Tikhonov's regularize method. This approach was first used in the Galerkin method, and thereafter, in 1995, Schneider proved that it also can be applied to the collocation method [24].

Now we can use integral representation (29) to evaluate the approximate value of at any point in ; that is, We use notation to represent the integral expression in (40); that is and at last, we have

#### 4. Algorithms of the IEWC Approach

In this section, the wavelet-based algorithms of the IEWC approach are presented in both of the general version and detailed version. To facilitate the implementation of the algorithm, the computation of matrix, which will be useful in performance of the algorithm, is also discussed in this section.

Algorithm 6 (boundary measure and wavelet (general)). One has the following.

Step 1. Choose the collocation points on , and evaluate matrixes and in (37) with the formula (38).

Step 2. Solve (37) with least square method to obtain coefficients .

Step 3. Choose a point in the domain needed to be transformed to new coordinate, and calculate the coefficients with the formula (41).

Step 4. Obtain the approximation using (41).

##### 4.1. Computation of the Matrix

In the Step  1 and Step  3 above, the main cost is the computation of In this subsection we will give a simple method to calculate it approximately. Let us first introduce some notations.

Definition 7. If a quadrature formula holds for any polynomial of degree less than or equal to , then it is called Gauss-type quadrature with scale function as its weight function, the dots , and the coefficients are called generalized Gauss-type quadrature dots and generalized Gauss-type quadrature weights.

Similar to the classical Gauss-type quadrature formula, we recall that is the support of the scale function ; therefore, we have the following.

Proposition 8. In formula (44), the necessary and sufficient condition for being generalized Gauss-type quadrature dots is being orthogonal to any polynomial with as the weight function; that is,

Proposition 9. In formula (44), if , the error of (44) is where .

A significant task is to determine and in formula (44). Let , , we obtain a nonlinear system [25]: First, we compute the right-hand side of (47) recursively. Suppose that , from the theory of wavelet, we can write from which we know and it can be computed recursively using [26] with , where reals are coefficients in . So far, the right-hand side in (47) is computed. Then, we solve the nonlinear system (47), so that coefficients and are obtained. For example, when , , we have and obtain Table 1.

Therefore, we have a very simple formula For , by use of spline function, we can prove that its error, , is Based on these results, we know that can be used to approximate and we will use it in our experiments.

Algorithm 10 (boundary measure and wavelet (detail)). Input and .Step 1Choose a resolution level and fix a rectangular domain covering .Calculate the index set by .Choose any collocation points .Step 2For to DoCompute matrix elements Obtain right-hand side from .End.Step 3Solve the linear systems and obtain .Step 4Choose the points in the domain to be transformed to new coordinate space.Step 5For to Do:Evaluate Calculate the new coordinate End.Output .

#### 5. Experiments

Experiments have been conducted to evaluate the performances of the new approach.

A couple of numerical experiments using the wavelet-based IEWC approach is presented in this section.

Experiment 1 (nonlinear-harmonic and its inverse). The harmonic transformation satisfies (5), which is recalled as follows: In this experiment, the equations of the boundary conditions , and the region are specified by the following forms: Geometric transform is viewed to map the location of input image to a location in the output image; it defines how the pixel values in the output image are to be arranged. The distorted image coordinates can be defined by the equations The primary idea is to find a mathematical model for the geometric distortion process, specifically, the equations and then apply the inverse process to find the restored image. To determine the necessary equations, we need to identify a set of points in the original image that matches another set of points in the distorted image. These sets of points are called tiepoints. In this experiment, the relationship between these tiepoints in two images is determined by the boundary condition, which is described in (60), and the proposed IEWC method can be used.
We use Figure 7 as an example, where Figure 7(a) is the original image. As the IEWC approach is applied to solve the harmonic transformation with the boundary condition (60), the original image is distorted and displayed in Figure 7(b). To achieve the restoration, then the CSIM method is used to perform the inverse transformation, and the restored image is illustrated in Figure 7(c). The details of the CSIM method can be found in our previous work [1].
To quantify and clarify the advantages of IEWC over the traditional approach, a comparison is given in Figure 8.
In the traditional approach, a harmonic distorted image is approximated by piecewise bilinear model [27], as shown in Figure 8. In the bilinear model, four points of each subregion are used to generate the equation: where , , are constants to be determined by solving the eight simultaneous equations. Because we have defined four tiepoints, thus we have eight equations, where , , , and are known. Now we can solve the eight equations for the eight unknowns , so that the necessary equations for the coordinate mapping can be obtained. The mapping equations , then are applied to all pairs needed. In practice, the values of , , , and are not likely to be integers. The simplest solution is the nearest neighbor method, where the pixel is assigned to the value of the closest pixel in the image. The restored image which is constructed by the classical piecewise bilinear method is presented in Figure 7(d), in which, 128 tiepoints are used on the boundary. This method does not necessarily provide optimal results.
To compare the above two restored images, we project both of them, as well as the original image and the distorted one, to the -axis and -axis. The results are shown in Figures 9 and 10, respectively. From these projections, it is clear that the results of IEWC approach are better than that of the traditional one.

Experiment 2 (nonlinear-harmonic and its inverse for a circle area). Consider In this experiment, the domain of the partial differential equation (PDE) is a circle area. Because the area is not a quadrilateral (four-sided polygon), it is difficult to use the piecewise bilinear method. It is shown that the program code of the IEWC approach is independent of the boundary form. The results are shown in Figure 11. To compare the results, the projections of them to the axis are shown in Figures 12 and 13.

#### 6. Conclusions

Usually the piecewise bilinear model can be used to all geometric degradation image regardless the different characteristics of images. We have found a new model (partial differential equation PDE) for the transformation in our previous work [1, 2] and aim to construct an efficient algorithm (IEWC) to solve the PDE in this paper.

In this paper, we have presented a wavelet-based approach for the harmonic transformation. Unlike the traditional methods, in the IEWC approach, the pixels needed to be transformed to new coordinates can be chosen arbitrarily. Meanwhile, the program code of our method is independent of the boundary, and we need only a set of the original coordinates of the pixels on the boundary of the image as well as their new coordinates in the transformed image. To make the algorithm more efficiency, Daubechies wavelet (scale) functions and a Gauss-type quadrature formula have been used. Different examples have been tested with the anticipated results.

Some further work is still under study, which is presented below. The tiepoints are unknown and should be guessed or calculated by some other methods. In this paper, a GA approach has been applied to extract the outer contours and find the tiepoints. In our further work, other GA, wavelet-based method, local search, immune approach, and so forth will be used to find their usage in this direction to construct a good interpolation method.

In addition, as discussed above, there are singularities along the boundary, which can be treated efficiently by wavelet. More recently, Donoho [9] has constructed a tool called curvelet to handle this kind of singularity, which is built from Meyer wavelet basis. In our further work, the curvelet will be utilized. In this way, the boundary measure will be approximated by the curvelet with the same accuracy as we use wavelet. Meantime, compared with wavelet or sinusoid basis, fewer terms will be computed when the curvelet will be used. It will make our algorithm more efficient.

Besides the further improvement of the algorithm, the proposed method will be combined with other techniques to solve more sophisticated problems. For example, when we restore the distorted image, some blurred pictures may occur. To obtain a clean image, which will be easier to be recognized by a pattern recognition system, the fusion technique will be applied.

#### Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

#### Acknowledgments

This work was financially supported by the Multi-Year Research of University of Macau under Grants no. MYRG205(Y1-L4)-FST11-TYY and no. MYRG187(Y1-L3)-FST11-TYY and by the National Natural Science Foundation of China under Grant no. 61273244. This research project was also supported by the Science and Technology Development Fund (FDCT) of Macau, under Contract nos. 100-2012-A3, 026-2013-A, and also by Guangxi Science and Technology Fund of China, under Contract no. 201203YB179.