Progressive Nonuniformity Correction for Aero-Optical Thermal Radiation Images via Bilateral Filtering and Bézier Surface Fitting

When a high-speed aircraft flies in the atmosphere, the imaging window is subjected to airflow friction, the kinetic energy of the airflow is converted into thermal energy, which makes the surface temperature of the imaging window rise unevenly, the imaging quality is significantly reduced, called aero-optical thermal radiation effects. The continuously increasing nonuniform thermal radiation bias field is not conducive to the precise identification of the target. To remove the bias field in degraded infrared images, this article proposes a progressive nonuniformity correction method. First, we establish a progressive thermal radiation effects correction model to estimate the thermal radiation bias field based on bilateral filtering and Bézier surface fitting. Then, to avoid overfitting the bias field, the degree of Bézier surface is reduced progressively during the iterative process. Finally, according to the properties of the heat transfer of the aero-optical thermal radiation effects, a gradient orientation prior is imposed for both the thermal radiation bias field and the latent clear image. Experiments on simulated degraded images and real degraded IR images show that our method can reduce the thermal radiation effects residual compared with the current aero-optical thermal radiation effects correction methods.


I. INTRODUCTION
W HEN high-speed aircraft flying in the atmosphere, the optical window of the infrared (IR) imaging system will be subject to high-speed airflow of aerodynamic heating, resulting in a nonuniform rise in window temperature. The essence of IR imaging detection is to detect the difference in radiation energy between the target and the background, the noise equivalent temperature difference of the current supersonic IR detection system is generally below 100mK, while the window under high-speed flight conditions, surge and other near-field high-temperature radiation sources are in the hundreds or even thousands of degrees Celsius, which is much higher than the temperature of the target background. The thermal radiation induced by heated optical window can significantly degrade the signal-to-noise ratio and even saturate the detector [1], which is called aero-optical thermal radiation effects [2], as shown in Fig. 1. The imaging results subjected to thermal radiation will be superimposed on the real IR image with a fixed pattern of thermal radiation noise, i.e., thermal radiation bias field, the quality of IR imaging is greatly reduced, and the target detection distance is significantly reduced [2], [3], [4], [5], [6], [7], [8], [9]. In addition, affected by high-speed airflow heating, the optical path outside the optical window will change, resulting in blurring, shaking and moving of the image. These are also referred to as aero-optical effects. This article focuses on the image degradation due to the aero-optical thermal radiation effects, the blurring, shaking and moving of the image can be restored by other image recovery techniques after the thermal radiation is removed.
To orrect the non-uniformity of IR images caused by aerooptical thermal radiation, the current physical methods are broadly classified into three categories: 1) changing the integration time of IR imaging system devices; 2) using heat-resistant materials and thickening the optical hood; and 3) using window cooling techniques. Based on the above three routes, the focal plane optimization technique and the fast infrared imaging This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ technique with adaptive variable frame rate variable integration time significantly improve the IR image quality by changing the integration time [10]. Anthony developed an optical window for harsh aerothermal environments using high purity transparent spinel [11]. Some heat-resistant materials are tested in [8], [12], [13] for thermal radiation suppression. HfO 2 film on the surface of CVD zinc sulfide window and yttrium oxide film on the surface of the sapphire window can delay the saturation of the imaging device and improve the resistance of the optical hood to thermal radiation interference [14]. Window cooling methods mainly include external jet cooling and internal cooling, and the cooling airflow of the external cooling scheme mixes with the airflow within the original boundary layer to form a shear/mixed layer complex flow, which may affect the transmission of target light [15], [16]. Nevertheless, window cooling techniques can effectively improve the aero-optical thermal environment or reduce the window temperature, and are still most effective physical means to suppress the window thermal radiation effects [17], [18]. However, the physical methods listed above are very expensive to implement, and the manufacturing process of optical windows with heat-resistant materials is complicated and the manufacturing period is long. An economically feasible option is to perform a nonuniformity correction for thermal radiation effects on the degraded image.
In the nature image processing field, the approaches for nonuniformity correction are mainly grayscale transformation methods [19], illumination-reflection based homomorphic filtering methods [20], and top-hat transformation based methods, etc. The top-hat transformation methods use morphological filtering for image opening operations, which can significantly destroy image edge information and are mostly used for binarization segmentation of simple images with details and cannot be used for rich detail images. In addition, the methods based on Retinex theory [21], [22], [23] can decompose the image into incident and reflected components by processing the extracted low-frequency illumination information, which also has some effect on aero-optical thermal radiation correction, but removing the illumination information will destroy the contrast of the image at the same time, which makes the dynamic range of the image significantly reduced. The methods mentioned above mainly focus on the non-uniformity of light reflected from image texture under natural illumination, which is affected by the image content information.
Many scholars have also researched in the field of nonuniformity correction of IR images [24], [25], [26], [27], [28], [29]. However, in the past decade, only a few scholars have studied nonuniformity correction casused by thermal radiation. Cao and Tisse proposed a correction method by gradient fitting [6]. Zheng et al. used the sparse property of the image gradient distribution and the segmental constancy of the image to correct the nonuniform intensity of medical images [30], and also had outstanding results in the correction of thermal radiation. Liu and Zhang proposed an nonuniformity correction method based on L 0 and L p regularization constraint [4], [31] for image gradients. Li et al. proposed a correction method based on multiscale bias field estimation [5] with L 1 constraint on the gradient of the corrected image and L 2 constraint on the gradient of the bias field. Shi et al. proposed low-frequency prior and sparse constraints in the gradient domain, combined with Butterworth filtering for low-frequency constraints on thermal radiation images [32]. Subsequently, Shi and Cheng estimated the thermal radiation bias field at multiple scales, and in combination with Chebyshev polynomials to fit the bias field [33], the multiscale information of the degraded image is effectively used. Chang proposed a nonuniformity correction method for thermal radiation effects based on deep multiscale residual networks [34], which is effective in correcting thermal radiation degraded images with fixed bias field scales, but it requires pre-acquisition of thermal radiation bias field radius and intensity for training. In general, although all of the above methods are effective in correcting the degraded IR images, all of them are non-progressive in nature, the corrected images are obtained by subtracting the bias field from the degraded images, and the degraded images remain constant. Therefore, the estimated thermal radiation bias field will directly affect the final correction result. If the thermal radiation bias field is not accurate, the corrected image will still have the residual thermal radiation. In order to accurately estimate the thermal radiation bias field and attenuate the thermal radiation residual, we propose progressive correction with bilateral filtering and Bézier surface fitting for aero-optical thermal radiation effects (BFBSF), the degraded image will be updated by progressive iterative correction. In the process of progressive correction, we use the output of the previous iteration as the degraded input of the next iteration of progressive correction, so that the correction result of the later iteration will contain the information of all previous iterations, the degraded image will gradually become clearer in each iteration, and finally approximate the real clear image. Compared with non-progressive correction, our method can obtain correction results with lower errors from the actual image and less thermal radiation residuals. A random sample of images from remote sensing images in the NWPU-RESISC45 dataset [35] and Middle Wave IR dataset images from [34] were used for thermal radiation simulation [see Fig. 2(a)] and compared with existing thermal radiation correction methods, and the results are shown in Fig. 2(b), where our correction results have the highest peak signal-to-noise ratio (PSNR) and structural similarity (SSIM).
The contributions of this article are summarized as follows: 1) A progressive nonuniformity correction model for thermal radiation is proposed. In contrast to the correction model with a single constant degraded image in previous methods, this article proposes a progressive correction model with degraded image changes, in which the image obtained from the previous iteration will be used as the degraded image for the next iteration, progressively reducing the thermal radiation residuals. To the best of our knowledge, this is the first model that uses progressive correction in the nonuniformity correction task. 2) According to the characteristics of thermal radiation bias field with smoothness and continuity, adaptive parametric bilateral filtering and Bézier surface with variable degree are used to fit the thermal radiation bias field, which improved the generalization ability of this model to the bias field compared with the traditional polynomial fitting.
3) The gradient orientation characteristics of degraded images, thermal radiation bias field and clear images are analyzed, and the gradient orientation regularization terms are proposed for the first time in the thermal radiation correction task to constrain the thermal radiation bias field and the corrected image. The estimation of the thermal radiation bias field is more accurate compared to other advanced methods.

II. PROGRESSIVE THERMAL RADIATION NONUNIFORMITY CORRECTION MODEL
In this article, degraded images are typically modeled using an additional bias field to clear images, and mathematically, a generic correction model for thermal radiation is expressed as where z and f denote the thermal radiation degraded image and the clear image, b denotes thermal radiation bias field, n denotes noise.
The goal of the progressive correction is to continuously separate the thermal radiation bias field b from the degraded image z and then remove it, finally obtain a clear image f , as shown in Fig. 3. Compared with the non-progressive method, the degraded image in the progressive correction will gradually approximate the real image, and the error of the correction result with the real image is lower. Subsequent experiments have demonstrated that progressive thermal radiation correction can remove thermal radiation residue while improving image quality. Based on this, a progressive nonuniform correction model for thermal radiation images is proposed in this article: (2) Noting (2) as BFBSF, as a single iteration correction, the result of the correction is derived from the following equation: In (2) and (3), i denotes the number of iterations,b i is the estimated bias field of thermal radiation for the ith iteration, f i−1 denotes the corrected image of the previous iteration and also serves as the degraded image of the current correction. At the first iteration,f 0 is initialized to the degraded image z.f i is the latent clear image corrected in current iteration, which will be used as the degraded image to enter the next iteration, so the input degraded image will be updated continuously during the subsequent iterations. R(f i ) and R(b i ) respectively denote the constraints on the iterative correction image and the thermal radiation bias field to ensure the correctness of the solution, which will be defined in detail in Sections II-B and II-C. η, α and β are penalty parameters. γ is weight parameter. Bézier m,n (b i ) indicates a Bézier surface fit tob i at m, n degree. Bézier surfaces have good smoothness and can reduce the influence of thermal radiation bias field b i by high frequency noise and image details. Compared with polynomial fitting, Bézier surfaces with Bernstein basis are simple to calculate and less prone to cause runge phenomenon at higher degree. During the iterative process with (3), the latent clear image is updated by progressively removing the thermal radiation bias field, so that the thermal radiation bias field b i converges gradually. When the thermal radiation bias field b i is fully converged, it is unable to continue to separate the thermal radiation bias field from the continuously updated degraded image (potential clear image), at which point f f last .
A. Bias Field Estimation 1) Using Surface Fitting to Solve b: The temperature field distribution on the surface of the optical hood is continuous due to aerodynamic heating [9], [36], [37], [38], [39], [40], [41], so the thermal radiation bias field has smoothness and continuity. In order to make the obtained bias field surface have the above properties and reduce the bias field affected by the details of degraded image, we use Bézier surface fitting to approximate the estimated bias field b i in (2).
First we downsample the estimated bias fieldb, solve the Bézier surface control point coordinates using (4), and then generate the upsampled thermal radiation bias field based on the control point coordinates. In particular, we will use a fitted approximation rather than interpolation, which mitigates the effect of individual pixels on the overall smoothness. The m × n degree Bézier surface equation is expressed as (4) where p(u, v) denotes the point from Bézier surface to be fitted, m, n represent the degree of surfaces,b down denotes downsampling of the estimated bias fieldb, and B i,n (t) denotes the ith basis function of the n degree Bernstein, i.e. B i, , u, v denote the parameterized coordinates of the points to be fitted, which are derived fromb down by the cumulative chord length parameterization method, to scale the x, y coordinate range ofb down to [0, 1]: where l ij = p i−1,j p ij , l * ij = p i,j−1 p ij , (i = 2, 3, . . . , r, j = 2, 3, . . . , s), are the chord lengths in x and y directions.
After parameterizingb down with (5), p(u, v) is known and B i,m (u), B j,n (v) can be derived by computing the basis functions, the problem is transformed into thesolving the control point Q (i,j) , for each coordinate component of where P is the set of p(u, v), (6) can be written as BQ x = P. The x-coordinate component of the control point is found by solving the above equation in the least squares sense, and the y, z components of the control point can be solved in the same way. After obtaining the control point Q, the thermal radiation bias field with the same scale as the degraded image can be generated using uniform parametric according to (4), denoted as b.
2) Using Bilateral Filtering to Solve forb: Since the solution ofb in (2) is not suitable, the solution space searching approach can be used to generate a set ofb based on the low-frequency information of the degraded image to form a solution space (B = {b 1 ,b 2 , . . . ,b n }), and search for theb i in the solution space B that minimizes the loss. Because bilateral filtering can retain or filter high frequency information depending on different parameters, the solution space B is generated from the degraded image using different filtering parameters based on fast bilateral filtering (BF) [42]. Letb = BF(z, θ 1 , θ 2 ), denotes the bilateral filtering of z with parameter θ 1 , θ 2 . Then, we estimate the optimal parameters and search for the optimal solution in the solution space, referring to (2), the loss function and searching approach are defined as A set of optimal parameters (θ 1 ,θ 2 ) can be found which make the loss function (7) obtain the minimum value. After bringing the optimal parameters into (9), the resultb is the solution that satisfies the first line of (2).

B. Gradient Orientation Regularization of Bias Field
In this subsection, the constraint on the bias field R(b) in (2) will be explained in detail. When the high-speed airflow heats the optical hood, if the optical hood is conical, hemispherical or parabolic in shape, the temperature distribution on the surface of the hood is not uniform and there is a heat center due to the physical properties of the optical hood and the aerodynamic heating effects [36], the temperature of the heat center will be higher than the surrounding temperature, so there is a center of thermal radiation on the IR image, and the gradient orientation of the pixels around the radiation center tends to the radiation center, as shown in Fig. 4(d). Since the gradient orientation of the clear image does not have a clear trend, we consider that there is an isotropic thermal radiation bias field in the thermal radiation degraded image with the gradient orientation shown in Fig. 4(f), which makes the gradient orientation of the degraded image have a tendency to converge toward the radiation center.
To make the gradient orientation of the bias field have the trend like Fig. 4(f), a regularization term based on the statistical gradient orientation is proposed: firstly, the image is divided into four blocks according to the radiation center, and the gradient orientation of the four blocks is made to tend to the radiation center by counting the pixel directions that tend to the radiation center in the four blocks and making them maximum. In addition, we introduce the L 1 parametric number for constraining the smoothness of the bias field. The bias field regularization term is defined as where s denotes the sth block of the four blocks, x s , y s denotes the pixel coordinates in the sth block, and |∇b| 1 denotes the L 1 norm of the gradient of b. f (b s , x s , y s , s) is denoted as a Boolean function of the gradient orientation, and is specified as which denotes the pixels in the four blocks with the gradient orientation towards the radiation center as 1 and the rest as 0. θ ∇ (b, x, y), (θ ∇ ∈ [0, 2π]) denotes the gradient orientation of the bias field b at (x, y). f (b s , x s , y s , s) will filter the pixels whose gradients are clustered toward the radiation center in the  sth block in the bias field b. Then, we add pixel filtered from all four blocks. The statistical process is shown in Fig. 5.

C. Gradient Orientation Regularization of Correction Image
The regularization term R(f ) in (2) will be described in detail in this subsection. For the correction images, we also constrain the gradient orientation. Compared with the degraded image, there is no significant trend in the gradient orientation of each block of the clear image, i.e., the gradient orientation distribution tends to be even among the four blocks, as shown in Fig. 6.
Since the edge information in the correction image can interfere with the gradient orientation statistics, for the correction image, we use a voting statistics approach, which is widely used in pedestrian detection [43] and face recognition [44]. The gradient orientation regularization of the corrected image is where hog(f s ) denotes the vector consisting of Histogram of Oriented Gradient (HOG), hog(f s ) denotes the mean vector of hog(f s ). The pixels are divided into cells of size 9 × 9 for voting, to determine the gradient orientation of the cell, and finally HOG statisticsare performed within each block, as shown in Fig. 7.  (7) and (8); Update estimated bias fieldb i by (9) withθ i 1 ,θ i 2 ; Update bias field b i by (4);  (7) and (8); Update estimated bias fieldb i by (9) withθ i 1 ,θ i 2 ; Update bias field b i by (4);

D. Progressive Solution Algorithm
Based on the BFBSF single correction in (3), the progressive BFBSF+ algorithm performs degraded image/latent clear image update by multiple iterations of (3). The progressive BFBSF shortens the algorithm execution time and maintains the stateof-the-art (SOTA) performance by freezing the adaptive filtering parameters on top of the progressive BFBSF+. The progressive BFBSF+ is shown as: The difference between the progressive BFBSF and the progressive BFBSF+ is the update of the filtering parametersθ 1 ,θ 2 . The experimental part is tested on 520 images from degraded dataset, showing that the adaptive filtering parameters fixation can reduce the algorithm running time while maintaining SOTA performance. The progressive BFBSF is shown as:

III. EXPERIMENT
In this section, firstly, comparisons experiment based on the simulated and real degraded images are presented. secondly, the selection of the regularization parameters in (7) will be explored.
In addition, the effect of the Bézier surface degree on the thermal radiation correction results and running time will be explored, and finally, the effectiveness of bilateral filter and surface fitting is tested.

1) Correction of Simulation Images:
To verify the performance of the proposed method. We compare other SOTA methods in this field using PSNR and SSIM. PSNR is defined simply by the mean square error (MSE). Two m × n monochrome images x and y, then their mean square error is defined as and the PSNR is defined as where MAX 2 x is the maximum possible pixel value of the image. SSIM is defined as where μ x , μ y , σ x , σ y , and σ xy are the local means, standard deviations, and cross-covariance for images. and the experimental results are shown in Fig. 8 and Table I. The comparison shows that Cao approach et al. [6] and Shi-19 approach [32] perform well in thermal radiation correction, but introduce some artifacts while correcting, as shown in aircraft image in Fig. 8. Li et al. [5] and Shi-22 [33] approach still has thermal radiation residual in the results, as shown in grassland image in Fig. 8. Zheng et al. [30] approach has reduced the image contrast after correction, as shown in street and residence images in Fig. 8. In contrast, the BFBSF approach reduces the thermal radiation residuals, so that both PSNR and SSIM are greatly improved. At the same time, the thermal radiation bias field fitted by the Bézier surface in this article has good smoothness, so corrected images can retain more gradient information while removing thermal radiation, and corrected images has fewer artifacts.
2) Correction of Real Images: We correct and compare the real thermal radiation degraded IR images with existing SOTA methods, as shown in Fig. 9. All IR images were taken by NASA's WB-57 aircraft [45], [46], [47], [48], IR images from the experiment are released with the code. By comparison, it can be found that BFBSF has better correction results on the nonuniformity correction of infrared thermal radiation, and the residual thermal radiation of corrected IR images is no longer obvious, moreover, there are less artifacts in corrected IR images. To quantify the correction results, we use coefficient of variation (CV) as a non-reference image quality evaluation index, which is defined as a quotient between standard deviation σ and mean value μ of the image x:    We select the smooth area of degraded images in Fig. 9, which is framed by the yellow box for comparison, as shown in Table II.
The degraded image has a higher CV than the corrected image.
Better restoration leads to the lower value of CV [5].

B. Parameters Selection
There are three regularization parameters η, α and β, a weight coefficient γ, and Bézier surface degree m, n exist in (2), (4) and (7), which are assumed to be independent of each other, and the selection of these parameters is analyzed and experimented in the following. In the experiments, simulated thermal radiation bias fields are added to the degraded dataset to generate simulated degraded images, followed by thermal radiation correction using the present method, and PSNR, SSIM are compared with the original images of the dataset to find the optimal parameter values.
1) Regularization Parameter Selection: Fixing two of the regularization parameters η, α and β in turn, the effect of the remaining parameters on the thermal radiation correction is tested on 520 images from the degraded dataset by taking the average PSNR, SSIM of corrected images and degraded images and as shown in Fig. 10.
A large parameter η can lead to poor correction results. The parameter α determines the level of clustering in the gradient orientation of the corrected image, and the parameter β determines the smoothing level of the estimated bias fieldb. Based on the highest PSNR and SSIM, the regularization parameters η, α and β are determined to be 0.2, 0.6 and 0.8.
2) Surface Parameters Selection: The weight coefficient γ and the surface degree m, n of the Bézier surface fitting have a large influence on the correction results. For the weight coefficient γ, if the value is too small or too large, it will affect  the accuracy of the thermal radiation bias field and lead to the PSNR of the corrected image decreasing, as shown in Fig. 10(d).
Combined with the experiment, we take the weight coefficient γ = 0.3 corresponding to the maximum PSNR to obtain the best correction result. Fig. 11 shows the effect of different surface degrees on PSNR and processing time tested on i7-11700 CPU. According to (2), the current iteration of latent clear imagef i is obtained by subtracting b i fromf i−1 in the previous iteration and will be used as the degraded image input for the next iteration. During the iteration, the thermal radiation estimation bias field b i will appear as texture information as the thermal radiation of the iterative image is gradually removed. We compare the processing time and PSNR of Bézier surface fitting with constant degree and with dynamic descending degree, as shown in Fig. 11.
If the higher-degree Bézier surface fitting is continued in subsequent iterations, overfitting will occur, resulting in a decrease in PSNR and a nonlinear increase in processing time. A dynamic decreasing degree surface fitting will avoid overfitting and significantly reduce the processing time. Therefore, we set the degree of the Bézier surface to 6 during the iteration and dynamically reduce the degree of the surface according to the number of iterations.
3) Selection of the Iteration Count: In the progressive thermal radiation correction model, the number of iterations is crucial for the final result. Taking the tennis court in Fig. 8 as degraded image, the Fig. 12 shows the effect of the iteration number on PSNR and SSIM. The PSNR and SSIM increase rapidly from the first to the 10th iterations, and the changing rate of PSNR and SSIM decreases significantly. At the number of iterations close to 20, the changing rate of PSNR and SSIM converge to 0. Therefore, the optimal number of iterations is taken as 20.

C. Effectiveness Study
Using the progressive BFBSF correction as a benchmark, the effectiveness of the regularization term, bilateral filtering, surface fitting, and the effect of the variable degree of the Bézier surface are tested separately, all experiments are performed on Matlab with i7-11700 CPU. The results are shown in Table III. The data fidelity term in the first line of (2) is retained, and the effectiveness of the remaining regularization terms is tested. After removing the R(b) term and the R(f ) term, the PSNR and SSIM of the experimental results are decreased. In addition, the effect of the regularization term on the processing time is not significant.
Since directly using Bézier surface fitting to degraded images can also filter the high-frequency information and conform to the smoothness and continuity of the thermal radiation bias field, direct surface fitting to degraded images without using bilateral filtering can still achieve good results. However, direct Bézier surface fitting cannot guarantee that the fitted sampling points do not contain detail information of the images, and when the fitted points overlap with image detail information, the obtained thermal radiation bias field will still contain some high-frequency information, leading to large deviations in the bias field calculation results. For the degraded image with richer details, bilateral filtering can well remove the high-frequency detail information and retain its low-frequency information. Therefore, it is necessary to perform bilateral filtering on degraded images.
The adaptive parameters of the bilateral filtering in BF-BSF need to be obtained by calculation. Corresponding to Algorithm 1, the progressive BFBSF+ performs the parameter estimation in each iteration and the correction effect has a higher PSNR and SSIM, but it comes at the cost of an extremely long execution time. In the progressive BFBSF, the bilateral filtering parameters are fixed after the first iteration, the PSNR and SSIM of corrected images are reduced in a small way, and the execution time decreases significantly, but the SOTA performance is still maintained in the simulated images of the degraded dataset, corresponding to Algorithm 2.
In our surface fitting effectiveness test, we find that if the surface fitting is abandoned, direct bilateral filtering of the degraded image also has some correction effect. As mentioned above, this is caused by the fact that bilateral filtering itself will try to eliminate high-frequency detail information in the image and retain low-frequency information, which is consistent with the low-frequency characteristics of the thermal radiation bias field. However, the thermal radiation bias field obtained only by bilateral filtering cannot satisfy the continuity and smoothness. The bias field after adding surface fitting is smooth and continuous, which is consistent with the thermal radiation model, so the PSNR and other values after adding surface fitting are further improved. Multi-scale and parallel (MSP) strategies were used in order to speed up the algorithm runtime. First, multi-scale search was used to optimize the searching process of (7) and reduce the search range of the solution space. In the first scale, we scale down the degraded image to 1/8 of the original image and adjust the search step to determine the approximate range of the optimal solution, which is noted as [θ 2 ] to obtain the solution space [θ (2) 1 , θ (2) 2 ]. In the third scale, the degraded image is reduced to 1/2 of the original image, and the search step is further reduced to search for the optimal solution in the range [θ (2) 1 , θ (2) 2 ]. Furthermore, by combining pipelined parallel techniques, we reduce the average correction time of the BFBSF algorithm to the millisecond level without degrading the PSNR.

IV. CONCLUSION
To address the residual problem of nonuniform correction of thermal radiation in IR images, we propose a method to progressively correct degraded images based on adaptive parametric bilateral filtering and Bézier surface fitting. The step of progressive correction makes the correction result gradually approximates the real result, which reduces the residual thermal radiation while improving the quality of IR images. Due to the smoothness and continuity of the thermal radiation bias field, we use adaptive parametric bilateral filtering and Bézier surface fitting in the iterative process to filter the high-frequency information of the degraded image, in order to separate the thermal radiation bias field. In addition, based on the characteristics of degraded IR images, we propose a gradient orientation prior, using the gradient orientation to constrain the corrected image and the thermal radiation bias field. In the experimental part, we test different parameter selections and their effects on the correction results, and investigate the effectiveness of surface fitting and bilateral filtering. Finally, we compare the calibration results of other state-of-the-art methods with ours. Our method achieves optimal results on 520 simulated degraded images from the degraded dataset, and also possesses good correction results on real thermal radiation images. The future improvement direction is to implement the algorithm in hardware to further reduce the algorithm running time, we are working on the following three aspects. 1) Build a multi-DSP, FPGA and ASIC combined architecture; 2) use C language to reconstruct the algorithm, using parallel processing will significantly reduce the algorithm execution time; 3) In real-time correction, the thermal radiation bias field does not change much in a short period of time, so the algorithm can be further optimized by considering the time dependence of the degraded images.