Integrations of both high resolution reconstruction and non-uniformity correction of infrared image sequence based on regularized maximum a posteriori

Abstract

During thermal imaging, it is vital importance to obtain high-performance images that non-uniformity noise in infrared focal plane array (IRFPA) should be eliminatined and the imaging spatial resolution should be improved as far as possible. Processing algorithms related to both of them have been hot topics, and attracted more and more attention of researchers. Considering that both high-resolution restoration algorithm of image sequences and scene-based non-uniformity correction (NUC) algorithm require multi-frame image sequences of target scene with micro-displacement, an integrated processing algorithm of high-resolution image reconstruction and NUC of infrared image sequences based on regularized maximum a posteriori (MAP) is proposed. Results of simulated and experimental thermal image suggested that this algorithm can suppress random noise and eliminate non-uniformity noise effectively, and high-resolution thermal imaging can be achieved.

Key words： infrared image; image sequences; motion estimation; non-uniformity correction (NUC); maximum a posteriori (MAP) restoration

Introduction

Infrared focal plane array (IRFPA) sensors are widely used in military, industry, medical imaging, environmental monitoring and other fields. There are many researches focusing on the improving resolution in thermal imaging to obtain high-quality images. Besides system random noise, there are many main factors to influence thermal imaging quality. One is fixed pattern noise (FPN) or spatial nonuniformity, which is key issue in IRFPA sensors; Others are the blur and distortion caused by optical aberration, atmospheric disturbance, movement and defocus during the imaging process, which can be described as point spread function (PSF); Moreover, because of the big pixel size and the small pixel number of IRFPA detector, not only a detector unit integral effect caused by spatial resolution is reduced, but also the discrete sampling imaging caused by spatial sampling rate is insufficient, produce under sampled imaging frequency aliasing effect.

At present, the non-uniformity correction (NUC) algorithms of IRFPA can be generally divided into two categories, namely reference source based two-point linear correction and scene-based adaptive correction. In recent years, the scene-based correction algorithm has been rapidly developed. More than ten NUC algorithms [1-4], such as artificial neural network, constant average statistics, Kalman filter algorithm, algebraic algorithm, etc, have been arisen, and these NUC algorithms generally require to estimate the inter-frame micro-displacement based on the image sequence with varying scene, for the purpose of effectively eliminating non-uniformity noise of IRFPA. On the other hand, the research of high-resolution reconstruction algorithm based on multiple undersampled images also has rapid development [5-8], of which the maximum a posteriori (MAP) based high-resolution reconstruction can regularize ill-posed problem [9,10], guarantee a unique solution and remove the noise effectively as well as. However, the MAP-based algorithm is been mostly used to aim at visible and infrared image restoration, and has not been related to the NUC. The idea of integrating NUC and multi-frame high-resolution restoration was firstly proposed by Armstrong et al. [11], who briefly described scene-based NUC techniques to pre-process the data for image registration and high image reconstruction. Later, a scene-based NUC and enhancement algorithm was proposed by Zhao et al. [12], the core of the proposed framework is a novel registration-based super-resolution method that is bootstrapped by statistical scene-based NUC method.

In this paper, the regularized MAP-based high-resolution image reconstruction is combined with NUC (RMAP-HR&NUC) to process the infrared image sequences polluted by non-uniformity noise. Thus, not only high-resolution image has been obtained, but also the non-uniformity has been corrected improving the image quality of the thermal imaging system.

Observation model

The image degradation process of thermal imaging system is shown in Fig. 1. BoldItalic denotes the discrete high-resolution infrared image with size M= M₁ × M₂, where M₁ = L₁ × N₁, and M₂ = L₂ × N₂, and the actually observed degradation images are denoted as {BoldItalic_k}, where k = 1,2, ..., K, and the size of BoldItalic_k is N= N₁ × N₂. L₁ and L₂ present resolution enhancement times of the row and column, respectively. Assuming that, the low-resolution image is obtained from the ideal high-resolution image BoldItalic by geometric space motion and distortion, linear optical system blurring, and the undersampling of infrared detector array, at the same time polluted by additive Gaussian noise. Because of the non-uniformity noise of IRFPA generally expressed by two parameters, namely gain and offset. Generally, after the pre-calibration correction in the system, the offset represents major part of dynamic FPN, in addition, because the low-frequency of image often varies slowly, the non-uniformity of each low-resolution image can be considered to be the same, and the degradation model of thermal imaging system containing non-uniformity noise can be expressed as

(1){\mathit{g}}_k={\mathit{D}}_k{\mathit{C}}_k{\mathit{F}}_k\mathit{f}+{\mathit{n}}_k+\mathit{b}={\mathit{H}}_k\mathit{f}+{\mathit{n}}_k+\mathit{b,}

where, the integrated degeneration matrix is BoldItalic_k = BoldItalic_kBoldItalic_kBoldItalic_k. BoldItalic_k denotes an M × M displacement matrix, which presents the relative motion between the kth frame image and the reference frame image, BoldItalic_k is the degraded matrix which represents optical blurring, with the size M × M, BoldItalic_k is the down-sampling matrix, with the size N × M, BoldItalic_k is the vector of N × 1, which presents the zero mean additive Gaussian noise in the kth frame image observed, and BoldItalic is the offset matrix, presenting the non-uniformity noise.

As the imaging system is generally linear space invariant, each low-resolution frame has the same degraded matrix BoldItalic and down-sampling matrix BoldItalic, therefore, the integrated degenerate matrix can be further expressed by BoldItalic_k = BoldItalic_k.

Fig.1 Imaging degradation model of infrared low-resolution image

Full size|PPT slide

Integrations of both high-resolution reconstruction and NUC based on regularized MAP

The regularized MAP high-resolution image reconstruction and NUC mean that, on the premise that the low-resolution infrared image sequences are known, the posterior probability of the high-resolution image and NUC offset matrix is maximized, namely,

(2)\langle \widehat{\mathit{f}},\widehat{\mathit{b}}\rangle =\arg \max \left[\ln \mathrm{Pr}\left(\mathit{f},\mathit{b}|{\mathit{g}}_{\mathrm{1}},{\mathit{g}}_{\mathrm{2}},...,{\mathit{g}}_K\right)\right].

According to Bayes conditional probability theory, taking into account the characteristic that the offset BoldItalic is unrelated with the high-resolution image BoldItalic, the MAP estimation of the high-resolution image BoldItalic and offset matrix BoldItalic are rewritten as

(3)\langle \widehat{\mathit{f}},\widehat{\mathit{b}}\rangle =\arg \max [\ln \mathrm{Pr}({\mathit{g}}_{\mathrm{1}},{\mathit{g}}_{\mathrm{2}},\mathrm{...},{\mathit{g}}_K|\mathit{f},\mathit{b})+\ln \mathrm{Pr}(\mathit{f})+\ln \mathrm{Pr}(\mathit{b})],

where, \ln \mathrm{Pr}({\mathit{g}}_{\mathrm{1}},{\mathit{g}}_{\mathrm{2}},\mathrm{...},{\mathit{g}}_K|\mathit{f},\mathit{b}) is the priori probability logarithm of maximum likelihood function, \ln \mathrm{Pr}(\mathit{f}) is the logarithm of a priori probability of high-resolution image BoldItalic, and \ln \mathrm{Pr}(\mathit{b}) is the logarithm of priori probability of offset parameter BoldItalic.

Assumed that, the degraded low-resolution image sequences are mutually independent, and image pixels are also statistically independent from each other. These image sequences have Gaussian noise with the same mean value of zero and a variance of {\sigma }_n^{\mathrm{2}}, then the probability-density function of the low-resolution image sequences is

(4)\mathrm{Pr}\left({\mathit{g}}_{\mathrm{1}},{\mathit{g}}_{\mathrm{2}},...,{\mathit{g}}_K|\mathit{f},\mathit{b}\right)={\displaystyle \prod _{k=\mathrm{1}}^K{\left({\displaystyle \frac{\mathrm{1}}{\sqrt{\mathrm{2}\mathrm{\pi }{\sigma }_n^{\mathrm{2}}}}}\right)}^N\exp \left(-{\displaystyle \frac{{\Vert {\stackrel{︵}{\mathit{g}}}_k-{\mathit{g}}_k\Vert }^{\mathrm{2}}}{\mathrm{2}{\sigma }_n^{\mathrm{2}}}}\right)}.

where \Vert ·\Vert denotes norm, {\stackrel{︵}{\mathit{g}}}_k={\mathit{H}}_k\mathit{f}+\mathit{b} refers to the simulated low-resolution images.

In order to solve ill-posed problems of the high-resolution reconstruction, the regularization method is used here to constrain the smoothness of the prior probability of high-resolution image BoldItalic, that is, to punish neighboring pixels with larger difference in the image, and the probability-density function is

(5)\mathrm{Pr}\left(\mathit{f}\right)=\exp \left(-{\displaystyle \frac{\mathrm{1}}{\rho }}{\Vert \mathit{Q}\mathit{f}\Vert }^{\mathrm{2}}\right),

where \rho is the control parameter, and BoldItalic is the linear high-pass filter operator, with two-dimensional Laplacian operator, etc. as an option, here we use the following operator:

(6)\mathit{Q}=\left[\begin{array}{ccccc}\mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{2}& -\mathrm{8}& \mathrm{2}& \mathrm{0}\\ \mathrm{1}& -\mathrm{8}& \mathrm{20}& -\mathrm{8}& \mathrm{1}\\ \mathrm{0}& \mathrm{2}& -\mathrm{8}& \mathrm{2}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}\end{array}\right].

Assuming that the non-uniformity offset parameter BoldItalic obeys the Gaussian distribution with a mean value of zero and a variance of {\sigma }^{\mathrm{2}}, then

(7)\mathrm{Pr}(\mathit{b})={\left({\displaystyle \frac{\mathrm{1}}{\sqrt{\mathrm{2}\mathrm{\pi }{\sigma }^{\mathrm{2}}}}}\right)}^N\exp \left\{-{\displaystyle \frac{{\Vert \mathit{b}\Vert }^{\mathrm{2}}}{\mathrm{2}{\sigma }^{\mathrm{2}}}}\right\}.

Taking Eqs. (4), (5), (7) into Eq. (3), the following equation is obtained

(8)\langle \widehat{\mathit{f}},\widehat{\mathit{b}}\rangle =\arg \max \left(-{\displaystyle \sum _{k=\mathrm{1}}^K{\displaystyle \frac{{\Vert {\stackrel{︵}{\mathit{g}}}_k-{\mathit{g}}_k\Vert }^{\mathrm{2}}}{\mathrm{2}{\sigma }_n^{\mathrm{2}}}}-{\displaystyle \frac{\mathrm{1}}{\rho }}{\Vert \mathit{Q}\mathit{f}\Vert }^{\mathrm{2}}-{\displaystyle \frac{{\Vert \mathit{b}\Vert }^{\mathrm{2}}}{\mathrm{2}{\sigma }^{\mathrm{2}}}}}\right).

Further assumed that {\sigma }_n^{\mathrm{2}} and {\sigma }^{\mathrm{2}} meet the following equation

(9){\sigma }^{\mathrm{2}}=\lambda {\sigma }_n^{\mathrm{2}},

where, λ denotes multiple factor.

Equation (8) can be interpretated to find the minimum of BoldItalic and BoldItalic with the sequence low-resolution infrared images BoldItalic_k are known beforehand:

(10)J(\mathit{f},\mathit{b})={\displaystyle \sum _{k=\mathrm{1}}^K{\displaystyle \frac{{\Vert {\mathit{H}}_k\mathit{f}+\mathit{b}-{\mathit{g}}_k\Vert }^{\mathrm{2}}}{\mathrm{2}}}+\alpha {\Vert \mathit{Q}\mathit{f}\Vert }^{\mathrm{2}}+{\displaystyle \frac{{\Vert \mathit{b}\Vert }^{\mathrm{2}}}{\mathrm{2}\lambda }}}\to \min ,

where, \alpha ={\sigma }_n^{\mathrm{2}}/\rho is the regularization parameter. Calculate the gradient of J(BoldItalic) to BoldItalic and BoldItalic:

(11){\nabla }_{f}J\left(\mathit{f},\mathit{b}\right)={\displaystyle \sum _{k=\mathrm{1}}^K{\mathit{H}}_k^{\mathrm{T}}\left({\mathit{H}}_k\mathit{f}+\mathit{b}-{\mathit{g}}_k\right)}+\mathrm{2}\alpha {\mathit{Q}}^{\mathrm{T}}\mathit{Q}\mathit{f},

(12){\nabla }_{b}J\left(\mathit{f},\mathit{b}\right)={\displaystyle \sum _{k=\mathrm{1}}^K\left({\mathit{H}}_k\mathit{f}+\mathit{b}-{\mathit{g}}_k\right)}+{\displaystyle \frac{\mathrm{1}}{\lambda }}\mathit{b}.

Through adopting the gradient based steepest descent iteration, the iteration equations for obtaining BoldItalic and BoldItalic are respectively shown as follows

(13){\mathit{f}}^{n+\mathrm{1}}={\mathit{f}}^n-\beta \left[{\displaystyle \sum _{k=\mathrm{1}}^k{\mathit{h}}_k^{\mathrm{t}}\left({\mathit{h}}_k{\mathit{f}}^n+{\mathit{b}}^n-{\mathit{g}}_k\right)}+\mathrm{2}\alpha {\mathit{q}}^{\mathrm{t}}\mathit{q}{\mathit{f}}^n\right],

(14){\mathit{b}}^{n+\mathrm{1}}={\mathit{b}}^n-\beta \left[{\displaystyle \sum _{k=\mathrm{1}}^k\left({\mathit{h}}_k{\mathit{f}}^n+{\mathit{b}}^n-{\mathit{g}}_k\right)}+{\displaystyle \frac{\mathrm{1}}{\lambda }}{\mathit{b}}^n\right],

where \beta is the iterative step for controlling the convergence properties and speed of algorithm.

In summary, the process of regularized MAP based high-resolution image reconstruction and NUC integration algorithm is as follows:

Step 1: Estimate global motion between the reference frame and the other frame of low-resolution image sequences, and map these displacements to the displacements of high resolution images according to restoration multiples (L₁ and L₂), thus the displacement matrix BoldItalic₁, BoldItalic₂, ..., BoldItalic_K is obtained. In this paper, employ gradient-based method [11,13] to obtain the global motion.

Step 2: After the cubic spline interpolation on the 1st frame BoldItalic₁ of the original low-resolution image, reduce the influence of the original non-uniformity using a 5×5 flat filter, and get {\mathit{f}}^{\mathrm{0}} as the initial estimate of high-resolution image; and set the initial offset matrix BoldItalic is 0.

Step 3: Enter into the large circulation, select \beta =\mathrm{1}, and use Eqs. (13) and (14) to implement iterative updates, calculate the gradient of the kth frame image, with the iterative process of Eqs. (11) and (12) as follows:

1)　Enter into the small circulation, introduce {\mathit{f}}^n into {\stackrel{︵}{\mathit{g}}}_k={\mathit{H}}_k\mathit{f}+\mathit{b} to get simulated low-resolution image {\widehat{\mathit{g}}}_k, where BoldItalic_k=BoldItalic_k. In the specific process, {\widehat{\mathit{f}}}^n is obtained from {\mathit{f}}^n which is compensated by BoldItalic_k, while {\widehat{\mathit{g}}}_k is obtained from {\widehat{\mathit{f}}}^n which is operated by matrix BoldItalic and BoldItalic. Where matrix BoldItalic is a Gaussian filter with size 5 × 5 to simulate the PSF of the system; and BoldItalic presents the down-sampling matrix.

2)　While k = K, complete the circulation, enter into Step 3), otherwise return to Step 1).

3)　Calculate the value of the constraint item, add to the above results, and get {\nabla }_{f}J(\mathit{f},\mathit{b}) and {\nabla }_{b}J(\mathit{f},\mathit{b}).

Step 4: When the difference between the high-resolution images obtained from two adjacent iterations is less than the preset threshold or the iteration time reaches the maximum number of iterations, terminate and exit from iterations, and output the high-resolution image BoldItalic. Otherwise, return to Step 3 and enter into the next round of circulation. The termination conditions are shown as follows

(15){\displaystyle \frac{\Vert {\mathit{f}}^{n+\mathrm{1}}-{\mathit{f}}^n\Vert }{\Vert {\mathit{f}}^n\Vert }}\le \delta ,

where \delta presents the threshold for terminating iterations.

Experimental results and discussion

The algorithm is used to process simulated and experimental infrared image sequences to verify its performance.

Restoration and correction processing experiments on simulated images

Fig.2 Experimental results on simulated video. (a) True high-resolution image, 240×320 pixel; (b) simulated frame-one clear low-resolution image, 120×160 pixel; (c) observed frame-one low-resolution image with μ = 0 and σ = 25, 120×160 pixel; (d) restored frame using the RMAP-HR&NUC algorithm (Q = 0.9737), 240×320 pixel; (e) low-resolution image of (b) using bilinear interpolation (Q = 0.9266), 240×320 pixel; (f) low-resolution image of (c) using bilinear interpolation (Q = 0.8809), 240×320 pixel

Full size|PPT slide

Figure 2 shows an example of high-resolution image reconstruction and NUC processing on simulated image sequence. First of all, a location in a large image (gray scale: 8-bit) is randomly selected, Given random displacements to simulate scene motion both in horizontal and vertical directions, and 15-frame images with size 240 × 320 are obtained as the original high-resolution image sequences using linear interpolation algorithm. Figure 2(a) illustrates the 1st high-resolution image frame. The down-sampling process is simulated by averaging the values of the four adjacent pixels per frame, 15-frame clear low-resolution image sequences with size 120×160 are achieved, and the 1st frame is shown as Fig. 2(b); and then, a pattern is randomly produced, which is conformed to Gaussian distribution (μ,σ). Adding it to the low-resolution image sequences as FPN (for the purpose of testing the adaptability of the algorithm under different non-uniformity pollution, three kinds of Gaussian noises are introduced to simulate non-uniformity for testing algorithm, namely ① μ = 0, σ = 5; ② μ = 0, σ = 15; ③ μ = 0, σ = 25, respectively), and introducing Gaussian noise (μ_n = 0, σ_n = 2) to image sequences. Figure 2(c) shows the 1st frame of the degraded image sequences while μ = 0, σ = 25. Figures 2(e) and 2(f) show the result of linear interpolation reconstruction of Figs. 2(b) and 2(c), respectively. The non-uniformity noise are around 25 grey levels while the total range is from 0 to 255 gray level, so set \alpha =\mathrm{0.01} and \lambda =\mathrm{100}, the proposed method is used to process the low-resolution degraded image sequence. Figure 2(d) shows the high-resolution image obtained through the reconstruction of the 15-frame low-resolution degraded image sequence. The figure also shows the restoration evaluation index value of image quality factor Q [14]. It is obvious that our algorithm removes non-uniformity effectively, and enables the reconstructed image closer to the original high-resolution image at the same time.

Further research is performed using the root mean square error (RMSE) [15] to evaluate the restored image BoldItalic and the offset parameter BoldItalic. Figure 3 shows the relative iteration numbers curves of the RMSE of BoldItalic and BoldItalic under different non-uniformity intensities. Through processing the degraded image sequence polluted by FPN, it can be found that: after about 30 iterations, the reconstructed high-resolution image BoldItalic is basically converged and stabilized; the offset parameter BoldItalic stops shocking after about 15 iterations, and gradually approaches the true offset value; while, the reconstructed image quality will decrease with the non-uniformity noise increases.

Fig.3 RMSE curves under different non-conformity condition. (a) RMSE curve of reconstructed image BoldItalic; (b) RMSE curve of offset parameter BoldItalic

Full size|PPT slide

Restoration and correction processing experiments on actual image

Figure 4 shows the test of high-resolution image reconstruction and NUC processing on the real infrared image sequence, where Fig. 4(a) shows the 1st frame of the actual 4-frame low-resolution infrared images with resolution 200 × 200, which shows that not only the resolution is low, but also the non-uniformity is serious. Through magnifying and reconstructing the 4-frame original infrared images, high-resolution image with resolution 400 × 400 is reconstructed (50 iterations). Figures 4(b) and 4(c) show the image reconstructed after the bilinear interpolation of the original image and reconstructed based on RMAP-HR&NUC algorithm respectively. It is obviously that, the image processed by the proposed algorithm is better than that processed by bilinear interpolation, which shows that our algorithm has a strong NUC and high-resolution image reconstruction capability.

Fig.4 Experiment on actual infrared image sequence. (a) Observed frame-one low-resolution infrared image, 200×200 pixel; (b) restored frame-one using bilinear interpolation algorithm, 400×400 pixel; (c) restored frame using RMAP-HR &NUC algorithm for 50 iterations, 400×400 pixel

Full size|PPT slide

Conclusions

This paper proposes the degradation model of thermal imaging system with non-uniformity noise. Considering that the inter-frame displacements need to be estimated both in super-resolution reconstruction and NUC of low-resolution infrared image sequences, a high-resolution reconstruction and NUC integrated mapping algorithm (RMAP-HR&NUC) is proposed based on regularized MAP, which can not only suppress random noise effectively, but also eliminate non-uniformity noise achieveing obvious image reconstruction results.

RMAP-HR&NUC algorithm approaches the optimal solution of non-uniformity offset correction matrix and restored image using the circulation iterative. Because of a larger processing load, it is difficult to perform hardware implementation. However, the post-processing for low-resolution infrared image sequence with serious non-uniformity (such as remote sensing image, reconnaissance image, etc.), can further improve the image quality for thermal imaging system, good adaptability and good result can be achieved. The algorithm is expected to become a more effective fast-processing technology in the furture.

References

1 Liu Z G, Hu X M, Lu J. An improved neural network non-uniformity correction for IRFPA. In: Proceedings of the Society for Photo-Instrumentation Engineers, 2009, 7383: 788330

2 Sui J, Jin W Q, Dong L Q. An adaptive nonuniformity correction algorithm for infrared line scanner based on local statistics. Chinese Optics Letters, 2007, 5(2): 74–76

3 Torres S N, Hayat M M. Kalman filtering for adaptive non-uniformity correction in infrared focal-plane arrays. J Opt Soc Am A Opt Image Sci Vis., 2003, 20(3):470–480

4 Ratliff B M, Hayat M M. An algebraic algorithm for non-uniformity correction in focal-plane arrays. . J Opt Soc Am A Opt Image Sci Vis. 2002, 19(9):1737–1747

5 Tsai R, Huang T. Multiframe image restoration and registration. Advances in Computer Vision and Image Processing, 1984, 1: 317–339

6 Su B H, Jin W Q. Super-resolution image resolution algorithm based on Poisson-Markov model. ACTA Electronica Sinica, 2003, 31(1): 41–44 (in Chinese)

7 Kim J Y, Park R H, Yang S. Super-resolution using POCS-based reconstruction with artifact reduction constraints. In: Proceedings of the Society for Photo-Instrumentation Engineers, 2005, 5960: 59605B

8 Sun G, Li Q H, Lu L. MAP algorithm to super-resolution of infrared images. In: Proceedings of the Society for Photo-Instrumentation Engineers, 2007, 6787: 67870K

9 Jonsson R. Regularization based super resolution imaging using FFT:s. In: Proceedings of the Society for Photo-Instrumentation Engineers, 2010, 5808: 122–131

10 Shen H F, Li P X, Zhang L P. Adaptive regularized MAP super-resolution reconstruction method. Geomatics and information science of Wuhan University, 2006, 31(11): 949–952 (in Chinese)

11 Armstrong E, Hayat M, Hardie R, Majeed M H. Non-uniformity correction for improved registration and high-resolution image reconstruction in IR imagery. In: Proceedings of the Society for Photo-Instrumentation Engineers, 1999, 3808: 150–161

12 Zhao W Y, Zhang C. Scene-based nonuniformity correction and enhancement: pixel statistics and subpixel motion. Journal of the Optical Society of America, 2008, 25(7): 1668–1681

PMID

13 Irani M, Peleg S. Improving resolution by image registration. CVGIP: graphical. Models and Image Process, 1991, 53(3): 231–239

14 Wang Z, Bovik A C. A universal image quality index. IEEE Signal Processing Letters, 2002, 9(3): 81–84

15 Torres S N, Vera E M, Rodrigo A R. Adaptive scene-based non-uniformity correction method for infrared-focal plane arrays. In: Proceedings of the Society for Photo-Instrumentation Engineers. Infrared Imaging Systems: Design, Analysis, Modeling, and Testing XIV. 2003, 5076: 130–139