The aim of this paper is to present a new stable method for smoothing and differentiating noisy data defined on a bounded domain with . The proposed method stems from the smoothing properties of the classical diffusion equation; the smoothed data are obtained by solving a diffusion equation with the noisy data imposed as the initial condition. We analyze the stability and convergence of the proposed method and we give optimal convergence rates. One of the main advantages of this method lies in multivariable problems, where some of the other approaches are not easily generalized. Moreover, this approach does not require strong smoothness assumptions on the underlying data, which makes it appealing for detecting data corners or edges. Numerical examples demonstrate the feasibility and robustness of the method even with the presence of a large amounts of noise.
It is very common in many scientific applications, such as biological imaging, image processing and computer vision, inverse source and identification problems, to have to differentiate functions specified by data. In many cases, the data is experimental and so substantial noise or imprecision may be present. Since differentiation is an ill-posed problem , naive numerical differentiation techniques will amplify the noise resulting in inaccurate and mostly useless results. A concrete example is given by the function
where is the noisy version of the true data , and is a small number representing the noise level. The computed derivative is then
Evidently, the maximum error in the data is , which is small, while the maximum error in the computed derivative can be arbitrarily large, depending on the value of the frequency .
The last example showcases the instability associated with the differentiation process: small errors in the data might induce large errors in the computed derivative. To overcome the instability issue of the differentiation operator, several regularization techniques have been developed. In , a perturbation based method is used for single variable functions. In , a mollification method is given. Smoothing splines is adopted in , and results for single variable functions are given. Other methods such as the Lanczos’ generalized derivative , and the differentiation by integration using Jacobi polynomials , are known for single variable functions only. For a summary of other popular methods, see . Despite its importance in real life applications, only few methods dealing with higher dimensional cases have reported [8,9, 10,11]. A side from limited applicability due to dimension restriction, the aforementioned approaches are relatively slow to suit real life applications, particularly with big data.
In this paper, the smoothing property of the classical diffusion equation is utilized to reconstruct the gradient of a multivariable function from given noise-corrupted data satisfying
where is a bounded domain. For the time being, we shall assume that ; see the Appendix for the general case. More precisely, the proposed method amounts to take the function as stable approximation to where solves the diffusion equation:
subject to the initial and boundary conditions:
where is the Laplacian operator. The time can be considered as a design parameter, and its rule is to control the smoothness of the desired approximation; in the context of regularization theory, represents the regularization parameter. The precise formulation of the suggested method as well as the motivation behind it will be given in more detail in the sequel.
Under reasonable assumptions on true data , we prove stability results and we propose optimal convergence rate for the proposed method. Particularly, we show that the method gives an optimal rate of convergence of order . Several examples and applications in computer vision demonstrating the validity and efficiency of the proposed approach are also provided.
The rest of this paper is composed of five sections. Section 2 recalls some preliminary and auxiliary results. In Sections 3, we state the precise formulation of the proposed method and give the motivation behind it, and then we prove the main convergence results in Section 4. We provide several numerical experiments and applications in Section 5. More detailed discussion on the general case is provided in the Appendix.
First, we recall some definitions and results from the theory of Sobolev spaces. Readers may refer to [12,13] for more detailed discussions. Throughout the rest of the paper, we assume that is a bounded domain in , , with sufficiently smooth boundary .
We use the notation to denote the usual Lebesgue space of square-integrable functions, which is a Hilbert space with the inner product
The induced norm will be denoted by .
Let be a multi-index, and set . For a positive integer , the Sobolev space is defined by
where the partial derivatives are taken in the distributional sense. It is a Hilbert space under the inner product
and a Banach space under the the corresponding induced norm. In particular, the norm on is given by
The Hilbert space is defined as the closure of in , where is the set of all functions that are infinitely differentiable on and compactly supported in . We have the characterization
where the boundary value is understood in the trace sense.
By , we mean the space of all restrictions on from functions in . It is a Hilbert space with the norm
From the Sturm-Liouville theory , the eigenvalues of the eigenvalue problem
form a nondecreasing sequence of positive numbers tending to infinity, and the corresponding (normalized) eigenfunctions form an orthonormal basis for the space . Moreover,
Associated with the eigenvalue problem (3), the Hilbert space is defined as follows:
The parameter characterizes the regularity properties of the space . We have , , , and , see .
3 Motivation of the proposed algorithm
To motivate the proposed regularization approach, let be the solution of the parabolic initial-boundary value problem given by equations (1) and (2). Let us set for , and think of as a parameter. Then, on separating variables, we have the following representation
and so, due to the exponential term, the Fourier coefficients of decay faster than those of . It is also evident that for small values of the parameter . Thus, from this brief analysis, we infer that is an approximation of that is smoother than the input data .
As a consequence of the aforementioned discussion, it would be reasonable to take , , as smooth approximations to for . Therefore, we suggest the following numerical differentiation algorithm:
Step 1. Compute the smoothed version of the noisy data by
Step 2. Take the gradient as an approximation of .
Few remarks are in order. First, from expansion (5), we see that the parameter plays the role of a smoothing parameter, that is, sufficiently small or large values of result in undersmoothing or oversmoothing, respectively. Second, we point out that the gradient in the second step can be computed exactly or using an appropriate difference scheme depending on the nature of the available data.
4 Stability and consistency analysis
Our goal next is to examine computed by the proposed algorithm as a stable approximation to . In particular, we study the stability and consistency of the aforementioned approach and prove a convergence result. Throughout the sequel, we assume , and for convenience, we also let
We start off with the following stability result.
Set and . Since , we have . Then, it follows from (4) and the uniqueness of the Fourier coefficients that
By using the fact that and the Parseval’s identity, we see that
which concludes the proof.□
Throughout the sequel, we shall assume that , and we set if , and if . We have the following consistency result:
If , then
Set , and so . Moreover, from (4), the definition of , and the uniqueness of the Fourier coefficients, we have
The last equality together with the fact that imply
Now, it is an easy matter to show that
for all , from which we see
which ends the proof.□
Suppose that for some , then there exists a constant K independent of t and such that
If we choose the parameter so that for some , then we obtain the convergence result
as . The convergence rate is optimal when ; in this case, we have
We obtain the fastest convergence when , in which case, and
provided that for some positive constant .
5 Numerical experiments
Since data acquisition results in only discrete data in many scientific applications, we assume and that the noisy data, which we denote by , is sampled at regular grid points using the formula
where is some Gaussian noise with mean 0, and . Let us further assume that (refer to the attached appendix for general treatment). We will designate the symbol to denote the relative noise level in the data measured in the norm, that is,
where is the usual Euclidean norm. We shall assess the quality of the recovered gradient by the relative error given by
Regarding the choice of the parameter , we employ the Morozov’s discrepancy principle (e.g., ), which amounts to choose such that
For such domain , the eigenvalues and eigenfunctions of (3) are
Therefore, the Fourier discrete sine transform of the data is expressed as follows:
Using the trapezoidal rule over the regular grid , the smoothed data given by (5) can be computed as follows:
where . Thus, the smoothed data can be computed by
where the symbol ( ) denotes the elementwise multiplication operator, and the hat ( ) stands for the Fourier discrete sine transform.
In the following examples, the derivatives of the smoothed data are computed via the central difference formula:
and similarly for the true and noisy data.
Now, we examine examples in several space dimensions. We point out that the computations are extremely quick owing to the fast Fourier transform (FFT) routine, which is available in many computer algebra systems.
The true data are given by the function
|0.2||0.190 (0.809)||0.085 (9.084)||0.043 (89.13)|
|0.02||0.045 (0.080)||0.030 (0.908)||0.020 (8.913)|
|0.002||0.007 (0.008)||0.012 (0.090)||0.008 (0.891)|
The results in Example 5.1 suggest a convergence rate of . Since for all , from Remark 4.1, we deduce the theoretic optimal order . The deterioration in the convergence rate is mainly due to discretization errors and the nonoptimal choice of the smoothing parameter .
We consider the recovery of the gradient of the function
|0.2||0.138 (1.671)||0.056 (16.64)||0.047 (24.97)|
|0.02||0.051 (0.167)||0.023 (1.664)||0.021 (2.497)|
|0.002||0.013 (0.016)||0.011 (0.166)||0.009 (0.249)|
The true data are given by the function
|0.2||0.061 (4.705)||0.038 (9.174)||0.038 (13.64)|
|0.02||0.030 (0.470)||0.021 (0.917)||0.020 (1.364)|
|0.002||0.015 (0.047)||0.011 (0.091)||0.009 (0.136)|
Now we present an application on edge detection (emphasizing) of digital images.
Image gradient is an essential ingredient in common image processing and computer vision applications, particularly in edge detection algorithms. Since the gradient is sensitive to noise, noisy images must be denoised or smoothed prior to the gradient estimation. In common edge detection algorithms, such as the Canny method, a Gaussian filter is typically used for this purpose. In the following demonstration, we employ our method to smooth out the noisy images, and then the image gradient is approximated using the central difference formula. Consequently, the magnitude of the image gradient at pixel is approximated by
where is the underlying image intensity function.
In this demonstration, we consider three noisy images that are obtained by adding a Poisson noise to the noise free images. The original, noisy, and smoothed images by the diffusion method are shown in Figure 4. To assess the quality of the smoothed (denoised) images, we computed the peak signal-to-noise ratio (PSNR) of the noisy and the smoothed images, given by
where and are the original and denoised images, respectively. The PSNR results are summarized in Table 4.
The gradient magnitudes of the noisy and smoothed images are plotted in Figure 5. For the sake of comparison, we also plot the edges of the noisy images obtained by the popular Canny method , which consists of five main steps: Gaussian filtering, measuring intensity gradient, edge thinning, double thresholding to determine strong and weak edges, and edge tracking via blob analysis to extract weak edges.
From Table 4, we see that the PSNR of the denoised images are significantly higher than those of noisy images, which indicates the denoised images by the proposed smoothing method have better qualities than the noisy images. Moreover, as it can be observed from Figures 4 and 5, the proposed technique can significantly enhance the quality of the recovered edges, and it seems that the smoothing stage has no bias toward certain direction; this is a major advantage since edge orientation is mostly unknown in real-life applications. It is evident from this experiment that diffusion smoothing gives plausible results, and it can be well suited for such type of applications.
|Image||Noisy image||Denoised image|
From the presented examples, we see that the proposed diffusion technique is very suitable for data smoothing and numerical differentiation, and it can excel in various scientific applications due to its speed (in particular, when utilizing the fast Fourier transform), supporting theory behind it, and further, it does not pose any dimension limitations; the method can be extended easily to problems in any dimension, contrary to the most existing methods in which they are limited to one-dimensional domains.
The smoothing property of the diffusion equation is utilized to estimate the gradient of functions specified by noisy data. We proved stability and consistency results and carried out the convergence analysis for the proposed method. Several examples in one and higher dimensional domains demonstrate the feasibility and robustness of the proposed approach.
The main advantage of the method lies in multivariable problems, where some of the other smoothing approaches do not easily generalize. Moreover, the method shows decent results under mild assumptions on the true data even in the presence of a large amount of noise, while these could be limitations for other methods. We believe the method can excel in various scientific applications due to its speed, relative ease of implementation, and the supporting theory behind it; we look forward to extending this approach and the theoretical results for higher derivatives estimations. Furthermore, appropriate a priori and a posteriori parameter choice strategies for choosing are under investigation.
The authors would like to thank the Scientific Research and Graduate Studies at Yarmouk University for supporting this paper, and they would like express their sincere gratitude to the anonymous referees for their helpful comments that helped us to improve the quality of the manuscript.
Conflict of interest: The authors state no conflict of interest.
If the true boundary data is known approximately by the function satisfying
then smoothed data can be defined as follows:
where the function satisfies (in the weak sense) the boundary-value problems:
To derive stability and convergence results, we set
where solves the boundary-value problems
From the definition of (e.g., ), we see that
where depends on and , provided . The last inequalities together with triangle inequality imply the main error result:
Moreover, if we let , and choose for some , then we obtain the convergence result
as . The convergence rate is optimal when , and in this case, we have
We obtain the fastest convergence when for some , in which case and
provided that for some positive constant .
 D. A. Murio, C. E. Mejia, and S. Zhan, Discrete mollification and automatic numerical differentiation, Comput. Math. Appl. 35 (1998), no. 5, 1–16, DOI: https://doi.org/10.1016/S0898-1221(98)00001-7. 10.1016/S0898-1221(98)00001-7Search in Google Scholar
 D. Y. Liu, O. Gibaru, and W. Perruquetti, Differentiation by integration with Jacobi polynomials, J. Comput. Appl. Math. 235 (2011), no. 9, 3015–3032, https://doi.org/10.1016/j.cam.2010.12.023. Search in Google Scholar
 I. Knowles and R. J. Renka, Methods for numerical differentiation of noisy data, Electron. J. Differ. Equ. Conf. 21 (2014), 235–246. Search in Google Scholar
 M. F. Al-Jamal, A. K. Alomari, and M. S. Gockenbach, Smoothing via elliptic operators with application to edge detection, Inverse Probl. Sci. Eng. 26 (2018), no. 5, 657–676, https://doi.org/10.1080/17415977.2017.1336552. Search in Google Scholar
 Z. Meng, Z. Zhao, D. Mei, and Y. Zhou, Numerical differentiation for two-dimensional functions by a Fourier extension method, Inverse Probl. Sci. Eng. 28 (2020), no. 1, 126–143, https://doi.org/10.1080/17415977.2019.1661410. Search in Google Scholar
 A. A. Yahya, J. Tan, B. Su, K. Liu, and A. N. Hadi, Image edge detection method based on anisotropic diffusion and total variation models, J. Eng. 2019 (2019), no. 2, 455–460, https://doi.org/10.1049/joe.2018.5345. Search in Google Scholar
 C. Brito-Loeza, R. Legarda-Saenz, and A. Martin-Gonzalez, A fast algorithm for a total variation based phase demodulation model, Numer. Methods Partial Differential Equations 36 (2020), no. 3, 617–636, DOI: https://doi.org/10.1002/num.22444. 10.1002/num.22444Search in Google Scholar
 R. C. McOwen, Partial Differential Equations: Methods and Applications, Prentice Hall, New Jersy, 2003. Search in Google Scholar
© 2022 Mohammad F. Al-Jamal et al., published by De Gruyter
This work is licensed under the Creative Commons Attribution 4.0 International License.