# Integral Equation-Wavelet Collocation Method for Geometric Transformation and Application to Image Processing

Authors: Lina Yang; Yuan Yan Tang; Xiang Chu Feng; Lu Sun
Journal: Abstract and Applied Analysis
 (2014)

---

## 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.


---

## Body

## 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].Figure 1
Example of the biquadratic geometric transformation.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(1)T:(x1,x2)⟶(u,v),
such that
(2)u=u(x1,x2),v=v(x1,x2).
Through this transformation, an image in Cartesian coordinates x1Ox2 is transformed into a new image in Cartesian coordinates uOv as shown in Figure 2. The properties of the transformation T are determined by functions u=u(x1,x2) and v=v(x1,x2); that is different functions can produce different kinds of geometric transformations.Figure 2
An image in Cartesian coordinatesx1Ox2 is transformed into a new image in Cartesian coordinates uOv by geometric transformation T.There are many models of the geometric transformations, which have been widely used in many disciplines [1–4]. 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 Figure1. 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
(3)[uv]=[[1-x1x1]00[1-x1x1]]×[[XP1XP4XP2XP3][YP1YP4YP2YP3]][1-x2x2],
where XPi and YPi  (i=1,2,3,4) 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 Figure3, 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.Figure 3
Restoration of image from the biquadratic distortion by inverse transformation.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 Figure4, 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.Figure 4
Harmonic distortion and restoration.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 u=f(x1,x2) and v=g(x1,x2). The harmonic transformation
(4)T:(x1,x2)⟶(u,v)
satisfies the partial differential equation:
(5)Δu(x1,x2)=0,(x1,x2)∈Ω,u|Γ=f(x1,x2),(x1,x2)∈Γ,Δv(x1,x2)=0,(x1,x2)∈Ωv|Γ=g(x1,x2),(x1,x2)∈Γ,
where
(6)Δ:∂2∂x12+∂2∂x22
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 Figure4, 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 Section2. 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. Twointegral 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 ei such that Ω≈∪iei. For example, region Ω in Figure 5 is divided into twenty two triangular elements based on twelve dots, which produce a lattice Π.Figure 5
Finite element method.Step
2 (discrete solution
u
). Use piecewise linear function uh to approximate the original solution u. Because uh is linear in each triangular element, therefore uh is fully determined by the values at lattice Π. To determine these values, replacing u with uh in thevariational form of the equation produces a set oflinear equations(7)Ku¯h=f¯h,
where the elements in vector u¯h are the values of uh at lattice Π, K is a known matrix, and f¯h is a right-hand side vector.Step
3 (consider the boundary conditions). Find the boundary dots on lattice Π, and thereafter, the values of uh on these dots are assigned to satisfy boundary condition f(x1,x2), which can be done by modifying the matrix K and the right-hand side f¯h in the equations Ku¯h=f¯h.Step
4 (solve equations). Solve the linear algebraic equations and finally obtain u|Π≈uh|Π=u¯h, where u|Π indicates the value of u at lattice Π.In our example, shown in Figure5, twelve values at lattice Π are to be determined, seven dots are on the boundary, and the correct values of uh at these seven boundary dots (dots 1 to 7) are given due to the boundary condition f(x1,x2). 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 solvingKu¯h=f¯h. 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 u|Π≈uh|Π=u¯h. 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 Algorithm1. 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 calledIntegral 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 Figure6. The basic idea of the method is briefly described in Figure 6.Figure 6
Diagram of the new approach.Phase 1. First, thepartial differential equation (PDE) (Laplace's equation) is changed into the form of integral equation and integral representation on boundary Γ, which are calledboundary integral equation (BIE) andboundary integral representation (BIR), respectively. There are two ways to do so, namely, direct method and indirect method [7, 8]. In this paper, theindirect method is utilized. Mathematically, the process for solving u can be written as follows:
(8)PDE:{Δu(x1,x2)=0,∀(x1,x2)∈Ωu|Γ=f(x1,x2),∀(x1,x2)∈Γ,⇓BoundaryIntegral(Indirect)Method⇓BIE:∫Γ‍ω(y)log⁡|x-y|dsy=f(x),∀(x1,x2)∈ΓBIR:u(x)=∫Γ‍ω(y)⋯dsy,∀(x1,x2)∈ω,
where ω(y) is a unknown function, which can be solved in BIE and will be used in BIE. dsy indicates the curvilinear integrate with respect to variable y=(y1,y2). The first shortcoming, which arises in the finite element method, can be overcome by this way.Similarly,v 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 u.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, theboundary 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 R2 rather than the special boundary Γ. They are calledplane integral equation (PIE) andplane 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:
(9)BIE:∫Γ‍ω(y)log⁡|x-y|dsy=f(x),∀(x1,x2)∈ΓBIR:u(x)=∫Γ‍ω(y)⋯dsy,∀(x1,x2)∈ω⇓BoundaryMeasureFormula(BMF)⇓PIE:∫R2‍ω(y)∥∂Ω∥log⁡|x-y|dsy=f(x),55555555555555555555555∀(x1,x2)∈ΓPIR:u(x)=∫R2‍ω(y)∥∂Ω∥⋯dsy,∀(x1,x2)∈ω,
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 f, which has a discontinuity across curve Γ; otherwise it is smooth. When a standard Fourier representation is applied to approximate f with the best m nonzero Fourier terms, f~mF, we have
(10)∥f-f~mF∥22≍m-1/2,m⟶∞.
When a wavelet representation is used to approximate f with the best m nonzero wavelet terms, f~mW, we can obtain
(11)∥f-f~mW∥22≍m-1,m⟶∞.
Note that if we compare the right-hand sides of the above two equations, that is, m-1/2 and m-1, 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:(12)PIE:∫R2‍ω(y)∥∂Ω∥log⁡|x-y|dsy=f(x),55555555555555555555555∀(x1,x2)∈ΓPIR:u(x)=∫R2‍ω(y)∥∂Ω∥⋯dsy,∀(x1,x2)∈ω⇓WaveletCollocation(WC)Technique⇓∑p,qh(p,q)j∫Ry2‍ϕpj(y1)ϕqj(y2)log⁡|xk-y|dy-f(xk)=0,∀(x1,x2)∈Γuj(x)=∑p,qh(p,q)j∫Ry2‍ϕpj(y1)ϕqj(y2)⋯dy,5555555555555555555555∀(x1,x2)∈ω⇓⇓Solutions.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 Figure2, 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:
(13)PDE:{Δu(x1,x2)=0,∀(x1,x2)∈Ωu|Γ=f(x1,x2),∀(x1,x2)∈Γ,⇓⇓BoundaryIntegral(Indirect)Method⇓⇓BIE:∫Γ‍ω(y)log⁡|x-y|dsy=f(x),∀(x1,x2)∈ΓBIR:u(x)=∫Γ‍ω(y)⋯dsy,∀(x1,x2)∈ω⇓⇓BoundaryMeasureFormula(BMF)⇓⇓PIE:∫R2‍ω(y)∥∂Ω∥log⁡|x-y|dsy=f(x),555555555555555555555555∀(x1,x2)∈ΓPIR:u(x)=∫R2‍ω(y)∥∂Ω∥⋯dsy,∀(x1,x2)∈ω⇓⇓WaveletCollocation(WC)Method⇓⇓∑p,qh(p,q)j∫Ry2‍ϕpj(y1)ϕqj(y2)log⁡|xk-y|dy-f(xk)=0,55555555555555555555555555555555∀(x1,x2)∈Γuj(x)=∑p,qh(p,q)j∫Ry2‍ϕpj(y1)ϕqj(y2)⋯dy,555555555555555555555∀(x1,x2)∈ω⇓⇓Solutions.


### 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,
(14)Δu=0,∀x=(x1,x2)∈Ω,u|Γ=f(x1,x2),∀x=(x1,x2)∈Γ,
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(15)∫Γ‍K(x,y)u(y)dsy=f(x),∀x=(x1,x2)∈Γ,
it is calledintegral equation of the first kind, where K(x,y) and f(x) are known functions and dsy indicates the integrate with variable y=(y1,y2). Otherwise, the equation with the form
(16)λu(x)-∫Γ‍K(x,y)u(y)dsy=f(x),∀x=(x1,x2)∈Γ
is called anintegral equation of the second kind, where λ is a known constant.Most of the researchers work on the integral equations of thesecond kind in both of the theories and applications, because some integral equations of thefirst 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 thefirst kind have been an increasingly popular approach to solve various boundary value problems [7]. In this paper, the boundary integral equations of thefirst 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 onGreen's formula,(17)∫Ω‍∇·FdΩ=∫Γ‍F·ndS,
where ∇=∂/∂x1+∂/∂x2, F=(f1(x),f2(x)), x=(x1,x2), and n 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. Step1. Solve the integral equation to find ω on Γ such that
(18)∫Γω(y)log⁡|x-y|dsy=-πf(x)+∫Γf(y)∂∂nlog⁡|x-y|dsy,5555555555555∀x=(x1,x2)∈Γ.Step2. ∀x∈Ω, calculate u(x) by formula
(19)u(x)=12π∫Γ(f(y)∂∂nlog⁡|x-y|-ω(y)log⁡|x-y|)dsy.(
2) Indirect Method. This method is based on thepotential theory, and the procedure of solving (5) is presented in the following.Step1. Solve the integral equation to find ω on Γ such that
(20)∫Γ‍ω(y)log⁡|x-y|dsy=f(x),∀x=(x1,x2)∈Γ.Step2. Calculate u(x) by formula
(21)u(x)=∫Γ‍ω(y)log⁡(|x-y|)dsy,∀x=(x1,x2)∈Ω.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
(22)diameter(Ω)<1.
As to the existence of the solution, [10] has proved that (20) is equivalent to
(23)∫Γ‍ω(y)log⁡|x-y|dsy=f(x)+c,∀x∈Γ,∫Γ‍ω(y)dsy=A,
and for arbitrary function f and constant A, (23) has unique solutions ω and c, 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 u(x) at any pixel x=(x1,x2) in the domain Ω. Note that the pixel x∈Ω 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 [11–16]. 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 Γ, theboundary 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 R2, whose boundary Γ can be presented by Lipschitz function F(x1,x2)=c, and χΩ denotes its characteristic function, which has value 1 if the point (x1,x2) is in domain Ω; otherwise is 0. The unit normal vector along Γ can be written as
(24)n=∇F|∇F|,
where ∇F is the gradient of F and |∇F| stands for its norm. They can generalize the vector-valued measure and the Radon measure, respectively, if F is only of bounded variation over domain Ω [17]. Hence, the boundary measure formula can be described below.Theorem 2.
For any integrable functionf defined on Γ, after extending f from Γ to R2, we have
(25)∫Γ‍fds=∫R2‍f∥∂Ω∥dx,
where ∥∂Ω∥=-∇χΩ·n is called boundary measure.Proof.
In fact, we can derive(26)∫Γ‍fds=∫Γ‍fn·nds=∫Ω‍∇f·ndx+∫Ω‍fdiv⁡ndx=∫R2‍∇f·χΩndx+∫R2‍fχΩdiv⁡ndx=∫R2‍∇f·χΩndx-∫R2‍∇(fχΩ)·ndx=-∫R2‍f∇χΩ·ndx=∫R2‍f∥∂Ω∥dx.
This establishes (25).Reference [18] has proved that the gradient of the characteristic functions ∇χΩ and n can be approximated by Daubechies scale or wavelet function in the Sobolev space H-1(Ω) or space H1(Ω), respectively. That ensures that
(27)f∥∂Ω∥=-f∇χΩ·n
can be approximated by Daubechies scale or wavelet function if f is integrable. We will use this formula in our new approach. It should be noted that ∥∂Ω∥=-∇χΩ·n 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 1/N, roughly speaking, requires N2 sinusoids in Fourier form but needs only about N 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:
(28)∫Ry2‍ω∥∂Ω∥log⁡|x-y|dy=f(x),∀x∈Γ,(29)u(x)=∫Ry2‍ω∥∂Ω∥log⁡|x-y|dy,∀x∈Ω.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 [19–22]. 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
(30)ϕpj(t)=2j/2ϕ(2jt-p).
Based on the boundary measure formula (25), the task now is to solve plane integral equation:
(31)∫Ry2‍ω∥∂Ω∥log⁡|x-y|dy=f(x),∀x∈Γ.
Recall that the support of the gradient of the characteristic function ∇χΩ is a compact domain in R2; in fact, it is a tubular neighborhood of Γ. Therefore, we need only finite number of scale functions to represent ω∥∂Ω∥=-ω∇χΩ·n. 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 u(x) for any x∈Ω. That is why we use integral equation of the first kind in this paper. As ω∥∂Ω∥ can be approximated by
(32)(ω∥∂Ω∥)j=∑(p,q)∈Λh(p,q)jϕpj(y1)ϕqj(y2),
substituting ω∥∂Ω∥ in (28) with (ω∥∂Ω∥)j produces the error representation e(x):
(33)e(x)=∫Ry2‍(ω∥∂Ω∥)jlog⁡|x-y|dy-f(x),∀x∈Γ.
Choose |Λ| collocation points {xk}⊂Γ, and let e(xk)=0,1≤k≤|Λ|, and then we have
(34)e(xk)=∫Ry2‍(ω∥∂Ω∥)jlog⁡|xk-y|dy-f(xk)=0,555555555555555555555555555551≤k≤|Λ|.
That is,
(35)∫Ry2‍(ω∥∂Ω∥)jlog⁡|xk-y|dy=f(xk),1≤k≤|Λ|.
That is,
(36)∑(p,q)∈Λh(p,q)j∫Ry2‍ϕpj(y1)ϕqj(y2)log⁡|xk-y|dy=f(xk).
Equation (36) is a set of linear algebraic equations, which can be rewritten in form of
(37)Kh¯=f¯,
from which, {h(p,q)j} is obtained, where matrix K with entries of |Λ|×|Λ| is
(38)Kk,(p,q)=∫Ry2ϕpj(y1)ϕqj(y2)log⁡|xk-y|dy,5555555551≤(p,q),k≤|Λ|,xk∈Γ,
and the right hand side vector f¯ is
(39)fk=f(xk),1≤k≤|Λ|.Remark 4.
Note that matrixK does not depend on Γ, except points {xk}k=1|Λ| that should be chosen on Γ. Meantime, the right-hand side in (37) does not need the representation of function f(x), except values {f(xk)}k=1|Λ| of the discrete pixels {xk}k=1|Λ| 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 matrixK 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, thediagonal preconditioning method [14] can be employed to reduce the condition number of matrix K. 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 u(x) at any point x in Ω; that is,
(40)u(x)=∫Γωlog⁡|x-y|dsy=∫Ry2ω∥∂ω∥log⁡|x-y|dy≈∑(p,q)∈Λ‍h(p,q)j∫Ry2ϕpj(y1)ϕqj(y2)log⁡|x-y|dy=uj(x),∀x∈Ω.
We use notation Q(p,q)(x) to represent the integral expression in (40); that is
(41)Q(p,q)(x)=∫Ry2‍ϕpj(y1)ϕqj(y2)log⁡|x-y|dy,∀x∈Ω,
and at last, we have
(42)uj(x)=∑(p,q)∈Λh(p,q)jQ(p,q),∀x∈Ω.


## 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,
(14)Δu=0,∀x=(x1,x2)∈Ω,u|Γ=f(x1,x2),∀x=(x1,x2)∈Γ,
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(15)∫Γ‍K(x,y)u(y)dsy=f(x),∀x=(x1,x2)∈Γ,
it is calledintegral equation of the first kind, where K(x,y) and f(x) are known functions and dsy indicates the integrate with variable y=(y1,y2). Otherwise, the equation with the form
(16)λu(x)-∫Γ‍K(x,y)u(y)dsy=f(x),∀x=(x1,x2)∈Γ
is called anintegral equation of the second kind, where λ is a known constant.Most of the researchers work on the integral equations of thesecond kind in both of the theories and applications, because some integral equations of thefirst 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 thefirst kind have been an increasingly popular approach to solve various boundary value problems [7]. In this paper, the boundary integral equations of thefirst 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 onGreen's formula,(17)∫Ω‍∇·FdΩ=∫Γ‍F·ndS,
where ∇=∂/∂x1+∂/∂x2, F=(f1(x),f2(x)), x=(x1,x2), and n 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. Step1. Solve the integral equation to find ω on Γ such that
(18)∫Γω(y)log⁡|x-y|dsy=-πf(x)+∫Γf(y)∂∂nlog⁡|x-y|dsy,5555555555555∀x=(x1,x2)∈Γ.Step2. ∀x∈Ω, calculate u(x) by formula
(19)u(x)=12π∫Γ(f(y)∂∂nlog⁡|x-y|-ω(y)log⁡|x-y|)dsy.(
2) Indirect Method. This method is based on thepotential theory, and the procedure of solving (5) is presented in the following.Step1. Solve the integral equation to find ω on Γ such that
(20)∫Γ‍ω(y)log⁡|x-y|dsy=f(x),∀x=(x1,x2)∈Γ.Step2. Calculate u(x) by formula
(21)u(x)=∫Γ‍ω(y)log⁡(|x-y|)dsy,∀x=(x1,x2)∈Ω.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
(22)diameter(Ω)<1.
As to the existence of the solution, [10] has proved that (20) is equivalent to
(23)∫Γ‍ω(y)log⁡|x-y|dsy=f(x)+c,∀x∈Γ,∫Γ‍ω(y)dsy=A,
and for arbitrary function f and constant A, (23) has unique solutions ω and c, 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 u(x) at any pixel x=(x1,x2) in the domain Ω. Note that the pixel x∈Ω 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 [11–16]. 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 Γ, theboundary 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 R2, whose boundary Γ can be presented by Lipschitz function F(x1,x2)=c, and χΩ denotes its characteristic function, which has value 1 if the point (x1,x2) is in domain Ω; otherwise is 0. The unit normal vector along Γ can be written as
(24)n=∇F|∇F|,
where ∇F is the gradient of F and |∇F| stands for its norm. They can generalize the vector-valued measure and the Radon measure, respectively, if F is only of bounded variation over domain Ω [17]. Hence, the boundary measure formula can be described below.Theorem 2.
For any integrable functionf defined on Γ, after extending f from Γ to R2, we have
(25)∫Γ‍fds=∫R2‍f∥∂Ω∥dx,
where ∥∂Ω∥=-∇χΩ·n is called boundary measure.Proof.
In fact, we can derive(26)∫Γ‍fds=∫Γ‍fn·nds=∫Ω‍∇f·ndx+∫Ω‍fdiv⁡ndx=∫R2‍∇f·χΩndx+∫R2‍fχΩdiv⁡ndx=∫R2‍∇f·χΩndx-∫R2‍∇(fχΩ)·ndx=-∫R2‍f∇χΩ·ndx=∫R2‍f∥∂Ω∥dx.
This establishes (25).Reference [18] has proved that the gradient of the characteristic functions ∇χΩ and n can be approximated by Daubechies scale or wavelet function in the Sobolev space H-1(Ω) or space H1(Ω), respectively. That ensures that
(27)f∥∂Ω∥=-f∇χΩ·n
can be approximated by Daubechies scale or wavelet function if f is integrable. We will use this formula in our new approach. It should be noted that ∥∂Ω∥=-∇χΩ·n 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 1/N, roughly speaking, requires N2 sinusoids in Fourier form but needs only about N 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:
(28)∫Ry2‍ω∥∂Ω∥log⁡|x-y|dy=f(x),∀x∈Γ,(29)u(x)=∫Ry2‍ω∥∂Ω∥log⁡|x-y|dy,∀x∈Ω.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 [19–22]. 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
(30)ϕpj(t)=2j/2ϕ(2jt-p).
Based on the boundary measure formula (25), the task now is to solve plane integral equation:
(31)∫Ry2‍ω∥∂Ω∥log⁡|x-y|dy=f(x),∀x∈Γ.
Recall that the support of the gradient of the characteristic function ∇χΩ is a compact domain in R2; in fact, it is a tubular neighborhood of Γ. Therefore, we need only finite number of scale functions to represent ω∥∂Ω∥=-ω∇χΩ·n. 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 u(x) for any x∈Ω. That is why we use integral equation of the first kind in this paper. As ω∥∂Ω∥ can be approximated by
(32)(ω∥∂Ω∥)j=∑(p,q)∈Λh(p,q)jϕpj(y1)ϕqj(y2),
substituting ω∥∂Ω∥ in (28) with (ω∥∂Ω∥)j produces the error representation e(x):
(33)e(x)=∫Ry2‍(ω∥∂Ω∥)jlog⁡|x-y|dy-f(x),∀x∈Γ.
Choose |Λ| collocation points {xk}⊂Γ, and let e(xk)=0,1≤k≤|Λ|, and then we have
(34)e(xk)=∫Ry2‍(ω∥∂Ω∥)jlog⁡|xk-y|dy-f(xk)=0,555555555555555555555555555551≤k≤|Λ|.
That is,
(35)∫Ry2‍(ω∥∂Ω∥)jlog⁡|xk-y|dy=f(xk),1≤k≤|Λ|.
That is,
(36)∑(p,q)∈Λh(p,q)j∫Ry2‍ϕpj(y1)ϕqj(y2)log⁡|xk-y|dy=f(xk).
Equation (36) is a set of linear algebraic equations, which can be rewritten in form of
(37)Kh¯=f¯,
from which, {h(p,q)j} is obtained, where matrix K with entries of |Λ|×|Λ| is
(38)Kk,(p,q)=∫Ry2ϕpj(y1)ϕqj(y2)log⁡|xk-y|dy,5555555551≤(p,q),k≤|Λ|,xk∈Γ,
and the right hand side vector f¯ is
(39)fk=f(xk),1≤k≤|Λ|.Remark 4.
Note that matrixK does not depend on Γ, except points {xk}k=1|Λ| that should be chosen on Γ. Meantime, the right-hand side in (37) does not need the representation of function f(x), except values {f(xk)}k=1|Λ| of the discrete pixels {xk}k=1|Λ| 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 matrixK 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, thediagonal preconditioning method [14] can be employed to reduce the condition number of matrix K. 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 u(x) at any point x in Ω; that is,
(40)u(x)=∫Γωlog⁡|x-y|dsy=∫Ry2ω∥∂ω∥log⁡|x-y|dy≈∑(p,q)∈Λ‍h(p,q)j∫Ry2ϕpj(y1)ϕqj(y2)log⁡|x-y|dy=uj(x),∀x∈Ω.
We use notation Q(p,q)(x) to represent the integral expression in (40); that is
(41)Q(p,q)(x)=∫Ry2‍ϕpj(y1)ϕqj(y2)log⁡|x-y|dy,∀x∈Ω,
and at last, we have
(42)uj(x)=∑(p,q)∈Λh(p,q)jQ(p,q),∀x∈Ω.


## 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 {xk}k=1|Λ| on Γ, and evaluate matrixes K and f¯ in (37) with the formula (38).Step
2. Solve (37) with least square method to obtain coefficients {h(p,q)j}(p,q)∈Λ.Step
3. Choose a point x in the domain needed to be transformed to new coordinate, and calculate the coefficients {Q(p,q)(x)}(p,q)∈Λ with the formula (41).Step
4. Obtain the approximation uj(x) using (41).


### 4.1. Computation of the Matrix

In the Step  1 and Step  3 above, the main cost is the computation of(43)∫Ry2‍f(x,y)ϕpj(y1)ϕqj(y2)dy.
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(44)∫‍ϕ(t)f(t)dt=∑k=0nAkf(tk)
holds for any polynomial of degree less than or equal to 2n+1, then it is called Gauss-type quadrature with scale function ϕ(t) as its weight function, the dots {tk}, and the coefficients {Ak} are called generalized Gauss-type quadrature dots and generalized Gauss-type quadrature weights.Similar to the classical Gauss-type quadrature formula, we recall that[0,2N-1] is the support of the scale function ϕ; therefore, we have the following.Proposition 8.
In formula (44), the necessary and sufficient condition for {tk} being generalized Gauss-type quadrature dots is wn(t)=Πk=0n(t-tk) being orthogonal to any polynomial Pn(t)(degree≤n) with ϕ(t) as the weight function; that is,
(45)∫02N-1‍ϕ(t)Pn(t)wn(t)dt=0.Proposition 9.
In formula (44), if f(t)∈C2n+2[0,2N-1], the error of (44) is
(46)R=∫02N-1ϕ(t)f(t)dt-∑k=0nAkf(tk)=f(2n+2)(ξ)ϕ(ξ)(2n+2)!∫02N-1wn2(t)dt,
where ξ∈[0,2N-1].A significant task is to determine{tk} and {Ak} in formula (44). Let f(t)=ti, 0≤i≤2n+1, we obtain a nonlinear system [25]:
(47)∑k=0nAktki=∫‍ϕ(t)tidt,i=0,…,2n+1.
First, we compute the right-hand side of (47) recursively. Suppose that 2n+1≤N-1, from the theory of wavelet, we can write
(48)tii!=∑k=-∞+∞Pi(k)ϕ(t-k),i=0,…,2n+1,
from which we know
(49)cm=∫-∞+∞‍tmm!ϕ(t)dt=Pm(0),
and it can be computed recursively using [26]
(50)c0=1ci=12i-1∑r=1iMrci-r
with Mr=(1/2)Σmhm(mr/r!), where reals hk are coefficients in ϕ(t)=2Σhkϕ(2t-k). So far, the right-hand side in (47) is computed. Then, we solve the nonlinear system (47), so that 2n+2 coefficients {tk} and {Ak} are obtained. For example, when n=1, N=3,4,5,6, we have
(51)A0≈0,A1≈1
and obtain Table 1.Table 1
N
x
1
=
c
3
8.17401
E
-
001
4
1.00539
5
1.19390
6
1.38216Therefore, we have a very simple formula(52)∫02N-1‍ϕ(t)f(t)dt≈f(c).
For f∈C4, by use of spline function, we can prove that its error, ɛ, is
(53)ɛ=O(h4),h=2N-12j.
Based on these results, we know that
(54)12jf(p+c2j,q+c2j)
can be used to approximate
(55)∫Ry2‍f(x,y)ϕpj(y1)ϕqj(y2)dy,
and we will use it in our experiments.Algorithm 10 (boundary measure and wavelet (detail)).
Input
{
x
¯
k
}
⊂
Γ and {f(x¯k)}.
Step
1
Choose a resolution levelj and fix a rectangular domain D covering {x¯k}.
Calculate the index setΛ={(p,q)} by ϕ(p,q)j∩D≠{0}.
Choose any|Λ| collocation points {xk}⊂{x¯k}.
Step
2
Fork=1 to |Λ| Do
Compute matrix elements(56)Kk,(p,q)(xk)=∫Ry2ϕpj(y1)ϕqj(y2)log⁡|xk-y|dy≈12log⁡[(xk1-p+c2j)2+(xk2-q+c2j)2].
Obtain right-hand side{f(xk)} from {f(x¯k)}.
End.
Step
3
Solve the linear systemsKh¯=f¯ and obtain h¯={h(p,q)j}.
Step
4
Choose the points{zk}k=1M in the domain to be transformed to new coordinate space.
Step
5
Fork=1 to M Do:
Evaluate(57)Q(p,q)(zk)=∫Ry2ϕpj(y1)ϕqj(y2)log⁡|zk-y|dy≈12log⁡[(zk1-p+c2j)2+(zk2-q+c2j)2].
Calculate the new coordinate(58)uj(zk)=∑(p,q)∈Λh(p,q)jQ(p,q)(zk).
End.
Output
{
u
j
(
z
k
)
}
k
=
1
M.


## 4.1. Computation of the Matrix

In the Step  1 and Step  3 above, the main cost is the computation of(43)∫Ry2‍f(x,y)ϕpj(y1)ϕqj(y2)dy.
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(44)∫‍ϕ(t)f(t)dt=∑k=0nAkf(tk)
holds for any polynomial of degree less than or equal to 2n+1, then it is called Gauss-type quadrature with scale function ϕ(t) as its weight function, the dots {tk}, and the coefficients {Ak} are called generalized Gauss-type quadrature dots and generalized Gauss-type quadrature weights.Similar to the classical Gauss-type quadrature formula, we recall that[0,2N-1] is the support of the scale function ϕ; therefore, we have the following.Proposition 8.
In formula (44), the necessary and sufficient condition for {tk} being generalized Gauss-type quadrature dots is wn(t)=Πk=0n(t-tk) being orthogonal to any polynomial Pn(t)(degree≤n) with ϕ(t) as the weight function; that is,
(45)∫02N-1‍ϕ(t)Pn(t)wn(t)dt=0.Proposition 9.
In formula (44), if f(t)∈C2n+2[0,2N-1], the error of (44) is
(46)R=∫02N-1ϕ(t)f(t)dt-∑k=0nAkf(tk)=f(2n+2)(ξ)ϕ(ξ)(2n+2)!∫02N-1wn2(t)dt,
where ξ∈[0,2N-1].A significant task is to determine{tk} and {Ak} in formula (44). Let f(t)=ti, 0≤i≤2n+1, we obtain a nonlinear system [25]:
(47)∑k=0nAktki=∫‍ϕ(t)tidt,i=0,…,2n+1.
First, we compute the right-hand side of (47) recursively. Suppose that 2n+1≤N-1, from the theory of wavelet, we can write
(48)tii!=∑k=-∞+∞Pi(k)ϕ(t-k),i=0,…,2n+1,
from which we know
(49)cm=∫-∞+∞‍tmm!ϕ(t)dt=Pm(0),
and it can be computed recursively using [26]
(50)c0=1ci=12i-1∑r=1iMrci-r
with Mr=(1/2)Σmhm(mr/r!), where reals hk are coefficients in ϕ(t)=2Σhkϕ(2t-k). So far, the right-hand side in (47) is computed. Then, we solve the nonlinear system (47), so that 2n+2 coefficients {tk} and {Ak} are obtained. For example, when n=1, N=3,4,5,6, we have
(51)A0≈0,A1≈1
and obtain Table 1.Table 1
N
x
1
=
c
3
8.17401
E
-
001
4
1.00539
5
1.19390
6
1.38216Therefore, we have a very simple formula(52)∫02N-1‍ϕ(t)f(t)dt≈f(c).
For f∈C4, by use of spline function, we can prove that its error, ɛ, is
(53)ɛ=O(h4),h=2N-12j.
Based on these results, we know that
(54)12jf(p+c2j,q+c2j)
can be used to approximate
(55)∫Ry2‍f(x,y)ϕpj(y1)ϕqj(y2)dy,
and we will use it in our experiments.Algorithm 10 (boundary measure and wavelet (detail)).
Input
{
x
¯
k
}
⊂
Γ and {f(x¯k)}.
Step
1
Choose a resolution levelj and fix a rectangular domain D covering {x¯k}.
Calculate the index setΛ={(p,q)} by ϕ(p,q)j∩D≠{0}.
Choose any|Λ| collocation points {xk}⊂{x¯k}.
Step
2
Fork=1 to |Λ| Do
Compute matrix elements(56)Kk,(p,q)(xk)=∫Ry2ϕpj(y1)ϕqj(y2)log⁡|xk-y|dy≈12log⁡[(xk1-p+c2j)2+(xk2-q+c2j)2].
Obtain right-hand side{f(xk)} from {f(x¯k)}.
End.
Step
3
Solve the linear systemsKh¯=f¯ and obtain h¯={h(p,q)j}.
Step
4
Choose the points{zk}k=1M in the domain to be transformed to new coordinate space.
Step
5
Fork=1 to M Do:
Evaluate(57)Q(p,q)(zk)=∫Ry2ϕpj(y1)ϕqj(y2)log⁡|zk-y|dy≈12log⁡[(zk1-p+c2j)2+(zk2-q+c2j)2].
Calculate the new coordinate(58)uj(zk)=∑(p,q)∈Λh(p,q)jQ(p,q)(zk).
End.
Output
{
u
j
(
z
k
)
}
k
=
1
M.


## 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 transformationT:(x,y)→(u,v) satisfies (5), which is recalled as follows:
(59)Δu(x1,x2)=0,(x1,x2)∈Ω,u|Γ=f(x1,x2),(x1,x2)∈Γ,Δv(x1,x2)=0,(x1,x2)∈Ω,v|Γ=g(x1,x2),(x1,x2)∈Γ.
In this experiment, the equations of the boundary conditions u|Γ=f(x1,x2), v|Γ=g(x1,x2) and the region Ω are specified by the following forms:
(60)f(x1,x2)=1.2x1,g(x1,x2)=x2+(x1-0.5)2,(61)Ω:0≤x1,x2≤0.5.
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
(62)u=u(x1,x2),v=v(x1,x2).
The primary idea is to find a mathematical model for the geometric distortion process, specifically, the equations u=u(x1,x2),v=v(x1,x2) 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 calledtiepoints. 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 Figure7 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 Figure8.
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:
(63)u=k1x1+k2x2+k3x1x2+k4,v=k5x1+k6x2+k7x1x2+k8,
where ki, i=1,2,…,8, are constants to be determined by solving the eight simultaneous equations. Because we have defined four tiepoints, thus we have eight equations, where x1, x2, u, and v are known. Now we can solve the eight equations for the eight unknowns ki, so that the necessary equations for the coordinate mapping can be obtained. The mapping equations u=u(x1,x2), v=v(x1,x2) then are applied to all (x1,x2) pairs needed. In practice, the values of x1, x2, u, and v 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 thex-axis and y-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.Experiment1: nonlinear harmonic and its inverse: (a) original image, (b) distorted image, (c) restored image by IEWC, and (d) restored image by bilinear method.
(a)
(b)
(c)
(d)Figure 8
A harmonic distorted image is approximated by piecewise bilinear model.Comparison in Experiment1: projection of images to the x-axis: (a) original image, (b) distorted image, (c) restored image by IEWC, and (d) restored image by the piecewise bilinear method.
(a)
Origx-axis
(b)
Changx-axis
(c)
Pdex-axis
(d)
Bi-128x-axisComparison in Experiment1: projection of images to the y-axis: (a) original image, (b) distorted image by harmonic transformation, (c) restored image by IEWC, and (d) restored image by the piecewise bilinear method.
(a)
Origy-axis
(b)
Changy-axis
(c)
Pdey-axis
(d)
Bi-128y-axisExperiment 2 (nonlinear-harmonic and its inverse for a circle area).
Consider(64)T2:(x1,x2)⟶(u,v),Ω:x12+x22≤0.42,f(x1,x2)=x1+x2,g(x1,x2)=x12-x22.
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.Experiment2: nonlinear harmonic and its inverse for a circle area.
(a)
(b)
(c)
(d)Comparison in Experiment2: projection of images to the x-axis: (a) original image, (b) distorted image, (c) restored image by IEWC, and (d) restored image by the piecewise bilinear method.
(a)
Origx-axis
(b)
Changx-axis
(c)
Pdex-axis
(d)
Bi-128x-axisComparison in Experiment2: projection of images to the y-axis: (a) original image, (b) distorted image, (c) restored image by IEWC, and (d) restored image by the piecewise bilinear method.
(a)
Origy-axis
(b)
Changy-axis
(c)
Pdey-axis
(d)
Bi-128y-axis


## 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 calledcurvelet 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.



---

*Source: 798080-2014-04-01.xml*