An Image Enhancement Algorithm Based on Fractional-Order Phase Stretch Transform and Relative Total Variation

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.


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. erefore, 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.
is 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. erefore, 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. e 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). is paper attempts to improve the treatment effect of the PST on the information of low and medium frequency. is method using phase-detection edges can simulate the diffraction process with an all-pass phase filter having a specific frequency and divergence dependency. e 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. e mathematical model of PST in the frequency domain is

Relate Work
where 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. L(u, v) is the frequency response of the local smoothing lowpass filter, K(u, v) � e jϕ(u,v) is a frequency-dependent nonlinear phase warping kernel function, and ϕ(u, v) 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 ϕ(u, v) 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 inversetangent 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 r � 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 e phase function in equation (4) was normalized to obtain ϕ N : 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 ϕ N (r, W, S) 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.  [11].

Fractional
In general, the p-order FRFT of the function x(t) can be expressed as X p (u) or F p x(t) as needed. F p x(t) can be seen as an operator F p acting on the signal x(t).
e basic definition of the FRFT was given below from the perspective of linear integral transformation. e p-order FRFT of the function x(t) in the time domain is defined as a linear integral operation: Among them, (7) can be further expressed as ; the linearity of FRFT given cannot indicate that this equation is unchanged, except for (p � 4n), because the kernel function is not only the function of (u, t), but also that of order p. At p � 1, α � π/2, A α � 1, and us, X 1 (u) is the common Fourier transform of x(t). Similarly, X − 1 (u) is the inverse transform of x(t) in the traditional Fourier transform. Because α � pπ/2 can only appear in the parameter position of the trigonometric function, the parameter p is defined by a period of 4, and it only needs to examine the interval All the above can be expressed as operators: where n, n ′ are the arbitrary integers.
e 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 0 ≤ p ≤ 1, the FRFT is the original function at p � 0, and it is the ordinary Fourier transform at p � 1; when p gradually changes from 0 to 1, its FRFT smoothly changes from the original function to the common Fourier transform. e FRFT can also be interpreted in another way. at is, the FRFT is defined as a rotation transform of a time-frequency plane, and the FRFT of p-order is a linear canonical transformation defined by a transformation matrix. Transform matrix is given as Among them, α � pπ/2. According to the definition of Radon transform, in a plane along different lines (the distance between the line and the origin is d, and the direction Computational Intelligence and Neuroscience angle is α), the line integral is made for f(x, y), and the obtained F(d, α) is the Radon transform of the function f. en, 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 x(t) rotating the angle of π/2 from the t 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 π/2 by Fourier transform; based on this, the FRFT is performed for the p 1 and p 2 -order on the function, that is, with the fractional-order operator F p 1 and F p 2 , the function is rotated by the α 1 angle (α 1 � p 1 π/2) and the α 2 angle (α 2 � p 2 π/2), respectively, and then mapped to p 1 and p 2 orders.

Discretization Method of Fractional Fourier Transform.
e 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. erefore, it is particularly important to study the discrete FRFT fast algorithm [12].
Ozaktas adopted the FRFT calculation method. is 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. e 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 Δx. is also indicates that the Wigner distribution of the signal is limited to the unit circle with Δx/2 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 a-order FRFTof the signal x(t) 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 g(t). en, the frequency domain bandwidth of g(t) becomes 2 times more than that of the signal x(t), so the sampling interval of g(t) should be 1/2Δx. e sampling interval of the original signal x(t) is 1/Δx. At this time, if the sampling interval of x(t) is to be changed into 1/2Δx, the signal x(t) needs to be interpolated twice to obtain the signal x(t) at the sampling interval 1/2Δx and then multiply by the intermediate signal g(t) at the sampling interval 1/2Δx.
Step 2: the signal g(t) is convolved with a chirp signal, because the bandwidth of g(t) is 2Δx, so according to the convolution theorem, the chirp signal can be represented by its band-limited form 2Δx, denoted as h(t): where H(Ω) is the Fourier transform of the convolved chirp signal. e convolution is written in discrete forms: is 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 X α (u) at the sampling interval of 1/2Δx. Because it is a mapping of N sampling points in the time domain to N sampling points in the fractional Fourier domain, the X α (u)sampling at the sampling interval 1/Δx can be obtained by performing a double extraction.
Let X α and x be the column vectors consisting of N sample points of X α (u)and x(t), respectively; then, the above calculation process can be written in matrix form: where D and J represent the extraction and interpolation operations and Λ and H are the corresponding chirp multiplication and chirp convolution operations. e method above expresses the FRFT as a convolution operation [13]. It is also possible to represent the FRFT in another form: · e j1/2(u− t) 2 csc α dt.  · e j1/2(m− n/2Δx) 2 csc α . (19) e 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 X α (u) sampling at a sampling interval of 1/Δx. Similarly, the above sampling process can be expressed in matrix form as

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 lowfrequency 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. e improved mathematical model is given as where 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, (u, v) is the frequency variable, L(u, v) is the frequency response of the local smoothing filter, K(u, v) � e jϕ(u,v) is a frequency-dependent nonlinear phase warping kernel function, and ϕ(u, v) is a nonlinear function of the frequency variable.

Comparison and Analysis.
It is assumed that the angle of the FRFT is α � X/Y, 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.
en, the edge detection test results under different values of X, y, s, and W are compared, as shown in Figure 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. e 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. e main principles of RTV are as follows.
To further enhance the contrast between texture and structure, especially for areas that are visually prominent, L and H were combined together to form a more efficient structural texture decomposition. e objective function is assumed as arg min Among them, (S p − I p ) 2 ensures that there is no significant deviation between the input and the result. e introduction of the new regularization term ( H x (p)/L x (p) + ε + H y (p)/L y (p) + ε ) 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.
Computational Intelligence and Neuroscience e objective function in equation (22) is nonconvex. erefore, 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, 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. e metric in the x direction was firstly discussed and then that in y direction. e penalty term is expanded as By regrouping the factor and the set element |(z x S) q |, it is derived as (24) e second line is an approximation because ε s for numerical stability was introduced.
rough factor rearrangement, the metric was decomposed into a quadratic term (z x S) 2 q and a nonlinear part u xq w xq , which are expressed as Equation (25) indicates that u x of each pixel actually merges adjacent gradient information in an isotropic spatial filtering manner; G σ is a Gaussian filter with a standard deviation of σ. e division operation in (25) is performed by elements; * is the convolution operator; and w x is only related to the pixel gradient.
Similarly, the ydirection penalty term can be expressed as Among them, the secondary y component partial derivatives (z y S) 2 q and u yq w yq are similar to the nonlinear part. ey are respectively expressed as rough these operations, equation (22) can be written in matrix form: where v s and v I are vector representations of S and I, respectively; C x and C y are Toeplitz matrices of the forward differential discrete gradient operators; and U x , U y , W x , and W y are diagonal matrices. eir diagonal values are is 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.

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 α � X/Y, 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.
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 α � X/Y; 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 fractionalorder PST, phase consistency, region growing, histogram equalization, Canny operator, and so forth, respectively. e specific test results are shown in Figures 8 and 9.
e 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    Computational Intelligence and Neuroscience 9 algorithm can add edge burr or noise while highlighting edge information.

Conclusions
is 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. en, 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.