Applications of the Moore-Penrose Inverse in Digital Image Restoration

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.

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.

Let 𝑇 be an 𝑟×𝑚 real matrix. Equations of the form 𝑇𝑥=𝑏,𝑇∈ℝ^𝑟×𝑚,𝑏∈ℝ^𝑟 occur in many pure and applied problems. It is known that when 𝑇 is singular, then its unique generalized inverse 𝑇^† (known as the Moore-Penrose 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. [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 the ginv 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 𝑟×𝑚 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 ℝ^𝑟×𝑚 the linear space of all 𝑟×𝑚 real matrices. For 𝑇∈ℝ^𝑟×𝑚, 𝑅(𝑇) will denote the range of 𝑇 and 𝑁(𝑇) the kernel of 𝑇. The generalized inverse 𝑇^† is the unique matrix that satisfies the following four Penrose equations:

𝑇𝑇^†=(𝑇𝑇^†)^∗,𝑇^†𝑇=(𝑇^†𝑇)^∗,𝑇𝑇^†𝑇=𝑇,𝑇^†𝑇𝑇^†=𝑇^†, (2.1) where 𝑇^∗ denotes the transpose matrix of 𝑇.

Let us consider the equation 𝑇𝑥=𝑏,𝑇∈ℝ^𝑟×𝑚,𝑏∈ℝ^𝑟, where 𝑇 is singular. If 𝑇 is an arbitrary matrix, then there may be none, one, or an infinite number of solutions, depending on whether 𝑏∈𝑅(𝑇) or not, and on the rank of 𝑇. However if 𝑏∉𝑅(𝑇), then the equation has no solution. Therefore, another point of view of this problem is the following: instead of trying to solve the equation ‖𝑇𝑥−𝑏‖=0, we are looking for a minimal norm vector 𝑢 that minimizes the norm ‖𝑇𝑢−𝑏‖. Note that this vector 𝑢 is unique. So, in this case we consider the equation 𝑇𝑥=𝑃_𝑅(𝑇)𝑏, where 𝑃_𝑅(𝑇) is the orthogonal projection on ℛ(𝑇). Since we are interested in the distance between 𝑇𝑥 and 𝑏, it is natural to make use of ‖𝑇‖₂ norm.

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

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

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

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

We will make use of this property (i.e., if 𝑇^† is the generalized inverse of 𝑇 then 𝑇^†𝑏=𝑢, where 𝑢 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 𝑇=∑^𝑛_𝑖=1𝑒_𝑖⊗𝑓_𝑖 is a rank-𝑛 operator then its generalized inverse is also a rank-𝑛 operator and for each 𝑥∈ℋ, it is defined by the relation 𝑇^†𝑥=^𝑛_𝑖=1𝜆_𝑖(𝑥)𝑒_𝑖, (2.2) where the functions 𝜆_𝑖 are the solution of an appropriately defined 𝑛×𝑛 linear system.

The proposed algorithm that implements the ideas of the above theorem, can be found in [2], and is calculating the corresponding 𝜆_𝑖's.

In the scientific area of image processing the analytical form of a linear degraded image is given by the following integral equation:

𝑥_out(𝑖,𝑗)=_𝐷𝑥_in(𝑢,𝑣)ℎ(𝑖,𝑗;𝑢,𝑣)𝑑𝑢𝑑𝑣, (2.3) where 𝑥_in(𝑢,𝑣) is the original image, 𝑥_out(𝑖,𝑗) represents the measured data from where the original image will be reconstructed, and ℎ(𝑖,𝑗;𝑢,𝑣) is a known Point Spread Function (PSF). 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 image 𝑥_in(𝑢,𝑣) 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 𝑟 rows and 𝑚 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:

𝑥_out(𝑖,𝑗)=^𝑟𝑚_𝑘=1𝑙=1𝑥_in(𝑢,𝑣)ℎ(𝑖,𝑗;𝑢,𝑣), (2.4) where 𝑖=1,2,…,𝑟 and 𝑗=1,2,…,𝑚.

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,

𝑥_out(𝑖,𝑗)=^𝑟𝑚_𝑘=1𝑙=1𝑥_in(𝑢,𝑣)ℎ(𝑖−𝑢,𝑗−𝑣)=𝑥_in(𝑖,𝑗)∗∗ℎ(𝑖,𝑗), (2.5) where ∗∗ indicates two-dimensional convolution.

In the formulation of (2.4), the noise can also be simulated by rewriting the equation as

𝑥_out(𝑖,𝑗)=^𝑟𝑚_𝑘=1𝑙=1𝑥_in(𝑢,𝑣)ℎ(𝑖,𝑗;𝑢,𝑣)+𝑛(𝑖,𝑗)=𝑥_in(𝑖,𝑗)∗∗ℎ(𝑖,𝑗)+𝑛(𝑖,𝑗), (2.6) where 𝑛(𝑖,𝑗) 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.

Given 𝑥_out (the digitized degraded X-ray image), then 𝑥_in is the deterministic original image that has to be recovered. The relation between these two components in matrix structure is the following:

𝐻𝑥_in=𝑥_out, (3.1) where 𝐻 represents a two-dimensional (𝑟×𝑚) priori knowledge matrix or it can be estimated from the degraded X-ray image using its Fourier spectrum [15]. The vector 𝑥_out is of 𝑟 entries, while the vector 𝑥_in is of 𝑚(=𝑟+𝑛−1) entries, where 𝑚>𝑟 and 𝑛 is the length of the blurring process in pixels. The problem consists of solving the underdetermined system of (3.1).

However, since there is an infinite number of exact solutions for 𝑥_in that satisfy the equation 𝐻𝑥_in=𝑥_out, 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 𝑟×𝑚 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,

min(‖𝑥_in−𝑥_out‖), (3.2) where 𝑥_in are the first 𝑟 elements of the unknown image 𝑥_in that has to be recovered subject to the constraint ‖𝐻𝑥_in−𝑥_out‖=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:

𝑊(𝑥_in)=‖𝑃𝑥_in−𝑥_out‖₂+𝜆‖𝐻𝑥_in−𝑥_out‖₂, (3.3) where 𝜆 (the Lagrange multiplier) must be large enough. The solution to this problem is the following:

𝑥_in=[𝜆𝐻^𝑇𝐻+𝑃^𝑇𝑃]⁻¹[𝜆𝐻+𝑃]^𝑇𝑥_out, (3.4) where 𝑃 is a projection 𝑟×𝑚 matrix, projecting on the vector 𝑥_out:

𝑃=[𝐼_𝑟|0⋯00⋯0⋮⋮⋮0⋯0], (3.5) and 𝐼_𝑟 is the 𝑟×𝑟 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

𝑘₁⋯𝑘_𝑛00000𝑘₁⋯𝑘_𝑛00000𝑘₁⋯𝑘_𝑛00⋮⋮⋮⋮⋮⋮⋮000⋯𝑘₁⋯𝑘_𝑛⋅𝑥_{in_1}𝑥_{in_2}𝑥_{in_3}⋮𝑥_{in_𝑚}=𝑥_{out_1}𝑥_{out_2}𝑥_{out_3}⋮𝑥_{out_𝑟}, (4.1) where the index 𝑛 indicates the linear motion blur in pixels. The element 𝑘₁,…,𝑘_𝑛 of the matrix are defined as 𝑘_𝑙=1/𝑛(1≤𝑙≤𝑛).

However, (4.1) can also be written in the pointwise form for 𝑖=1,…,𝑟,

𝑥_out1(𝑖)=𝑛^𝑛−1_ℎ=0𝑥_in(𝑖+ℎ), (4.2) which describes an underdetermined system of 𝑟 simultaneous equations and 𝑚=𝑟+𝑛−1 unknowns. The objective is to calculate the original column per column data of the image.

For this reason, given each column [𝑥_{out_1},𝑥_{out_2},𝑥_{out_3},…,𝑥_{out_𝑟}]^𝑇 of a degraded blurred image 𝑥_out, (4.1) results the corresponding column [𝑥_{in_1},𝑥_{in_2},𝑥_{in_3},…,𝑥_{in_𝑚}]^𝑇 of the original image.

As we have seen, the matrix 𝐻 is an 𝑟×𝑚 matrix, and the rank of 𝐻 is less or equal to 𝑚. 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 𝑥_out∈ℛ(𝐻). The Moore-Penrose Inverse is clearly suitable for our case, since we can have a minimum norm solution for every 𝑥_out∈ℛ(𝐻), and a minimal norm least squares solution for every 𝑥_out∉ℛ(𝐻).

The proposed algorithm has been tested on a simulated blurred image produced by applying the matrix 𝐻 on the original image. This can be represented as

𝑥_out1(𝑖,𝑗)=𝑛^𝑛−1_ℎ=0𝑥_in(𝑖,𝑗+ℎ), (4.3) where 𝑖=1,…,𝑟,𝑗=1,…,𝑚 for 𝑚=𝑟+𝑛−1, and 𝑛 is the linear motion blur in pixels.

Following the above, and the analysis given in Section 3, there is an infinite number of exact solutions for 𝑥_in that satisfy 𝐻𝑥_in=𝑥_out, but from Proposition 2.2, only one of them minimizes the norm ‖𝐻𝑥_in−𝑥_out‖.

We will denote this unique vector by 𝑥_in. So, 𝑥_in can be easily found from the equation:

𝑥_in=𝐻^†𝑥_out. (4.4)

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. Figure 1, Original Image, shows such a deterministic original X-ray image. The image has been divided into 𝑟 = 900 rows and 𝑚 = 1200 columns. In Figure 1, Degraded Image, we present the degraded X-ray image for 𝑛=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.

(a)

(b)

(c)

(d)

Figure 1 

Restoration of a simulated degraded image for length of the blurring process 𝑛=30.

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:

ISNR=10log₁₀{∑_𝑖,𝑗[𝑥_in(𝑖,𝑗)−𝑥_out(𝑖,𝑗)]²∑_𝑖,𝑗[𝑥_in(𝑖,𝑗)−𝑥_in(𝑖,𝑗)]²}, (5.1) where 𝑥_in and 𝑥_out represent the original deterministic image and degraded image respectively, and 𝑥_in is the corresponding restored image. Figure 2 shows the corresponding ISNR and computational time values.

(a)

(b)

Figure 2 

ISNR and computational time versus number of pixels in the blurring process (𝑛=1,…,60).

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 of 𝑛. It is clear that the computational load is less by using the Moore-Penrose inverse matrix.

In Table 1, 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+𝑎), where 𝑎=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.

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.

(a)

(b)

Figure 3 

Restoration of a simulated degraded (𝑛=40) and noisy (salt and pepper) image.

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. Figure 3, 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 𝑛.

(a)

(b)

Figure 4 

ISNR and computational time of a degraded and noisy image versus number of pixels in the blurring process.

In order to justify the time efficiency of the proposed method in Table 2 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.

The second set of tests aimed at accenting the reconstruction error between the original image 𝑥_in and the reconstructed image 𝑥_in for various values of linear motion blur, 𝑛. The calculated quantity is the normalized reconstruction error given by

1𝐸=∑^𝑟_𝑖=1∑^𝑚_𝑗=1[𝑥_in(𝑖,𝑗)]^{2𝑟𝑚𝑖=1𝑗=1}[𝑥_in(𝑖,𝑗)−𝑥_in(𝑖,𝑗)]², (5.2) using the generalized inverse reconstructed method.

Figure 5 shows the reconstruction error by increasing the number of pixels in the blurring process 𝑛=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 rank 𝑟×𝑚 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.