# Shannon Wavelet Precision Integration Method for Pathologic Onion Image Segmentation Based on Homotopy Perturbation Technology

Authors: Haihua Wang; Shu-Li Mei
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

Image segmentation variational method is good at processing the images with blurry and complicated contours, which is useful in quality identification of pathologic picture of onion. An adaptive Shannon wavelet precise integration method (WPIM) on digital image segmentation was proposed based on the image processing variational model to improve the processing speed and eliminate the artifacts of the images. First, taking full advantage of the interpolation property of the Shannon wavelet function, a multiscale Shannon wavelet interpolation scheme was constructed based on the homotopy perturbation method (HPM). The image pixels of the Burkholderia cepacia (ex-Burkholder) infected onions were taken as the collocation points of the WPIM. Then, with this scheme, the image segmentation model (C-V model) can be discretized into a system of nonlinear ODEs and solved by the half-analytical scheme combining the HPM and the precision integration method. At last, the numerical precision and efficiency of WPIM were discussed and compared with other common segmentation methods such as OSTU method and Sobel operator. The results show that the contour curve of the segmentation object obtained by the new method has many excellent properties such as closed and clear topological structure and the artifacts can be eliminated.


---

## Body

## 1. Introduction

Rot is the primary factor of the onion storage losses. The damages from the pests, soak and over-nitrogen often lead to the onion infected with pathogen. Cross-infection of onions accelerates rotting in the package. Therefore, grading and classification of the onions are necessary in the postharvest processing to decrease the losses of onions.The manual identification is time-consuming and unreliable. Image measurement technology [1] is the new method for grading the quality of agricultural products. Image segmentation is an important part of the image measurement. Accurate segmentation techniques can be used for rotting detection and grading of onions. The classical segementation methods such as Sobel, Canny, quadtree, and OTSU algorithm, take the gradient of the image as the feature descriptor directly in image segmentation. But the object boundary and target pixels obtained by these methods are often unclosed, which is makes it difficult to analyze geometric characteristics of the target for connection and econometric analysis of segmentation results. In the recent years, the wavelet precise integration method [2–10] has been developed to solve the nonlinear PDEs for image processing, which can improve the efficiency and precision of the image processing effectively. WPIM is an image segmentation variational method based on the C-V model, with which the segmentation results with closed object contour can be obtained. In addition, the continuity of wavelet function is helpful to eliminate the artifact introduced by the difference method. There is no doubt that these merits are necessary in the image measurement and product quality grading.In the recent decades, many wavelets which have compact support, smoothness, and other properties have been constructed. Cattani studied the properties of the Shannon wavelet function, which possesses many advantages such as orthogonality, continuity, and differentiability [11]. It also has the advantage over the Hermite DAF in that it is an interpolating function, producing matrix equations that have the potential to be relatively sparse. In addition, the second order approximation of a C2-function, based on Shannon wavelet functions, is given [12]. The approximation is compared with the wavelet reconstruction formula and the error of approximation is explicitly computed [13, 14]. Furthermore, Shannon wavelet has been used to solve the fractional calculus problems in the recent years [15–17]. A perceived disadvantage of the Shannon scaling function is that it tends to zero quite slowly as |x|→∞. A direct consequence of this is that, when calculating the derivatives, a large number of the nodal values will contribute significantly. However, comparing with other wavelets, the Shannon scaling function is one that possesses more excellent numerical properties such as the interpolation, orthogonality, and smoothness. Therefore, Shannon wavelet is employed to construct the multiscale interpolator in our research.The purpose of this research is to construct a multiscale two-dimension wavelet interpolation operator based on the homotopy perturbation method [18–20], with which the image segmentation PDE (C-V model) can be discretized into a system of nonlinear ODEs. Combining the precise integration method with homotopy perturbation method [21, 22] for solving nonlinear problems, a fast wavelet numerical algorithm for P-M and C-V model in image processing can be obtained. With this new image segmentation method, the pathologic onion classification and grading can be processed precisely and efficiency.


## 2. Chan-Vese Model

In order to solve Mumford-Shah model with Euler-Lagrange method, a simplified model was deduced by Chan and Vese, in which Euclid length was employed instead of Hausdorff length. This simplified model can also be called Chan-Vese model, which can be expressed as follows:(1)ECV(c1,c2,C)=λ1∫Ω1(I0-c1)2dxdy+λ2∫Ω2(I0-c2)2dxdy+ν|C|,(2)ci=meanΩi(u0)=∫Ωiu0(x,y)dxdyArea(Ωi),i=1,2,
where λ1 and λ2 are positive constants and c1 and c2 are the average gray level values inside (Ω1) and outside (Ω2) of the object contour, respectively. I0 denotes the image to process, |C| is the length of the object contour, and ν is the weight parameter. According to the level set method, the contour curve of the objects should be embedded into the level set function as follows:
(3)C={(x,y)∣(x,y)∈Ω,ϕ(x,y)=0},Ω1={(x,y)∣(x,y)∈Ω,ϕ(x,y)>0},Ω2={(x,y)∣(x,y)∈Ω,ϕ(x,y)<0}.
Then, the level set-based C-V model can be rewritten as follows:
(4)E(c1,c2,ϕ)=λ1∫Ω|I0-c1|2H(ϕ)dxdy+λ2∫Ω|I0-c2|2(1-H(ϕ))dxdy+ν∫Ω|H(ϕ)|dxdy,H(ϕ)={1,ϕ≥00,ϕ<0,δε=επ(ε2+ϕ2).
Using the variational method, the PDEs with respect to the variable ϕ can be obtained as follows:
(5)∂ϕ∂t=δε(ϕ)[νdiv⁡(∇ϕ|∇ϕ|)-λ1|I0-c1|2+λ2|I0-c2|2].
Obviously, div⁡(∇ϕ/|∇ϕ|) is the curvature of the level set function ϕ, and δε(ϕ) is used to constrain the growth of the level set function.The solution of (5) is the level set function ϕ(x,y,t) at time t. The zero level set is the object contour curve, which can be obtained by solving the equation ϕ(x,y,t)=0. We focus on multiscale wavelet image segmentation discrete format and the corresponding numerical solution method.


## 3. HPM-Based Wavelet Interpolation Operator Construction Schemes

Let the definition domain of the image be(xmin⁡,xmax⁡)×(ymin⁡,ymax⁡); the discretization points can be defined as (xk1j,yk2j), where j is a scale parameter and k1 and k2 are position parameters. So
(6)xk1j=xmin⁡+k1xmax⁡-xmin⁡2j,yk1j=ymin⁡+k2ymax⁡-ymin⁡2j,j,k1,k2∈ℤ.In addition,wk1,k2j(m,n)(x,y) denotes the multiscale wavelet function and the corresponding mth and nth derivatives with respect to x and y, respectively. The level set function ϕ(x,y,t) and the corresponding derivative function can be discretized as follows:
(7)ϕJ(m,n)(x,y,t)=∑k01=01∑k02=01ϕ(xk010,yk020)wk01,k020(m,n)(x,y)+∑j=0J-1∑k11=02j-1∑k12=02j-1[w2k11+1,2k12+1j+1(m,n)(x,y)αj,k11,k121(t)w2k11+1,2k12j+1(m,n)(x,y)ffffffffffffff+αj,k11,k122(t)w2k11,2k12+1j+1(m,n)(x,y)ffffffffffffff+αj,k11,k123(t)w2k11+1,2k12+1j+1(m,n)(x,y)],
where j and J are constants, which denote the wavelet scale number and the maximum of the scale number, respectively. αj,k11,k121, αj,k11,k122, and αj,k11,k123 are the wavelet coefficients at the points (xk1j,yk2j). According to the interpolation wavelet transform theory, the wavelet coefficients can be written as
(8)αj,k1,k21=ϕ(xj+1,2k1+1,yj+1,2k2)-Ijϕ(xj+1,2k1+1,yj+1,2k2),αj,k1,k22=ϕ(xj+1,2k1,yj+1,2k2+1)-Ijϕ(xj+1,2k1,yj+1,2k2+1),αj,k1,k23=ϕ(xj+1,2k1+1,yj+1,2k2+1)-Ijϕ(xj+1,2k1+1,yj+1,2k2+1),
where Ij denotes the multilevel interpolation operator. In order to obtain the multilevel interpolation operator, it is necessary to express the wavelet coefficients αj,k1,k21,αj,k1,k22,αj,k1,k23 as a weighted sum of u in all of the collocation points in the J-level. Therefore, we should give the definition of the restriction operator as follows:
(9)Rk1,k2,m1,m2l,l,j,j={1,xk1l=xm1j,yk2l=ym2j0,otherwise.
Using the restriction operator, u(x2k1+1j+1,y2k2j+1),u(x2k1j+1,y2k2+1j+1) and u(x2k1+1j+1,y2k2+1j+1) can be rewritten as
(10)ϕ(x2k1+1j+1,y2k2j+1)=∑n1=02J∑n2=02JR2k1+1,2k2,n1,n2j+1,j+1,J,Jϕ(xn1J,yn2J),ϕ(x2k1j+1,y2k2+1j+1)=∑n1=02J∑n2=02JR2k1,2k2+1,n1,n2j+1,j+1,J,Jϕ(xn1J,yn2J),ϕ(x2k1+1j+1,y2k2+1j+1)=∑n1=02J∑n2=02JR2k1+1,2k2+1,n1,n2j+1,j+1,J,Jϕ(xn1J,yn2J).
Introducing the extension operators C1, C2, and C3 and substituting (10) into (8), the wavelet coefficients can be rewritten as(11)αj,k1,k21=∑n1=02J∑n2=02JR2k1+1,2k2,n1,n2j+1,j+1,J,Jϕ(xn1J,yn2J)-[∑n1=02J∑n2=02J∑k01=02j0∑k02=02j0Rk01,k02,n1,n2j0,j0,J,Jϕ(xn1J,yn2J)wk01,k02j0(x2k1+1j+1,y2k2j+1)ddddd+∑j1=j0j-1∑n1=02J∑n2=02J∑k11=02j1∑k12=02j1(C1k11,k12,n1,n2j1,j1,J,Jw2k11+1,2k12j1+1(x2k1+1j+1,y2k2j+1)ϕ(xn1J,yn2J)ddddddddddddddddddddddddd+C2k11,k12,n1,n2j1,j1,J,Jw2k11,2k12+1j1+1(x2k1+1j+1,y2k2j+1)u(xn1J,yn2J)ddddddddddddddddddddddddd+C3k11,k12,n1,n2j1,j1,J,Jw2k11+1,2k12+1j1+1(x2k1+1j+1,y2k2j+1)ϕ(xn1J,yn2J))∑n2=02J∑k01=02j0∑k02=02j0Rk01,k02,n1,n2j0,j0,J,Jϕ(xn1J,yn2J)wk01,k02j0(x2k1+1j+1,y2k2j+1)]=∑n1=02J∑n2=02JC1k1,k2,n1,n2j,j,J,Jϕ(xn1J,yn2J).αj,k1,k22 and αj,k1,k23 are similar to αj,k1,k21. From the above equation, the extension operator can be obtained as(12)C1k1,k2,n1,n2j,j,J,J=R2k1+1,2k2,n1,n2j+1,j+1,J,J-[∑k01=02j0∑k02=02j0Rk01,k02,n1,n2j0,j0,J,Jϕ(xn1J,yn2J)wk01,k02j0(x2k1+1j+1,y2k2j+1)ddddd+∑j1=j0j-1∑n2=02J∑k11=02j1(C1k11,k12,n1,n2j1,j1,J,Jw2k11+1,2k12j1+1(x2k1+1j+1,y2k2j+1)ϕ(xn1J,yn2J)dddddddddddddddddd+C2k11,k12,n1,n2j1,j1,J,Jw2k11,2k12+1j1+1(x2k1+1j+1,y2k2j+1)ϕ(xn1J,yn2J)dddddddddddddddddd+C3k11,k12,n1,n2j1,j1,J,Jw2k11+1,2k12+1j1+1(x2k1+1j+1,y2k2j+1)ϕ(xn1J,yn2J))∑k01=02j0∑k02=02j0Rk01,k02,n1,n2j0,j0,J,Jϕ(xn1J,yn2J)wk01,k02j0(x2k1+1j+1,y2k2j+1)].C2 and C3 can be obtained with the same method. Therefore, the calculation time complexity of the wavelet transform coefficients αj,k11,k121, αj,k11,k122, and αj,k11,k123 is O((1/3)42J-1).Substitutingαj,k11,k121, αj,k11,k122, and αj,k11,k123, C1, C2, and C3 into (2), the multilevel wavelet interpolation operator can be obtained as
(13)In1,n2(x,y)=∑k01=02j0∑k02=02j0Rk01,k02,n1,n2j0,j0,J,Jwk01,k02j0(x,y)+∑j=j0J-1∑k1=02j∑k2=02j(C1k1,k2,n1,n2j,j,J,Jw2k1+1,2k2j+1(x,y)ϕ(xn1J,yn2J)ffffffffff+C2k1,k2,n1,n2j,j,J,Jw2k1,2k2+1j+1(x,y)ϕ(xn1J,yn2J)fffffffifff+C3k1,k2,n1,n2j,j,J,Jw2k1+1,2k2+1j+1(x,y)ϕ(xn1J,yn2J)).
Then, (7) can be rewritten as
(14)ϕJ(m,n)(x,y,t)=∑n12J∑n22JIn1,n2(x,y)ϕ(xn1J,yn2J).
Substituting (14) into (5), the multilevel wavelet discretization scheme of PERONA-MALIK model can be obtained.The purpose of constructing the multilevel wavelet collocation method is to decrease the amount of the collocation points and then improve the efficiency of the algorithm. But the efficiency will be eliminated if the computation complexity of the multilevel wavelet interpolation operator is too high. It is easy to understand that the interpolation wavelet coefficient is the error between the interpolation result and the exact result at the same collocation point. And so the wavelet coefficient must be the function of the parametert. In other words, the wavelet coefficient should vary with the time parameter t. Then, the interpolation operator can be viewed as a nonlinear problem. HPM is efficient and effective tool to solve nonlinear problem. Aiming to improve the efficiency of the multilevel wavelet interpolation operators, HPM would be employed to construct a novel interpolation operator in this section.For convenience,ϕ and its derivative in (5) should be rewritten as
(15)∂ϕ∂t=F(t,x,y,ϕ,∂ϕ∂x,∂ϕ∂y,∂2ϕ∂x2,∂2ϕ∂x∂y,∂2ϕ∂y2)(t>0)ϕ(x,y,0)=ϕ0(x,y),dϕJ(x,y,t)dt=F[t,x,y,ϕJ(x,y,t),ϕJ(1,0)(x,y,t),fffddffddddffϕJ(0,1)(x,y,t),ϕJ(2,0)(x,y,t),fffddffddddffϕJ(1,1)(x,y,t),ϕJ(0,2)(x,y,t)],
respectively, where
(16)ϕJ(x,y,t)=∑k01=01∑k02=01ϕ(xk010,yk020)wk01,k020(x,y)+∑j=0J-1∑k11=02j-1∑k12=02j-1[αj,k11,k121w2k11+1,2k12j+1(x,y)ddddddddddddiddd+αj,k11,k122w2k11,2k12+1j+1(x,y)dddddddddddddddd+αj,k11,k123w2k11+1,2k12+1j+1(x,y)].
The value of ϕJ(x,y,tn) at tn is denoted as ϕn, and
(17)F[tn,x,y,ϕJ(x,y,tn),ϕJ(1,0)(x,y,tn),ϕJ(0,1)(x,y,tn),fϕJ(2,0)(x,y,tn),ϕJ(1,1)(x,y,tn),ϕJ(0,2)(x,y,tn)]
is denoted as Fn. And then a linear homotopy function can be constructed as
(18)ϕJ(x,y,t)=(1-ε)Fn+εFn+1.
It is easy to identify the homotopy parameter as
(19)ε(t)=t-tntn+1-tnt∈[tn,tn+1]∴ε∈[0,1].
According to the perturbation theory, the solution of (18) can be expressed as the power series expansion of ε as follows:
(20)ϕJ=ϕ0J+εϕ1J+ε2ϕ2J+⋯.
Substituting (20) into (18) and rearranging based on powers of ε-terms, we have
(21)ε0:ϕ0J=Fn,ε1:ϕ1J=Fn+1-Fn,⋮
According to HPM, we obtain the wavelet coefficients αj,k1,k21(tn+1), αj,k1,k22(tn+1), αj,k1,k23(tn+1) at tn as follows:(22)αj,k1,k21=ϕ(x2k1+1j+1,y2k2j+1)-Ijϕ(x2k1+1j+1,y2k2j+1)=ϕ(x2k1+1j+1,y2k2j+1)-[∑k01=01∑k02=01ϕ(xk010,yk020)wk01,k020(x2k1+1j+1,y2k2j+1)ddddd+∑j1=0j-1∑k11=02j1∑k12=02j1(αj1,k11,k121w2k11+1,2k12j1+1(x2k1+1j+1,y2k2j+1)dddgggggggggiggdd+αj1,k11,k122w2k11,2k12+1j1+1(x2k1+1j+1,y2k2j+1)ddgggggggggggddd+αj1,k11,k123w2k11+1,2k12+1j1+1(x2k1+1j+1,y2k2j+1))∑k01=01∑k02=01ϕ(xk010,yk020)],αj,k1,k22=ϕ(x2k1j+1,y2k2+1j+1)-Ijϕ(x2k1j+1,y2k2+1j+1)=ϕ(x2k1j+1,y2k2+1j+1)-[∑k01=01∑k02=01ϕ(xk010,yk020)wk01,k020(x2k1j+1,y2k2+1j+1)ddddd+∑j1=0j-1∑k11=02j1∑k12=02j1(αj1,k11,k121w2k11+1,2k12j1+1(x2k1j+1,y2k2+1j+1)ddggggggggggddd+αj1,k11,k122w2k11,2k12+1j1+1(x2k1j+1,y2k2+1j+1)dddggggggggggdd+αj1,k11,k123w2k11+1,2k12+1j1+1(x2k1j+1,y2k2+1j+1))∑k01=01∑k02=01ϕ(xk010,yk020)],αj,k1,k23=ϕ(x2k1+1j+1,y2k2+1j+1)-Ijϕ(x2k1+1j+1,y2k2+1j+1)=ϕ(x2k1+1j+1,y2k2+1j+1)-[∑k01=01∑k02=01ϕ(xk010,yk020)wk01,k020(x2k1+1j+1,y2k2+1j+1)ddddf+∑j1=0j-1∑k11=02j1∑k12=02j1(αj1,k11,k121w2k11+1,2k12j1+1(x2k1+1j+1,y2k2+1j+1)ffddddddddddddddff+αj1,k11,k122w2k11,2k12+1j1+1(x2k1+1j+1,y2k2+1j+1)fffdddddddddddddff+αj1,k11,k123w2k11+1,2k12+1j1+1(x2k1+1j+1,y2k2+1j+1))∑k01=01∑k02=01ϕ(xk010,yk020)].Obviously, the calculation time complexity of the wavelet transform coefficients αj,k1,k21, αj,k1,k22, and αj,k1,k23 is O(4J), which is decreased greatly than that in (8) which is O((1/3)42J-1).Substituting the wavelet transform efficient (22) into (16), we obtain
(23)ϕJ(x,y,tn+1)=ϕJ(x,y,tn)+Δt2[ϕ0J(0,2)(x,y,tn+1))F(tn,x,y,ϕJ(x,y,tn),ϕJ(1,0)(x,y,tn),ddddgfddϕJ(0,1)(x,y,tn),ϕJ(2,0)(x,y,tn),ddddffddϕJ(1,1)(x,y,tn),ϕJ(0,2)(x,y,tn))dddffgdd+F(tn+1,x,y,ϕ0J(x,y,tn+1),dddfgggffddϕ0J(1,0)(x,y,tn+1),ϕ0J(0,1)(x,y,tn+1),dddfgggffddϕ0J(2,0)(x,y,tn+1),ϕ0J(1,1)(x,y,tn+1),dddggfgffddϕ0J(0,2)(x,y,tn+1))]
and the derivative function
(24)ϕJ(m,n)(x,y)=∑k01=01∑k02=01ϕ(xk010,yk020)wk01,k020(m,n)(x,y)+∑j=0J-1∑k11=02j-1∑k12=02j-1[αj,k11,k121w2k11+1,2k12j+1(m,n)(x,y)ddddddddddd+αj,k11,k122w2k11,2k12+1j+1(m,n)(x,y)ddddddddddd+αj,k11,k123w2k11+1,2k12+1j+1(m,n)(x,y)].
Obviously, the computation complexity is decreased greatly comparing with (14).


### 3.1. The Multiscale Interpolation Wavelet Approximation of the C-V Model

There are many ways to solve partial differential equations and the most typical method is the difference method. This method uses the flat function to describe image, approximately the surface function. But, it is easy to cause artifacts phenomenon, affecting the accuracy of image segmentation. Wavelet function has both smooth and compactly supported characteristics. Besides, performance of multiscale analysis can be used to construct the multiscale adaptive interpolation operator for solving nonlinear partial differential equations. The wavelet approximation of the level set function and its derivative with respect tox and y, respectively, can be expressed as follows:
(25)ϕJ(m,n)(x,y,t)=∑k01=01∑k02=01ϕ(xk010,yk020)wk01,k020(m,n)(x,y)+∑j=0J-1∑k11=02j-1∑k12=02j-1[αj,k11,k121(t)w2k11+1,2k12j+1(m,n)(x,y)dddddddddddfd+αj,k11,k122(t)w2k11,2k12+1j+1(m,n)(x,y)ddddddddddddd+αj,k11,k123(t)w2k11+1,2k12+1j+1(m,n)(x,y)],
where wk1,k2j(m,n)(x,y) is the wavelet function and its m- and n-order derivative with respect to x and y, respectively. αj,k11,k121(t), αj,k11,k122(t), and αj,k11,k123(t) are wavelet transform coefficients. We convert (2) to wavelet multiscale discrete format of the level set function by stead of (25):
(26)dϕJ(x,y,t)dt=δε(ϕJ(x,y,t))×[νdiv⁡(∇ϕJ(x,y,t)|∇ϕJ(x,y,t)|)-λ1|I0-c1|2+λ2|I0-c2|2].It is the most direct way for dynamic adaptation to only retain the distribution points of corresponding wavelet coefficients that satisfy the condition(27)min⁡(|αj,k11,k121(t)|,|αj,k11,k122(t)|,|αj,k11,k123(t)|)≥ε.Time domain numerical integration of partial differential equations is an iterative process; therefore, some points, which are possible important next step, need to be kept to enable the algorithm to track singularities of solutions. So adjacent points of distribution points also should be kept. The adjacent region can be delineated as follows:(28)|s-j|≤M,|xks-xij|≤εx,|xks-xij|≤εy,
where s, j are numbers of different scale wavelet and k,i,M∈Z, εx,εy are constant.


### 3.2. Nonlinear Discrete Ordinary Differential Equations

Becauseν is a small parameter in ordinary differential equation (26), value of νdiv⁡(∇ϕJ(x,y,t)/|∇ϕJ(x,y,t)|) is low. Equation (26) can be converted into
(29)dϕJ(x,y,t)dt=δε(ϕJ(x,y,t))×[νdiv⁡(∇ϕJ(x,y,t0)|∇ϕJ(x,y,t0)|)-λ1|I0-c1|2+λ2|I0-c2|2].
The solution of (29) is
(30)ϕJ(x,y,t)=(εm2π)1/3+(m2πε5)-1/3,
where
(31)m=νdiv⁡(∇ϕJ(x,y,t0)|∇ϕJ(x,y,t0)|)-λ1|I0-c1|2+λ2|I0-c2|2.We can get the solution of ordinary differential equations (26) by iterative solution (30).


## 3.1. The Multiscale Interpolation Wavelet Approximation of the C-V Model

There are many ways to solve partial differential equations and the most typical method is the difference method. This method uses the flat function to describe image, approximately the surface function. But, it is easy to cause artifacts phenomenon, affecting the accuracy of image segmentation. Wavelet function has both smooth and compactly supported characteristics. Besides, performance of multiscale analysis can be used to construct the multiscale adaptive interpolation operator for solving nonlinear partial differential equations. The wavelet approximation of the level set function and its derivative with respect tox and y, respectively, can be expressed as follows:
(25)ϕJ(m,n)(x,y,t)=∑k01=01∑k02=01ϕ(xk010,yk020)wk01,k020(m,n)(x,y)+∑j=0J-1∑k11=02j-1∑k12=02j-1[αj,k11,k121(t)w2k11+1,2k12j+1(m,n)(x,y)dddddddddddfd+αj,k11,k122(t)w2k11,2k12+1j+1(m,n)(x,y)ddddddddddddd+αj,k11,k123(t)w2k11+1,2k12+1j+1(m,n)(x,y)],
where wk1,k2j(m,n)(x,y) is the wavelet function and its m- and n-order derivative with respect to x and y, respectively. αj,k11,k121(t), αj,k11,k122(t), and αj,k11,k123(t) are wavelet transform coefficients. We convert (2) to wavelet multiscale discrete format of the level set function by stead of (25):
(26)dϕJ(x,y,t)dt=δε(ϕJ(x,y,t))×[νdiv⁡(∇ϕJ(x,y,t)|∇ϕJ(x,y,t)|)-λ1|I0-c1|2+λ2|I0-c2|2].It is the most direct way for dynamic adaptation to only retain the distribution points of corresponding wavelet coefficients that satisfy the condition(27)min⁡(|αj,k11,k121(t)|,|αj,k11,k122(t)|,|αj,k11,k123(t)|)≥ε.Time domain numerical integration of partial differential equations is an iterative process; therefore, some points, which are possible important next step, need to be kept to enable the algorithm to track singularities of solutions. So adjacent points of distribution points also should be kept. The adjacent region can be delineated as follows:(28)|s-j|≤M,|xks-xij|≤εx,|xks-xij|≤εy,
where s, j are numbers of different scale wavelet and k,i,M∈Z, εx,εy are constant.


## 3.2. Nonlinear Discrete Ordinary Differential Equations

Becauseν is a small parameter in ordinary differential equation (26), value of νdiv⁡(∇ϕJ(x,y,t)/|∇ϕJ(x,y,t)|) is low. Equation (26) can be converted into
(29)dϕJ(x,y,t)dt=δε(ϕJ(x,y,t))×[νdiv⁡(∇ϕJ(x,y,t0)|∇ϕJ(x,y,t0)|)-λ1|I0-c1|2+λ2|I0-c2|2].
The solution of (29) is
(30)ϕJ(x,y,t)=(εm2π)1/3+(m2πε5)-1/3,
where
(31)m=νdiv⁡(∇ϕJ(x,y,t0)|∇ϕJ(x,y,t0)|)-λ1|I0-c1|2+λ2|I0-c2|2.We can get the solution of ordinary differential equations (26) by iterative solution (30).


## 4. Experiment and Discussion of Onion Infected Region Segmentation

Figure1(a) is a 256×302 image of an onion infected by sour skin virus. We noticed that the onion has a water-soaked appearance. Compared with the background, the grayscale difference between the water-soaked appearance and the healthy part is smaller. So it is beneficial to compare performance of the different algorithms. Figure 1(b) is an ideal segmentation results. Segmentation target is infected regions of onion; however infected part is often not uniform. In the image, the difference of gradient is less than 1 at upper left part due to slight infection, so the algorithm is difficult to precise segmentation based on global threshold. Therefore, the best one of the different segmentation algorithms can distinguish the majority of the infected region without serious oversegmentation. Then, it is easy to identify and measure the infected portions by using a priori knowledge.Burkholderia cepacia (ex Burkholder) infected onion and the target segmentation.
(a)
Original image
(b)
Segmentation target


### 4.1. Comparison among Different Segmentation Methods

The common image segmentation methods, including watershed algorithm, Sobel operator and Canny edge detection algorithm, Otsu algorithm, and an effective and common quad-tree decomposition algorithm, were selected for comparison. Shannon wavelet was employed to construct the wavelet interpolation operator. The representation of Shannon wavelet is based upon approximating the Dirac delta function as a band-limited function and is given by(32)w(x)=sin(πx)πx.Consider a one-dimensional functionf(x), x∈[a,b]. A discrete point sequence of the variable x is defined as
(33)xn=a+b-a2j·n,j∈ℤ,n=0,1,2,…,2j
and the corresponding discrete point sequence of the scaling function ϕ(x) can be defined as
(34)wj,n(x)=wj(x-xn)=sin(2jπ/(b-a))(x-xn)(2jπ/(b-a))(x-xn).The first and second order derivatives ofϕj(x-xn) at the discrete point xk are
(35)ϕj′(xk-xn)={0,k=n2jcos⁡[π(k-n)](k-n)(b-a),k≠n,ϕj′′(xk-xn)={-π23((b-a)/2j)2,k=n-2cos⁡[π(k-n)]((b-a)/2j)2(k-n)2,k≠n.
The corresponding 2-dimension weight function can be represented as the tensor product form of the above equations.The experimental procedure is described as follows:(1)
convert the infested onion image to grayscale (Figure2(a)) and solve for the gradient map (Figure 2(b));(2)
use the grayscale image to test the Sobel operator, Canny operator, Qtdecomp algorithms, and wavelet precise integration method; Use the gradient mat to test the OSTU method and watershed algorithm;(3)
applying the watershed method to segment the image which has been processed by OSTU method in order to avoid the oversegmentation from the watershed algorithm.Comparison of various segmentation methods.
(a)
Grayscale image
(b)
Gradient image
(c)
OSTU
(d)
Watershed 1
(e)
Watershed 2
(f)
Qtdecomp
(g)
Sobel
(h)
Canny
(i)
Wavelet precise integrationThe segmentation results are shown in Figure2. The watershed algorithm segmentation result shows serious oversegmentation (Figure 2(d)) and cannot recognize the infected part. Although OSTU method separated part infected region of the onion, the partition boundary is discontinuity and is difficult to measure infection specific gravity. To avoid oversegmentation, OSTU was overlapped with the watershed segmentation. The result is shown in Figure 2(e) in which the partition boundary is clear but is unable to distinguish the virus infected part. The Sobel operator recognition on the part of the infection is also not good (Figure 2(g)). Canny operator and Qtdecomp algorithm identified the area of infection, but the boundary points of segmentation region are disorder and cannot be measured. The precise integration method presented in this paper can identify the infected area clearly. So, it is helpful to onion evaluation and classification.In fact, it is impossible to project image segmentation by a single algorithm. The important reason that partial differential equations are effective for image segmentation is that the method integrated many image segmentation principles to the model of partial differential equations. In this paper, the C-V image segmentation model is a global convex optimization variational model, which is established on image piecewise smooth (c1 and c2 are the average gray values inside (Ω1) and outside (Ω2) of the object contour, respectively). To ensure the accuracy of image segmentation, the curvature of the image, the border gradient, and level set function evolution were taken into account in image segmentation. It means that the global convex optimization model of image segmentation has been built based on integration of a variety of image segmentation theories and has obvious advantages. In addition, the method of iterative solution of the self-adaptive method can also be integrated into the segmentation process to ensure the accuracy of segmentation method further. However, the speed of the algorithm will be affected. Therefore, it is important to find efficient and accurate numerical method.


### 4.2. Efficiency Comparison of Multiscale Adaptive Wavelet Numerical Method and Difference Method

The C-V model was used for 256 × 302 images segmentation and divided difference method was used to disperse partial differential equations. So discrete 7312 (256×302) ordinary differential equations are huge solving workload. But the adaptive wavelet precise integration method can reduce the scale to 9576 equations. It can improve solution efficiency greatly due to less workload and low memory demand. Of course, the use of adaptive wavelet precise integration method for solving the number of distribution points will dynamically change as the solution process. In addition, as shown in Figure 3 distribution points are relatively dense within the ellipse ring and another location was sparse. The evident grayscale difference between the infected and the healthy parts lead to this special points distribution. Furthermore, distribution points also exhibit regular matrix form, which result from block solving method of wavelet transform to improve the efficiency. The matrix-like distribution is from the boundary effect among the different blocks. The interval wavelet can effectively reduce the range effect, but it will also increase the computation work of the wavelet transform.Figure 3
Adaptive wavelet collocation points on level set.In this paper, difference method was tested in MATLAB. The wavelet interpolation operator was implemented with VC programming and other parts with MATLAB programming. On the same computer, difference method takes 0.3 seconds, the adaptive wavelet precise integration method takes 0.18 seconds. The results also show that the wavelet transform of the iterative process reduces the overall computational efficiency of the algorithm.


## 4.1. Comparison among Different Segmentation Methods

The common image segmentation methods, including watershed algorithm, Sobel operator and Canny edge detection algorithm, Otsu algorithm, and an effective and common quad-tree decomposition algorithm, were selected for comparison. Shannon wavelet was employed to construct the wavelet interpolation operator. The representation of Shannon wavelet is based upon approximating the Dirac delta function as a band-limited function and is given by(32)w(x)=sin(πx)πx.Consider a one-dimensional functionf(x), x∈[a,b]. A discrete point sequence of the variable x is defined as
(33)xn=a+b-a2j·n,j∈ℤ,n=0,1,2,…,2j
and the corresponding discrete point sequence of the scaling function ϕ(x) can be defined as
(34)wj,n(x)=wj(x-xn)=sin(2jπ/(b-a))(x-xn)(2jπ/(b-a))(x-xn).The first and second order derivatives ofϕj(x-xn) at the discrete point xk are
(35)ϕj′(xk-xn)={0,k=n2jcos⁡[π(k-n)](k-n)(b-a),k≠n,ϕj′′(xk-xn)={-π23((b-a)/2j)2,k=n-2cos⁡[π(k-n)]((b-a)/2j)2(k-n)2,k≠n.
The corresponding 2-dimension weight function can be represented as the tensor product form of the above equations.The experimental procedure is described as follows:(1)
convert the infested onion image to grayscale (Figure2(a)) and solve for the gradient map (Figure 2(b));(2)
use the grayscale image to test the Sobel operator, Canny operator, Qtdecomp algorithms, and wavelet precise integration method; Use the gradient mat to test the OSTU method and watershed algorithm;(3)
applying the watershed method to segment the image which has been processed by OSTU method in order to avoid the oversegmentation from the watershed algorithm.Comparison of various segmentation methods.
(a)
Grayscale image
(b)
Gradient image
(c)
OSTU
(d)
Watershed 1
(e)
Watershed 2
(f)
Qtdecomp
(g)
Sobel
(h)
Canny
(i)
Wavelet precise integrationThe segmentation results are shown in Figure2. The watershed algorithm segmentation result shows serious oversegmentation (Figure 2(d)) and cannot recognize the infected part. Although OSTU method separated part infected region of the onion, the partition boundary is discontinuity and is difficult to measure infection specific gravity. To avoid oversegmentation, OSTU was overlapped with the watershed segmentation. The result is shown in Figure 2(e) in which the partition boundary is clear but is unable to distinguish the virus infected part. The Sobel operator recognition on the part of the infection is also not good (Figure 2(g)). Canny operator and Qtdecomp algorithm identified the area of infection, but the boundary points of segmentation region are disorder and cannot be measured. The precise integration method presented in this paper can identify the infected area clearly. So, it is helpful to onion evaluation and classification.In fact, it is impossible to project image segmentation by a single algorithm. The important reason that partial differential equations are effective for image segmentation is that the method integrated many image segmentation principles to the model of partial differential equations. In this paper, the C-V image segmentation model is a global convex optimization variational model, which is established on image piecewise smooth (c1 and c2 are the average gray values inside (Ω1) and outside (Ω2) of the object contour, respectively). To ensure the accuracy of image segmentation, the curvature of the image, the border gradient, and level set function evolution were taken into account in image segmentation. It means that the global convex optimization model of image segmentation has been built based on integration of a variety of image segmentation theories and has obvious advantages. In addition, the method of iterative solution of the self-adaptive method can also be integrated into the segmentation process to ensure the accuracy of segmentation method further. However, the speed of the algorithm will be affected. Therefore, it is important to find efficient and accurate numerical method.


## 4.2. Efficiency Comparison of Multiscale Adaptive Wavelet Numerical Method and Difference Method

The C-V model was used for 256 × 302 images segmentation and divided difference method was used to disperse partial differential equations. So discrete 7312 (256×302) ordinary differential equations are huge solving workload. But the adaptive wavelet precise integration method can reduce the scale to 9576 equations. It can improve solution efficiency greatly due to less workload and low memory demand. Of course, the use of adaptive wavelet precise integration method for solving the number of distribution points will dynamically change as the solution process. In addition, as shown in Figure 3 distribution points are relatively dense within the ellipse ring and another location was sparse. The evident grayscale difference between the infected and the healthy parts lead to this special points distribution. Furthermore, distribution points also exhibit regular matrix form, which result from block solving method of wavelet transform to improve the efficiency. The matrix-like distribution is from the boundary effect among the different blocks. The interval wavelet can effectively reduce the range effect, but it will also increase the computation work of the wavelet transform.Figure 3
Adaptive wavelet collocation points on level set.In this paper, difference method was tested in MATLAB. The wavelet interpolation operator was implemented with VC programming and other parts with MATLAB programming. On the same computer, difference method takes 0.3 seconds, the adaptive wavelet precise integration method takes 0.18 seconds. The results also show that the wavelet transform of the iterative process reduces the overall computational efficiency of the algorithm.


## 5. Conclusions

Shannon wavelet precise integration method is a new image segmentation method based on the C-V model, which was used to construct adaptive wavelet interpolation operator due to multiscale characteristics of wavelet transform combined with the time precise integration technology. The method makes full use of the multiscale characteristics and the high precise performance of precise integration method. Compared to the gradient method and wavelet transform method of image segmentation, object boundary obtained by WPIM segmentation method is clear and closed; compared to the watershed method, the WPIM method avoids over- and undersegmentation problems and is very suitable for measurement of image segmentation, such as onion quality assessment.The adaptive interpolation operator in the Shannon wavelet precision integration method can reduce the amount of the collocation points and improve the calculation efficiency. As the interpolation operator contains a wavelet transform process, the corresponding algorithm needs to do a wavelet transform between each two iteration time steps. So the cost of the wavelet transformation is an important part of the calculation amount of the algorithm. Compactly supported orthogonal wavelet function can be expected to solve the problem.



---

*Source: 601841-2014-03-19.xml*