Mode decomposition evolution equations

II.A High order PDE based low-pass filters

Motivated by a number of physical phenomena, such as pattern formation in alloys, glasses, polymer, combustion and biological systems, Wei introduced some of the first family of high order nonlinear PDEs for image processing in 1999⁷⁴

where u(r, t) is the image function. Here d_q(u(r), |∇u(r)|, t) and e(u(r), |∇u(r)|, t) are edge sensitive diffusion coefficients and enhancement operator respectively. Equation (1) is a generalization of the Perona-Malik equation,⁵⁴ where the latter is recovered at q = 0 and e(u(r), |∇u(r)|, t) = 0. The hyper diffusion coefficients d_q(u(r), |∇u(r)|, t) were chosen as

where the values of constant d_q0 depend on the noise level, and σ₀ and σ₁ were chosen as the local statistical variance of u and ∇u

The notation $\overline{X(\mathbf{r})}$ denotes the local average of X(r) centered at position r. The statistical measure based on the variance is important for discriminating image edges from noise. By using this measure, one can avoid the use of preprocessing, i.e., the convolution of the noise image with a smooth mask. In general, both Neumann type or periodic boundary conditions can be imposed. However, if arbitrarily high order PDEs are devised, periodic boundary conditions are better choices, as are chosen for the application with Fourier pseudospectral methods in this paper. The numerical application of Eq. (1) to image denoising and restoration was demonstrated and the performance of this equation was compared with that of the Perona-Malik equation⁵⁴ by Wei⁷⁴ and many other researchers.^{22, 23, 44} The well-posedness of the generalized Perona-Malik equation was analyzed in terms of the existence and uniqueness of the solution.^{7, 24, 25, 37, 82} It was argued that the mathematical properties of the generalized Perona-Malik equation differ from those of other high order PDEs.³⁷ The stability of Eq. (1) follows from appropriate choice of diffusion coefficients d_q(u(r), |∇u(r)|, t), e.g., the sign of d_q should be (−1)^q. However, other choices of the signs were also explored in forward-and-backward diffusion processes to achieve simultaneous image denoising and enhancement.²²

II.B Nonlinear PDE based high-pass filters

By 2000, there had been a large number of publications that explored PDE based image processing approaches. However, most of the research was focused on the use of even-order partial derivatives as a means to smooth images, i.e., as low-pass filters. These PDE approaches were inefficient for image edge detection, particularly for images with a large amount of texture. This is due to the fact that edge detection is a high-pass filtering operation, while the diffusion process is inherently a low-pass filtering process. Wei and Jia addressed this issue by introducing a pair of weakly coupled nonlinear evolution equations in 2002⁷⁷

where u(r, t) and v(r, t) are scalar fields with Neumann boundary conditions imposed, ε_u and ε_v are coupling strengths. F₁ and F₂ are general nonlinear functions which can be chosen as the Perona-Malik operator

with d_j(|∇u_j|) = d_j0 exp[−|∇u_j|²/(2σ²)], u₁ = u and u₂ = v. The initial values for both scalar fields were chosen to be the same image of interest, i.e., u(r, 0) = v(r, 0) = I(r). In the theory of nonlinear dynamics, Eqs. (4) and (5) constitute a synchronization system. To attain an appropriate image edge contrast, the original image must evolve under dramatically different time scales, i.e., either d₁₀ ≫ d₂₀ or d₂₀ ≫ d₁₀. The coupling strengths ε₁ and ε₂ were set to be relatively small (i.e., ε_u ≅ ε_v ~ 0) so that the rate of change of u or v was dominated by the diffusion process. The coupling term plays the role of relative fidelity. This becomes clear when d₁₀ ~ 0. The image edge was defined as the difference of two dynamical systems

To ensure a normal performance of diffusion operator, one of the diffusion coefficients d₁₀ or d₂₀ must be of similar amplitude as that used in the Perona-Malik dynamics. If we choose d₂₀ to be a normal one, the requirement of d₂₀ ≫ d₁₀ implies that d₁₀ ≈ 0. As such, w(r, t) in Eq. (7) behaves like a band pass filter. In the extreme case, if we set d₁₀ = 0, we have a PDE based high-pass filter

This setting is often used in our practical applications.^{64, 65} The PDE based high-pass filters have been shown to be very robust and efficient. They provide 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.⁷⁷

II.C Mode decomposition evolution equations (MoDEEs)

The high-pass PDE filter discussed above does not discriminate different high frequency modes or components. The goal of the proposed project is to devise a PDE based system for systematic separation of all the mode components, including high order modes. This requires a generalization of the PDE based high-pass PDE filter. To this end, we first cast Eqs. (4) and (5) into a matrix form

where all quantities are same as those defined in the last section. Considering the fact that w ~ u − v, we can arrive at a more compact form

Since ε_u ≅ ε_v ~ 0, for simplicity, we drop the coupling strengths (ε_u + ε₁) in the evolution equation of w

where we have relabled variables v and w as v₁ and w₁, respectively. Note that the PDE for “edge” w in Eq. (11) is very different from the other two genuine PDEs for u and v. The solution of w is easily obtained by integrating the corresponding PDE once only, provided that u and v are given. Equation system (11) is hereafter referred to as the first MoDEE system.

For problems involving high frequency modes, we need to construct the high order MoDEEs. The matrix form (11) can be easily generalized to include fourth order derivatives

where differential equations of u and v₂ employ the fourth order nonlinear low-pass filters given in Eq. (1). Here w₂ is more sensitive to high frequency components than w₁ is. Therefore, w₂ contains the information related to the difference between high order modes. As discussed earlier, it is computationally cheaper to solve edge functions w₂ than to solve u and v₂. We have utilized the fourth order operator designed by Wei⁷⁴ in the above formulation. The initial values of u and v₂ are set to the original image or signal of interest. The initial value of w₂ can be set to zero. Equation system (12) is hereafter referred to as the second MoDEE system.

If the task is merely to detect image edges, one may just use the first MoDEE system or the high-pass filter of Wei and Jia to obtain desirable results as demonstrated in our earlier work.^{64, 77} However, if one wants to separate high frequency modes, the second MoDEE system is a better candidate than the first one.

For highly oscillatory signals, images and functions, in order to separate all the high frequency modes, which is desirable in many applications, we need the MoDEEs containing even higher order PDEs. To this end, we construct an nth order system

Equation (13) can be solved to generate nth edge functions w_n. Similar to the first MoDEE system, the nth MoDEE system works when one sets d_uj+1 ~0, j = 0, · · ·, n − 1, ε_u ~ 0, and ε_n ~ 0. Since all PDEs are designated as smoothing operators, all odd diffusion coefficients, such as d_u1, d_u3, · · ·, d_n1, d_n3, · · · are to be positive and all even diffusion coefficients are to be negative. In order to separate high order modes, it is practical to set all d_uj+1 = 0 and d_nj+1 = 0, except for d_un and d_nn.

III.A Linear MoDEEs

In general, we need to solve systems of coupled, nonlinear, high order PDEs as shown in Eq. (13), which can be technically demanding. However, for many practical problems, we can design simplified MoDEE models which also work well. Based on the physical understanding, various simplifications of the MoDEE systems are presented to reduce the computational complexity. Numerical algorithms for solving the high order MoDEEs are also discussed in the following paragraphs.

The first approximation is to set ε_j = 0 in the MoDEEs. This approximation leads to the decoupling of all MoDEEs. The uncoupled MoDEEs are much easier to solve. The solutions of decoupled PDEs are also easy to interpret, which makes the selection of mode parameters easy.

Another useful approximation is to eliminate nonlinearity and positional dependence of all the diffusion coefficients. The nonlinear models work better for most application as evident from our own experiments. However, for certain class of problems with smooth functions, the nonlinear edge sensitive diffusion coefficients do not play the critical role. Signal processing is a typical example. Therefore, the linear MoDEE systems are computationally favored. As an example, we can obtain a linear form of the nth MoDEE system from Eq. (13):

The linear MoDEE systems can be easily solved with the Fourier pseudospectral methods.^{65, 84}

III.B Algebraic MoDEEs

To further simplify the MoDEE systems, we can make use of the algebraic relation given in the original formulation of Wei and Jia.⁷⁷ The use of algebraic relations saves the computational cost of integrating the edge equations. We first consider a set of linear equations modified from Eq. (14)

These linear equations offer copies of smoothed image functions. After solving Eq. (15) for u and v_n, we determine the mode functions by the simple form of Wei and Jia as that given in Eq. (7)

By the appropriate selection of parameters, w_n(r, t) delivers desirable mode functions. It behaves as a band-pass filter.

Another simple choice is to determine the mode functions by the form given in Eq. (8)

where I(r) is an initial value. As such, we do not need to solve the PDE for u. It behaves as a high-pass filter. We expect that w_n in Eq. (17) performs similar to w_n in Eq. (16) does when the parameters are appropriately selected.

Obviously, there are still many other ways to construct mode vectors. Essentially, for a given application, the selection of appropriate v_n for the construction of jth order mode component can be formulated as an optimization problem. It is beyond the scope of the present work to discuss other alternative MoDEEs.

Similar to the IFD algorithm,^{42, 70} we extract the highest frequency component first, then residual is used as the initial value for u and v_n to extract the second highest frequency component. We continue this procedure until all components or intrinsic mode functions are separated. Note that this procedure is inherently nonlinear, even if the linear MoDEEs are used. Unlike the IFD algorithm, no (internal) iteration is used — to extract an intrinsic mode, we only solve the MoDEE once.

We also noted that one can extract the lowest intrinsic mode function first. Then the residual is used as the initial value in the MoDEE to extract second lowest intrinsic mode function. One continues this procedure until all the intrinsic modes are systematically separated one by one.

Using the PDE transform algorithm discussed above, intrinsic mode functions are systematically extracted from residues, a perfect reconstruction is straightforward by summing all the modes and final residue I_n, i.e.,

In practical computations, either the periodic boundary condition or the Neumann boundary condition can be used in all the proposed MoDEEs. When Fourier spectral methods is used for integration in this paper, periodic boundary condition is applied.

III.C Numerical methods for high-order MoDEEs

The proposed MoDEEs involve many high order PDEs which have to be solved in an appropriate manner to avoid the accumulation of numerical errors. Essentially, there are two important issues in the solution of high order evolution PDEs. The first issue concerns the accuracy of approximating high order derivatives. The rule of thumb here is that to approximate high order derivatives, high order numerical methods are required.^{72, 87} High order finite difference methods based on the Lagrange polynomials, Fourier pseudospectral methods and local spectral methods are suitable approaches for approximating high order derivatives. The second issue is the stability constraint in solving high order evolution PDEs. When the explicit (forward) Euler type of time discretization is used to solve an nth order evolution PDE, the step size Δt of time is constrained by Δt ~ (Δx)ⁿ, where Δt is the spatial grid spacing. Therefore, it is desirable to use alternative direction implicit (ADI) type of implicit methods for the time discretization. ADI type of methods has been previously developed by us to solve high order geometric flow equations.^{5, 18}

To solve the decoupled linear MoDEEs (15), we use the Fourier pseudospectral method. Spectral methods have been a popular choice for the numerical solution of various wave problems in recent years. As global methods, the Fourier spectral methods usually are much more efficient than local methods (e.g., finite difference and finite element methods) for certain classes of nonlinear problems. We applied Fourier pseudospectral method to the solution of the Navier-Stokes equation and other PDEs.^{65, 84} This approach is unconditionally stable and particularly efficient for solving the high oder MoDEEs proposed in Eq. (15). A windowed Fourier pseudospectral method has been developed for hyperbolic conservation laws, i.e., Euler equations.⁶⁵

IV.A Intrinsic mode decomposition

In this section, the MoDEE is applied to and validated on several benchmark testing cases to demonstrate the flexibility, efficiency, and accuracy of the method. In this paper, we use Fourier spectral methods to numerically integrate PDEs in one-step.

We first illustrate the efficiency and accuracy of the MoDEE algorithm III.B using signal f(x) = sin(x)+ sin(1.2x) + cos(2x) + sin(12x), as shown in Figure 1. In the figure, x-axis is in the domain of [0, 2π]. The total signal f(x) is shown in solid black curve, and the two modes with close frequencies, sin(1.2x) and sin(x), are chosen to be plotted with red dashed and blue dotted lines, respectively. The MoDEE decomposes the original signal into four modes which agree very well with the analytical results, and no visible difference can be observed when the 32nd order PDE is used in the MoDEE. Note that mode 3 of sin(1.2x) and mode 4 of sin(x) are closely embedded in f(x) and are thus difficult to be separated using many methods other than Fourier transform or wavelet methods. To investigate the MoDEE method in detail, various high order PDEs are used to separate the same signal. In Table 1, the MoDEE method illustrates its advances in separating out various modes clearly. When the highest order of the PDE used in the MoDEE is increased from 8 to 32, the L₂ error of various modes separated numerically decreases by orders of magnitude. The signal here is perfect candidate for the application of the Fourier method, and the MoDEE gives satisfactory decomposition accurately.

In Figures 2 and 3, plots are shown to compare the effect of high order PDEs on mode separation as illustrated as in Table 1. Mode 1 of sin(12x) is plotted in the Figure 2. Red circles in three panels show the result using the MoDEE method with different highest order of PDEs, and black curves show the analytical results of sin(12x) as reference. The highest order of PDEs as in ∇²ⁿ in Eq. (15) are 8, 16, and 32 respectively. For this mode, the use of up to the 16th order PDE is enough for the convergence. Similar plots of the decomposition of the third mode sin(1.2x) from the fourth one sin(x) are shown in Figure 3. A higher order PDE of up to 32nd order is required for the separation of closely overlapping modes 3 and 4. The results illustrate how important it is and when it is necessary to include very high order PDEs for signal processing. In addition, it is shown that the MoDEE using very high order PDEs is still very stable and does not diverge.

Lastly, we illustrate the spectral accuracy achieved by the MoDEE using signal composed of two modes of cos(mπx) and cos((m − 1) πx), where m is the maximum frequency supported by the mesh size (which is 0.05 in this example, and thus m = 20), i.e. two grid points for each oscillation period. In Figure 4, f(x) = cos(mπx) + cos((m − 1) πx) is decomposed into two frequency modes, cos(mπx) and cos((m − 1) πx), using the MoDEE. The comparison of numerical results are shown in Figures 4(b) and 4(c) respectively. Black curve shows the exact values, and red circles show the numerical results. Highly accurate MoDEE decomposition is observed in the figures where the 32nd order PDE is employed in the MoDEE algorithm as in Section III.A.

IV.B Image denoising, edge detection and enhancement

Image processing is becoming increasingly important in many areas of research in physical, mathematical and biological sciences.^{16, 29, 30, 40, 56} In particular, edge detection is a key issue in pattern recognition, computer vision, target tracking and image processing.^{2, 3, 9, 50, 60, 62} Closely related to edge detection is denoising. PDE methods provide powerful tools image processing. The MoDEE is designed for both accurately decomposing modes with different frequency (or “frequency-modes”) with spectral accuracy and effectively decomposing modes with different functions (of “functional-modes”) for image enhancement like edge detection, denoising, segmentation, etc.

The 512×512 grey-scale Barbara image as shown in Figure 5(a) is used to test the MoDEE algorithm. The same image has been extensively used in the literature due to its fine details associated with multiple edges. In Figure 5(b) random white Gaussian noise with signal to noise ratio (SNR) of 9.8 dB was added to the original image. For a quantitative view of the magnitude of SNR and effect of PDEs on the image enhancement, values of grey scale along a randomly chosen horizontal line is plotted before and after image processing. In Figure 5(a), such a horizontal line is specified by a dark grey line crossing the forehead of the Barbara image. In Figure 5(c), the black curve shows the grey scale values along the chosen horizontal line in the original Barbara image, and the red curve shows the values of the same horizontal line in the noise-added Barbara image 5(b). The added noises overlap closely with the edges of the original image, such that any denoising would smear the fine details of edges.

The MoDEE algorithm in III.A is applied to denoise the image while preserving the edge details as best as possible by taking difference between two PDEs with different diffusion constant and inclusion of high order PDEs. In Figure 6, values of grey scale along the horizontal line Figure 5(a) are plotted. Black solid curves show the original grey scale without noise. Blue curves show the grey scale of the same horizontal line after image processing using the MoDEE algorithm with evolution time Δt = 10, d_u2 = 0.05 and d_n2 = 0.4. Red curves show the grey scale using the same MoDEE algorithm with different value of d_u2 = 0. It is clear from the image, especially in those intervals as magnified and shown in Figures 6(b), 6(c) and 6(d), that the blue curves capture more details of various image edges with better accuracy than those red curves. The better accuracy illustrates the advantage of employing the difference between two high-pass filters in image processing as is employed in the MoDEE algorithms. For visual observation and comparison, denoised images using the two aforementioned MoDEE equations and parameters are shown in Figures 6(e) and 6(f). SNR is improved from 16.9 dB in Figure 6(e) to 17.4 dB in 6(e). Though the improvement in SNR alone is not very significant, the visual enhancement is more significant due to the improvement in capturing more fine edge details.

As given by Eq. (18), the MoDEE algorithm allows straightforward reconstruction of the original image through various modes and residue. In terms of image processing, the MoDEE algorithm automatically extracts noise, denoised image and image edges at the same time. The denoised image is shown in Figure 6(f), and the residual is therefore shown in Figure 7(a) which corresponds to the image edge. In addition, as aforementioned, the MoDEE is a general scheme allowing feature extraction and secondary processing of the decomposed functional modes. The residual in Figure 7(a) which contains image edges can be further post-processed by applying MoDEE algorithm once more. To achieve good performance, both the residue and edge mode from previous MoDEE application are utilized in the second MoDEE application. The reference-assisted “denoised” mode of the image in Figure 7(a) is shown in Figure 7(b). Comparing the two figures, it can be observed that edges in Figure 7(a) are finer and sharper, representing good edge extraction. On the other hand, edges in Figure 7(b) have varying density and look more natural to human eyes with more physical features captured in detail, e.g. the shade on the elbow of the left hand of Barbara. The edges in Figure 7(b) can be seen as an edge “representation” rather than edge extraction, where stronger and more meaningful edges are shown as denser lines while more subtle edges are shown as finer lines. In many applications this representation may prove to be very useful.

V.A Magnetic resonance images

MRI is built upon the physical theory of nuclear magnetic resonance (NMR): a nucleon such as proton possesses spin which has value of multiples of 1/2, and those nuclei with unpaired spins tend to align with the orientation of the externally applied magnetic field. After alignment, high frequency pulses, usually in the range of radio frequency (RF) of MHz, are emitted into a slice plane perpendicular to the external field to excite the aligned spins. RF waves are then turned off such that excited nucleus undergo a relaxation process to resume its original distribution. Different types of biological structures and tissues have different characteristic relaxation time which can be measured by MRI hardware and software to construct associated spatial images. Using MRI technique, each pulse sequence exploits some specific physical or chemical property of the protons of small, mobile molecules like water and lipids. Therefore, one could depict structural and functional information from living tissue at the sub-millimeter scale.

Though the medical imaging techniques have advanced tremendously in terms of spatial resolution, acquisition speed, and signal-to-noise ratio, medical images thus created still need to be carefully enhanced to account for factors like signal intensity inhomogeneities (i.e. bias fields), noise, and other artifacts. Noise, as one example, in MRI enters the data samples in k-space and competes with the NMR signal due to random fluctuations in the receiving coil electronics and in the patient body. In the case of function MRI (fMRI)⁸ for brain mapping and diffusion tensor imaging (DTI)^{73, 78} for studying neural fibers, the amount of random thermal noise entering MRI data in the time domain in the acquisitions limits the performance and usefulness of quantitative MRI diagnostics such as voxel-based tissue classification, extraction of organ shape or tissue boundaries, estimation of physiological parameters like tissue perfusion and contrast agent permeability from dynamic imaging.^{31, 44} Despite of much progress made in post-processing methods for noise reduction and image enhancement, it remains a challenge to find robust and interacting ways for noise removal, boundary preservation applicable to the different MRI acquisition techniques.

In Figure 8, the MoDEE algorithm using PDE (9) is used to detect edges for blur medical images obtained using MRI. Figure 8(a) is obtained using cardiac magnetic resonance imaging (cMRI) for a 21-year old woman referred for clinical management of idiopathic dilated cardiomyopathy (the image is authored by the European society of cardiology). By applying the MoDEE algorithm with evolution time Δt = 10, d_u2 = 0.2 and d_n2 = 0.52, one obtains satisfactory results with enhanced edge detection in Figure 8(b) which shows a dilated left ventricle, without evidence of tissue abnormalities (e.g. scars, patchy fibrosis).

V.B Magnetic resonance angiography images

Magnetic resonance angiography (MRA) image stands for a group of important techniques based on MRI to image blood vessels. MRA is often used to evaluate the arteries of the neck and brain, the thoracic and abdominal aorta, the renal arteries and the legs to diagnose symptoms like stenosis (abnormal narrowing) and occlusion or aneurysms (vessel wall dilatations, at risk of rupture), etc. An advantage of MRA compared to invasive catheter angiography is the non-invasive character of the examination. In particular, compared to computed tomography and angiography and catheter angiography, MRA does not expose patient to any ionizing radiation. The greatest drawbacks of the method are its comparatively high cost and its somewhat limited spatial resolution. In addition, the length of time the scans take can also be an issue, with computed tomography being far quicker. Therefore, image processing technique is important for post-processing the images generated by MRA to save cost, improve resolution and extract/highlight important features like fine details of edges of arteries.

The image shown in Figure 9(a) is downloaded from the website of Magnetic Resonance and Image Analysis Research Centre, University of Liverpool. The image is used to show the vessels known as the circle of will is in the brain. To test the ability of edge detection and noise removal, original image is added with noise. For a realistic setting the total noise added to the MRA image 9(a) is composed of two parts, η(x, y) = σ (x, y) + γ (x, y), where σ (x, y) is white Gaussian noise with maximum amplitude 25 and zero mean, and γ (x, y) is a periodically oscillating cosine noise given by

The noisy image is shown in Figure 9(b). To remove both noises closely embedded in the original image, we apply the MoDEE algorithm in Section III.A including 20th order PDE with d_u2 = 0.01 and d_n2 = 0.47. The resulting enhanced image is shown in Figure 9(c), in which both types of noises have been reduced and edges are clearer for medical diagnosis.

V.C X-ray computed tomography images

Computed tomography (CT) is a technology of using X-rays to image biological and human body. Compared to visible spectra imaging, X-rays have wavelengths between nanometer and picometer scales, which are energetic enough to penetrate biological tissues. Denser structures like the bones are more efficient at absorbing X-rays compared to softer tissues. This leads to the noninvasive medical imaging. CT is a popular technique for reconstructing the X-ray images with delicate hardware design as well as computer software. In CT scanning, a localized X-ray source and corresponding detector are put opposite to each other and rotate around the human body. The 2D slice image is then synthesized from the X-ray signals around all the directions. In the last few decades, CT technique has been widely used in medical imaging of soft tissues, hard bonds, blood vessels, etc.

In Figures 10(a) and 10(c), two different CT abdomen images are shown corresponding to different patient with slightly different shapes and edges of the organ in the abdomen. The MoDEE algorithm using PDE (9) is applied to detect edges in these two original CT images. The same set of MoDEE coefficients are used with d_u2 = 0.22 and d_n2 = 0.5. In particular, figure 10(a) shows the CT image of abdomen of a 68-year old male, image authored by Maclennan, Radiologist, Royal Alexandra Hospital, Paisley, United Kingdom. Edge enhanced image using the MoDEE algorithm is shown in Figure 10(b) which shows with better clarity (e.g. finer edges and tiny circles) the calcific density mass with flat upper margin in neck of gallbladder. The flat margin suggests that whatever is causing the mass is layering under the effect of gravity. Figure 10(c) shows a similar cross-sectional image of abdomen, with image authored by Department of Health and Human Services, using computerized axial tomography (CAT) scanning technique. The MoDEE captures the small differences which are reflected in the final edge-extracted image 10(d). The shape of the tiny circles in the liver are clearly highlighted.

As a short conclusion, through the applications and images shown in this section, it is illustrated that the MoDEE algorithm provides a robust and powerful tool for various image processing and enhancement tasks of interest to both clinical and research scientists equipped with various types of medical imaging techniques.