Table of Contents Author Guidelines Submit a Manuscript
Mathematical Problems in Engineering
Volume 2011 (2011), Article ID 380807, 20 pages
http://dx.doi.org/10.1155/2011/380807
Research Article

Efficient and Effective Total Variation Image Super-Resolution: A Preconditioned Operator Splitting Approach

1School of Computer Science and Technology, Nanjing University of Science and Technology, Nanjing 210094, China
2Department of Information and Computing Science, Guangxi University of Technology, Liuzhou 545006, China
3Department of Applied Mathematics, Nanjing University of Science and Technology, Nanjing 210094, China

Received 23 August 2010; Revised 30 November 2010; Accepted 4 January 2011

Academic Editor: J. J. Judice

Copyright © 2011 Li-Li Huang et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

Super-resolution is a fusion process for reconstructing a high-resolution image from a set of low-resolution images. This paper proposes a novel approach to image super-resolution based on total variation (TV) regularization. We applied the Douglas-Rachford splitting technique to the constrained TV-based variational SR model which is separated into three subproblems that are easy to solve. Then, we derive an efficient and effective iterative scheme, which includes a fast iterative shrinkage/thresholding algorithm for denoising problem, a very simple noniterative algorithm for fusion part, and linear equation systems for deblurring process. Moreover, to speed up convergence, we provide an accelerated scheme based on precondition design of initial guess and forward-backward splitting technique which yields linear systems of equations with a nice structure. The proposed algorithm shares a remarkable simplicity together with a proven global rate of convergence which is significantly better than currently known lagged diffusivity fixed point iteration algorithm and fast decoupling algorithm by exploiting the alternating minimizing approach. Experimental results are presented to illustrate the effectiveness of the proposed algorithm.

1. Introduction

Multiframe image super-resolution (SR) is one of the promising techniques in image processing community since it enables us to obtain an image with a resolution that exceeds the hardware limitation, for example, the number of pixels in a charge-coupled device (CCD). Super-resolution is the process of combining a sequence of low-resolution (LR) noisy blurred images to produce a high-resolution (HR) image of sequence. It overcomes the inherent resolution limitation by bringing together the additional information from each image.

Generally, SR techniques can be divided into two broad categories: frequency domain methods and spatial domain methods. Most of the earlier SR work was developed in frequency domain using discrete Fourier transform (DFT), such as the work of Tsai and Huang [1], Kim et al. [2, 3], where high-frequency information is extracted from low-frequency data in the given LR frames. Many other popular frequency domain methods were proposed in Discrete cosine transform (DCT) domain and wavelet domain [47]. Although the frequency domain methods are intuitively simple and computationally cheap, they are extremely sensitive to model error [8], limiting their use. In a spatial domain alternative, Patti et al. [9] and Stark and Oskoui [10] proposed projection onto convex sets (POCSs) algorithm. A related method, the iterative back projection, was developed in [11, 12]. Although projection-based algorithms are usually robust to noise and allow some modeling flexibility, they are also known for their low rate of convergence. The hybrid maximum likelihood (ML)-maximum a posteriori (MAP)-POCS method was proposed in [13]. Based on these basic reconstruction methods, researchers have produced many extended algorithms, such as nonlocal-means (NLMs) based approach [14], multidimensional kernel regression-based approach [15], the joint formulation of reconstruction and registration [1619], algorithms for multi-spectral and color [20, 21], hyperspectral [22], and compressed [23, 24] sequence.

The spatial domain methods discussed so far are generally confronted with the problem of slow convergence and expensive computation. To apply the SR algorithm to practical situations, many novel and powerful algebraic techniques have been proposed to reduce the computation complexity. For example, authors of [2529] proposed efficient preconditioners to accelerate convergence of a conjugate gradient minimization algorithm. Xiao et al. [30] proposed an efficiency SR reconstruction algorithm employing the Armijo rule to identify the step length instead of the exact line search and replaced numerical approximation of the gradient of the MAP object function by analytic approximation. For pure translational motion and common space invariant blurring model, Elad and Hel-Or [31] proposed a novel SR algorithm that separates fusion and deblurring. The fusion method is a very simple noniterative algorithm, while preserving its optimality in ML sense. Farsiu et al. [32] proposed an efficient two-stage method for minimizing a novel framework combining a robust 𝑙1 norm fidelity term and a bilateral prior, leading to an initial Median Shift-And-Add operation on Bayer-filtered LR data followed by a deblurring and interpolating stage. Huang et al. [33] proposed a fast decoupling algorithm by exploiting the alternating minimization approach.

In this paper, we propose a general framework for multiple shifted and linear space-invariant noisy blurred LR image frames which subsume several well-known SR models. The proposed model combines the TV regularization to formulate the SR image reconstruction as an optimization problem.

The purpose of this paper is to study an efficient TV-based SR reconstruction algorithm. There are two major contributions in this paper. As the first contribution, we propose an efficient algorithm that takes full advantage of the problem structures; that is, geometrical motion matrices, blur matrix, and the first-order finite-difference matrix all have block-circulant-circulant-block (BCCB) structure under periodic boundary condition. As such, we propose to compute the minimizer of our SR model by applying Douglas-Rachford splitting (DBS) techniques, respectively, alternating direction methods of multipliers (ADMM), which separate the SR model into three subproblems that are easy to solve. As the second contribution, we provide an accelerated scheme based on precondition design of initial value and forward-backward splitting (FBS) to speed up convergence. Our method can separate the SR treatment into measurements fusion, denoising, and deblurring. The fusion part is shown to be a very simple noniterative algorithm. The denoising problem can be solved by a linear time shrinkage operation. Fast Fourier transform (FFT) is employed to solve the deblurring problem. Finally, experimental results are presented to illustrate the effectiveness of the proposed algorithm.

The outline of this paper is as follows. In Section 2, we present the image observation model of the SR problem, then propose a TV-based SR model. In Section 3, we present an efficient SR reconstruction algorithm. Experimental results are provided in Section 4. Finally, concluding remarks are given in Section 5.

2. Problem Formulation

2.1. Observation Model

The image observation model is employed to relate the desired referenced HR image to all the observed LR images. Consider the desired HR image of size 𝐿1𝑁1×𝐿2𝑁2 written in lexicographical notation as the vector 𝑧=[𝑧1,𝑧2,,𝑧𝑁]𝑇, where 𝑁=𝐿1𝑁1×𝐿2𝑁2. Let the parameters 𝐿1 and 𝐿2 be the subsampling factors in the horizontal and vertical directions, respectively; each observed LR image has the size 𝑁1×𝑁2. Thus, the LR image can be represented as 𝑦𝑘=[𝑦𝑘1,𝑦𝑘2,,𝑦𝑘𝑀]𝑇, for 𝑘=1,2,,𝐾 and 𝑀=𝑁1×𝑁2. A popular model assumes that LR images {𝑦𝑘}𝐾𝑘=1 are generated from HR image 𝑧 through a sequence of operations that includes (i) geometrical motions 𝑀𝑘, (ii) a linear space-invariant blur 𝐵, (iii) a subsampling step represented by 𝑆, and finally (iv) an additive white Gaussian noise 𝑛𝑘 with zero mean that represents both measurements noise and model mismatch [32]. All these are linear operators, represented by a matrix multiplying the image they operate on. We assume hereafter that 𝐵 and 𝑆 are identical for all images in the sequence. This model leads to the following set of equations, where all images are ordered lexicographically:𝑦𝑘=𝑆𝐵𝑀𝑘𝑧+𝑛𝑘=𝑊𝑘𝑧+𝑛𝑘for𝑘=1,2,,𝐾,(2.1) where 𝑊𝑘=𝑆𝐵𝑀𝑘 represents the imaging system.

The recovery of 𝑧 from {𝑦𝑘}𝐾𝑘=1 is thus an inverse problem, combining motion compensation, denoising, deblurring, scaling-up operation, and fusion of the different images, which all merged to one. The quality of the desired SR image depends on the assumption that 𝑆,𝐵, and 𝑀𝑘 are known, or the accuracy in estimating the degraded operators. Throughout this paper, we assume that 𝑆,𝐵, and 𝑀𝑘 are known. The decimation 𝑆 is dependent on the resolution scale factor that we aim to achieve, and, as such, it is easily fixed. In this work, we shall assume that this resolution factor is an integer 𝑠1 on both axes. In most cases, the blur 𝐵 refers to the camera point spread function (PSF), and, therefore, it is also accessible. Even if this is not the case, the blurring is typically dependent on few parameters, and those, in the worst case, can be manually set. To be identical with the work of Elad and Farsiu [31, 32], we focus on the simplest of the motion models, namely, the translational model. Reference [34] detailed the several reasons for this. We believe that an in-depth study of this simple case allows much insight to be gained about the problems inherent to SR image reconstruction.

2.2. TV Regularization-Based SR Model

In general, SR is an ill-posed problem either because the information contained in the observed LR images is not sufficient or because it has great sensitivity to the noise. Procedures adopted to stabilize the inversion of ill-posed problem are called regularization. In such stabilization scheme, we reconstruct the original HR image by finding the minimizer of some appropriate functional𝐸(𝑧)=Φ(𝑧)+𝜇Ψ(𝑧),(2.2) where Φ(𝑧) is a regularization term which includes prior information about the original image, Ψ(𝑧) denotes the data fitting term depending on the given LR image {𝑦𝑘}𝐾𝑘=1, and 𝜇 is a positive parameter which controls the tradeoff between the two terms for minimization.

In general, the data fitting term can be deduced from an MAP estimation. If 𝑧 is corrupted by additive white Gaussian noise, the MAP estimation will lead to the data fitting term Ψ(𝑧)=(1/2)𝐾𝑘=1𝑊𝑘𝑧𝑦𝑘2𝐿2(Ω𝐿), where Ω𝐿Ω and Ω denote a bounded and open domain of continuous LR image and HR image in 𝑅2, respectively. For regularization term, the popular choice is total variation seminorm 𝑧TV=Ω|𝑧|d𝑥, which was first proposed for image denoising [35], because TV norm can better preserve sharp edges or object boundaries that are usually the most important features to recover.

As stated previously, to invert the degradation process in (2.1), we can formulate a TV regularization model which requires solving the variational problem:min𝑧BV(Ω)𝐸(𝑧)=𝑧TV+𝜇2𝐾𝑘=1𝑊𝑘𝑧𝑦𝑘2𝐿2Ω𝐿,(2.3) where BV(Ω) is the space of functions of bounded variation. Note that the fitting term in (2.3) is strictly convex and coercive and the TV regularity term is also convex (though not strictly so) and lower semicontinuous. So the objective function 𝐸(𝑧) is globally strictly convex and possesses a unique minimizer. In terms of optimization, these are desirable properties.

3. An Operator Splitting Approach to TV-Based Super-Resolution

To solve the desired HR image of (2.3), commonly used method is the gradient descent method [1932, 36]. Although this approach is simple, the nonlinearly and poor conditioning of the problem make this approach very slow. A more efficient class of solvers are those based on a linearized gradient method which solves the associated Euler-Lagrange equation via a lagged diffusivity fixed-point iteration [28, 29]. In each iteration of the linearized gradient method, a linear system needs to be solved, which becomes more and more difficult as 𝐵 becomes more ill-conditioned. Another group of algorithms is based on the well-known variable-splitting and penalty techniques in optimization. These ideas have gained wide application in image processing, such as works in [3741]. Recently, Huang et al. [33] modified the SR model (2.3) by adding a quadratic term to get a simpler alternating minimization algorithm. The drawback of this method is the same as the lagged diffusivity fixed point method.

To efficiently solve the SR problem (2.3), in this section, we will show that the operator splitting method can be used to divide the problem (2.3) into subproblems that can be solved in sequence, and each of them permits a closed form solution. Among the current splitting methods, the most prominent splitting schemes are forward-backward splitting, double-backward splitting, Peaceman-Rachford splitting, and Douglas-Rachford splitting. In this paper, we will focus on the forward-backward splitting and Douglas-Rachford splitting. One may refer to [42, 43] and the references therein for more details.

3.1. FBS and DRS

Let 𝐻 be a real Hilbert space, and let 𝐴,𝐵𝐻2𝐻 be two set-valued operators. We assume that 𝐴 and 𝐵 are maximal monotone; that is, their resolvents 𝐽𝐴=(𝐼+𝐴)1 and 𝐽𝐵=(𝐼+𝐵)1 exist and are firmly nonexpansive. The problem which we will describe as a fundamental problem can be written in the form of a common zero inclusion problem0𝐴(̂𝑥)+𝐵(̂𝑥).(3.1)

The idea of the forward-backward splitting algorithm is that, for any constant 𝛾>0, we have0𝐴(̂𝑥)+𝐵(̂𝑥)̂𝑥𝛾𝐵(̂𝑥)̂𝑥+𝛾𝐴(̂𝑥)̂𝑥𝐽𝛾𝐴(𝐼𝛾𝐵)(̂𝑥).(3.2) This leads to the following result.

Theorem 3.1 (FBS [43, Theorem 2.3.17]). Suppose that 𝐴𝐻2𝐻 is maximal monotone and 𝐵𝐻𝐻 is a monotone operator such that 𝜂𝐵 is firmly nonexpansive for some 𝜂>0. Furthermore, assume that a solution of (3.1) exists. Then, for every start element 𝑥0 and step size 𝛾(0,2𝜂), the forward-backward splitting algorithm 𝑥𝑛+1𝐽𝛾𝐴(𝐼𝛾𝐵)(𝑥𝑛)(3.3) converges weakly to an element of the set of solutions (𝐴+𝐵)1(0).

We now describe the Douglas-Rachford splitting scheme, which does exhibit general convergence, at least when used with a constant step size in finite-dimensional spaces. To introduce it, we can rewrite the fixed point relation (3.2) as follows:0𝐴(̂𝑥)+𝐵(̂𝑥)̂𝑥𝐽𝛾𝐴(𝐼𝛾𝐵)(̂𝑥)̂𝑥+𝛾𝐵(̂𝑥)𝐽𝛾𝐴(𝐼𝛾𝐵)(̂𝑥)+𝛾𝐵(̂𝑥)̂𝑥𝐽𝛾𝐵𝐽𝛾𝐴(𝐼𝛾𝐵)(̂𝑥)+𝛾𝐵(̂𝑥).(3.4) This leads to the following result.

Theorem 3.2 (DBS [43, Theorem 2.3.21]). Let 𝐴,𝐵𝐻2𝐻 be maximal monotone, and assume that a solution of (3.1) exists. Then, for any elements 𝑦0 and 𝑥0 and any step size 𝛾>0, the following Douglas-Rachford splitting algorithm converges weakly to an element ̂𝑦: 𝑦𝑛+1=𝐽𝛾A(2𝑥𝑛𝑦𝑛)+𝑦𝑛𝑥𝑛,𝑥𝑛+1=𝐽𝛾𝐵𝑦𝑛+1.(3.5) Furthermore, it holds that ̂𝑥=𝐽𝛾𝐵(̂𝑦) satisfies 0𝐴(̂𝑥)+𝐵(̂𝑥). If 𝐻 is finite-dimensional, then the sequence {𝑥𝑛} converges to ̂𝑥.

3.2. Proposed Algorithm

In this subsection, we will apply the DRS to dual problem of the minimization functional (2.3). We first rewrite the energy functional (2.3) in a discrete formmin𝑧𝑅𝐿1𝑁1×𝐿2𝑁2𝐸(𝑧)=𝜙(𝐷𝑧)Φ(𝑧)+𝜇Ψ(𝑧)(3.6) with Ψ(𝑧)=(1/2)𝐾𝑘=1𝑊𝑘𝑧𝑦𝑘2, 𝜙(𝐷𝑧)=𝑖,𝑗(𝐷𝑧)𝑖,𝑗 is the discrete total variation of 𝑧. Here, denotes Euclidean norm, and 𝐷 is given by (𝐷𝑧)𝑖,𝑗=𝐷+𝑥𝑧𝑖,𝑗;𝐷+𝑦𝑧𝑖,𝑗=𝑧𝑖,𝑗+1𝑧𝑖,𝑗;𝑧𝑖+1,𝑗𝑧𝑖,𝑗,for1𝑖𝐿1𝑁1,1𝑗𝐿2𝑁2,(3.7) where 𝐷+𝑥 and 𝐷+𝑦 are forward difference operators with periodic boundary condition 𝑧𝑖,𝐿2𝑁2+1=𝑧𝑖,1 and 𝑧𝐿1𝑁1+1,𝑗=𝑧1,𝑗. Therefore, 𝐷 is a BCCB matrix.

Consider the equivalent problem of (3.6)min𝑧𝑅𝐿1𝑁1×𝐿2𝑁2𝑤𝑅2𝐿1𝑁1×𝐿2𝑁2{𝐸(𝑧)=𝜙(𝑤)+𝜇Ψ(𝑧)},s.t.𝑤=𝐷𝑧.(3.8) The Lagrangian for problem (3.8) is𝐿(𝑤,𝑧,𝜆)=𝜙(𝑤)+𝜇Ψ(𝑧)𝜆,𝑤𝐷𝑧,(3.9) where the dual variable 𝜆𝑅2𝐿1𝑁1×𝐿2𝑁2 can be thought of as a vector of Lagrange multipliers. Therefore, the dual problem of (3.8) ismax𝜆inf𝑤,𝑧𝐿(𝑤,𝑧,𝜆)=min𝜆𝜙(𝜆)+𝜇Ψ𝜇1𝐷𝜆,(3.10) where 𝜙(resp., Ψ) denotes conjugate function of 𝜙 (resp., Ψ).

Define operators 𝐴 and 𝐵 by𝐴(𝜆)=𝜕𝜙(𝜆),𝐵(𝜆)=𝜇𝜕Ψ𝜇1𝐷𝜆.(3.11) By Fermat's rule, solving the dual problem (3.10) is equivalent to finding 𝜆 such that0𝐴(𝜆)+𝐵(𝜆).(3.12)

By formally applying DRS (3.5) to (3.12) with 𝛼 as the step size, respectively, the ADMM iterations [4451] are given by 𝑤𝑛+1=argmin𝑤𝜙(𝑤)𝜆𝑛,𝑤+𝛼2𝑤𝐷𝑧𝑛2,(3.13)𝑧𝑛+1=argmin𝑧𝜇Ψ(𝑧)+𝜆𝑛,𝐷𝑧+𝛼2𝑤𝑛+1𝐷𝑧2,(3.14)𝜆𝑛+1=𝜆𝑛+𝛼𝐷𝑧𝑛+1𝑤𝑛+1.(3.15) As pointed out by Setzer in [43, 49], we note that this iteration scheme coincides with DRS algorithm (3.5) with 𝑦𝑛=𝛾(𝛼𝜆𝑛+𝐷𝑧𝑛), 𝑥𝑛=𝛾𝛼𝜆𝑛, and 𝛾=𝛼. Since operators 𝜙 and Ψ are proper, lower semicontinuous, and convex, the operators 𝐴 and 𝐵 are maximal monotone see [43]. According to Theorem 3.2, the sequence {𝑧𝑛} converges to solution of (3.6). We also note that the above ADMM algorithm coincides with the alternating split Bregman algorithm proposed by Goldstein and Osher [40] with 𝜕𝜙(𝑤𝑛+1)=𝜆𝑛+1 and 𝜕Ψ(𝑧𝑛+1)=(1/𝜇)𝜆𝑛+1.

3.2.1. 𝑤-Subproblem

It is not difficult to show that the minimization of (3.13) with respect to 𝑤 is equivalent to solving 𝐿1𝑁1×𝐿2𝑁2 two-dimensional problem of the formmin𝑤𝑖,𝑗𝑤𝑖,𝑗+𝛼2𝑤𝑖,𝑗(𝐷𝑧𝑛)𝑖,𝑗+𝜆𝑛𝑖,𝑗𝛼2,(3.16) for which the unique minimizer is given by the following two-dimensional shrinkage formula:𝑤𝑛+1𝑖,𝑗=max(𝐷𝑧𝑛)𝑖,𝑗+𝜆𝑛𝑖,𝑗𝛼1𝛼,0(𝐷𝑧𝑛)𝑖,𝑗+𝜆𝑛𝑖,𝑗/𝛼(𝐷𝑧𝑛)𝑖,𝑗+𝜆𝑛𝑖,𝑗/𝛼,(3.17) where the convention 0(0/0)=0 is followed. Here, shrinkage formula (3.17) which serves as nonlinear low-pass filter to restored HR image is tantamount to denoising treatment.

3.2.2. 𝑧-Subproblem

Subproblem (3.14) is quadratic in 𝑧, and the minimizer 𝑧𝑛+1 is given by the normal equations𝐷𝑇𝐷+𝜇𝛼𝑘𝑊𝑇𝑘𝑊𝑘𝑧𝑛+1=𝐷𝑇𝑤𝑛+1𝜆𝑛𝛼+𝜇𝛼𝑘𝑊𝑇𝑘𝑦𝑘,(3.18) where 𝑊𝑇𝑘 is the transpose operator of 𝑊𝑘.

Note that 𝐷𝑇𝐷 in (3.18) is BCCB matrix and can be diagonalized by FFT. Moreover, with periodic boundary condition, both 𝐵 and 𝑀𝑘 are BCCB matrices. Exploiting the fact that the product order of two BCCB matrices can commute, we get that 𝑀𝑇𝑘𝐵𝑇=𝐵𝑇𝑀𝑇𝑘, 𝐵𝑀𝑘=𝑀𝑘𝐵. Then, 𝑘𝑊𝑇𝑘𝑊𝑘=𝐵𝑇𝑅𝐵 (𝑅=𝑘𝑀𝑘𝑇𝑆𝑇𝑆𝑀𝑘 is a diagonal matrix [31]). However, 𝑘𝑊𝑇𝑘𝑊𝑘=𝐵𝑇𝑅𝐵 does not have BCCB structure. Therefore, it does not allow us to apply FFT implementation to (3.18) directly as done in [37, 38]. The quadratic term Ψ(𝑧)=(1/2)𝐾𝑘=1𝑊𝑘𝑧𝑦𝑘2 in (3.14) that couples the variable 𝑧 by the matrix 𝑊𝑘 makes the algorithm computationally expensive. There are some techniques to overcome these problems. In [50], the authors introduced the three constrains: 𝑤1=𝐵𝑧, 𝑤2=𝐷𝑧, 𝑤3=𝑧 and used the alternating split Bregman technique to maximally decouple the variables. In [51], the linear system was solved noniteratively by using Sherman-Morrison-Woodbury (SMW) matrix inversion formula and FFT to diagonalize the Hessian matrix of the energy functional. In this paper, we use the forward-backward splitting method (3.3) to efficiently solve the 𝑧-subproblem (3.14), which can decouple the variable 𝑧 and the constraint matrix 𝑊𝑘 and make the Hessian matrix of the energy functional have BCCB structure.

Let 𝑓1(𝑧)=(𝛼/2𝜇)𝑤𝑛+1𝐷𝑧2+𝜆𝑛/𝜇,𝐷𝑧 and 𝑓2(𝑧)=(1/2)𝐾𝑘=1𝑊𝑘𝑧𝑦𝑘2. Then, 𝑓2(𝑧𝑛)=𝐵𝑇(𝑅𝐵𝑧𝑛𝑠) with 𝑠=𝑘𝑀𝑘𝑇𝑆𝑇𝑦𝑘. Define operators 𝐴=𝜕𝑓1 and 𝐵=𝑓2; the FBS method (3.3) with 𝛾 as the step size applied to (3.14) leads to the iterative scheme𝑧𝑛+1=argmin𝑧𝛼2𝑤𝑛+1𝐷𝑧2+𝜆𝑛,𝐷𝑧+𝜇2𝛾𝑧𝑧𝑛+𝛾𝑓2(𝑧𝑛)2.(3.19) The minimizer 𝑧𝑛+1 is given by the normal equations𝐷𝑇𝐷+𝜇𝛾𝛼𝐼𝑧𝑛+1=𝐷𝑇𝑤𝑛+1𝜆𝑛𝛼+𝜇𝛾𝛼𝑧𝑛𝛾𝑓2(𝑧𝑛).(3.20) Under the periodic boundary condition, (3.20) can be solved by two FFTs, which simultaneously performs the LR measurements fusion and deblurring treatment. In the next subsection, the fusion part will be shown to be a simple noniterative.

Notice that 𝐵=𝑓2 is Lipschitz continuous with Lipschitz constant 𝛽=𝑅𝐵2, where 𝐵 (resp., 𝑅) is the matrix norm of 𝐵 (resp., 𝑅). From [41,Theorem 2.3.19], (1/𝛽)𝐵 is firmly nonexpansive. Then, the convergence of 𝑧𝑛 is ensured by Theorem 3.1, while 𝛾]0,2/𝛽[.

3.2.3. Optimal Initial Guess

Firstly, borrowed the idea of [31], we take initial data 𝑧0 from maximum likelihood (ML) estimation; that is, 𝑧0 is the minimizer of the optimization problem as follows:̂𝑧=argmin𝑧𝑘𝑆𝐵𝑀𝑘𝑧𝑦𝑘2.(3.21) The gradient descent algorithm suggests the following iterative equation for the solution of (3.21):̂𝑧𝑛+1=̂𝑧𝑛Δ𝑡𝑘𝑀𝑇𝑘𝐵𝑇𝑆𝑇𝑆𝐵𝑀𝑘̂𝑧𝑛𝑦𝑘=̂𝑧𝑛Δ𝑡𝐵𝑇𝑘𝑀𝑇𝑘𝑆𝑇𝑆𝑀𝑘𝐵̂𝑧𝑛𝑦𝑘,(3.22) where Δ𝑡 denotes step size. Let us define the blurred super-resolution image by ̂𝑢𝑛=𝐵̂𝑧𝑛. Multiplying both sides of (3.22) with 𝐵, we get̂𝑢𝑛+1=̂𝑢𝑛+Δ𝑡𝐵𝐵𝑇𝑠𝑅̂𝑢𝑛,(3.23) where 𝑠=𝑘𝑀𝑇𝑘𝑆𝑇𝑦𝑘 and 𝑅=𝑘𝑀𝑇𝑘𝑆𝑇𝑆𝑀𝑘. According to [52], the steady state solution of (3.23) is given bŷ𝑢=𝑅1𝑠.(3.24) We note that 𝑀𝑇𝑘𝑆𝑇 copies the values from the LR grid to the HR grid after proper shifting and zero filling and 𝑆𝑀𝑘 copies a selected set of pixels in HR grid back on the LR grid (Figure 1 illustrates the effect of upsampling and downsampling matrices 𝑆𝑇 and 𝑆). Neither of these two operations changes the pixel value. It is easy to show that 𝑅 is a diagonal matrix. Each diagonal entry in 𝑅 corresponds to one pixel in the super-resolution image. Its value is a nonnegative integer, counting the number of measurements contributing to it. The fusion image 𝑠 is simply the addition of the measurements after proper zero-filling interpolation and motion compensation. Thus, ̂𝑢=𝑅1𝑠 is none other than the pixel-wise average of the measurement. Therefore, the noise of ̂𝑢 is reduced due to the averaging.

380807.fig.001
Figure 1: Effect of 𝑆𝑇 and 𝑆 on image (resolution enhancement factor is three).

Because ̂𝑢=𝐵̂𝑧, Wiener filter was applied to ̂𝑢. Then, the restoration image ̂𝑧 is taken as the initial data 𝑧0.

As stated previously, precondition design procedure of the initial data 𝑧0 is summed up in Algorithm 1 as follows:

alg1
Algorithm 1

3.2.4. Algorithm Description

To sum up the above arguments, the complete resulting algorithm is summarized in Algorithm 2 as follows.

alg2
Algorithm 2

3.2.5. Some Complexity Notes

It is clear that the complexity of the proposed algorithm mainly includes three parts. The calculation in (3.17) and (3.15) have linear-time complexity of order 𝑂(𝑁2) for an 𝑁×𝑁 image. Hence, the 𝑤-subproblem (3.16) and 𝜆 can be solved quickly. The solution of the 𝑧-subproblem (3.20) requires two FFTs (including one inverse FFT), which has a total complexity in the order of 𝑂(𝑁2log(𝑁2))=𝑂(𝑁2log(𝑁)).

4. Experimental Results

In this section, we present some experimental examples to demonstrate the performance of our method. We use the three 256×256 test images (“Lena” Figure 2(b), “Cameraman” Figure 4(b), and “Fingerprint” Figure 8(b)) for the synthetic test. A sequence of LR frames of 64×64 pixels from the original image is generated as follows. First, the original image was shifted by one pixel in the vertical direction. Then, to simulate the effect of camera PSF, this shifted image was convolved with a Gaussian blur kernel. The resulting image was subsampled by the factor 4 in each direction. The same approach with different motion vectors in the vertical and horizontal directions was used to produce eight LR images from the original scene. The resulting LR frames were corrupted with white Gaussian noise. All experiments were performed under Windows XP and MATLAB v7.1 running on a desktop with an Intel Core Dual Processor 3.00 GHz and 4.00 GB of memory.

fig2
Figure 2: Results of different resolution enhancement methods applied to “Lena” image degraded by Gaussian blur with support size 3×3, standard deviation 0.5, and Gaussian noise with variance 5. (a) One LR image; (b) original image; reconstructed image (c) using LDFP method, (d) using AMM method, and (e) using the proposed method (𝜇=130).

For the objective comparison between the original HR and SR reconstructed images, we measure the peak signal-to-noise ratio (PSNR) and the relative error (ReErr) defined as PSNR=10log10𝐿1𝑁1×𝐿2𝑁2max{𝑧}2𝑧𝑧2,ReErr=𝑧𝑧2𝑧2,(4.1) where 𝑧 and 𝑧 are the original and the SR reconstructed images, respectively, and 𝐿1𝑁1×𝐿2𝑁2 represents the image size.

We compare the proposed method (operator splitting method, OSM) with the lagged diffusivity fixed point iteration (LDFP) [28, 29] and alternating minimization method (AMM) [33]. In all the tests, we set the initial guess 𝑧0=𝐵1𝑅1𝑠. The choice of parameters in three methods all base on the tradeoff between reconstruction effect and computing time. In the proposed method, the value of 𝛾 and 𝛼 are fixed to be 1.5 and 4, respectively. The stopping criterion of all the methods is that the relative difference (ReDiff) between the successive iterative of the SR reconstructed image should satisfy the following inequality: ReDi=𝑧𝑛+1𝑧𝑛𝑧𝑛+1104.(4.2)

In the first test, we apply Gaussian kernel with window size 3×3, standard deviation 0.5, and different noise level (noise variance 𝜎2=5,20). One of LR frames is presented in Figures 25, 8(a), respectively. The corresponding reconstructed images of the three methods are shown in Figures 2, 4, 8(c)8(e) and Figures 3, 5(b)5(d), respectively. Our second test uses Gaussian kernel with support size 15×15, standard deviation 1.7, and noise variance 5. One of LR images is presented in Figures 6, 7(a), respectively. The corresponding reconstructed images by the three methods are shown in Figures 6, 7(b)7(d). We can see theses SR reconstructed images by different methods are very similar in real visual perception. In Table 1, we compare their reconstruction performances in PSNRs and ReErrs. On one hand, we see from the table that both PSNRs and ReErrs of the reconstructed images by the proposed method are better than those by the LDFP and AMM method. On the other hand, it is clear from Table 1 that the proposed method is more efficient (in iterations and computation times) than the other two methods.

tab1
Table 1: The PSNR, ReErr, number of iterations, and computational times of the reconstructed images using three methods. The numbers in the bracket in the “Iterations” row refer to the total number of inner conjugate gradient iterations.
fig3
Figure 3: Results of different resolution enhancement methods applied to “Lena” image degraded by Gaussian blur with support size 3×3, standard deviation 0.5, and Gaussian noise with variance 20. (a) One LR image; reconstructed image (b) using LDFP method, (c) using AMM method, and (d) using the proposed method (𝜇=130).
fig4
Figure 4: Results of different resolution enhancement methods applied to “Cameraman” image degraded by Gaussian blur with support size 3×3, standard deviation 0.5, and Gaussian noise with variance 5. (a) One LR image; (b) original image; reconstructed image (c) using LDFP method, (d) using AMM method, and (e) using the proposed method (𝜇=130).
fig5
Figure 5: Results of different resolution enhancement methods applied to “Cameraman” image degraded by Gaussian blur with support size 3×3, standard deviation 0.5, and Gaussian noise with variance 20. (a) One LR image; reconstructed image (b) using LDFP method, (c) using AMM method, and (d) using the proposed method (𝜇=130).
fig6
Figure 6: Results of different resolution enhancement methods applied to “Lena” image degraded by Gaussian blur with support size 15×15, standard deviation 1.7, and Gaussian noise with variance 5. (a) One LR image; reconstructed image (b) using LDFP method, (c) using AMM method, and (d) using the proposed method (𝜇=500).
fig7
Figure 7: Results of different resolution enhancement methods applied to “Cameraman” image degraded by Gaussian blur with support size 15×15, standard deviation 1.7, and Gaussian noise with variance 5. (a) one LR image; reconstructed image (b) using LDFP method, (c) using AMM method, and (d) using the proposed method (𝜇=500).
fig8
Figure 8: Results of different resolution enhancement methods applied to “Fingerprint” image degraded by Gaussian blur with support size 3×3, standard deviation 0.5, and Gaussian noise with variance 5. (a) One LR image, (b) original image; reconstructed image (c) using LDFP method, (d) using AMM method, and (e) using the proposed method (𝜇=130).

In Figure 9, we show that the convergence of the proposed method is faster than LDFP and AMM methods. The 𝑥-axis refers to the number of iterations. The 𝑦-axis in Figures 9(a), 9(b) refers to the PSNR and the 𝑦-axis in Figures 9(c), 9(d) refers to the relative difference between the successive iteration of the reconstructed image. These figures show that the proposed method can provide good quality of reconstructed images in an efficient manner.

fig9
Figure 9: Convergence performance of three methods in the Gaussian blur with support size 3×3, standard deviation 0.5, and Gaussian noise with variance 5 case: measured by PSNR value of the SR reconstructed (a) “Lena” image, (b) “Cameraman” image, measured by ReDiff value of the SR reconstructed (c) “Lena” image, (d) “Cameraman” image.

5. Conclusion

This paper proposes a general framework for multiple shifted and linear space-invariant blurred LR image frames which subsume several well-known SR models. The proposed model combines total variation (TV) regularization to formulate the SR image reconstruction as an optimization problem. Then, we propose an efficient algorithm that takes full advantage of the problem structures. As such, we propose to compute the minimizer of our SR model by applying DRS techniques (resp., ADMM) which separated the SR model into three subproblems that can be easily solved. Moreover, to speed up convergence, we provide an accelerated scheme based on precondition design of initial value and FBS. The proposed algorithm reduces the computational complexity. The good performance of the proposed explicit algorithm has been tested for synthetic data sets of several images degraded with Gaussian blur and contaminated with Gaussian white noise. Numerical results indicate that the algorithm recovers well edges and small features not appearing in the original degraded images. The experimental results indicate that the proposed algorithm has considerable effectiveness in terms of both objective measurements and visual evaluation.

Acknowledgments

This work was supported in part by the Natural Science Foundation of China under Grant no. 60802039, by the Doctoral Foundation of Ministry of Education of China under Grant no. 20070288050, NUST Research Funding under Grant no. 2010ZDJH07, and also sponsored by Natural Science Foundation of Jiangsu under Grant no. BK2010488 and “Qing Lan Project” of Jiangsu Province. The authors would like to express their gratitude to the anonymous referees for making helpful and constructive suggestions.

References

  1. R. Tsai and T. Huang, “Multi-frame image restoration and registration,” Advances in Computer Vision and Image Processing, vol. 1, no. 2, pp. 317–339, 1984. View at Google Scholar
  2. S. P. Kim, N. K. Bose, and H. M. Valenzuela, “Recursive reconstruction of high resolution image from noisy undersampled multiframes,” IEEE Transactions on Acoustics, Speech, and Signal Processing, vol. 38, no. 6, pp. 1013–1027, 1990. View at Publisher · View at Google Scholar · View at Scopus
  3. S. P. Kim and W. Y. Su, “Recursive high-resolution reconstruction of blurred multiframe images,” IEEE Transactions on Image Processing, vol. 2, no. 4, pp. 534–539, 1993. View at Publisher · View at Google Scholar · View at Scopus
  4. S. Rhee and M. G. Kang, “Discrete cosine transform based regularized high-resolution image reconstruction algorithm,” Optical Engineering, vol. 38, no. 8, pp. 1348–1356, 1999. View at Publisher · View at Google Scholar · View at Scopus
  5. R. H. Chan, T. F. Chan, L. Shen, and Z. Shen, “Wavelet algorithms for high-resolution image reconstruction,” SIAM Journal on Scientific Computing, vol. 24, no. 4, pp. 1408–1432, 2003. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  6. M. K. Ng, C. K. Sze, and S. P. Yung, “Wavelet algorithms for deblurring models,” International Journal of Imaging Systems and Technology, vol. 14, no. 3, pp. 113–121, 2004. View at Publisher · View at Google Scholar · View at Scopus
  7. N. Nguyen and P. Milanfar, “A wavelet-based interpolation-restoration method for superresolution (wavelet superresolution),” Circuits, Systems, and Signal Processing, vol. 19, no. 4, pp. 321–338, 2000. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at Scopus
  8. S. Borman and R. L. Stevenson, “Super-resolution from image sequences-a review,” in Proceedings of the Midwest Symposium on Circuits and Systems, pp. 9–12, Notre Dame, Ind, USA, August 1998.
  9. A. J. Patti, M. I. Sezan, and A. M. Tekalp, “Superresolution video reconstruction with arbitrary sampling lattices and nonzero aperture time,” IEEE Transactions on Image Processing, vol. 6, no. 8, pp. 1064–1076, 1997. View at Publisher · View at Google Scholar · View at Scopus
  10. H. Stark and P. Oskoui, “High-resolution image recovery from image-plane arrays, using convex projections,” Journal of the Optical Society of America A, vol. 6, no. 11, pp. 1715–1726, 1989. View at Publisher · View at Google Scholar · View at Scopus
  11. S. Peleg, D. Keren, and L. Schweitzer, “Improving image resolution using subpixel motion,” CVGIP: Graphical Models and Image Processing, vol. 54, pp. 181–186, 1992. View at Publisher · View at Google Scholar
  12. M. Irani and S. Peleg, “Improving resolution by image registration,” CVGIP: Graphical Models and Image Processing, vol. 53, no. 3, pp. 231–239, 1991. View at Publisher · View at Google Scholar · View at Scopus
  13. M. Elad and A. Feuer, “Restoration of a single superresolution image from several blurred, noisy, and undersampled measured images,” IEEE Transactions on Image Processing, vol. 6, no. 12, pp. 1646–1658, 1997. View at Publisher · View at Google Scholar · View at Scopus
  14. M. Protter, M. Elad, H. Takeda, and P. Milanfar, “Generalizing the nonlocal-means to super-resolution reconstruction,” IEEE Transactions on Image Processing, vol. 18, no. 1, pp. 36–51, 2009. View at Publisher · View at Google Scholar · View at Scopus
  15. H. Takeda, P. Milanfar, M. Protter, and M. Elad, “Super-resolution without explicit subpixel motion estimation,” IEEE Transactions on Image Processing, vol. 18, no. 9, pp. 1958–1975, 2009. View at Publisher · View at Google Scholar · View at Scopus
  16. J. Chung, E. Haber, and J. Nagy, “Numerical methods for coupled super-resolution,” Inverse Problems, vol. 22, no. 4, pp. 1261–1272, 2006. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  17. R. C. Hardie, K. J. Barnard, and E. E. Armstrong, “Joint MAP registration and high-resolution image estimation using a sequence of undersampled images,” IEEE Transactions on Image Processing, vol. 6, no. 12, pp. 1621–1633, 1997. View at Publisher · View at Google Scholar · View at Scopus
  18. N. A. Woods, N. P. Galatsanos, and A. K. Katsaggelos, “Stochastic methods for joint registration, restoration, and interpolation of multiple undersampled images,” IEEE Transactions on Image Processing, vol. 15, no. 1, pp. 201–213, 2006. View at Publisher · View at Google Scholar
  19. H. Shen, L. Zhang, B. Huang, and P. Li, “A MAP approach to joint motion estimation, segmentation, and super resolution,” IEEE Transactions on Image Processing, vol. 16, no. 2, pp. 479–490, 2007. View at Publisher · View at Google Scholar
  20. R. Sasahara, H. Hasegawa, I. Yamada, and K. Sakaniwa, “A color super-resolution with multiple nonsmooth constraints by hybrid steepest descent method,” in Proceedings of the IEEE International Conference on Image Processing (ICIP '05), vol. 1, pp. 857–860, Genova, Italy, September 2005. View at Publisher · View at Google Scholar · View at Scopus
  21. S. Farsiu, M. Elad, and P. Milanfar, “Multiframe demosaicing and super-resolution of color images,” IEEE Transactions on Image Processing, vol. 15, no. 1, pp. 141–159, 2006. View at Publisher · View at Google Scholar · View at Scopus
  22. T. Akgun, Y. Altunbasak, and R. M. Mersereau, “Super-resolution reconstruction of hyperspectral images,” IEEE Transactions on Image Processing, vol. 14, no. 11, pp. 1860–1875, 2005. View at Publisher · View at Google Scholar · View at Scopus
  23. C. A. Segall, A. K. Katsaggelos, R. Molina, and J. Mateos, “Bayesian resolution enhancement of compressed video,” IEEE Transactions on Image Processing, vol. 13, no. 7, pp. 898–910, 2004. View at Publisher · View at Google Scholar · View at Scopus
  24. C. A. Segall, R. Molina, and A. K. Katsaggelos, “High-resolution images from low-resolution compressed video,” IEEE Signal Processing Magazine, vol. 20, no. 3, pp. 37–48, 2003. View at Publisher · View at Google Scholar · View at Scopus
  25. N. Nguyen, P. Milanfar, and G. Golub, “A computationally efficient super-resolution image reconstruction algorithm,” IEEE Transactions on Image Processing, vol. 10, no. 4, pp. 573–583, 2001. View at Google Scholar
  26. A. C. Yau, N. K. Bose, and M. K. Ng, “An efficient algorithm for superresolution in medium field imaging,” Multidimensional Systems and Signal Processing, vol. 18, no. 2-3, pp. 173–188, 2007. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  27. N. K. Bose, M. K. Ng, and A. C. Yau, “A fast algorithm for image super-resolution from blurred observations,” EURASIP Journal on Applied Signal Processing, vol. 2006, 14 pages, 2006. View at Publisher · View at Google Scholar · View at Scopus
  28. M. K. Ng, H. Shen, E. Y. Lam, and L. Zhang, “A total variation regularization based super-resolution reconstruction algorithm for digital video,” EURASIP Journal on Advances in Signal Processing, vol. 2007, Article ID 74585, 16 pages, 2007. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at Scopus
  29. T. F. Chan, M. K. Ng, A. C. Yau, and A. M. Yip, “Superresolution image reconstruction using fast inpainting algorithms,” Applied and Computational Harmonic Analysis, vol. 23, no. 1, pp. 3–24, 2007. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  30. C. B. Xiao, J. Yu, and Y. Xue, “A high-efficiency super-resolution reconstruction algorithm from image/video sequences,” in Proceedings of the 3rd IEEE International Conference on Signal Image Technologies and Internet Based Systems (SITIS '07), pp. 573–580, Shanghai, China, December 2007. View at Publisher · View at Google Scholar · View at Scopus
  31. M. Elad and Y. Hel-Or, “A fast super-resolution reconstruction algorithm for pure translational motion and common space-invariant blur,” IEEE Transactions on Image Processing, vol. 10, no. 8, pp. 1187–1193, 2001. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at Scopus
  32. S. Farsiu, M. D. Robinson, M. Elad, and P. Milanfar, “Fast and robust multiframe super resolution,” IEEE Transactions on Image Processing, vol. 13, no. 10, pp. 1327–1344, 2004. View at Publisher · View at Google Scholar · View at Scopus
  33. L. L. Huang, L. Xiao, Z. H. Wei, and J. Zhang, “A fast decoupling algorithm for image super-resolution reconstruction of space-invariant system,” Acta Automatica Sinica, vol. 36, no. 2, pp. 229–236, 2010 (Chinese). View at Publisher · View at Google Scholar · View at Scopus
  34. S. Farsiu, D. Robinson, M. Elad, and P. Milanfar, “Advances and challenges in super-resolution,” International Journal of Imaging Systems and Technology, vol. 14, no. 2, pp. 47–57, 2004. View at Publisher · View at Google Scholar · View at Scopus
  35. L. I. Rudin, S. Osher, and E. Fatemi, “Nonlinear total variation based noise removal algorithms,” Physica D, vol. 60, no. 1–4, pp. 259–268, 1992. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at Scopus
  36. L. Zhang, H. Zhang, H. Shen, and P. Li, “A super-resolution reconstruction algorithm for surveillance images,” Signal Processing, vol. 90, no. 3, pp. 848–859, 2010. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at Scopus
  37. Y. Wang, W. Yin, and Y. Zhang, “A fast algorithm for image deblurring with total variation regularization,” Tech. Rep. 07-10, Department of Computational and Applied Mathematics, Rice University, Houston, Tex, USA, 2007. View at Google Scholar
  38. Y. Wang, J. Yang, W. Yin, and Y. Zhang, “A new alternating minimization algorithm for total variation image reconstruction,” SIAM Journal on Imaging Sciences, vol. 1, no. 3, pp. 248–272, 2008. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  39. Y. Huang, M. K. Ng, and Y.-W. Wen, “A fast total variation minimization method for image restoration,” Multiscale Modeling & Simulation, vol. 7, no. 2, pp. 774–795, 2008. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  40. T. Goldstein and S. Osher, “The split Bregman method for L1-regularized problems,” SIAM Journal on Imaging Sciences, vol. 2, no. 2, pp. 323–343, 2009. View at Publisher · View at Google Scholar
  41. L. Xiao, L. L. Huang, and Z. H. Wei, “A weberized total variation regularization-based image multiplicative noise removal algorithm,” EURASIP Journal on Advances in Signal Processing, vol. 2010, Article ID 490384, 2010. View at Publisher · View at Google Scholar · View at Scopus
  42. H. G. Chen, Forward backward splitting techniques: theory and applications, Ph.D. thesis, University of Washington, 1994.
  43. S. Setzer, Splitting methods in image processing, Ph.D. thesis, University of Mannheim, 2009.
  44. J. Yang, Y. Zhang, and W. Yin, “A fast alternating direction method for TVL1-L2 signal reconstruction from partial Fourier data,” IEEE Journal on Selected Topics in Signal Processing, vol. 4, no. 2, pp. 288–297, 2010. View at Publisher · View at Google Scholar
  45. P. L. Combettes and J. C. Pesquet, “A Douglas-Rachford splitting approach to nonsmooth convex variational signal recovery,” IEEE Journal on Selected Topics in Signal Processing, vol. 1, no. 4, pp. 564–574, 2007. View at Publisher · View at Google Scholar · View at Scopus
  46. G. Steidl and T. Teuber, “Removing multiplicative noise by douglas-rachford splitting methods,” Journal of Mathematical Imaging and Vision, vol. 36, no. 2, pp. 168–184, 2010. View at Publisher · View at Google Scholar · View at Scopus
  47. J. Eckstein and D. P. Bertsekas, “On the Douglas-Rachford splitting method and the proximal point algorithm for maximal monotone operators,” Mathematical Programming, vol. 55, no. 3, pp. 293–318, 1992. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  48. E. Esser, “Applications of Lagrangian-based alternating direction methods and connections to split Bregman,” Tech. Rep. 09-31, UCLA CAM, 2009. View at Google Scholar
  49. S. Setzer, “Split Bregman algorithm, Douglas-Rachford splitting, and frame shrinkage,” in Proceedings of the 2nd International Conference on Scale Space Methods and Variational Methods in Computer Vision, vol. 5567 of Lecture Notes in Computer Science, pp. 464–476, Springer, 2009.
  50. S. Setzer, G. Steidl, and T. Teuber, “Deblurring Poissonian images by split Bregman techniques,” Journal of Visual Communication and Image Representation, vol. 21, no. 3, pp. 193–199, 2010. View at Publisher · View at Google Scholar · View at Scopus
  51. M. V. Afonso, J. M. Bioucas-Dias, and M. A. T. Figueiredo, “Fast image recovery using variable splitting and constrained optimization,” IEEE Transactions on Image Processing, vol. 19, no. 9, pp. 2345–2356, 2010. View at Publisher · View at Google Scholar
  52. D. Bertsekas, Nonlinear Programming, Athena Scientific, Belmont, Mass, USA, 1995.