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

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 r × m matrices provides the computation framework of the method. The method has been tested by reconstructing an image represented by an r ×m 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.


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 highresolution 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 r × m 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 T be an r × m real matrix.Equations of the form Tx b, T ∈ R r×m , b ∈ R r occur in many pure and applied problems.It is known that when T is a singular square matrix, then its unique generalized inverse is defined.In the case when T is a real r × m matrix, Penrose showed that there is a unique matrix satisfying the four Penrose equations, called the generalized inverse of T , noted by T † .
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.

Preliminaries and Notation
We shall denote by R r×m the linear space of all r × m real matrices.For T ∈ R r×m , R T will denote the range of T and N T the kernel of T. The generalized inverse T † is the unique matrix that satisfies the following four Penrose equations: where T * denotes the conjugate transpose matrix of T , and T T denotes the transpose matrix of T .
Let us consider the equation Tx b, T ∈ R r×m , b ∈ R r , where T is singular.If b / ∈ R T , then the equation has no solution.Therefore, instead of trying to solve the equation Tx−b 0, we are looking for a minimal norm vector u that minimizes the norm Tu − b .Note that this vector u is unique.So, in this case we consider the equation Tx P R T b, where P R T is the orthogonal projection on R T .Since we are interested in the distance between Tx and b, it is natural to make use of the 2-norm.
The following two propositions can be found in 10 .
Proposition 2.1.Let T ∈ R r×m and b ∈ R r , b / ∈ R T .Then, for u ∈ R m , the following are equivalent: This set of solutions is closed and convex; therefore, it has a unique vector with minimal norm.In the literature cf. 10 , B is known as the set of the least square solutions.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 τ u, v that is used in order to convert an r × r pixel image s x, y from the spatial domain to some other domain in which the image exhibits more readily reducible features is given in the following equation: where u and v are the coordinates in the transform domain and g x, y, u, v denote the general basis function used by the transform.Similarly, the inverse transform is given as where h x, y, u, v represents the inverse of the basis function g x, y, u, v .The two-dimensional version of the function g x, y, u, v 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 x. τ u, y 1 √ r r x 1 s x, y g x, u .

2.4
Now we transform across y τ u, v 1 √ r r y 1 τ u, y g y, u 2.5 and substitute 2.4 : s x, y g x, u g y, u .

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: in which T, S, G, and H are the matrix equivalents of τ, s, g and h 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 S HGSG T H T .

2.8
In order for the transform to be reversible we need H to be the inverse of G and H T to be the inverse of G T , that is, HG G T H T I. Given that G is orthogonal it is trivial to show that this is satisfied when H G T .Given H is merely the transpose of G the inverse function for g x, y, u, v h x, y, u, v is also separable.

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: where x in a , b is the original image, x out a, b represents the measured data from where the original image will be reconstructed and h a, b 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: where i 1, 2, . . ., r and j 1, 2, . . ., m.
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 twodimensional horizontal and vertical convolution, that is, Given a grey value vector x out the digitized degraded image , then x in is the deterministic original image that has to be recovered.The relation between these two components in matrix structure is the following: where H represents a two-dimensional r × m 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., n 30 .The vector x out is of r entries, while the vector x in is of m r n − 1 entries, where m > r and n 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 x in that satisfy the equation Kx in x out , where K 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 r × m matrices.The method retains a restored signal whose norm is smaller than any other solution.

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 nonstationary 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, where x in are the first r elements of the unknown image x in that has to be recovered subject to the constraint Following Proposition 2.2, there is only one solution that minimizes the norm Hx in − x out from an infinite number of exact solutions for x in that satisfy the equation Hx in x out .This unique vector is denoted by x in and can be easily found from: 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 : x in 1 x in 2 x in 3 . . .
x out 1 x out 2 x out 3 . . .
where the index n indicates the linear motion blur in pixels.The elements k 1 , . . ., k n of the matrix are defined as: Similar kernel to that of H, that is, H T can be constructed in order to simulate the vertical motion.The product of HH T 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 n 30.The objective is to calculate the inverse matrix of the blurring kernel HH T and then applied back simple multiplication in the spectral domain to the degraded blurred image x out .That is presented in the following section.

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 s t is given by where ω 2πf is the angular frequency.Since the set of exponentials forms an orthogonal basis the signal can be reconstructed from the projection values: s ω exp iωt dω.

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 x out x in H.

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: x in x out H † .

5.4
The reconstructed image is the inverse Fourier transform of x in .By using our method it not only the advantage of fast recovery but also provides us with an operator H † 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 as rearranging the above equation leads to thus, F m, n can be written in matrix form as where K S y, n * denotes the conjugate transpose of the forward kernel K S y, n .Using the same principles but writing 5.5 in a form where the increasing indexes correspond to higher-frequency coefficients we obtain

5.8
The Fourier coefficients can be seen as the projection coefficients of the image S XY onto a set of complex exponential basis functions that lead to the basis matrix: The approximation of an image S XY in the least square sense can be expressed in terms of the projection matrix P kl : as where T and −1 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 k l 30, k l 100 and k l 450, respectively.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.

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 A, B be bounded operators on H with closed range.Then AB † B † A † if and only if the following three conditions hold.
i The range of AB is closed.
ii A † A commutes with BB * .
iii BB † commutes with A * A.
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 H the underlying Hilbert space, B H the set of all bounded linear operators acting on H and Lat T denotes the set of all closed subspaces of the underlying Hilbert space H invariant under T .

Proposition 5.2. Let A, T ∈ B H be two operators such that A is invertible and T has closed range. Then
5.12 Remark 5.3.In the case when A and B 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 H the corresponding matrix of the picture has dimension less or equal than the rank of the FT matrix W.
For the case of the DFT in matrix form we have that W W T , a symmetric matrix and the inverse DFT matrix is W −1 W * .
The DFT matrix is a r × r matrix, the image has a corresponding matrix H with dimensions r × m, and the generalized inverse of the image H † is a m × r matrix.
In order to satisfy Proposition 5.2, the range of H, R H , must have dimension k less than n rank H k ≤ n .If this condition is satisfied, we have the following theorem, which is an easy consequence of the above discussion: Theorem 5.4.One has

5.13
When k ≤ n, this condition always holds.

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 n of the blurring process.
In most applications noise is unavoidable and a real observation is thus often modeled by Hx in ns x out , 6.1 provided that the noise ns 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 x in k, l h i, j; k, l ns i, j , 6.2 where i 1, 2, . . ., r and j 1, 2, . . ., m.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.
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 : where x in and x out represent the original deterministic image and degraded image respectively, and x 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 k l 10, . . ., 450 .
The second set of tests aimed at accenting the reconstruction error between the original image S XY and the reconstructed image S XY for various values of Fourier coefficients.The calculated quantity is the normalized reconstruction error given by S XY − S XY 2 .6.4 Figure 7 shows the reconstruction error by increasing the moments to k l 450.

Conclusions
In this study, we introduce a fast computational method based on the calculation of the Moore-Penrose inverse of r × m 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.

Proposition 2 . 2 .
Let T ∈ R r×m and b ∈ R r , b / ∈ R T , and the equation Tx b.Then, if T † is the generalized inverse of T , one has that T † b u, where u is the minimal norm solution defined above.

Figure 1 :
Figure 1: Simulation of an a Original image and b Degraded image.

n 30 Figure 2 :
Figure 2: Spectral representations of the degraded image presented in Figure 1 b .

Figure 3 :
Figure 3: Reconstructed images for a k l 30 b k l 100, and c k l 450.

Figure 4 :Figure 5 :
Figure 4: Computational time versus length n of the blurring process in the a horizontal and b bidirectional blurring processes.

Figure 7 :
Figure 7: Reconstruction Error for a bidirectional motion-blurred image.