## 1.

## Introduction

The spatial noise caused by nonuniformity of individual detector elements in the infrared focal plane array (IRFPA) can limit the overall performance of imaging systems. In general, the responsivity of individual detector elements is assumed to be linear. Therefore, it is possible to correct the nonuniformity using the well-known two-point correction (TPC) method, where two blackbodies at different temperatures are employed.^{1} Since the nonuniformity normally drifts in time^{2} and the correction capability of the TPC is degraded with repeated operations,^{3} we need to recalibrate the IRFPA. In practice, there are two difficulties in using the TPC method to compensate the residual nonuniformity (RNU) resulting from the temporal drifit or repeated operations: (1) we have to maintain two distinct heat sources in imaging systems; and (2) real-time video operation is interrupted during the correction process.^{4}

Various scene-based nonuniformity correction (SBNUC) algorithms have been proposed to solve these problems. In general, SBNUC schemes can be broadly divided into two categories: constant statistics (CS) methods^{4}5.^{–}^{6} and least mean square (LMS) methods.^{7}8.9.^{–}^{10} The original CS method assumes that the temporal mean and standard deviation of each pixel are constant over time and space.^{5} The performance of the original CS method is reliable as long as the assumption is valid. As pointed out in Refs. 4 and 6, however, thousands of image frames are required to hold the assumption. Zhang and Zhao^{4} proposed a local constant statistics (LCS) method, which assumes that the temporal statistics are constant in a local region around each pixel. The LCS method improved the correction performance for the same number of input frames.^{4} Later, Zuo et al.^{6} generalized the LCS method by introducing a new constraint called multiscale CS.

In the LMS methods, the correction parameters are learned by minimizing the LMS error between corrected images and desired output images. The minimization is performed frame-by-frame using the stochastic gradient descent (SGD) technique.^{11} Since the desired output images are not available, spatially low-pass filtered input images are used as the desired ones.^{7} However, the performance of the LMS method is degraded at strong edge points as reported in Refs. 8 to 10. Thus, several methods improve the estimation accuracy of correction parameters by suppressing the influence of strong edge points in the minimization process. Vera and Torres^{8} adaptively adjust the learning rate which is a fixed parameter in the original SGD techinque, according to the local spatial standard deviations of input images. In this way, the influence of strong edge regions, where the local standard deviation is large, is reduced in the minimization process. This method is further improved in Ref. 9 to remove burn-in artifacts caused by temporally slowly-varying image regions. An approach proposed in Ref. 10 updates correction parameters only when sufficient change occurs between consecutive images.

Lately, different error functions are proposed as alternatives to the LMS error in the original LMS method.^{12}13.^{–}^{14} Interframe registration-based LMS method^{12}^{,}^{13} first registers the previous corrected image with the current input image by assuming that only slight translational motion exists. Then, this technique minimizes the LMS error between the previously shifted image and the currently corrected image. Vera et al.^{14} minimize the total variation of corrected images to obtain the correction parameters.

Although these previous SBNUC algorithms may reduce the RNU, they have a common problem: a relatively large number of image frames are required to acquire the correction parameters.^{15} Two recently proposed SBNUC algorithms perform NUC using several image frames^{16} or even only two image frames.^{15} These methods can achieve good performance when the relative motion in successive image frames is a small translation along vertical or horizontal directions. However, large displacements between consecutive image frames can occur in some applications.^{17}^{,}^{18} To deal with the motion constraint of the previous approaches,^{15}^{,}^{16} we propose a new SBNUC method that estimates the correction parameters using several image frames. In the proposed method, we utilize the prior information on the parameters regarding the responsivity and the true scene irradiance. There is no restriction on the motion in successive image frames unless it is static.

The rest of this article is organized as follows. In Sec. 2, the proposed SBNUC method is detailed. In Sec. 3, the performance of the proposed method is evaluated. Conclusion is drawn in Sec. 4.

## 2.

## Proposed Method

In this section, we first formulate an optimization problem for correcting the RNU followed by its numerical solution.

## 2.1.

### Formulation

Let us assume that the characteristic of each detector element in the IRFPA is linear.^{5}^{,}^{7} Then, the acquired signal $y(i,j,t)$ for the $(i,j)$’th detector element at time $t$ is given by

^{15}

^{,}

^{19}

^{,}

^{20}Given the image observation model [Eq. (1)], our objective is to estimate $x(i,j,t)$ and $o(i,j,t)$. This problem can be solved by minimizing the proposed energy function, which consists of three terms:

## (2)

$$E(x,o)=\phantom{\rule{0ex}{0ex}}\sum _{t\in \mathrm{\Omega}}\sum _{i,j\in I}D(x,o)+{\lambda}_{o}{C}_{1}(o)+{\lambda}_{x}{C}_{2}(x),$$Solving Eq. (3) alone is an underconstrained problem in which the number of unknowns is greater than that of equations. Thus, a regularization approach is taken to estimate $x(i,j,t)$ and $o(i,j,t)$ in this work. The regularization term for the offset $o(i,j,t)$ is defined as follows:

${C}_{1}(o)$ is derived from the observation that the offset changes very slowly in time.^{7}^{,}^{9} In other words, the offset remains almost constant for several consecutive image frames.^{20} This regularization term favors the offset with small changes along the time axis. If we correctly estimate $o(i,j,t)$ for the given image frames, their temporal variation is negligible, which means that ${C}_{1}(o)$ is very close to zero.

The last term ${C}_{2}(x)$ in Eq. (2) is introduced to regularize the scene irradiance $x(i,j,t)$. In general, $x(i,j,t)$ is smooth in the spatial domain. This fact is implictly used in the original LMS method, where the desired image is the spatially low-pass filtered input image.^{7} Thus, it is natural for us to require $x(i,j,t)$ to be smooth in the spatial domain. Since the degree of the spatial smoothness can be measured via the image gradient, ${C}_{2}(x)$ is given by

## (5)

$${C}_{2}(x)={w}_{i}[y(i,j,t)]{\left[\frac{\partial x(i,j,t)}{\partial i}\right]}^{2}\phantom{\rule{0ex}{0ex}}+{w}_{j}[y(i,j,t)]{\left[\frac{\partial x(i,j,t)}{\partial j}\right]}^{2}.$$The smoothness term ${C}_{2}(x)$ is proportional to the magnitude of the spatial intensity change of the scene irradiance. Therefore, the smoother the scene irradiance is, the smaller the value of ${C}_{2}(x)$ is. However, the smoothness constraint is not appropriate at edge points as pointed out in Sec. 1. Since the spatial variation of $x(i,j,t)$ is normally greater than that of RNU,^{2} the large spatial variation in the input image $y(i,j,t)$ is mainly due to the edge points of $x(i,j,t)$. We adaptively adjust the effect of the smoothness constraint according to the gradient of $y(i,j,t)$. The weighting factors ${w}_{i}[y(i,j,t)]$ and ${w}_{j}[y(i,j,t)]$ in Eq. (5) are defined as follows

The exponent $\gamma $ controls the sensitivity to the spatial gradients of $y(i,j,t)$ and $\epsilon $ is a small constant that prevents dividing by zero. Since the weighting factors are inversely proportional to the spatial gradients of $y(i,j,t)$, the smoothness constraint has little effect on edge regions. These weighting factors are the same as smoothness weights for the image smoothing operator in Ref. 21.

## 2.2.

### Numerical Solution

The proposed energy function $E(x,o)$ in Eq. (2), a function of two variables, is nonconvex. We minimize the energy function by solving two convex subproblems in an alternating way with initial estimates ${x}^{(0)}=y$ and ${o}^{(0)}=0$:

The above process is repeated until there is no significant change in the estimates ${x}^{(n)}$ and ${o}^{(n)}$.

We compute the partial derivative of Eq. (2) with respect to $x$ in order to solve Eq. (8). First, we represent Eq. (2) in matrix notation as follows:

## (10)

$$\sum _{t\in \mathrm{\Omega}}{[\mathbf{y}(t)-\mathbf{x}(t)-\mathbf{o}(t)]}^{\mathbf{T}}[\mathbf{y}(t)-\mathbf{x}(t)-\mathbf{o}(t)]\phantom{\rule{0ex}{0ex}}+{\lambda}_{x}[\mathbf{x}{(t)}^{\mathbf{T}}{\mathbf{D}}_{\mathbf{i}}^{\mathbf{T}}{\mathbf{W}}_{\mathbf{i}}(t){\mathbf{D}}_{\mathbf{i}}\mathbf{x}(t)+\mathbf{x}{(t)}^{\mathbf{T}}{\mathbf{D}}_{\mathbf{j}}^{\mathbf{T}}{\mathbf{W}}_{\mathbf{j}}(t){\mathbf{D}}_{\mathbf{j}}\mathbf{x}(t)],$$## (11)

$$\frac{\partial E}{\partial \mathbf{x}}=\sum _{t\in \mathrm{\Omega}}[\mathbf{A}(t)\mathbf{x}(t)-\mathbf{b}(t)],$$## (12)

$$\mathbf{A}(t)={\lambda}_{x}[{\mathbf{D}}_{\mathbf{i}}^{\mathbf{T}}{\mathbf{W}}_{\mathbf{i}}(t){\mathbf{D}}_{\mathbf{i}}+{\mathbf{D}}_{\mathbf{j}}^{\mathbf{T}}{\mathbf{W}}_{\mathbf{j}}(t){\mathbf{D}}_{\mathbf{j}}],$$Therefore, we solve a large system of linear equations [i.e., $\mathbf{A}(t)\mathbf{x}(t)=\mathbf{b}(t)$] for each $\mathbf{x}(t)$. The conjugate-gradient (CG) method is used to solve the linear equations in this work since the matrix $\mathbf{A}$ is sparse, symmetric, and positive definite.^{22}

We rewrite Eq. (2) in different matrix notation from that represented in Eq. (10) to solve Eq. (9):

## (14)

$$\sum _{i,j\in I}{[\mathbf{y}(i,j)-\mathbf{x}(i,j)-\mathbf{o}(i,j)]}^{\mathbf{T}}[\mathbf{y}(i,j)-\mathbf{x}(i,j)-\mathbf{o}(i,j)]\phantom{\rule{0ex}{0ex}}+{\lambda}_{o}[\mathbf{o}{(i,j)}^{\mathbf{T}}{\mathbf{D}}_{\mathbf{t}}^{\mathbf{T}}{\mathbf{D}}_{\mathbf{t}}\mathbf{o}(i,j)],$$## (15)

$$\frac{\partial E}{\partial \mathbf{o}}=\sum _{i,j\in I}[{\mathbf{A}}^{\prime}\mathbf{o}(i,j)-{\mathbf{b}}^{\prime}(i,j)],$$## (16)

$${\mathbf{A}}^{\prime}=\mathbf{I}+{\lambda}_{o}{\mathbf{D}}_{\mathbf{t}}^{\mathbf{T}}{\mathbf{D}}_{\mathbf{t}},$$$I$ denotes identity matrix in Eq. (16). The CG method is used here again to obtain the values of the offset for each $\mathbf{o}(i,j)$.

## 3.

## Simulation Results

To our best knowledge, no study has been reported on correcting the RNU with several image frames that have large displacements. Therefore, no comparison is made with existing SBNUC methods in this work. The regularization parameters are empirically set to ${\lambda}_{o}=6$, ${\lambda}_{x}=0.1$ for all experiments in this article, and the value of $\gamma $ is determined to be 0.2. At first, we test the convergence of the proposed method with eight synthetic images as shown in Fig. 1. The eight images are generated by adding the artificial offset^{9}^{,}^{15} to calibrated infrared images caputured by a $320\times 256$ InSb focal plane array camera operating in the 3 to 5 *μ*m range. The RNU is generally composed of two patterns, the low-frequency one and white noise-like one as reported in Ref. 23. However, only the white noise-like pattern is usually prominent to observers together with natural scenes. This is due to the masking effect of the human visual system, which attenuates contrast sensitivity at low-spatial frequencies.^{24} Therefore, the artificial offset is generated as realizations of independent identically distributed Gaussian random variables.^{9}

We plot the proposed energy function in Eq. (2) against the number of iterations. As shown in Fig. 2, the value of the energy function drops quickly. We obtain good results with 11 iterations in our expeirments. In Fig. 3, we present images corrected by the proposed NUC method. Close-up views of some parts of the images are depicted in Fig. 3(b) to help the reader observe the visual quality improvement. The proposed method produces acceptable results no matter how complex the spatial distribution is, as shown in Fig. 3.

We investigate the effect of the number of input images on the estimated scene irradiance using the eight synthetic images. Table 1 shows the peak signal-to-noise ratio results for different numbers of input images. As the number of input images increases, we can obtain more accurate results. This can be described by the fact that the information gained from consecutive image frames leads to high-quality NUC results and enhanced temporal consistency in the offset. Note that, however, raising the number of input images increases processing time. Thus, selecting a proper number of input images demands trade-off between the computational complexity and the image quality.

## Table 1

Peak signal-to-noise ratio results of the proposed method with various numbers of input images.

Image # | Input | 2 images | 4 images | 6 images | 8 images |
---|---|---|---|---|---|

1 | 27.65 | 29.61 | 31.4 | 32.19 | 32.6 |

2 | 27.64 | 29.52 | 31.37 | 32.17 | 32.6 |

3 | 27.62 | — | 31.34 | 32.19 | 32.65 |

4 | 27.59 | — | 31.27 | 32.17 | 32.64 |

5 | 27.59 | — | — | 32.12 | 32.62 |

6 | 27.57 | — | — | 32.01 | 32.55 |

7 | 27.54 | — | — | — | 32.48 |

8 | 27.52 | — | — | — | 32.38 |

Average | 27.59 | 29.57 | 31.35 | 32.14 | 32.57 |

We also perform an experiment on two sets of real infrared images as shown in Fig. 4. We have collected the two sets of images using a $640\times 512$ InSb focal plane array camera operating in the 3 to 5 *μ*m range. One set of images in Fig. 4(a) shows drastic intensity changes among them due to atmospheric effects. The other set of images in Fig. 4(b) has relatively large motion. Objective results for the proposed method are provided in Table 2. We employ a roughness metric^{8}^{,}^{9}^{,}^{15} which is defined by

^{9}The correction results of the proposed method are depicted in Figs. 5 and 6. Similar to the simulated nonuniformity case, our method consistently suppresses the RNU as shown in Figs. 5(b) and 6(b). We also present the difference images between the input and the corrected in Fig. 7 to visualize the RNU corrected by the proposed method.

## Table 2

Roughness results for real images.

Image # | Fig. 4(a) | Fig. 4(b) | ||
---|---|---|---|---|

Input | Proposed | Input | Proposed | |

1 | 0.2211 | 0.0717 | 0.2183 | 0.0626 |

2 | 0.1290 | 0.0798 | 0.3981 | 0.0618 |

3 | 0.4618 | 0.0648 | 0.4338 | 0.0465 |

4 | 0.7389 | 0.0484 | — | — |

5 | 0.3975 | 0.0685 | — | — |

We have implemented the proposed method using C language. The simulation is performed on a PC with an Intel i7 3.40-GHz CPU and 4-GB memory. Our optimization procedure takes 5.6 and 3.8 s for the set of images in Figs. 4(a) and 4(b), respectively.

## 4.

## Conclusion

We presented a regularization approach to SBNUC with several image frames. Our work formulated the SBNUC process as the energy minization problem that incorporates the slowly varying nature of the detector responsivity and the smoothness constraint for the scene irradiance. In the proposed method, no assumption was made about the motion among input images except that the motion is static. Therefore, the proposed method can be used in applications where only several image frames are available and large displacements exist among the given images. Simulation results on both synthetic and real infrared images demonstrated that the proposed method can reduce the RNU. In future works, we plan to apply more advanced numerical techniques to reduce the computational complexity of the proposed method.

## References

## Biography

**Jun-Hyung Kim** received his BS and PhD degrees in electronic engineering from Korea University in 2006 and 2012, respectively. He has worked for the Agency for Defense Development since 2012. His current research interests are in the area of image processing and infrared imaging system. He is a member of SPIE.

**Jieun Kim** received her BS degree in electrical engineering from Busan National University in 2002, and her MS degree in electrical engineering from KAIST in 2004. She has worked for the Agency for Defense Development since 2004. Her current interests include digital image processing, target detection, and tracking.

**Sohyun Kim** is currently a research member at the Agency for Defense Development in Korea, and has over 10 years of experience in developing electro-optic systems. Her experience includes target detection algorithm design and real-time implementation for video tracker for infrared images. She holds a BS in physics from Sogang University and an MS in information and communications from Gwangju Institute of Science and Technology.

**Joohyoung Lee** received his BS and MS degress in electronics engineering from Dankuk University, Republic of Korea, in 1990 and 1992, respectively. Since 1992, he has been a principal researcher in the Electro-Optics Laboratory at the Agency for Defense Development. His research interests include analog and digital signal processing for IRST, low-noise electronics, IRST system test and evaluation for small target detection.

**Boohwan Lee** received his BS, MS, and PhD degrees in electrical engineering and computer science from Kyungpook National University, Daegu, Republic of Korea in 1991, 1993, and 2006, respectively. He has worked as a principal researcher for the Agency for Defense Development since 1993. His current interests include digital image processing, target detection, and tracking.