Super-resolution for imagery from integrated microgrid polarimeters

doi:10.1364/OE.19.012937

« Show journal navigation

Super-resolution for imagery from integrated microgrid polarimeters

Russell C. Hardie, Daniel A. LeMaster, and Bradley M. Ratliff »View Author Affiliations

Optics Express, Vol. 19, Issue 14, pp. 12937-12960 (2011)
http://dx.doi.org/10.1364/OE.19.012937

View Full Text Article

Acrobat PDF (1874 KB)

Journal Search
Article Lookup

Article Tools

Citations

Alert me when this article is cited
Export Citation/Save

Full Text
Article Info
References (16)
Cited By
Figures (15)
Metrics

Abstract

Imagery from microgrid polarimeters is obtained by using a mosaic of pixel-wise micropolarizers on a focal plane array (FPA). Each distinct polarization image is obtained by subsampling the full FPA image. Thus, the effective pixel pitch for each polarization channel is increased and the sampling frequency is decreased. As a result, aliasing artifacts from such undersampling can corrupt the true polarization content of the scene. Here we present the first multi-channel multi-frame super-resolution (SR) algorithms designed specifically for the problem of image restoration in microgrid polarization imagers. These SR algorithms can be used to address aliasing and other degradations, without sacrificing field of view or compromising optical resolution with an anti-aliasing filter. The new SR methods are designed to exploit correlation between the polarimetric channels. One of the new SR algorithms uses a form of regularized least squares and has an iterative solution. The other is based on the faster adaptive Wiener filter SR method. We demonstrate that the new multi-channel SR algorithms are capable of providing significant enhancement of polarimetric imagery and that they outperform their independent channel counterparts.

1. Introduction

Passive polarimetric imaging has emerged in recent years as a useful sensing modality for detection and recognition [1

J. S. Tyo, D. L. Goldstein, D. B. Chenault, and J. A. Shaw, “Review of passive imaging polarimetry for remote sensing applications,” Appl. Opt. 45(22), 5453–5469 (2006). [CrossRef] [PubMed]

]. The added information provided by polarization can aid in clutter reduction and object characterization. Imagery from microgrid polarimeters is obtained by incorporating a mosaic of pixel-wise micropolarizers on a focal plane array (FPA). Such imagers are akin to Bayer pattern color cameras, except the mosaic pattern is for polarization state, rather than color. With a typical microgrid polarimetric sensor, four polarimetric state images are acquired using a 2 × 2 pattern of polarizers. In this manner, each polarization image is obtained by subsampling the full FPA image by 2 in each dimension. Thus, the effective pixel pitch for each channel is twice that of the native FPA (1/2 the sampling frequency). Because of the desire for a wide field of view and other factors, many imaging systems are equipped with optics that allow some undersampling for a given FPA. By employing the microgrid with the same optics and FPA, undersampling and resulting aliasing artifacts are increased [2

B. M. Ratliff, C. F. LaCasse, and J. S. Tyo, “Interpolation strategies for reducing IFOV artifacts in microgrid polarimeter imagery,” Opt. Express 17(11), 9112–9125 (2009). [CrossRef] [PubMed]

, 3

J. S. Tyo, C. F. LaCasse, and B. M. Ratliff, “Total elimination of sampling errors in polarization imagery obtained with integrated microgrid polarimeters,” Opt. Lett. 34(20), 3187–3189 (2009). [CrossRef] [PubMed]

Multi-frame super-resolution (SR) enhancement methods have been widely demonstrated to improve the effective sampling rate of undersampled imagers [4

S. C. Park, M. K. Park, and M. G. Kang, “Super-resolution image reconstruction: a technical overview,” IEEE Signal Process. Mag. 20(3), 21–36 (2003). [CrossRef]

]. By using multiple low-resolution (LR) frames from video, these SR algorithms exploit sub-pixel motion between frames to effectively enhance the sampling rate for the native FPA. By doing so, these SR methods effectively trade temporal resolution for spatial resolution. In addition to enhancing the sampling rate and reducing aliasing, the SR algorithms typically perform restoration to reduce noise and blur. Generally speaking, iterative image reconstruction based SR methods are among the best performing algorithms, but come with a relatively high computational cost [4

S. C. Park, M. K. Park, and M. G. Kang, “Super-resolution image reconstruction: a technical overview,” IEEE Signal Process. Mag. 20(3), 21–36 (2003). [CrossRef]

]. The simplest SR algorithms, in a computational sense, are the interpolation-restoration approaches [4

S. C. Park, M. K. Park, and M. G. Kang, “Super-resolution image reconstruction: a technical overview,” IEEE Signal Process. Mag. 20(3), 21–36 (2003). [CrossRef]

]. Here we shall consider these two classes of SR algorithms as applied to polarimetric imagery. Because of their ability to help overcome undersampling, we believe multi-frame SR algorithms are a powerful match for microgrid polarimeters. To the best of our knowledge, the only publications to date exploring this connection were the result of preliminary work led by two of the current coauthors [5

B. M. Ratliff, J. S. Tyo, W. T. Black, and C. F. LaCasse, “Exploiting motion-based redundancy to enhance microgrid polarimeter imagery,” Proc. SPIE 7461, 74610K (2009). [CrossRef]

, 6

D. A. Lemaster, “Resolution enhancement by image fusion for microgrid polarization imagers,” in IEEE Aerospace Conference, Big Sky, MT (2010).

]. We believe the current paper represents the first comprehensive treatment of the subject. In [5

B. M. Ratliff, J. S. Tyo, W. T. Black, and C. F. LaCasse, “Exploiting motion-based redundancy to enhance microgrid polarimeter imagery,” Proc. SPIE 7461, 74610K (2009). [CrossRef]

], Ratliff et al used the iterative multi-frame SR algorithm in [7

R. C. Hardie, K. J. Barnard, J. G. Bognar, E. E. Armstrong, and E. A. Watson, “High resolution image reconstruction from a sequence of rotated and translated frames and its application to an infrared imaging system,” Opt. Eng. 37(1), 247–260 (1998). [CrossRef]

] applied to the four demosaiced polarimetric channels independently. LeMaster explored a simple and fast multi-frame SR algorithm in [6

D. A. Lemaster, “Resolution enhancement by image fusion for microgrid polarization imagers,” in IEEE Aerospace Conference, Big Sky, MT (2010).

] that registered and interlaced samples from multiple frames, also using independent channel processing. The work in [5

B. M. Ratliff, J. S. Tyo, W. T. Black, and C. F. LaCasse, “Exploiting motion-based redundancy to enhance microgrid polarimeter imagery,” Proc. SPIE 7461, 74610K (2009). [CrossRef]

, 6

D. A. Lemaster, “Resolution enhancement by image fusion for microgrid polarization imagers,” in IEEE Aerospace Conference, Big Sky, MT (2010).

] demonstrated that SR can be used to enhance the resolution of polarimetric data from a microgrid imager, but neither fully exploited the redundant information available in the microgrid array.

Here we present two new multi-channel multi-frame SR algorithms specifically designed for processing imagery from microgrid polarimeters. The new SR methods are designed to provide improved performance by exploiting correlation between the polarimetric channels and by using a detailed observation model that includes an accurate point spread function (PSF). One of the two new SR algorithms presented here uses a form of regularized least squares (RLS) related to that in [7

] and has an iterative solution. The proposed RLS algorithm is novel in that here we jointly estimate all of the polarimetric images and use a new regularization function that exploits correlation among the polarimetric channels. The other proposed algorithm is a new type of adaptive wiener filter (AWF) SR method [8

R. Hardie, “A fast image super-resolution algorithm using an adaptive Wiener filter,” IEEE Trans. Image Process. 16(12), 2953–2964 (2007). [CrossRef] [PubMed]

] designed for microgrid image data. The microgrid AWF method leads to a faster implementation that is more readily suitable to real-time processing with appropriate hardware than the RLS method. We study the performance of the two new microgrid SR methods using imagery from a long-wave infrared (LWIR) microgrid polarimeter. We also present results using simulated microgrid data derived from a visible polarimetric imaging system. These simulated data are used to allow for a new quantitative performance analysis. Our results demonstrate that there can be significant benefits to using SR with microgrid polarimeters. We also show how multichannel SR processing that exploits the correlation between the polarimetric channels yields improved performance over independent channel processing.

The remainder of this paper is organized as follows. In Section 2, background on polarimetric imagery is briefly reviewed and the observation models are presented. The new microgrid RLS SR algorithm is presented in Section 3 and the new microgrid AWF SR method is presented in Section 4. Experimental results are presented in Section 5 to illustrate the efficacy of the proposed methods. Finally, we offer conclusions in Section 6.

2. Observation model

The sensor we are focusing on is equipped with four polarization analyzers oriented at angles θ ₁ = 0°, θ ₂ = 45°, θ ₃ = 90°, and θ ₄ = 135°. The microgrid configuration for acquiring these four polarimetric images is shown in Fig. 1. As shown, the microgrid image can be deinterlaced to produce the four individual polarimetric images. The choice of angles is designed to make the computation of Stokes vector images convenient [1

J. S. Tyo, D. L. Goldstein, D. B. Chenault, and J. A. Shaw, “Review of passive imaging polarimetry for remote sensing applications,” Appl. Opt. 45(22), 5453–5469 (2006). [CrossRef] [PubMed]

, 9

R. M. A. Azzam and N. M. Bashara, Ellipsometry and Polarized Light (North-Holland, 1977).

]. In general, let the ideally-sampled high-resolution (HR) polarimetric image, with P channels and N pixels per channel, be represented in lexicographical notation as p = [p ₁, p ₂,..., p_NP ] ^T . Breaking this down into individual polarimetric channels, we get

p = {[p_{θ_{1}}^{T}, p_{θ_{2}}^{T}, \dots, p_{θ_{P}}^{T}]}^{T}

, where p _{θ_i} = [p_{θ_i} _,1, p_{θ_i} ,₂,..., p_{θ_i,N} ] ^T . For such a four-angle sensor, the first three Stokes images [1

J. S. Tyo, D. L. Goldstein, D. B. Chenault, and J. A. Shaw, “Review of passive imaging polarimetry for remote sensing applications,” Appl. Opt. 45(22), 5453–5469 (2006). [CrossRef] [PubMed]

] can be estimated as

\begin{array}{l} s_{0, j} & = & (p_{0, j} + p_{45, j} + p_{90, j} + p_{135, j}) / 2 \\ s_{1, j} & = & p_{0, j} - p_{90, j} \\ s_{2, j} & = & p_{45, j} - p_{135, j} \end{array} .

(1)

The s ₀ term is the total light intensity, and s ₁ and s ₂ present linear polarization content. Note that since we have four polarimetric states being measured and three Stokes images relating to linear polarization, there is some redundancy. That is, ideally we should have s _0, _j = p _0, _j + p _90, _j = p _45, _j + p _135, _j . Later we shall seek to exploit this redundancy in our SR algorithms. Another useful statistic, related to the Stokes images, is the degree of linear polarization (DoLP) [1

J. S. Tyo, D. L. Goldstein, D. B. Chenault, and J. A. Shaw, “Review of passive imaging polarimetry for remote sensing applications,” Appl. Opt. 45(22), 5453–5469 (2006). [CrossRef] [PubMed]

] given by

{DoLP}_{j} = \frac{\sqrt{s_{1, j}^{2} + s_{2, j}^{2}}}{s_{0, j}} .

(2)

Fig. 1 Microgrid pattern for the polarimetric imager considered here. Four polarization images are acquired using a 2 × 2 repeating pattern of micropolarizers on the FPA.

For multi-frame SR, let the number of observed LR microgrid frames be K. Let the LR observed pixels from frame k be expressed in lexicographical notation as g(k), where k = 1, 2,...,K. Breaking down each observed frame by polarimetric channel we get

g (k) = {[g_{θ_{1}}^{T} (k), g_{θ_{2}}^{T} (k), \dots, g_{θ_{P}}^{T} (k)]}^{T}

, where g _{θ_i} (k) = [g_{θ_i} ,₁ (k), g_{θ_i} ,₂ (k),..., g_{θ_i} , _M (k)] ^T . Let the complete set of all observed pixels in lexicographical notation be specified by g = [g ^T (1), g ^T (2),...,g ^T (K)] ^T = [g ₁, g ₂,..., g_KMP ] ^T .

A block diagram illustrating the observation model relating p to the observed frame g(k) is presented in Fig. 2. The ideal polarimatric data in p is assumed to be sampled above the Nyquist rate. The motion between frames is parameterized with the vector s(k). Next, the model includes a linear shift invariant blurring using the discrete PSF for the optical system. This is followed by down-sampling by a factor of L. Finally, the down-sampled images are mosaiced into the pattern in Fig. 1 and noise is introduced to obtain g(k). This observation model can be mathematically expressed such that each observed pixel is a weighted sum of the ideal pixels in p plus noise. The particular weights depend on the system PSF, the motion between frames, and the microgrid layout. In particular, the LR pixels can be expressed in terms of ideal HR pixels as follows

g_{θ, i} (k) = \sum_{j = 1}^{N} h_{θ, i, j} (k) p_{θ, j} + n_{θ, i} (k),

(3)

where h_θ _,i,j (k) is the contribution of p_θ _,j to g_θ _,i (k) and n_θ _,i (k) is the noise term for g_θ _,i (k). This model can be express compactly for all frames jointly as

g = Hp + n,

(4)

where H is a KMP × NP matrix and n is a KMP × 1 vector of noise samples.

Fig. 2 Block diagram illustrating the observation model relating p to the observed frames g(k), for k = 1,2,...,K.

The observation model in Eqs. (3) and (4) will be used for the RLS SR algorithm. However, to aid in the development of the fast AWF algorithm here, we consider an alternate observation model based on that in [8

R. Hardie, “A fast image super-resolution algorithm using an adaptive Wiener filter,” IEEE Trans. Image Process. 16(12), 2953–2964 (2007). [CrossRef] [PubMed]

]. The alternate observation model is shown in Fig. 3. Here the motion model is commuted with the PSF and integrated into what becomes a nonuniform sampling operator. This alternative observation model is equivalent to the model in Fig. 2 for translational motion. It is also equivalent for rotational motion when the PSF is circularly symmetric [8

R. Hardie, “A fast image super-resolution algorithm using an adaptive Wiener filter,” IEEE Trans. Image Process. 16(12), 2953–2964 (2007). [CrossRef] [PubMed]

]. The power of this model is that it leads to the computationally simple interpolation-restoration SR algorithms [8

R. Hardie, “A fast image super-resolution algorithm using an adaptive Wiener filter,” IEEE Trans. Image Process. 16(12), 2953–2964 (2007). [CrossRef] [PubMed]

]. Based on this model, nonuniform interpolation can be used to recover noisy, but uniformly spaced, samples of f (x,y,θ_i ) from g. Then, restoration is used to estimate the corresponding samples of p(x,y, θ_i ).

Fig. 3 Alternative observation model valid for translational motion (and rotational motion when the PSF is circularly symmetric). Here the motion model is commuted with the PSF and integrated into what becomes a nonuniform sampling operator.

Let us now consider the modeling of the system PSF. We model the optical transfer function (OTF) with three components

H (u, v) = H_{dif} (u, v) H_{abr} (u, v) H_{det} (u, v),

(5)

where u and v are the horizontal and vertical spatial frequencies in cycles per millimeter, H _dif(u, v) is from the diffraction-limited optics, H _abr(u, v) models optical aberrations, and detector integration is included as H _det(u, v). The OTF from diffraction-limited optics with a circular pupil function is given by [10

J. Goodman, Introduction to Fourier Optics (McGraw-Hill, 1968).

]

H_{dif} (u, v) = {\begin{matrix} \frac{2}{π} [{cos}^{- 1} (ρ / ρ_{c}) - (ρ / ρ_{c}) \sqrt{1 - {(ρ / ρ_{c})}^{2}}] & for ρ < ρ_{c} \\ 0 & else \end{matrix},

(6)

where

ρ = \sqrt{u^{2} + v^{2}}

and

ρ_{c} = \frac{1}{λ 𝒩}

, and 𝒩 is the F/# of the optics. Since a purely diffraction-limited optical system is perhaps overly idealistic, we also include the following aberration OTF [11

R. R. Shannon, “Aberrations and their effect on images,” Proc. SPIE 531, 27–37 (1985).

]

H_{abr} (u, v) = {\begin{matrix} 1 - {(W_{RMS} / 0.18)}^{2} {(1 - 4 (ρ / ρ_{c} - 0.5))}^{2} & for ρ < ρ_{c} \\ 0 & else \end{matrix},

(7)

where W_RMS is the root mean square wavefront error [11

R. R. Shannon, “Aberrations and their effect on images,” Proc. SPIE 531, 27–37 (1985).

]. The detector component of the OTF is found as the Fourier transform of a mask with the geometry of the active area of an individual detector in the FPA. Finally, the system PSF is the inverse Fourier transform of the overall OTF in Eq. (5).

The microgrid polarimetric imager used here has the following specifications. It has a spectral bandwidth of λ = 7.8 – 9.8 μm. We use a wavelength of λ = 9 μm for our PSF model. The system uses F/2 optics and has a full FPA pixel pitch of 25 μm with approximately 100% fill factor rectangular detectors. Note that for this microgrid imager, the pixel pitch for pixels of like polarization is actually 50 μm. Thus, the effective active area for the like-polarization sub-array detectors is modeled as spanning half the 50 μm like-polarization pixel pitch in each dimension. This equates to a 25% fill factor (taking into account the skipping of detectors in the FPA as shown in Fig. 1). For abberation modeling, we use W_RMS = 1/14 [11

R. R. Shannon, “Aberrations and their effect on images,” Proc. SPIE 531, 27–37 (1985).

]. Cross sections of the relevant 2-D modulation transfer functions (MTFs) are shown in Fig. 4(a) for this system. The continuous PSF is shown in Fig. 4(b). The discrete impulse invariant PSF, for use in the observation model operating on the HR image, can be found by sampling the continuous PSF with a sampling period of 50/L μm.

Fig. 4 (a) Cross sections of theoretical 2-D MTF and its components for the microgrid polarimetric imager used here. (b) Theoretical continuous PSF for L = 4.

Note that the Nyquist or folding frequency (i.e., 1/2 sampling frequency) for a single channel is 10 cycles/mm, as shown in green in Fig. 4(a). This is well below the 55.56 cycles/mm cut-off frequency of the optics. Any scene frequency content above the Nyquist frequency is “folded” into lower frequencies, creating aliasing artifacts. Thus, this system is likely to exhibit significant aliasing without SR. Another point to note from the overall MTF is that it has a low pass blurring effect. However, the blurring cannot be combatted with a linear shift-invariant filter without first (or jointly) addressing the undersampling issue. Otherwise, the aliasing artifacts will simply be amplified. Based on the overall system MTF shown in Fig. 4(a), we believe that L = 4 is a practical choice for SR processing for this sensor. The effective sampling frequency is then 80 cycles/mm, with a Nyquist frequency of 40 cycles/mm as shown in yellow in Fig. 4(a). In addition to aliasing reduction, we also hope to restore the imagery from the blurring effects of the overall MTF.

3. RLS SR for microgrid polarimeters

Based on the observation model depicted in Fig. 2 and detailed in Eq. (3), we propose an RLS estimator for p. The output of the RLS algorithm is given by

\hat{p} = \underset{p}{argmin} C (p) .

(8)

The cost function to be minimized is given by

\begin{matrix} C (p) & = & \frac{1}{σ_{n}^{2}} \sum_{θ} \sum_{k = 1}^{K} \sum_{i = 1}^{M} {(g_{θ, i} (k) - \sum_{j = 1}^{N} h_{θ, i, j} (k) p_{θ, j})}^{2} + \sum_{m = 0}^{2} \frac{1}{σ_{s_{m}}^{2}} \sum_{i = 1}^{N} {(\sum_{j = 1}^{N} w_{i, j} s_{m, j})}^{2} . \\ + \frac{1}{σ_{r}^{2}} \sum_{j = 1}^{N} {(p_{0, j} + p_{90, j} - p_{45, j} - p_{135, j})}^{2} \end{matrix}

(9)

The first term in Eq. (9) is the main term that is minimized when p projected through observation model matches the observed data g. Since this term alone may be insufficient to provide a unique solution, the second and third terms act as a regularization function. The second term acts as a smoothness constraint on the Stokes images. For this smoothness term we use a discrete Laplacian operator given by

w_{i, j} = {\begin{matrix} 1 & i = j \\ - 1 / 4 & j is cardinal neighbor of i \end{matrix} .

(10)

Finally, the third term in Eq. (9) exploits the redundancy in the four polarimetric measurements. This term is minimized when p _0, _j + p _90, _j = p _45, _j + p _135, _j . For an ideal noise-free sensor, with no cross-talk between polarization channels, this equality should be true [1

J. S. Tyo, D. L. Goldstein, D. B. Chenault, and J. A. Shaw, “Review of passive imaging polarimetry for remote sensing applications,” Appl. Opt. 45(22), 5453–5469 (2006). [CrossRef] [PubMed]

] and this cost function term would go to zero. However, in practice, noise and non-ideal sensor characteristics make this a small but non-zero term. Note that each term in Eq. (9) has a corresponding scaling parameter. The value of these constants govern the relative impact of each term on the cost function, and ultimately the final estimate. Specifically, the first term is scaled by

1 / σ_{n}^{2}

, which would typically be related to the expected noise variance in the observation model. With high noise level, we would expect a larger error in the first cost function term, even for the true image as the estimate. Thus, employing a larger

σ_{n}^{2}

reduces the cost associated with this expected higher error. The second term has scaling parameters of

1 / σ_{s_{m}}^{2}

, for m = 0, 1, 2. These reflect our a priori assumptions about the variance of the Laplacian of the Stokes images as defined using Eq. (10). Finally, the third term is scaled by

1 / σ_{r}^{2}

, which reflects our expectation for the variance of the difference between p _0, _j + p _90, _j and p _45, _j + p _135, _j . Note that it is only the relative values of these scaling constants that matter for the algorithm.

While this RLS approach is based on that in [7

], what makes Eq. (9) novel are the regularization terms. The RLS algorithm in [7

] uses only one regularization term, a smoothness term for the single channel data. Here we have multichannel data with physical constraints, like the redundancy exploited in the third term in Eq. (9). Note also, that the smoothness terms are applied here to the Stokes images, and not the original polarimetric channels. This is done for several reasons. First, the Stokes images may be a final product for such an imaging system. By having smoothness terms applied here, we better ensure that this product will be free from high frequency artifacts. Furthermore, {s _1, _j } and {s _2, _j } are known to generally have relatively low energy (i.e., low variance). Thus, applying a term favoring small high-frequency energy in these channels is sensible. The Stokes image {s _0, _j } is the total intensity image. This should contain most of the energy. One may choose to use the smoothness constraint on this term minimally by employing a high σ_s _₀. The parameters σ_s _₁ and σ_s _₂ may be smaller, because of the expected lower energy in these images. Finally, because the Stokes terms interrelate the different polarization channels, it allows us to exploit their distinct spatial sampling positions better than if we applied a smoothness constraint on each channel independently. The approach of applying smoothness regularization to the Stokes images in polarimetric data may be viewed as akin to using luminance and chrominance coordinates in color processing [12

T. Gotoh and M. Okutomi, “Direct super-resolution and registration using raw CFA images,” in IEEE Conference on Computer Vision and Pattern Recognition , vol. 2, pp. 600–607 (Los Alamitos, CA, 2004).

–14

L. Condat, “A generic variational approach for demosaicking from an arbitrary color filter array,” in Proceedings of IEEE ICIP , pp. 1625–1628 (2009).

]. It is interesting to note that the RLS estimate can be viewed as a maximum a posteriori (MAP) estimator in the case of independent and identically distributed Gaussian noise [15

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 Trans. Image Process. 6(12), 1621–1633 (1997). [CrossRef] [PubMed]

]. In the MAP framework, the quadratic nature of the regularization terms in Eq. (9) can be seen as resulting from a Gaussian prior pdf.

To minimize the cost function in Eq. (9), one may use a gradient descent or conjugate gradient algorithm similar to that presented in [7

]. Here we develop the somewhat simpler gradient descent approach. We initialize the RLS estimation process by registering all K frames to a common reference. A convenient choice of reference is the most recent observed frame. Note that we need to register the frames to subpixel accuracy. We have found that we get better registration performance by deinterlacing the channels, applying the registration algorithm in [7

], and then averaging the registration parameters obtained from the four channels. With this and the PSF from Section 2, we have a complete observation model. Next, we take the reference image and use single frame interpolation to obtain an initial estimate of p, denoted p̂ ⁰. To improve this estimate using gradient descent, we need the gradient of Eq. (9) with respect to p, denoted ∇ _p C(p). This gradient is composed of the partial derivative of the cost function with respect to each of the HR samples in p. This can be expressed as

\nabla_{p} C (p) = {[\nabla_{p_{0}} C {(p)}^{T}, \nabla_{p_{45}} C {(p)}^{T}, \nabla_{p_{90}} C {(p)}^{T}, \nabla_{p_{135}} C {(p)}^{T}]}^{T},

(11)

where

\nabla_{p θ} C (p) = {[\frac{\partial C (p)}{\partial p_{θ, 1}}, \frac{\partial C (p)}{\partial p_{θ, 2}}, \dots, \frac{\partial C (p)}{\partial p_{θ, N}},]}^{T} .

(12)

The partial derivatives making up the gradient can be computed from Eq. (9), yielding

\frac{\partial C (p)}{\partial p_{0, q}} = H (0, q) + S (0, q) + S (1, q) + R (q),

(13)

\frac{\partial C (p)}{\partial p_{45, q}} = H (45, q) + S (0, q) + S (2, q) - R (q),

(14)

\frac{\partial C (p)}{\partial p_{90, q}} = H (90, q) + S (0, q) - S (1, q) + R (q),

(15)

\frac{\partial C (p)}{\partial p_{135, q}} = H (135, q) + S (0, q) - S (2, q) - R (q),

(16)

where the observation model gradient term is given by

H (θ, q) = \frac{2}{σ_{n}^{2}} \sum_{k = 1}^{K} \sum_{i = 1}^{M} h_{θ, i, q} (k) e_{θ, i} (k),

(17)

and the model error component is

e_{θ, i} (k) = \sum_{j = 1}^{N} h_{θ, i, j} (k) p_{θ, j} - g_{θ, i} (k) .

(18)

The Stokes gradient term is given by

S (m, q) = \frac{2}{σ_{s_{m}}^{2}} \sum_{i = 1}^{N} w_{i, q} (\sum_{j = 1}^{N} w_{i, j} s_{m, j}) .

(19)

Finally, the redundancy gradient term is

R (q) = \frac{2}{σ_{r}^{2}} (p_{0, q} + p_{90, q} - p_{45, q} - p_{135, q}) .

(20)

Given the gradient, the gradient descent updates are computed as follows

{\hat{p}}^{n} = {\hat{p}}^{n - 1} + ɛ^{n} \nabla_{p} C (p) |_{p = {\hat{p}}^{n - 1}},

(21)

where ɛⁿ is the step size for iteration n. The iterations may be performed a predetermined number of times, or may be stopped with a halting criterion, such as ||p̂ ⁿ – p̂ ⁿ ⁻¹|| < T, where T is a threshold.

The optimum step size at each iteration can be found by optimizing C(p̂ ⁿ ) as a function of ɛⁿ , given p̂ ⁿ ⁻¹. In particular, we set the derivative of C(p̂ ⁿ ) with respect to ɛ ⁿ equal to zero, and solve for ɛⁿ . This calculation results in an optimum step size of the form

ɛ^{n} = \frac{N_{1} + N_{2} + N_{3}}{D_{1} + D_{2} + D_{3}} .

(22)

The individual terms in the optimum step size are given by

N_{1} = \frac{1}{σ_{n}^{2}} \sum_{θ} \sum_{k = 1}^{K} \sum_{i = 1}^{M} (\sum_{j = 1}^{N} h_{θ, i, j} (k) γ_{θ, j}) e_{θ, i} (k),

(23)

N_{2} = \sum_{m = 0}^{2} \frac{1}{σ_{s_{m}}^{2}} \sum_{i = 1}^{N} (\sum_{j = 1}^{N} w_{i, j} s_{m, j}) (\sum_{j = 1}^{N} w_{i, j} S_{m, j}),

(24)

N_{3} = \frac{1}{σ_{r}^{2}} \sum_{j = 1}^{N} (p_{0, j} + p_{90, j} - p_{45, j} - p_{135, j}) (γ_{0, j} + γ_{90, j} - γ_{45, j} - γ_{135, j}),

(25)

D_{1} = \frac{1}{σ_{n}^{2}} \sum_{θ} \sum_{k = 1}^{K} \sum_{i = 1}^{M} {(\sum_{j = 1}^{N} h_{θ, i, j} (k) γ_{θ, j})}^{2},

(26)

D_{2} = \sum_{m = 0}^{2} \frac{1}{σ_{s_{m}}^{2}} \sum_{i = 1}^{N} {(\sum_{j = 1}^{N} w_{i, j} S_{m, j})}^{2},

(27)

D_{3} = \frac{1}{σ_{r}^{2}} \sum_{j = 1}^{N} {(γ_{0, j} + γ_{90, j} - γ_{45, j} - λ_{135, j})}^{2},

(28)

where

γ_{θ, j} = \frac{\partial C (p)}{\partial p_{θ, j}} |_{p = {\hat{p}}^{n}}

(29)

and

\begin{array}{l} S_{0, j} & = & (γ_{0, j} + γ_{45, j} + γ_{90, j} + γ_{135, j}) / 2 \\ S_{1, j} & = & γ_{0, j} - γ_{90, j} \\ S_{2, j} & = & γ_{45, j} - γ_{135, j} \end{array} .

(30)

Computing the optimum step size in this fashion tends to be faster than performing a search. A final point to note about the RLS microgrid SR method is that we can easily treat “bad” pixels by simply setting the error, e_θ _,i (k) as defined in Eq. (18), to zero for any θ, i, and k for which g_θ _,i (k) is known to be a bad pixel value. This is convenient as it obviates the need for a separate bad pixel correction algorithm, even when using a single frame. This is how we address bad pixels for the RLS SR results presented in Section 5. Other bad pixel replacement strategies that exploit microgrid redundancy have been proposed in [16

B. M. Ratliff, J. S. Tyo, J. K. Boger, W. T. Black, D. L. Bowers, and M. P. Fetrow, “Dead pixel replacement in LWIR microgrid polarimeters,” Opt. Express 15(12), 7596–7609 (2007). [CrossRef] [PubMed]

4. AWF SR for microgrid polarimeters

Unlike the RLS, the AWF SR method [8

R. Hardie, “A fast image super-resolution algorithm using an adaptive Wiener filter,” IEEE Trans. Image Process. 16(12), 2953–2964 (2007). [CrossRef] [PubMed]

] is non-iterative and has a potentially faster implementation than the RLS. Here we propose a new type of the AWF SR algorithm [8

R. Hardie, “A fast image super-resolution algorithm using an adaptive Wiener filter,” IEEE Trans. Image Process. 16(12), 2953–2964 (2007). [CrossRef] [PubMed]

] for use with a microgrid imager. As in [8

R. Hardie, “A fast image super-resolution algorithm using an adaptive Wiener filter,” IEEE Trans. Image Process. 16(12), 2953–2964 (2007). [CrossRef] [PubMed]

], the microgrid AWF uses a moving observation window on the HR grid as shown in Fig. 5. The observation window is shown for L = 4 and is populated with pixels from a single frame for simplicity. With multiple frames of data with motion between them, the observation window will contain additional sets of pixels from each polarization state, positioned according to the motion. In general, let the observation window span W_x HR pixel spacings in the horizontal direction and W_y HR pixel spacings in the vertical direction. All of the LR pixels of any polarization that lie within the span of this observation window are placed into an observation vector g _i = [g_i _,1, g_i _,2,..., g_{i,J_i} ] ^T , where i is the window positional index and J_i is the number of LR pixels within the window span. The samples within the observation window will be used to estimate the HR pixels within the generally smaller D_x × D_y sample estimation window at the center of the observation window, as shown in Fig. 5. Thus, the AWF SR method is effectively performing nonuniform interpolation and restoration simultaneously. One of the main differences here, versus the AWF in [8

R. Hardie, “A fast image super-resolution algorithm using an adaptive Wiener filter,” IEEE Trans. Image Process. 16(12), 2953–2964 (2007). [CrossRef] [PubMed]

], is that the observation window contains LR pixels from multiple polarizations. Furthermore, at each HR pixel location, we must estimate all of the polarization states, not just a single value.

Fig. 5 Observation window for the microgrid AWF algorithm for polarimetric data with L = 4. For simplicity, the window is shown populated with pixels from only a single frame.

The HR pixel estimates for the samples within the estimation window are obtained using a weighted sum of the LR pixels in the observation window. This is expressed as

{\hat{d}}_{i} = W_{i}^{T} g_{i},

(31)

where d̂ _i = [d̂_i _,1, d̂_i _,2,..., d̂_{i,D_xD_yP} ] ^T and W _i is a J_i × D_xD_yP matrix of weights. Each column of W _i contains weights used to estimate the value of one particular HR pixel at one polarization state inside the estimation window. The observation window moves across the HR grid in a raster scan fashion stepping by increments of D_x and D_y in the horizontal and vertical directions, respectively (non-overlapping estimation sub-windows).

We wish to employ the minimum mean squared error (MSE) filter weights in Eq. (31). The additional challenge we face when processing microgrid data is that we must not only consider the spatial arrangement of samples in a given observation window, but also their polarization state. Given a correlation model for the observed samples, the minimum MSE weights [8

R. Hardie, “A fast image super-resolution algorithm using an adaptive Wiener filter,” IEEE Trans. Image Process. 16(12), 2953–2964 (2007). [CrossRef] [PubMed]

] are given by

W_{i} = R_{i}^{- 1} P_{i},

(32)

where

R_{i} = E {g_{i} g_{i}^{T}}

is the autocorrelation matrix for the observation vector and

P_{i} = E {g_{i} d_{i}^{T}}

is the cross-correlation between the true desired vector d _i and the observation vector. The goal now is to provide the statistics in R _i and P _i . As in [8

R. Hardie, “A fast image super-resolution algorithm using an adaptive Wiener filter,” IEEE Trans. Image Process. 16(12), 2953–2964 (2007). [CrossRef] [PubMed]

], let f _i be the noise-free version of g _i , such that g _i = f _i + n _i , where n _i is the random noise vector. We shall assume that the noise vector is zero-mean with independent and identically distributed elements of variance

σ_{n}^{2}

. In this case, the autocorrelation matrix for the observation vector is given by

R_{i} = E {g_{i} g_{i}^{T}} = E {f_{i} f_{i}^{T}} + σ_{n}^{2} I

(33)

and cross-correlation matrix is

P_{i} = E {g_{i} d_{i}^{T}} = E {f_{i} d_{i}^{T}} .

(34)

To fill the matrices in Eqs. (33) and (34), we need a statistical model to provide the needed correlations. Consider a correlation model, specifically designed for microgrid data, that incorporates knowledge of the spatial separation between samples, and the polarization state of each. For the desired data, we denote this correlation function as r_{d_θd_ϕ} (x,y). The variables θ and ϕ represent the polarization angles of the two samples involved, and x,y represents the spatial separation. Here we propose the following model

r_{d_{θ} d ϕ} (x, y) = σ_{θ, ϕ}^{2} ρ_{θ, ϕ}^{\sqrt{x^{2} + y^{2}}},

(35)

where

σ_{θ, ϕ}^{2}

controls correlation between samples at the same spatial location and ρ_θ,ϕ controls the decay of the correlation with separation distance. We have obtained useful results using ρ_θ,ϕ = ρ for all θ, ϕ, and

σ_{θ, ϕ}^{2} = {\begin{matrix} σ_{p}^{2} & θ = ϕ \\ α σ_{p}^{2} & otherwise \end{matrix} .

(36)

Note that if α = 0, we assume no correlation between samples of different polarizations regardless of their spatial proximity. This gives rise to an algorithm essentially equivalent to deinterlacing the microgrid data and applying independent AWF SR [8

R. Hardie, “A fast image super-resolution algorithm using an adaptive Wiener filter,” IEEE Trans. Image Process. 16(12), 2953–2964 (2007). [CrossRef] [PubMed]

] to each channel. In most cases, this is not the best approach, as there is often significant correction across polarization channels. By letting α = 1, this is equivalent to applying a single AWF to the mosaic microgrid data as if it were a standard FPA and each pixel represents the same physical quantity. In that case, spatial separation is the only factor altering the correlation model. This is also not likely the best approach as samples from the same polarization generally do exhibit higher correlation than those of different polarizations. Thus, one would want to select α in the range 0 ≤ α ≤ 1 based on expected scene content. A relatively high α is appropriate when we have weakly polarized data, like we may see in LWIR. A lower α may be used for more highly polarized data with a high average DoLP.

Given r_{d_θd_ϕ} (x,y), it can be shown that the cross-correlation function between p(x,y, θ) and f (x,y, ϕ), as shown in Fig. 3, can be expressed as

r_{d_{θ} f_{ϕ}} (x, y) = r_{d_{θ} d_{ϕ}} (x, y) * h (x, y) .

(37)

The cross-correlation between f (x,y, θ) and f (x,y,ϕ) is given by

r_{f_{θ} f_{ϕ}} (x, y) = r_{d_{θ} d_{ϕ}} (x, y) * h (x, y) * h (- x, - y) .

(38)

By registering the LR frames, the spatial locations and corresponding polarizations of all the LR pixels that make up g _i , which are the same for f _i , are known. Thus, the horizontal and vertical displacements between these samples can be computed. Evaluating Eq. (38) using the displacements and the corresponding polarizations allows us to fill the correlation matrix

E {f_{i} f_{i}^{T}}

. Next, R _i can be found from Eq. (33), given the noise variance. Similarly, we compute the displacements between the samples in the estimation window and the LR pixels in the observation window. With these displacements and the corresponding polarizations, we evaluate Eq. (37) and complete P _i as given in Eq. (34). Finally, the weights are found using Eq. (32).

As described in [8

R. Hardie, “A fast image super-resolution algorithm using an adaptive Wiener filter,” IEEE Trans. Image Process. 16(12), 2953–2964 (2007). [CrossRef] [PubMed]

], there are circumstances where one only needs to compute one set of weights for processing all the observation windows for an entire output frame. This allows for very fast processing. This fortunate case occurs when we have global translational motion and the statistical parameters are assumed to be constant across the image. With global translational motion, the same spatial pattern of LR samples is observed in each observation window, so long as the observation and estimation window sizes are an integer multiple of L. If we have bad pixels, we can simply exclude those from the observation vector g _i , so that they do not contribute to the output. However, this requires that a distinct set of weights be computed around any bad pixels, as it changes the spatial distribution of samples in the observation window. If processing speed is critical, a bad pixel replacement algorithm can be employed prior to AWF filtering [16

]. This faster approach is how we address bad pixels for the AWF SR results presented in Section 5.

An example of some polarimetric AWF filter weights is shown in Fig. 6. The brightest pixel corresponds to a weight of 1.41 and the darkest is −0.17. The remaining weights are scaled linearly. For reference, note that the perimeter weights are very close to zero. These weights are computed for K = 1 frame, P = 4 polarimetric channels, and an upsampling factor of L = 2. Thus, we assume we have a single microgrid frame and we form an estimate of each polarization at each FPA location. This can be viewed as demosaicing with restoration. This simple case is shown to make interpreting the weights more straightforward than with mulitple frames and/or more upsampling. In computing these weights, we assume the noise variance is

σ_{n}^{2} = 1

and

σ_{p}^{2} = 100

. The observation window size is W_x = W_y = 10 and the estimation window size is D_x = D_y = 2. The weights correspond to estimating the upper left pixel in the estimation window. This location is highlighted with a plus sign in Fig. 6. Each of the four columns in Fig. 6 corresponds to a different polarization sample being estimated for the one spatial location, and each row corresponds to using a different α in computing the weights. In particular, the columns from left to right correspond to polarization angles of θ = 0°, 45°, 90° and 135°. The rows from top to bottom correspond to α = 0, 0.7, 1.0. Note that when α = 0 (top row), the weights are only nonzero for pixels on the microgrid of the same polarization as the output. On the other extreme, when α = 1 (bottom row), the weights treat all pixels on the migrogrid array as equivalent, ignoring their polarization differences. Thus, the output for all polarizations is formed with the same weights, yielding the same output. This would be appropriate when the DoLP is assumed to be zero, and all polarization channels should be equal. The case where α = 0.7 (middle row) yields an interesting result, where samples from all polarizations are combined to form each output, especially when estimating at a location where that polarization is not measured (columns 2–4). The advantage of using samples with polarizations other than the one being estimated is that they offer spatial locations not represented by pixels of like polarization. Spatially close pixels of a different polarization may be more correlated to the sample we are estimating than ones farther away that have the same polarization. Thus, using α, we control the weight given to samples with polarizations different than the polarization being estimated relative to those of like polarization. When multiple frames are present, the picture is more intricate as samples of various polarizations are distributed potentially nonuniformly throughout the observation window. However, the principle is the same where α balances spatial correlation with polarization correlation.

Fig. 6 AWF filter weights for K = 1 frame, P = 4 channels, W_x = W_y = 10 observation window, D_x = D_y = 2 estimation window and L = 2 upsampling (single frame demosaicing with restoration). The weights correspond to estimating the pixel location with the red plus sign. The columns from left to right correspond to polarization angles of θ = 0°, 45°, 90° and 135°, respectively. The rows from top to bottom correspond to α = 0, 0.7, 1.0.

We have found that in some cases we can improve the microgrid AWF results by preprocessing the microgrid to compensate for disparate local statistics of the individual channels. In particular, if we subtract the local means and divide by the local standard deviations of the individual channels, they tend to exhibit greater correlation. This allows us to use a larger α and better exploit the spatial sampling diversity provided by the microgrid as shown in Fig. 1. To better describe this method, let a(i, j) represent a sample from the microgrid array, where i = 1,2,...,M is the 2 × 2 super-pixel index and j is the polarization sample index within the super-pixel. Let j = 1,2,3,4 correspond to θ = 0°, 45°, 90⁰, 135°, respectively. Local means, denoted ā(i, j), are estimated by deinterlacing the channels and applying a Gaussian low-pass filter with standard deviation of σ pixels. Local standard deviations, denoted ã(i, j), are estimated by subtracting the local mean and then using a Gaussian low-pass filter on the squared difference, and taking the square root of the result. In terms of these variables, the preprocessing step is expressed as

b (i, j) = \frac{a (i, j) - ā (i, j)}{ã (i, j)} .

(39)

After SR processing, the local means and standard deviations are interpolated, shifted appropriately to match the HR grid, and reintroduced to the SR result to restore the estimated channels to their proper dynamic range. Even without SR processing, we have found that a local statistics fusion (LSF) method can be useful as a stand-alone demosaicing approach. In that case, the output is given by

b (i, j, k) = {\begin{matrix} a (i, j) & j = k \\ ã_{j} (i, k) \frac{a (i, j) - ā (i, j)}{ã (i, j)} + ā_{j} (i, k) & otherwise \end{matrix} .

(40)

Here b(i, j,k) is the output for super-pixel i = 1,2,...,M, position j = 1,2,3,4 within the super-pixel, and with polarization index k = 1,2,3,4. The variables ā_j (i, k) and ã_j (i, k) are the mean and standard deviation, respectively, for super-pixel i and polarization k repositioned through interpolation to align with spatial position j. Note that the low-pass nature of the statistic images make them more suited to interpolation than the raw data. This method works best when there is high correlation between the polarimetric channels, but perhaps a scale and bias difference.

5. Experimental results

In this section, we present a number of experimental results to compare and contrast the performance of the AWF and RLS microgrid SR methods with each other and with other single frame methods. The first set of results are from the LWIR microgrid polarimetric imager described in Section 2. We also consider data from a visible polarimetric sensor, where we simulate microgrid sampling. These data allow for quantitative error analysis.

5.1. LWIR microgrid polarimetric data

The first set of results are shown in Fig. 7 and continue in Fig. 8. A single raw LWIR microgrid image is shown in Fig. 7(a). The sensor has a number of bad pixels including a stripe of bad pixels near column 120. The bad pixel map is shown as a red contour in Fig. 7(a). The remaining images in Figs. 7 and 8 are various estimates of p ₀ for L = 4. In particular, the result of using single-frame channel-independent bicubic interpolation is shown in Fig. 7(b). The output using the single-frame frequency-domain method by Tyo et al [3

], which is built upon the work in [2

B. M. Ratliff, C. F. LaCasse, and J. S. Tyo, “Interpolation strategies for reducing IFOV artifacts in microgrid polarimeter imagery,” Opt. Express 17(11), 9112–9125 (2009). [CrossRef] [PubMed]

], is shown in Fig. 7(c). The output for the single-frame LSF method from Eq. (40) with σ = 1 is shown in Fig. 7(d). The microgrid AWF SR method using only a single frame with ρ = α = 0.7 and σ = 1 produces the result shown in Fig. 7(e). The microgrid RLS SR method using only a single frame with

σ_{n}^{2} = 1

σ_{s_{0}}^{2} = 100

, and

σ_{s_{1}}^{2} = σ_{s_{2}}^{2} = σ_{r}^{2} = 10

produces the result shown in Fig. 7(f). The selection of the cost function parameters has been informed by quantitative analysis using simulated data as well as subjective evaluation. Note that the bicubic image is quite blurred and shows aliasing artifacts. These artifacts are most pronounced on the starboard wing and vertical stabilizer of the aircraft. Here the leading edges of these surfaces appear jagged. The Tyo and LSF methods do appear to reduce some of the more obvious aliasing artifacts. They do this by exploiting the sampling diversity of all polarizations on the microgrid array. The AWF and RLS outputs using a single frame appear to show further improvements as they also provide some deconvolution of the system PSF.

Fig. 7 LWIR microgrid imagery results. (a) Raw microgrid image with bad pixel mask shown in red. Estimates of p ₀ using an upsampling factor of L = 4 are shown for (b) bicubic interpolation of polarimetric channel g ₀(1) (c) Tyo et al method from [3

] (d) LSF method (e) single frame AWF with correlation model parameters of ρ = α = 0.7 (f) single frame RLS.

Fig. 8 LWIR microgrid imagery results continued. (a) Single frame AWF with α = 0.0 (b) single frame AWF with α = 1.0 (c) independent-channel single-frame RLS (d) independent-channel 20 frame RLS (e) 20 frame microgrid AWF SR with ρ = α = 0.7 (f) 20 frame microgrid RLS SR.

To see the impact of α on the microgrid AWF SR method, we show the single frame result using α = 0.0 and α = 1.0 in Figs. 8(a) and 8(b), respectively. Note that for α = 0.0, the polarimetric channels are treated independently [8

R. Hardie, “A fast image super-resolution algorithm using an adaptive Wiener filter,” IEEE Trans. Image Process. 16(12), 2953–2964 (2007). [CrossRef] [PubMed]

] and the result is similar to that of bicubic interpolation. Using α = 1.0, the AWF treats all polarimetric channels as being the same. While this method can sometimes produce a reasonable {s _0, _j } image, the {s _1, _j } and {s _2, _j } Stokes images are zero or have only low spatial frequency content from the local statistics processing in Eq. (39). Note that cross-hatching artifacts can be seen in Fig. 8(b), especially near the cockpit window at the front of the aircraft. These types of artifacts tend to become significantly more pronounced with α = 1.0 when more input frames are used. The single-frame RLS algorithm operating on the polarimetric channels independently [5

B. M. Ratliff, J. S. Tyo, W. T. Black, and C. F. LaCasse, “Exploiting motion-based redundancy to enhance microgrid polarimeter imagery,” Proc. SPIE 7461, 74610K (2009). [CrossRef]

] is shown in Fig. 8(c). Again, the output looks similar to bicubic interpolation as only samples of one polarization are used. Using the independent-channel RLS with 20 input frames gives rise to the output shown in Fig. 8(d). Here the results are much better. Given a large number of input frames, one can often get a good result using only samples of polarization matching that of the output. The 20 frame microgrid AWF output with ρ = α = 0.7 and σ = 9 is presented in Fig. 8(e). Finally, the 20 frame microgrid RLS output with

σ_{n}^{2} = 5

σ_{s_{0}}^{2} = 100

, and

σ_{s_{1}}^{2} = σ_{s_{2}}^{2} = σ_{r}^{2} = 10

are shown in Fig. 8(f). The noise variance term is increased here to add robustness to possible registration errors. These multiframe results appear to show the best detail on the roofs of the buildings in the foreground. Many diagonal lines that appear in these images look aliased in most of the other images.

To illustrate the typical convergence behavior of the microgrid RLS algorithm, the cost function and its components are shown in Fig. 9 for the single frame LWIR result in Fig. 7(e). Note that the total cost decreases smoothly with gradient-descent iteration. The observation model (or data) term also declines with iteration. The Stokes regularization terms actually increase with iteration. This is because the image estimate becomes sharper than the interpolated starting image as the iterations progress. These terms help keep the final estimate from becoming overly noisy. Note that the constants in the cost function can be modified to give more or less weight to any of these components in order to shape the characteristics of the final solution.

Fig. 9 Microgrid RLS SR algorithm cost function and its components for the single frame LWIR result in Fig. 7(e).

The DoLP images calculated from various estimates of p are shown in Fig. 10. Since DoLP involves all of the polarimetric channels, these results reveal characteristics of the estimates in all channels. DoLP also tends to expose sampling errors rather dramatically because it is made up of image differences, as can be seen in Eq. (2). The DoLP images for bicubic interpolation, the Tyo method, and the LSF method are shown in Figs. 10(a)–10(c), respectively. The DoLP image for the independent channel RLS method [5

B. M. Ratliff, J. S. Tyo, W. T. Black, and C. F. LaCasse, “Exploiting motion-based redundancy to enhance microgrid polarimeter imagery,” Proc. SPIE 7461, 74610K (2009). [CrossRef]

] using 20 frames is shown in Fig. 10(d). Finally, the DoLP images for 20 frame microgrid AWF SR and 20 frame microgrid RLS SR are shown in Figs. 10(e) and 10(f), respectively. Based on subjective analysis, we believe that the multiframe microgrid AWF and RLS methods appear to produce the sharpest results with the least aliasing.

Fig. 10 DoLP images for LWIR microgrid data. (a) Bicubic interpolation (b) Tyo et al (c) LSF method (d) independent-channel 20 frame RLS (e) 20 frame microgrid AWF SR with ρ = α = 0.7 (f) 20 frame microgrid RLS SR.

5.2. Visible polarimetric data with simulated microgrid sampling

To allow for quantitative analysis with truth imagery, we use a visible imaging system equipped with a linear polarization filter. The camera is an Imaging Source DMK 21BU04. This is a 640 × 480 8-bit grayscale camera with a Sony ICX098BL CCD sensor with 5.6μm detectors. The camera is fitted with a Computar 5mm F/1.4 lens. Full frames at each of the four polarization angles are acquired of an outdoor scene containing two vehicles. The polarization filter is rotated in 45° increments between acquisition to provide the four full-frame polarization images. We define these data to be the true p. These data are then put through the observation model with L = 4 and a noise variance of 1. We use a set of randomly generated shifts and the PSF in Fig. 4 to create 20 simulated microgrid images. The Stokes intensity image results for these data are shown in Fig. 11. In particular, one of the simulated microgrid images is shown in Fig. 11(a). The corresponding true {s _0, _j } Stokes image is shown in Fig. 11(b). Bicubic interpolation yields the result presented in Fig. 11(c). Note that the bicubic interpolation result looks blurred relative to the true image. Since {s _0, _j } is an average of the four polarimetric channels, some of the aliasing artifacts that are more clearly seen in the individual polarimetric channels are reduced here. The output of the Tyo method is shown in Fig. 11(d) and the 20 frame microgrid AWF output with ρ = α = 0.7 and σ = 9 is shown in Fig. 11(e). Finally, the 20 frame microgrid RLS output with

σ_{n}^{2} = 1

σ_{s_{0}}^{2} = 100

, and

σ_{s_{1}}^{2} = σ_{s_{2}}^{2} = σ_{r}^{2} = 10

is shown in Fig. 11(f).

Cross-sections of the true image {s _0, _j } from Fig. 11 as well as several of the corresponding estimate images along Row 234 and about Column 67 are shown in Fig. 12. These cross-sections are centered about a dark line in the sidewalk concrete, providing us with approximate line spread functions. These provide some insight into the resolution characteristics of the various image estimates. Here it appears that the microgrid RLS provides the narrowest line spread, followed closely by the microgrid AWF. The single frame methods clearly have much broader line spread functions. This is to be expected as these methods do not include any deblurring. It is interesting to note that the Tyo method does appear to have a more favorable line spread function than bicubic interpolation.

Fig. 11 Visible polarimetric data with simulated microgrid sampling. (a) Simulated microgrid frame (b) true {s _0, _j } image for an upsampling factor of L = 4 (c) estimate using bicubic interpolation (d) Tyo et al method (e) 20 frame microgrid AWF SR with ρ = α = 0.7 (f) 20 frame microgrid RLS SR.

Fig. 12 Cross-sections of the true image {s _0, _j } and corresponding estimates along Row 234 and about Column 67 from Fig. 11. These plots represent approximate line spread functions and provide some insight into the resolution of the various estimates.

The corresponding DoLP images for the key methods operating on the simulated microgrid data are provided in Fig. 13. The true DoLP image is shown in Fig. 13(a). The DoLP generated with bicubic interpolation is shown in Fig. 13(b). Here, the aliasing artifacts are exposed, especially along edges in the scene content. The DoLP images for Tyo and LSF are shown in Figs. 13(c) and 13(d), respectively. The DoLP image for the 20 frame microgrid AWF is shown in Fig. 13(e) and the DoLP image for the 20 frame microgrid RLS is shown in Fig. 13(f). Based on subjective analysis, we believe that the multiframe AWF and RLS methods appear to provide the best estimates of the true Stokes intensity image and DoLP image.

Fig. 13 DoLP images for simulated microgrid data. (a) True DoLP image (b) bicubic interpolation (c) Tyo et al method (d) LSF method (e) 20 frame microgrid AWF SR with ρ = α = 0.7 (f) 20 frame microgrid RLS SR.

The MSEs between the true Stokes image and the various estimates are shown in Fig. 14 as a function of the number of input frames. In particular, the MSEs for {s _0, _j }, {s _1, _j }, and {s _2, _j } are shown in Figs. 14(a)–14(c), respectively. Note that the Tyo and LSF methods do outperform bicubic interpolation, as expected. Furthermore, the single-frame microgrid AWF with α = 0.7 and the microgrid RLS methods do better yet. We believe this is in part due to the fact that these methods incorporate knowledge of the PSF and perform some deblurring. As more input frames are used, the microgrid AWF and RLS SR estimates improve, with diminishing returns as the sampling grid gets well populated by numerous shifted frames. Note also that the α = 0.7 result outperforms the independent channel AWF (α = 0.0). The multichannel RLS also outperforms its independent channel counterpart [5

B. M. Ratliff, J. S. Tyo, W. T. Black, and C. F. LaCasse, “Exploiting motion-based redundancy to enhance microgrid polarimeter imagery,” Proc. SPIE 7461, 74610K (2009). [CrossRef]

]. It is also interesting to note that the AWF with α = 1.0 provides a relatively low MSE on {s _0, _j }, but does a poor job with {s _1, _j } and {s _2, _j } that does not improve with input frames. This is because with α = 1.0, the {s _1, _j } and {s _2, _j } estimates are nearly zero, save for the impact of the local statistics processing beginning with Eq. (39). These results suggest that when estimating {s _0, _j } with a relatively small number of frames, it is beneficial to assume more correlation between channels (i.e., use a higher α) to more fully exploit the spatial sampling diversity offered by the different channels. However, this comes at the expense of the polarization content in {s _1, _j } and {s _2, _j }. We also see that when enough frames are available, the benefit of multichannel processing is somewhat reduced, due to the abundance of spatial samples of each polarization.

Fig. 14 MSE versus number of input frames for the Stokes images from the simulated microgrid imagery (a) {s _0, _j } (b) {s _1, _j } (c) {s _2, _j }.

On whole, the MSE results show that the RLS outperforms the AWF. We believe this is in part due to the finite window size of the AWF as well as the natural way in which the RLS incorporates the redundancy constraint in the cost function. However, the AWF offers a significant computational savings and yields a very competitive result with similar MSE. To illustrate the computational savings, consider the following run times using MATLAB on a PC with an Intel Core i7 CPU with 3.07GHz processor. For the 20 frame AWF, the run time was 0.31s per output frame. For the 20 frame RLS, the run time was 82.46s per output frame. For reference, the run times per output frame were 0.093s, 0.071s, and 0.133s for bicubic, Tyo, and LSF, respectively.

Finally, the MSE results in Fig. 15 show how the α parameter impacts the MSE performance of the microgrid AWF filter on the three Stokes images. These results are for the 20 frame microgrid AWF output with ρ = 0.7 and σ = 9. Note that the MSE on {s _0, _j } is minimized with an α close to, but distinctly less than, one. The MSE on the other Stokes images is obtained with an α near 0.5. Again, this shows us that if one assumes too much correlation between channels (i.e., high α), the pure polarimetric content in {s _1, _j } and {s _2, _j } is degraded. On the other hand, assuming too little correlation predominantly hurts our estimate of {s _0, _j }. We have found a good compromise in the value used in our earlier results of α = 0.7.

Fig. 15 MSE versus α for the 20 frame microgrid AWF output with ρ = 0.7.

6. Conclusions

Microgrid imaging systems are likely to suffer from undersampling to a greater extent than their single-channel counterparts. We have presented two new SR algorithms designed for the restoration of microgrid polarimetric data in the presence of aliasing, blur and noise. These are the microgrid RLS and microgrid AWF algorithms. We have also presented the simple LSF demosacing algorithm based on local statistics. Each in their own way, the new microgrid SR methods exploit the sampling diversity offered by the microgrid pattern and correlation between different polarization channels. Unlike simpler independent-channel methods, the new microgrid algorithms use samples from all polarizations to estimate each output at each polarization. This allows the new microgrid algorithms to outperform their independent channel counterparts. By using multiple frames, when available, we further exploit sampling diversity provided by relative motion between the scene and the FPA. Even with a single frame, the AWF and RLS provided lower MSE than the other single frame methods tested here. As the number of input frames increases, the algorithms perform even better. We did observe that independent-channel SR processing can do reasonably well when a large number of input frames are available. However, the microgrid approach using all the polarizations is critical to combatting aliasing with a relatively small number of input frames. Operating with a small number of input frames is desirable in many applications because we are less likely to encounter scene motion, motion parallax, or other changes that make registration more difficult. Finally, the microgrid AWF SR method has a low computational complexity for translational motion. It is much faster than the RLS in that case, but enjoys a very similar performance, as seen in the subjective and quantitative results presented here.

Acknowledgments

This work was sponsored under AFRL contract FA8650-10-7028. The authors would like to thank Dr. Kenneth Barnard at AFRL for assisting in getting this project underway. We would also like to thank Ernie Atkins at DRS Sensors and Tracking for providing the LWIR data. Finally, thanks to Charles LaCasse for providing the MATLAB implementation of the Tyo et al polarimetric demosaicing algorithm used here.

References and links

1.	J. S. Tyo, D. L. Goldstein, D. B. Chenault, and J. A. Shaw, “Review of passive imaging polarimetry for remote sensing applications,” Appl. Opt. 45(22), 5453–5469 (2006). [CrossRef] [PubMed]
2.	B. M. Ratliff, C. F. LaCasse, and J. S. Tyo, “Interpolation strategies for reducing IFOV artifacts in microgrid polarimeter imagery,” Opt. Express 17(11), 9112–9125 (2009). [CrossRef] [PubMed]
3.	J. S. Tyo, C. F. LaCasse, and B. M. Ratliff, “Total elimination of sampling errors in polarization imagery obtained with integrated microgrid polarimeters,” Opt. Lett. 34(20), 3187–3189 (2009). [CrossRef] [PubMed]
4.	S. C. Park, M. K. Park, and M. G. Kang, “Super-resolution image reconstruction: a technical overview,” IEEE Signal Process. Mag. 20(3), 21–36 (2003). [CrossRef]
5.	B. M. Ratliff, J. S. Tyo, W. T. Black, and C. F. LaCasse, “Exploiting motion-based redundancy to enhance microgrid polarimeter imagery,” Proc. SPIE 7461, 74610K (2009). [CrossRef]
6.	D. A. Lemaster, “Resolution enhancement by image fusion for microgrid polarization imagers,” in IEEE Aerospace Conference, Big Sky, MT (2010).
7.	R. C. Hardie, K. J. Barnard, J. G. Bognar, E. E. Armstrong, and E. A. Watson, “High resolution image reconstruction from a sequence of rotated and translated frames and its application to an infrared imaging system,” Opt. Eng. 37(1), 247–260 (1998). [CrossRef]
8.	R. Hardie, “A fast image super-resolution algorithm using an adaptive Wiener filter,” IEEE Trans. Image Process. 16(12), 2953–2964 (2007). [CrossRef] [PubMed]
9.	R. M. A. Azzam and N. M. Bashara, Ellipsometry and Polarized Light (North-Holland, 1977).
10.	J. Goodman, Introduction to Fourier Optics (McGraw-Hill, 1968).
11.	R. R. Shannon, “Aberrations and their effect on images,” Proc. SPIE 531, 27–37 (1985).
12.	T. Gotoh and M. Okutomi, “Direct super-resolution and registration using raw CFA images,” in IEEE Conference on Computer Vision and Pattern Recognition , vol. 2, pp. 600–607 (Los Alamitos, CA, 2004).
13.	S. Farsiu, M. Elad, and P. Milanfar, “Multi-frame demosaicing and super-resolution of color images,” IEEE Trans. Image Process. 15, 141–159 (2006). [CrossRef] [PubMed]
14.	L. Condat, “A generic variational approach for demosaicking from an arbitrary color filter array,” in Proceedings of IEEE ICIP , pp. 1625–1628 (2009).
15.	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 Trans. Image Process. 6(12), 1621–1633 (1997). [CrossRef] [PubMed]
16.	B. M. Ratliff, J. S. Tyo, J. K. Boger, W. T. Black, D. L. Bowers, and M. P. Fetrow, “Dead pixel replacement in LWIR microgrid polarimeters,” Opt. Express 15(12), 7596–7609 (2007). [CrossRef] [PubMed]

OCIS Codes
(100.6640) Image processing : Superresolution
(120.2130) Instrumentation, measurement, and metrology : Ellipsometry and polarimetry
(120.5410) Instrumentation, measurement, and metrology : Polarimetry
(110.5405) Imaging systems : Polarimetric imaging

ToC Category:
Imaging Systems

History
Original Manuscript: May 6, 2011
Revised Manuscript: June 1, 2011
Manuscript Accepted: June 2, 2011
Published: June 20, 2011

Citation
Russell C. Hardie, Daniel A. LeMaster, and Bradley M. Ratliff, "Super-resolution for imagery from integrated microgrid polarimeters," Opt. Express 19, 12937-12960 (2011)
http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-19-14-12937

Sort: Author | Year | Journal | Reset

References

J. S. Tyo, D. L. Goldstein, D. B. Chenault, and J. A. Shaw, “Review of passive imaging polarimetry for remote sensing applications,” Appl. Opt. 45(22), 5453–5469 (2006). [CrossRef] [PubMed]
B. M. Ratliff, C. F. LaCasse, and J. S. Tyo, “Interpolation strategies for reducing IFOV artifacts in microgrid polarimeter imagery,” Opt. Express 17(11), 9112–9125 (2009). [CrossRef] [PubMed]
J. S. Tyo, C. F. LaCasse, and B. M. Ratliff, “Total elimination of sampling errors in polarization imagery obtained with integrated microgrid polarimeters,” Opt. Lett. 34(20), 3187–3189 (2009). [CrossRef] [PubMed]
S. C. Park, M. K. Park, and M. G. Kang, “Super-resolution image reconstruction: a technical overview,” IEEE Signal Process. Mag. 20(3), 21–36 (2003). [CrossRef]
B. M. Ratliff, J. S. Tyo, W. T. Black, and C. F. LaCasse, “Exploiting motion-based redundancy to enhance microgrid polarimeter imagery,” Proc. SPIE 7461, 74610K (2009). [CrossRef]
D. A. Lemaster, “Resolution enhancement by image fusion for microgrid polarization imagers,” in IEEE Aerospace Conference, Big Sky, MT (2010).
R. C. Hardie, K. J. Barnard, J. G. Bognar, E. E. Armstrong, and E. A. Watson, “High resolution image reconstruction from a sequence of rotated and translated frames and its application to an infrared imaging system,” Opt. Eng. 37(1), 247–260 (1998). [CrossRef]
R. Hardie, “A fast image super-resolution algorithm using an adaptive Wiener filter,” IEEE Trans. Image Process. 16(12), 2953–2964 (2007). [CrossRef] [PubMed]
R. M. A. Azzam and N. M. Bashara, Ellipsometry and Polarized Light (North-Holland, 1977).
J. Goodman, Introduction to Fourier Optics (McGraw-Hill, 1968).
R. R. Shannon, “Aberrations and their effect on images,” Proc. SPIE 531, 27–37 (1985).
T. Gotoh and M. Okutomi, “Direct super-resolution and registration using raw CFA images,” in IEEE Conference on Computer Vision and Pattern Recognition , vol. 2, pp. 600–607 (Los Alamitos, CA, 2004).
S. Farsiu, M. Elad, and P. Milanfar, “Multi-frame demosaicing and super-resolution of color images,” IEEE Trans. Image Process. 15, 141–159 (2006). [CrossRef] [PubMed]
L. Condat, “A generic variational approach for demosaicking from an arbitrary color filter array,” in Proceedings of IEEE ICIP , pp. 1625–1628 (2009).
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 Trans. Image Process. 6(12), 1621–1633 (1997). [CrossRef] [PubMed]
B. M. Ratliff, J. S. Tyo, J. K. Boger, W. T. Black, D. L. Bowers, and M. P. Fetrow, “Dead pixel replacement in LWIR microgrid polarimeters,” Opt. Express 15(12), 7596–7609 (2007). [CrossRef] [PubMed]

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.

Click here to see a list of articles that cite this paper

Figures


Fig. 1	Fig. 2	Fig. 3


Fig. 4	Fig. 5	Fig. 6


Fig. 7	Fig. 8	Fig. 9


Fig. 10	Fig. 11	Fig. 12


Fig. 13	Fig. 14	Fig. 15