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

Authors: Spiros Chountasis; Vasilios N. Katsikis; Dimitrios Pappas
Journal: Mathematical Problems in Engineering
 (2010)

---

## 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 ofr×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.


---

## Body

## 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 allr×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. LetT be an r×m real matrix. Equations of the form Tx=b,  T∈ℝr×m,  b∈ℝ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.


## 2. Preliminaries and Notation

We shall denote byℝr×m the linear space of all r×m real matrices. For T∈ℝ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: (2.1)TT†=(TT†)*,T†T=(T†T)*,TT†T=T,T†TT†=T†,
where T* denotes the conjugate transpose matrix of T, and TT denotes the transpose matrix of T.Let us consider the equationTx=b,  T∈ℝr×m,  b∈ℝ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=PR(T)b, where PR(T) is the orthogonal projection on ℛ(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.
LetT∈ℝr×m and b∈ℝr,  b∉R(T). Then, for u∈ℝm, the following are equivalent:  (i)
Tu=PR(T)b,(ii)
∥Tu-b∥≤∥Tx-b∥, for all  x∈ℝm,(iii)
T*Tu=T*b.Letℬ={u∈ℝm∣T*Tu=T*b}. 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.
LetT∈ℝr×m and b∈ℝ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.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:(2.2)τ(u,v)=1r∑x=1r∑y=1rs(x,y)g(x,y,u,v),
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(2.3)s(x,y)=1r∑u=1r∑v=1rτ(u,v)h(x,y,u,v),
where h(x,y,u,v) represents the inverse of the basis function g(x,y,u,v).The two-dimensional version of the functiong(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.(2.4)τ′(u,y)=1r∑x=1rs(x,y)g(x,u).
Now we transform across y(2.5)τ(u,v)=1r∑y=1rτ′(u,y)g(y,u)
and substitute  (2.4):(2.6)τ(u,v)=1r∑x=1r∑y=1rs(x,y)g(x,u)g(y,u).
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)T=GSGT,S=HTHT,
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(2.8)S=HGSGTHT.
In order for the transform to be reversible we need H to be the inverse of G and HT to be the inverse of GT, that is, HG=GTHT=I.Given thatG is orthogonal it is trivial to show that this is satisfied when H=GT. 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.


## 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:(3.1)xout(a,b)=∫∫Dxin(a′,b′)h(a-a′,b-b′)da′db′,
where xin(a′,b′) is the original image, xout(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:(3.2)xout(i,j)=1rm∑k=1r∑l=1mxin(k,l)h(i,j;k,l),
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 two-dimensional (horizontal and vertical) convolution, that is,(3.3)xout(i,j)=1rm∑k=1r∑l=1mxin(k,l)h(i-k,j-l).
Given a grey value vector xout (the digitized degraded image), then xin is the deterministic original image that has to be recovered. The relation between these two components in matrix structure is the following:(3.4)xout=Hxin,
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 xout is of r entries, while the vector xin 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 xin that satisfy the equation Kxin=xout, 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 ofr×m 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. Figure1(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, (4.1)min∥x̂in-xout∥,
where x̂in are the first r elements of the unknown image xin that has to be recovered subject to the constraint(4.2)∥Hxin-xout∥2=0.
Following  Proposition 2.2, there is only one solution that minimizes the norm ∥Hxin-xout∥ from an infinite number of exact solutions for xin that satisfy the equation Hxin=xout.Simulation of an (a) Original image and (b) Degraded image.
(a)(b)Figure 2
Spectral representations of the degraded image presented in Figure1(b).This unique vector is denoted byx̂in and can be easily found from:(4.3)xin=H†xout.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):(4.4)[k1…kn00000k1…kn00000k1…kn00⋮⋮⋮⋮⋮⋮⋮000…k1…kn]·[xin_1xin_2xin_3⋮xin_M]=[xout_1xout_2xout_3⋮xout_N],
where the index n indicates the linear motion blur in pixels. The elements k1,…,kn of the matrix are defined as: ki=1/n(1≤i≤n). Similar kernel to that of H, that is, HT can be constructed in order to simulate the vertical motion. The product of HHT 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 HHT and then applied back (simple multiplication in the spectral domain) to the degraded blurred image xout. 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 signals(t) is given by (5.1)s̃(ω)=12π∫-∞+∞s(t)exp(-iωt)dt,
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: (5.2)s(t)=12π∫-∞+∞s̃(ω)exp(iωt)dω.
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(5.3)x̃out=x̃inH̃.Figure2 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: (5.4)x̃in=x̃outH̃†.
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(5.5)F(m,n)=1XY∑x=1X∑y=1YSXYexp(-2πi((x-1)(m-1)X+(y-1)(n-1)Y)),
rearranging the above equation leads to (5.6)F(m,n)=1XY∑x=1Xexp((-2πi(x-1)(m-1)X)∑y=1YSXYexp(-2πi(y-1)(n-1)Y)),
thus, F(m,n) can be written in matrix form as (5.7)F(m,n)=KS(x,m)SXYKS(y,n)*,
where KS(y,n)* denotes the conjugate transpose of the forward kernel KS(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)F(m,n)=1XY∑x=1X∑y=1YSXY·exp[-2πi((x-1)(m-(k-1)/2-1)X+(y-(l-1)/2-1)(n-1)Y)].The Fourier coefficients can be seen as the projection coefficients of the imageSXY onto a set of complex exponential basis functions that lead to the basis matrix: (5.9)Bkl(m,n)=1kexp[-2πi(m-1)(n-(l-1)/2-1)k].
The approximation of an image SXY in the least square sense can be expressed in terms of the projection matrix Pkl: (5.10)Pkl=(BXk)TSXYBYl
as (5.11)SXY′=BXk(BXkTBXkT)-1Pkl(BYlTBYl)-1BYlT=(BXk)-Pkl(BYl)†,
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.Reconstructed images for (a)k=l=30 (b) k=l=100, and (c) k=l=450.
(a)(b)(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.
LetA,B be bounded operators on ℋ with closed range. Then (AB)†=B†A† if and only if the following three conditions hold.(i)
The range ofAB 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ℋ the underlying Hilbert space, ℬ(ℋ) the set of all bounded linear operators acting on ℋ and LatT denotes the set of all closed subspaces of the underlying Hilbert space ℋ invariant under T.Proposition 5.2.
LetA,T∈ℬ(ℋ) be two operators such that A is invertible and T has closed range. Then
(5.12)(AT)†=T†A-1iffℛ(T)∈Lat(A*A).Remark 5.3.
In the case whenA 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 ofH (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 thatW=WT, (a symmetric matrix) and the inverse DFT matrix is W-1=W*.The DFT matrix is ar×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  Proposition5.2,  the range of H, ℛ(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)(WH)†=H†W-1iffℛ(H)∈Lat(W*W)=In.
When k≤n, this condition always holds.


## 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.
LetA,B be bounded operators on ℋ with closed range. Then (AB)†=B†A† if and only if the following three conditions hold.(i)
The range ofAB 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ℋ the underlying Hilbert space, ℬ(ℋ) the set of all bounded linear operators acting on ℋ and LatT denotes the set of all closed subspaces of the underlying Hilbert space ℋ invariant under T.Proposition 5.2.
LetA,T∈ℬ(ℋ) be two operators such that A is invertible and T has closed range. Then
(5.12)(AT)†=T†A-1iffℛ(T)∈Lat(A*A).Remark 5.3.
In the case whenA 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 ofH (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 thatW=WT, (a symmetric matrix) and the inverse DFT matrix is W-1=W*.The DFT matrix is ar×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  Proposition5.2,  the range of H, ℛ(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)(WH)†=H†W-1iffℛ(H)∈Lat(W*W)=In.
When k≤n, 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. Figures4(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.Computational time versus lengthn of the blurring process in the (a) horizontal and (b) bidirectional blurring processes.
(a)(b)In most applications noise is unavoidable and a real observation is thus often modeled by(6.1)Hxin+ns=xout,
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(6.2)xout(i,j)=∑k=1r∑l=1mxin(k,l)h(i,j;k,l)+ns(i,j),
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.(a) Motion-blurred and noisy image, (b) Generalized inverse reconstructed and filtered image.
(a)(b)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]:(6.3)ISNR=10log10{∑i,j[xin(i,j)-xout(i,j)]2∑i,j[xin(i,j)-x̃in(i,j)]2},
where xin and xout  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).Figure 6
ISNR versus number of moments.The second set of tests aimed at accenting the reconstruction error between the original imageSXY and the reconstructed image SXY′ for various values of Fourier coefficients. The calculated quantity is the normalized reconstruction error given by (6.4)E=1∑x=1X∑y=1Y[SXY]2∑x=1X∑y=1Y[SXY-SXY′]2.
Figure 7 shows the reconstruction error by increasing the moments to k=l=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 ofr×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.



---

*Source: 750352-2010-06-23.xml*