Introduction

In the process of digital registration, the image is recorded into a two-dimensional digital array using elements of a rectangular photomatrix. In this case, whatever technology is used, the spatial resolution of the resulting image is determined by the size of a single photocell of the matrix. And while there has been tremendous progress in high-resolution digital cameras, there is still a need to improve the resolution of digital images. At the same time, a certain limit has already been reached in improving the quality of images by purely technological methods, since with a decrease in the size of elements of photo matrices, the cost of their production increases significantly [1,2].

That is why purely digital methods of image processing to improve their spatial resolution are of great interest. Such a need arises, in particular, in a situation where it is not possible to obtain a high-resolution image, but it is possible to obtain an excessive number of low-quality images, which must then be processed in an optimal way.

The results of this work are important for solving applied problems of automatic scanning, in most problems of digital image processing, in areas related to stereo reconstruction of optical and thermal images, in the development of modern night vision devises, for processing surveillance data from video cameras, and in other fields of science and technology.

And although theoretical developments in this area have existed for several decades, only the rapid progress in the development of modern computer technology has made it possible to successfully apply computational algorithms for solving problems of improving the quality of digital images (their use was previously complicated by the enormous computational costs that computers could not provide several decades ago).

SDM-2021: All-Russian conference, August 24-27, 2021, Novosibirsk, Russia reznik@iae.nsk.su (A. L. Reznik); soloviev@iae.nsk.su (A. A. Soloviev); torgov@iae.nsk.su (A. V. Torgov)

It should also be noted that the choice of optimal processing algorithms (with their help the quality of processed digital images will be improved), in each specific case, depends on two components: first, on a priori knowledge of the statistical characteristics of images, and secondly, on the required specific output image parameters that must be achieved as a result of applying the developed processing algorithm. This explains the presence of a large number of studies in the field of creating mathematical and software-algorithmic methods for improving the quality of matrix images. Each of these studies has a specific purpose and is being developed to address one of the many challenges in digital imaging.

Improving the resolution of digital images based on the calculation of the image with the least variance

As we mentioned above, the possibility of increasing the resolution of images requires reducing the size of the elements of the photomatrix. However, there are a number of technical difficulties here. A possible way out of this situation is multiple acquisition of the same image by a "badresolution" photo sensor, which changes its position during the shooting. In this case, using algorithms that effectively process the results of such a subpixel scan, it is possible to obtain an image with a higher spatial resolution.

Improving the spatial resolution of images -one-dimensional case

Let us consider a one-dimensional registration scheme (see Figure 1). Here the number 𝑁 corresponding to the dimension of the restored vector 𝑋 = (𝑥 1 , 𝑥 2 , . . . , 𝑥 𝑁 ) is a multiple of 𝑙 -the number of resolution elements located into the integrating aperture, that is, the field size (in one-dimensional in the case, the size of the interval) of scanning is an integer number of times larger than the size of the aperture: 𝑁 = 𝑛 × 𝑙.

The described registration mode leads to an underdetermined problem, when the observed data are insufficient for accurate reconstruction of the signal (image), and its statistically valid estimate must be constructed. Various authors [3,4] from the field of signal and image processing have developed mathematical models, computational schemes and algorithms for solving such problems, which allow constructing rather effective schemes for solving specific applied problems. For example, in [5], such an estimate is obtained by algebraic methods by means of pseudo-inversion of matrices. The resulting solution has a number of advantages, but also has certain disadvantages. One of them is that the solution obtained with the help of pseudo-inversion is generally unbalanced (the sum of the weight coefficients with which the samples-observations are included in the solution is different for different elements of the generated image signal), and this leads to an increase in the variance the restored field, which is not always acceptable.

In contrast to the classical approach, which leads to a solution with the minimum norm, in this work we are looking for a solution with the minimum variance (energy), which is not the same in the general case. The solution with the minimum norm is the "least bright" signal corresponding to the system of observations, while the solution with the minimum variance selects the "smoothest" solution from all images that satisfy the observation system.

Let us write the system of equations corresponding to the observation vector 𝐼 = (𝐼 1 , 𝐼 2 , . . . , 𝐼 (𝑛−1)𝑙+1 ):

⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ 𝑥 1 + 𝑥 2 + • • • + 𝑥 𝑙 = 𝐼 1 , 𝑥 2 + 𝑥 3 + • • • + 𝑥 𝑙+1 = 𝐼 2 , . . . 𝑥 (𝑛−1)𝑙+1 + 𝑥 (𝑛−1)𝑙+2 + • • • + 𝑥 𝑛𝑙 = 𝐼 (𝑛−1)𝑙+1 . (1)

It is easy to see that in this case the mean value of the signal

⟨𝑥⟩ = 𝑥 1 + 𝑥 2 + • • • + 𝑛 𝑁 𝑁 = = (𝑥 1 + • • • + 𝑥 𝑙 ) + (𝑥 𝑙+1 + • • • + 𝑥 2𝑙 ) + (𝑥 (𝑛−1)𝑙+1 + • • • + 𝑥 𝑛𝑙 ) 𝑁 = = 𝐼 1 + 𝐼 𝑙+1 + 𝐼 2𝑙+1 + • • • + 𝐼 (𝑛−1)𝑙+1 𝑁(2)

is a constant expressed in terms of the elements of the observation vector 𝐼 and independent of the variables 𝑥 𝑖 . Let's write the expression for the variance

𝐷 𝑥 = 1 𝑁 − 1 𝑁 ∑︁ 𝑖=1 (𝑥 𝑖 − ⟨𝑥⟩) 2 = 1 𝑁 − 1 𝑁 ∑︁ 𝑖=1 (︀ 𝑥 2 𝑖 − 2⟨𝑥⟩𝑥 𝑖 + ⟨𝑥⟩ 2 )︀ = = 1 𝑁 − 1 (︃ 𝑁 ∑︁ 𝑖=1 𝑥 2 𝑖 )︃ − 1 𝑁 − 1 (𝐼 1 + 𝐼 𝑙+1 + 𝐼 2𝑙+1 + • • • + 𝐼 (𝑛−1)𝑙+1 ) 2 𝑁 = = 1 𝑁 − 1 ||𝑥|| 2 − const ⇒ min .(3)

To find this solution, we will do the following. Let us fix free variables 𝑥 1 , 𝑥 2 , . . . , 𝑥 𝑙−1 and express through them the remaining variables 𝑥 𝑙 , 𝑥 𝑙+1 , . . . , 𝑥 𝑁 :

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ 𝑥 𝑙 = 𝐼 1 − 𝑥 1 − • • • − 𝑥 𝑙−1 , 𝑥 𝑙+1 = 𝑥 1 + (𝐼 2 − 𝐼 1 ), 𝑥 𝑙+2 = 𝑥 2 + (𝐼 3 − 𝐼 2 ), . . . 𝑥 2𝑙−1 = 𝑥 𝑙−1 + (𝐼 𝑙 − 𝐼 𝑙−1 ), 𝑥 2𝑙 = 𝐼 1 − 𝑥 1 − • • • − 𝑥 𝑙−1 + (𝐼 𝑙+1 − 𝐼 𝑙 ), 𝑥 2𝑙+1 = 𝑥 1 + (𝐼 2 − 𝐼 1 ) + (𝐼 𝑙+2 − 𝐼 𝑙+1 ), 𝑥 2𝑙+2 = 𝑥 2 + (𝐼 3 − 𝐼 2 ) + (𝐼 𝑙+3 − 𝐼 𝑙+2 ), . . . 𝑥 3𝑙−1 = 𝑥 𝑙−1 + (𝐼 𝑙 − 𝐼 𝑙−1 ) + (𝐼 2𝑙 − 𝐼 2𝑙−1 ), 𝑥 3𝑙 = 𝐼 1 − 𝑥 1 − • • • − 𝑥 𝑙−1 + (𝐼 𝑙+1 − 𝐼 𝑙 ) + (𝐼 2𝑙+1 − 𝐼 2𝑙 ), . . . 𝑥 (𝑛−1)𝑙+1 = 𝑥 1 + (𝐼 2 − 𝐼 1 ) + (𝐼 𝑙+2 − 𝐼 𝑙+1 ) + • • • + (𝐼 (𝑛−2)𝑙+2 − 𝐼 (𝑛−2)𝑙+1 ), 𝑥 (𝑛−1)𝑙+2 = 𝑥 2 + (𝐼 3 − 𝐼 2 ) + (𝐼 𝑙+3 − 𝐼 𝑙+2 ) + • • • + (𝐼 (𝑛−2)𝑙+3 − 𝐼 (𝑛−2)𝑙+2 ), . . . 𝑥 𝑛𝑙−1 = 𝑥 𝑙−1 + (𝐼 𝑙 − 𝐼 𝑙−1 ) + (𝐼 2𝑙 − 𝐼 2𝑙−1 ) + • • • + (𝐼 (𝑛−1)𝑙 − 𝐼 (𝑛−1)𝑙−1 ), 𝑥 𝑛𝑙 = 𝐼 1 −𝑥 1 −. . .−𝑥 𝑙−1 +(𝐼 𝑙+1 − 𝐼 𝑙 )+(𝐼 2𝑙+1 − 𝐼 2𝑙 )+. . .+(𝐼 (𝑛−1)𝑙+1 − 𝐼 (𝑛−1)𝑙 ).(4)

Substituting ( 4) into (2) and equating to zero the partial derivatives of the resulting expression with respect to the variables 𝑥 1 , 𝑥 2 , . . . , 𝑥 𝑙−1 , we obtain a set of (𝑙 − 1) relations

𝑛𝑥 1 + 𝑛𝑥 2 + • • • + 2𝑛𝑥 𝑖 + • • • + 𝑛𝑥 𝑙−1 = = 𝑛𝐼 − [︀(

and, after simple transformations,

𝑛𝑥 𝑖 = 𝑛𝑥 1 −{(𝑛−1) [(𝐼 𝑖+1 −𝐼 𝑖 )−(𝐼 2 −𝐼 1 )]+(𝑛−2) [(𝐼 𝑙+𝑖+1 −𝐼 𝑙+𝑖 ) − (𝐼 𝑙+2 −𝐼 𝑙+1 )]+ +. . .+ 1× [︀ (𝐼 (𝑛−2)𝑙+𝑖+1 −𝐼 (𝑛−2)𝑙+𝑖 )−(𝐼 (𝑛−2)𝑙+2 −𝐼 (𝑛−2)𝑙+1 ) ]︀}︀ , 𝑖 = 2, . . . , 𝑙−1.(6)

Successively applying the last relation (6) to the variables 𝑥 2 , . . . , 𝑥 𝑙−1 and then substituting the obtained expressions into the first of the equations of system (5), we obtain a solution for the element 𝑥 1 :

𝑥 1 = [︂ 2𝑛 − 1 + (𝑙 − 2)(𝑛 − 1) 𝑛𝑙 ]︂ 𝐼 1 + [︂ − 𝑛 − 1 𝑛 ]︂ 𝐼 2 + [︂ 1 𝑛𝑙 ]︂ 𝐼 (𝑛−1)𝑙+1 + + 𝑛−2 ∑︁ 𝑖=1 {︂[︂ (𝑛 − 𝑖 − 1)𝑙 + 1 𝑛𝑙 ]︂ 𝐼 𝑖𝑙+1 + [︂ − 𝑛 − 𝑖 − 1 𝑛 ]︂ 𝐼 𝑖𝑙+2 }︂ .(7)

Further, from ( 6) we find a solution for the elements 𝑥 2 , . . . , 𝑥 𝑙−1 :

𝑥 𝑖 = 𝑥 1 − 1 𝑛 𝑛−2 ∑︁ 𝑠=0 (𝑛 − 𝑠 − 1) [(𝐼 𝑠𝑙+𝑖+1 − 𝐼 𝑠𝑙+𝑖 ) − (𝐼 𝑠𝑙+2 − 𝐼 𝑠𝑙+1 )] , 𝑖 = 2, . . . , 𝑙 − 1. (8)

And, finally, we successively obtain a solution for all other elements of the photomatrix:

𝑥 𝑗 = 𝐼 𝑗−𝑙+1 − 𝑗−1 ∑︁ 𝑞=𝑗−𝑙+1

𝑥 𝑞 , 𝑗 = 𝑙, . . . , 𝑁.

(9)

Improving the spatial resolution of images -two-dimensional case

The problem of reconstructing a two-dimensional field with a minimum dispersion from a set of low-resolution images shifted relative to each other is formulated as follows (Figure 2). We need to construct an estimate for image elements 𝑋 = (𝑥 𝑚𝑛 ), which would satisfy the conditions:

1) 𝑖+𝑙−1 ∑︁ 𝑝=𝑖 𝑗+𝑙−1 ∑︁ 𝑞=𝑗 𝑥 𝑝𝑞 = 𝐼 𝑖𝑗 , 𝑖 = 1, . . . , 𝑀 − 𝑙 + 1, 𝑗 = 1, . . . , 𝑀 − 𝑙 + 1, 2) 𝑀 ∑︁ 𝑝=1 𝑁 ∑︁ 𝑞=1 (𝑥 𝑝𝑞 − ⟨𝑥⟩) 2 ⇒ min, where ⟨𝑥⟩ = 1 𝑀 𝑁 𝑀 ∑︁ 𝑝=1 𝑁 ∑︁ 𝑞=1 𝑥 𝑝𝑞 . (10)

In the discrete case to restore the two-dimensional field 𝑋 = (𝑥 𝑚𝑛 ), 𝑚 = 1, . . . , 𝑀 , 𝑛 = 1, . . . , 𝑁 according to the observed data, which is a matrix 𝐼 = (𝐼 𝑖𝑗 ), 𝑖 = 1, . . . , 𝑀 − 𝑙 + 1, 𝑗 = 1, . . . , 𝑁 − 𝑙 + 1, the solution with the minimum variance is found by analogy with the onedimensional case with the only difference that the variance of the reconstructed field is expressed in terms of its bordering elements (i.e., elements 𝑥 𝑖𝑗 for which 𝑖 ≤ 𝑙 − 1 or 𝑗 ≤ 𝑙 − 1). Thus, the dimension of the problem being solved decreases from 𝑀 × 𝑁 to (𝑙 − 1) × (𝑀 + 𝑁 − 𝑙 + 1), and the filter matrix, which is responsible for the formation of the restored field and does not depend on the observation results, is calculated in advance only once. If the field size is a multiple of the linear size 𝑙 of the integrating aperture, i.e. when 𝑀 = 𝑚 × 𝑙 and 𝑁 = 𝑛 × 𝑙, the solution is found using factorization by two-fold sequential application of the one-dimensional procedure described above. For this, new vectors are formed 𝑦 𝑗 = (𝑦 𝑗1 , 𝑦 𝑗1 , . . . , 𝑦 𝑗𝑚 ), 𝑗 = 1, . . . , 𝑁 −𝑙+1:

𝑦 𝑗𝑘 = 𝑙 ∑︁ 𝑟=1 𝑥 𝑘,𝑗+𝑟−1 , 𝑘 = 1, . . . , 𝑀(11)

and then, for each of 𝑗 = 1, . . . , 𝑁 − 𝑙 + 1, a one-dimensional problem is solved (according to the already described algorithm):

1)

𝑙 ∑︁ 𝑖=1 𝑦 𝑗,𝑘+𝑖−1 = 𝐼 𝑘𝑗 , 𝑘 = 1, . . . , 𝑀 − 𝑙 + 1, 2) 𝑀 ∑︁ 𝑖=1 (𝑦 𝑖𝑗 − ⟨𝑦 𝑗 ⟩) 2 ⇒ min, where ⟨𝑦 𝑗 ⟩ = 1 𝑀 𝑀 ∑︁ 𝑖=1 𝑀 ∑︁ 𝑖=1 𝑦 𝑗𝑖 .(12)

In conclusion, when all the elements of the solution (12) are found by the "columns", i.e. all variables (𝑗 = 1, . . . , 𝑁 − 𝑙 + 1, 𝑘 = 1, . . . , 𝑀 ) are defined, the initial field is restored as a result of the independent solution of 𝑀 one-dimensional problems for each 𝑘 = 1, . . . , 𝑀 :

A great advantage of the developed algorithms is that the filter matrix for specific parameters of restoration can be calculated in advance, which significantly increases the speed of restoration and, in fact (using modern computing equipment), allows the restoration of images in "real time".

Results

The proposed methods have been tested on a large number of real and artificially generated control digital images. In the overwhelming majority of cases, a high quality of recovery has been demonstrated [6,7]. Examples of image reconstruction with minimal variance are shown in Figure 3.

Improving the resolution of video sequences

To restore the spatial resolution of video sequences, the same algorithms are used that are described in Section 2 with the only difference that each of the frames of the video sequence (we use standard video with 24 frames per second) must be processed independently, which leads to a significantly longer calculation time. This time can be reduced by the preliminary calculation of the filter matrix for fixed video parameters, however, the whole video processing process can be quite long for high-resolution video.

When restoring video sequences, the key problem is not the speed of calculations, but the precise positioning of video cameras (it is necessary to achieve the same conditions for subpixel shift, as in the case of processing two-dimensional digital images, as shown in Figure 2).

Simulation on test video sequences showed that with correct positioning, these algorithms can be successfully used for video sequences with obtaining video files of higher spatial resolution. Figure 4 shows examples of image restoration with 9 cameras installed in a 3 × 3 array.

Conclusion

The most important criterion in evaluating the developed algorithms is to determine how the images with the minimum variance correspond to real images. Studies have shown that most real images are quite "smooth", which makes the developed algorithms very effective in improving the spatial resolution of such images. An important distinctive feature of the developed methods is that when performing calculations on special computers, the matrix templates necessary for carrying out all the required mathematical operations can be calculated in advance, which significantly reduces the computational capacity of the algorithms.

The proposed methods have been tested on a large number of real and artificially generated control digital images and video sequences. The high quality of restoration has been demonstrated in the overwhelming majority of cases.