Next Article in Journal
Advances in Hydrogen, Carbon Dioxide, and Hydrocarbon Gas Sensor Technology Using GaN and ZnO-Based Devices
Next Article in Special Issue
Fundamentals of in Situ Digital Camera Methodology for Water Quality Monitoring of Coast and Ocean
Previous Article in Journal
Immobilization of HRP in Mesoporous Silica and Its Application for the Construction of Polyaniline Modified Hydrogen Peroxide Biosensor
Previous Article in Special Issue
CMOS Image Sensor with a Built-in Lane Detector
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

A High Resolution Color Image Restoration Algorithm for Thin TOMBO Imaging Systems

School of Electrical, Electronic and Computer Engineering, The University of Western Australia, 35 Stirling Highway, Crawley, WA 6009, Australia
*
Author to whom correspondence should be addressed.
Sensors 2009, 9(6), 4649-4668; https://doi.org/10.3390/s90604649
Submission received: 20 May 2009 / Revised: 5 June 2009 / Accepted: 5 June 2009 / Published: 15 June 2009
(This article belongs to the Special Issue Image Sensors 2009)

Abstract

:
In this paper, we present a blind image restoration algorithm to reconstruct a high resolution (HR) color image from multiple, low resolution (LR), degraded and noisy images captured by thin (< 1mm) TOMBO imaging systems. The proposed algorithm is an extension of our grayscale algorithm reported in [1] to the case of color images. In this color extension, each Point Spread Function (PSF) of each captured image is assumed to be different from one color component to another and from one imaging unit to the other. For the task of image restoration, we use all spectral information in each captured image to restore each output pixel in the reconstructed HR image, i.e., we use the most efficient global category of point operations. First, the composite RGB color components of each captured image are extracted. A blind estimation technique is then applied to estimate the spectra of each color component and its associated blurring PSF. The estimation process is formed in a way that minimizes significantly the interchannel cross-correlations and additive noise. The estimated PSFs together with advanced interpolation techniques are then combined to compensate for blur and reconstruct a HR color image of the original scene. Finally, a histogram normalization process adjusts the balance between image color components, brightness and contrast. Simulated and experimental results reveal that the proposed algorithm is capable of restoring HR color images from degraded, LR and noisy observations even at low Signal-to-Noise Energy ratios (SNERs). The proposed algorithm uses FFT and only two fundamental image restoration constraints, making it suitable for silicon integration with the TOMBO imager.

Graphical Abstract

1. Introduction

TOMBO (Thin Observation Module by Bound Optics) imaging systems are a new generation of paper-thin imaging products, which integrate micro-optics, photo-sensing elements and processing-circuitry all on one single silicon chip (Figure 1). These paper-thin imagers exhibit a thickness of less than a millimeter because they do not image the scene through a single imaging lens but through an array of microlenses [2]. The concept of a TOMBO imager was proposed and demonstrated by Tanida et al. [38]. The structure of a TOMBO imager is shown in Figure 1. It consists of an array of imaging units, each comprises a microlens that is associated with a subset of photo-sensitive pixel array. Individual imaging units are optically isolated by an opaque wall to prevent crosstalk (Figure 1). As a result, each individual imaging unit visualizes part of the scene. The output of the TOMBO imager is thus a mosaic of low resolution (LR) unit images. To reconstruct a high resolution (HR) image from the acquired set of LR images, Tanida et al. first proposed a Back-Projection (BP) method [6], which requires complete knowledge of the imaging system point spread function (PSF). The HR image of the scene is obtained by multiplying the captured LR images by the inverse (pseudo-inverse) of the known PSF. Tanida et al. proposed a second image reconstruction method (the “pixel rearrange method”) [7], which computes the cross-correlation peaks between captured unit images to arrange and align unit image pixels in the HR image of the scene. A de-shading pre-processing step is introduced to compensate for the shading introduced by the separation walls (Figure 1).
We have previously proposed a novel spectral-based image restoration algorithm that require neither prior information about the imaging system nor the original scene [1]. Furthermore, the proposed algorithm alleviates the need for de-shading and rearrangement processing. In this paper, we extend this algorithm to color images. We examine the difference in characteristics between grayscale and color images to develop a model for the color TOMBO imager. Previous work on color TOMBO imagers directly applied grayscale HR reconstruction algorithms to color images without considering the cross-correlation between color channels, and thus resulted in color artifacts [913]. In this paper, we exploit the global category of point operations for image restoration (see Figure 2) [14]. In this process, each pixel of the restored image is obtained by using information (pixels) from all captured LR images [1520].
The proposed spectral-based color image restoration method averages out all LR captured images, making the color channels globally independent of each other. Compared to previously reported color restoration techniques [9], this proposed algorithm uses FFT and only two fundamental image restoration constraints, which makes it suitable for silicon integration with a TOMBO imager. The paper is organized as follows: in section 2, we develop a model for a color TOMBO imager and derive the mathematical formulation of a novel blind color image restoration method. Section 3 details the color image restoration algorithm. Results are discussed in section 4. Finally, a conclusion is given in section 5.

2. Image Restoration Method for TOMBO Color Imaging Systems

In this section, we extend the grayscale image restoration algorithm reported in [1] to color TOMBO imagers. In the restoration process, we consider the global point operations based on multiple images. By using this category of point operations, we exploit all available information in the mosaic of simultaneously captured color images (see Figure 2). In addition, the global category is reported to have the ability to remove significant additive noise [1520].

2.1. System Model

Consider a TOMBO color system with (μ × μ) captured color images as shown in Figure 1. Each captured color image represents a blurred, LR and noisy observation of an original HR scene. The mathematical model (Figure 3) for the system can be described by
[ g i , j ( x , y , R ) g i , j ( x , y , G ) g i , j ( x , y , B ) ] = { [ h i , j ( x , y , R ) h i , j ( x , y , G R ) h i , j ( x , y , B R ) h i , j ( x , y , G R ) h i , j ( x , y , G ) h i , j ( x , y , B G ) h i , j ( x , y , B R ) h i , j ( x , y , B G ) h i , j ( x , y , B ) ] [ f ( x , y , R ) f ( x , y , G ) f ( x , y , B ) ] + [ v i , j ( x , y , R ) v i , j ( x , y , G ) v i , j ( x , y , B ) ] } D
  • gi,j(x,y,ϑ),ϑ ∈ {R, G, B} represents the blurred, LR and noisy color component for the ith,jth captured unit image with resolution (M × N) pixels per color
  • hi,j(x,y,ϑ) is an (l × l) PSF that represents the overall channel blur affecting gi,j(x,y,ϑ) unit image for the color component ϑ, also called the intrachannel. We assume here that the blur is different for each color of each unit image
  • hi,j(x, y, GR),hi,j(x, y, BR),hi,j(x, y, BG) are (l × l) PSFs representing the overall mutual relation between red-green, red-blue and green-blue respectively.
  • “* *” represents the 2-D convolution operator w.r.t x, y
  • f(x, y, ϑ) is the ϑ color component of the original scene with resolution (M × N) > (M × N) per color component
  • vi,j(x, y, ϑ) is the additive 2-D, zero mean white Gaussian noise per color component that affect the unit image gi,j(x,y,ϑ)
  • D is the down-sampling operator representing the LR process
Our main goal is to develop a restoration method that is able to reconstruct a HR, noiseless, color image of the original scene using only the (μ × μ) LR, blurred and noisy TOMBO color images gi,j (x, y, ϑ). The proposed method will have the following characteristics: (i) it does not require prior information about the imaging system nor the original scene, and thus performs blind image restoration, (ii) it can remove blur and additive noise from the HR color image, (iii) it exploits all available information contained in the captured LR images, and (iv) it requires minimal computational complexity.

2.2. Formulation of the Restoration Method

The general model of the color TOMBO imaging system is described by Equation (2). Let us consider all signal terms involved in the captured red color component of a unit image (i, j):
g i , j ( x , y , R ) = [ h i , j ( x , y , R ) f ( x , y , R ) + h i , j ( x , y , G R ) f ( x , y , G ) + h i , j ( x , y , B R ) f ( x , y , B ) + v i , j ( x , y , R ) ] D , i , j = 1 , , μ
By applying the two dimensional -transform to the model in Equation (2), we have
G i , j ( z 1 , z 2 , R ) = 1 D 2 k = 0 D 1 l = 0 D 1 H i , j ( z 1 1 D e j 2 π k D , z 2 1 D e j 2 π l D , R ) F ( z 1 1 D e j 2 π k D , z 2 1 D e j 2 π l D , R ) + 1 D 2 k = 0 D 1 l = 0 D 1 H i , j ( z 1 1 D e j 2 π k D , z 2 1 D e j 2 π l D , G R ) F ( z 1 1 D e j 2 π k D , z 2 1 D e j 2 π l D , G R ) + 1 D 2 k = 0 D 1 l = 0 D 1 H i , j ( z 1 1 D e j 2 π k D , z 2 1 D e j 2 π l D , B R ) F ( z 1 1 D e j 2 π k D , z 2 1 D e j 2 π l D , B R ) + 1 D 2 k = 0 D 1 l = 0 D 1 V i , j ( z 1 1 D e j 2 π k D , z 2 1 D e j 2 π l D , R ) ,
where z 1 = e j 2 π f 1, z 2 = e j 2 π f 2andGi,j(z1, z2, R), is the -transform of gi,j(x, y, R).
The system model in Equation (3) is partitioned into the following terms:
G i , j ( z 1 , z 2 , R ) = T i , j a ( z 1 1 D , z 2 1 D , R ) + T i , j b ( z 1 1 D , z 2 1 D , R ) + T i , j c ( z 1 1 D , z 2 1 D , R ) + T i , j d ( z 1 1 D , z 2 1 D , G R ) + T i , j e ( z 1 1 D , z 2 1 D , B R )
where,
  • T i , j a ( z 1 1 D , z 2 1 D , R ) = H i , j ( z 1 1 D , z 2 1 D , R ) F ( z 1 1 D , z 2 1 D , R ) + V i , j ( z 1 1 D , z 2 1 D , R ) represents the image of interest plus the noise term (defined in-frequency band useful terms),
  • T i , j b ( z 1 1 D , z 2 1 D , R ) = k = 1 D 1 l = 1 D 1 H i , j ( z 1 1 D e j 2 π k D , z 2 1 D e j 2 π l D , R ) F ( z 1 1 D e j 2 π k D , z 2 1 D e j 2 π l D , R ), are the aliasing out of frequency band image terms,
  • T i , j c ( z 1 1 D , z 2 1 D , R ) = k = 1 D 1 l = 1 D 1 V i , j ( z 1 1 D e j 2 π k D , z 2 1 D e j 2 π l D , R ), are the aliasing out-of-frequency band noise terms.
  • T i , j d ( z 1 1 D , z 2 1 D , G R ) = k = 1 D 1 l = 1 D 1 H i , j ( z 1 1 D e j 2 π k D , z 2 1 D e j 2 π l D , G R ) F ( z 1 1 D e j 2 π k D , z 2 1 D e j 2 π l D , G ), are the GR overall cross-talk terms.
  • T i , j e ( z 1 1 D , z 2 1 D , B R ) = k = 1 D 1 l = 1 D 1 H i , j ( z 1 1 D e j 2 π k D , z 2 1 D e j 2 π l D , B R ) F ( z 1 1 D e j 2 π k D , z 2 1 D e j 2 π l D , B ), are the BR overall cross-talk terms.
In the above terms, the constant 1 D 2 is not shown for simplification.
By multiplying both sides of Equation (4) with the complex conjugate of the red (i,j) blurring PSF, H i , j ( z 1 1 D , z 2 1 D , R ), i.e., H i , j ( z 1 1 D , z 2 1 D , R ) and similarly by F ( z 1 1 D , z 2 1 D , R ) and applying the ensemble average operator, E{}, we have,
E { G i , j H i , j ( z 1 1 D , z 2 1 D , R ) } = E { T i , j a H i , j ( z 1 1 D , z 2 1 D , R ) } + E { T i , j b H i , j ( z 1 1 D , z 2 1 D , R ) } + E { T i , j c H i , j ( z 1 1 D , z 2 1 D , R ) } + E { T i , j d H i , j ( z 1 1 D , z 2 1 D , G R ) } + E { T i , j e H i , j ( z 1 1 D , z 2 1 D , B R ) } E { G i , j F ( z 1 1 D , z 2 1 D , R ) } = E { T i , j a F ( z 1 1 D , z 2 1 D , R ) } + E { T i , j b F ( z 1 1 D , z 2 1 D , R ) } + E { T i , j c F ( z 1 1 D , z 2 1 D , R ) } + E { T i , j d F ( z 1 1 D , z 2 1 D , G R ) } + E { T i , j e F ( z 1 1 D , z 2 1 D , B R ) }
which leads to:
C G i , j H i , j ( z 1 1 D , z 2 1 D , R ) = F ( z 1 1 D , z 2 1 D , R ) C H i , j H i , j ( z 1 1 D , z 2 1 D , R ) + C V i , j H i , j ( z 1 1 D , z 2 1 D , R ) + C T i , j b H i , j ( z 1 1 D , z 2 1 D , R ) + C T i , j c H i , j ( z 1 1 D , z 2 1 D , R ) + C T i , j d H i , j ( z 1 1 D , z 2 1 D , G R ) + C T i , j e H i , j ( z 1 1 D , z 2 1 D , B R ) C G i , j F ( z 1 1 D , z 2 1 D , R ) = H i , j ( z 1 1 D , z 2 1 D , R ) C F F ( z 1 1 D , z 2 1 D , R ) + C V i , j F ( z 1 1 D , z 2 1 D , R ) + C T i , j b F ( z 1 1 D , z 2 1 D , R ) + C T i , j c F ( z 1 1 D , z 2 1 D , R ) + C T i , j d F ( z 1 1 D , z 2 1 D , G R ) + C T i , j e F ( z 1 1 D , z 2 1 D , B R )
where cross-spectra CXY* (z1, z2, R) = E {X(z1, z2, R)Y* (z1, z2, R)} over z1, z2.
Since the frequency bandwidth of the blurring PSF function Hi,j(z1, z2, R) is limited [15, 19, 21], the terms T i , j b ( z 1 1 D , z 2 1 D , R ) and T i , j c ( z 1 1 D , z 2 1 D , R ) are not in the same frequency band where H i , j ( z 1 1 D , z 2 1 D , R ) is. Therefore, the last two cross-spectral terms of Equation (6) can be discarded, leading to:
C G i , j H i , j ( z 1 1 D , z 2 1 D , R ) = F ( z 1 1 D , z 2 1 D , R ) C H i , j H i , j ( z 1 1 D , z 2 1 D , R ) + C V i , j H i , j ( z 1 1 D , z 2 1 D , R ) + C T i , j d H i , j ( z 1 1 D , z 2 1 D , G R ) + C T i , j e H i , j ( z 1 1 D , z 2 1 D , B R ) C G i , j F ( z 1 1 D , z 2 1 D , R ) = H i , j ( z 1 1 D , z 2 1 D , R ) C F F ( z 1 1 D , z 2 1 D , R ) + C V i , j F ( z 1 1 D , z 2 1 D , R ) + C T i , j d F ( z 1 1 D , z 2 1 D , G R ) + C T i , j e F ( z 1 1 D , z 2 1 D , B R )
To minimize the effects due to interchannel cross-correlation terms (cross-talks) C T i , j d H i , j ( z 1 1 D , z 2 1 D , R ) and C T i , j e H i , j ( z 1 1 D , z 2 1 D , R ) in Equation (7), a typical approach would be to apply a decorrelation stretch technique to the captured LR images [14, 22]. However, such technique uses principal component analysis, which is computationally expensive and thus not suited for silicon integration with a TOMBO imager. An alternative solution is to minimize the cross-correlation terms by using ergodicity, one of the wide-sense stationary properties of signals [23]. Let us consider the cross-correlations between the color channels, which depend on the reflection characteristics of the imaged objects. If the object exhibits strong reflection at the region of correlated wavelengths, then the cross-correlation between the channels could be strong. However, since the characteristics of one object is almost independent of those of other objects, the averaged interchannel cross-correlations over the entire image can be very weak. As a result, the RGB components of a color image can be regarded as being locally correlated but globally independent of each other [23]. In this paper, we exploit this property by using a global category of point operations for image restoration. In this process, each pixel of the restored image is obtained by using information (pixels) from all captured LR images (Figure 2). This is carried out by computing the cross-spectra between the different color channels. Since the cross-spectra are nothing but the Fourier transform of the cross-correlations, our algorithm essentially averages out the interchannel cross-correlations over the entire image. As a result, the averaged interchannel cross-correlations become very weak [23]. Our spectral-based restoration method has the effect of further minimizing the interchannel cross-correlations because it averages out over all captured LR images. As a result, the cross-spectral terms of Equation 7, C T i , j d H i , j ( z 1 1 D , z 2 1 D , G R ) the C T i , j e H i , j ( z 1 1 D , z 2 1 D , B R ), C T i , j d F ( z 1 1 D , z 2 1 D , G R ), C T i , j e F ( z 1 1 D , z 2 1 D , B R ) can be neglected, leading to:
C G i , j H i , j ( z 1 1 D , z 2 1 D , R ) = F ( z 1 1 D , z 2 1 D , R ) C H i , j H i , j ( z 1 1 D , z 2 1 D , R ) + C V i , j H i , j ( z 1 1 D , z 2 1 D , R ) C G i , j F ( z 1 1 D , z 2 1 D , R ) = H i , j ( z 1 1 D , z 2 1 D , R ) C F F ( z 1 1 D , z 2 1 D , R ) + C V i , j F ( z 1 1 D , z 2 1 D , R )
where C V i , j H i , j ( z 1 1 D , z 2 1 D , R ) and C V i , j F ( z 1 1 D , z 2 1 D , R ) are the residual errors representing the in-band cross spectral terms between original images, in the PSFs and in the independent additive noise respectively. These errors can be assumed 2-D, independent and identically distributed (i.i.d.) signals under some regularity conditions [24]. By extending the previous analysis to the two other color components G and B, (8) can be generalized:
[ C G i , j H i , j ( z 1 1 D , z 2 1 D , R ) C G i , j H i , j ( z 1 1 D , z 2 1 D , G ) C G i , j H i , j ( z 1 1 D , z 2 1 D , B ) ] [ C H i , j H i , j ( z 1 1 D , z 2 1 D , R ) 0 0 0 C H i , j H i , j ( z 1 1 D , z 2 1 D , G ) 0 0 0 C H i , j H i , j ( z 1 1 D , z 2 1 D , B ) ] [ F ( z 1 1 D , z 2 1 D , R ) F ( z 1 1 D , z 2 1 D , G ) F ( z 1 1 D , z 2 1 D , B ) ] + [ C V i , j H i , j ( z 1 1 D , z 2 1 D , R ) C V i , j H i , j ( z 1 1 D , z 2 1 D , G ) C V i , j H i , j ( z 1 1 D , z 2 1 D , B ) ]
and similarly
[ C G i , j F ( z 1 1 D , z 2 1 D , R ) C G i , j F ( z 1 1 D , z 2 1 D , G ) C G i , j F ( z 1 1 D , z 2 1 D , B ) ] [ C F F ( z 1 1 D , z 2 1 D , R ) 0 0 0 C F F ( z 1 1 D , z 2 1 D , G ) 0 0 0 C F F ( z 1 1 D , z 2 1 D , B ) ] [ H i , j ( z 1 1 D , z 2 1 D , R ) H i , j ( z 1 1 D , z 2 1 D , G ) H i , j ( z 1 1 D , z 2 1 D , B ) ] + [ C V i , j F ( z 1 1 D , z 2 1 D , R ) C V i , j F ( z 1 1 D , z 2 1 D , G ) C V i , j F ( z 1 1 D , z 2 1 D , B ) ]

2.3. Restoration Process

From Equations (8), (9), and (10), if some prior information or constraints about either the original or the blur is given, it is then possible to restore the original image. This can be done independently for each color component or jointly for all color components using Equations (9), and (10). For example, if the PSF associated with each captured gi,j(x,y, R) unit image, H i , j ( z 1 1 D , z 2 1 D , R ), is pre-determined or can be estimated, then it becomes possible to restore the original image using,
F ( z 1 1 D , z 2 1 D , R ) = C G i , j H i , j ( z 1 , z 2 , R ) C V i , j H i , j ( z 1 1 D , z 2 1 D , R ) C H i , j H i , j ( z 1 1 D , z 2 1 D , R )
The above equation is valid under some constraints [24, 25]. One constraint is that the number of zeros in Hi,j(z1, z2, R) is much less than the size of the image. Because our method sums all PSFs spectra when implementing Equation (11), this will minimize the probability of having a division by zero. In general, if zeros exist but only in small numbers compared to the image size, a tolerance value can be added to the denominator to avoid division by zero. Other problems associated with image and PSF restoration methods are (i) the length £ and the PSF are unknown and most critically (ii) the impact of residual error terms, which can significantly deteriorate the restoration process [26, 27]. However, since a TOMBO imaging system can capture multiple observations of the same scene, it is possible to reduce the effects of the error terms significantly by using averaged cross-spectral techniques [25]. The use of global point operations will also minimize the additive noise [1520]. To clarify this point, consider an imaging system with (μ × μ) captured images, the averaged spectral and cross-spectral techniques can be then applied to provide results similar to Equation (8) using spectral estimates instead of the true estimates. In mathematical form, F ( z 1 1 D , z 2 1 D , R ) and similarly H i , j ( z 1 1 D , z 2 1 D , R ) can be estimated using the averaged cross-spectra defined by
i = 1 μ j = 1 μ C ^ G i , j H i , j ( z 1 , z 2 , R ) = F ( z 1 1 D , z 2 1 D , R ) i = 1 μ j = 1 μ C ^ H i . j H i , j ( z 1 1 D , z 2 1 D , R ) + i = 1 μ j = 1 μ C ^ V i , j H i , j ( z 1 1 D , z 2 1 D , R )
where XY*(z1,z2,R) = X(z1,z2, R) Y*(z1, z2, R) is an estimate of the cross-spectra between X(z1, z2, R) and Y(z1,z2, R). For sufficient number of μ images, the last summation (error term of Equation (12)), is nothing but the mean value of an i.i.d. signal which has a zero mean [24, 25]. In this situation and after employing interpolation techniques, we have:
F ^ ( z 1 , z 2 , R ) i = 1 μ j = 1 μ C ^ G i , j H i , j ( z 1 , z 2 , R ) i = 1 μ j = 1 μ C ^ H i , j H i , j ( z 1 , z 2 , R )
Under some constraints [17, 19, 20, 22, 25-28], the above estimates can be used to restore the original image.

3. Color Image Restoration Algorithm

In this section, we employ the work formulated above to develop a color image restoration algorithm. Using Equation (13), it is possible to restore the original image using iterative techniques. During the restoration process, the algorithm will only impose two fundamental image restoration constraints (positivity and support region):
  • For restoring the image
    f ^ ( x , y , R ) = { | f ^ ( x , y , R ) | , < x , y > [ L M , L N ] F M , | f ^ ( x , y , R ) | F , < x , y > [ L M , L N ] 0 otherwise
    where, FM is the mean value of a mesh of pixels in the region < x, y > ∈ [LM, LN] surrounding symmetrically the pixel index (x, y), F∞ is a threshold of a large value pre-defined for the case when the summation of the PSF spectra in Equation 13 is close to zero.
  • For estimating the PSFs
    h ^ i , j ( x , y , R ) = { | h ^ i , j ( x , y , R ) | 0 , < x , y , [ , ] otherwise
where, LD is the interpolation or up-sampling factor needed to restore the HR image, (LM × LN) is the size of the restored image which can be greater or equal to the size of the original image (M × N), and (ℓ × ℓ) is the size of the estimated PSF.
Practical considerations for the iterative algorithm implementation can be summarized as follows:
  • Pixel amplitudes that reach values greater than 255 are scaled using the following histogram normalization,
    f ^ scaled ( x , y , R ) = a × ( f ^ ( x , y , R ) f min R ) ( f max , R f min , R ) + b
    where a and b are usually 255 and 0 respectively (but other values can be also used to adjust contrast and brightness), fmax, R and fmin, R are the maximum and minimum pixel values of the color component R. For improved performance, fmm,R and fmax,R are usually chosen as the 5% and 95% levels in the histogram distribution respectively (confidence interval).
  • The mean value of the input image(s) and the output image is to be maintained (note that there are twice as many green pixels as red/blue pixel for the Bayer filter)
  • To resolve the problem of having zeros or nulls in the spectra, the following equation for the interpolated f(x,y, R) is used:
    F ( z 1 , z 2 , R ) i = 1 μ j = 1 μ C ^ G i , j H i , j ( z 1 , z 2 , R ) α 1 i = 1 μ j = 1 μ C ^ H i , j H i , j ( z 1 , z 2 , R ) α 2
    followed by:
    F ^ ( n + 1 ) ( z 1 , z 2 , R ) = β F ^ ( n ) ( z 1 , z 2 , R ) + ( 1 β ) F ^ ( z 1 , z 2 , R )
    where α1,α2 are two positive numbers representing the extent of additive noise in the residual terms, β is a recursive stability factor controls the amount of information needed from posterior image estimates (0 < β < 1 ), n + 1 is the current iteration number. Typical values of β range between 0.5 to 0.9. Like in adaptive systems, the value of β can be adjusted to avoid such that impulsive-like outputs.
  • For initialization, one of the images is used as an initial estimate of the HR image. The up-sampling and interpolation process is done by zero-padding in the spatial domain between the image samples. Afterwards the FFT is applied. In the Fourier domain, a single spectrum is then taken out of the repetitive spectra using a low pass filter with cut-off frequency ( π L ) and zeroing the rest of the spectrum. Finally, inverse fast Fourier transform (IFFT) is applied to inverse back to the image domain. It is essential that the zero-padding be done such that the zero frequency components remain the same. In addition, zero-padding should be applied to both positive and negative frequencies. Unlike existing techniques that use lower order functions for interpolation (cubic interpolation used in [2]), our method uses the more efficient sinc function.
  • We use the 2-D fast fourier transform (FFT) to estimate spectra and cross spectra needed for the algorithm
The proposed HR color image restoration algorithm is detailed in Table 1.

3.1. Convergence Analysis

From Equation (18) and after considering the situation where n = 0,1,2,… ,no (note that also changes at each iteration), after some mathematical manipulation it can be shown that
F ^ ( n 0 + 1 ) ( z 1 , z 2 , R ) = β ( n 0 + 1 ) F ^ ( 0 ) ( z 1 , z 2 , R ) + j = 0 n 0 β j ( 1 β ) F ( n 0 j ) ( z 1 , z 2 , R )
For a value of 0 < β < 1 and when the number of iteration no → ∞ or large enough, the above equation can be approximated to:
F ^ ( n 0 + 1 ) ( z 1 , z 2 , R ) = j 0 n 0 β j ( 1 β ) F ( n 0 j ) ( z 1 , z 2 , R )
Based on spectral estimation principals [24, 25], and by using the image estimator provided in Equation (17), it can be proven that:
E { F ^ ( z 1 , z 2 , R ) } = F ( z 1 , z 2 , R )
Thus
E { F ^ ( n 0 + 1 ) ( z 1 , z 2 , R ) } = j = 0 n 0 β j ( 1 β ) E { F ( n 0 j ) ( z 1 , z 2 , R ) }
E { F ^ ( n 0 + 1 ) ( z 1 , z 2 , R ) } = ( 1 β ) F ( z 1 , z 2 , R ) j = 0 n 0 β j = F ( z 1 , z 2 , R ) ( 1 β ) ( 1 β ( n 0 + 1 ) ) ( 1 β ) = F ( z 1 , z 2 , R ) ( 1 β ( n 0 + 1 ) ) | n 0 = F ( z 1 , z 2 , R )
From the above analysis, it can be clearly seen that for any value of 0 < β < 1, the algorithm will converge to an unbiased estimator of the original image.

4. Results and Discussion

To evaluate the performance of the proposed color image restoration algorithm, we consider the following tasks:
  • Restore a HR image from multiple blurred, LR and noisy “simulated” TOMBO color images.
  • Restore a HR image from multiple blurred, LR and noisy images and compare the results with the previous method [9, 10].
  • Restore a HR image from multiple blurred, noisy “real” TOMBO color images.

4.1. Examples of Simulated Images

In this section, we test the performance of our proposed algorithm in restoring a HR image from simulated TOMBO images of “Lena” [29]. The algorithm performance is compared with the pixel rearrange method developed in [2] for TOMBO imagers. Simulation parameters for the generated LR, blurred and noisy images are given in Table 2 for noiseless and noisy cases. The simulated images are generated f(ollowing the pro)cedure described in [1]. In all simulations, the SNER is defined as SNER = 10 log ( Signal Energy Noise Enegry ).
Figure 4 shows that our restoration algorithm can restore a HR color image from the simulated LR, blurred TOMBO color images in the absence of additive noise. Our method performs better than the pixel rearrange method because the pixel rearrange method can not align captured images. In Figure 5, we tested the algorithm with the simulation parameters given in Table 2, but with additive noise. Our algorithm appears almost insensitive to additive noise, whereas a significant amount of noise can still be seen in the image restored by of the pixel rearrange method. From the two simulation results, we can see that our proposed algorithm converges rapidly.

4.2. Comparison with Existing Image Restoration Methods

In this section, we compare our algorithm with the advanced restoration methods developed by Sina in [9] and Sroubek in [10]. We generated 12 blurred, LR, noiseless and noisy images following the procedure explained in [9]. Results for the two cases are shown in Figure 6, and the PSNR [9] for each restoration is summarized in Table 3.
At high SNRs, our algorithm does not perform as well as Sina's, which enforces more constraints including (i) a penalty term to enforce similarities between the raw data and the HR estimate (data fidelity penalty term), (ii) a penalty term to encourage sharp edges in the luminance component of the HR image (spatial luminance penalty term), (iii) a penalty term to encourage smoothness in the chrominance component of the HR image (spatial chrominance penalty term), and (iv) A penalty term to encourage homogeneity of the edge location and orientation in different color bands (inter-color dependencies penalty term). In Sroubek's method, the following constraints are enforced: (i) a fidelity term, (ii) a smoothing term using variational integrals, (iii) a consistency term that binds the different volatile PSFs, (iv) a smoothing term to overcome the higher nullity of integer-factors, and (v) an anisotropic term for edge preservation. Although our algorithm only imposes two fundamental constraints, its performance is visually satisfactory, as seen in Figure 6. Our method achieved a PSNR of 16.348 dB in contrast to a PSNR of 21.986 dB and 18.250 dB for Sina's and Sroubek's methods, respectively. A higher PSNR, however, does not always correlate to subjective quality [30, 31]. For instance, the Shift-And-Add restored image (in Sina et al. [9] page 148) has a PSNR of 17.17 dB, but is clearly of poor quality compared to the one we restored. At low SNRs (10.46 dB), our algorithm visually outperforms both aforementioned methods (Figure 6). However, their PSNRs (14.13 dB for Sina's and 13.78 dB for Sroubek's) are higher than the 10.89 dB achieved by our algorithm. This can be explained by the fact that both Sina's and Sroubek's methods minimize regularized energy functions.
Furthermore, we have observed that Sina's method experience some limitations when applied to images blurred using semi-gaussian PSFs. This could be due to the gaussian kernel used in Sina's approach. It is also important to point out that Sina's method uses multiple frames of an image with different displacements (i.e. diversity is guaranteed, see section VI in [9]), while our method uses multiple simultaneous observations of the same scene captured by the TOMBO imager. The diversity between all captured images cannot be guaranteed for the case of our TOMBO imager.

4.3. Examples of Real Images

In this section, we investigate the performance of our proposed algorithm with real captured TOMBO color images. In this example, the captured images are that of a “teddy bear” picture, with a plan object located at 200 mm from the sensor module. Each unit imager has 60 × 60 pixels and each pixel is 6.25 μm × 6.25 μm. The microlens array has the following characteristics: 1.3 mm focal length, 0.5 mm diameter of aperture, and a 0.5 mm pitch for the microlens array. The TOMBO imager integrates a color filter array with a Bayer pattern. Demosaicing is achieved on-chip. To test the performance of the restoration algorithms at lower SNRs, a zero mean white Gaussian noise is added manually to the captured LR images. System parameters for this real example are given in Table 4.
Figure 7 demonstrates that our proposed algorithm is able to restore a HR color image of the “teddy bear” picture using real captured TOMBO LR observations. In addition, our algorithm can significantly minimize the additive noise. It outperforms the conventional pixel rearrange method [33, 34]. Furthermore, our algorithm can restore the HR color image within only 15 iterations.

5. Conclusions

A blind color image restoration method is proposed for the reconstruction of HR color images using multiple LR, degraded and noisy color images captured by TOMBO imagers. The proposed spectral-based method only imposes two fundamental image restoration constraints (positivity and support region) to (i) correct the blur that affects the captured LR images, (ii) minimize the interchannel cross-correlations between RGB color components, (iii) significantly reduce the impact of additive noise, and (iv) reconstruct a HR color image. The computation complexity of the algorithm is low compared with existing techniques because it uses FFT for spectral estimation.
The proposed restoration algorithm has a rapid convergence rate of 10 to 20 iterations. Results show that the proposed algorithm is capable of restoring a HR image from the degraded LR color TOMBO images even when the SNER is as low as 5 dB. The proposed algorithm uses FFT and only two fundamental image restoration constraints, which makes it suitable for silicon integration with the TOMBO imager.

Acknowledgments

This work is supported by the Australian Research Council's Discovery Project DP0664909. The authors would like to thank Prof. Tanida and Dr. Ryoichi Horisaki for providing us with images of their experimental TOMBO imager. Thanks also go to Prof. Peyman Milanfar and Prof. Sina Farsiu for providing us with the lighthouse images and the MDSP resolution enhancement software; Dr Sroubek Filip and Dr Barbara Zitova for providing us with the BSR superresolution and blind deconvolution GUI. The authors would also like to thank Ms. Sabine Betts from Microelectronics Research Group (MRG) for her support.

References and Notes

  1. El-Sallam, A.; Boussaid, F. Spectral-based blind image restoration method for thin TOMBO imagers. Sensors 2008, 8, 6108–6124. [Google Scholar]
  2. Tanida, J.; Kumagai, T.; Yamada, K.; Miyatake, S.; Ishida, K.; Morimoto, T.; Kondou, N.; Miyazaki, D.; Ichioka, Y. Thin observation module by bound optics (TOMBO): Concept and experimental verification. Appl. Opt. 2001, 40, 1806–1813. [Google Scholar]
  3. Yamada, K.; Tanida, J.; Yu, W.; Miyatake, S.; Ishida, K.; Miyazaki, D. Fabrication of diffractive microlens array for opto-electronic hybrid information system. Proc. Diffract. Opt. 1999, 22, 52–53. [Google Scholar]
  4. Tanida, J.; Kumagai, T.; Yamada, K.; Miyatake, S.; Ishida, K.; Morimoto, T.; Kondou, N.; Miyazaki, D.; Ichioka, Y. Thin observation module by bound optics-TOMBO: an optoelectronic image capturing system. Proc. SPIE Opt. Comput. 2000, 4086, 1030–1036. [Google Scholar]
  5. Tanida, J.; Yamada, K. TOMBO: thin observation module by bound optics. Proceedings of the 15th Annual Meeting of the IEEE in Lasers and Electro-Optics, Glasgow, Scotland, November 10-14, 2002; 1, pp. 233–234.
  6. Kitamura, Y.; Shogenji, R.; Yamada, K.; Miyatake, S.; Miyamoto, M.; Morimoto, T.; Masaki, Y.; Kondou, N.; Miyazaki, D.; Tanida, J.; Ichioka, Y. Reconstruction of a high-resolution image on a compound-eye image-capturing system. Appl. Opt. 2004, 43, 1719–1727. [Google Scholar]
  7. Nitta, K.; Shogenji, R.; Miyatake, S.; Tanida, J. Image reconstruction for thin observation module by bound optics by using the iterative backprojection method. Appl. Opt. 2006, 45, 2893–2900. [Google Scholar]
  8. Yamada, K.; Ishida, K.; Shougenji, R.; Tanida, J. Development of three dimensional endoscope by Thin Observation by Bound Optics (TOMBO). Proceedings of World Automation Congress WAC'06, Budapest, Hungary, July 24-26, 2006; pp. 1–4.
  9. Farsiu, S.; Elad, M.; Milanfar, P. Multiframe demosaicing and super-resolution of color images. IEEE Trans. Image Process. 2006, 15, 141–159. [Google Scholar]
  10. Filip, S.; Flusser, J. Multichannel blind deconvolution of spatially misaligned images. IEEE Trans. Image Process. 2005, 14, 874–883. [Google Scholar]
  11. Yu, H.; Yap, K.; Li, C; Chau, L. A new color image regularization scheme for blind image decon-volution. Proceedings of the IEEE ICASSP'08, Las Vegas, Nevada, USA, March 30 - April 4, 2008.
  12. He, Y.; Yap, K.H.; Chen, L.; Chau, L.P. A novel hybrid model framework to blind color image deconvolution. IEEE Trans. Syst. Man Cybern. Part A: Syst. Humans 2008, 38, 862–880. [Google Scholar]
  13. Vega, M.; Molina, R.; Katsaggelos, A.K. A Bayesian super-resolution approach to demosaicing of blurred images. EURASIP 2006, 25072:1–25072:12. [Google Scholar]
  14. Ohta, Y.I.; Kanade, T.; Sakai, T. Color information for region segmentation. J. Comput. Graph. Image Process. 1980, 13, 222–241. [Google Scholar]
  15. Gonzalez, R.C.; Woods, R.E. Digital Image Processing, 3rd ed.; Prentice Hall: Englewood Cliffs, NJ, USA, August 2007. [Google Scholar]
  16. Pratt, W.K. Digital Image Processing, 4th ed.; Wiley: Hoboken, NJ, USA, 2007. [Google Scholar]
  17. Lagendijk, R.L.; Biemond, J. Iterative Identification and Restoration of Images; Springer-Verlag Inc.: New York, NY, USA, 2007. [Google Scholar]
  18. Jain, A.K. Fundamentals of Digital Image Processing; Prentice-Hall: Englewood Cliffs, NJ, USA, 1989. [Google Scholar]
  19. Banham, M.R.; Katsaggelos, A.K. Digital image restoration. IEEE Signal Process. Mag. 1997, 14, 24–41. [Google Scholar]
  20. Zhao, W.; Pope, A. Image restoration under significant additive noise. IEEE Signal Process. Lett. 2007, 14, 401–404. [Google Scholar]
  21. Ohyama, N.; Yachida, M.; Badique, E.; Tsujiuchi, J.; Honda, T. Least squares filter for color image restoration. J. Opt. Soc. Am. A 1988, 5, 19–24. [Google Scholar]
  22. Zing, X.Y.; Chen, Y.W.; Nakao, Z. Classification of remedy sensed images using independent component analysis and spatial consistency. J. Adv. Comput. Intell. Intell. Inform. 2004, 8, 216–217. [Google Scholar]
  23. Boo, K.J.; Bose, N.K. Multispectral image restoration with multisensors. Proceedings of International Conference on Image Processing, Lausanne, Switzerland, September 16-19, 1996; 3, pp. 995–998.
  24. Brillinger, D.R. Time Series: Data Analysis and Theory; Holden-Day, Inc.: San Francisco, CA, USA, 1981. [Google Scholar]
  25. Kay, S.M. Modern Spectral Estimation: Theory and Application; Prentice Hall: Englewood Cliffs, NJ, USA, 1988. [Google Scholar]
  26. Choi, K.; Schulz, T.J. Signal-processing approaches for image-resolution restoration for TOMBO imagery. Appl. Opt. 2008, 47, B104–B116. [Google Scholar]
  27. Kundur, D.; Hatzinakos, D. Blind image deconvolution revisited. IEEE Signal Process. Mag. 1996, 13, 61–63. [Google Scholar]
  28. Chen, L.; Yap, K.; He, Y. Efficient recursive multichannel blind image restoration. EURASIP J. Adv. Signal Process. 2007, 2007, 1–10. [Google Scholar]
  29. Munson, D.C. A note on lena. IEEE Trans. Image Process. 1996, 5, 1–2. [Google Scholar]
  30. Huynh-Thu, Q.; Ghanbari, M. Scope of validity of PSNR in image/video quality assessment. Electron. Lett. 2008, 44, 800–801. [Google Scholar]
  31. Girod, A. What's wrong with mean-squared error? In Digital images and human vision; MIT Press: Cambridge, MA, USA, 1993; pp. 207–220. [Google Scholar]
  32. Shogenji, R.; Kitamura, Y.; Yamada, K.; Miyatake, S.; Tanida, J. Color imaging with an integrated compound imaging system. Opt. Expr. Opt. Soc. Am. 2003, 11, 2109–2117. [Google Scholar]
  33. Tanida, J.; Shogenji, R.; Kitamura, Y.; Yamada, K.; Miyamoto, M.; Miyatake, S. Multispectral imaging using compact compound optics. Opt. Expr. Opt. Soc. Am. 2004, 12, 1643–1655. [Google Scholar]
  34. Kanaev, A.V.; Ackerman, J.R.; Fleet, E.F.; Scibner, D.A. TOMBO sensors with scene-independent superresolution processing. Opt. lett. 2007, 32, 2855–2857. [Google Scholar]
Figure 1. The architecture of a color TOMBO imaging system.
Figure 1. The architecture of a color TOMBO imaging system.
Sensors 09 04649f1
Figure 2. Point operations categories.
Figure 2. Point operations categories.
Sensors 09 04649f2
Figure 3. System model for the color TOMBO system.
Figure 3. System model for the color TOMBO system.
Sensors 09 04649f3
Figure 4. Simulation results, 6 × 6 images, D = 4, L = 4, ℓ = 7, no additive noise.
Figure 4. Simulation results, 6 × 6 images, D = 4, L = 4, ℓ = 7, no additive noise.
Sensors 09 04649f4
Figure 5. Simulation results, 6 × 6 images, D = 4, L = 4, ℓ = 7, SNER = 4.968 dB.
Figure 5. Simulation results, 6 × 6 images, D = 4, L = 4, ℓ = 7, SNER = 4.968 dB.
Sensors 09 04649f5
Figure 6. Simulation results for 12 blurred LR images of the lighthouse, in noiseless, and noisy, SNER=10.4616 dB, D = 4, L = 4, ℓ = 5.
Figure 6. Simulation results for 12 blurred LR images of the lighthouse, in noiseless, and noisy, SNER=10.4616 dB, D = 4, L = 4, ℓ = 5.
Sensors 09 04649f6
Figure 7. Experimental results, 4 × 4 unit images, D = 4, L = 4, ℓ = 5, SNER = 15.774 dB.
Figure 7. Experimental results, 4 × 4 unit images, D = 4, L = 4, ℓ = 5, SNER = 15.774 dB.
Sensors 09 04649f7
Table 1. Proposed color image restoration algorithm.
Table 1. Proposed color image restoration algorithm.
Step 1: Set the values of L, ℓ, α1, α2,β, FM, F
Step 2: Select the color to be restored ϑ ε {R, G, B}
Step 3: Iteration n = 1, initialize (x, y, ϑ), hence (f1, f2,ϑ) = FFT {(x, y, ϑ)}
Step4:For i,j = 1,2, …, μ estimate
C ^ G i , j F ( f 1 , f 2 , ϑ ) = G i , j ( f 1 , f 2 , ϑ ) F ^ ( f 1 , f 2 , ϑ )
H i , j ( f 1 , f 2 , ϑ ) = C ^ G i , j F ( f 1 , f 2 , ϑ ) α 1 F ^ ( f 1 , f 2 , ϑ ) F ^ ( f 1 , f 2 , ϑ ) + α 2 h i , j ( x , y , ϑ ) = IFFT { H i , j ( f 1 , f 2 , ϑ ) }
Step 5: Impose PSF constraints to get the accurate estimates
h ^ i , j ( x , y , ϑ ) = { | h i , j ( x , y , ϑ ) | , 0 < x , y > [ , ] otherwise H ^ i , j ( f 1 , f 2 , ϑ ) = FFT { h ^ i , j ( x , y , ϑ ) }
Step 6: Estimate the biased image spectra
F ( f 1 , f 2 , ϑ ) = i = 1 μ j = 1 μ C ^ G i , j H i , j ( f 1 , f 2 , ϑ ) α 1 i = 1 μ j = 1 μ C ^ H i , j H i , j ( f 1 , f 2 , ϑ ) + α 2 f ( x , y , ϑ ) = IFFT { F ( f 1 , f 2 , ϑ ) }
Step 7: Impose the image constraints
f ( x , y , ϑ ) = { | f ( x , y , ϑ ) | < x , y > [ L M , L N ] F M , | f ( x , y , ϑ ) | F , < x , y > [ L M , L N ] 0 otherwise
then estimate the original image by updating the estimates using
f ^ ( x , y , ϑ ) = β f ^ ( x , y , ϑ ) + ( 1 β ) f ( x , y , ϑ )
Step 8: Scale the estimated images pixels or use histogram normalization to find a and b, then adjust the image using
f ^ ( x , y , ϑ ) = a × ( f ^ ( x , y , ϑ ) f ^ min , ϑ ) ( f ^ max , ϑ f ^ min , ϑ ) + b
Step 9: Repeat from Step 4 until convergence, then repeat for another color ϑ
Table 2. Input parameters for simulated TOMBO images.
Table 2. Input parameters for simulated TOMBO images.
μ × μM × NSNERLM ×LNα1α2β# Iterations
Figure 46 × 660 × 60-7240 × 2400.1100.120
Figure 56 × 660 × 604.968 dB7240 × 2400.1100.120
Table 3. PSNR values (in dB) for the restoration methods.
Table 3. PSNR values (in dB) for the restoration methods.
Sina [9]Sroubek [10]This Work
Noiseless21.98618.25016.348
Noisy14.1313.7810.89
Table 4. Input parameters for real captured TOMBO images, D = 4.
Table 4. Input parameters for real captured TOMBO images, D = 4.
μ × μM × NSNER↑LLM× LNα11α2β# of Iterations
Figure 74 ×460 × 015.774 dB54240400.0010.0010.925

Share and Cite

MDPI and ACS Style

El-Sallam, A.A.; Boussaid, F. A High Resolution Color Image Restoration Algorithm for Thin TOMBO Imaging Systems. Sensors 2009, 9, 4649-4668. https://doi.org/10.3390/s90604649

AMA Style

El-Sallam AA, Boussaid F. A High Resolution Color Image Restoration Algorithm for Thin TOMBO Imaging Systems. Sensors. 2009; 9(6):4649-4668. https://doi.org/10.3390/s90604649

Chicago/Turabian Style

El-Sallam, Amar A., and Farid Boussaid. 2009. "A High Resolution Color Image Restoration Algorithm for Thin TOMBO Imaging Systems" Sensors 9, no. 6: 4649-4668. https://doi.org/10.3390/s90604649

APA Style

El-Sallam, A. A., & Boussaid, F. (2009). A High Resolution Color Image Restoration Algorithm for Thin TOMBO Imaging Systems. Sensors, 9(6), 4649-4668. https://doi.org/10.3390/s90604649

Article Metrics

Back to TopTop