`Mathematical Problems in EngineeringVolumeÂ 2010Â (2010), Article IDÂ 750352, 14 pageshttp://dx.doi.org/10.1155/2010/750352`
Research Article

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

1Hellenic Transmission System Operator, 18545 Piraeus, Greece
2General Department of Mathematics, Technological Education Institute of Piraeus, 12244 Athens, Greece
3Department of Statistics, Athens University of Economics and Business, 76 Patission Str., 10434 Athens, Greece

Received 17 October 2009; Accepted 19 April 2010

Copyright Â© 2010 Spiros Chountasis et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### 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: 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 , 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: 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 as 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 . Now we transform across and substituteâ€‰â€‰(2.4): 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 , 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 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: where is the original image, 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: where and

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, Given a grey value vector (the digitized degraded image), then is the deterministic original image that has to be recovered. The relation between these two components in matrix structure is the following: 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., ). The vector is of entries, while the vector is of 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 that satisfy the equation , 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, where are the first elements of the unknown image that has to be recovered subject to the constraint Following Proposition 2.2, there is only one solution that minimizes the norm from an infinite number of exact solutions for that satisfy the equation .

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 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): where the index indicates the linear motion blur in pixels. The elements of the matrix are defined as: . 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 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 . 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 where is the angular frequency. Since the set of exponentials forms an orthogonal basis the signal can be reconstructed from the projection values: 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

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: The reconstructed image is the inverse Fourier transform of . 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 as rearranging the above equation leads to thus, can be written in matrix form as 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

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: The approximation of an image in the least square sense can be expressed in terms of the projection matrix : as 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 , and , respectively.

Figure 3: Reconstructed images for (a) (b) , and (c) .

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

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 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.

Figure 4: Computational time versus length 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 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 where and .

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 . 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.

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]: where and represent the original deterministic image and degraded image respectively, and 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 .

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 Figure 7 shows the reconstruction error by increasing the moments to .

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.

#### References

1. K. R. Castleman, Digital Processing, Prentice-Hall, Eglewood Cliffs, NJ, USA, 1996.
2. M. El-Sayed Wahed, â€śImage enhancement using second generation wavelet super resolution,â€ť International Journal of Physical Sciences, vol. 2, no. 6, pp. 149â€“158, 2007.
3. H. J. Trussell and S. Fogel, â€śIdentification and restoration of spatially variant motion blurs in sequential images,â€ť IEEE Transactions on Image Processing, vol. 1, pp. 123â€“126, 1992.
4. S. Chountasis, V. N. Katsikis, and D. Pappas, â€śApplications of the Moore-Penrose inverse in digital image restoration,â€ť Mathematical Problems in Engineering, vol. 2009, Article ID 170724, 12 pages, 2009.
5. M. Hillebrand and C. H. Müller, â€śOutlier robust corner-preserving methods for reconstructing noisy images,â€ť The Annals of Statistics, vol. 35, no. 1, pp. 132â€“165, 2007.
6. A. K. Katsaggelos, J. Biemond, R. M. Mersereau, and R. W. Schafer, â€śA general formulation of constrained iterative restoration algorithms,â€ť in Proceedings of the IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP '85), pp. 700â€“703, Tampa, Fla, USA, 1985.
7. R. W. Schafer, R. M. Mersereau, and M. A. Richards, â€śConstrained iterative restoration algorithms,â€ť Proceedings of the IEEE, vol. 69, no. 4, pp. 432â€“450, 1981.
8. A. Ben-Israel and T. N. E. Greville, Generalized Inverses: Theory and Application, Springer, New York, NY, USA, 2nd edition, 2002.
9. V. N. Katsikis and D. Pappas, â€śFast computing of the Moore-Penrose inverse matrix,â€ť Electronic Journal of Linear Algebra, vol. 17, pp. 637â€“650, 2008.
10. C. Groetsch, Generalized Inverses of Linear Operators, Marcel-Dekker, New York, NY, USA, 1973.
11. R. Gonzalez and P. Wintz, Digital Processing, Addison-Wesley, New York, NY, USA, 2nd edition, 1987.
12. S. A. Dudani, K. J. Breeding, and R. B. McGhee, â€śAircraft identification by moment invariants,â€ť IEEE Transactions on Computers, vol. 26, no. 1, pp. 39â€“46, 1977.
13. P. Milanfar, W. C. Karl, and A. S. Willsky, â€śA moment-based variational approach to tomographic reconstruction,â€ť IEEE Transactions on Image Processing, vol. 5, no. 3, pp. 459â€“470, 1996.
14. R. Mukundan, S. H. Ong, and P. A. Lee, â€śImage analysis by Tchebichef moments,â€ť IEEE Transactions on Image Processing, vol. 10, no. 9, pp. 1357â€“1364, 2001.
15. T. B. Nguyen and B. J. Oommen, â€śMoment-preserving piecewise linear approximations of signals and images,â€ť IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 19, no. 1, pp. 84â€“91, 1997.
16. M. R. Teague, â€śImage analysis via the general theory of moments,â€ť Journal of the Optical Society of America, vol. 70, no. 8, pp. 920â€“930, 1980.
17. C.-H. Teh and R. T. Chin, â€śOn image analysis by the methods of moments,â€ť IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 10, no. 4, pp. 496â€“513, 1988.
18. D. C. Schuurman and D. W. Capson, â€śVideo-rate eigenspace methods for position tracking and remote monitoring,â€ť in Proceedings of the 5th IEEE Southwest Symposium on Image Analysis and Interpretation, pp. 45â€“49, 2002.
19. M. Pawlak, â€śOn the reconstruction aspects of moment descriptors,â€ť IEEE Transactions on Information Theory, vol. 38, no. 6, pp. 1698â€“1708, 1992.
20. R. Bouldin, â€śThe pseudo-inverse of a product,â€ť SIAM Journal on Applied Mathematics, vol. 24, no. 4, pp. 489â€“495, 1973.
21. S. Karanasios and D. Pappas, â€śGeneralized inverses and special type operator algebras,â€ť Facta Universitatis (NIS), Series Mathematics and Informatics, no. 21, pp. 41â€“48, 2006.
22. M. R. Banham and A. K. Katsaggelos, â€śDigital image restoration,â€ť IEEE Signal Processing Magazine, vol. 14, pp. 24â€“41, 1997.