Applications of the Moore-Penrose Inverse in Digital Image Restoration
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.
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. ). 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 ) 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:
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 , 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 .
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. ), 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 ]). Let be a Hilbert space. If is a rank- operator then its generalized inverse is also a rank- operator and for each , it is defined by the relation 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 , 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:
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). 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 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:
In this paper we adopt the use of a shift invariant model for the blurring process as in . Therefore, the analytically expression for the degraded system is given by a two-dimensional (horizontal and vertical) convolution, that is,
where indicates two-dimensional convolution.
In the formulation of (2.4), the noise can also be simulated by rewriting the equation as
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 (the digitized degraded X-ray 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 priori knowledge matrix or it can be estimated from the degraded X-ray image using its Fourier spectrum . 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.1).
However, since there is an infinite number of exact solutions for that satisfy the equation , 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,
where are the first elements of the unknown image that has to be recovered subject to the constraint 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:
where (the Lagrange multiplier) must be large enough. The solution to this problem is the following:
where is a projection matrix, projecting on the vector :
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
where the index indicates the linear motion blur in pixels. The element of the matrix are defined as .
However, (4.1) can also be written in the pointwise form for
which describes an underdetermined system of simultaneous equations and unknowns. The objective is to calculate the original column per column data of the image.
For this reason, given each column of a degraded blurred image , (4.1) results the corresponding column 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 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 . The Moore-Penrose Inverse is clearly suitable for our case, since we can have a minimum norm solution for every , and a minimal norm least squares solution for every .
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
where for , and is the linear motion blur in pixels.
We will denote this unique vector by . So, can be easily found from the equation:
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 . 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.
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:
where and represent the original deterministic image and degraded image respectively, and is the corresponding restored image. Figure 2 shows the corresponding ISNR and computational time values.
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 , where . 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.
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 . 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.
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 . Accordingly, Figure 4 provides a graphical representation for the ISNR and the computational time, of the reconstructed and filtered image for different values of .
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.
The second set of tests aimed at accenting the reconstruction error between the original image and the reconstructed image for various values of linear motion blur, . The calculated quantity is the normalized reconstruction error given by
using the generalized inverse reconstructed method.
Figure 5 shows the reconstruction error by increasing the number of pixels in the blurring process . The solid line represents the degraded image case and the dotted line the noisy image.
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.
M. R. Banham and A. K. Katsaggelos, “Digital image restoration,” IEEE Signal Processing Magazine, vol. 14, no. 2, pp. 24–41, 1997.View at: Google Scholar
K. R. Castleman, Digital Processing, Prentice-Hall, Eglewood Cliffs, NJ, USA, 1996.
H. J. Trussell and S. Fogel, “Identification and restoration of spatially variant motion blurs in sequential images,” IEEE Transactions of Image Processing, vol. 1, no. 1, pp. 123–126, 1992.View at: Google Scholar
D. L. Tull and A. K. Katsaggelos, “Iterative restoration of fast-moving objects in dynamic image sequences,” Optical Engineering, vol. 35, no. 12, pp. 3460–3469, 1996.View at: Google Scholar
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.View at: Google Scholar
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.View at: Google Scholar
A. K. Katsaggelos, J. Biemond, R. M. Mersereau, and R. W. Schafer, “A general formulation of constrained iterative restoration algorithms,” in Proceedings of IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP '85), pp. 700–703, Tampa, Fla, USA, 1985.View at: Google Scholar
C. W. Groetsch, Generalized Inverses of Linear Operators: Representation and Approximation, Monographs and Textbooks in Pure and Applied Mathematics, no. 3, Marcel Dekker, New York, NY, USA, 1977.View at: MathSciNet
R. Gonzalez and P. Wintz, Digital Processing, Addison-Wesley, Reading, Mass, USA, 2nd edition, 1987.
M. Sondhi, “Restoration: the removal of spatially invariant degradations,” Proceedings of the IEEE, vol. 60, no. 7, pp. 842–853, 1972.View at: Google Scholar