# Applications of the Moore-Penrose Inverse in Digital Image Restoration

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

---

## Abstract

This paper presents a fast computational method that finds application
in a broad scientific field such as digital image restoration. The proposed method provides a new approach to the problem of image reconstruction by using the Moore-Penrose inverse. The resolution of the reconstructed image remains at a very high level but the main advantage of the method was found on the computational load that has been decreased considerably
compared to the classic techniques.


---

## Body

## 1. Introduction

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, who was apparently unaware of Moore's work. These two definitions are equivalent (as it was pointed by Rao in 1956) and since then, 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 singular, then its unique generalized inverse T† (known as the Moore-Penrose 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. [1]). 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 method (see theginv function in [2]) in order to estimate the Moore-Penrose inverse matrix. The computational effort required for the ginv function in order to obtain 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.The field of image restoration has seen a tremendous growth in interest over the last two decades [3–8]. The recovery of an original image from degraded observations is of crucial importance and finds application in several scientific areas including medical imaging and diagnosis, military surveillance, satellite and astronomical imaging, and remote sensing. A number of various algorithms have been proposed and intensively studied for achieving a fast-recovered and high-resolution reconstructed images (see, e.g., [9–11]). In this paper the proposed method provide us a fast computational algorithm for the calculation of the Moore-Penrose inverse of full rank r×m matrices in order to have a fast and accurate digital image restoration. This method is compared with the well-known and also commonly used algorithm of Lagrange multipliers.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, also the required computational time provides another criterion of comparison. In this work, it is clear that the proposed method is marginally better for the former, while for the latter criterion, it is substantially less.


## 2. Preliminaries and Notation

We will 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 transpose matrix of T.Let us consider the equationTx=b,T∈ℝr×m,b∈ℝr, where T is singular. If T is an arbitrary matrix, then there may be none, one, or an infinite number of solutions, depending on whether b∈R(T) or not, and on the rank of T. However if b∉R(T), then the equation has no solution. Therefore, another point of view of this problem is the following: 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 ∥T∥2 norm.The following two propositions can be found in [12].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∥,∀x∈ℝm,(iii)
T*Tu=T*b.LetB={u∈ℝm∣T*Tu=T*b}. This set of solutions is closed and convex, therefore, it has a unique vector with minimal norm. Note that in literature (cf. [12]), B is known as the set of the least square solutions.Proposition 2.2.
LetT∈ℝr×m and b∈ℝr,b∉R(T), and 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 will make use of this property (i.e., ifT† is the generalized inverse of T then T†b=u, where u is the minimal norm solution) for the construction of an alternative method in image processing inverse problems.The algorithm used in this paper for the fast computing of the Moore-Penrose inverse is based on the following theorem.Theorem 2.3 (see [13, Theorem 3.2]).
Letℋ be a Hilbert space. If T=∑i=1nei⊗fi is a rank-n operator then its generalized inverse is also a rank-n operator and for each x∈ℋ, it is defined by the relation
(2.2)T†x=∑i=1nλi(x)ei,
where the functions λi are the solution of an appropriately defined n×n linear system.The proposed algorithm that implements the ideas of the above theorem, can be found in [2], and is calculating the corresponding λi's.In the scientific area of image processing the analytical form of a linear degraded image is given by the following integral equation:(2.3)xout(i,j)=∬Dxin(u,v)h(i,j;u,v)dudv,
where xin(u,v) is the original image, xout(i,j) represents the measured data from where the original image will be reconstructed, and h(i,j;u,v) is a known Point Spread Function (PSF). The PSF depends on the measurement imaging system and is defined as the output of the system for an input point source.A digital image described in a two-dimensional discrete space is derived from an analogue imagexin(u,v) in a two-dimensional continuous space through a sampling process that is frequently referred to as digitization. The two-dimensional continuous image is divided into r rows and m columns. The intersection of a row and a column is termed a pixel. The discrete model for the above linear degradation of an image can be formed by the following summation:(2.4)xout(i,j)=∑k=1r∑l=1mxin(u,v)h(i,j;u,v),
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 [14]. Therefore, the analytically expression for the degraded system is given by a two-dimensional (horizontal and vertical) convolution, that is,(2.5)xout(i,j)=∑k=1r∑l=1mxin(u,v)h(i-u,j-v)=xin(i,j)**h(i,j),
where ** indicates two-dimensional convolution.In the formulation of (2.4), the noise can also be simulated by rewriting the equation as(2.6)xout(i,j)=∑k=1r∑l=1mxin(u,v)h(i,j;u,v)+n(i,j)=xin(i,j)**h(i,j)+n(i,j),
where n(i,j) is an additive noise introduced by the system.However, in this paper the noise is image related which means that the noise has been added to the image.


## 3. Restoration of a Blurry Image

This work introduces a new technique for the removal of blur in an image caused by the uniform linear motion. The method assumes that the linear motion corresponds to a discrete number of pixels and is aligned with the horizontal or vertical sampling.Givenxout (the digitized degraded X-ray 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.1)Hxin=xout,
where H represents a two-dimensional (r×m) priori knowledge matrix or it can be estimated from the degraded X-ray image using its Fourier spectrum [15]. 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.1).However, since there is an infinite number of exact solutions forxin that satisfy the equation Hxin=xout, an additional criterion that finds a sharp restored vector is required. Our work provides a new criterion for restoration of a blurred image including a fast computational method in order to calculate the Moore-Penrose inverse of full rank r×m matrices. The method retains a restored signal whose norm is smaller than any other solution. The computational load for the method is compared with the already known methods.The criterion for restoration of a blurred image that we are using is the minimum distance of the measured data, that is,(3.2)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 ∥Hxin-xout∥=0. In fact, zero is not always attained, but following Proposition 2.1(ii), the norm is minimized.The corresponding optimization problem can be written, alternatively, by minimizing the Lagrange multipliers expression:(3.3)W(xin)=∥Pxin-xout∥2+λ∥Hxin-xout∥2,
where λ (the Lagrange multiplier) must be large enough. The solution to this problem is the following:(3.4)x̃in=[λHTH+PTP]-1[λH+P]Txout,
where P is a projection r×m matrix, projecting on the vector xout:(3.5)P=[Ir|0⋯00⋯0⋮⋮⋮0⋯0],
and Ir is the r×r identity matrix.In general, the PSF varies independently with respect to both (horizontal and vertical) directions, because the degradation of a PSF may depend on its location in the image. An example of this kind of behavior is an optical system that suffers strong geometric aberrations. However, in most of the studies, the PSF is accurately written as a function of the horizontal and vertical displacements independently of the location within the field of view.


## 4. The Generalized Inverse Approach

A blurred image, in our case an X-ray 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(4.1)[k1⋯kn00000k1⋯kn00000k1⋯kn00⋮⋮⋮⋮⋮⋮⋮000⋯k1⋯kn]·[xin_1xin_2xin_3⋮xin_m]=[xout_1xout_2xout_3⋮xout_r],
where the index n indicates the linear motion blur in pixels. The element k1,…,kn of the matrix are defined as kl=1/n(1≤l≤n).However, (4.1) can also be written in the pointwise form for i=1,…,r,(4.2)xout(i)=1n∑h=0n-1xin(i+h),
which describes an underdetermined system of r simultaneous equations and m=r+n-1 unknowns. The objective is to calculate the original column per column data of the image.For this reason, given each column[xout_1,xout_2,xout_3,…,xout_r]T of a degraded blurred image xout, (4.1) results the corresponding column [xin_1,xin_2,xin_3,…,xin_m]T of the original image.As we have seen, the matrixH is an r×m matrix, and the rank of H is less or equal to m. Therefore, the linear system of equations is underdetermined. The proper generalized inverse for this case is a left inverse, which is also called a {1,2,4} inverse, in the sense that it needs to satisfy only the three of the four Penrose equations. A left inverse gives the minimum norm solution of this underdetermined linear system, for every xout∈ℛ(H). The Moore-Penrose Inverse is clearly suitable for our case, since we can have a minimum norm solution for every xout∈ℛ(H), and a minimal norm least squares solution for every xout∉ℛ(H).The proposed algorithm has been tested on a simulated blurred image produced by applying the matrixH on the original image. This can be represented as(4.3)xout(i,j)=1n∑h=0n-1xin(i,j+h),
where i=1,…,r,j=1,…,m for m=r+n-1, and n is the linear motion blur in pixels.Following the above, and the analysis given in Section3, there is an infinite number of exact solutions for xin that satisfy Hxin=xout, but from Proposition 2.2, only one of them minimizes the norm ∥Hxin-xout∥.We will denote this unique vector byx̃in. So, x̃in can be easily found from the equation:(4.4)x̃in=H†xout.The following section presents comparative results that highlight the performance of the generalized inverse.


## 5. Experimental Results

In this section we apply the proposed method on an X-ray image and present numerical results that compare the proposed method with the Lagrangian.The numerical tasks have been performed using Matlab programming language. Specifically, the Matlab 7.4 (R2007b) environment was used on an Intel(R) Pentium(R) Dual CPU T2310 @ 1.46 GHz 1.47 GHz 32-bit system with 2 GB of RAM memory running on the Windows Vista Home Premium Operating System.


### 5.1. Applications on a Degraded Image

The X-ray photography provides a crucial method of diagnostic by using the image analysis. Figure1, Original Image, shows such a deterministic original X-ray image. The image has been divided into r = 900 rows and m = 1200 columns. In Figure 1, Degraded Image, we present the degraded X-ray image for n=30. Finally, from Figure 1, Lagrange Reconstructed Image, and Generalized Inverse Reconstructed Image, it is clearly seen that the details of the original image have been recovered. These figures demonstrate two different methods of reconstruction, the Lagrange multipliers and the generalized inverse, respectively.Restoration of a simulated degraded image for length of the blurring processn=30.
(a)(b)(c)(d)Clearly, the difference between the two methods is insignificant, and can hardly be seen by the human eye. For this reason, the improvement in signal-to-noise ratio (ISNR) has been chosen in order to compare the reconstructed images obtained by the proposed algorithm and the Lagrange multiplier algorithm.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:(5.1)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 2 shows the corresponding ISNR and computational time values.ISNR and computational time versus number of pixels in the blurring process (n=1,…,60).
(a)(b)As far as the ISNR concerns, the difference is small, of the scale of 1 dB. The required computational time provides another criterion of comparison for a range ofn. It is clear that the computational load is less by using the Moore-Penrose inverse matrix.In Table1, we present the numerical results for the calculation of the ISNR as well as the computational time for various random matrices with dimensions 1000×(1000+a), where a=5,10,…,50. It is clear that although the generalized inverse reconstructed method provides better results for the ISNR case the main advantage of the proposed method is its time efficiency.Table 1
ISNR and computational time results for 10 random matrices.
aGinvLagrangeGinvLagrangeISNRISNRcomputation timecomputation time50.35340.35876.42109.2040100.34850.36357.078010.0970150.34840.36188.794010.6780200.34750.35689.199011.6580250.34570.36989.776012.0000300.35370.364310.181012.5540350.35460.365110.771012.5320400.35240.362311.123013.1120450.36420.366011.899014.3230500.35590.377812.140014.9670


### 5.2. Applications on a Degraded and Noisy Image

Noise may be introduced into an image in a number of different ways. In (2.6) the noise has been introduced in an additive way. Here, 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 3, degraded image, we present the original image while a motion blurred and a salt and pepper noise has been added to it.Restoration of a simulated degraded (n=40) and noisy (salt and pepper) image.
(a)(b)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 modifies a pixel value, given the original value of the pixel and the values that surrounding it. Figure3, Generalized Inverse Reconstructed Inverse, demonstrates the result of applying a low-pass rotationally symmetric Gaussian filter of standard deviation equal to 45. Accordingly, Figure 4 provides a graphical representation for the ISNR and the computational time, of the reconstructed and filtered image for different values of n.ISNR and computational time of a degraded and noisy image versus number of pixels in the blurring process.
(a)(b)In order to justify the time efficiency of the proposed method in Table2 we present, for the computation time, the deviation as well as the slope in the cases of degraded image and degraded and noisy image.Table 2
Deviation and slope.
DeviationDeviationSlopeSlopedegradeddegraded and noisydegradeddegraded and noisyLagrange0.82080.1090.04760.0015generalized inverse0.43210.12340.0251−0.0029The second set of tests aimed at accenting the reconstruction error between the original imagexin and the reconstructed image x̃in for various values of linear motion blur, n. The calculated quantity is the normalized reconstruction error given by(5.2)E=1∑i=1r∑j=1m[xin(i,j)]2∑i=1r∑j=1m[xin(i,j)-x̃in(i,j)]2,
using the generalized inverse reconstructed method.Figure5 shows the reconstruction error by increasing the number of pixels in the blurring process n=1,…,60. The solid line represents the degraded image case and the dotted line the noisy image.Figure 5
Reconstruction error of a degraded and noisy image versus number of pixels in the blurring process.


## 5.1. Applications on a Degraded Image

The X-ray photography provides a crucial method of diagnostic by using the image analysis. Figure1, Original Image, shows such a deterministic original X-ray image. The image has been divided into r = 900 rows and m = 1200 columns. In Figure 1, Degraded Image, we present the degraded X-ray image for n=30. Finally, from Figure 1, Lagrange Reconstructed Image, and Generalized Inverse Reconstructed Image, it is clearly seen that the details of the original image have been recovered. These figures demonstrate two different methods of reconstruction, the Lagrange multipliers and the generalized inverse, respectively.Restoration of a simulated degraded image for length of the blurring processn=30.
(a)(b)(c)(d)Clearly, the difference between the two methods is insignificant, and can hardly be seen by the human eye. For this reason, the improvement in signal-to-noise ratio (ISNR) has been chosen in order to compare the reconstructed images obtained by the proposed algorithm and the Lagrange multiplier algorithm.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:(5.1)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 2 shows the corresponding ISNR and computational time values.ISNR and computational time versus number of pixels in the blurring process (n=1,…,60).
(a)(b)As far as the ISNR concerns, the difference is small, of the scale of 1 dB. The required computational time provides another criterion of comparison for a range ofn. It is clear that the computational load is less by using the Moore-Penrose inverse matrix.In Table1, we present the numerical results for the calculation of the ISNR as well as the computational time for various random matrices with dimensions 1000×(1000+a), where a=5,10,…,50. It is clear that although the generalized inverse reconstructed method provides better results for the ISNR case the main advantage of the proposed method is its time efficiency.Table 1
ISNR and computational time results for 10 random matrices.
aGinvLagrangeGinvLagrangeISNRISNRcomputation timecomputation time50.35340.35876.42109.2040100.34850.36357.078010.0970150.34840.36188.794010.6780200.34750.35689.199011.6580250.34570.36989.776012.0000300.35370.364310.181012.5540350.35460.365110.771012.5320400.35240.362311.123013.1120450.36420.366011.899014.3230500.35590.377812.140014.9670


## 5.2. Applications on a Degraded and Noisy Image

Noise may be introduced into an image in a number of different ways. In (2.6) the noise has been introduced in an additive way. Here, 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 3, degraded image, we present the original image while a motion blurred and a salt and pepper noise has been added to it.Restoration of a simulated degraded (n=40) and noisy (salt and pepper) image.
(a)(b)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 modifies a pixel value, given the original value of the pixel and the values that surrounding it. Figure3, Generalized Inverse Reconstructed Inverse, demonstrates the result of applying a low-pass rotationally symmetric Gaussian filter of standard deviation equal to 45. Accordingly, Figure 4 provides a graphical representation for the ISNR and the computational time, of the reconstructed and filtered image for different values of n.ISNR and computational time of a degraded and noisy image versus number of pixels in the blurring process.
(a)(b)In order to justify the time efficiency of the proposed method in Table2 we present, for the computation time, the deviation as well as the slope in the cases of degraded image and degraded and noisy image.Table 2
Deviation and slope.
DeviationDeviationSlopeSlopedegradeddegraded and noisydegradeddegraded and noisyLagrange0.82080.1090.04760.0015generalized inverse0.43210.12340.0251−0.0029The second set of tests aimed at accenting the reconstruction error between the original imagexin and the reconstructed image x̃in for various values of linear motion blur, n. The calculated quantity is the normalized reconstruction error given by(5.2)E=1∑i=1r∑j=1m[xin(i,j)]2∑i=1r∑j=1m[xin(i,j)-x̃in(i,j)]2,
using the generalized inverse reconstructed method.Figure5 shows the reconstruction error by increasing the number of pixels in the blurring process n=1,…,60. The solid line represents the degraded image case and the dotted line the noisy image.Figure 5
Reconstruction error of a degraded and noisy image versus number of pixels in the blurring process.


## 6. Conclusions

In this study, we introduce a novel computational method based on the calculation of the Moore-Penrose inverse of full rankr×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, we present the results by comparing our method and that of the Lagrange multipliers algorithm, a well-established and intensively studied algorithm used for fast-recovered and high-resolution reconstructed images. The efficiency of the generalized inverse is evidenced by the presented simulation results. 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.



---

*Source: 170724-2009-11-05.xml*