Abstract
The main purpose of image enhancement technology is to improve the quality of the image to better assist those activities of daily life that are widely dependent on it like healthcare, industries, education, and surveillance. Due to the influence of complex environments, there are risks of insufficient detail and low contrast in some images. Existing enhancement algorithms are prone to overexposure and improper detail processing. This paper attempts to improve the treatment effect of Phase Stretch Transform (PST) on the information of low and medium frequencies. For this purpose, an image enhancement algorithm on the basis of fractional-order PST and relative total variation (FOPSTRTV) is developed to address the task. In this algorithm, the noise in the original image is removed by low-pass filtering, the edges of images are extracted by fractional-order PST, and then the images are fused with extracted edges through RTV. Finally, extensive experiments were used to verify the effect of the proposed algorithm with different datasets.
1. Introduction
Image sharpening means enhancing, that is, highlighting edges and texture features, while suppressing those unimportant planner areas, planner areas, and parts with constant grey levels. In addition, the original image is fused with enhanced edge and texture images (the simplest fusion is a pixel-by-pixel weighted average) to ensure a smooth transition between its grey levels. Therefore, in the enhanced image, the edges and key textures are highlighted, and other areas are also sharpened.
In the process of image enhancement, the causes of image degradation have not been analyzed, and the processed image is not necessarily close to the original image. In the literature [1], the K-means clustering algorithm was used to divide the image into several grey-scale intervals and then equate them separately. This algorithm has a better-enhanced effect on the images for grayscale distribution at both ends and more pixels in the grey-scale regions. Literature [2] proposed an algorithm to enhance the adaptive image based on a double histogram equation, which has to keep the average brightness of the output image close to the original image and avoid the tendency to magnify. Literature [3] proposed an image enhancement algorithm based on wavelet analysis and Retinex algorithm. Wavelet transform decomposes the image into multiscale images, removes noise from images with different frequencies, and uses a Retinex algorithm to enhance image details. Therefore, the combination of the two methods can improve the overall visual effect of the image and better highlight the details of the image. In the literature [4, 5], a better enhancement effect was achieved by enhancing the image contrast and the singular matrix of the image defined in the wavelet domain. In [6], an improved algorithm is proposed, which combines the adaptive total variation model and the impact filter model. Literature [7] used the minimum error method to automatically calculate the segmentation point of the piecewise linear grey-scale transformation, which improved the visual effect of the image. In [8], the improved histogram equalization method and the contrast enhancement algorithm were combined to improve the local information contrast enhancement of low illumination images. An adaptive smoothing and enhancement method for images was proposed in [9]; when filtering, the filter weights of the processed pixels are dynamically determined; after smoothing the inside of the image region, the edge of the region in the image is also sharpened and enhanced, which effectively solves the contradiction between image smoothing and enhancement.
The above literature all involves the image enhancement algorithms based on time domain and frequency domain variation. But it has been rarely found in domestic and foreign literature about the image enhancement for phase information changes. In 2015, M.H. Asghari and B. Jalali proposed a digital image transformation inspired by physical phenomena, called Phase Stretch Transform (PST). This paper attempts to improve the treatment effect of the PST on the information of low and medium frequency.
2. Relate Work
2.1. Basic PST
In 2015, M.H. Asghari and B. Jalali proposed a digital image transformation inspired by physical phenomena, called Phase Stretch Transform, which can simulate the propagation process of electromagnetic waves in the diffraction medium with the warping dispersive dielectric function. This method using phase-detection edges can simulate the diffraction process with an all-pass phase filter having a specific frequency and divergence dependency. The output phase prototype of the filter shows the change in image intensity value in the spatial domain, and the edge detection can be achieved after the thresholding and morphological processing of the phase prototype.
Figure 1 shows the process of digital image edge detection based on PST. A local filter kernel function first smoothed the original image, and then the phase operation of the nonlinear frequency function was performed in the frequency domain, called PST. Finally, edge detection was achieved through postprocessing such as thresholding and morphological filtering.
The mathematical model of PST in the frequency domain iswhere A(m, n) is the angle image, “∠” is an angle operation, B(m, n) is an original input image, FFT2 and IFFT2 are a two-dimensional fast Fourier transform and an inverse transform, respectively, and (u, v) is the frequency variable. is the frequency response of the local smoothing low-pass filter, is a frequency-dependent nonlinear phase warping kernel function, and is a nonlinear function of the frequency variable.
Although any phase kernel function can be considered in the PST, the derivative of the kernel phase function is a linear or sublinear function of the frequency variable according to the results of the literature [10]. For the phase kernel function prototype, a simple example is the inverse-tangent function of the “S” type. For the sake of simplicity, if it is further required that the phase warping operation is isotropic in the plain of the frequency domain, the degree of warping is only related to the polar radius r in the planar polar coordinate system of the o-uv frequency, but independent of the polar angle θ, that is, assuming that the nuclear phase prototype of the PST is circularly symmetric about the frequency variable, the PST core phase shall be obtained:where r is the polar radius in the planar polar o-uv polar coordinate system of frequency domain; θ is the polar angle, and its relationship with the uv frequency variable is , . If it is required that the derivative of φ polar(r) with respect to r is an inverse tangent of the S-type, then
φ polar(r): note that the uv frequency plane after the Fourier transform of the image is a finite region, then φpolar(r) can be solved according to
The phase function in equation (4) was normalized to obtain :
For the phase function in equation (5), the phase tensile strength parameter S and the warping parameter W in the nonlinear distortion stretching transformation were added to obtain the final kernel phase function with strength parameter S and warping parameter in the PST transformation:where tan-1(•) represents the inverse-tangent function, ln(•) is the natural logarithm, and rmax represents the maximum frequency polar radium of the uv frequency plane.
A PST having a kernel phase shaped as shown in equation (6) is applied to the spectrum of the digital image, forming a phase image A(m, n). For the application of the edge detection, subsequence operations can be implemented such as translation, thresholding, and morphological processing.
2.2. Fractional Order
2.2.1. Basic Theory of Fractional Order
The Fractional Fourier Transform (FrFT) has more new unique features based on retaining the original properties of the traditional Fourier transform, so it can be considered to be a generalized Fourier transform [11].
In general, the -order FRFT of the function can be expressed as or as needed. can be seen as an operator acting on the signal .
The basic definition of the FRFT was given below from the perspective of linear integral transformation. The -order FRFT of the function in the time domain is defined as a linear integral operation:
Among them, and is called as the kernel function of the FRFT, where , , , and is an integer.
By variable substitution and , equation (7) can be further expressed as
In equations (1)–(8), ; the linearity of FRFT given cannot indicate that this equation is unchanged, except for , because the kernel function is not only the function of , but also that of order . At = 1, , , and
Thus, is the common Fourier transform of . Similarly, is the inverse transform of in the traditional Fourier transform. Because can only appear in the parameter position of the trigonometric function, the parameter is defined by a period of 4, and it only needs to examine the interval . At , ; at , .
All the above can be expressed as operators:where are the arbitrary integers.
The additivity of fractional order is a very important property of the FRFT, which can be expressed as
By using Gauss integrals to give direct integral representations, the operation shall be simpler, and then
In summary, the FRFT can be explained as follows: only considering the interval , the FRFT is the original function at , and it is the ordinary Fourier transform at ; when gradually changes from 0 to 1, its FRFT smoothly changes from the original function to the common Fourier transform.
The FRFT can also be interpreted in another way. That is, the FRFT is defined as a rotation transform of a time-frequency plane, and the FRFT of -order is a linear canonical transformation defined by a transformation matrix. Transform matrix is given as
Among them, . According to the definition of Radon transform, in a plane along different lines (the distance between the line and the origin is , and the direction angle is ), the line integral is made for f(x, y), and the obtained is the Radon transform of the function . Then, the matrix can be regarded as the two-dimensional rotation matrix on the time-frequency plane.
As shown in Figure 2, the Fourier transform can be considered as the representation of the function rotating the angle of from the axis to the ω axis on the time-frequency plane; that is, the original function is mapped from the time domain to the frequency domain of the angle by Fourier transform; based on this, the FRFT is performed for the and -order on the function, that is, with the fractional-order operator and , the function is rotated by the angle () and the angle (), respectively, and then mapped to and orders.
2.2.2. Discretization Method of Fractional Fourier Transform
The fast Fourier transform greatly promotes the development of the Fourier transform. Similarly, the fast FRFT algorithm will also rapidly develop the FRFT in the field of signal processing. Therefore, it is particularly important to study the discrete FRFT fast algorithm [12].
Ozaktas adopted the FRFT calculation method. This method performs N-point sampling on the time domain of the original function and also maps it to N sample points in the fractional Fourier domain. The computational complexity of the algorithm is O(N log N). Before using this method to calculate the FRFT, the original signal must be dimensionally normalized. Afterwards, in both the time and frequency domains, the representation of the signals is dimensionless and the length of support is equal to . This also indicates that the Wigner distribution of the signal is limited to the unit circle with as the radius and the origin of the time-frequency plane as the centre. In order to obtain an efficient calculation method, the calculation of the FRFT is decomposed into a convolutional form. According to the definition of FRFT, the -order FRFT of the signal can be written as
It can be seen from equations (1)–(13) that the calculation of the FRFT can be decomposed into three steps. Step 1: the signal x(t) was multiplied by a chirp function, to obtain the intermediate result, which is recorded as . Then, the frequency domain bandwidth of becomes 2 times more than that of the signal x(t), so the sampling interval of should be . The sampling interval of the original signal x(t) is . At this time, if the sampling interval of x(t) is to be changed into , the signal x(t) needs to be interpolated twice to obtain the signal x(t) at the sampling interval and then multiply by the intermediate signal at the sampling interval . Step 2: the signal is convolved with a chirp signal, because the bandwidth of is , so according to the convolution theorem, the chirp signal can be represented by its band-limited form , denoted as h(t): where is the Fourier transform of the convolved chirp signal. The convolution is written in discrete forms: This convolution formula can be calculated using the Fast Fourier Transform algorithm (FFT). Step 3: multiply another chirp signal, to obtain the 2N sample points of at the sampling interval of . Because it is a mapping of N sampling points in the time domain to N sampling points in the fractional Fourier domain, the sampling at the sampling interval can be obtained by performing a double extraction.
Let and x be the column vectors consisting of N sample points of and x(t), respectively; then, the above calculation process can be written in matrix form:where and represent the extraction and interpolation operations and and are the corresponding chirp multiplication and chirp convolution operations.
The method above expresses the FRFT as a convolution operation [13]. It is also possible to represent the FRFT in another form:
Sampling equations (2)–(13), it is calculated as
The summation of equation (18) can be expressed as a convolutional form of the signal, which was calculated using FFT. Finally, a 2x extraction was performed to obtain the sampling at a sampling interval of . Similarly, the above sampling process can be expressed in matrix form as
3. Theoretical Analysis
3.1. Fractional-Order-Based PST
3.1.1. PST Improvement Based on Fractional Order
It has been proven that the PST phase helps in detecting edges in the image as well as dramatic changes in intensity values [14]. However, edge detection for low image variations is less than ideal.
In image enhancement, the fractional differential can enhance the high-frequency component while enhancing the mid-low frequency component. To a certain extent, it can also nonlinearly preserve the DC component of the image, making the image texture details clearer and overcoming the defect that the integer-order differential can weaken the low-frequency information. With the efforts of scholars, the fractional differential enhancement has achieved certain results. Yang Zhuzhong et al. [15] used the G-L fractional differential to construct the Tiansi differential operator for image enhancement and stated that the differential order between 0.4 and 0.6 has a better enhancement effect; when the order is too large, the noise is also enhanced. In order to pass more low-frequency information, Wu Ruifang et al. [16] improved the Tiansi operator template. Wang Weixing et al. [17] modified the template to 8 different directions for greatly enhancing the edge information. Zhang Yu et al. [18] proposed an adaptive fractional-order enhancement algorithm, which applies an exponential function to determine the fractional order. Chen Qingli et al. [19] put forward the fractional-order enhancement algorithm based on Caputo's definition, which lays a good foundation for our study.
To this end, the paper combines PST and fractional order to improve the stretching operator of PST. The improved mathematical model is given aswhere A(m, n) is the angle image, “∠” is the angle operation, B(m, n) is the original input image, FRFT(∗) and FRFT(-∗), respectively, represent the fractional Fourier transform and Inverse transform, is the frequency variable, is the frequency response of the local smoothing filter, is a frequency-dependent nonlinear phase warping kernel function, and is a nonlinear function of the frequency variable.
3.2. Comparison and Analysis
It is assumed that the angle of the FRFT is , where X is the vector in the X-axis direction and Y is the vector in the Y-axis direction; the tensile strength parameter of PST is S and the warping parameter is W.
Then, the edge detection test results under different values of X, y, s, and W are compared, as shown in Figure 3.
(a)
(b)
(c)
(d)
3.3. Relative Total Variation (RTV) Analysis
To achieve the image enhancement effect, the extracted image from the edge and the original image need to be superimposed. The superimposed image highlights the edge features and achieves the enhanced result, but it may have noise, edge burrs, and so on. For this reason, a relative total variation (RTV) should be processed on the image to fuse its main and background images. The main principles of RTV are as follows.
To further enhance the contrast between texture and structure, especially for areas that are visually prominent, and were combined together to form a more efficient structural texture decomposition.
The objective function is assumed as
Among them, ensures that there is no significant deviation between the input and the result. The introduction of the new regularization term can achieve the effect of removing image texture, which is called relative total variation (RTV). λ in formula (22) is a weight; ε is a small positive number, avoiding being divided by zero.
The objective function in equation (22) is nonconvex. Therefore, its solution cannot be obtained in a common manner. Based on the quadratic measure penalty, the objective function was proposed, which is an effective solution method to make linear optimization (Szeliski, 2006; Lischinski et al., 2006; and Krishnan and Szeliski, 2011).
Our approach is to decompose RTV metrics into nonlinear terms and quadratic terms. Interestingly, the problem of the nonlinear part can be transformed into a series of linear equations by means of an iterative reweighted least squares method.
The metric in the direction was firstly discussed and then that in direction. The penalty term is expanded as
By regrouping the factor and the set element , it is derived as
The second line is an approximation because for numerical stability was introduced. Through factor rearrangement, the metric was decomposed into a quadratic term and a nonlinear part , which are expressed as
Equation (25) indicates that of each pixel actually merges adjacent gradient information in an isotropic spatial filtering manner; is a Gaussian filter with a standard deviation of . The division operation in (25) is performed by elements; ∗ is the convolution operator; and is only related to the pixel gradient.
Similarly, the direction penalty term can be expressed as
Among them, the secondary component partial derivatives and are similar to the nonlinear part. They are respectively expressed as
Through these operations, equation (22) can be written in matrix form:where and are vector representations of and , respectively; and are Toeplitz matrices of the forward differential discrete gradient operators; and , , , and are diagonal matrices. Their diagonal values are
This makes it possible to achieve a special iterative optimization process. Due to the decomposition of the nonlinear part and the quadratic part, a numerically stable approximation is naturally obtained, which was found to be very effective in estimating fast structural and texture images in experiments.
4. Comparison of Test Results
In this article, the test is divided into two parts. In the first part, the test data is based on tensile strength parameter S and the warping parameter W of the fractional PST takes different values. Let the angle of the FRFT be , where X is the vector in the X-axis direction, and the vector in which Y is in the Y-axis direction takes different values; the smoothness parameter Lam and the texture element parameter Sig of the RTV take different values. Taking multiple sets of images as test objects, the specific test results obtained are shown in Figures 4–7 and Table 1.
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
For the above test results, it is found that the tensile strength parameter S and the warping parameter W of the fractional PST are close to 0.5 and 10, respectively, and the angle of the FRFT is ; when the values of X and Y are closer to 1, more edge information can be obtained, but more noise is added; when the smoothness parameter Lam and the texture size parameter Sig of the RTV are larger, the smoothing effect will be better, but the important edge information of the image will be blurred.
In the second part of the test, a comparative analysis was performed on image enhancement based on the fractional-order PST, phase consistency, region growing, histogram equalization, Canny operator, and so forth, respectively. The specific test results are shown in Figures 8 and 9.
(a)
(b)
(c)
(d)
The results of the above test show that some image enhancement algorithms such as the region growing-based image enhancement algorithm do not highlight the edge information while enhancing the high-frequency information of the image, and the resulting image will mask some important information; the image enhancement algorithms based on the canny operator and so forth may blur small edge information while presenting the edge portions of the image; the phase consistency-based edge enhancement algorithm can add edge burr or noise while highlighting edge information.
5. Conclusions
This article proposed a novel image enhancement algorithm, called FOPSTRTV to improve the quality of the image. Using the proposed algorithm to process the image, the noise in the original image is first eliminated by low-pass filtering. Then, the edge of the image is extracted through the differential order PST because the differential order PST processing can help to get the phase information and extract the image edge information well; in contrast, the edge extraction of the ordinary image is based on the time domain and frequency domain information. Finally, RTV can be used to fuse the image, giving a good visual effect. Extensive experiments with several representative algorithms demonstrated that FOPSTRTV is competitive and promising.
Data Availability
The data are available online at https://github.com/sunwww168/PSTRTV.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.