Digital Image Reconstruction in the Spectral Domain Utilizing the Moore-Penrose Inverse

Abstract

The field of image restoration has seen a tremendous growth in interest over the last two decades. The recovery of an original image from degraded observations is a crucial method and finds application in several scientific areas including medical imaging and diagnosis, military surveillance, satellite and astronomical imaging, and remote sensing. The proposed approach presented in this work employs Fourier coefficients for moment-based image analysis. The main contributions of the presented technique, are that the image is first analyzed in orthogonal basis matrix formulation increasing the selectivity on image components, and then transmitted in the spectral domain. After the transmission has taken place, at the receiving end the image is transformed back and reconstructed from a set of its geometrical moments. The calculation of the Moore-Penrose inverse of 𝑟×𝑚 matrices provides the computation framework of the method. The method has been tested by reconstructing an image represented by an 𝑟×𝑚 matrix after the removal of blur caused by uniform linear motions. The noise during the transmission process is another issue that is considered in the current work.

1. Introduction

To recover a sharp image from its blurry observation is the problem known as image deblurring. It frequently arises in imaging sciences and technologies and is crucial for allowing to detect important features and patterns [1–3]. A number of various algorithms have been proposed and intensively studied for achieving a fast recovered and high-resolution reconstructed images [4–7].

The main objective of this article is the development of such an algorithm that allows us to restore a blurred image, based on a new fast computational method that calculates the Moore-Penrose inverse of all 𝑟×𝑚 matrices. The proposed method has been applied on the spectral domain of the image. Emphasis has been put on the spectral domain due to its properties for fast recovery and transmission of the image. When approximating an image, an often used method is to divide it into blocks. The proposed approach is based on the consideration that an image is an element of a vector space, which can be expressed as a linear combination of the elements of any orthogonal basis of this space. Accordingly, the image is first analyzed in orthogonal basis matrix formulation and then reconstructed from a set of its geometrical moments. The advantage of representing and recovered any image by choosing a few coefficients is the faster transmission of the image as well as the increased robustness of the method when the image is subject to various attacks that can be introduced during the transmission of the data, including additive noise, cropping and compression.

In this study, the transform that has been chosen is the Fourier transform (FT) because of its ability to transform the intensity information of an image into frequency-related information. That transformation has the advantage of describing global features of an image by the low-frequency components of the transform. Moreover, it is possible to selectively ignore certain components with minimal effect on the final image.

In this work, the generalized inverse of a singular square matrix and of a rectangular matrix plays a crucial role.

The notion of the generalized inverse of a (square or rectangular) matrix was first introduced by Moore in 1920, and again by Penrose in 1955. These two definitions are equivalent, and the generalized inverse of a matrix is also called the Moore-Penrose inverse. Let 𝑇 be an 𝑟×𝑚 real matrix. Equations of the form 𝑇𝑥=𝑏,  𝑇∈ℝ^𝑟×𝑚,  𝑏∈ℝ^𝑟 occur in many pure and applied problems. It is known that when 𝑇 is a singular square matrix, then its unique generalized inverse is defined. In the case when 𝑇 is a real 𝑟×𝑚 matrix, Penrose showed that there is a unique matrix satisfying the four Penrose equations, called the generalized inverse of 𝑇, noted by 𝑇^†.

There are several methods for computing the Moore-Penrose inverse matrix (cf. [8]). One of the most commonly used methods is the Singular Value Decomposition (SVD) method. This method is very accurate but also time-intensive since it requires a large amount of computational resources, especially in the case of large matrices.

In the present work, we apply a very fast and reliable algorithm, in order to estimate the Moore-Penrose inverse matrix, presented in [9], and applied in [4]. The computational effort required for the computation of the generalized inverse is substantially lower, particularly for large matrices, compared to those provided by the SVD method. In addition, we obtain reliable and very accurate approximations.

2. Preliminaries and Notation

We shall denote by ℝ^𝑟×𝑚 the linear space of all 𝑟×𝑚 real matrices. For 𝑇∈ℝ^𝑟×𝑚, 𝑅(𝑇) will denote the range of 𝑇 and 𝑁(𝑇) the kernel of 𝑇. The generalized inverse 𝑇^† is the unique matrix that satisfies the following four Penrose equations: 𝑇𝑇^†=(𝑇𝑇^†)^∗,𝑇^†𝑇=(𝑇^†𝑇)^∗,𝑇𝑇^†𝑇=𝑇,𝑇^†𝑇𝑇^†=𝑇^†, (2.1) where 𝑇^∗ denotes the conjugate transpose matrix of 𝑇, and 𝑇^𝑇 denotes the transpose matrix of 𝑇.

Let us consider the equation 𝑇𝑥=𝑏,  𝑇∈ℝ^𝑟×𝑚,  𝑏∈ℝ^𝑟, where 𝑇 is singular. If 𝑏∉𝑅(𝑇), then the equation has no solution. Therefore, instead of trying to solve the equation ‖𝑇𝑥−𝑏‖=0, we are looking for a minimal norm vector 𝑢 that minimizes the norm ‖𝑇𝑢−𝑏‖. Note that this vector 𝑢 is unique. So, in this case we consider the equation 𝑇𝑥=𝑃_𝑅(𝑇)𝑏, where 𝑃_𝑅(𝑇) is the orthogonal projection on ℛ(𝑇). Since we are interested in the distance between 𝑇𝑥 and 𝑏, it is natural to make use of the 2-norm.

The following two propositions can be found in [10].

Proposition 2.1. Let 𝑇∈ℝ^𝑟×𝑚 and 𝑏∈ℝ^𝑟,  𝑏∉𝑅(𝑇). Then, for 𝑢∈ℝ^𝑚, the following are equivalent: (i)𝑇𝑢=𝑃_𝑅(𝑇)𝑏, (ii)‖𝑇𝑢−𝑏‖≤‖𝑇𝑥−𝑏‖, for all 𝑥∈ℝ^𝑚, (iii)𝑇^∗𝑇𝑢=𝑇^∗𝑏.

Let ℬ={𝑢∈ℝ^𝑚∣𝑇^∗𝑇𝑢=𝑇^∗𝑏}. This set of solutions is closed and convex; therefore, it has a unique vector with minimal norm. In the literature (cf. [10]), ℬ is known as the set of the least square solutions.

Proposition 2.2. Let 𝑇∈ℝ^𝑟×𝑚 and 𝑏∈ℝ^𝑟, 𝑏∉𝑅(𝑇), and the equation 𝑇𝑥=𝑏. Then, if 𝑇^† is the generalized inverse of 𝑇, one has that 𝑇^†𝑏=𝑢, where 𝑢 is the minimal norm solution defined above.

We shall make use of this property for the construction of an alternative method in image processing inverse problems.

The general pointwise definition of the transform 𝜏(𝑢,𝑣) that is used in order to convert an 𝑟×𝑟 pixel image 𝑠(𝑥,𝑦) from the spatial domain to some other domain in which the image exhibits more readily reducible features is given in the following equation:1𝜏(𝑢,𝑣)=𝑟^𝑟𝑟_𝑥=1𝑦=1𝑠(𝑥,𝑦)𝑔(𝑥,𝑦,𝑢,𝑣), (2.2) where 𝑢 and 𝑣 are the coordinates in the transform domain and 𝑔(𝑥,𝑦,𝑢,𝑣) denote the general basis function used by the transform. Similarly, the inverse transform is given as1𝑠(𝑥,𝑦)=𝑟^𝑟𝑟_𝑢=1𝑣=1𝜏(𝑢,𝑣)ℎ(𝑥,𝑦,𝑢,𝑣), (2.3) where ℎ(𝑥,𝑦,𝑢,𝑣) represents the inverse of the basis function 𝑔(𝑥,𝑦,𝑢,𝑣).

The two-dimensional version of the function 𝑔(𝑥,𝑦,𝑢,𝑣) in (2.2) can typically be derived as a series of one-dimensional functions. Such functions are referred to as being “separable”, we can derive the separable two-dimensional functions as follows: The transform is performed across 𝑥.𝜏^′1(𝑢,𝑦)=√𝑟^𝑟_𝑥=1𝑠(𝑥,𝑦)𝑔(𝑥,𝑢). (2.4) Now we transform across 𝑦1𝜏(𝑢,𝑣)=√𝑟^𝑟_𝑦=1𝜏^′(𝑢,𝑦)𝑔(𝑦,𝑢) (2.5) and substitute  (2.4):1𝜏(𝑢,𝑣)=𝑟^𝑟𝑟_𝑥=1𝑦=1𝑠(𝑥,𝑦)𝑔(𝑥,𝑢)𝑔(𝑦,𝑢). (2.6) We can use an identical approach in order to write (2.2) and its inverse (2.3) in matrix form, using the standard orthonormal basis:𝑇=𝐺𝑆𝐺^𝑇,𝑆=𝐻𝑇𝐻^𝑇, (2.7) in which 𝑇,𝑆,𝐺, and 𝐻 are the matrix equivalents of 𝜏,𝑠,𝑔 and ℎ respectively. This is due to our use of orthogonal basis functions, meaning the basis function is its own inverse. Therefore, it is easy to see that the complete process to perform the transform, and then invert it is thus𝑆=𝐻𝐺𝑆𝐺^𝑇𝐻^𝑇. (2.8) In order for the transform to be reversible we need 𝐻 to be the inverse of 𝐺 and 𝐻^𝑇 to be the inverse of 𝐺^𝑇, that is, 𝐻𝐺=𝐺^𝑇𝐻^𝑇=𝐼.

Given that 𝐺 is orthogonal it is trivial to show that this is satisfied when 𝐻=𝐺^𝑇. Given 𝐻 is merely the transpose of 𝐺 the inverse function for 𝑔(𝑥,𝑦,𝑢,𝑣)ℎ(𝑥,𝑦,𝑢,𝑣) is also separable.

3. A Model for Linear-Degraded Image (Motion Blur)

In the scientific area of image processing the analytically form of a Linear-Degraded image is given by the following integral equation:𝑥_out(𝑎,𝑏)=_𝐷𝑥_in(𝑎^′,𝑏^′)ℎ(𝑎−𝑎^′,𝑏−𝑏^′)𝑑𝑎^′𝑑𝑏^′, (3.1) where 𝑥_in(𝑎^′,𝑏^′) is the original image, 𝑥_out(𝑎,𝑏) represents the measured data from where the original image will be reconstructed and ℎ(𝑎,𝑏) is a known Point Spread Function (PSF) that depends on the measurement system.

The discrete model for the above linear degradation of an image can be formed by the following summation:𝑥_out1(𝑖,𝑗)=√𝑟𝑚^𝑟𝑚_𝑘=1𝑙=1𝑥_in(𝑘,𝑙)ℎ(𝑖,𝑗;𝑘,𝑙), (3.2) where 𝑖=1,2,…,𝑟 and 𝑗=1,2,…,𝑚.

In this paper we adopt the use of a shift invariant model for the blurring process as in [11]. Therefore, the analytically expression for the degraded system is given by a two-dimensional (horizontal and vertical) convolution, that is,𝑥_out1(𝑖,𝑗)=√𝑟𝑚^𝑟𝑚_𝑘=1𝑙=1𝑥_in(𝑘,𝑙)ℎ(𝑖−𝑘,𝑗−𝑙). (3.3) Given a grey value vector 𝑥_out (the digitized degraded image), then 𝑥_in is the deterministic original image that has to be recovered. The relation between these two components in matrix structure is the following:𝑥_out=𝐻𝑥_in, (3.4) where 𝐻 represents a two-dimensional 𝑟×𝑚 a priori knowledge matrix or it can be estimated from the degraded image using its Fourier spectrum. The latter can be seen in Figure 2 where the distance between two successive lines indicates the length of the blurring process in pixels (i.e., 𝑛=30). The vector 𝑥_out is of 𝑟 entries, while the vector 𝑥_in is of 𝑚(=𝑟+𝑛−1) entries, where 𝑚>𝑟 and 𝑛 is the length of the blurring process in pixels. The problem consists of solving the underdetermined system of (3.4). However, since there is an infinite number of exact solutions for 𝑥_in that satisfy the equation 𝐾𝑥_in=𝑥_out, where 𝐾 is any two-dimensional matrix, an additional criterion that finds a sharp restored image is required.

Our work provides a new criterion for restoration of a blurred image including a novel fast computational method in order to calculate the Moore-Penrose inverse of 𝑟×𝑚 matrices. The method retains a restored signal whose norm is smaller than any other solution.

4. Deblurring in the Spatial Domain

Images can be viewed as non-stationary two-dimensional signals with edges, textures, and deterministic objects at different locations. Figure 1(a) shown such an image. Although non-stationary signals are, in general, characterized by their local features rather than their global ones, it is possible to recover images by introducing global constrains on either its spatial or spectral resolution. The criterion for restoration of a blurred image that we are using is the minimum distance of the measured data, that is, min‖𝑥_in−𝑥_out‖, (4.1) where 𝑥_in are the first 𝑟 elements of the unknown image 𝑥_in that has to be recovered subject to the constraint‖𝐻𝑥_in−𝑥_out‖²=0. (4.2) Following Proposition 2.2, there is only one solution that minimizes the norm ‖𝐻𝑥_in−𝑥_out‖ from an infinite number of exact solutions for 𝑥_in that satisfy the equation 𝐻𝑥_in=𝑥_out.

(a)

(b)

Figure 1 

Simulation of an (a) Original image and (b) Degraded image.

Figure 2 

Spectral representations of the degraded image presented in Figure 1(b).

This unique vector is denoted by 𝑥_in and can be easily found from:𝑥_in=𝐻^†𝑥_out. (4.3)

A blurred image that has been degraded by a uniform linear motion in the horizontal direction, usually results of camera panning or fast object motion can be expressed as in (3.4):𝑘₁…𝑘_𝑛00000𝑘₁…𝑘_𝑛00000𝑘₁…𝑘_𝑛00⋮⋮⋮⋮⋮⋮⋮000…𝑘₁…𝑘_𝑛⋅𝑥_{in_1}𝑥_{in_2}𝑥_{in_3}⋮𝑥_{in_𝑀}=𝑥_{out_1}𝑥_{out_2}𝑥_{out_3}⋮𝑥_{out_𝑁}, (4.4) where the index 𝑛 indicates the linear motion blur in pixels. The elements 𝑘₁,…,𝑘_𝑛 of the matrix are defined as: 𝑘_𝑖=1/𝑛(1≤𝑖≤𝑛). Similar kernel to that of 𝐻, that is, 𝐻^𝑇 can be constructed in order to simulate the vertical motion. The product of 𝐻𝐻^𝑇 in the spectral domain produces a bidirectional blur. By simulating and applied that product kernel to the original image we obtained the degraded version that is presented in Figure 1(b), for 𝑛=30. The objective is to calculate the inverse matrix of the blurring kernel 𝐻𝐻^𝑇 and then applied back (simple multiplication in the spectral domain) to the degraded blurred image 𝑥_out. That is presented in the following section.

5. Deblurring in the Spectral Domain: Application of the Moments on Image Reconstruction

In view of the importance of the frequency domain, the FT has become one of the most widely used signal analysis tool across many disciplines of science and engineering. The FT is generated by projecting the signal on to a set of basis functions, each of which is a sinusoid with a unique frequency. The FT of a time signal 𝑠(𝑡) is given by ̃1𝑠(𝜔)=√2𝜋^+∞−∞𝑠(𝑡)exp(−𝑖𝜔𝑡)𝑑𝑡, (5.1) where 𝜔=2𝜋𝑓 is the angular frequency. Since the set of exponentials forms an orthogonal basis the signal can be reconstructed from the projection values: 1𝑠(𝑡)=√2𝜋^+∞−∞𝑠(̃𝜔)exp(𝑖𝜔𝑡)𝑑𝜔. (5.2) Following the property of the FT that the convolution in the spatial domain is translated into simple algebraic product in the spectral domain (3.4) can be written in the form𝑥_out=𝑥_in𝐻. (5.3)

Figure 2 shows the spectral representation of the degraded image obtained using the above equation.

In order to obtain back the original image, (5.3) is solved in the Fourier space: 𝑥_in=𝑥_out𝐻^†. (5.4) The reconstructed image is the inverse Fourier transform of 𝑥_in. By using our method it not only the advantage of fast recovery but also provides us with an operator 𝐻^† that exists even for not full rank nonsquare matrices. In this section the whole process of deblurring and restoring the original image is done in the spectral domain. It provides us the ability of fast recovering and algorithmic simplicity.

Moments are particularly popular due to their compact description, their capability to select differing levels of detail and their known performance attributes (see [3, 12–17]). It is a well-recognised property of moments that they can be used to reconstruct the original function, that is, none of the original image information is lost in the projection of the image on to the moment basis functions, assuming an infinite number of moments are calculated. Another property for the reconstruction of a band-limited image using its moments is that while derivatives give information on the high frequencies of a signal, moments provide information on its low frequencies. It is known that the higher-order moments capture increasingly higher frequencies within a function and in the case of an image the higher frequencies represent the detail of the image. This is also consistent with work on other types of reconstruction, such as eigenanalysis where it has been found that increasing numbers of eigenvectors are required to capture image detail [18] and again exceed the number required for recognition. Describing images with moments instead of other more commonly used image features means that global properties of the image are used rather than local properties. Moments provide information on its low-frequency of an image. Applying the Fourier coefficients a low-pass approximation of the original image is obtained. It is well known that any image can be reconstructed from its moments in the least-squares sense. Discrete orthogonal moments provide a more accurate description of image features by evaluating the moment components directly in the image coordinate space.

The reconstruction of an image from its moments is not necessarily unique. Thus, all possible methods must impose extra constraints in order to make its moments uniquely solve the reconstruction problem.

The most common reconstruction method of an image from some of its moments is based on the least squares approximation of the image using orthogonal polynomials [14, 19]. In this paper the constraint that introduced is related to the bandwidth of the image and provides a more general reconstruction method. We must keep in mind that this constraint is a global, for a local one a joint bilinear distribution such as Wigner or wavelet must be used.

In a discrete Fourier domain the two-dimensional Fourier coefficients are defined as1𝐹(𝑚,𝑛)=√𝑋𝑌^𝑋𝑌_𝑥=1𝑦=1𝑆_𝑋𝑌exp(−2𝜋𝑖((𝑥−1)(𝑚−1)𝑋+(𝑦−1)(𝑛−1)𝑌)), (5.5) rearranging the above equation leads to 1𝐹(𝑚,𝑛)=√𝑋𝑌^𝑋_𝑥=1exp((−2𝜋𝑖(𝑥−1)(𝑚−1)𝑋)^𝑌_𝑦=1𝑆_𝑋𝑌exp(−2𝜋𝑖(𝑦−1)(𝑛−1)𝑌)), (5.6) thus, 𝐹(𝑚,𝑛) can be written in matrix form as 𝐹(𝑚,𝑛)=𝐾_𝑆(𝑥,𝑚)𝑆_𝑋𝑌𝐾_𝑆(𝑦,𝑛)^∗, (5.7) where 𝐾_𝑆(𝑦,𝑛)^∗ denotes the conjugate transpose of the forward kernel 𝐾_𝑆(𝑦,𝑛).

Using the same principles but writing (5.5) in a form where the increasing indexes correspond to higher-frequency coefficients we obtain 1𝐹(𝑚,𝑛)=√𝑋𝑌^𝑋𝑌_𝑥=1𝑦=1𝑆_𝑋𝑌⋅exp[−2𝜋𝑖((𝑥−1)(𝑚−(𝑘−1)/2−1)𝑋+(𝑦−(𝑙−1)/2−1)(𝑛−1)𝑌)]. (5.8)

The Fourier coefficients can be seen as the projection coefficients of the image 𝑆_𝑋𝑌 onto a set of complex exponential basis functions that lead to the basis matrix: 𝐵_𝑘𝑙1(𝑚,𝑛)=√𝑘exp[−2𝜋𝑖(𝑚−1)(𝑛−(𝑙−1)/2−1)𝑘]. (5.9) The approximation of an image 𝑆_𝑋𝑌 in the least square sense can be expressed in terms of the projection matrix 𝑃_𝑘𝑙: 𝑃_𝑘𝑙=(𝐵_𝑋𝑘)^𝑇𝑆_𝑋𝑌𝐵_𝑌𝑙 (5.10) as 𝑆^′_𝑋𝑌=𝐵_𝑋𝑘(𝐵^𝑇_𝑋𝑘𝐵^𝑇_𝑋𝑘)⁻¹𝑃_𝑘𝑙(𝐵^𝑇_𝑌𝑙𝐵_𝑌𝑙)⁻¹𝐵^𝑇_𝑌𝑙=(𝐵_𝑋𝑘)⁻𝑃_𝑘𝑙(𝐵_𝑌𝑙)^†, (5.11) where ()^𝑇 and ()⁻¹ denote the transpose and the inverse of the given matrix. The operations ()⁻ and ()^† stand for the left and right Moore-Penrose pseudoinverses and are unique. Among the multiple inverse solutions it chooses the one with minimum norm. When considering image reconstruction from moments, the number of moments required for accurate reconstruction will be related to the frequencies present within the original image. For a given image size it would appear that there should be a finite limit to the frequencies that are present in the image and for a binary image that limiting frequency will be relatively low. As the higher-order moments approach this frequency the reconstruction will become more accurate. Figures 3(a), 3(b) and 3(c) present the reconstructed image for the cases of 𝑘=𝑙=30, 𝑘=𝑙=100 and 𝑘=𝑙=450, respectively.

(a)

(b)

(c)

Figure 3 

Reconstructed images for (a) $k=l=30$ (b) $k=l=100$, and (c) $k=l=450$.

From the reconstruction point of view the basis matrix is applied to both original image and blurring kernel transforming these into spectral domain. After the inversion of the blurring kernel, its product with the degraded image is applied to inverted basis functions for the reconstruction of the original image.

5.1. The Fourier Transform and the Reverse Order Law

A natural question to ask about the FT and the Moore-Penrose inverse, is weather the reverse order law for their inverses holds or not. In general, the reverse order law for generalized inverses holds under certain conditions. The following proposition is a restatement of a part of R. Bouldin's theorem [20] that holds for operators and matrices.

Proposition 5.1. Let 𝐴,𝐵 be bounded operators on ℋ with closed range. Then (𝐴𝐵)^†=𝐵^†𝐴^† if and only if the following three conditions hold.(i)The range of 𝐴𝐵 is closed.(ii)𝐴^†𝐴 commutes with 𝐵𝐵^∗.(iii)𝐵𝐵^† commutes with 𝐴^∗𝐴.

A corollary of the above theorem is the following proposition that can be found in [21] and we will use in our case.

Let us denote by ℋ the underlying Hilbert space, ℬ(ℋ) the set of all bounded linear operators acting on ℋ and 𝐿𝑎𝑡𝑇 denotes the set of all closed subspaces of the underlying Hilbert space ℋ invariant under 𝑇.

Proposition 5.2. Let 𝐴,𝑇∈ℬ(ℋ) be two operators such that 𝐴 is invertible and 𝑇 has closed range. Then (𝐴𝑇)^†=𝑇^†𝐴⁻¹iffℛ(𝑇)∈𝐿𝑎𝑡(𝐴^∗𝐴). (5.12)

Remark 5.3. In the case when 𝐴 and 𝐵 are matrices, their range is always closed, so the first condition is always satisfied.

As it is shown in this case, equality holds if the space spanned by the orthonormal basis of the range of 𝐻 (the corresponding matrix of the picture) has dimension less or equal than the rank of the FT matrix 𝑊.

For the case of the DFT in matrix form we have that 𝑊=𝑊^𝑇, (a symmetric matrix) and the inverse DFT matrix is 𝑊⁻¹=𝑊^∗.

The DFT matrix is a 𝑟×𝑟 matrix, the image has a corresponding matrix 𝐻 with dimensions 𝑟×𝑚, and the generalized inverse of the image 𝐻^† is a 𝑚×𝑟 matrix.

In order to satisfy Proposition 5.2, the range of 𝐻, ℛ(𝐻), must have dimension 𝑘 less than 𝑛 (rank (𝐻) = 𝑘≤𝑛).

If this condition is satisfied, we have the following theorem, which is an easy consequence of the above discussion:

Theorem 5.4. One has (𝑊𝐻)^†=𝐻^†𝑊⁻¹iffℛ(𝐻)∈𝐿𝑎𝑡(𝑊^∗𝑊)=𝐼_𝑛. (5.13) When 𝑘≤𝑛, this condition always holds.

6. Error-Time Analysis

Several tests were performed in order to evaluate the effectiveness of the method. The first set of tests aimed at the reconstruction time for different motion blur transformations. Figures 4(a) and 4(b) show the computational time for the reconstruction of the original image shown in Figure 1(a). The computational time has been calculated for both horizontal and bidirectional (horizontal and vertical) degradation of the image and were plotted against the length 𝑛 of the blurring process.

(a)

(b)

Figure 4 

Computational time versus length $n$ of the blurring process in the (a) horizontal and (b) bidirectional blurring processes.

In most applications noise is unavoidable and a real observation is thus often modeled by 𝐻𝑥_in+𝑛𝑠=𝑥_out, (6.1) provided that the noise 𝑛𝑠 is additive, although multiplicative noise can be handled similarly. In the formulation of (3.3) the noise can also be simulated by rewriting the equation as𝑥_out(𝑖,𝑗)=^𝑟𝑚_𝑘=1𝑙=1𝑥_in(𝑘,𝑙)ℎ(𝑖,𝑗;𝑘,𝑙)+𝑛𝑠(𝑖,𝑗), (6.2) where 𝑖=1,2,…,𝑟 and 𝑗=1,2,…,𝑚.

However, in this article, we simulate a noise model where a number of pixels are corrupted and randomly take on a value of white and black (salt and pepper noise) with noise density equal to 0.02. The image that we receive from a faulty transmission line can contain this form of corruption. In Figure 5(a) we present the original image while a motion-blurred and a salt and pepper noise has been added to it. Image processing and analysis are based on filtering the content of the images in a certain way. The filtering process is basically an algorithm that modify a pixel value, given the original value of the pixel and the values of the surrounding it. Figure 5(b) demonstrates the result of applying a low-pass Gaussian filter on a generalized inverse reconstructed image.

(a)

(b)

Figure 5 

(a) Motion-blurred and noisy image, (b) Generalized inverse reconstructed and filtered image.

The Improvement in Signal-to-Noise Ratio (ISNR) is a criterion that has been used extensively for the purpose of objectively testing the performance of image processing algorithms [22]:ISNR=10log₁₀{∑_𝑖,𝑗[𝑥_in(𝑖,𝑗)−𝑥_out(𝑖,𝑗)]²∑_𝑖,𝑗[𝑥_in(𝑖,𝑗)−𝑥_in(𝑖,𝑗)]²}, (6.3) where 𝑥_in and 𝑥_out represent the original deterministic image and degraded image respectively, and 𝑥_in is the corresponding restored image. Figure 6 shows the corresponding ISNR values of the proposed generalized inverse reconstruction algorithm for a number of Fourier coefficients (𝑘=𝑙=10,…,450).

Figure 6 

ISNR versus number of moments.

The second set of tests aimed at accenting the reconstruction error between the original image 𝑆_𝑋𝑌 and the reconstructed image 𝑆^′_𝑋𝑌 for various values of Fourier coefficients. The calculated quantity is the normalized reconstruction error given by 1𝐸=∑^𝑋_𝑥=1∑^𝑌_𝑦=1[𝑆_𝑋𝑌]^{2𝑋𝑌𝑥=1𝑦=1}[𝑆_𝑋𝑌−𝑆^′_𝑋𝑌]². (6.4) Figure 7 shows the reconstruction error by increasing the moments to 𝑘=𝑙=450.

Figure 7 

Reconstruction Error for a bidirectional motion-blurred image.

7. Conclusions

In this study, we introduce a fast computational method based on the calculation of the Moore-Penrose inverse of 𝑟×𝑚 matrix, with particular focus on problems arising in image processing. We are motivated by the problem of restoring blurry and noisy images via well-developed mathematical methods and techniques based on the inverse procedures in order to obtain an approximation of the original image. By using the proposed algorithm, the resolution of the reconstructed image remains at a very high level, although the main advantage of the method was found on the computational load that has been decreased considerably compared to the other methods and techniques. In this paper the results presented where fully demonstrated in the spectral domain. Orthogonal moments have demonstrated significant energy compaction properties that are desirable in the field of image processing, especially in feature and object recognition. The advantage of representing and recovered any image by choosing a few Fourier coefficients, is the faster transmission of the image as well as the increased robustness when the image is subject to various attacks that can be introduced during the transmission of the data, including additive noise. The results of this work are well established by simulating data. Besides digital image restoration, our work on generalized inverse matrices may also find applications in other scientific fields where a fast computation of the inverse data is needed.

The proposed method can be used in any kind of matrix, square or rectangular, full rank or not; so the dimensions and the nature of the image do not play any role in this application.