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

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.

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 𝐶²-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 |𝑥|→∞. 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: 𝐸^CV𝑐₁,𝑐₂,𝐶=𝜆_1Ω1𝐼₀−𝑐₁²𝑑𝑥𝑑𝑦+𝜆_2Ω2𝐼₀−𝑐₂²|𝐶|,𝑐𝑑𝑥𝑑𝑦+𝜈_𝑖=mean_Ω𝑖𝑢₀=∫_Ω𝑖𝑢₀𝑥,𝑦𝑑𝑥𝑑𝑦AreaΩ_𝑖,𝑖=1,2, (1)(2) where 𝜆₁ and 𝜆₂ are positive constants and 𝑐₁ and 𝑐₂ are the average gray level values inside (Ω₁) and outside (Ω₂) of the object contour, respectively. 𝐼₀ denotes the image to process, |𝐶| 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: ∣,Ω𝐶=𝑥,𝑦𝑥,𝑦∈Ω,𝜙𝑥,𝑦=0₁=∣,Ω𝑥,𝑦𝑥,𝑦∈Ω,𝜙𝑥,𝑦>0₂=𝑥,𝑦∣𝑥,𝑦∈Ω,𝜙𝑥,𝑦<0. (3) Then, the level set-based C-V model can be rewritten as follows: 𝐸𝑐₁,𝑐₂,𝜙=𝜆_1Ω𝐼₀−𝑐₁²𝐻𝜙𝑑𝑥𝑑𝑦+𝜆_2Ω𝐼₀−𝑐₂²1−𝐻𝜙𝑑𝑥𝑑𝑦+𝜈_Ω𝐻𝜙𝐻𝜙=𝛿𝑑𝑥𝑑𝑦,1,𝜙≥00,𝜙<0,_𝜀=𝜀𝜋𝜀²+𝜙². (4) Using the variational method, the PDEs with respect to the variable 𝜙 can be obtained as follows: 𝜕𝜙𝜕𝑡=𝛿_𝜀𝜙𝜈div∇𝜙∇𝜙−𝜆₁𝐼₀−𝑐₁²+𝜆₂𝐼₀−𝑐₂². (5) 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 𝜙(𝑥,𝑦,𝑡) at time 𝑡. The zero level set is the object contour curve, which can be obtained by solving the equation 𝜙(𝑥,𝑦,𝑡)=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 (𝑥_min,𝑥_max)×(𝑦_min,𝑦_max); the discretization points can be defined as (𝑥^𝑗_𝑘1,𝑦^𝑗_𝑘2), where 𝑗 is a scale parameter and 𝑘₁ and 𝑘₂ are position parameters. So 𝑥^𝑗_𝑘1=𝑥_min+𝑘₁𝑥_max−𝑥_min2^𝑗,𝑦^𝑗_𝑘1=𝑦_min+𝑘₂𝑦_max−𝑦_min2^𝑗,𝑗,𝑘₁,𝑘₂∈ℤ. (6)

In addition, 𝑤_𝑘^{𝑗(𝑚,𝑛)}_1,𝑘2(𝑥,𝑦) denotes the multiscale wavelet function and the corresponding 𝑚th and 𝑛th derivatives with respect to 𝑥 and 𝑦, respectively. The level set function 𝜙(𝑥,𝑦,𝑡) and the corresponding derivative function can be discretized as follows: 𝜙^{𝐽(𝑚,𝑛)}=𝑥,𝑦,𝑡_{1𝑘011=0𝑘02=0}𝜙𝑥⁰_𝑘01,𝑦⁰_𝑘02𝑤_𝑘^0(𝑚,𝑛)_01,𝑘02+𝑥,𝑦_{𝐽−12𝑗=0𝑗−1𝑘112=0𝑗−1𝑘12=0}𝛼¹_{𝑗,𝑘11,𝑘12}(𝑡)𝑤^{𝑗+1(𝑚,𝑛)}_{2𝑘11+1,2𝑘12}𝑥,𝑦+𝛼²_{𝑗,𝑘11,𝑘12}(𝑡)𝑤^{𝑗+1(𝑚,𝑛)}_{2𝑘11,2𝑘12+1}𝑥,𝑦+𝛼³_{𝑗,𝑘11,𝑘12}(𝑡)𝑤^{𝑗+1(𝑚,𝑛)}_{2𝑘11+1,2𝑘12+1}𝑥,𝑦, (7) where 𝑗 and 𝐽 are constants, which denote the wavelet scale number and the maximum of the scale number, respectively. 𝛼¹_{𝑗,𝑘11,𝑘12}, 𝛼²_{𝑗,𝑘11,𝑘12}, and 𝛼³_{𝑗,𝑘11,𝑘12} are the wavelet coefficients at the points (𝑥^𝑗_𝑘1,𝑦^𝑗_𝑘2). According to the interpolation wavelet transform theory, the wavelet coefficients can be written as 𝛼¹_{𝑗,𝑘1,𝑘2}=𝜙𝑥_{𝑗+1,2𝑘1+1},𝑦_{𝑗+1,2𝑘2}−𝐼_𝑗𝜙𝑥_{𝑗+1,2𝑘1+1},𝑦_{𝑗+1,2𝑘2},𝛼²_{𝑗,𝑘1,𝑘2}=𝜙𝑥_{𝑗+1,2𝑘1},𝑦_{𝑗+1,2𝑘2+1}−𝐼_𝑗𝜙𝑥_{𝑗+1,2𝑘1},𝑦_{𝑗+1,2𝑘2+1},𝛼³_{𝑗,𝑘1,𝑘2}=𝜙𝑥_{𝑗+1,2𝑘1+1},𝑦_{𝑗+1,2𝑘2+1}−𝐼_𝑗𝜙𝑥_{𝑗+1,2𝑘1+1},𝑦_{𝑗+1,2𝑘2+1}, (8) where 𝐼_𝑗 denotes the multilevel interpolation operator. In order to obtain the multilevel interpolation operator, it is necessary to express the wavelet coefficients 𝛼¹_{𝑗,𝑘1,𝑘2},𝛼²_{𝑗,𝑘1,𝑘2},𝛼³_{𝑗,𝑘1,𝑘2} as a weighted sum of 𝑢 in all of the collocation points in the 𝐽-level. Therefore, we should give the definition of the restriction operator as follows: 𝑅_𝑘^{𝑙,𝑙,𝑗,𝑗}_{1,𝑘2,𝑚1,𝑚2}=1,𝑥^𝑙_𝑘1=𝑥^𝑗_𝑚1,𝑦^𝑙_𝑘2=𝑦^𝑗_𝑚20,otherwise. (9) Using the restriction operator, 𝑢(𝑥^𝑗+1_2𝑘1+1,𝑦^𝑗+1_2𝑘2),𝑢(𝑥^𝑗+1_2𝑘1,𝑦^𝑗+1_2𝑘2+1) and 𝑢(𝑥^𝑗+1_2𝑘1+1,𝑦^𝑗+1_2𝑘2+1) can be rewritten as 𝜙𝑥^𝑗+1_2𝑘1+1,𝑦^𝑗+1_2𝑘2=^2𝐽_𝑛1²₌₀^𝐽_𝑛2=0𝑅^{𝑗+1,𝑗+1,𝐽,𝐽}_{2𝑘1+1,2𝑘2,𝑛1,𝑛2}𝜙𝑥^𝐽_𝑛1,𝑦^𝐽_𝑛2,𝜙𝑥^𝑗+1_2𝑘1,𝑦^𝑗+1_2𝑘2+1=^2𝐽_𝑛1²₌₀^𝐽_𝑛2=0𝑅^{𝑗+1,𝑗+1,𝐽,𝐽}_{2𝑘1,2𝑘2+1,𝑛1,𝑛2}𝜙𝑥^𝐽_𝑛1,𝑦^𝐽_𝑛2,𝜙𝑥^𝑗+1_2𝑘1+1,𝑦^𝑗+1_2𝑘2+1=^2𝐽_𝑛1²₌₀^𝐽_𝑛2=0𝑅^{𝑗+1,𝑗+1,𝐽,𝐽}_{2𝑘1+1,2𝑘2+1,𝑛1,𝑛2}𝜙𝑥^𝐽_𝑛1,𝑦^𝐽_𝑛2. (10) Introducing the extension operators 𝐶1, 𝐶2, and 𝐶3 and substituting (10) into (8), the wavelet coefficients can be rewritten as𝛼¹_{𝑗,𝑘1,𝑘2}=^2𝐽_𝑛1²₌₀^𝐽_𝑛2=0𝑅^{𝑗+1,𝑗+1,𝐽,𝐽}_{2𝑘1+1,2𝑘2,𝑛1,𝑛2}𝜙𝑥^𝐽_𝑛1,𝑦^𝐽_𝑛2−^2𝐽_𝑛1²₌₀^𝐽_𝑛2²₌₀^𝑗0_𝑘01²₌₀^𝑗0_𝑘02=0𝑅^𝑗0,𝑗0_𝑘^,𝐽,𝐽_{01,𝑘02,𝑛1,𝑛2}𝜙𝑥^𝐽_𝑛1,𝑦^𝐽_𝑛2𝑤^𝑗0_{𝑘01,𝑘02}𝑥^𝑗+1_2𝑘1+1,𝑦^𝑗+1_2𝑘2+^𝑗−1_𝑗1=𝑗0^2𝐽_𝑛1²₌₀^𝐽_𝑛2²₌₀^𝑗1_𝑘11²₌₀^𝑗1_𝑘12=0𝐶1^𝑗1,𝑗1_𝑘^,𝐽,𝐽_{11,𝑘12,𝑛1,𝑛2}𝑤^𝑗1+1_{2𝑘11+1,2𝑘12}𝑥^𝑗+1_2𝑘1+1,𝑦^𝑗+1_2𝑘2𝜙𝑥^𝐽_𝑛1,𝑦^𝐽_𝑛2+𝐶2^𝑗1,𝑗1_𝑘^,𝐽,𝐽_{11,𝑘12,𝑛1,𝑛2}𝑤^𝑗1+1_{2𝑘11,2𝑘12+1}𝑥^𝑗+1_2𝑘1+1,𝑦^𝑗+1_2𝑘2𝑢𝑥^𝐽_𝑛1,𝑦^𝐽_𝑛2+𝐶3^𝑗1,𝑗1_𝑘^,𝐽,𝐽_{11,𝑘12,𝑛1,𝑛2}𝑤^𝑗1+1_{2𝑘11+1,2𝑘12+1}𝑥^𝑗+1_2𝑘1+1,𝑦^𝑗+1_2𝑘2𝜙𝑥^𝐽_𝑛1,𝑦^𝐽_𝑛2=^2𝐽_𝑛1²₌₀^𝐽_𝑛2=0𝐶1_𝑘^{𝑗,𝑗,𝐽,𝐽}_{1,𝑘2,𝑛1,𝑛2}𝜙𝑥^𝐽_𝑛1,𝑦^𝐽_𝑛2. (11)𝛼²_{𝑗,𝑘1,𝑘2} and 𝛼³_{𝑗,𝑘1,𝑘2} are similar to 𝛼¹_{𝑗,𝑘1,𝑘2}. From the above equation, the extension operator can be obtained as𝐶1_𝑘^{𝑗,𝑗,𝐽,𝐽}_{1,𝑘2,𝑛1,𝑛2}=𝑅^{𝑗+1,𝑗+1,𝐽,𝐽}_{2𝑘1+1,2𝑘2,𝑛1,𝑛2}−^2𝑗0_𝑘01²₌₀^𝑗0_𝑘02=0𝑅^𝑗0,𝑗0_𝑘^,𝐽,𝐽_{01,𝑘02,𝑛1,𝑛2}𝜙𝑥^𝐽_𝑛1,𝑦^𝐽_𝑛2𝑤^𝑗0_{𝑘01,𝑘02}𝑥^𝑗+1_2𝑘1+1,𝑦^𝑗+1_2𝑘2+^𝑗−1_𝑗1=𝑗0^2𝐽_𝑛2²₌₀^𝑗1_𝑘11=0𝐶1^𝑗1,𝑗1_𝑘^,𝐽,𝐽_{11,𝑘12,𝑛1,𝑛2}𝑤^𝑗1+1_{2𝑘11+1,2𝑘12}𝑥^𝑗+1_2𝑘1+1,𝑦^𝑗+1_2𝑘2𝜙𝑥^𝐽_𝑛1,𝑦^𝐽_𝑛2+𝐶2^𝑗1,𝑗1_𝑘^,𝐽,𝐽_{11,𝑘12,𝑛1,𝑛2}𝑤^𝑗1+1_{2𝑘11,2𝑘12+1}𝑥^𝑗+1_2𝑘1+1,𝑦^𝑗+1_2𝑘2𝜙𝑥^𝐽_𝑛1,𝑦^𝐽_𝑛2+𝐶3^𝑗1,𝑗1_𝑘^,𝐽,𝐽_{11,𝑘12,𝑛1,𝑛2}𝑤^𝑗1+1_{2𝑘11+1,2𝑘12+1}𝑥^𝑗+1_2𝑘1+1,𝑦^𝑗+1_2𝑘2𝜙𝑥^𝐽_𝑛1,𝑦^𝐽_𝑛2. (12)𝐶2 and 𝐶3 can be obtained with the same method. Therefore, the calculation time complexity of the wavelet transform coefficients 𝛼¹_{𝑗,𝑘11,𝑘12}, 𝛼²_{𝑗,𝑘11,𝑘12}, and 𝛼³_{𝑗,𝑘11,𝑘12} is 𝑂((1/3)4^2𝐽−1).

Substituting 𝛼¹_{𝑗,𝑘11,𝑘12}, 𝛼²_{𝑗,𝑘11,𝑘12}, and 𝛼³_{𝑗,𝑘11,𝑘12}, 𝐶1, 𝐶2, and 𝐶3 into (2), the multilevel wavelet interpolation operator can be obtained as 𝐼_𝑛1,𝑛2=𝑥,𝑦_{2𝑗0𝑘012=0𝑗0𝑘02=0}𝑅^𝑗0,𝑗0_𝑘^,𝐽,𝐽_{01,𝑘02,𝑛1,𝑛2}𝑤^𝑗0_{𝑘01,𝑘02}+𝑥,𝑦_{𝐽−1𝑗=𝑗02𝑗𝑘1=02𝑗𝑘2=0}𝐶1_𝑘^{𝑗,𝑗,𝐽,𝐽}_{1,𝑘2,𝑛1,𝑛2}𝑤^𝑗+1_{2𝑘1+1,2𝑘2}𝑥,𝑦𝜙𝑥^𝐽_𝑛1,𝑦^𝐽_𝑛2+𝐶2_𝑘^{𝑗,𝑗,𝐽,𝐽}_{1,𝑘2,𝑛1,𝑛2}𝑤^𝑗+1_{2𝑘1,2𝑘2+1}𝑥,𝑦𝜙𝑥^𝐽_𝑛1,𝑦^𝐽_𝑛2+𝐶3_𝑘^{𝑗,𝑗,𝐽,𝐽}_{1,𝑘2,𝑛1,𝑛2}𝑤^𝑗+1_{2𝑘1+1,2𝑘2+1}𝑥,𝑦𝜙𝑥^𝐽_𝑛1,𝑦^𝐽_𝑛2. (13) Then, (7) can be rewritten as 𝜙^{𝐽(𝑚,𝑛)}𝑥,𝑦,𝑡=^2𝐽_𝑛1^2𝐽_𝑛2𝐼_𝑛1,𝑛2𝑥,𝑦𝜙𝑥^𝐽_𝑛1,𝑦^𝐽_𝑛2. (14) 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 parameter 𝑡. In other words, the wavelet coefficient should vary with the time parameter 𝑡. 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 𝜕𝜙𝜕𝑡=𝐹𝑡,𝑥,𝑦,𝜙,𝜕𝜙,𝜕𝑥𝜕𝜙,𝜕𝜕𝑦²𝜙𝜕𝑥²,𝜕²𝜙,𝜕𝜕𝑥𝜕𝑦²𝜙𝜕𝑦²()𝜙𝑡>0𝑥,𝑦,0=𝜙₀𝑥,𝑦,𝑑𝜙^𝐽𝑥,𝑦,𝑡𝑑𝑡=𝐹𝑡,𝑥,𝑦,𝜙^𝐽𝑥,𝑦,𝑡,𝜙^𝐽(1,0),𝜙𝑥,𝑦,𝑡_𝐽(0,1)𝑥,𝑦,𝑡,𝜙^𝐽(2,0),𝜙𝑥,𝑦,𝑡_𝐽(1,1)𝑥,𝑦,𝑡,𝜙^𝐽(0,2)𝑥,𝑦,𝑡, (15) respectively, where 𝜙^𝐽=𝑥,𝑦,𝑡_{1𝑘011=0𝑘02=0}𝜙𝑥⁰_𝑘01,𝑦⁰_𝑘02𝑤⁰_{𝑘01,𝑘02}+𝑥,𝑦_{𝐽−12𝑗=0𝑗−1𝑘112=0𝑗−1𝑘12=0}𝛼¹_{𝑗,𝑘11,𝑘12}𝑤^𝑗+1_{2𝑘11+1,2𝑘12}𝑥,𝑦+𝛼²_{𝑗,𝑘11,𝑘12}𝑤^𝑗+1_{2𝑘11,2𝑘12+1}𝑥,𝑦+𝛼³_{𝑗,𝑘11,𝑘12}𝑤^𝑗+1_{2𝑘11+1,2𝑘12+1}𝑥,𝑦. (16) The value of 𝜙^𝐽(𝑥,𝑦,𝑡_𝑛) at 𝑡_𝑛 is denoted as 𝜙_𝑛, and 𝐹𝑡_𝑛,𝑥,𝑦,𝜙^𝐽𝑥,𝑦,𝑡_𝑛,𝜙^𝐽(1,0)𝑥,𝑦,𝑡_𝑛,𝜙^𝐽(0,1)𝑥,𝑦,𝑡_𝑛,𝜙^𝐽(2,0)𝑥,𝑦,𝑡_𝑛,𝜙^𝐽(1,1)𝑥,𝑦,𝑡_𝑛,𝜙^𝐽(0,2)𝑥,𝑦,𝑡_𝑛 (17) is denoted as 𝐹_𝑛. And then a linear homotopy function can be constructed as 𝜙^𝐽𝑥,𝑦,𝑡=(1−𝜀)𝐹_𝑛+𝜀𝐹_𝑛+1. (18) It is easy to identify the homotopy parameter as 𝜀(𝑡)=𝑡−𝑡_𝑛𝑡_𝑛+1−𝑡_𝑛𝑡∈𝑡_𝑛,𝑡_𝑛+1∴𝜀∈[0,1]. (19) According to the perturbation theory, the solution of (18) can be expressed as the power series expansion of 𝜀 as follows: 𝜙^𝐽=𝜙^𝐽₀+𝜀𝜙^𝐽₁+𝜀²𝜙^𝐽₂+⋯. (20) Substituting (20) into (18) and rearranging based on powers of 𝜀-terms, we have 𝜀⁰:𝜙^𝐽₀=𝐹_𝑛,𝜀¹:𝜙^𝐽₁=𝐹_𝑛+1−𝐹_𝑛,⋮ (21) According to HPM, we obtain the wavelet coefficients 𝛼¹_{𝑗,𝑘1,𝑘2}(𝑡_𝑛+1), 𝛼²_{𝑗,𝑘1,𝑘2}(𝑡_𝑛+1), 𝛼³_{𝑗,𝑘1,𝑘2}(𝑡_𝑛+1) at 𝑡_𝑛 as follows:𝛼¹_{𝑗,𝑘1,𝑘2}=𝜙𝑥^𝑗+1_2𝑘1+1,𝑦^𝑗+1_2𝑘2−𝐼_𝑗𝜙𝑥^𝑗+1_2𝑘1+1,𝑦^𝑗+1_2𝑘2=𝜙𝑥^𝑗+1_2𝑘1+1,𝑦^𝑗+1_2𝑘2−¹_𝑘01¹_=0𝑘02=0𝜙𝑥⁰_𝑘01,𝑦⁰_𝑘02𝑤⁰_{𝑘01,𝑘02}𝑥^𝑗+1_2𝑘1+1,𝑦^𝑗+1_2𝑘2+^𝑗−1_𝑗1²₌₀^𝑗1_𝑘11²₌₀^𝑗1_𝑘12=0𝛼¹_{𝑗1,𝑘11,𝑘12}𝑤^𝑗1+1_{2𝑘11+1,2𝑘12}𝑥^𝑗+1_2𝑘1+1,𝑦^𝑗+1_2𝑘2+𝛼²_{𝑗1,𝑘11,𝑘12}𝑤^𝑗1+1_{2𝑘11,2𝑘12+1}𝑥^𝑗+1_2𝑘1+1,𝑦^𝑗+1_2𝑘2+𝛼³_{𝑗1,𝑘11,𝑘12}𝑤^𝑗1+1_{2𝑘11+1,2𝑘12+1}𝑥^𝑗+1_2𝑘1+1,𝑦^𝑗+1_2𝑘2,𝛼²_{𝑗,𝑘1,𝑘2}=𝜙𝑥^𝑗+1_2𝑘1,𝑦^𝑗+1_2𝑘2+1−𝐼_𝑗𝜙𝑥^𝑗+1_2𝑘1,𝑦^𝑗+1_2𝑘2+1=𝜙𝑥^𝑗+1_2𝑘1,𝑦^𝑗+1_2𝑘2+1−¹_𝑘01¹_=0𝑘02=0𝜙𝑥⁰_𝑘01,𝑦⁰_𝑘02𝑤⁰_{𝑘01,𝑘02}𝑥^𝑗+1_2𝑘1,𝑦^𝑗+1_2𝑘2+1+^𝑗−1_𝑗1²₌₀^𝑗1_𝑘11²₌₀^𝑗1_𝑘12=0𝛼¹_{𝑗1,𝑘11,𝑘12}𝑤^𝑗1+1_{2𝑘11+1,2𝑘12}𝑥^𝑗+1_2𝑘1,𝑦^𝑗+1_2𝑘2+1+𝛼²_{𝑗1,𝑘11,𝑘12}𝑤^𝑗1+1_{2𝑘11,2𝑘12+1}𝑥^𝑗+1_2𝑘1,𝑦^𝑗+1_2𝑘2+1+𝛼³_{𝑗1,𝑘11,𝑘12}𝑤^𝑗1+1_{2𝑘11+1,2𝑘12+1}𝑥^𝑗+1_2𝑘1,𝑦^𝑗+1_2𝑘2+1,𝛼³_{𝑗,𝑘1,𝑘2}=𝜙𝑥^𝑗+1_2𝑘1+1,𝑦^𝑗+1_2𝑘2+1−𝐼_𝑗𝜙𝑥^𝑗+1_2𝑘1+1,𝑦^𝑗+1_2𝑘2+1=𝜙𝑥^𝑗+1_2𝑘1+1,𝑦^𝑗+1_2𝑘2+1−¹_𝑘01¹_=0𝑘02=0𝜙𝑥⁰_𝑘01,𝑦⁰_𝑘02𝑤⁰_{𝑘01,𝑘02}𝑥^𝑗+1_2𝑘1+1,𝑦^𝑗+1_2𝑘2+1+^𝑗−1_𝑗1²₌₀^𝑗1_𝑘11²₌₀^𝑗1_𝑘12=0𝛼¹_{𝑗1,𝑘11,𝑘12}𝑤^𝑗1+1_{2𝑘11+1,2𝑘12}𝑥^𝑗+1_2𝑘1+1,𝑦^𝑗+1_2𝑘2+1+𝛼²_{𝑗1,𝑘11,𝑘12}𝑤^𝑗1+1_{2𝑘11,2𝑘12+1}𝑥^𝑗+1_2𝑘1+1,𝑦^𝑗+1_2𝑘2+1+𝛼³_{𝑗1,𝑘11,𝑘12}𝑤^𝑗1+1_{2𝑘11+1,2𝑘12+1}𝑥^𝑗+1_2𝑘1+1,𝑦^𝑗+1_2𝑘2+1. (22)Obviously, the calculation time complexity of the wavelet transform coefficients 𝛼¹_{𝑗,𝑘1,𝑘2}, 𝛼²_{𝑗,𝑘1,𝑘2}, and 𝛼³_{𝑗,𝑘1,𝑘2} is 𝑂(4^𝐽), which is decreased greatly than that in (8) which is 𝑂((1/3)4^2𝐽−1).

Substituting the wavelet transform efficient (22) into (16), we obtain 𝜙^𝐽𝑥,𝑦,𝑡_𝑛+1=𝜙^𝐽𝑥,𝑦,𝑡_𝑛+Δ𝑡2𝐹𝑡_𝑛,𝑥,𝑦,𝜙^𝐽𝑥,𝑦,𝑡_𝑛,𝜙^𝐽(1,0)𝑥,𝑦,𝑡_𝑛,𝜙^𝐽(0,1)𝑥,𝑦,𝑡_𝑛,𝜙^𝐽(2,0)𝑥,𝑦,𝑡_𝑛,𝜙^𝐽(1,1)𝑥,𝑦,𝑡_𝑛,𝜙^𝐽(0,2)𝑥,𝑦,𝑡_𝑛+𝐹𝑡_𝑛+1,𝑥,𝑦,𝜙^𝐽₀𝑥,𝑦,𝑡_𝑛+1,𝜙₀^𝐽(1,0)𝑥,𝑦,𝑡_𝑛+1,𝜙₀^𝐽(0,1)𝑥,𝑦,𝑡_𝑛+1,𝜙₀^𝐽(2,0)𝑥,𝑦,𝑡_𝑛+1,𝜙₀^𝐽(1,1)𝑥,𝑦,𝑡_𝑛+1,𝜙₀^𝐽(0,2)𝑥,𝑦,𝑡_𝑛+1 (23) and the derivative function 𝜙^{𝐽(𝑚,𝑛)}=𝑥,𝑦_{1𝑘011=0𝑘02=0}𝜙𝑥⁰_𝑘01,𝑦⁰_𝑘02𝑤_𝑘^0(𝑚,𝑛)_01,𝑘02+𝑥,𝑦_{𝐽−12𝑗=0𝑗−1𝑘112=0𝑗−1𝑘12=0}𝛼¹_{𝑗,𝑘11,𝑘12}𝑤^{𝑗+1(𝑚,𝑛)}_{2𝑘11+1,2𝑘12}𝑥,𝑦+𝛼²_{𝑗,𝑘11,𝑘12}𝑤^{𝑗+1(𝑚,𝑛)}_{2𝑘11,2𝑘12+1}𝑥,𝑦+𝛼³_{𝑗,𝑘11,𝑘12}𝑤^{𝑗+1(𝑚,𝑛)}_{2𝑘11+1,2𝑘12+1}𝑥,𝑦. (24) 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 to 𝑥 and 𝑦, respectively, can be expressed as follows: 𝜙^{𝐽(𝑚,𝑛)}=𝑥,𝑦,𝑡_{1𝑘011=0𝑘02=0}𝜙𝑥⁰_𝑘01,𝑦⁰_𝑘02𝑤_𝑘^0(𝑚,𝑛)_01,𝑘02+𝑥,𝑦_{𝐽−12𝑗=0𝑗−1𝑘112=0𝑗−1𝑘12=0}𝛼¹_{𝑗,𝑘11,𝑘12}(𝑡)𝑤^{𝑗+1(𝑚,𝑛)}_{2𝑘11+1,2𝑘12}𝑥,𝑦+𝛼²_{𝑗,𝑘11,𝑘12}(𝑡)𝑤^{𝑗+1(𝑚,𝑛)}_{2𝑘11,2𝑘12+1}𝑥,𝑦+𝛼³_{𝑗,𝑘11,𝑘12}(𝑡)𝑤^{𝑗+1(𝑚,𝑛)}_{2𝑘11+1,2𝑘12+1}𝑥,𝑦, (25) where 𝑤_𝑘^{𝑗(𝑚,𝑛)}_1,𝑘2(𝑥,𝑦) is the wavelet function and its 𝑚- and 𝑛-order derivative with respect to 𝑥 and 𝑦, respectively. 𝛼¹_{𝑗,𝑘11,𝑘12}(𝑡), 𝛼²_{𝑗,𝑘11,𝑘12}(𝑡), and 𝛼³_{𝑗,𝑘11,𝑘12}(𝑡) are wavelet transform coefficients. We convert (2) to wavelet multiscale discrete format of the level set function by stead of (25): 𝑑𝜙^𝐽𝑥,𝑦,𝑡𝑑𝑡=𝛿_𝜀𝜙^𝐽×𝑥,𝑦,𝑡𝜈div∇𝜙^𝐽𝑥,𝑦,𝑡∇𝜙^𝐽𝑥,𝑦,𝑡−𝜆₁𝐼₀−𝑐₁²+𝜆₂𝐼₀−𝑐₂². (26)

It is the most direct way for dynamic adaptation to only retain the distribution points of corresponding wavelet coefficients that satisfy the condition min𝛼¹_{𝑗,𝑘11,𝑘12}(𝑡),𝛼²_{𝑗,𝑘11,𝑘12}(𝑡),𝛼³_{𝑗,𝑘11,𝑘12}(𝑡)≥𝜀. (27)

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) where 𝑠, 𝑗 are numbers of different scale wavelet and 𝑘,𝑖,𝑀∈𝑍, 𝜀_𝑥,𝜀_𝑦 are constant.

3.2. Nonlinear Discrete Ordinary Differential Equations

Because 𝜈 is a small parameter in ordinary differential equation (26), value of 𝜈div(∇𝜙^𝐽(𝑥,𝑦,𝑡)/|∇𝜙^𝐽(𝑥,𝑦,𝑡)|) is low. Equation (26) can be converted into 𝑑𝜙^𝐽𝑥,𝑦,𝑡𝑑𝑡=𝛿_𝜀𝜙^𝐽×𝑥,𝑦,𝑡𝜈div∇𝜙^𝐽𝑥,𝑦,𝑡₀∇𝜙^𝐽𝑥,𝑦,𝑡₀−𝜆₁𝐼₀−𝑐₁²+𝜆₂𝐼₀−𝑐₂². (29) The solution of (29) is 𝜙^𝐽𝑥,𝑦,𝑡=𝜀𝑚2𝜋^1/3+𝑚2𝜋𝜀^5−1/3, (30) where 𝑚=𝜈div∇𝜙^𝐽𝑥,𝑦,𝑡₀∇𝜙^𝐽𝑥,𝑦,𝑡₀−𝜆₁𝐼₀−𝑐₁²+𝜆₂𝐼₀−𝑐₂². (31)

We can get the solution of ordinary differential equations (26) by iterative solution (30).

4. Experiment and Discussion of Onion Infected Region Segmentation

Figure 1(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.

(a) Original image

(b) Segmentation target

Figure 1 

Burkholderia cepacia (ex Burkholder) infected onion and the target segmentation.

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 𝑤(𝑥)=sin(𝜋𝑥).𝜋𝑥 (32)

Consider a one-dimensional function 𝑓(𝑥), 𝑥∈[𝑎,𝑏]. A discrete point sequence of the variable 𝑥 is defined as 𝑥_𝑛=𝑎+𝑏−𝑎2^𝑗⋅𝑛,𝑗∈ℤ,𝑛=0,1,2,…,2^𝑗 (33) and the corresponding discrete point sequence of the scaling function 𝜙(𝑥) can be defined as 𝑤_𝑗,𝑛(𝑥)=𝑤_𝑗𝑥−𝑥_𝑛=sin2^𝑗𝜋/(𝑏−𝑎)𝑥−𝑥_𝑛2^𝑗𝜋/(𝑏−𝑎)𝑥−𝑥_𝑛. (34)

The first and second order derivatives of 𝜙_𝑗(𝑥−𝑥_𝑛) at the discrete point 𝑥_𝑘 are 𝜙^′_𝑗𝑥_𝑘−𝑥_𝑛=20,𝑘=𝑛_𝑗cos[𝜋(𝑘−𝑛)]()𝜙𝑘−𝑛)(𝑏−𝑎,𝑘≠𝑛,_′′𝑗𝑥_𝑘−𝑥_𝑛=−𝜋²3(𝑏−𝑎)/2^𝑗2−[𝜋(,𝑘=𝑛2cos𝑘−𝑛)](𝑏−𝑎)/2^𝑗2(𝑘−𝑛)²,𝑘≠𝑛. (35) 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 (Figure 2(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.

(a) Grayscale image

(b) Gradient image

(c) OSTU

(d) Watershed 1

(e) Watershed 2

(f) Qtdecomp

(g) Sobel

(h) Canny

(i) Wavelet precise integration

Figure 2 

Comparison of various segmentation methods.

The segmentation results are shown in Figure 2. 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 (𝑐₁ and 𝑐₂ are the average gray values inside (Ω₁) and outside (Ω₂) 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.

Conflict of Interests

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

Acknowledgments

The authors appreciate the funding support from National Natural Science Foundation of China (Award nos. 41171184 and 41171337). The software tools were provided by the Food and Fiber Sensing Lab of University of Georgia and the Computer Center of China Agricultural University.