Image Reconstruction Methods for MATLAB Users - A Moore-Penrose Inverse Approach

1. Introduction

In the last decades the Moore-Penrose pseudoinverse has found a wide range of applications in many areas of Science and became a useful tool for different scientists dealing with optimization problems, data analysis, solutions of linear integral equations, etc. At first we will present a review of some of the basic results on the so-called Moore-Penrose pseudoinverse of matrices, a concept that generalizes the usual notion of inverse of a square matrix, but that is also applicable to singular square matrices or even to non-square matrices.

The notion of the generalized inverse of a (square or rectangular) matrix was first introduced by H. Moore in 1920, and again by R. 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 $A$be a $r\times m$real matrix. Equations of the form $Ax=b,A\in {\mathbb{R}}^{r\times m},b\in {\mathbb{R}}^r$occur in many pure and applied problems. It is known that when $T$is singular, then its unique generalized inverse $A^{†}$ (known as the Moore- Penrose inverse) is defined. In the case when $A$is a real $r\times m$matrix, Penrose showed that there is a unique matrix satisfying the four Penrose equations, called the generalized inverse of$A$, noted by$A^{†}$.

An important question for applications is to find a general and algorithmically simple way to compute$A^{†}$. There are several methods for computing the Moore-Penrose inverse matrix (cf. [2]). The most common approach uses the Singular Values Decomposition (SVD). 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. Therefore, many other methods can be used for the numerical computation of various types of generalized inverses, see [16]; [25]; [30]. For more on the Moore-Penrose inverse, generalized inverses in general and their applications, there are many excellent texbooks on this subject, see [2]; [30]; [4].

The Moore-Penrose pseudoinverse is a useful concept in dealing with optimization problems, as the determination of a “least squares” solution of linear systems. A typical application of the Moore-Penrose inverse is its use in Image and signal Processing and Image restoration.

The field of image restoration has seen a tremendous growth in interest over the last two decades, see [1]; [5]; [6]; [14]; [28]; [29]. 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 [10]; [15]; [22].

The presented method in this article is based on the use of the Moore-Penrose generalized inverse of a matrix and provides us a fast computational algorithm for a fast and accurate digital image restoration. This article is an extension of the work presented in [7]; [8].

2. Theoretical background

2.1. The Moore-Penrose inverse

We shall denote by ${\mathbb{R}}^{r\times m}$ the algebra of all $r\times m$real matrices. For$T\in {\mathbb{R}}^{r\times m}$, $R\left(T\right)$will denote the range of $T$and $N\left(T\right)$the kernel of $T.$The generalized inverse $T^{†}$ is the unique matrix that satisfies the following four Penrose equations:

$$T{T}^{†}=(T{T}^{†}{)}^{\mathrm{*}},\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }{T}^{†}T=(T^{†}T{)}^{\mathrm{*}},\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }T{T}^{†}T=T,\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }{T}^{†}T{T}^{†}=T^{†},$$E1

where $T^{\mathrm{*}}$ denotes the transpose matrix of $T.$

Let us consider the equation$Tx=b,T\in {\mathbb{R}}^{r\times m},b\in {\mathbb{R}}^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\in R\left(T\right)$or not, and on the rank of$T$. But if$b\notin R\left(T\right)$, then the equation has no solution. Therefore, another point of view of this problem is the following: instead of trying to solve the equation$\parallel Tx-b\parallel =0$, we are looking for a minimal norm vector $u$that minimizes the norm$\parallel Tu-b\parallel$. Note that this vector $u$is unique. So, in this case we consider the equation$Tx=P_{R\left(T\right)}b$, where $P_{R\left(T\right)}$ is the orthogonal projection on$\mathcal{R}\left(T\right)$. Since we are interested in the distance between $Tx$and$b$, it is natural to make use of $\parallel T{\parallel }_2$ norm.

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

Proposition 0.1 Let $T\in {\mathbb{R}}^{r\times m}$and$b\in {\mathbb{R}}^r,b\notin R\left(T\right)$. Then, for$u\in {\mathbb{R}}^m$, the following are equivalent:

1. $$Tu=P_{R\left(T\right)}b$$E2
2. $$\parallel Tu-b\parallel \le \parallel Tx-b\parallel ,\forall x\in {\mathbb{R}}^m$$E3
3. $$T^{\mathrm{*}}Tu=T^{\mathrm{*}}b$$E4

Let$\mathbb{B}=\{u\in {\mathbb{R}}^m|T^{\mathrm{*}}Tu=T^{\mathrm{*}}b\}$. This set of solutions is closed and convex; it therefore has a unique vector $u_0$ with minimal norm. In fact, $\mathbb{B}$is an affine manifold; it is of the form$u_0+\mathcal{N}\left(T\right)$. In the literature (c.f. [12]), $\mathbb{B}$is known as the set of the least square solutions. Proposition 0.2 Let $T\in {\mathbb{R}}^{r\times m}$and$b\in {\mathbb{R}}^r,b\notin R\left(T\right)$, and the equation$Tx=b$. Then, if $T^{†}$ is the generalized inverse of$T$, we have 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.

2.2. Image restoration problems

The general pointwise definition of the transform $\tau (u,v)$that is used in order to convert an $r\times 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:

$$\tau (u,v)=\frac{1}{r}\sum _{x=1}^r\mathrm{‍}\sum _{y=1}^r\mathrm{‍}s(x,y)g(x,y,u,v)$$E5

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:

$$s(x,y)=\frac{1}{r}\sum _{u=1}^r\mathrm{‍}\sum _{v=1}^r\mathrm{‍}\tau (u,v)h(x,y,u,v)$$E6

where $h(x,y,u,v)$represents the inverse of the basis function$g(x,y,u,v)$.

The two dimensional version of the function $g(x,y,u,v)$in Equation (1) 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 been performed across $x$

$$\tau \mathrm{\text{'}}(u,y)=\frac{1}{\sqrt{r}}\sum _{x=1}^r\mathrm{‍}s(x,y)g(x,u)$$E7

Moreover we transform across$y$:

$$\tau (u,v)=\frac{1}{\sqrt{r}}\sum _{y=1}^r\mathrm{‍}\tau \mathrm{\text{'}}(u,y)g(y,u)$$E8

and using Equation (3) we have

$$\tau (u,v)=\frac{1}{r}\sum _{x=1}^r\mathrm{‍}\sum _{y=1}^r\mathrm{‍}s(x,y)g(x,u)g(y,u)$$E9

We can use an identical approach in order to write Equation (1) and its inverse (Equation 2) in matrix form, using the standard orthonormal basis:

$$T=GS{G}^T,\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }S=HT{H}^T$$E10

in which $T,S,G$and $H$are the matrix equivalents of $\tau ,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:

$$S=HGS{G}^T{H}^T$$E11

In order for the transform to be reversible we need $H$to be the inverse of $G$and $H^T$ to be the inverse of$G^T$, i.e.,$HG=G^T{H}^T=I$.

Given that G is orthogonal it is trivial to show that this is satisfied when$H=G^T$. 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.

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

$$x_{out}(i,j)=\int \mathrm{‍}{\int }_D\mathrm{‍}{x}_{in}(u,v)h(i,j;u,v)dudv$$E12

where $x_{in}(u,v)$ is the original image, $x_{out}(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 image $x_{in}(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

$$x_{out}(i,j)=\sum _{u=1}^r\mathrm{‍}\sum _{v=1}^m\mathrm{‍}{x}_{in}(u,v)h(i,j;u,v)$$E13

where $i=\mathrm{1,2},\dots ,r$and $j=\mathrm{1,2},\dots ,m.$

In this work 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 i.e.,

$$x_{out}(i,j)=\sum _{u=1}^r\mathrm{‍}\sum _{v=1}^m\mathrm{‍}{x}_{in}(u,v)h(i-u,j-v)=x_{in}(i,j)\mathrm{*}\mathrm{*}h(i,j)$$E14

where $\mathrm{*}\mathrm{*}$ indicates two dimensional convolution.

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

$$x_{out}(i,j)=\sum _{u=1}^r\mathrm{‍}\sum _{v=1}^m\mathrm{‍}{x}_{in}(u,v)h(i,j;u,v)+n(i,j)=x_{in}(i,j)\mathrm{*}\mathrm{*}h(i,j)+n(i,j)$$E15

where $n(i,j)$is an additive noise introduced by the system.

However, in this work the noise is image related which means that the noise has been added to the image.

2.3. The Fourier Transform, the Haar basis and the moments in image reconstruction problems

Moments are particularly popular due to their compact description, their capability to select differing levels of detail and their known performance attributes (see [3]; [9];[17]; [18]; [19]; [20]; [26]; [27]; [28]. It is a well-recognised property of moments that they can be used to reconstruct the original function, i.e., 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 ([23], ) 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 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 ([19]; [21]). 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.

2.3.1. The Fourier Basis

In view of the importance of the frequency domain, the Fourier Transform (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 $s\left(t\right)$is given by

$$\tilde{s}\left(\omega \right)=\frac{1}{\sqrt{2\pi }}{\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}\mathrm{‍}s\left(t\right)exp(-i\omega t)dt$$E16

where $\omega =2\pi f$is the angular frequency. Since the set of exponentials forms an orthogonal basis the signal can be reconstructed from the projection values

$$s\left(t\right)=\frac{1}{\sqrt{2\pi }}{\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}\mathrm{‍}\tilde{s}\left(\omega \right)exp\left(i\omega t\right)d\omega$$E17

Following the property of the FT that the convolution in the spatial domain is translated into simple algebraic product in the spectral domain Equation (8) can be written in the form

$${\tilde{x}}_{out}={\tilde{x}}_{in}\tilde{H}$$E18

In a discrete Fourier domain the two-dimensional Fourier coefficients are defined as

$$F(m,n)=\frac{1}{\sqrt{XY}}\sum _{x=1}^X\mathrm{‍}\sum _{y=1}^Y\mathrm{‍}{S}_{XY}exp(-2\pi i(\frac{(x-1)(m-1)}{X}+\frac{(y-1)(n-1)}{Y}\left)\right)$$E19

rearranging the above equation leads to

$$F(m,n)=\frac{1}{\sqrt{XY}}\sum _{x=1}^X\mathrm{‍}exp(-2\pi i\frac{(x-1)(m-1)}{X})\sum _{y=1}^Y\mathrm{‍}{S}_{XY}exp(-2\pi i\frac{(y-1)(n-1)}{Y}))$$E20

thus, $F(m,n)$can be written in matrix form as:

$$F(m,n)=K_S(x,m)S_{XY}{K}_S(y,n{)}^{\mathrm{*}}$$E21

where $K_S(y,n{)}^{\mathrm{*}}$ denotes the conjugate transpose of the forward kernel $K_S(y,n).$

Using the same principles but writing Equation (12) in a form where the increasing indexes correspond to higher frequency coefficients we obtain

$$F(m,n)=\frac{1}{\sqrt{XY}}\sum _{x=1}^X\mathrm{‍}\sum _{y=1}^Y\mathrm{‍}{S}_{XY}\cdot \mathrm{e}\mathrm{x}\mathrm{p}[-2\pi i(\frac{(x-1)(m-\frac{(k-1)}{2}-1)}{X}+\frac{(y-\frac{(l-1)}{2}-1)(n-1)}{Y}\left)\right]$$E22

The Fourier coefficients can be seen as the projection coefficients of the image $S_{XY}$ onto a set of complex exponential basis functions that lead to the basis matrix:

$$B_{kl}(m,n)=\frac{1}{\sqrt{k}}exp[-2\pi i\frac{(m-1)(n-\frac{(l-1)}{2}-1)}{k}]$$E23

The approximation of an image $S_{XY}$ in the least square sense, can be expressed in terms of the projection matrix $P_{kl}$ :

$$P_{kl}=(B_{Xk}{)}^T{S}_{XY}{B}_{Yl}$$E24

as

$${S\mathrm{\text{'}}}_{XY}=B_{Xk}\left(B_{Xk}^T{B}_{Xk}^T{)}^{-1}{P}_{kl}\right(B_{Yl}^T{B}_{Yl}{)}^{-1}{B}_{Yl}^T=\left(B_{Xk}{)}^{-}{P}_{kl}\right(B_{Yl}{)}^{†}$$E25

where $({)}^T$ and $({)}^{-1}$ denote the transpose and the inverse of the given matrix. The operations $({)}^{-}$ and $({)}^{†}$ stand for the left and right inverses, both are equal to the Moore-Penrose inverse, 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.

2.3.2. The Haar basis

The reconstruction of an image from its moments is not necessarily unique. Thus, all possible methods must impose extra constraints in order to its moments uniquely solve the reconstruction problem. In this method the constraint that introduced is related to the number of coefficients and the spatial resolution of the image. The Haar basis is unique among the functions we have examined as it actually defines what is referred to as a ’wavelet’. Wavelet functions are a class of functions in which a ’mother’ function is translated and scaled to produce the full set of values required for the full basis set. Limiting the resolution of an image means eliminating those regions of smaller size than a given one. The Haar coefficients are obtained from the projection of the image onto the discrete Haar functions $B_{k,l}\left(m\right)$ for $k$a power of 2, and are defined as

$$B_{k,l}\left(m\right)=\frac{1}{\sqrt{k}},$$E26

in the case$l=1$, and for $l>1$

$$B_{k,l}\left(m\right)=\left\{\begin{array}{llll}+\sqrt{\frac{q}{k}},& if& p\le m<p+\frac{k}{2q}& \\ -\sqrt{\frac{q}{k}},& if& p+\frac{k}{2q}\le m\le p+\frac{k}{q}& \\ 0,& otherwise& & \end{array}\right.$$E27

with $q=2^{\left[lo{g}_2\right(l-1\left)\right]}$and$p=\frac{k(l-1-q)}{q}+1$, where $[.]$ stands for the function$\mathrm{f}ix\left(x\right)$, which rounds the elements of $x$to the nearest integer towards zero.

3. Restoration of a blurry image in the spatial domain

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$x_{out}$, then $x_{in}$ is the deterministic original image that has to be recovered. The relation between these two components in matrix structure is the following :

$$H{x}_{in}=x_{out},$$E28

where $H$represents a two dimensional $(r\times m)$ priori knowledge matrix or it can be estimated from the degraded X-ray image using its Fourier spectrum ([24]). The vector$x_{out}$, is of $r$entries, while the vector $x_{in}$ 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 equations (Eq. 13).

However, since there is an infinite number of exact solutions for $x_{in}$ that satisfy the equation$H{x}_{in}=x_{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 $r\times 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, i.e.,

$$\mathrm{m}\mathrm{i}\mathrm{n}\left(‖x_{in}^{\mathrm{*}}-x_{out}‖\right),$$E29

where $x_{in}^{\mathrm{*}}$ are the first $r$elements of the unknown image $x_{in}$ that has to be recovered subject to the constraint$‖H{x}_{in}-x_{out}‖=0$. In fact, zero is not always attained, but following Proposition 0.1(ii) the norm is minimized.

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.

3.1. The generalized inverse approach

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 follows, as desribed in Eq. (13):

$$\left[\begin{array}{lllllll}{k}_1& \dots & k_n& 0& 0& 0& 0\\ 0& k_1& \dots & k_n& 0& 0& 0\\ 0& 0& k_1& \dots & k_n& 0& 0\\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮\\ 0& 0& 0& \dots & k_1& \dots & k_n\\ & & & & & & \end{array}\right]\cdot \left[\begin{array}{l}{x}_{in\_1}\\ x_{in\_2}\\ x_{in\_3}\\ ⋮\\ x_{in\_m}\\ \end{array}\right]=\left[\begin{array}{l}{x}_{out\_1}\\ x_{out\_2}\\ x_{out\_3}\\ ⋮\\ x_{out\_r}\\ \end{array}\right]$$E30

where the index $n$indicates the linear motion blur in pixels. The element $k_1,\dots ,k_n$ of the matrix are defined as:$k_l=1/n\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }(1\le l\le n)$.

Equation (3) can also be written in the pointwise form for$i=1,\dots ,r$,

$$x_{out}\left(i\right)=\frac{1}{n}\sum _{h=0}^{n-1}\mathrm{‍}{x}_{in}(i+h)$$E31

that 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 $[x_{out\_1},x_{out\_2},x_{out\_3},\dots x_{out\_r}{]}^T$ of a degraded blurred image$x_{out}$, Eq. (3) results the corresponding column

$[x_{in\_1},x_{in\_2},x_{in\_3},\dots ,x_{in\_m}{]}^T$of the original image.

As we have seen, the matrix $H$is a $r\times 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$x_{out}\in \mathcal{R}\left(H\right)$. The Moore-Penrose Inverse is clearly suitable for our case, since we can have a minimum norm solution for every$x_{out}\in \mathcal{R}\left(H\right)$, and a minimal norm least squares solution for every$x_{out}\notin \mathcal{R}\left(H\right)$.

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

$$x_{out}(i,j)=\frac{1}{n}\sum _{h=0}^{n-1}\mathrm{‍}{x}_{in}(i,j+h)$$E32

where $i=1,\dots ,r\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }j=1,\dots ,m$for$m=r+n-1$, and $n$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 $x_{in}$ that satisfy the equation$H{x}_{in}=x_{out}$, but from proposition 2.2, only one of them minimizes the norm$\parallel H{x}_{in}-x_{out}\parallel$.

We shall denote this unique vector by${\widehat{x}}_{in}$. So, ${\widehat{x}}_{in}$can be easily found from the equation :

$${\widehat{x}}_{in}=H^{†}{x}_{out}$$E33

The following section presents results that highlight the performance of the generalized inverse.

4. Experimental results

In this section we apply the proposed method on an boat picture and present the numerical results.

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.

4.1. Recovery from a degraded image

Figure 1(a) provides the original boat picture. In Figure 1(b), we present the degraded boat picture where the length of the blurring process is equal to$n=60$. Finally, in Figure 1(c) we present the reconstructed image using the Moore- Penrose inverse approach. As we can see, it is clearly seen that the details of the original image have been recovered.

Figure 1.

a) Original Image (b) Blurred image for a length of the blurring process n = 60 (c) Restoration of a simulated degraded image with a length of the blurring process n = 60.

The Improvement in Signal to Noise Ratio (ISNR) has been chosen in order to present the reconstructed images obtained by the proposed algorithm. It provides a criterion that has been used extensively for the purpose of objectively testing the performance of image processing algorithms expressed as:

$$ISNR=10{\mathrm{l}\mathrm{o}\mathrm{g}}_{10}\left\{\frac{\sum _{i,j}\mathrm{‍}{\left[x_{in}(i,j)-x_{out}(i,j)\right]}^2}{\sum _{i,j}\mathrm{‍}{\left[x_{in}(i,j)-{\tilde{x}}_{in}(i,j)\right]}^2}\right\},$$E34

where $x_{in}$ and $x_{out}$ represent the original deterministic image and degraded image respectively, and ${\widehat{x}}_{in}$ is the corresponding restored image. Figure 2(a) shows the corresponding ISNR values. for increasing the number of pixels in the blurring process$n=1,\dots ,60$.

The second set of tests aimed at accenting the reconstruction error between the original image $x_{in}$ and the reconstructed image ${\widehat{x}}_{in}$ for various values of linear motion blur,$n$. The calculated quantity is the normalized reconstruction error given by

$$E=\frac{1}{\sqrt{\sum _{i=1}^r\mathrm{‍}\sum _{j=1}^m\mathrm{‍}\left[x_{in}\right(i,j){]}^2}}\sqrt{\sum _{i=1}^r\mathrm{‍}\sum _{j=1}^m\mathrm{‍}\left[x_{in}\right(i,j)-{\widehat{x}}_{in}(i,j){]}^2}$$E35

using the generalized inverse reconstructed method.

Figure 2(b) shows the reconstruction error by increasing the number of pixels in the blurring process $n=1,\dots ,60.$

Figure 2.

a) ISNR and (b) Reconstruction Error calculations vs number of pixels in the blurring process$(n=1,\dots ,60)$.

4.2. Recovery from a degraded and noisy image

Noise may be introduced into an image in a number of different ways. In Equation (10) 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(b), we present the original boat image while a motion blurred and a salt and pepper noise has been added to it.

Figure 3.

a) Noisy Image (b) Blurred and noisy (salt and pepper) image for length of the blurring process n = 60 (c) Restoration of a simulated degraded (n = 60) 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. Accordingly, Figure 4(a) provides a graphical representation for the ISNR of the reconstructed and filtered image for different values of$n$. Moreover, Figure 4(b) shows the reconstruction error by increasing the number of pixels in the blurring process$n=1,\dots ,60$.

Figure 4.

a) ISNR and (b) Reconstruction Error calculations for a noisy and blurred image vs number of pixels in the blurring process$(n=1,\dots ,60)$.

5. Deblurring in the spatial and spectral domain: Application of the Haar and Fourier moments on image reconstruction.

As mentioned before, images can be viewed as non-stationary two-dimensional signals with edges, textures, and deterministic objects at different locations. 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 objective is to calculate the inverse matrix of the blurring kernel H and then applied back (simple multiplication in the spectral domain) to the degraded blurred image xout. Figure 5 shows the spectral representation of the degraded image obtained using Equation (11).

In order to obtain back the original image, Equation (13) is solved in the Fourier space

$${\tilde{x}}_{in}={\tilde{x}}_{out}{\tilde{H}}^{†}$$E36

The reconstructed image is the inverse Fourier transform of${\tilde{x}}_{in}$. By using our method not only we have the advantage of fast recovery but also provide us with an operator ${\tilde{H}}^{†}$ that exists even for not full rank non square matrices. In this section the whole process of deblurring and restoring the original image is done in the spatial domain by using the Haar basis moments and in the spectral domain by applied the Fourier basis moments on the image. It provides us the ability of fast recovering and algorithmic simplicity. The former, obtained by using directly the original image and analysed that on its moments. The method is robust in the presence of noise, as can be seen on the results. In the latter, 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. The method provides almost the same robustness for the case of degradation and noise presence as for the spatial moment analysis case.

Figure 5.

Spectral representations of the degraded image for n=30.

Figures $6\left(a\right),6\left(b\right)$ and $6\left(c\right)$ present the reconstructed image using the Fourier basis,for the cases of $k=l=30,k=l=100$and$k=l=450$, respectively.

Figure 6.

Fourier based moment reconstructed images for (a) k = l = 30 (b) k = l = 100 and (c) k = l = 450.

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.

Figures $7\left(a\right),7\left(b\right)$ and $7\left(c\right)$ present the reconstructed image using the Haar basis,for the cases of $k=l=30,k=l=100$and$k=l=450$, respectively.

Figure 7.

Haar based moment reconstructed images for (a) k = l = 30 (b) k = l = 100 and (c) k = l = 450.

Figure 8.

a) ISNR and (b) Reconstruction Error calculations for a noisy and blurred image vs number of pixels in the blurring process$(n=1,\dots ,60)$. The blue and red lines indicate the usage of Fourier and Haar based moment analysis of the image, respectively.

Figure 9.

a) ISNR and (b) Reconstruction Error calculations for a noisy and blurred image vs number of moments$(k=l=1,\dots ,200)$. The blue and red lines indicate the usage of Fourier and Haar based moment analysis of the image, respectively.

Figures 8(a) and 8(b) show the ISNR and the Reconstruction Error accordingly, for various lengths of the blurring processes. Graphical representations on these Figures correspond to moment values $k=l=450$(blue line for the Fourier moment and red line for the Haar moment case). The image is corrupted with white and black (salt and pepper) noise with noise density equal to 0.02. After the moment analysis took place a low pass rotationally symmetric Gaussian filter of standard deviation equal to 45 were applied. Finally, on Figures 9(a) and 9(b) we present the ISNR and the Reconstruction Error respectively, for a number of moments, $k=l=1,\dots ,200$and keeping the number of blurring process at a high level equal to$n=60$. Similarly, to the previous cases the value of the black and white noise density is equal to the 0.02 and a low-pass Gaussian filter was used for the filtering process.

6. Conclusions

In this study, we introduced a novel computational method based on the calculation of the Moore-Penrose inverse of full rank $r\times 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. The efficiency of the generalized inverse is evidenced by the presented simulation results. In this chapter the results presented were demonstrated in the spatial and spectral domain of the image. 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 Haar coefficients (spatial domain) or Fourier coefficients (spectral domain), 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 so the dimensions and the nature of the image do not play any role in this application

Appendix

In this section we provide the interested readers with the Matlab codes used in this article.

The following Matlab functions where used to calculate the Fourier and the Haar basis coefficients, and the blurring matrix of the images used.

Function that calculates the Fourier Basis Coefficients (FBC) of an image.

Function that calculates the Haar Basis Coefficients (HBC) of an image.

Function that calculates the blurring matrix of an image.