## Support-assisted optical superresolution of low-resolution image sequences: the one-dimensional problem

Optics Express, Vol. 17, Issue 25, pp. 23213-23233 (2009)

http://dx.doi.org/10.1364/OE.17.023213

Acrobat PDF (296 KB)

### Abstract

We analyze the problem of optical superresolution (OSR) of a one-dimensional (1D) incoherent spatial signal from undersampled data when the support of the signal is known in advance. The present paper corrects and extends our previous work on the calculation of Fisher information (FI) and the associated Cramer-Rao lower bound (CRB) on the minimum error for estimating the signal intensity distribution and its Fourier components at spatial frequencies lying beyond the optical band edge. The faint-signal and bright-signal limits emerge from a unified noise analysis in which we include both additive noise of detection and shot noise of photon counting via an approximate Gaussian statistical distribution. For a large space-bandwidth product, we derive analytical approximations to the exact expressions for FI and CRB in the faint-signal limit and use them to argue why achieving any significant amount of unbiased bandwidth extension in the presence of noise is a uniquely challenging proposition. Unlike previous theoretical work on the subject of support-assisted bandwidth extension, our approach is not restricted to specific forms of the system transfer functions, and provides a unified analysis of both digital and optical superresolution of undersampled data.

© 2009 Optical Society of America

## 1. Introduction

1. S. Prasad, “Digital superresolution and the generalized sampling theorem,” J. Opt. Soc. Am. A **24**, 311–325 (
2007) [CrossRef]

*a priori*knowledge in addition to the image data in order to perform true bandwidth extension, or optical superresolution (OSR), by means of data inversion. In fact, most algorithms that claim to achieve OSR incorporate prior knowledge of one form or another, often implicitly, into their data inversion protocol. Enforcing a finite-support constraint in the physical domain is the essential basis of certain OSR algorithms [3

3. R. Gerchberg, “Superresolution through error energy reduction,” Opt. Acta **21**, 709–721 (
1974). [CrossRef]

4. A. Papoulis, “A new algorithm in spectral analysis and band-limited extrapolation,” IEEE Trans. Circuits Syst. **CAS-22**, 735–742 (
1975). [CrossRef]

5. S. Plevritis and A. Macovski, “Spectral extrapolation of spatially bounded images,” IEEE Trans. Medical Imaging **14**, 487–497 (
1995). [CrossRef]

6. W. Richardson, “Bayesian-based iterative method of image restoration,” J. Opt. Soc. Am. **62**, 55–59 (
1972). [CrossRef]

7. L. Lucy, “An iterative technique for the rectification of observed distributions,” Astron. J. **79**, 745–754 (
1974) [CrossRef]

8. L. Shepp and Y. Vardi, “Maximum likelihood reconstruction in positron emission tomography,” IEEE Trans. Medical Imaging **1**, 113–122 (
1982). [CrossRef]

9. T. Holmes, “Maximum-likelihood image restoration adapted for noncoherent optical imaging,” J. Opt. Soc. Am. A **5**, 666–673 (
1988). [CrossRef]

10. T. Holmes and Y.-H. Liu, “Richardson-Lucy/maximum likelihood image restoration algorithm for fluorescence microscopy: further testing,” Appl. Opt. **28**, 4930–4938 (
1989). [CrossRef] [PubMed]

11. J.-A. Conchello, “Superresolution and convergence properties of the expectation-maximization algorithm for maximum-likelihood deconvolution of incoherent images,” J. Opt. Soc. Am. A **15**, 2609–2619 (
1998). [CrossRef]

12. E. Boukouvala and A. Lettington, “Restoration of astronomical images by an iterative superresolving algorithm,” Astron. Astrophys. **399**, 807–811 (
2003). [CrossRef]

14. R. Narayan and R. Nityananda,“Maximum entropy image restoration in astronomy,” Ann. Rev. Astron. Astrophys. **24**, 127–170 (
1986). [CrossRef]

16. D. Fried, “Analysis of the CLEAN algorithm and implications for superresolution,” J. Opt. Soc. Am. A **12**, 853–860 (
1995). [CrossRef]

17. P. Magain, F. Courbin, and S. Sohy, “Deconvolution with correct sampling,” Astrophys. J. **494**, 472–477 (
1998). [CrossRef]

18. F. Pijpers, “Unbiased image reconstruction as an inverse problem,” Mon. Not. Roy. Astron. Soc. **307**, 659–668 (
1999). [CrossRef]

20. J. Starck, E. Pantin, and F. Murtagh, “Deconvolution in astronomy: A review,” Publ. Astron. Soc. Pacific **114**, 1051–1069 (
2002). [CrossRef]

*object*, is known in advance. This work was an attempt to generalize a large body of previous theoretical literature [21

21. C. Rushforth and R. Harris, “Restoration, resolution, and noise,” J. Opt. Soc. Am. **58**, 539–545 (
1968). [CrossRef]

24. M. Bertero and E. R. Pike, “Resolution in diffraction-limited imaging, a singular value analysis: I. The case of coherent illumination,” Opt. Acta **29**, 727–746 (
1982). [CrossRef]

25. M. Bertero and C. De Mol, “Superresolution by data inversion,” Progress in Optics **XXXVI**, 129–178 (
1996). [CrossRef]

26. C. Matson and D. Tyler, “Primary and secondary superresolution by data inversion,” Opt. Express **14**, 456–473 (
2006). [CrossRef] [PubMed]

5. S. Plevritis and A. Macovski, “Spectral extrapolation of spatially bounded images,” IEEE Trans. Medical Imaging **14**, 487–497 (
1995). [CrossRef]

26. C. Matson and D. Tyler, “Primary and secondary superresolution by data inversion,” Opt. Express **14**, 456–473 (
2006). [CrossRef] [PubMed]

## 2. Imaging model

*h*

_{0}(

*x*), a continuously defined image,

*g*(

*x*), may be expressed in terms of the object intensity function,

*f*(

*x*), as the convolution

*p*(

*x*) is the pixel sensitivity function,

*x*=

_{k}*kδw*is the position coordinate of the center of the

*k*th pixel in terms of the pixel pitch,

*δw*, and

*n*is the additive read-out noise at the

_{k}*k*th pixel. A complete noise model for the sensor must also include the Poisson noise of counting for which the first term on the right-hand side (RHS) must be regarded as the mean count rate at the

*k*th pixel.

1. S. Prasad, “Digital superresolution and the generalized sampling theorem,” J. Opt. Soc. Am. A **24**, 311–325 (
2007) [CrossRef]

27. D. Robinson and P. Milanfar, “Statistical performance analysis of superresolution,” IEEE Trans. Image Process. **15**, 1413–1428 (
2006). [CrossRef] [PubMed]

*H*

_{0}(

*u*), and the detector-pixel transfer function (PTF),

*P*(

*u*), we may express the intensity at the

*k*th pixel of the

*m*th image as [2]

*t*is the sub-pixel shift of that image and

_{m}*n*

^{(m)}

*is the additive noise at its*

_{k}*k*th pixel.

*P*(

*u*)=

*Cδ w*sinc(

*uδw*), where the sinc function is defined as sinc(

*x*)≡sin

*πx*/(

*πx*) and

*C*is the quantum efficiency of the sensor pixels. We shall assume

*C*=1 in the rest of the paper. For a clear aperture of optical bandwidth

*B*

_{0}, the incoherent OTF vanishes outside the spatial-frequency range |

*u*|<

*B*

_{0}where its functional form is (1-|

*u*|/

*B*

_{0}). Expression (3) may thus be rewritten as the following integral over the normalized frequency

*u*̄=

*u*/

*B*

_{0}:

*δx*=1/(2

_{c}*B*

_{0}), specifically

*x*̄

*=*

_{k}*x*/

_{k}*δx*and

_{c}*t*̄

*=*

_{m}*t*/

_{m}*δx*. The symbol

_{c}*χ*denotes the ratio,

*δw*/

*δxc*, of detector-sampling and critical optical-sampling intervals. For

*χ*<1, the detector oversamples the image while for

*χ*>1 it undersamples the image.

## 3. Finite object support and the object spectrum

*f*(

*x*) is supported compactly on the interval [-

*L*,

*L*], it has a Fourier-series expansion,

*F*(

*u*), in terms of its spatial-frequency samples,

*F*

*=*

_{ℓ}*F*(

*u*=

*ℓ/2L*), as

*B*

_{0},

*B*

_{0}], is aliased into the measurement band and thus in principle contained in the acquired data. The finiteness of the support is directly responsible for this aliasing. Note, however, that because the width of the sinc function in Eq. (7) is of order 1/(2

*L*), the maximum bandwidth extension that can be expected is also of that order. This is consistent with previous results on unbiased, primary OSR [28

28. P. Sementilli, B. Hunt, and M. Nader, “Analysis of of the limit of superresolution in incoherent imaging,” J. Opt. Soc. Am. A **10**, 2265–2276 (
1993). [CrossRef]

25. M. Bertero and C. De Mol, “Superresolution by data inversion,” Progress in Optics **XXXVI**, 129–178 (
1996). [CrossRef]

*Q*=2

*LB*

_{0}is the space-bandwidth product (SBP) and

*u*̄ has been simply relabeled as

*u*.

*ℓ*, of the index,

_{max}*i.e*., to the range -

*ℓ*≤

_{max}*ℓ*≤

*ℓ*. The ratio

_{max}*S*=

*ℓ*/

_{max}*Q*, when it exceeds 1, is the sought OSR factor. In the absence of noise, the (2

*SQ*+1) object spectral components contained in the truncated sum can be solved for by acquiring at least (2

*SQ*+1)/N independent LR images, each with

*N*pixels.

## 4. Fisher information and estimation error for fourier samples

*L*,

*L*]. A simply calculable lower bound on this error is the CRB, determined by the diagonal elements of the inverse of the FI matrix [29].

*Fℓ*, sometimes called its pseudo-variance, is given by the corresponding diagonal element of the inverse of a complex Hermitian FI matrix. The

*ℓℓ*′ element of the FI matrix is defined as a statistical self-average in two equivalent ways,

*P*=

*P*({

*g*

^{(m)}

*}|{*

_{k}*F*,

_{ℓ}*F**

*}) denotes the PDF of the image data set, {*

_{ℓ}*g*

^{(m)}

*k*}, corresponding to the given object. The reality of the object implies

*F**

*=*

_{ℓ}*F*-

*ℓ*. It is therefore sufficient to consider only non-negative

*ℓ*,

*ℓ*′ values when considering the matrix elements of

**J**.

*biased*estimation of SR frequencies, so our analysis based on using simply the inverse of the FI matrix to provide a lower bound on the variance of any

*unbiased*estimation must be regarded as merely suggestive of the actual sensitivity limits on the performance of most practical OSR algorithms. Limits on a biased estimation of the SR frequencies may also be calculated in the FI formalism, but such limits depend critically on the gradient of the estimator bias relative to the parameters being estimated by the specific algorithm and, unlike unbiased-estimation CRBs, are thus algorithm-dependent. It is worth noting that the absence or presence of estimator bias in OSR has been utilized recently to sharpen the discussion of OSR in terms of primary and secondary OSR [26

26. C. Matson and D. Tyler, “Primary and secondary superresolution by data inversion,” Opt. Express **14**, 456–473 (
2006). [CrossRef] [PubMed]

*k*th pixel of the

*m*th image given by the first term on the RHS of Eq. (8) and variance

*σ*

^{2}

*. In this case, as we argue in Appendix A, the average (10) is evaluated easily in terms of the coefficients,*

_{D}*c*

^{(m)}

*k*ℓ,

*c*

^{(m)}-

*is the column vector (*

_{k}*c*

^{(m)}

_{k0},

*c*

^{(m)}

_{k1},…)

*and*

^{T}*c*

^{(m)†}

_{-k}its Hermitian transpose.

*σ*≫1 or the mean count rate 〈

_{D}*g*

^{(m)}

*〉≫1, we may do this approximately via a Gaussian PDF with mean 〈*

_{k}*g*

^{(m)}

*〉 and variance (*

_{k}*σ*

^{2}

*+〈*

_{D}*g*

^{(m)}

*〉). Under such a Gaussian approximation to the exact PDF for the mixed-noise model, which we assume to be valid throughout this paper, the FI (12) is modified by bringing*

_{k}*σ*

^{2}

*inside the*

_{D}*m*,

*k*double sum and replacing it by the sum noise variance,

*σ*

^{2}

*+〈*

_{D}*g*

^{(m)}

*〉, pixel-by-pixel,*

_{k}## 5. The large-Q limit

*ℓ*,

*ℓ*′) plane of matrix-element labels, namely 0≤

*ℓ*,

*ℓ*′<

*Q*;

*Q*≤

*ℓ*,

*ℓ*′≤

*ℓ*; 0≤

_{max}*ℓ*<

*Q*≤

*ℓ*′≤

*ℓ*; and 0≤

_{max}*ℓ*′<

*Q*≤′≤

*ℓ*. The off-diagonal submatrices, which are Hermitian adjoints of each other, represent the cross-sensitivity of the estimation of the SR frequencies on the DR frequencies and vice-versa.

_{max}*Q*≫1, the aliasing sinc function inside the integral (9) is sharply peaked at

*u*=

*ℓ/Q*with a width of order 1/

*Q*. Hence for the spatial frequencies well inside the passband for which |

*Q*-

*ℓ*|≫1, we may approximate the integral by replacing the relatively slowly varying factor in the integrand, namely (1-|

*u*|)sinc(

*χu*/2), by its value at

*u*=

*ℓ/Q*and extending the limits of integration about

*v*≡(

*u*-

*ℓ*/

*Q*)=0 to ±∞. Thus, we have

*δ*function, we may express the integral in Eq. (14) as

*t*̃

*=*

_{m}*t*̄

*/*

_{m}*χ*is the sub-pixel shift of the mth image as a fraction of an image-pixel width. The integral (16) evaluates as (1/

*Q*)Θ(1-

*χ*|

*k*-

*t*̃

*|/*

_{m}*Q*), where Θ denotes the Heaviside unit step function that takes the value 1 when its argument is positive and 0 otherwise. Putting all these simplifications together, we have the following approximate expression for the coefficient

*c*

^{(m)}

*valid for large*

_{k}*Q*:

*k*-

*t*̄

*|<*

_{m}*Q*/

*χ*. This constrains the sum for each image m to be over only those pixels

*k*over which the image signal is supported.

### 5.1. The large-additive-noise limit

*i.e*. when

*σ*

^{2}

*≫〈*

_{D}*g*

^{(m)}

*k*〉 for all

*m*,

*k*, so the denominator of the summand can be replaced simply by

*σ*

^{2}

*. The double sum then reduces to a simple geometric-series inner sum over*

_{D}*k*running from -

*Q*/

*χ*to

*Q/χ*, which equals (2

*Q*/

*χ*)

*ℓ′, followed by a trivial outer sum over*

_{δℓ}*m*, which equals

*M*, the total number of LR images. This leads to the following approximately diagonal form for the diagonal FI submatrix corresponding to frequency samples inside the passband:

*Q*, the elements of the off-diagonal submatrices are essentially zero, at least well in their interior. As a result, regardless of the form of the remaining, diagonal FI submatrix that involves only the SR frequencies, the inverse of the overall FI matrix, when similarly sub-divided into four submatrices, consists of diagonal submatrices that are simply inverses of the corresponding diagonal submatrices of the overall FI matrix. In other words, because of the approximate vanishing of the cross-sensitivity-submatrix elements, the CRBs on the errors of an unbiased estimation of the DR frequencies well within the passband or of the SR frequencies lying well outside are independent of each other for large SBP.

*u*=

_{ℓ}*ℓ/Q*, where the detector PTF (

*χ*/2)sinc(

*χuℓ*/2) vanishes, which can happen for

*χ*>2. This divergence may be traced to expression (14) which fails to provide an accurate approximation of the exact expression (9) near such frequencies. The finite width of the support-based aliasing sinc function cannot be neglected in Eq. (9) even in the large-

*Q*limit, so that the exact CRB is in fact finite for such frequencies. The prior knowledge of the object support injects information at the zero crossings of the PTF within the optical passband, thus enabling DSR throughout that band. Without prior information, such object frequencies would be filtered out perfectly by the sensor PTF and could not be estimated.

*Q*

^{2}unlike the 1/

*Q*scaling of expression (19). Thus, since

*Q*≫1, it has much smaller elements, and its inverse, which determines the CRBs, has much larger elements than those given by Eq. (20). The CRBs for these SR frequencies are, in fact, considerably larger than the simple scaling argument suggests, since when

*ℓ*,

*ℓ*′≫

*Q*, the

*ℓ*,

*ℓ*′ dependence of expression (21) may be approximated as being of a factorized form, (-1)

*ℓ*-

*ℓ*′/(

*ℓ*

*ℓ*′). Such matrices are highly ill-conditioned, since their rank is small (if the factorization were exact for all

*ℓ*,

*ℓ*′, then the rank would be 1, independent of the dimension of the submatrix). Physically, this is due to the fact that the image data contained within the diffraction passband become increasingly insensitive to – and thus increasingly lacking in information about – increasingly larger SR frequencies. We expect the rise of the CRB for these SR frequencies to become dramatically steep as

*ℓ*is increased well beyond

_{max}*Q*.

## 6. Numerical results and discussion of OSR in the spatial-frequency plane

*w*is a fractional width parameter and

*K*may be interpreted as the integrated source flux whenever

*w*is small compared to 1. For this object, we now illustrate several essential features of the variation of the CRB for a joint, unbiased estimation of its spatial frequencies for four different values of

*Q*, namely 5, 10, 15, and 20. On the plots the spatial-frequency samples are labeled by their index

*ℓ*, which takes integer values between 0 and

*ℓ*. In our numerical calculations of the exact expression (13) for the FI and its inverse, the set of jointly estimated samples was enlarged incrementally by adding the next higher-frequency sample pair at a time,

_{max}*i.e*., by incrementing

*ℓ*by 1. Three values of

_{max}*ℓ*, namely

_{max}*Q*,

*Q*+1, and

*Q*+2, were considered corresponding to the inclusion of zero, one, and two pairs of samples outside the optical bandwidth, respectively. Incrementing

*ℓ*by one unit corresponds to a fractional bandwidth extension of amount 1/

_{max}*Q*, which is a 20% extension for

*Q*=5 but only a 5% extension for

*Q*=20.

*w*=1) from a single (

*M*=1) critically sampled (

*χ*=1) LR image for the four values of

*Q*. The peak SNR, defined as the ratio of the maximum object pixel intensity to the standard deviation,

*σ*, of the additive noise, was taken to be very small at 0.001, so the additive noise was dominant and the CRBs, according to Eq. (13), essentially independent of the object brightness distribution. Here the CRBs were normalized by dividing out a factor of

_{D}*σ*

^{2}

*before plotting them. A number of observations may be made. First, the CRBs for estimating the DR frequencies remain essentially unchanged for moderate to large values of*

_{D}*Q*(Figs. 1(c) and (d)) even as more and more of the SR frequencies are jointly estimated. This is in consistency with our earlier result that at least for large

*Q*the cross-sensitivity of the DR and SR frequencies in their joint estimation is expected to be small. Furthermore, the approximate result (20), indicated by the dashed line segments on each plot, is quite accurate even for

*Q*=15, although for smaller

*Q*values the greater cross-sensitivity of DR and SR frequencies causes DR-frequency estimation to become more noisy as more SR frequencies are jointly estimated. Second, the normalized CRB, which may be regarded as the inverse of the maximum SNR with which the object spatial power spectrum may be estimated, rises dramatically, in accordance with the approximate analysis of Sec. 5, as the optical band edge is approached from below and then exceeded. Finally, with increasing

*Q*the CRB for estimating the DR frequencies increases linearly with

*Q*, when expressed as a function of the scaled sample index,

*ℓ/Q*, but the SR-frequency-estimation variance has dramatically larger CRBs.

^{3}. In this bright-source limit, the estimation problem is dominated by the shot (Poisson) noise of counting, and the CRBs were normalized by dividing them by the square of the mean source intensity (which is proportional to the strength of the dc component of the source). The source width parameter was sufficiently large so that the source intensity had only a weak variation over its support, and the general behavior of the CRB in the present bright-source limit is, as expected, much the same as that we saw for the uniform-additive-noise considered in the previous figures. However, the normalized CRB value at each spatial frequency is about 3–4 orders smaller because of the large source brightness.

^{3}). In all of our figures,

*L*=1 and

*Q*=20, so the optical system bandwidth,

*B*

_{0}=

*Q*/(2

*L*), equals 10. The spatial spectrum of the Gaussian source (22) has a characteristic width of order 2/(

*πwQ*), when expressed as a fraction of the optical system bandwidth. The support-based sinc sampling function, on the other hand, has a chacteristic frequency width of order 2/

*Q*. For

*w*=1, the source spectral width is then less than the sampling-function bandwidth, which itself is a small fraction of the optical system bandwidth for

*Q*≫1. Thus, there is little signal at the SR frequencies that can get coupled into the optical bandwidth, leading to essentially no support-induced gain in information at the DR or SR frequencies for such wide Gaussian sources. By contrast, for an extremely narrow source,

*w*=0.05, the source spectral profile has a significant value at the optical band edge. In spite of the narrow bandwidth of the spectral sampling function (for large

*Q*), the relatively slow variation of the object spectrum near the band edge implies a relatively strong aliasing of the SR frequencies of the source into the optical band. However, this makes it difficult to disambiguate and thus estimate the SR frequencies individually. The optimal condition holds when the source spectral profile is neither too small nor too slowly varying at the band edge so the support-based aliasing yields good sensitivity of the data just below the band edge to the individual SR frequencies lying just above. This condition is met for

*w*=0.1 for which we find that the CRBs for an unbiased estimation of the frequency samples on the whole are smaller by a factor of 2–3 compared to those for other source widths. The two limits of very narrow and very wide sources lead to little source-based information gain when

*Q*is not small.

*M*, of sub-pixel-shifted LR images used to estimate the SR frequencies. An important motivation for our present work has been to advance an integrated treatment of DSR and OSR when images are in general under-sampled by the detector pixels. For the case of critically sampled images, the use of more than one LR image can affect only the relative variance of estimation noise by reducing it by a factor of

*M*corresponding to the standard noise reduction due to multiple independent measurements of a specific quantity. This is seen in Figs. 5 for the same additive-noise-dominated estimation CRBs as covered in Fig. 1, except that

*M*=5. We verified the same factor-

*M*reduction factor for critically sampled images in the bright-source limit too.

*χ*>1), the availability of multiple sub-pixel-shifted LR images should provide for nontrivial DSR as well as noise reduction. We show this in the next set of figures where

*χ*equals 1 (Fig. 6(a)), 2 (Fig. 6(b)), and 5 (Fig. 6(c)), while keeping

*M*=5, in the bright-source limit. For

*M*larger than

*χ*, not only is DSR possible but the normalized CRBs are also reduced by the noise-averaging factor

*M/χ*. This is readily seen for

*χ*=1 by comparing Fig. 6(a) with Fig. 2(d), particularly at the lower spatial frequencies. By contrast, as the last figure (Fig. 6(d)) in this set shows, having fewer images than the under-sampling factor,

*i.e*.,

*M*<

*χ*, in a joint estimation of the DR and SR frequencies dramatically increases the noise of estimating even those DR frequencies that would be normally resolved by the detector. The same, although not shown here, is true in the faint-source limit as well. These considerations underscore the fact that the number of independent spatial degrees of freedom for a support-limited, band-limited imaging environment is given by the total number of critical samples, spaced consecutively by distance 1/(2

*B*) over the spatial support, 2

*L*, of the object, namely by 2

*Q*, as has been widely recognized before [21

21. C. Rushforth and R. Harris, “Restoration, resolution, and noise,” J. Opt. Soc. Am. **58**, 539–545 (
1968). [CrossRef]

25. M. Bertero and C. De Mol, “Superresolution by data inversion,” Progress in Optics **XXXVI**, 129–178 (
1996). [CrossRef]

*ℓ*=8 and 16 signal reduced information at spatial frequencies that are multiples of 2/

*χ*(times

*B*) where the detector PTF vanishes identically. Without any support-based convolutions in the spatial-frequency domain, these frequencies would have infinite CRBs (i.e., no data sensitivity) and these cusps would be replaced by true divergences, as we noted earlier in Sec. 5.1. This is an instance of information at certain spatial frequencies, even those within the optical passband, that without any prior knowledge would be simply inaccessible to measurements [1

1. S. Prasad, “Digital superresolution and the generalized sampling theorem,” J. Opt. Soc. Am. A **24**, 311–325 (
2007) [CrossRef]

*χ*>2, the CRB plots exhibit cusp-like behavior at those spatial frequencies whose index,

*ℓ*, is an integral multiple of

*ℓ*

_{0}≡2

*Q*/

*χ*. For

*χ*=2, as in Fig. 6(b), a PTF zero occurs precisely at the optical band edge. The vanishing of both the PTF and OTF there yields dramatically higher CRBs for frequency samples near the band edge, as a comparison with Fig. 6(a) easily shows.

*Q*the minimum noise variance of adding a single pair of SR spatial frequencies to the reconstructed image,

*i.e*., increasing

*ℓ*by 1, is essentially a constant factor, of order 10

_{max}^{1.5}to 10

^{2}, larger than before. This represents a logarithmic scaling of the number of estimated SR frequencies with image SNR, whether in the additive or Poisson-noise-dominated regime, a conclusion that has also been reached by previous researchers [28

28. P. Sementilli, B. Hunt, and M. Nader, “Analysis of of the limit of superresolution in incoherent imaging,” J. Opt. Soc. Am. A **10**, 2265–2276 (
1993). [CrossRef]

**XXXVI**, 129–178 (
1996). [CrossRef]

**14**, 456–473 (
2006). [CrossRef] [PubMed]

## 7. OSR in the image domain

*ℓ*/(2

_{max}*L*) is the largest spatial frequency considered in the object. Note that

*ℓ*=

_{max}*Q*refers to frequency at the optical band edge of the imager. It is assumed that the object pixel indices

*k*,

*k*′ run between -

*ℓ*and

_{max}*ℓ*-1 as well, corresponding to successive pixels that are

_{max}*L*/

*ℓ*apart. Note that this is consistent with the Shannon-Nyquist sampling theorem, since for a signal of bandwidth

_{max}*ℓ*/(2

_{max}*L*) there are in all 2

*ℓ*independent, critical, physical-domain samples successively separated by distance

_{max}*L/ℓ*over the [-

_{max}*L,L*) support of the signal. Thus, as increasingly higher degree of bandwidth extension, or OSR, is sought and achieved, it is possible to estimate physical-domain pixels that are increasingly finer.

*f*

_{0}is taken to be a sufficiently large positive number so the object intensity is non-negative everywhere over the support, [-

*L*,

*L*]. The spatial frequencies,

*ν*

_{1},

*ν*

_{2}, and

*ν*

_{3}, were so chosen that at least one of them was just above the optical band edge,

*B*

_{0}=

*Q*/(2

*L*). The object and its critically sampled mean image are shown in Fig. 7. Note the large loss of spatial resolution in the image.

^{3}corresponding to the bright-source limit, while

*Q*was fixed at 20. For

*L*=1,

*i.e*., for a support size of 2 units, the largest of the three sinusoidal frequencies, 8, 9, and 11, we chose to illustrate our results lies just above the band edge at 10. A variety of values for the detector undersampling factor

*χ*and the number

*M*of LR images were considered to give us a fuller picture of how the theoretical prediction of estimation fidelity fares against the physically expected results. From the figures, we arrive at the following important conclusions. First, the CRB is 2–4 orders of magnitude lower near the support boundaries than deep inside the support interior, at least when some OSR (

*ℓ*>

_{max}*Q*) is sought. This is expected since the knowledge of support provides the most information about the object near its boundaries, rather than in its interior, as the object is required, with no error, to have zero pixel intensity everywhere outside the support. This space-variant character has also been noted by other investigators in both one and two dimensional studies [5

5. S. Plevritis and A. Macovski, “Spectral extrapolation of spatially bounded images,” IEEE Trans. Medical Imaging **14**, 487–497 (
1995). [CrossRef]

**14**, 456–473 (
2006). [CrossRef] [PubMed]

*χ*=2,

*M*=2, to critically sampling the object, but that misses its spatial frequency 11 that lies just beyond the optical band width. On the other hand, for

*M*=3, the subpixel shift is a third of a pixel from one image to the next, which carries significant information about the SR frequency 11, unlike the coarser sub-pixel shift of 1/2. This is seen in dramatically reduced CRBs when

*M*is changed from 2 to 3. The improvement continues, but a far slower pace as

*M*is increased further, in a manner consistent with the 1/

*M*reduction of the squared noise-to-signal ratio with multiple, independently performed measurements (Fig. 8 (d)), as noted earlier.

*M*is sufficiently large. This quantifies the cost overhead of DSR, which requires a sufficiently large LR image sequence before the fidelity of OSR can begin to approach that for a single critically sampled image.

*ℓ*<

_{max}*Q*, as true for the bottom two curves on each sub-plot, the CRBs are essentially uniform over the entire image except close to the support boundaries. Mathematically, this is so because

**J**

^{-1}is a diagonal matrix in this case, as seen from Eq. (20), and thus the CRB of the

*k*image pixel,

*C*

^{(ph)}

*, as given by Eq. (25), becomes independent of the pixel index*

_{k}*k*.

## 8. Conclusions

**14**, 456–473 (
2006). [CrossRef] [PubMed]

*a priori*.

*F*of the present paper, becomes more involved. These generalizations will be discussed in future papers.

_{ℓ}### A. Gaussian Approximation to the Mixed Gaussian-Poisson Noise PDF and Approximate FI

*X*is distributed according to a Poisson distribution with mean 〈

*X*〉 and

*N*according to a zero-mean Gaussian PDF with variance

*σ*

^{2}. Its characteristic function,

*Q*(λ), defined as the expected value of exp(

_{Y}*iλY*), is the product of the characteristic functions of

*X*and

*N*evaluated at λ, namely exp{〈

*X*〉[exp(

*iλ*)-1]} and exp(-

*λ*

^{2}

*σ*

^{2}/2), respectively. Expanding the exponent of the product of these two functions in powers of λ yields the following expression for

*QY*(

*λ*):

*λ*|≤

*O*(1)/(

*σ*

^{2}+〈

*X*〉)

^{1/2}. For either

*σ*or 〈

*X*〉 large compared to 1, this implies that only values of |

*λ*| small compared to 1 are important. For such small |

*λ*|, the exponent of the second exponential factor in Eq. (28) is of order 〈

*X*〉

*λ*

^{3}≤

*O*(1)〈

*X*〉/(

*σ*

^{2}+〈

*X*〉)

^{3/2}, which is small compared to 1 provided either 〈

*X*〉≫1 or

*σ*≫1. In either case, the second exponential factor may be set approximately equal to 1 over the range of values of

*λ*over which the first exponential factor is significant. For a given 〈

*X*〉, the PDF of

*Y*, which is the inverse Fourier transform of

*QY*(

*λ*), may thus be approximated as

*σ*≫1 or 〈

*X*〉≫1.

*P*(

_{Y}*y*|〈

*X*〉) with respect to 〈

*X*〉 generates a quantity whose expected value is the FI relative to 〈

*X*〉. With approximation (29), which is valid if

*σ*

^{2}+〈

*X*〉≫1, this set of steps yields the following FI:

*y*-〈

*X*〉)〉=0.

*X*〉 is a differentiable function of a vector of parameters,

*λ*̱=(

*λ*

_{1},…,

*λ*), then the elements of the FI matrix with respect to those parameters are given by a generalization of Eq. (30) that follows immediately from the chain rule of differentiation,

_{n}*M*independent measurement channels for which 〈

*X*〉 is a vector with each component depending functionally on the parameter vector

*λ*̱. In this case we need simply append a label to the mean value, say 〈

*X*〉, to indicate a specific measurement channel and then sum expression (31) over all channels,

_{m}*m*=1,…,

*M*. For the case of a linear functional dependence of 〈

*X*〉 on

_{m}*λ*̱, the desired expression (13) is obtained in this way. By setting the mean values equal to 0 in that expression, we recover the case of additive noise alone, namely Eq. (11).

### B. FI and CRB in the Large-Q, Additive-Noise-Dominated Limit

*σ*

^{2}

*≫max(*

_{D}*g*

^{(m)}

*), the FI matrix elements are given by Eq. (11). Substituting expression (9) into expression (11), first for*

_{k}*ℓ*and then in its complex-conjugated form for

*ℓ*′ using a different integration variable, say

*u*′, yields a double frequency integral for

*Jℓℓ*′,

*S*(

*u*-

*u*′),

*t*̃

*=*

_{m}*m/M*. By separating out the summand into a product of an

*m*-dependent factor and a

*k*-dependent factor, we may separate the double sum into a product of two single geometric-series sums, which may be evaluated simply. The following expression for the double sum (33) thus results:

*ix*)=2

*i*exp(-

*ix*/2) sin(

*x*/2) to obtain the second equality from the first.

*Q*→∞, the ratio of the two sines in the final expression in Eq. (34) is highly peaked at

*u*=

*u*′ within a width of order 1/

*Q*, so for the purposes of evaluating the frequency integrals, it may be replaced in terms of the Dirac delta function by (2

*M/χ*)

*δ*(

*u*-

*u*′). Because of the vanishing of the denominator of this ratio at |

*u*-

*u*′|=2

*Mn*/

*χ, n*=1,2,…, however, this replacement is not quite accurate, but whenever

*M*≥

*χ*, this is of no consequence since |

*u*-

*u*′| cannot exceed 2 inside the bandlimited double-frequency integral (32). Physically speaking, assuming

*M*≥

*χ*is sensible since we must acquire at least as many pixels of information as there are parameters to estimate. For a down-sampling factor

*χ*, each LR image has only of order 2

*Q/χ*meaningful pixels over which the image is supported, so we must acquire at least

*M*=

*χ*images to estimate (2

*Q*+1) or more (for finite OSR) spatial-frequency samples of the object brightness distribution.

*δ*function, we may perform one of the two integrals in Eq. (32), say the

*u*′ integral, trivially so the following expression for the FI matrix elements is obtained in the limit of large

*Q*:

*ℓ*and

*ℓ*ℓ′.

*B.1.*|

*ℓ*|,|

*ℓ*′|<

*Q*

*u*=

*ℓ/Q*and

*u*=

*ℓ*′/

*Q*within a width of order 1/

*Q*. But since two dissimilar values of

*ℓ*and

*ℓ*′ differ by at least 1, this implies that for

*Q*≫1 the integral (35) will be vanishingly small unless

*ℓ*=

*ℓ*′. The narrow width of these sinc functions enables us to perform the integral approximately by replacing all slowly varying functions of u in the integrand by their values at

*u*=

*ℓ/Q*and extending the integration limits to ±∞. Use of the well known integral,

*Jℓℓ*′.

*B.2*. |

*ℓ*|,|

*ℓ*′|>

*Q*

*ℓ*,

*ℓ*′>

*Q*, we may simplify the product of the latter two sinc functions in Eq. (35) as

*Q*, its contribution to the integral (35) is of order 1/

*Q*when compared to the contribution of the first term. We may thus ignore the rapidly oscillating term, and the expression (35) reduces to the form,

*ℓ*-

*ℓ*′)

*π*=(-1)

*-*

^{ℓ}*′.*

^{ℓ}*B.3*. |

*ℓ*|<

*Q*, |

*ℓ*′|>

*Q*

*or*|

*ℓ*|>

*Q*, |ℓ′|<

*Q*

*ℓ*|<

*Q*but |

*ℓ*′|>

*Q*, then the sinc(

*Qu*-

*ℓ*) factor inside the integral (35) is peaked within the range of integral (35) and may be replaced by (1/

*Q*)

*δ*(

*u*-

*ℓ*/

*Q*) in the large

*Q*limit. This enables us to evaluate the integral as

*ℓ*-

*ℓ*′)=0 for all integral

*ℓ*,

*ℓ*′. The other sub-case, |

*ℓ*|>

*Q*,|

*ℓ*′|<

*Q*, is symmetrical, and the same result (39) holds once again.

*Q*becomes increasingly large and

*ℓ*,

*ℓ*′ take values well in the interior of their ranges considered here.

## Acknowledgments

## References and links

1. | S. Prasad, “Digital superresolution and the generalized sampling theorem,” J. Opt. Soc. Am. A |

2. | S. Prasad, “Digital and optical superresolution of low-resolution image sequences,” Unconventional Imaging Conference, 2007 Annual SPIE Meeting, San Diego, CA, Aug 25–29, 2007: Proc. SPIE |

3. | R. Gerchberg, “Superresolution through error energy reduction,” Opt. Acta |

4. | A. Papoulis, “A new algorithm in spectral analysis and band-limited extrapolation,” IEEE Trans. Circuits Syst. |

5. | S. Plevritis and A. Macovski, “Spectral extrapolation of spatially bounded images,” IEEE Trans. Medical Imaging |

6. | W. Richardson, “Bayesian-based iterative method of image restoration,” J. Opt. Soc. Am. |

7. | L. Lucy, “An iterative technique for the rectification of observed distributions,” Astron. J. |

8. | L. Shepp and Y. Vardi, “Maximum likelihood reconstruction in positron emission tomography,” IEEE Trans. Medical Imaging |

9. | T. Holmes, “Maximum-likelihood image restoration adapted for noncoherent optical imaging,” J. Opt. Soc. Am. A |

10. | T. Holmes and Y.-H. Liu, “Richardson-Lucy/maximum likelihood image restoration algorithm for fluorescence microscopy: further testing,” Appl. Opt. |

11. | J.-A. Conchello, “Superresolution and convergence properties of the expectation-maximization algorithm for maximum-likelihood deconvolution of incoherent images,” J. Opt. Soc. Am. A |

12. | E. Boukouvala and A. Lettington, “Restoration of astronomical images by an iterative superresolving algorithm,” Astron. Astrophys. |

13. | B. R. Frieden, “Image enhancement and restoration,” in |

14. | R. Narayan and R. Nityananda,“Maximum entropy image restoration in astronomy,” Ann. Rev. Astron. Astrophys. |

15. | J. Hogboom, “Aperture synthesis with a non-regular distribution of interferometer baselines,” Astron. Astrophys. Supp. |

16. | D. Fried, “Analysis of the CLEAN algorithm and implications for superresolution,” J. Opt. Soc. Am. A |

17. | P. Magain, F. Courbin, and S. Sohy, “Deconvolution with correct sampling,” Astrophys. J. |

18. | F. Pijpers, “Unbiased image reconstruction as an inverse problem,” Mon. Not. Roy. Astron. Soc. |

19. | R. Puetter and R. Hier, “Pixon sub-diffraction space imaging,” in Proc. SPIE |

20. | J. Starck, E. Pantin, and F. Murtagh, “Deconvolution in astronomy: A review,” Publ. Astron. Soc. Pacific |

21. | C. Rushforth and R. Harris, “Restoration, resolution, and noise,” J. Opt. Soc. Am. |

22. | G. Toraldo Di Francia, “Degrees of freedom of an image,” J. Opt. Soc. Am. |

23. | B. R. Frieden, “Evaluation, design, and extrapolation methods for optical signals based on the use of the prolate functions,” in |

24. | M. Bertero and E. R. Pike, “Resolution in diffraction-limited imaging, a singular value analysis: I. The case of coherent illumination,” Opt. Acta |

25. | M. Bertero and C. De Mol, “Superresolution by data inversion,” Progress in Optics |

26. | C. Matson and D. Tyler, “Primary and secondary superresolution by data inversion,” Opt. Express |

27. | D. Robinson and P. Milanfar, “Statistical performance analysis of superresolution,” IEEE Trans. Image Process. |

28. | P. Sementilli, B. Hunt, and M. Nader, “Analysis of of the limit of superresolution in incoherent imaging,” J. Opt. Soc. Am. A |

29. | H. Van Trees, |

**OCIS Codes**

(100.2000) Image processing : Digital image processing

(100.3010) Image processing : Image reconstruction techniques

(100.3190) Image processing : Inverse problems

(100.6640) Image processing : Superresolution

**ToC Category:**

Image Processing

**History**

Original Manuscript: September 30, 2009

Revised Manuscript: November 17, 2009

Manuscript Accepted: November 18, 2009

Published: December 3, 2009

**Citation**

Sudhakar Prasad and Xuan Luo, "Support-assisted optical superresolution of low-resolution image sequences: the one-dimensional problem," Opt. Express **17**, 23213-23233 (2009)

http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-17-25-23213

Sort: Year | Journal | Reset

### References

- S. Prasad, "Digital superresolution and the generalized sampling theorem," J. Opt. Soc. Am. A 24, 311-325 (2007) [CrossRef]
- S. Prasad, "Digital and optical superresolution of low-resolution image sequences," Proc. SPIE 6712, 67120E 1-11 (2007).
- R. Gerchberg, "Superresolution through error energy reduction," Opt. Acta 21, 709-721 (1974). [CrossRef]
- A. Papoulis, "A new algorithm in spectral analysis and band-limited extrapolation," IEEE Trans. Circuits Syst. CAS-22, 735-742 (1975). [CrossRef]
- S. Plevritis and A. Macovski, "Spectral extrapolation of spatially bounded images," IEEE Trans. Medical Imaging 14, 487-497 (1995). [CrossRef]
- W. Richardson, "Bayesian-based iterative method of image restoration," J. Opt. Soc. Am. 62, 55-59 (1972). [CrossRef]
- L. Lucy, "An iterative technique for the rectification of observed distributions," Astron. J. 79, 745-754 (1974) [CrossRef]
- L. Shepp and Y. Vardi, "Maximum likelihood reconstruction in positron emission tomography," IEEE Trans. Medical Imaging 1, 113-122 (1982). [CrossRef]
- T. Holmes, "Maximum-likelihood image restoration adapted for noncoherent optical imaging," J. Opt. Soc. Am. A 5, 666-673 (1988). [CrossRef]
- T. Holmes and Y.-H. Liu, "Richardson-Lucy/maximum likelihood image restoration algorithm for fluorescence microscopy: further testing," Appl. Opt. 28, 4930-4938 (1989). [CrossRef] [PubMed]
- J.-A. Conchello, "Superresolution and convergence properties of the expectation-maximization algorithm for maximum-likelihood deconvolution of incoherent images," J. Opt. Soc. Am. A 15, 2609-2619 (1998). [CrossRef]
- E. Boukouvala and A. Lettington, "Restoration of astronomical images by an iterative superresolving algorithm," Astron. Astrophys. 399, 807-811 (2003). [CrossRef]
- B. R. Frieden, "Image enhancement and restoration," in Picture Processing and Digital Filtering, T. Huang, ed., (Springer-Verlag, NY, 1975), 177-248.
- R. Narayan and R. Nityananda,"Maximum entropy image restoration in astronomy," Ann. Rev. Astron. Astrophys. 24, 127-170 (1986). [CrossRef]
- J. Hogboom, "Aperture synthesis with a non-regular distribution of interferometer baselines," Astron. Astrophys. Supp. 15, 417-426 (1974).
- D. Fried, "Analysis of the CLEAN algorithm and implications for superresolution," J. Opt. Soc. Am. A 12, 853-860 (1995). [CrossRef]
- P. Magain, F. Courbin, and S. Sohy, "Deconvolution with correct sampling," Astrophys. J. 494, 472-477 (1998). [CrossRef]
- F. Pijpers, "Unbiased image reconstruction as an inverse problem," Mon. Not. Roy. Astron. Soc. 307, 659-668 (1999). [CrossRef]
- R. Puetter and R. Hier, "Pixon sub-diffraction space imaging," Proc. SPIE 7094, 709405-709405-12 (2008).
- J. Starck, E. Pantin, and F. Murtagh, "Deconvolution in astronomy: A review," Publ. Astron. Soc. Pacific 114, 1051-1069 (2002). [CrossRef]
- C. Rushforth and R. Harris, "Restoration, resolution, and noise," J. Opt. Soc. Am. 58, 539-545 (1968). [CrossRef]
- G. Toraldo Di Francia, "Degrees of freedom of an image," J. Opt. Soc. Am. 59799-805 (1969).
- B. R. Frieden, "Evaluation, design, and extrapolation methods for optical signals based on the use of the prolate functions," in Progress in Optics, E. Wolf, ed., (North-Holland, Amsterdam, 1971), Vol. 9, 313-407.
- M. Bertero and E. R. Pike, "Resolution in diffraction-limited imaging, a singular value analysis: I. The case of coherent illumination," Opt. Acta 29, 727-746 (1982). [CrossRef]
- M. Bertero and C. De Mol, "Superresolution by data inversion," Progress in Optics 36, 129-178 (1996). [CrossRef]
- C. Matson and D. Tyler, "Primary and secondary superresolution by data inversion," Opt. Express 14, 456-473 (2006). [CrossRef] [PubMed]
- D. Robinson and P. Milanfar, "Statistical performance analysis of superresolution," IEEE Trans. Image Process. 15, 1413-1428 (2006). [CrossRef] [PubMed]
- P. Sementilli, B. Hunt, and M. Nader, "Analysis of of the limit of superresolution in incoherent imaging," J. Opt. Soc. Am. A 10, 2265-2276 (1993). [CrossRef]
- H. Van Trees, Detection, Estimation, and Modulation Theory (Wiley, New York, 1968).

## Cited By |
Alert me when this paper is cited |

OSA is able to provide readers links to articles that cite this paper by participating in CrossRef's Cited-By Linking service. CrossRef includes content from more than 3000 publishers and societies. In addition to listing OSA journal articles that cite this paper, citing articles from other participating publishers will also be listed.

« Previous Article | Next Article »

OSA is a member of CrossRef.