Evolution-Operator-Based Single-Step Method for Image Processing

This work proposes an evolution-operator-based single-time-step method for image and signal processing. The key component of the proposed method is a local spectral evolution kernel (LSEK) that analytically integrates a class of evolution partial differential equations (PDEs). From the point of view PDEs, the LSEK provides the analytical solution in a single time step, and is of spectral accuracy, free of instability constraint. From the point of image/signal processing, the LSEK gives rise to a family of lowpass filters. These filters contain controllable time delay and amplitude scaling. The new evolution operator-based method is constructed by pointwise adaptation of anisotropy to the coefficients of the LSEK. The Perona-Malik-type of anisotropic diffusion schemes is incorporated in the LSEK for image denoising. A forward-backward diffusion process is adopted to the LSEK for image deblurring or sharpening. A coupled PDE system is modified for image edge detection. The resulting image edge is utilized for image enhancement. Extensive computer experiments are carried out to demonstrate the performance of the proposed method. The major advantages of the proposed method are its single-step solution and readiness for multidimensional data analysis.


INTRODUCTION
Image denoising, restoration, edge detection, and enhancement are fundamental problems in image processing, computer vision, and artificial intelligence [1]. The solution to these problems is crucial to automatic control, robotics, imaging, target tracking, and telecommunication. A variety of methods have been proposed to tackle these problems in the past few decades. Among these methods, partialdifferential-equation-(PDE-) based methods have recently received much attention. Witkin [2] introduced a parabolic evolution equation (i.e., the heat equation) for image denoising. The basic idea behind Witkin's approach is to evolve an original image under a diffusion process so as to effectively remove noise in the image surface. Such process is formally equivalent to the standard Gaussian filter. His approach has been intensively studied in terms of the scalespace process [3][4][5]. Based on a discrete version of the diffusion equation where Wei and Zhao showed that [6] the effect of the diffusion after S iterations can be expressed as a convolution where the filter W(k, S) can be explicitly given by Therefore, the result of the multiple-step time evolution described in discrete diffusion equation (1) can be exactly A two-dimensional (2D) generalization of (3) is straightforward. Wei and Zhao also showed that W(k, S) is a generalization of the Hanning filter [7] and applied it to chaos synchronization [6] and trend estimation [8].
Like the Gaussian filter, the diffusion process not only removes noise, but also smears image edges and leads to poor visual effects. To address such an issue, Perona and Malik proposed anisotropic diffusion process [9] in which the diffusion coefficient is replaced by a function of position gradient ∂u(r, t) ∂t = ∇ · d |∇u(r, t)| ∇u(r, t) , u(r, 0) = I(r), where I(r) is the original image and d(|∇u|) is a generalized diffusion coefficient which is so designed that its values are very small at the edge of an image. Perona and Malik suggested the Gaussian d(x) = exp[−x 2 /2σ 2 ] and the Lorentz d(x) = [1 + (x/σ) 2 ] −1 , where σ is a constant to be tuned for a particular application. The Perona-Malik equation has recently stimulated much interest in image processing and applied mathematical communities [5,[10][11][12][13][14][15][16][17][18][19]. It is commonly believed that the Perona-Malik equation facilitates a new and potentially more effective algorithm for noise removing, image restoration, edge detection, and image enhancement. Catté et al. [11] pointed out the possibility of generating additional maximum and minimum which do not belong to the initial image data. The problem was analyzed in detail by Kichenassamy [12]. Further studies [5,11,13,19] showed that the anisotropic diffusion process may break down when the gradient generated by noise is comparable to image edges and features. A preconvolution with a smoothing function, such as a Gaussian filter was proposed as a regularization procedure to alleviate the instability. Such a preprocess in fact degrades image quality. One of the present authors [20] proposed a statistical method which discriminates noise from image edges. The idea is to choose σ in the Gaussian or Lorentz as a local statistical variance of |∇u|, σ(r, t) = |∇u − ∇u| 2 , where X(r) denotes the local average of X(r) centered at point r. The area of the local average can be selected in particular application. This simple approach works well in restoration of noisy images. More sophisticated approaches, including variational methods, curvature, and active contours, have been extensively explored in the literature [21][22][23][24][25][26].
Another problem with the original Perona-Malik equation concerns its efficiency in noisy removing and image enhancement. This problem was addressed by Wei [20] by introducing high-order edge-controlled diffusion operators. A fourth-order variation-based PDE was proposed by Chan et al. [27] for image restoration. You and Kaveh [28] introduced a fourth-order curvature diminishing diffusion equation. The high-order Perona-Malik equation [20] takes the form ∂u(r, t) ∂t = q ∇ · d q u(r, t), ∇u(r, t) ∇∇ 2q u(r, t) + e u(r, t), ∇u(r, t) .
A linear and discrete version of the fourth-order diffusion operator was tested for time series analysis [6]. Clearly, the motivations behind the generalized Perona-Malik equation (10) include high-order phenomenological kinetic equations for pattern formation in alloys, glasses, polymer, combustion, and biological systems, and hyperdiffusion used to stabilize the direct numerical simulation of turbulence flow by damping small-scale oscillations around the Nyquist frequency. Image sharpening can be achieved by a balance between forward diffusion and backward diffusion in (10) via appropriately choosing the signs of d 1 and d 2 . Bertozzi and Greer [29][30][31] have carried out mathematical analysis of the generalized Perona-Malik equation (10). The latter was found to be efficient for treating medical magnetic resonance images [32]. Gilboa et al. [33] constructed an efficient image sharpening scheme by using (10) with a triple-well potential for the kinetic production e(u(r, t), |∇u(r, t)|). Witelski and Bowen [34] proposed alternating-direction implicit (ADI) schemes for the integration of (9). The other problem in the Perona-Malik equation is its inefficiency for image edge detection, particularly for edge detection of images with a large amount of texture. This is due to the fact that edge detection is naturally a highpass filtering operation, while the diffusion process is inherently a lowpass filtering process. Wei and Jia [35] address this problem by introducing a pair of weakly coupled nonlinear equations: where ε 1 and ε 2 are the coupling strengths. Here, F 1 and F 2 are nonlinear functions and can both be chosen as the Perona-Malik operator provided that the common initial value is a digital image field u(r, 0) = v(r, 0) = I(r). This coupled PDE approach has been shown [35] to be very robust and efficient, and it provides superior results in image edge detection compared to those obtained by using other existing approaches, such as the Sobel, Prewitt, and Canny operators, and by anisotropic diffusion. Furthermore, the time in the evolution equation is an extra parameter for images and often does not have much justification for its necessity. In particular, there is no rigorous and common rule to determine the period of the time evolution. An automatically controlled generalized Perona-Malik operator was proposed [20] by introducing a timedependent factor where e 0 and t 0 are positive constants. Since each term in the generalized Perona-Malik equation may have a different time dependence, this actually turns the time into an additional dimension for optimizing image properties. Such an idea was adopted and the related technique was further explored by Gilboa et al. [36]. They proposed two time-dependent diffusion coefficients for (7) in the form of with k(t) = 1/( + α 1 t) and where α 1 , α 2 are parameters for controlling the diffusion rate, is a constant, and k is the gradient threshold. Controlled evolution is an important concept for Perona-Malik type of evolution equations. In fact, two types of controls have been introduced in the literature. One control, originally proposed by Wei [20], is to control the magnitude of an individual operator as an explicit function of time, such as (15), so that its effect on the image can be optimized. The other is to control the exit time of the time evolution based on the characteristic of the evolving image. The second type of control was proposed by Wei and Jia [35] to stop the time evolution of their couple edge detection equations based on the difference in variance of two evolving images. Certainly, more research work on the automatic control of time evolution is required in order to make PDE-based image processing a really competitive approach against a host of other methods, such as wavelets, Wiener filter, neural networks, and so forth.
One advantage of PDE-based image processing methods is being elegant and mathematically rigorous. Another advantage is their readiness for systematic generalization to the processing of three-dimensional (3D) and higherdimensional data. Geometric selectivity and parameter controllability also endow these methods great flexibility and potential for a variety of applications in medical imaging, astronomic imaging, pattern recognition, and computer vision. Nevertheless, there are obstacles which prevent PDE-based image processing methods from being used in practical applications, that is, industrial codes. First of all, multiple iterations involved in solving nonlinear PDEs are in general much more expensive than the single-step operation of most classic image processing methods. This lack of efficiency is very severe for the processing of large 3D medical image data. Moreover, the numerical solution of Perona-Malik type of nonlinear PDEs is nontrivial. Numerical instability exists when the time step is too large in explicit schemes. For implicit schemes, stability may not be a problem, but other undesired features might be introduced to the solution.
Many efficient solution techniques were proposed [37][38][39][40] for the linear heat equation. However, when the emphasis of the diffusion process is changed from the linear filtering to the nonlinear one, poor computational efficiency becomes a main obstruction in the practical applications of the nonlinear diffusion equation. In time, explicit schemes like forward Euler (first-order), multilevel explicit finite-difference schemes are most widely used. However, these explicit schemes require very small time steps in order to be stable. It is doomed that the whole diffusion procedure is rather time consuming. Consequently, accelerated explicit schemes are desirable. A fast explicit numerical scheme (λ-resolution) was presented to approximate the solution of the linear diffusion filtering [41]. More sophisticated and absolutely stable approaches are semi-implicit schemes [11,18]. Weickert et al. [18] proposed an additive operator splitting (AOS) scheme for nonlinear diffusion process. In contrast to traditional multiplicative splitting such as the ADI and locally one-dimensional splitting, all axes are treated in the same manner in AOS scheme. However, this scheme needs to compute matrix inversions which are usually a source of computational errors and require much computer memory. In space, since the pixel structure of digital images provides a natural discretization on a fixed rectangular grid, it is not surprising that finite-difference-based methods, such as secondorder, fourth-order, and sixth-order central finite-difference schemes, are commonly used for the task of discretization. Alternatives to the finite-difference methods include finiteelement methods [42][43][44], finite-volume schemes [45,46], pseudospectral approaches [47], wavelets [20,47,48], multigrid methods [49], fast level-set methods [50], high-order ENO [51], lattice Boltzmann methods [52], and stochastic simulations [53].
Recently, a local spectral (LS) method, also called the discrete singular convolution (DSC) [54], was proposed as an efficient numerical tool for solving PDEs. The mathematical foundations of the DSC algorithm are the theory of distribution and wavelet analysis. The DSC algorithm has found to be useful in a number of scientific and engineering applications [55][56][57][58][59][60][61][62][63]. More recently, a family of local spectral evolution kernels (LSEKs) [64] have been constructed for analytically 4 International Journal of Biomedical Imaging integrating a class of evolution PDEs: where A, B, and C are functions of time, and Re(A) ≥ 0. The proposed local spectral method has the same level of accuracy as that of spectral methods in space, and is analytic in time. The local spectral evolution kernels are constructed by using Hermite functions. However, unlike the standard global Hermite spectral methods [65], the LSEKs adopt uniform grids and have banded matrices, just like other DSC kernels. In fact, the length of the computational stencil, and thus the accuracy of the LSEK, can be controlled according to the needs of the problem of interest, which is another feature of the DSC inherited from wavelet analysis. Moreover, the local property endows the LSEK method with sufficient flexibility to handle moderately complex geometry and complex boundary conditions. In principle, Fourier spectral methods can also solve many evolution equations analytically. However, they are restricted to periodic boundary conditions and lead to distortion near image boundaries. The objective of the present paper is to introduce an evolution-operator-based single-step image processing method. The LSEK systematically incorporates anisotropic diffusion coefficients at each pixel. The new method arrives at a desired evolution time in a single step. Moreover, the proposed method is free of instability constraint and is of spectral accuracy for integrating evolution PDEs (18). Furthermore, by an appropriate design of coefficient C(t), one can achieve the automatic termination of the time integration. In this paper, we explore the use of the proposed single-step method for several important tasks, including image deblurring, denoising, edge detection, and enhancement. Image denoising is accomplished by incorporating anisotropic diffusion type of coefficients in the LSEK, while image deblurring is realized by an appropriate forward-backward diffusion process. The coupled PDE system (11) is simplified and utilized as a model for constructing new schemes for image edge detection, and the resulting image edge is employed for image enhancement. It is believed that the proposed evolution-operatorbased method is one of the fastest schemes for image processing. In particular, the proposed PDE-based method has great potential for processing large-scale multidimensional data and for data mining.
This paper is organized as follows. In Section 2, the theory of local spectral evolution kernel (LSEK) is briefly reviewed. The frequency response of the LSEK is studied. The efficiency of the LSEK is examined by 1D and 2D diffusion processes that are exactly solvable. Numerical results are compared with those obtained by using other standard methods. A LSEK-based quasinonlinear method is proposed for image processing. Section 3 is devoted to the applications of the proposed method. Computer experiments are conducted to demonstrate the performance and efficiency of the proposed single-step method for a number of image processing tasks. This paper ends with a conclusion.

THEORY AND ALGORITHM
In this section, we give a brief review on the local spectral method and its evolution kernels. Discrete Fourier analysis is carried out to study the filter properties of the LSEK. Accuracy of the the kernel and efficiency of the scheme are examined. A new evolution-operator-based method is introduced for image processing.

Local spectral evolution kernels
The discrete singular convolution (DSC) is a general framework for constructing local spectral methods. It is an efficient approach for the numerical realization of singular convolution where η(x) is an element of the space of test functions, and T(t − x) a singular kernel. Interesting examples include singular kernels of delta type and Hilbert (Able) type. The latter plays an important role in the theory of analytic functions, processing of analytical signals, theory of linear responses, and Radon transform. The delta-type kernels are of the form where δ(x) is the delta distribution. We focus on these singular kernels since they are key elements in the approximation theory and the numerical solution of differential equations. Because of their singular nature, their approximation is necessary for numerical computations. A variety of candidates are available for the approximation in the literature. Among these examples, Shannon's delta kernel is of particular interest. It essentially gives rise to classic Fourier spectral methods. Some of local spectral kernels [54] are constructed by the regularized Shannon kernel where h is the grid spacing. The delta distribution can also be approximated by classical polynomials. Korevaar [66] proposed a Hermite function approximation to the delta distribution. Independently, Hoffman et al. [67] derived the same approximation and an expression for its arbitrary derivatives where the Hermite function h n is defined as Here H n (x) is the classic Hermite polynomial.
Yuhui Sun et al.

5
In the DSC algorithm, a function and its derivatives can be approximated by (24) where 2M + 1 is the computational bandwidth, or effective kernel support. This approach has been extensively validated by solving nonlinear PDEs [58,61], vibration analysis [59,60], Maxwell's equations [55,57], incompressible and compressible Navier-Stokes equations [62,63], and image edge detection [56]. It is also used in this work for computing derivatives.
What is crucial for the present work is an analytical integrator of the PDEs given in (18). Their formal solution can be expressed as where K(x, x , t, t ) is an evolution operator, defined as Here the standard inner product The operator can be converted into an algebraic expression by using a Fourier projection operator and then carrying out the integration over the wavenumber. The resulting kernel is diagonal in the position representation K(x, x , t, t ) = K(x − x , t, t ) and can be cast in a local spectral representation One therefore obtains [64] Here σ t t is given by Expression (28) is called a local spectral evolution kernel (LSEK). The case of A = 0 can be recovered by setting σ t t = σ. Another simplification occurs when A, B, and C are independent of time, that is, For an initial function f (x, t ), the solution f (x, t) of (18) at an arbitrary time t can be analytically given by or for a grid point x j , This implementation is the same as that for other local spectral kernels, given in (24). Equation (31), together with (28), implies that the evolution operation has been converted into an algebraic operation and the local spectral solution (31) to PDEs given in (18) is analytical. This approach is therefore free of stability constraint and is in principle of spectral accuracy when the support of the stencil is sufficiently large. For image and data analysis, it is necessary to deal with higher dimensions. If the operator in (18) its solution on an arbitrary grid point (x1 j1 , . . . , xi ji , . . . , xq jq , t) can easily be constructed by tensor product where definitions of quantities with subscripts are self-evident.
The local spectral method also provides a freedom to choose a desired level of accuracy so that the computing time is reduced. This can be achieved by an appropriate choice of M h , σ, and M. Both forward and backward time evolutions can be resolved by using the LSEK as long as Re(σ t t ) ≥ 0 for a given t. A detailed analysis and validation of the LSEK including an efficient algorithm for solving complex nonlinear PDEs via time splitting formulation are out of the scope of the present work and have been accounted elsewhere [64].

Filter properties
To achieve a better understanding of the LSEK, we study its filter properties by discrete Fourier analysis. Such a study allows us to optimally choose the DSC parameters, that is, the computational bandwidth M h and the scaled Gaussian window size r = σ/h. The Fourier analysis is a classical technique for characterizing the Fourier resolution of an interpolation or differentiation scheme applied to a class of compactly supported periodic functions and for examining the frequency response of a filter.
Note that if the function f (x, t) in (30) is regarded as a two-parameter family of signal, the convolution kernel 6 International Journal of Biomedical Imaging K h,σ (x − x k , t, t ) is a linear time-invariant system and is practically used as a finite impulse response (FIR) filter. This filter is stable and inevitable for C < ∞. On the one hand, when B = C = 0 and A(t) is a constant, (18) reduces to the heat equation. As the LESK provides analytical solution to the heat equation, one might assume that the analytic integrator is a Gaussian filter as the action of heat equation is known to be equivalent to the Gaussian filter [2]. However, this is not true. In fact, it looks like a wavelet filter in the coordinate domain, see Figure 1. Mathematically, it is called a kernel of Dirichlet type in contrast to the kernel of positive type, for example, a Gaussian kernel. On the other hand, when B = 0, the kernel K h,σ (x, t, t ) is obviously symmetric and its Fourier response is of linear phase and real, since it contains only even terms of the Hermite functions. As it is also the evolution kernel for the parabolic PDE (when A = 0), it must be a lowpass filter in general for t > t as shown Figure 2. It is noted that the LSEK lowpass is almost ideal for small t, and appropriate M h and r values. In general, large M h and r values lead to a better approximation to the ideal lowpass filter. Nevertheless, since they are constructed from polynomial approximation, they are Chebyshev type of filters as illustrated with appropriate M h and r values, see Figure 3. Moreover, for a given set of other parameters, when t − t is getting larger, that is, long time evolution, the LSEK filter response becomes more and more focused in the low-frequency part, corresponding to large amount of diffusion.
In general, the parameter B gives rise to a translation in the coordinate, which leads to an extra phase factor in the frequency response, see Figure 4. Obviously, in signal analysis, it is the time-shift factor that is useful for the circuit design. In other words, the drafting term B(t)(∂/∂x) in (18) is a timeshift operator from the point view of signal processing. Since B(t) is a function of time, it in fact can be used for nonlinear time shifting. As shown in Figure 5, the parameter C provides an overall scaling factor (either exponential growth or decay) to the amplitude of the signal, which can be useful in practical applications.

Numerical test
We illustrate the proposed local spectral method for the diffusion equation in this section. The diffusion equation is a special case of the class of PDEs of (18) (B = 0 and C = 0). The most important fact is that the diffusion equation has extensive applications in image processing. The accuracy of the proposed local spectral method is obtained by comparing its results with exact solutions. Two advantages of the present local spectral approach (accurate and fast) have been fully demonstrated since both forward and backward time evolutions can be resolved by using the LSEK as long as Re(σ t t ) ≥ 0 for a given t. When A t t is positive in (29), it is obviously a forward propagation process. In this case, the numerical solution of the forward diffusion equation at any time can be obtained by (31) in only one time step. However, when time A t t is negative, the problem becomes a backward propagating process, and we still can use (31) to get the solution as long as we ensure Re(σ t t ) ≥ 0. However, the backward diffusion is an unsteady problem. In order to obtain a reliable numerical solution, we have to terminate the backward diffusion Yuhui Sun et al.
in a limited time. As was mentioned earlier, the accuracy of the proposed local spectral method is controllable by adjusting its support parameter M and the Hermite parameters M h and σ. In this numerical example, we choose the Hermite parameters and the computational bandwidth as M h = 128, σ = 3.7Δ, and M = 45 to illustrate the spectral order of accuracy of the method. Unless stated otherwise, the above set of parameters is used in this section. In Example 2 (2D diffusion process), we will compare the performance of the present local spectral method with that of Fourier spectral method which is realized by using the fast Fourier transformation (FFT).

Example 1.
We consider a one-dimensional linear diffusion process given by where α is the diffusion coefficient. With an initial deltafunction distribution localized at x 0 , the analytic solution is The computation is conducted using a sufficiently large interval [−10, 10] of coordinate space to ensure that boundary reflection is negligible. Fifty grid points are used for this interval. The diffusion coefficient α is set to be 0.5 in all computations. Since it is difficult to accurately represent a delta function on a grid as the initial condition, we choose the initial value as the Gaussian of (34) at t = 1.0 for the forward diffusion process and the Gaussian of (34) at t = 2.0 8 International Journal of Biomedical Imaging as the initial condition for the backward diffusion process. The numerical solutions at time t = 1.0, 3.0, 5.0, 7.0, and 9.0 are given in Figure 6. Two kinds of numerical error measures, that is, L ∞ and L 2 are used to evaluate the quality of the present local spectral method. The computational errors are listed in Table 1. It can be observed from Table 1 that the present local spectral method gives highly accurate numerical results. On the other hand, it can be seen clearly that the results of the backward diffusion are not as good as those of the forward diffusion. The loss of accuracy is due to the fact that the backward diffusion is itself an unsteady pro-cess. Moreover, since no time evolution scheme is needed in the present method, there is no accumulation of time discretization error. This property guarantees the accuracy and the speed of the present method and shows great advantage over many other numerical algorithms.

Example 2.
We consider a two-dimensional linear diffusion equation It is a good test problem because it admits an exact solution of the form where a is a positive parameter and we will use a = 0.7 in this test. To compare with the FFT method, we impose periodic boundary conditions in both xand y-directions. The computation domain is taken to be [0, 10] × [0, 10] with equal spatial discretization. Although the grid can be set to be any integer values for the LSEK method, for the purpose of comparing with the FFT method in which the grid is required to be the power of 2, we choose N x = N y = 32 in our computation. The forward diffusion is initialized to be a Gaussian wave at t = 0.1, where the point (x 0 , y 0 ) is the center of the 2D square domain. The detailed LSEK method used for the 2D problem is given in (33). Our interest is to explore the efficiency of the present LS method. Therefore, we will not only compare the numerical errors of the LSEK method and the FFT method but also the computing time of both methods. They are given Yuhui Sun et al.

9
in Table 2. It is clear from Table 2 that the present LSEK method is as accurate and fast as the FFT method when they are solving the diffusion equations. By comparing the present local spectral method and the Fourier spectral method, we reach conclusions which are summarized as follows.
(i) The accuracy of the proposed LSEK method is controlled by its computational bandwidth M, Hermite parameters M h and σ. Under choice of M = 45, M h =128, and σ=3.7Δ, the LS method can provide the same level of accuracy as that of Fourier spectral method.
(ii) The global nature of the Fourier spectral method makes it difficult for practical applications to problems of complicated boundary condition and complex geometry. However, the present LSEK is a local method which can be used for treating complicated boundary conditions and geometries.
(iii) For the Fourier spectral method, the number of grid points is required to be the power of 2, but there is no such a limitation for our LS method.
(iv) The present LSEK method transforms the diffusion equation into algebraic calculations. This implies that the LSEK method can speed up the computation dramatically, and thereby reduce the computational cost. Specifically, the CPU cost of the LSEK scales as O(MN), where M and N are half of computational bandwidth and the number of grid points to deal with. However, asymptotically, the LSEK method requires less CPU time than the Fourier pseudospectral method does as the latter scales as O (N log N).
In practical applications, especially in image processing, the aforementioned precision is entirely unnecessary. A three-percent error does not produce much difference in visual perception. Therefore, in the following image experiments, we choose M h = 36, σ = 1.7Δ, and M = 8. In fact, even less costing parameters can be used for certain tasks in image processing, such as image deblurring and denoising where the image quality is very low to start with.

Adaptation of anisotropy
Obviously, the Perona-Malik type of anisotropic diffusion equations cannot be solved directly with the LSEK in a single step, albeit the local spectral method given in (24) can be used to solve these equations iteratively. Fortunately, a close examination on the solution scheme, (31) or (33), and the LSEK (28), one realizes that the solution scheme is of collocation type and the kernel allows a localized modification of coefficients A, B, and C, that is, This flexibility of selecting local coefficients endows us with a new single-step evolution-operator-based method for image processing and data analysis.
In Perona-Malik equation (7), the break up of the operator leads to a gradient controlled diffusion term and a gradient controlled drafting term. The local strength of the both  terms can be easily computed and incorporated in the LSEK expression.
Since the diffusion coefficients computed from Perona-Malik type of equations may vary from point to point, possible N 2 sets of LSEK weights may be computed. This requires much computation if N is about 512. On the other hand, since image resolution is normally limited to 8 bits, it is obviously unnecessary to use more than 256 sets of LSEK weights. We therefore classify the diffusion coefficients into a total of N c groups. For each group, we calculate a set of LSEK weights according to the mean value of the group. In most practical applications, taking N c ∼ 20 is enough and it costs less than 1 second to compute these weights. In the above argument, we have assumed that the drafting coefficient behaves similar to the diffusion coefficient at each point. If this is not the case, the discussed classification method can be easily modified. Thus, the basic procedure for the proposed evolutionoperator-based image processing method is the following.
where K h,σ,Ac,Bc (k x h, t) is obtained from (28) by setting t = 0, A t t = A c t, and B t t = B c t. Although the proposed method has its root in nonlinear PDEs, it does not exactly solve the anisotropic diffusion equation. Instead, it is constructed by the pointwise adaptation of anisotropy to the evolution kernel for the isotropic diffusion equation. Its major advantages are its single-step operation, freeness of numerical instability, readiness for multidimensional problems. Moreover, by appropriately choosing B and C, systematic time shift and amplitude scaling can be easily accomplished in a single step.

APPLICATIONS
In this section, we apply the evolution-operator-based singlestep approach for image deblurring, denoising, image edge detection, and image enhancement. The performance of the proposed method is illustrated by a large number of standard test images and is compared with that of other existing methods.

Image deblurring
Image blur refers to serious degradation to image qualities by either the environment of an imaging object, or the process of imaging, or the processing of an image [1]. Backgrounds of imaging objects, such as fog, air pollution, texture to a fingerprint residual and illumination, are common environmental sources of image blur. Motions of the imaging object and/or camera, or out-of-focus lens generate blur in the process of imaging. Jiang and Wang [68,69] discussed image blurring due to reconstruction techniques. Very often, image processing may lead to blur during a variety of mathematical operations, such as Gaussian convolution, interpolation, or image registration. Deblurring is a basic image processing task and the subject is so extensively studied that it is virtually impossible to give a comprehensive review to include all important contributions and developments in the literature.
Essentially, deblurring methods can be classified as linear and nonlinear, global or local, time-dependent or timeindependent, a priori knowledge based or blind, and so forth. Linear filters [70][71][72][73] do not make use of local information of an image under processing, and are usually simple and low cost. Many linear filers, such as Wiener filter, Hanning filter, and Fourier filtering, work effectively when there is a priori knowledge about the source of the blur. Moreover, linear filters are often utilized to obtain initial information for other nonlinear filters. Nonlinear filters make use of the local information of an image to gain a better effect of deblurring. Some nonlinear filters assume partially a priori knowledge of an image, such as Kalman filter and Lucy-Richardson algorithm. The former makes optimal use of imprecise data on a linear (or nearly linear) system with Gaussian errors to continuously update the best estimate of system's current state, while the latter maximizes the likelihood that the resulting image is an instance of the blurred image, assuming the Poisson noise statistics. Many blind deblurring algorithms can be used effectively when no information about the source or sta-tus of the distortion (blurring and noise) is known. They utilize the local information or Fourier spectrum of an image intensively and achieve the goal of deblurring by iterations. Obviously, Perona-Malik type of PDE-based image deblurring algorithms are time-dependent, nonlinear, blind filters. Wang et al. proposed an efficient blind deblurring algorithm for spiral computed tomography (CT) images [74]. While wavelet transform is a linear algorithm in nature [75], its application to image processing can be made nonlinear by appropriately incorporating image information into the selection of truncation parameters [76][77][78][79]. Recently, some timeindependent total-variation (TV-) based regularization functionals have attracted much attention because of their ability of recovering image edges during deblurring, see Rudin et al. [14], Li and Santosa [80], Dobson and Vogel [81], and Chan et al. [82]. These methods are often formulated in terms of Euler-Lagrange equations and are associated with ill-posed inverse problems. Ideally, TV-based methods should work as blind algorithms. Nevertheless, the main difficulty associated with TV methods, as indicated by many researchers, is the determination of the unknown regularization parameter(s), which balances the fidelity and smoothness of a TV solution. Consequently, these methods are most efficient when some a priori knowledge about the image is given.
A common used degradation process is often presented in terms of a blur operation and additive noise where η stands for a white additive Gaussian noise and where R is a linear operator representing the blur (usually a convolution). Given degradation version of image I 0 , the restoration problem is then to obtain I knowing (41). In this study, we will concern ourselves only with two aspects of image restoration which are deblurring and denoising. We also make use of the sharpening property of the backward diffusion to restore a blurry image. The effect of noise amplification (an inherent byproduct of backward diffusion process) will be minimized by terminating the diffusion process in a limited time or by choosing an appropriate diffusion coefficient. We consider the anisotropic diffusion equations, as a model. Here n is a unit vector, outward normal to the boundary ∂Ω. The diffusion coefficient A is a function of the smoothed gradient magnitude s = G σ * ∇I evaluated at t = 0. We further create the LSEK analogy of (42), which is given by (40) with B = C = 0. Here, we design a smoothed gradient-magnitude-dependent diffusion coefficient for deblurring: Yuhui Sun et al. where is a small strictly positive real number for avoiding zero denominator and k is a threshold parameter like that in the Peronal-Malik model. Here k can be a constant threshold or be adjusted locally by gradient. In the latter case, the parameter k varies gradually with the image, and edge sharpening is accomplished by different thresholds in different locations. As an illustration we present in Figure 7 a plot of the suggested diffusion coefficient A with k = 80. It is seen that the diffusion coefficient considered in this work is strictly increasing and negative. Negative diffusion coefficient is for edge sharpening and increasing tendency ensures less backward diffusion for large gradient magnitudes.
We use the parrot image for our demonstration. The blurry parrot image is obtained by convolving the initial parrot image with a Gaussian kernel (σ = 2). As for another Gaussian kernel used for the calculation of the smoothed gradient G σ * ∇I of original image, we set σ = 3. We stop the diffusion process at time t = 1.0. Figure 8 shows the results of deblurring. For a more clear observation, we zoomed Figure 8 out near one parrot's head and the zoomed images are shown in Figure 9. We found that edges at different locations have been improved. For example, the edges near parrot's beak and those around the eye are well sharpened. Figure 10 shows the deblurring process for a building. The blurring and deblurring processes are similar to those described in the parrot image. In this case, we are interested in the features that characterize the architecture, such as edges and lines.
In fact, one can interpret the above deblurring process as an adaptive resolution of the linear diffusion equation which is conceptually different from the conventional nonlinear diffusion filtering in which the nonlinearity is spatial.

Image denoising
Noise is one of the most common sources of image degradation [83]. It is often referred to as the rapid changes over spatial extension in image intensities or color plan. Either imaging acquisition process or naturally occurring phenomena can lead to noise contamination. Typically, noise is assumed to be additive and satisfies a mathematical model, such as the Poisson, or the Gaussian, or Laplacian. Some noise is multiplicative, for example, the speckle noise. The image denoising problem refers to the process of recovering an image contaminated by noise [84] and is one of the most elementary operations in image processing.
The earlier algorithms include (linear) spatial filtering approach [73] and frequency cut-off filters in the Fourier domain, which are first motivated by the primitive ideas from signal processing. These filters reduce the noise by removing the high-frequency components. However, they also blur images. Stochastic modeling [70][71][72] originated from estimation theory has also been introduced for noise removing. In the past two decades, a number of new techniques have been introduced to improve the quality of image denoising. Statistical approaches such as maximum-likelihood estimation and adaptive (Wiener) filters are developed. Local adaptive filters are combined with impulse removal filters in the transform domain to remove not only white and mixed noise, but also their mixtures [85]. Wavelet theory [75] has been used for denoising by many researchers [78,[86][87][88][89][90][91]. The essential idea in wavelet schemes is to decompose the noisy image into a number of frequency bands, and remove the noise in appropriate frequency bands. Wavelet approaches usually perform well with a viable strategy to discriminate noise from image in each wavelet subdomain. PDE-based denoising methods were proposed to preserve image edges during the process of noise reduction as briefly reviewed in the introduction.
If impressive results are the main reasons for using nonlinear diffusion filtering in image processing, poor efficiency is the main reason for less widely spread application of nonlinear diffusion filtering in industrial codes. Finitedifference schemes are extensively used to discretize nonlinear anisotropic diffusion equations and they make the whole filtering procedure quite time consuming. We resolve this problem by using the LSEK. The efficiency of the proposed LSEK for diffusion equations has been already validated in Section 2. In the present work, we explore the performance of the LSEK image denoising. We take two different LSEK approaches.
In our first denoising approach, we choose a timeindependent diffusion coefficient of the Gaussian suggested by Perona and Malik: where the variance σ is based on a local statistical estimate of (8). By using this variance, no additional kernel smoothing is needed. The success of this approach has been fully demonstrated in [20]. In this study, we used the proposed evolution-operator-based single-step method to obtain a fast denoising process. A Gaussian noise contaminated Lena image has been used to test the performance. It is found that the PSNR of the noisy Lena image is 21.03 db. After denoising, the quality of the Lena image has been improved and the resulting PSNR is 29.81 db. Both noisy and denoised Lena images were presented in Figure 11. We also examine our scheme for an MRI image of brain. An improvement of 6.30 db in the PSNR is obtained for this medical image, see Figure 12.
In the second approach, we use (40) directly for image denoising in a single step. From (29) and (28), it is seen that the diffusion coefficient A t t = t t A(τ)dτ controls the local diffusion. Therefore, one can tailor A t t for a given purpose. In this work, we propose the following gradient-dependent A t 0 : where s = |G σ * ∇I| is the smoothed gradient magnitude, β and are two parameters controlling the tendency of A t 0 . We fix β = 16 and vary between 0 and 0.1. Figure 13 shows Yuhui Sun et al. the change of A t 0 with sets of = 0, 0.002, 0.01, and 0.1 at the time t = 100. Obviously, = 0 almost corresponds to the linear diffusion process. Therefore, we set to be larger than zero in order to preserve edges of images. In our test, we take = 0.01. With β and fixed, we plot the relation between A t 0 and the time t in Figure 14. One can observe that A t 0 basically does not change any more as the time t is large enough. This implies that the denoising process is stable with respect to the stopping time. Thus, the stopping time can be chosen to be quite large without worrying about possible oversmoothing.
We classify A t 0 by equally dividing the values of A t 0 into N c = 20 groups. We therefore generate 20 sets of LSEK coefficients accordingly. We still use the last two noisy images for our denoising demonstration, see Figure 15. The PSNR for the restored Lena image is 29.41 db and that for the restored brain image is 25.22 db. It is can be seen that the performance of this denoising procedure with a time-dependent diffusion coefficient is almost the same as that of the time-independent one. As for the stopping time, we can set it to be 5, 10, or a larger number, it is not very sensitive.
Obviously, more efficient diffusion coefficients can be designed to remove blur and noise from images, but they are beyond the scope of the present work. What we want to convey to the reader is just an efficient single-step procedure for image deblurring and denoising.

Image edge detection
Edges characterize object boundaries and are therefore crucial for image segmentation and registration, and for the  identification of objects in scenes. In recent years, considerable research has been carried out in the development and application of efficient algorithms for the detection of the image edges. A well-known approach in the edge detection is to associate edges with zero crossing of second-order derivative. The resulting operator in use is called Laplacian of Gaussian (LoG) operator [92]. An alternative approach is the method of curve fitting [93,94], in which an approximated function is built based on the least-square method. Therefore, a derivative on the fitting curve yields the approximate differentiation for each pixel directly. Compared with the finite-difference method, curve fitting approach results in better performance in a noisy environment. They are, however, computationally more expensive. Recently, Canny [95] has formulated edge detection as an optimization problem and defined an optimal filter which can be efficiently approximated by the first derivative of the Gaussian function. With a recursive procedure [96], the Canny detector provided an efficient way for image noise filtering and edge detection. Other edges detection methods include differentiation-based methods that use the logarithmic image processing (LIP) models [97], contrast-based methods [98], and relaxation labeling techniques [99]. In fact, these delicate techniques can achieve better performance in one way or another for common tasks. However, for images with large amount of texture, these traditional edge detection techniques are no longer reliable because of severe edge distortion. The latter is due to the wrong prediction of high-frequency edge components by low-accuracy methods. Most recently, Wei and Jia [35] proposed a novel synchronization-based realistic edge detection scheme that has demonstrated a success in detecting edges of high-texture images. In this scheme, the residual of synchronization, defined by the difference in the synchronization of two dynamic systems, gives rise to image edges as if obtained by a highpass filter. However, solving two dynamic systems makes it computationally expensive. In this work, we propose a simplified version of this synchronization scheme.
The essential idea of this new edge detector is quite simple. A smoothed version (SV) of the original image I is obtained by modifying the forward nonlinear diffusion operator (42). In such a nonlinear operator, the diffusion coefficient at image edges is significantly large, with contrast to the original use of nonlinear diffusion operators. As a result, image values near the edges are effectively suppressed. This smoothing operation, from the point of view of image processing, is a lowpass filter. Whereas, image edge detection should be a highpass filtering process. Therefore, edge detection version (EDV) denoted by I EDV (x, y) is then expressed as the difference between the original image and the smoothed image: It is noted that image edges are effectively emphasized in I EDV (x, y) due to the edge suppressive nonlinear diffusion operator. Here C 1 and C 2 are two constants, either 1 or −1.
We emphasize that the sign of C 1 and C 2 must be the same.
If the edge of the object has a white color on a black background, we take C 1 = C 2 = 1. On the contrary, if the edge of the object has a black color on a white background, we choose C 1 = C 2 = −1. For instance, the tiger image is composed of the black object (a tiger) on the white background, C 1 and C 2 are, therefore, set to be −1. While the pinecone image is composed of the white object on the black background, we take C 1 = C 2 = 1. Computational cost is very low in the present approach since only a single step is needed in running the nonlinear diffusion equation with our LSEK, see Table 3. Nevertheless, one can further speed up the edge detection by using a linear diffusion operator, for example, running the forward diffusion (42) till t = 3 by our LSEK, with a constant diffusion coefficient A = 1. In order to assess the computational complexity of the LSEK method, we compared in Table 4 the number of basic operations for four different methods in a single time step. These methods include the adaptive resolution scheme [41], the Perona-Malik nonlinear diffusion equation [9] resolved by the conventional explicit scheme, the AOS scheme [18] without the step of regularization, the LSEK for a linear diffusion process. Table 4 provides a quantitative analysis of computational complexity for these schemes. Multiplication or division and the addition or subtraction are calculated separately, as given in [18,41]. In Table 4, N denotes the number of pixels to treat, m is the dimension of the problem considered, and W is the half of the computational bandwidth. Note that the first two methods require multiple time steps, while the last two methods can be integrated in a single time step. However, the reliability of the implicit AOS scheme for a large time step needs to be validated.
We next discuss an efficient edge extraction scheme. As it is well known, thresholding is one of the most common approaches. Usually, thresholding may be viewed as an operation that tests a function T of the form where p(x, y) denotes some local property of the point. Fortunately, experience has shown that the threshold T is connected to the estimated variance which is a measure of contrast. The choice of variance can be where σ 2 denotes the local variance. Here, I EDVi,j represents local average and is defined as where K and L are the support of subimage in xand ydirections.
The edge detector is susceptible to streaking if it uses only a single threshold. Streaking is the breaking up of an edge contour caused by fluctuation around the value of the threshold. Therefore, a more effective threshold scheme uses two thresholds T high and T low with T high ≈ 2T low . In this case, streaking occurs only when it fluctuates above the high threshold and below the low threshold. Therefore, the probability of streaking is greatly reduced. The probability of isolated false edge points is also reduced because the strength of such points must be above the high threshold. We note that if the threshold T low is too small, there will be some false edges. Similarly, if T high is too large, some portion of real contour may be missing. In our scheme, the ratio of the variance of I EDV to that of original image is used to determine T low , from which T high can be fixed, that is, where σ EDI and σ I represent the variance of the edge detection image and that of the original image, respectively. In addition, nonmaximum suppression is performed to thin the edges in the process of edge extraction. For more detailed description, the reader is referred to [95].
To demonstrate the capabilities of the edge detector described above, we carry out experiments on various synthetic and real images, and compare the results with the output of four other detectors: Sobel detector, Prewitt detector, Canny detector, and anisotropic diffusion scheme (Perona-Malik model).
A variety of images are employed for the edge detection. We first consider the Barbara image, which is a challenging test. Figures 16(a), 16(b), 16(c), 16(d), 16(e), and 16(f) show the Barbara image and edges of the Barbara image detected by the Sobel detector, the Prewitt detector, the Canny detector, the anisotropic diffusion scheme, and the present LSEK edge detection method. It is seen that the present LSEK edge detection method is the best. Facial edges are accurately extracted and image texture is clearly detected. However, none of the other four detectors gives correct edge response to the textured part of the image. The regularity of the thin lines in the original image has been entirely distorted.
We next consider a few standard texture images (baboon, tiger, Goldhill, bridge, boats, and pinecone), see Figure 17. In general, the processing of texture images is a difficult task for most existing edge detection methods because the involvement of high-frequency components. For example, most existing methods are not able to distinguish baboon's bristles from the background. Moreover, it is a challenge to fully detect the thin ropes on the boats by many other existing methods. Similar fine texture, such as that in the grass and trees can be found in other four images. The ability of resolving such details is crucial to edge detection and pattern recognition. It is seen from Figure 18 that the proposed LSEK does an excellent job in the edge detection of these high-texture images. The detected image edges provide even clearer features than the original images do. For example, the pattern on the house roof and window in Goldhill image is not very clear in the original images. This feature is visible from the detected edges.
Having validated the proposed method for the edge detection of texture images, we finally apply it to a group of medical images as shown in Figure 19. The edge detection of X-ray images can be difficult because of small gradient in image resulted from small variation in X-ray attenuation, see Lung01, Lung02, and Lung04. Large amount of vessels can lead to high-texture images as shown in Lung03. Figure 20 indicates that the proposed method is able to reveal a variety of features in these images.
Based on the above discussion, it is evident that the present LSEK edge detection method is one of the best edge detectors in performance. Moreover, the present method requires only a single-time-step operation. It is these two advantages, high efficiency and excellent performance, that qualify the proposed LSEK edge detection method as one of the best approaches for image edge detection.

Digital image enhancement
In general, image enhancement is a process designed to increase the usefulness of the image. When images are enhanced for human viewers, the objective is to improve perceptual aspects, including image quality, intelligibility, or visual appearance. One important class of image enhancement is to modify the contrast and/or dynamic range of the image. Another class of image enhancement is to reduce the degradation of the image. In this work, we focus on the former class of issues, that is, to improve the contrast of mammograms.
A large number of investigations have been made to enhance mammography feature while reducing noise. Gordon and Rangayyan [100] used an adaptive neighborhood image processing technique to enhance the contrast of features relevant to mammography. However, this method enhanced the contrast of mammography features as well as noise. Dhawan et al. [101][102][103] developed an adaptive neighborhood-based image processing technique following Gordon's work. They utilized low-level analysis and knowledge in desired feature design to improve the contrast of specific features.
Filtering as an oldest technique has also been extensively used in medical imaging. For example, Tahoces et al. [104] developed a method for the enhancement of chest and breast radiographies by an automatic spatial filtering. In their work, a linear combination of an original image and two smoothed images were used. The process was completed by a nonlinear contrast stretching.
The difference image is a commonly used enhancement technique and is the basis of most existing computer-aided diagnosis systems. The difference image is usually obtained by subtracting a suppressed image from the original one to remove the background [105].
As we have discussed in the last subsection, the edge of an image is obtained based on the difference between 18 International Journal of Biomedical Imaging where the parameter C 5 is used to control the maxima in scale-space and C 3 is a magnification parameter. Basically, C 5 is an element of (1, 2) and C 3 is an element of (1, 2). In the above equation, the smoothed version of the image is still produced by solving the linear forward diffusion equation (42) which can be rapidly accomplished by the LSEK. The diffusion coefficient is also fixed to be A = 1 and the whole diffusion process is terminated at time t = 10.
We first consider four mammography images collected from the FTP site ftp://peipa.essex.ac.uk/ipa/pix/mias. The original images from the FTP site all have the size of 1024 × 1024 pixels. However, the images used in our study are cropped to the useful portion, as shown in Figure 21. The above image enhancement scheme employs two important parameters C 5 and C 3 , which are set to 1.95 and 1.05, respectively. Figure 22 shows good results for almost all test images, that is, abnormalities on all mammography images can be seen more clearly in the enhanced images than in their original ones. We next consider four CT images of a mouse as shown in Figure 23. The same parameters used in the last example are utilized for these images. Results shown in Figure 24 indicate that the proposed method is very efficient for medical image enhancement. Nevertheless, it is worth mentioning that more attractive image enhancement results may be obtained if one chooses the diffusion equation coefficient adaptively.

CONCLUSION
Partial-differential-equation-(PDE-) based approaches have great potential for the processing of images and multidimensional data because they can be made to be systematic and automatic for large-scale high-dimensional image data.
Yuhui Sun et al. Nevertheless, research in this direction is often hindered by problems in parameter optimization and extra computing time. The present work proposes an evolution-operatorbased single-time-step method for image and multidimensional data processing. By utilizing a local spectral evolution kernel (LSEK) that analytically integrates a class PDEs, we show that a number of image processing operations, such as image deblurring, denoising, edge detection, and enhancement, can be effectively carried out in a single step of time integration. Filter properties of the LSEK are studied by using the Fourier analysis. Local information about the image function and its gradient is analyzed and incorporated in the time evolution of the anisotropic-type equation. In image denoising, a local-variance-based method is employed and is compared to the regularization type of approaches. Automatic exit schemes are proposed to stop image evolution in the denoising. We have achieved about 8 dB in the peak signal-to-noise improvement. In image edge detection,

22
International Journal of Biomedical Imaging

24
International Journal of Biomedical Imaging the proposed method utilized a simplified version of the previous coupled PDE algorithm. The present method is compared with a number of standard approaches, including Sobel, Prewitt, Canny detectors, and anisotropic diffusion operators. High-texture images and X-ray lung images are employed to demonstrate the proposed method. It is found that the proposed scheme provides some of the best performances for the edge detection of high-texture images, which is a challenge to other existing methods. Image enhancement is carried out by the combination of the image edges and the original image. Significant enhancement is achieved due to reliable image textures.

ACKNOWLEDGMENT
This work was supported in part by NSF Grant IIS-0430987 and by IRGP Grant 71-4834.