Super-Resolution Reconstruction of Depth Image Based on Kriging Interpolation

Abstract

The super-resolution of depth images is a research hotspot. In this study, the classical Kriging algorithm is applied to the spatial interpolation of depth images, together with the fractional-order differential method for edge recognition, to realise the super-resolution reconstruction of depth images. The resulting interpolation model improves the edge performance of Kriging interpolation by harnessing the superior characteristics of fractional-order differential edge recognition and effectively solving the edge blurring problem in super-resolution interpolation of depth images. Experimental results show that, compared with the classical algorithms, the super-resolution reconstruction based on Kriging interpolation is greatly improved in terms of visual effects and the peak signal-to-noise ratio of the depth image. In particular, edge recognition based on fractional-order differentiation solves the image blurring problem at the edges of the depth images. Inspection of the point clouds of the depth images shows that the output of the proposed interpolation model has obvious fractal characteristics.

1. Introduction

With the development of computer vision and image processing, deep image processing has become a hot research topic in recent years. Deep image processing has both theoretical significance and important practical applications. For example, in the field of virtual reality video games, where the player’s pose movements are acquired through somatosensory devices and subjected to human–computer interaction and modelling, depth images acquired in real time must be processed to improve the quality of the models constructed from depth data. Another topic of importance is the super-resolution reconstruction of depth images, which is critical for applications in vision fields such as virtual reality video games, biomedicine, and augmented reality, as these applications require highly accurate depth image data. High-precision super-resolution reconstruction helps generate highly realistic three-dimensional (3D) object surface models. In addition, artificial intelligence and autonomous driving both require highly accurate depth image data. Developments in computer vision and image processing promise to make the exploration of depth-image-based super-resolution algorithms a research hotspot.

Recently, image super-resolution using deep learning methods has achieved remarkable success in producing high-quality and visually pleasing results [1,2,3]. However, these methods usually require a large amount of training data and computational resources, making them impractical for real-time applications. To overcome this limitation, some researchers have explored the use of non-deep learning-based techniques, such as interpolation-based methods [4,5]. While these methods are computationally efficient, they may not be suitable for processing complex depth images.

Recently, a few studies have proposed the use of fractional-order differentiation in interpolation-based methods for depth image super-resolution [6,7]. These methods have been shown to be effective in preserving texture information and edge features as well as improving noise reduction performance. In comparison with traditional interpolation-based methods, these methods offer better image reconstruction and greater robustness.

Depth image super-resolution is a signal processing technique for obtaining high-resolution images from single or multiple frames of low-resolution images [8] and has promising practical applications [9]. Image super-resolution techniques can be broadly classified into four categories: interpolation-based [10], reconstruction-based [11], enhancement-based [12], and learning-based [1,2] methods. Fractional-order calculus has been widely applied in cybernetics, probability theory, and signal processing [13]. One study [14] proposed a numerical implementation of fractional-order calculus in the field of signal analysis and its initial application to two-dimensional (2D) image processing; the same study also proposed the Tiansi mask operator based on fractional-order calculus, which achieved good results in the field of image enhancement and edge detection. Other studies [15,16,17] have proposed image denoising algorithms based on a small probability model and gradient-adaptive fractional-order calculus, achieving good results in image denoising and enhancement.

Spatial interpolation is a technique used to obtain attribute information at the unsampled points (interpolation points) by studying the attribute information of the sampled points. Using a mathematical model, the mapping relationship between the attribute domain and the spatial domain is analysed on the basis of known information from a certain number of sampled points [18]. Spatial interpolation techniques are widely used in tasks such as precipitation prediction, soil level classification, and groundwater simulation [19]. Numerous interpolation methods are in common use, such as inverse distance weighted interpolation [20], Kriging interpolation [21,22], minimum curvature interpolation, the Sheppard method, natural neighbour interpolation [20], nearest-neighbour interpolation, multiple regression interpolation, the radial basis function method, the linear interpolation triangular network method, and the moving average method.

The classical super-resolution algorithm for images introduces the Kriging interpolation algorithm by exploiting the fact that a depth image is also a distance image. In this study, the fractional-order calculus operator is introduced into the super-resolution reconstruction process of depth images, and the fractional-order differential mask is constructed to obtain high-frequency information such as image edges for originally low-resolution images. In addition, fractional-order edge matching is used to improve the accuracy of the Kriging interpolation process, and the fractional-order differential kernel is used to construct and improve the Kriging interpolation diffusion function. The experimental results show that the algorithm effectively improves the super-resolution reconstruction and enhances the subjective and objective outputs of the reconstruction process compared with the traditional algorithm.

2. Depth Image super-resolution Algorithm

Super-resolution of depth images requires image enlargement by interpolation of images. The main theoretical basis of such methods is that there is a purely sampling relationship between low-resolution images and high-resolution images. In the up-sampling process, only the pixel values of the low-resolution image need to be fitted to obtain a continuous image, and then sampling can be performed to obtain the target-resolution image. This method benefits from conceptual simplicity and low computational difficulty and is widely used for image super-resolution in industry. However, due to theoretical limitations, the super-resolution effect is not satisfactory when the images contain noise or blur, i.e., when the relationship between the high- and low-resolution images is not purely a sampling relationship.

2.1. Classical Image Interpolation Methods

OpenCV2′s resize function includes a number of commonly used interpolation algorithms: nearest-neighbour interpolation (INTER_NEAREST) [23], bilinear interpolation (INTER_LINEAR) [24,25], resampling interpolation (INTER_AREA) [26], cubic interpolation (INTER_CUBIC) [27], and LANCZOS4 cubic interpolation (INTER_LANCZOS4) [28]. Resampling interpolation is based on resampling using the pixel region relationship. Therefore, resampling is preferred in image extraction (decimation). Meanwhile, when the image is zoomed, INTER_NEAREST is preferred. These five classical interpolation algorithms are widely used in industry. However, recent years have seen the emergence of several new interpolation algorithms, including the fractal theory-based algorithm [29], the wavelet algorithm [30], the contour template algorithm [31], the partial differential equation-based algorithm [32], the topography-based inverse distance algorithm [33], and the Kriging algorithm [34]. These algorithms have been shown to have superior performance in various applications, such as image super-resolution, terrain modelling, and geological analysis, and therefore merit consideration for future research.

This study evaluates the effectiveness of the five classical interpolation algorithms mentioned at the start of the preceding paragraph with respect to the super-resolution of depth images. Using experimental data obtained from the depth image dataset in [35], the depth images and point clouds are first 50% down-sampled and then 200% up-sampled to examine the effectiveness of the interpolation algorithms. The super-resolution results are shown in Figure 1, Figure 2, Figure 3, Figure 4, Figure 5 and Figure 6 and

Table 1. Comparison of classical interpolation algorithms.

Interpolation Method	Nearest Neighbor Interpolation	Bilinear Interpolation	Resampling Interpolaton	Cubic Interpolation	Lanczos Interpolation
PSNR (dB)	29.795	34.154	33.570	34.515	33.654

, the images show good clarity and completeness after scaling.

The five classical interpolation algorithms are tested experimentally with respect to their ability to meet the requirements of image appearance. In terms of image edge completeness and clarity, the cubic interpolation and Lanczos interpolation algorithms achieve the highest values. However, after conversion of the depth images into point clouds, it is evident that the cubic interpolation and Lanczos interpolation algorithms introduce a large amount of smoothed data to optimise the appearance of the image. From the point cloud or depth distance perspectives, these data can be seen as a source of additional noise. Overall, the above experimental results show that the classical image interpolation algorithms, except the nearest-neighbour interpolation algorithm, inevitably introduce noise that cannot be eliminated. If the images need to be converted to point clouds and reconstructed in 3D, the noise introduced will greatly affect the final result and be difficult to eliminate. If high-quality, super-resolution depth images are to be obtained, improvements are needed, starting with the interpolation algorithm.

2.2. Interpolation Algorithm Based on Kriging Interpolation

A depth image is itself a distance image and is similar to a geological/topographic image after conversion to a point cloud. Considering the concept of inverse distance used in adaptive weighting, an algorithm suitable for geology/topography is introduced. The algorithm is based on inverse distance interpolation and Kriging interpolation, both of which take into account the spatial relationship between the point to be interpolated and the valid points around them [20,22].

Kriging interpolation is a method for the unbiased optimal estimation of regionalised variables in a finite region (used to estimate values that are spatially correlated, such as air quality, where values at closely spaced locations are similar). ‘Unbiased’ means that the expectation of the difference between the estimated and actual values is equal to zero, and ‘optimal’ means that the variance between the estimated and actual values is minimised. Based on this feature, Kriging interpolation is much more effective than other interpolation methods. In simple terms, Kriging interpolation does the following. Given some coordinates and the real values at those coordinates, which we call sampling points, these points are fed into the Kriging algorithm to estimate the values at other unknown locations [22].

Kriging interpolation analyses the regional variation of variables to achieve an optimal linear unbiased estimation [20,22]. The method consists of two steps: (1) the variational function model is established using the values of the spatial attributes of the regional variables, which are used to express the structural and stochastic characteristics of those attributes; (2) on the basis of the analysis of the spatial attributes, the domain range is determined and the Kriging equation is solved to estimate the values of the attributes with interpolation points.

$${\widehat{\mathrm{Z}}}_{\mathrm{o}}=\sum _{\mathrm{i}=1}^{\mathrm{n}}{\mathsf{\lambda }}_{\mathrm{i}}{\mathrm{z}}_{\mathrm{i}}$$

(1)

where ${\mathrm{z}}_{\mathrm{i}}$ is the estimated value at the point $\left(\mathrm{x},\mathrm{y}\right)$ and ${\mathsf{\lambda }}_{\mathrm{i}}$ is the weight coefficient, which is a weighted sum of the data at all known points in the space and is used to estimate the value of unknown points. However, the weight coefficients are not the inverse of the distance but are the optimal set of coefficients that can satisfy the minimum difference between the valuation ${\hat{\mathrm{z}}}_{\mathrm{o}}$ of the point $\left(\mathrm{x},\mathrm{y}\right)$ and the true value of ${\mathrm{z}}_{\mathrm{o}}$, i.e.,

$$\underset{{\mathsf{\lambda }}_{\mathrm{i}}}{\min }\mathrm{Var}\left({\hat{\mathrm{z}}}_{\mathrm{o}}-{\mathrm{z}}_{\mathrm{o}}\right)$$

(2)

Additionally, the unbiased condition must be satisfied:

$$\mathrm{E}\left({\hat{\mathrm{z}}}_{\mathrm{o}}-{\mathrm{z}}_{\mathrm{o}}\right)=0$$

(3)

The main difference between different Kriging interpolation methods is the underlying assumptions. In this paper, only the assumptions and applications of ordinary Kriging interpolation are presented.

Ordinary Kriging interpolation assumes that the spatial properties $\mathrm{z}$ are homogeneous. All points in space $\left(\mathrm{x},\mathrm{y}\right)$ have the same expectation c with variance ${\mathsf{\sigma }}^2$. That is, for any point $\left(\mathrm{x},\mathrm{y}\right)$, we have

$$\mathrm{E}\left[\mathrm{z}\left(\mathrm{x},\mathrm{y}\right)\right]=\mathrm{E}\left[\mathrm{z}\right]=\mathrm{c}$$

(4)

$$\mathrm{Var}\left[\mathrm{z}\left(\mathrm{x},\mathrm{y}\right)\right]{=\mathsf{\sigma }}^2$$

(5)

In other words, the value at any point z(x,y) is composed of the regional mean c and the random deviation at that point R(x,y):

$$\mathrm{z}\left(\mathrm{x},\mathrm{y}\right)=\mathrm{E}\left[\mathrm{z}\left(\mathrm{x},\mathrm{y}\right)\right]+\mathrm{R}\left(\mathrm{x},\mathrm{y}\right)=\mathrm{c}+\mathrm{R}\left(\mathrm{x},\mathrm{y}\right)$$

(6)

where 𝑅(x,y) denotes the deviations at the point (x,y), and their variances are all constants.

$$\mathrm{Var}\left[\mathrm{z}\left(\mathrm{x},\mathrm{y}\right)\right]{=\mathsf{\sigma }}^2$$

(7)

First, we analyse the unbiased estimation condition:

$$\mathrm{E}\left({\hat{\mathrm{z}}}_{\mathrm{o}}-{\mathrm{z}}_{\mathrm{o}}\right)=0$$

(8)

Substitute

$${\widehat{\mathrm{Z}}}_{\mathrm{o}}=\sum _{\mathrm{i}=1}^{\mathrm{n}}{\mathsf{\lambda }}_{\mathrm{i}}{\mathrm{z}}_{\mathrm{i}}$$

(9)

We obtain

$$\mathrm{E}\left({\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{n}}{\mathsf{\lambda }}_{\mathrm{i}}{\mathrm{z}}_{\mathrm{i}}{-\mathrm{z}}_{\mathrm{o}}\right)=0$$

(10)

Furthermore, because $\mathrm{E}\left[\mathrm{z}\right]=\mathrm{c}$ for any z, one of the constraints on ${\mathsf{\lambda }}_{\mathrm{i}}$ can be obtained.

$${\displaystyle \sum }_{\mathrm{I}=1}^{\mathrm{n}}{\mathsf{\lambda }}_{\mathrm{i}}=1$$

(11)

The analysis estimation error $\mathrm{Var}\left({\hat{\mathrm{z}}}_{\mathrm{o}}-{\mathrm{z}}_{\mathrm{o}}\right)$ is denoted here by J:

$$\mathrm{J}=\mathrm{Var}\left({\hat{\mathrm{z}}}_{\mathrm{o}}-{\mathrm{z}}_{\mathrm{o}}\right)$$

(12)

Furthermore, we have

$$\begin{gathered} \mathrm{J}=\mathrm{Var}\left({\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{n}}{\mathsf{\lambda }}_{\mathrm{i}}{\mathrm{z}}_{\mathrm{i}}{-\mathrm{z}}_{\mathrm{o}}\right) \\ =\mathrm{Var}\left({\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{n}}{\mathsf{\lambda }}_{\mathrm{i}}{\mathrm{z}}_{\mathrm{i}}-2\mathrm{Cov}\left({\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{n}}{\mathsf{\lambda }}_{\mathrm{i}}{\mathrm{z}}_{\mathrm{i}}{,\mathrm{z}}_{\mathrm{o}}\right){-\mathrm{z}}_{\mathrm{o}}\right)+\mathrm{Cov}\left({\mathrm{z}}_{\mathrm{o}}{,\mathrm{z}}_{\mathrm{o}}\right) \\ ={\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{n}}{\displaystyle \sum }_{\mathrm{j}=0}^{\mathrm{n}}{\mathsf{\lambda }}_{\mathrm{i}}{\mathsf{\lambda }}_{\mathrm{j}}\mathrm{Cov}\left({\mathrm{z}}_{\mathrm{i}}{,\mathrm{z}}_{\mathrm{j}}\right)-2{\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{n}}{\mathsf{\lambda }}_{\mathrm{i}}\mathrm{Cov}\left({\mathrm{z}}_{\mathrm{i}}{,\mathrm{z}}_{\mathrm{o}}\right)+\mathrm{Cov}\left({\mathrm{z}}_{\mathrm{o}}{,\mathrm{z}}_{\mathrm{o}}\right) \end{gathered}$$

(13)

By definition, ${\mathrm{C}}_{\mathrm{ij}}$=$\mathrm{Cov}\left({\mathrm{z}}_{\mathrm{i}},{\mathrm{z}}_{\mathrm{j}}\right)$=$\mathrm{Cov}\left({\mathrm{R}}_{\mathrm{i}},{\mathrm{R}}_{\mathrm{j}}\right)$, where ${\mathrm{R}}_{\mathrm{i}}={\mathrm{z}}_{\mathrm{i}}-\mathrm{c}$ denotes the deviation of the attribute value at point $\left({\mathrm{x}}_{\mathrm{i}},{\mathrm{y}}_{\mathrm{i}}\right)$ from the regional average attribute value

$$\mathrm{J}={\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{n}}{\displaystyle \sum }_{\mathrm{j}=0}^{\mathrm{n}}{\mathsf{\lambda }}_{\mathrm{i}}{\mathsf{\lambda }}_{\mathrm{j}}{\mathrm{C}}_{\mathrm{ij}}-2{\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{n}}{\mathsf{\lambda }}_{\mathrm{i}}{\mathrm{C}}_{\mathrm{io}}{+\mathrm{C}}_{\mathrm{oo}}$$

(14)

Defining the semivariance function ${\mathrm{r}}_{\mathrm{ij}}={\mathsf{\sigma }}^2-{\mathrm{C}}_{\mathrm{ij}}$, which is substituted into J, we have the following:

$$\begin{gathered} \mathrm{J}={\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{n}}{\displaystyle \sum }_{\mathrm{j}=0}^{\mathrm{n}}{\mathsf{\lambda }}_{\mathrm{i}}{\mathsf{\lambda }}_{\mathrm{j}}\left({\mathsf{\sigma }}^2{-\mathrm{r}}_{\mathrm{ij}}\right)-2{\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{n}}{\mathsf{\lambda }}_{\mathrm{i}}\left({\mathsf{\sigma }}^2{-\mathrm{r}}_{\mathrm{io}}\right)+\left({\mathsf{\sigma }}^2-{\mathrm{r}}_{\mathrm{oo}}\right) \\ ={\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{n}}{\displaystyle \sum }_{\mathrm{j}=0}^{\mathrm{n}}{\mathsf{\lambda }}_{\mathrm{i}}{\mathsf{\lambda }}_{\mathrm{j}}\left({\mathsf{\sigma }}^2\right)-{\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{n}}{\displaystyle \sum }_{\mathrm{j}=0}^{\mathrm{n}}{\mathsf{\lambda }}_{\mathrm{i}}{\mathsf{\lambda }}_{\mathrm{j}}\left({\mathrm{r}}_{\mathrm{ij}}\right)-2{\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{n}}{\mathsf{\lambda }}_{\mathrm{i}}\left({\mathsf{\sigma }}^2\right)+2{\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{n}}{\mathsf{\lambda }}_{\mathrm{i}}\left({\mathrm{r}}_{\mathrm{io}}\right)+\left({\mathsf{\sigma }}^2{-\mathrm{r}}_{\mathrm{oo}}\right) \end{gathered}$$

(15)

Given that ${{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}{\mathsf{\lambda }}_{\mathrm{i}}=1$, after substitution we have

$$\begin{gathered} {\mathrm{J}=\mathsf{\sigma }}^2-{\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{n}}{\displaystyle \sum }_{\mathrm{j}=0}^{\mathrm{n}}{\mathsf{\lambda }}_{\mathrm{i}}{\mathsf{\lambda }}_{\mathrm{j}}\left({\mathrm{r}}_{\mathrm{ij}}\right){-2\mathsf{\sigma }}^2+2{\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{n}}{\mathsf{\lambda }}_{\mathrm{i}}\left({\mathrm{r}}_{\mathrm{io}}\right)+\left({\mathsf{\sigma }}^2{-\mathrm{r}}_{\mathrm{oo}}\right) \\ =2{\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{n}}{\mathsf{\lambda }}_{\mathrm{i}}\left({\mathrm{r}}_{\mathrm{io}}\right)-{\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{n}}{\displaystyle \sum }_{\mathrm{j}=0}^{\mathrm{n}}{\mathsf{\lambda }}_{\mathrm{i}}{\mathsf{\lambda }}_{\mathrm{j}}\left({\mathrm{r}}_{\mathrm{ij}}\right){-\mathrm{r}}_{\mathrm{oo}} \end{gathered}$$

(16)

The goal is to find a set of ${\mathsf{\lambda }}_{\mathrm{i}}$ that minimises $\mathrm{J}$. As $\mathrm{J}$ is a function of ${\mathsf{\lambda }}_{\mathrm{i}}$, we can directly calculate the partial derivative of $\mathrm{J}$ with respect to ${\mathsf{\lambda }}_{\mathrm{i}}$, and find the zero-valued point. That is, we solve the following for ${\mathsf{\lambda }}_{\mathrm{i}}$:

$$\frac{\partial \mathrm{J}}{\partial {\mathsf{\lambda }}_{\mathrm{i}}}=0;\mathrm{i}=1,2,3,\dots \mathrm{n}$$

(17)

It is also necessary to satisfy ${\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{n}}{\mathsf{\lambda }}_{\mathrm{i}}=1$. Additionally, the optimisation problem with this constraint is solved using the Lagrange multiplier method by first constructing a new objective function.

$$\mathrm{J}+2\varphi \left({\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{n}}{\mathsf{\lambda }}_{\mathrm{i}}-1\right)$$

(18)

where $\varphi$ is the Lagrangian multiplier. We solve for the set of parameters $\varphi ,{\mathsf{\lambda }}_1,{\mathsf{\lambda }}_2,\cdots ,{\mathsf{\lambda }}_{\mathrm{n}}$ that minimises this cost function and thereby minimises J under the constraint ${\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{n}}{\mathsf{\lambda }}_{\mathrm{i}}=1$, which yields

$$\{\begin{array}{l}\frac{\partial \left(\mathrm{J}+2\varphi \left({{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}{\mathsf{\lambda }}_{\mathrm{i}}-1\right)\right)}{\partial {\mathsf{\lambda }}_{\mathrm{i}}}=0;\mathrm{i}=1,2,3,\dots \mathrm{n}\\ \frac{\partial \left(\mathrm{J}+2\varphi \left({{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}{\mathsf{\lambda }}_{\mathrm{i}}-1\right)\right)}{\partial \varphi }=0\end{array}$$

(19)

Substituting into J, we obtain

$$\{\begin{array}{c}\frac{\partial \left(2{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}{\mathsf{\lambda }}_{\mathrm{i}}\left({\mathrm{r}}_{\mathrm{io}}\right)-{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}{{\displaystyle \sum }}_{\mathrm{j}=0}^{\mathrm{n}}{\mathsf{\lambda }}_{\mathrm{i}}{\mathsf{\lambda }}_{\mathrm{j}}\left({\mathrm{r}}_{\mathrm{ij}}\right){-\mathrm{r}}_{\mathrm{oo}}+2\varphi \left({{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}{\mathsf{\lambda }}_{\mathrm{i}}-1\right)\right)}{\partial {\mathsf{\lambda }}_{\mathrm{i}}}=0;\\ \frac{\partial \left(2{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}{\mathsf{\lambda }}_{\mathrm{i}}\left({\mathrm{r}}_{\mathrm{io}}\right)-{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}{{\displaystyle \sum }}_{\mathrm{j}=0}^{\mathrm{n}}{\mathsf{\lambda }}_{\mathrm{i}}{\mathsf{\lambda }}_{\mathrm{j}}\left({\mathrm{r}}_{\mathrm{ij}}\right){-\mathrm{r}}_{\mathrm{oo}}+2\varphi \left({{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}{\mathsf{\lambda }}_{\mathrm{i}}-1\right)\right)}{\partial \varphi }=0\\ (\mathrm{i}=1,2,3,\text{…}\mathrm{n})\\ \end{array}$$

(20)

$$\{\begin{array}{c}\frac{\partial \left(2{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}{\mathsf{\lambda }}_{\mathrm{i}}\left({\mathrm{r}}_{\mathrm{io}}\right)-{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}{{\displaystyle \sum }}_{\mathrm{j}=0}^{\mathrm{n}}{\mathsf{\lambda }}_{\mathrm{i}}{\mathsf{\lambda }}_{\mathrm{j}}\left({\mathrm{r}}_{\mathrm{ij}}\right){-\mathrm{r}}_{\mathrm{oo}}+2\varphi \left({{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}{\mathsf{\lambda }}_{\mathrm{i}}-1\right)\right)}{\partial {\mathsf{\lambda }}_{\mathrm{i}}}=0;\\ \frac{\partial \left(2{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}{\mathsf{\lambda }}_{\mathrm{i}}\left({\mathrm{r}}_{\mathrm{io}}\right)-{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}{{\displaystyle \sum }}_{\mathrm{j}=0}^{\mathrm{n}}{\mathsf{\lambda }}_{\mathrm{i}}{\mathsf{\lambda }}_{\mathrm{j}}\left({\mathrm{r}}_{\mathrm{ij}}\right){-\mathrm{r}}_{\mathrm{oo}}+2\varphi \left({{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}{\mathsf{\lambda }}_{\mathrm{i}}-1\right)\right)}{\partial \varphi }=0\\ \mathrm{i}=1,2,3,\cdots \mathrm{n}\end{array}$$

(21)

For the reason that ${\mathrm{C}}_{\mathrm{ij}}=\mathrm{Cov}\left({\mathrm{z}}_{\mathrm{i}},{\mathrm{z}}_{\mathrm{j}}\right)={\mathrm{C}}_{\mathrm{ji}}$, the same semivariance function ${\mathrm{r}}_{\mathrm{ij}}={\mathrm{r}}_{\mathrm{ij}}$, and further simplification shows that

$$\{\begin{array}{l}{\mathrm{r}}_{\mathrm{io}}-{\displaystyle {\displaystyle \sum }_{\mathrm{j}=1}^{\mathrm{n}}}{\mathrm{r}}_{\mathrm{ij}}{\mathsf{\lambda }}_{\mathrm{j}}+\varphi =0;\mathrm{i}=1,2,3,\cdots \mathrm{n}\\ {\displaystyle {\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{n}}}{\mathsf{\lambda }}_{\mathrm{i}}=1\end{array}$$

(22)

The equations to be solved for the above weighting factor ${\mathsf{\lambda }}_{\mathrm{i}}$ are first written in the form of a linear system of equations, as follows.

$$\{\begin{array}{c}{\mathrm{r}}_{11}{\mathsf{\lambda }}_1{+\mathrm{r}}_{12}{\mathsf{\lambda }}_2+\cdots {\mathrm{r}}_{1\mathrm{n}}{\mathsf{\lambda }}_{\mathrm{n}}-\varphi {=\mathrm{r}}_{1\mathrm{o}}\\ {\mathrm{r}}_{21}{\mathsf{\lambda }}_1{+\mathrm{r}}_{22}{\mathsf{\lambda }}_2+\cdots {\mathrm{r}}_{2\mathrm{n}}{\mathsf{\lambda }}_{\mathrm{n}}-\varphi {=\mathrm{r}}_{2\mathrm{o}}\\ ⋮\\ {\mathrm{r}}_{\mathrm{n}1}{\mathsf{\lambda }}_1{+\mathrm{r}}_{\mathrm{n}2}{\mathsf{\lambda }}_2+\cdots {\mathrm{r}}_{\mathrm{nn}}{\mathsf{\lambda }}_{\mathrm{n}}-\varphi {=\mathrm{r}}_{\mathrm{no}}\\ {\mathsf{\lambda }}_1{+\mathsf{\lambda }}_2+\cdots {+\mathsf{\lambda }}_{\mathrm{n}}=1\end{array}$$

(23)

The system of equations is written in matrix form, and the inverse matrix solution procedure is performed.

By the definition of semi-covariance, ${\mathrm{r}}_{\mathrm{ij}}={\mathsf{\sigma }}^2-{\mathrm{C}}_{\mathrm{ij}}$ and its equivalent form ${\mathrm{r}}_{\mathrm{ij}}=\frac{1}{2}\mathrm{E}\left[{\left({\mathrm{z}}_{\mathrm{i}}-{\mathrm{z}}_{\mathrm{j}}\right)}^2\right]$.

By the previous definition (hypothetical conditions), we have $\mathrm{Z}\left(\mathrm{x},\mathrm{y}\right)=\mathrm{c}+\mathrm{R}\left(\mathrm{x},\mathrm{y}\right)$, so we then have ${\mathrm{R}}_{\mathrm{i}}={\mathrm{z}}_{\mathrm{i}}-\mathrm{c}$. Further, we have ${\mathrm{z}}_{\mathrm{i}}-{\mathrm{z}}_{\mathrm{j}}={\mathrm{R}}_{\mathrm{i}}-{\mathrm{R}}_{\mathrm{j}}$ under these conditions:

$$\begin{gathered} {\mathrm{r}}_{\mathrm{ij}}=\frac{1}{2}\mathrm{E}\left[{\left({\mathrm{R}}_{\mathrm{i}}{-\mathrm{R}}_{\mathrm{j}}\right)}^2\right] \\ =\frac{1}{2}\mathrm{E}\left[{\mathrm{R}}_{\mathrm{i}}{}^2{-2\mathrm{R}}_{\mathrm{ij}}{+\mathrm{R}}_{\mathrm{j}}{}^2\right] \\ =\frac{1}{2}{\mathrm{E}[\mathrm{R}}_{\mathrm{i}}{}^2]+\frac{1}{2}{\mathrm{E}[\mathrm{R}}_{\mathrm{j}}{}^2]-\mathrm{E}\left[{\mathrm{R}}_{\mathrm{i}}{\mathrm{R}}_{\mathrm{j}}\right] \end{gathered}$$

(24)

This follows from the fact that

$${\mathrm{E}[\mathrm{R}}_{\mathrm{i}}{}^2{] = \mathrm{E}[\mathrm{R}}_{\mathrm{j}}{}^2{] = \mathrm{E}[(\mathrm{z}}_{\mathrm{i}}{-\mathrm{c})}^2{] = \mathrm{Var}(\mathrm{z}}_{\mathrm{i}}{) = \mathsf{\sigma }}^2$$

(25)

$$\mathrm{E}\left[{\mathrm{R}}_{\mathrm{i}}{\mathrm{R}}_{\mathrm{j}}\right]=\mathrm{E}\left[\left({\mathrm{z}}_{\mathrm{i}}-\mathrm{c}\right)\left({\mathrm{z}}_{\mathrm{j}}-\mathrm{c}\right)\right]=\mathrm{Cov}\left({\mathrm{z}}_{\mathrm{i}}{,\mathrm{z}}_{\mathrm{j}}\right){=\mathrm{C}}_{\mathrm{ij}}$$

(26)

After substitution, we obtain

$${\mathrm{r}}_{\mathrm{ij}} = \frac{1}{2}{\mathsf{\sigma }}^2+\frac{1}{2}{\mathsf{\sigma }}^2{-\mathrm{C}}_{\mathrm{ij}}{ = \mathsf{\sigma }}^2{-\mathrm{C}}_{\mathrm{ij}}$$

(27)

The above proof confirms that ${\mathrm{r}}_{\mathrm{ij}}=\frac{1}{2}\mathrm{E}\left[{\left({\mathrm{z}}_{\mathrm{i}}-{\mathrm{z}}_{\mathrm{j}}\right)}^2\right]$. The following procedure is then used to calculate ${\mathrm{r}}_{\mathrm{ij}}$. Here, we continue to use the principle of similarity of attributes in space. This spatial similarity is expressed by ${\mathrm{r}}_{\mathrm{ij}}=\frac{1}{2}{\left({\mathrm{z}}_{\mathrm{i}}-{\mathrm{z}}_{\mathrm{j}}\right)}^2$ in terms of distance, defining the geometric distance between ${\mathrm{z}}_{\mathrm{i}}$ and ${\mathrm{z}}_{\mathrm{j}}$ as

$${\mathrm{d}}_{\mathrm{ij}}{=\mathrm{d}(\mathrm{z}}_{\mathrm{i}}{,\mathrm{z}}_{\mathrm{j}}{)=\mathrm{d}((\mathrm{x}}_{\mathrm{i}}{,\mathrm{y}}_{\mathrm{i}}{),(\mathrm{x}}_{\mathrm{j}}{,\mathrm{y}}_{\mathrm{j}}\left)\right)=\sqrt{{\left({\mathrm{x}}_{\mathrm{i}}{-\mathrm{x}}_{\mathrm{j}}\right)}^2+{\left({\mathrm{y}}_{\mathrm{i}}{-\mathrm{y}}_{\mathrm{j}}\right)}^2}$$

(28)

Kriging interpolation assumes that there is a functional relationship between ${\mathrm{r}}_{\mathrm{ij}}$ and ${\mathrm{d}}_{\mathrm{ij}}$, which can be linear, quadratic, exponential or logarithmic. To identify this relationship, we perform the following steps on the observed data set $\left\{\mathrm{z}({\mathrm{x}}_1,{\mathrm{y}}_1),\mathrm{z}({\mathrm{x}}_2,{\mathrm{y}}_2),\cdots ,\mathrm{z}({\mathrm{x}}_{\mathrm{n}-1},{\mathrm{y}}_{\mathrm{n}-1}),\mathrm{z}({\mathrm{x}}_{\mathrm{n}},{\mathrm{y}}_{\mathrm{n}}),\right\}$.

First, we calculate the distance between any two points

$${\mathrm{d}}_{\mathrm{ij}}=\sqrt{{\left({\mathrm{x}}_{\mathrm{i}}{-\mathrm{x}}_{\mathrm{j}}\right)}^2+{\left({\mathrm{y}}_{\mathrm{i}}{-\mathrm{y}}_{\mathrm{j}}\right)}^2}$$

(29)

and semi-covariance

$${\mathrm{r}}_{\mathrm{ij}}=\frac{1}{2}{\left({\mathrm{z}}_{\mathrm{i}}{-\mathrm{z}}_{\mathrm{j}}\right)}^2$$

(30)

It is possible to obtain ${\mathrm{n}}^2$ (d_ij, r_ij) data pairs by plotting all the d and r values in a scatter plot to find an optimal curve to fit d and r and obtain a functional relationship.

$$\mathrm{r}=\mathrm{r}\left(\mathrm{d}\right)$$

(31)

Furthermore, for any two points $\mathrm{z}({\mathrm{x}}_{\mathrm{i}},{\mathrm{y}}_{\mathrm{i}}),\mathrm{z}({\mathrm{x}}_{\mathrm{j}},{\mathrm{y}}_{\mathrm{j}})$, we first calculate their distance, d_ij. Next, according to the functional relationship, we obtain the semi-covariance r_ij. Considering the computational complexity, simple Kriging is used for the interpolation of depth images. The ordinary Kriging formula is as follows.

$${\widehat{\mathrm{Z}}}_{\mathrm{o}}=\sum _{\mathrm{i}=1}^{\mathrm{n}}{\mathsf{\lambda }}_{\mathrm{i}}{\mathrm{z}}_{\mathrm{i}}$$

(32)

The simple Kriging formula is as follows.

$${\widehat{\mathrm{Z}}}_{\mathrm{o}}-\mathrm{c}=\sum _{\mathrm{i}=1}^{\mathrm{n}}{\mathsf{\lambda }}_{\mathrm{i}}{\mathrm{z}}_{\mathrm{i}}-\mathrm{c}$$

(33)

where c is the mathematical expectation of the attribute value, i.e., $\mathrm{E}\left[\mathrm{z}\right]=\mathrm{c}$. That is, in ordinary Kriging interpolation, the attribute value of an unknown point is considered to be a weighted sum of the attribute values of known points.

The simple Kriging Formula (33) completes the derivation of the estimation of the values of depth images by interpolation, according to which different peripheral depth values can be selected for interpolation. Based on this theoretical derivation, a depth image interpolation algorithm is constructed to super-resolve depth images by Kriging interpolation. The method first down-samples the image by 50% and then up-samples it by 200%, where the sampling algorithm uses Kriging interpolation. The super-resolution results are shown in Figure 7 which shows that the Kriging interpolation algorithm is affected by edges and voids. It can also be seen from the output of the Kriging interpolation algorithm that rapid changes at edges can lead to highly discrepant interpolation results from the actual situation.

3. Fractional-Order Differential Edge Recognition

Image enhancement processing is an important branch of digital image processing. Under real-world environmental conditions, the visual effectiveness of image capture is poor, which necessitates image enhancement techniques. Such techniques include highlighting certain characteristics of the target object in the image and extracting the characteristic parameters of the target object from the digital image. The goal is to improve the recognition, tracking, and understanding of the target [36]. The main objective of image enhancement processing is to highlight the parts of interest in an image and attenuate or remove unwanted information so that the useful information is enhanced to obtain a more useful image, or the image can be made more suitable for analysis and processing by humans or machines. The main features of an image (i.e., the main information) include edges, texture details, and other features. The enhancement of image features plays an important role in improving the visual effectiveness of an image and also lays a good foundation for subsequent image processing (e.g., image recognition, image segmentation, or super resolution). Fractional-order differentiation significantly enhances the high- and medium-frequency components of the signal while nonlinearly strengthening the medium-frequency components of the signal to a certain extent as well as nonlinearly preserving the low-frequency and DC components of the signal. Therefore, fractional-order differentiation is beneficial in terms of both significantly enhancing the image edge and contour information and strengthening the texture details in the smooth regions of the image. Due to this property of fractional-order differentiation, the fractional calculus theory in image processing operations has received increasing attention [13,14,15,16,17].

It is generally believed that the main benefit of applying methods based on fractional-order calculus to digital images is to add a degree of freedom, i.e., the fractional-order v (0 < [v] < 1). By choosing the appropriate order and constructing the corresponding mask operator, better image enhancement (v > 0) and denoising (v < 0) can be achieved. Huang G et al. derived fractional-order calculus and verified the application of the fractional-order integral in image processing [37]. They also proposed a fractional-order differential enhancement and fractional-order integral denoising model based on the Grünwald–Letnikov (G-L) definition of 2D images.

Introducing the G-L definition of the fractional-order calculus, the $-v$ order integral of the G-L expression [17] is given by the following equation:

$${}_{\mathrm{a}}{}^{\mathrm{G}}\mathrm{D}_{\mathrm{x}}^{-\mathsf{\nu }}=\underset{\mathrm{h}\to 0}{\lim }{\mathrm{h}}^{\mathsf{\nu }}{{\displaystyle \sum }}_{\mathrm{n}=0}^{\frac{\mathrm{t}-\mathrm{a}}{\mathrm{h}}}\left(\begin{array}{c}-\mathrm{v}\\ \mathrm{j}\end{array}\right)\mathrm{f}\left(\mathrm{x}-\mathrm{jh}\right),\mathrm{v}\in \mathrm{R}$$

(34)

The fractional-order integral can be approximately expanded by defining the one-dimensional (1D) signal f(t) of the duration interval as [a,t]. According to the unit h = 1 equipartition, we have $\mathrm{m}=\left[\frac{\mathrm{t}-\mathrm{a}}{\mathrm{h}}\right]=\left[\mathrm{t}-\mathrm{a}\right]$, implying that there is a 1D signal of v (v ≥ 0). The equivalent expression of the fractional-order differentiation of order $v$ is

$$\frac{{\mathrm{d}}^{\mathrm{v}}\mathrm{t}}{{\mathrm{dt}}^{\mathrm{v}}}\approx \mathrm{f}\left(\mathrm{t}\right)+{\left(-1\right)}^1\left(\mathrm{v}\right)\mathrm{f}\left(\mathrm{t}-1\right)+{\left(-1\right)}^2\left(\frac{\mathrm{v}(\mathrm{v}+1)}{2}\right)\mathrm{f}\left(\mathrm{t}-2\right)+\dots +{\left(-1\right)}^{\mathrm{j}}\frac{\mathsf{\Gamma }\left(\mathsf{\nu }+1\right)}{\mathsf{\Gamma }\left(\mathrm{j}+1\right)\mathsf{\Gamma }\left(\mathsf{\nu }-\mathrm{j}+1\right)}\mathrm{f}\left(\mathrm{t}-\mathrm{j}\right)$$

(35)

Define the 2D image signal $\mathrm{I}\left(\mathrm{x},\mathrm{y}\right)$, and assume that the fractional-order derivatives of $\mathrm{I}\left(\mathrm{x},\mathrm{y}\right)$ are separable for the x-axis and y-axis directions under certain conditions. Furthermore, we use the separability of the Fourier transform to extend the fractional-order calculus from one dimension to two dimensions. By equating the 2D image signal $\mathrm{I}\left(\mathrm{x},\mathrm{y}\right)$ according to the unit time with h = 1, the fractional-order calculus of the x-axis and y-axis can be obtained, and v is the fractional-order calculus order (v ≥ 0).

From Equation (35), the approximate solution of partial fractional-order calculus under the G-L definition can be derived, enabling us to solve the numerical expression of the fractional-order calculus operator along the x-axis and y-axis directions, namely

$$\frac{{\mathrm{d}}^{\mathrm{v}}\mathrm{I}(\mathrm{x},\mathrm{y})}{{\mathrm{dx}}^{\mathrm{v}}}\approx \mathrm{I}(\mathrm{x},\mathrm{y})+{\left(-1\right)}^1\left(\mathrm{v}\right)\mathrm{I}(\mathrm{x}-1,\mathrm{y})+{\left(-1\right)}^2\left(\frac{\mathrm{v}(\mathrm{v}-1)}{2}\right)\mathrm{I}(\mathrm{x}-2,\mathrm{y})+\text{…}+{\left(-1\right)}^{\mathrm{j}}\frac{\mathsf{\Gamma }\left(\mathsf{\nu }+1\right)}{\mathsf{\Gamma }\left(\mathrm{j}+1\right)\mathsf{\Gamma }\left(\mathsf{\nu }-\mathrm{j}+1\right)}\mathrm{I}(\mathrm{x}-\mathrm{j},\mathrm{y})$$

(36)

$$\frac{{\mathrm{d}}^{\mathrm{v}}\mathrm{I}(\mathrm{x},\mathrm{y})}{{\mathrm{dy}}^{\mathrm{v}}}\approx \mathrm{I}(\mathrm{x},\mathrm{y})+{\left(-1\right)}^1\left(\mathrm{v}\right)\mathrm{I}(\mathrm{x},\mathrm{y}-1)+{\left(-1\right)}^2\left(\frac{\mathrm{v}(\mathrm{v}-1)}{2}\right)\mathrm{I}(\mathrm{x},\mathrm{y}-2)+\text{…}+{\left(-1\right)}^{\mathrm{j}}\frac{\mathsf{\Gamma }\left(\mathsf{\nu }+1\right)}{\mathsf{\Gamma }\left(\mathrm{j}+1\right)\mathsf{\Gamma }\left(\mathsf{\nu }-\mathrm{j}+1\right)}\mathrm{I}(\mathrm{x},\mathrm{y}-\mathrm{j})$$

(37)

That is,

$$\frac{{\mathrm{d}}^{\mathrm{v}}\mathrm{I}(\mathrm{x},\mathrm{y})}{{\mathrm{dx}}^{\mathrm{v}}}=\underset{\mathrm{N}\to \infty }{\lim }\left[{{\displaystyle \sum }}_{\mathrm{m}=0}^{\mathrm{N}-1}{\left(-1\right)}^{\mathrm{j}}\frac{\mathsf{\Gamma }\left(\mathsf{\nu }+1\right)}{\mathsf{\Gamma }\left(\mathrm{j}+1\right)\mathsf{\Gamma }\left(\mathsf{\nu }-\mathrm{j}+1\right)}\mathrm{I}(\mathrm{x}-\mathrm{j},\mathrm{y})\right]$$

(38)

$$\frac{{\mathrm{d}}^{\mathrm{v}}\mathrm{I}(\mathrm{x},\mathrm{y})}{{\mathrm{dy}}^{\mathrm{v}}}=\underset{\mathrm{N}\to \infty }{\lim }\left[{{\displaystyle \sum }}_{\mathrm{m}=0}^{\mathrm{N}-1}{\left(-1\right)}^{\mathrm{j}}\frac{\mathsf{\Gamma }\left(\mathsf{\nu }+1\right)}{\mathsf{\Gamma }\left(\mathrm{j}+1\right)\mathsf{\Gamma }\left(\mathsf{\nu }-\mathrm{j}+1\right)}\mathrm{I}(\mathrm{x},\mathrm{y}-\mathrm{j})\right]$$

(39)

From Equations (38) and (39), we can obtain v (v ≥ 0). The coefficients of the fractional-order differential operator of order 𝑅 are as follows.

$${\mathrm{R}=(-1)}^{\mathrm{j}}\left(\begin{array}{c}\mathrm{v}\\ \mathrm{j}\end{array}\right){=(-1)}^{\mathrm{j}}\frac{\mathsf{\Gamma }(\mathsf{\nu }+1)}{\mathsf{\Gamma }(\mathrm{j}+1)\mathsf{\Gamma }(\mathsf{\nu }-\mathrm{j}+1)}$$

(40)

Assuming a mask size of $3\times 3$, i.e., setting N = 3, the approximate solutions in both axis directions can be obtained from Equations (38) and (39).

$$\frac{{\mathrm{d}}^{\mathrm{v}}\mathrm{I}(\mathrm{x},\mathrm{y})}{{\mathrm{dx}}^{\mathrm{v}}}\approx {\mathrm{I}(\mathrm{x},\mathrm{y})+(-1)}^1{\left(\mathrm{v}\right)\mathrm{I}(\mathrm{x}-1,\mathrm{y})+(-1)}^2(\frac{\mathrm{v}(\mathrm{v}-1)}{2}\left)\mathrm{I}\right(\mathrm{x}-2,\mathrm{y})$$

(41)

$$\frac{{\mathrm{d}}^{\mathrm{v}}\mathrm{I}(\mathrm{x},\mathrm{y})}{{\mathrm{dy}}^{\mathrm{v}}}\approx {\mathrm{I}(\mathrm{x},\mathrm{y})+(-1)}^1{\left(\mathrm{v}\right)\mathrm{I}(\mathrm{x},\mathrm{y}-1)+(-1)}^2(\frac{\mathrm{v}(\mathrm{v}-1)}{2}\left)\mathrm{I}\right(\mathrm{x},\mathrm{y}-2)$$

(42)

Considering that the other six directions in the image have similar definitions, the approximate solution of the fractional-order differential in those six directions can be derived. Thus, a fractional-order differential operator filter with 8 directions can be obtained, and the operator is rotationally invariant.

From Equations (40)–(42), we can obtain the fractional-order differential mask and a 5 × 5 fractional-order differential mask template v (v ≥ 0). The coefficients of the fractional-order differential convolution mask are

$${\mathrm{a}}_{*}^{\mathsf{\nu }}{=[\mathrm{R}}_0^{\mathrm{v}}{,\mathrm{R}}_1^{\mathrm{v}}{,\mathrm{R}}_2^{\mathrm{v}}],(\mathrm{j}=0,1,2)$$

(43)

where

$${\mathrm{R}}_0^{\mathrm{v}}=1$$

(44)

$${\mathrm{R}}_1^{\mathrm{v}}=-\mathrm{v}$$

(45)

$${\mathrm{R}}_2^{\mathrm{v}}=\frac{\mathrm{v}(\mathrm{v}-1)}{2}$$

(46)

The definition mask template has 8 directions, where by definition the x-axis positive and negative directions have $a_0^v$ and $a_{180}^v$. In the anticlockwise direction, there are $a_{45}^v$, $a_{135}^v$, $a_{225}^v,$ and $a_{315}^v$, and the positive and negative directions of the y-axis have $a_{90}^v$ and $a_{270}^v$. The constructed masks are as follows.

Where the coefficients are as follows

$$a_0^v{=\mathrm{R}}_0^{\mathrm{v}}=1$$

(47)

$$a_1^v{=\mathrm{R}}_1^{\mathrm{v}}=-\mathrm{v}$$

(48)

$$a_2^v{=\mathrm{R}}_2^{\mathrm{v}}=\frac{\mathrm{v}(\mathrm{v}-1)}{2}$$

(49)

In any image, there is a certain similarity between neighbouring pixels. To reduce the unnecessary spatial and temporal complexity, the feature information of the local neighbourhood of the target pixel should be fully utilised, considering that the pixels closer to the target have higher similarity to the target. Given 8 directions of 3 × 3 masks a ∗ v, for a 5 × 5 masks of image points I(x,y), the convolution calculation is performed as follows:

$$\mathrm{I}{\left(\mathrm{x},\mathrm{y}\right)}_{*}^{\mathsf{\nu }}=\mathrm{I}\left(\mathrm{x},\mathrm{y}\right){*\mathrm{a}}_{*}^{\mathsf{\nu }}$$

(50)

The result of the convolution calculation is linearly weighted by the sum of the convolutions in each direction to obtain the final result:

$$\begin{gathered} \mathrm{I}{\left(\mathrm{x},\mathrm{y}\right)}_{*}^{\mathsf{\nu }}=\frac{\mathrm{I}{\left(\mathrm{x},\mathrm{y}\right)}_0^{\mathsf{\nu }}}{\mathrm{sum}(\mathrm{x},\mathrm{y})}\times \mathrm{I}{\left(\mathrm{x},\mathrm{y}\right)}_0^{-\mathsf{\nu }}+\frac{\mathrm{I}{\left(\mathrm{x},\mathrm{y}\right)}_{45}^{\mathsf{\nu }}}{\mathrm{sum}(\mathrm{x},\mathrm{y})}\times \mathrm{I}{\left(\mathrm{x},\mathrm{y}\right)}_{45}^{\mathsf{\nu }}+\frac{\mathrm{I}{\left(\mathrm{x},\mathrm{y}\right)}_{90}^{\mathrm{v}}}{\mathrm{sum}(\mathrm{x},\mathrm{y})}\times \mathrm{I}{\left(\mathrm{x},\mathrm{y}\right)}_{90}^{\mathsf{\nu }} \\ +\frac{\mathrm{I}{\left(\mathrm{x},\mathrm{y}\right)}_{135}^{\mathsf{\nu }}}{\mathrm{sum}(\mathrm{x},\mathrm{y})}\times \mathrm{I}{\left(\mathrm{x},\mathrm{y}\right)}_{135}^{\mathsf{\nu }}+\frac{\mathrm{I}{\left(\mathrm{x},\mathrm{y}\right)}_{180}^{\mathsf{\nu }}}{\mathrm{sum}(\mathrm{x},\mathrm{y})}\times \mathrm{I}{\left(\mathrm{x},\mathrm{y}\right)}_{180}^{\mathsf{\nu }}+\frac{\mathrm{I}{\left(\mathrm{x},\mathrm{y}\right)}_{225}^{\mathsf{\nu }}}{\mathrm{sum}(\mathrm{x},\mathrm{y})}\times \mathrm{I}{\left(\mathrm{x},\mathrm{y}\right)}_{225}^{\mathsf{\nu }} \\ +\frac{\mathrm{I}{\left(\mathrm{x},\mathrm{y}\right)}_{270}^{\mathsf{\nu }}}{\mathrm{sum}(\mathrm{x},\mathrm{y})}\times \mathrm{I}{\left(\mathrm{x},\mathrm{y}\right)}_{270}^{\mathsf{\nu }}+\frac{\mathrm{I}{\left(\mathrm{x},\mathrm{y}\right)}_{315}^{\mathsf{\nu }}}{\mathrm{sum}(\mathrm{x},\mathrm{y})}\times \mathrm{I}{\left(\mathrm{x},\mathrm{y}\right)}_{315}^{\mathsf{\nu }} \end{gathered}$$

(51)

where

$$\begin{gathered} \mathrm{sum}\left(\mathrm{x},\mathrm{y}\right)=\mathrm{I}{\left(\mathrm{x},\mathrm{y}\right)}_0^{\mathsf{\nu }}+\mathrm{I}{\left(\mathrm{x},\mathrm{y}\right)}_{45}^{\mathsf{\nu }}+\mathrm{I}{\left(\mathrm{x},\mathrm{y}\right)}_{90}^{\mathsf{\nu }}+\mathrm{I}{\left(\mathrm{x},\mathrm{y}\right)}_{135}^{\mathsf{\nu }} \\ +\mathrm{I}{\left(\mathrm{x},\mathrm{y}\right)}_{180}^{\mathsf{\nu }}+\mathrm{I}{\left(\mathrm{x},\mathrm{y}\right)}_{225}^{\mathsf{\nu }}+\mathrm{I}{\left(\mathrm{x},\mathrm{y}\right)}_{270}^{\mathsf{\nu }}+\mathrm{I}{\left(\mathrm{x},\mathrm{y}\right)}_{315}^{\mathsf{\nu }} \end{gathered}$$

(52)

After the fractional-order differential mask convolution calculation, there is a continuous process of enlargement or reduction. Hence, it is necessary to construct the fractional-order differential normalisation factor q to unify the final enlargement and reduction factors. The factor is calculated as follows:

$$\mathrm{q}={\displaystyle \sum }_{\mathsf{\theta }=0}^{360}{\mathrm{I}(\mathrm{x},\mathrm{y})}_{\mathsf{\theta }}^{\mathrm{v}}(\mathsf{\theta }=0,45,90,135,180,215,270,315)$$

(53)

Importing q into Equation (54), the new $\mathrm{I}{\left(\mathrm{x},\mathrm{y}\right)}_{\mathrm{m}}^{\mathsf{\nu }}$ after fractional-order differential filtering of order ν using a 5 × 5 mask is obtained as follows.

$$\mathrm{I}{\left(\mathrm{x},\mathrm{y}\right)}_{\mathrm{m}}^{\mathsf{\nu }}=\frac{\mathrm{I}\left(\mathrm{x},\mathrm{y}\right){*\mathrm{a}}_{*}^{\mathsf{\nu }}}{\mathrm{q}}$$

(54)

In 2D images, both edges and noise are discontinuous points with local characteristics, and both have rapid changes in the pixel values of the corresponding domain, so they are high-frequency signals. The edges have the characteristics of orderliness, directionality, and structure, and they have a larger energy and range relative to the noise. Edge detection models based on fractional-order differentiation have shown good results. The authors of [13] recommended the use of v = 0.6 in fractional-order differentiation based on their experiments, which enabled them to achieve richer edge information compared with the classical integer-order edge detection operator. From the depth image and the corresponding point cloud map after the fractional-order differential mask convolution calculation in Figure 8. It can be seen from Figure 9 that the fractional-order differential method can identify the edges accurately.

4. Kriging Interpolation Super-Resolution Algorithm Based on Fractional-Order Differentiation

Depth images consist of ordered distance data distributed spatially, and the orderliness makes them conducive to parallel computation, a technology that may realise efficient filtering and convolution computation in the future.

Depth images also contain peripheral correlation, and the high-frequency points of edges can be identified by fractional-order differentiation. The high-frequency points are interpolated using the spatial local neighbourhood, while the non-high-frequency points are interpolated by Kriging interpolation. The flowchart in Figure 10 depicts the fractional-order differential Kriging interpolation super-resolution model, which is divided into the following processes: acquiring the depth image, down-sampling the depth image, performing fractional-order differential edge detection, interpolating the down-sampled depth map using the fractional-order differential Kriging interpolation super-resolution model, and obtaining the high-resolution depth image.

Table 2. Fractional-order differential Kriging interpolation super-resolution (200%).

Datasets [35] Number	Fractional-Order Differential Kriging Interpolation		Nearest Neighbor Interpolation
	$\mathit{v}$	PSNR	PSNR
00033	0.6	33.457	29.795
03236	0.6	37.230	27.738
03528	0.6	38.598	29.868
04797	0.6	47.182	26.580
05989	0.6	39.888	29.468
08343	0.6	38.391	25.008

shows the comparison between the kriging interpolation algorithm based on fractional order differentiation and the nearest neighbor interpolation algorithm. It can be seen that the PSNR of kriging interpolation based on fractional order differentiation is better than the nearest neighbor interpolation algorithm.

Using open-source datasets [35], a comparison of nearest-neighbour interpolation and Kriging interpolation is performed. In the Kriging interpolation method, the images are first up-sampled by 200% and then down-sampled by 50%. As shown in Table 2, using peak signal-to-noise ratio (PSNR) as the evaluation metric, Kriging interpolation outperforms nearest-neighbour interpolation in terms of reflecting texture details. Specifically, the super-resolution method based on fractional-order differential Kriging interpolation improves the PSNR by 3–18 dB, a significant enhancement. The super-resolved depth images are converted into point cloud maps for comparative analysis, and it can be seen that most of the edges are correctly recovered. These experimental results show that the fractional-order differential Kriging interpolation algorithm performs very well in the super-resolution of depth images and has significance for further research and application.

The super-resolution results are shown in Figure 11 and Table 2, it can be seen that the edge blurring problem can be solved by edge detection, but the adoption of the simple nearest-neighbour interpolation scheme also causes unwanted jagged edge effects. Subsequently, the super-resolution model based on Kriging interpolation is verified for larger-scale magnification factors using 25% down-sampling and 400% up-sampling. The data in

Table 3. Super-resolution results of the method based on fractional-order differential Kriging interpolation (400% up-sampling) and the nearest-neighbour interpolation method.

Datasets [35] Number	Fractional-Order Differential Kriging Interpolation		Nearest-Neighbour Interpolation
	$\mathit{v}$	PSNR	PSNR
00033	0.6	28.828	25.385
03236	0.6	32.290	23.383
03528	0.6	33.677	25.825
04797	0.6	40.595	22.456
05989	0.6	34.459	25.503
08343	0.6	33.819	21.051

show that the super-resolution method based on Kriging interpolation via fractional-order differential edge detection achieves a greater improvement in super-resolution than that of nearest-neighbour interpolation.

The super-resolution results are shown in Figure 12, Figure 13, Figure 14, Figure 15, Figure 16, Figure 17 and Figure 18 and Table 3. The results of Figure 12 are compared in side view in Figure 13 as an example. As can be seen from Figure 13, the results of the interpolation methods differ in the details. The edge detection method using fractional-order differentiation is implemented to solve the problem of blurred edge points. In terms of texture details, the Kriging interpolation yields more details than the nearest-neighbour interpolation. In particular, the Kriging interpolation yields results with more texture details and a certain level of randomness and fractal characteristics, indicating its superiority over inverse distance interpolation in terms of detail representation.

Figure 14, Figure 15, Figure 16, Figure 17 and Figure 18 show the super-resolution results after 25% down-sampling and 400% up-sampling, and it can be seen that the depth images after super-resolution meet the requirements for effective reconstruction and recognition. Moreover, for complex and variable data, such as the multiple bouquets in Figure 18 Kriging interpolation considerably outperforms nearest-neighbour interpolation.

5. Discussion

This study has shown that most of the classical image interpolation algorithms are not suitable for super-resolution of complex depth images with voids and rapid edge changes. Through comparisons of the point clouds of interpolated depth images, it is found that only nearest-neighbour interpolation of the classical image interpolation algorithms is applicable to depth images. The Kriging interpolation model, especially the variant based on spatial distance, can obtain better results that reflect the spatial characteristics to a greater extent. The introduction of fractional-order differential edge recognition successfully addresses the edge blurring problem of Kriging interpolation. The experimental results demonstrate the superiority of the improved Kriging interpolation based on fractional-order differentiation in terms of the super-resolution of depth images, which preserve texture information and edge properties. Objective evaluations show that the Kriging interpolation with fractional-order differential edge detection significantly improves the interpolation effect compared with the classical interpolation algorithm, with a PSNR improvement of 3–18 dB. Furthermore, the images produced by Kriging interpolation have more fractal characteristics, as would be expected for super-resolution based on interpolation. Overall, the Kriging interpolation super-resolution method based on fractional-order differential edge detection has better image reconstruction ability and stronger robustness than the classical interpolation algorithms.

While our proposed method shows promising results, there are still some limitations that should be addressed in future work. For instance, our current approach assumes that the depth images are noise-free, which may not be realistic in practise. Future work could explore the integration of denoising techniques with our approach to improve its robustness. Additionally, while we have demonstrated the effectiveness of our approach on synthetic depth images, its performance on real-world data should be further evaluated.

Overall, our work provides a valuable contribution to the field of depth image super-resolution and has the potential to be applied in a range of practical applications, such as 3D reconstruction, augmented reality, and medical imaging.