Automatic Inhomogeneous Background Correction for Spatial Target Detection Image Based on Partition Processing

Abstract

High-resolution imaging with wide field of view (FoV) ground-based telescopes is often affected by skylight background and noise due to the detector, resulting in an inhomogeneous background. In this paper, we propose an improved method for spatial image non-uniformity correction based on partition processing. First, an evaluation metric is introduced to evaluate the partition size and automatically iterate a suitable partition value for different scenarios based on the different operating conditions of the telescope. Then, we iteratively calculate the mean and variance in each partitioned region to extract the background of each partitioned region. Finally, after applying bilinear interpolation to the background extracted from each region, the inhomogeneous background is obtained and removed from the original image. The experiments on the simulated and real images show that the proposed method can effectively remove the inhomogeneous background of spatial images and meet the requirements of the real-time processing of high-resolution images under long exposure conditions.

1. Introduction

Following the growing number of space debris, space target detection becomes a complex and challenging subject. Ground-based photoelectric telescopes are the most commonly used equipment in space target detection; several organizations have carried out observations with ground-based photoelectric telescopes and published their results [1,2,3,4,5]. Meanwhile, the removal of the inhomogeneous background from their captured images is an important prerequisite for improving space target detection capabilities. Image degradation not only affects image quality and the image signal-to-noise ratio (SNR), but it also has a serious impact on the subsequent image segmentation and target detection. Therefore, the inhomogeneous background correction of spatial images becomes a necessary pre-processing step. In addition, for applications where small telescopes search for space targets, remote unattended automatic observation with multi-station deployment can be achieved if the equipment has the characteristics of a lower cost, automatic observation, and easy deployment [6].

Methods of geometric modification, image smoothing and filtering, image morphological operations, and image arithmetic operations are used to process spatial images taken by ground-based telescopes. Each of these operations enhances or highlights the region of interest for later viewing, processing, etc. However, image processing operations can cause some damage to the original information of the image, which is known as the process of image degradation.

Although ground-based telescope systems can use adaptive optics or phase modulation to improve image quality, these systems are usually expensive and complex [7,8,9]. Even more effective image enhancement techniques are applied, such as deconvolution and image processing. It still cannot be achieved in a way that satisfies the real-time processing requirements [10,11,12]. A lot of research has been published in the field of image processing on the issue of inhomogeneous background correction. A simple method is the flat-field correction under uniform illumination [13,14]. Although it can solve the problem in a fixed scene, the correction must be performed again when the observation environment changes. For ground-based telescopes, obtaining a flat-field image under uniform illumination is more problematic. Another way to solve the inhomogeneous background is image processing. Depending on the number of images processed, the methods can be classified as multi-frame image methods and single-frame image methods. The widely used multi-frame image method extracts the inhomogeneous background from multiple images of the same scene. Yuan et al. [15] proposed an improved radial gradient correction method with a better radial uniformity, achieving more accurate details in the corrected images. Litvinov et al. [16] proposed a method to simultaneously estimate the camera radiometric response, camera gain, and vignetting using multiple consecutive images. However, a limitation of the multi-frame processing methods is adapting to the drastic changes in the scene, such as observing a fast-moving target in low earth orbit (LEO).

Single-frame inhomogeneous background correction methods are more flexible because they do not require the accumulation of multi-frame information. The single-frame inhomogeneous background correction studies mainly concentrate on the field of computed tomography, which does not require high real-time processing [17,18]. In the field of spatial target detection, Zhang et al. [19] used maximum likelihood estimation for correction. Partitioning the image into segments is currently the most commonly used method for the real-time correction of an inhomogeneous background in a single-frame image. However, the segmentation threshold is mostly a fixed empirical value, which cannot be adjusted correspondingly over the frames.

This paper presents a single-frame method for automatically removing the inhomogeneous background from the images taken with telescopes. It is more suitable for images captured by small telescopes of an integrated design and meets the real-time processing requirements within 2.5 s at 4 Kpx × 4 Kpx resolution. By introducing a split-size evaluation index, the segmentation value can be adjusted over different scenes. Experiments with simulated images and captured images in this paper show that the method has a better quality and performance compared to the processing method with fixed segmentation thresholds.

2. Equipment and Methods

We first describe the data capture device information in Section 2.1. Next, the workflow of the method is described. The complete workflow of the method is illustrated in Figure 1. The method can be divided into 3 steps, which are the automatic partition, the background-estimating, and the background-removing step. The automatic partition step determines the split value and partitions the image into several regions. The background-estimating step extracts the background within each region. The background-removing step calculates the full image background and then removes it from the original image. The steps will be explained in detail in Section 2.2, Section 2.3 and Section 2.4, respectively.

2.1. Experimental Equipment

The image acquisition equipment for the experiment is shown in Figure 2. The system consists of an optical system, a camera, a focusing mechanism, a motion control system, and an integrated control box. We chose a scientific CMOS detector with a resolution of 4096 px × 4096 px and a pixel size of 9 × 9 μm. The detector was cooled to −40 °C to improve imaging and increase the telescope’s observing power. The telescope’s optical system is designed as a transmissive imaging system, which reduces the size and weight of the optical system. The optical system has an aperture of 150 mm and a field of view of 6.5° × 6.5°. For the integrated design of the telescope, the software system, including the image processing and motion control modules, runs on NVIDIA Jetson AGX Xavier. At the same time, the real-time processing of the telescope is challenged by the high resolution of the images and the low computing power of the processing devices.

The method in this paper is designed to meet the real-time processing requirements of these devices, and we use images taken with each of these devices in a city with a strong skylight background and in a countryside with a weak skylight background. To facilitate the presentation of the results of our method, we store the captured images and then process them offline. At the same time, the images were not compressed in any way to preserve their original information.

2.2. Determination of Region Split Size

The main causes of image degradation are the atmospheric scattering noise and the CMOS detector noise [20]. Mohamed et al. [21] and Sonnenschein et al. [22] pointed out that the noise generated by atmospheric scattering is mainly the Poisson noise. Pain et al. [23] and Nguyen et al. [24] illustrated the sources of CMOS noise, where the dark current noise follows a Poisson distribution, and the readout noise follows a Gaussian distribution. An inhomogeneous image is a non-exposure distorted image. In other words, a large number of pixels are non-saturated, and only the space target or star is the pixel under the maximum grey level. An imaging system with optical system interference can be expressed as:

$${\mathrm{I}}^{\prime }=I+V\times 𝒮+N$$

(1)

where $I$ and ${\mathrm{I}}^{\prime }$ are the pure image and the image with the inhomogeneous background. $V$ is the change caused by the optical system, $𝒮$ is the inhomogeneous background caused by stray light or the detector, and $N$ is the noise of the detector. Assuming that the changes caused by the optical system are slow and uniform, stray light does not cause abnormal exposures in some images. The corrected image can be obtained by the following equation:

$$I={\mathrm{I}}^{\prime }-\mathfrak{B}$$

(2)

where $\mathfrak{B}$ is an inhomogeneous background. Our goal is to determine a method $\mathfrak{F}$ to extract information about the inhomogeneous background from the image.

$$\mathfrak{B}=\mathfrak{F}\left({\mathrm{I}}^{\prime }\right)$$

(3)

Although the distribution of the abovementioned noise is inhomogeneous, the variation in the inhomogeneous background within a certain small image area is very flat. Therefore, we split the image into several small area blocks $R_{block}$. A larger block size can reduce the processing time but it cannot accurately estimate the background of the original image. A smaller split block size can estimate the background of the original image more accurately but requires a correspondingly longer processing time.

Our method determines the split size in step 1. First, a small initial split size is set to split the image $I^{\prime }$. This initial split size is set differently for different scenes and different image resolutions. In this article, the initial split size is set to 10 px. Then, the mean value of pixels within each segmentation block is calculated to obtain the evaluation image $S(m,n)$ of size $m\times n$, which is used to evaluate the segmentation quality. Since the processing is low-pass filtering, the smoothness between each region is a characterization of the changing trend of the high-frequency part of the image. Therefore, an effective smoothness-evaluating operator must perform high-pass filtering on the image. One way to high-pass filter an image is to determine its second-order derivative. For a two-dimensional image, the Laplace operator ${\nabla }^{2}S(m,n)$ can be defined as follows:

$${\nabla }^{2}S(m,n)=\frac{{\partial }^{2}S(m,n)}{\partial m^2}+\frac{{\partial }^{2}S(m,n)}{\partial n^2}$$

(4)

We note that the second-order derivatives of the Laplace operator in the $m$ and $n$ directions may have opposite signs and tend to cancel each other out. This phenomenon occurs frequently on telescope-taken images, meaning the Laplace operator may be unstable. We overcome this problem by modifying the Laplace operator ${{\nabla }^2}_{M}S(m,n)$ as:

$${{\nabla }^2}_{M}S(m,n)=\left|\frac{{\partial }^{2}S}{\partial m^2}\right|+\left|\frac{{\partial }^{2}S}{\partial n^2}\right|$$

(5)

The modified Laplace operator is always larger or equal to the pre-correction Laplace operator. The discretization of a Laplace operator is usually approximated as a 3 × 3 operator. Therefore, the discrete approximation of the modified Laplace (ML) operator can be obtained as:

$$ML\left(m,n\right)=\left|2S\left(m,n\right)-S\left(m-1,n\right)-S\left(m+1,n\right)\right|+\left|2S\left(m,n\right)-S\left(m,n-1\right)-S\left(m,n+1\right)\right|$$

(6)

$ML\left(m,n\right)$ evaluates the undulation between small regions of the discriminated image [25]. Further, the mean ML value between all segmented regions in the evaluation image $S\left(m,n\right)$ is regarded as the evaluation value $E_1$ of the original image at the current segmentation size.

$${\mathrm{E}}_1=\frac{1}{m\times n}\sum _{m,n}\left[\left|2S\left(m,n\right)-S\left(m-1,n\right)-S\left(m+1,n\right)\right|+\left|2S\left(m,n\right)-S\left(m,n-1\right)-S\left(m,n+1\right)\right|\right]$$

(7)

We iteratively scale up the size of the segmented regions with the preset step size then update the evaluation value until the kth evaluation value $E_k$ satisfies the following termination condition.

$$\frac{E_k-E_{k-1}}{E_{k-1}}<0.01$$

(8)

If this condition is satisfied, the rate of change in the evaluation value at this time is small. We consider that the split size is suitable for the current frame.

2.3. Region Background Extraction

In Section 2.2, the mean value of the region is simply taken as the background of the block, but the calculation of the block background should exclude the bright stars and faint points to obtain an accurate background. For star maps taken by telescopes, most of the image elements belong to the background grayscale, and the range is concentrated. At the same time, stars and faint points are presented with high and low grayscale values, respectively. To extract the regional background of the image accurately, the dark and bright spots should be removed when processing the inhomogeneity correction. The removal is performed as follows. First, by calculating the mean value $m_{block}^1$ and the standard deviation ${\sigma }_{block}^1$ in a block.

$$m_{block}^1=\frac{\sum _{(x,y)\in R_{block}}{f}_{block}(x,y)}{N_{block}}$$

(9)

$${\sigma }_{block}^1=\sqrt{\frac{\sum _{(x,y)\in R_{block}}{(f_{block}\left(x,y\right)-m_{block}^1)}^2}{N_{block}}}$$

(10)

where $f_{block}(x,y)$ means the gray value at $(x,y)$ in the region block $R_{block}$, and $N_{block}$ represents the number of pixels in the $R_{block}$. Then, we calculate the segmentation thresholds $T_{block}^{1}L$ and $T_{block}^{1}H$ for the removal:

$$T_{block}^{1}L=m_{block}^1-3\times {\sigma }_{block}^1$$

(11)

$$T_{block}^{1}H=m_{block}^1+3\times {\sigma }_{block}^1$$

(12)

Next, the images are culled according to the following rules:

$$g_{block}^1\left(x,y\right)=\left\{\begin{array}{c}0,f_{block}\left(x,y\right)\le T_{block}^{1}L\\ f_{block}\left(x,y\right),T_{block}^{1}L<f_{block}\left(x,y\right)<T_{block}^{1}H\\ 0,f_{block}(x,y)\ge T_{block}^{1}H\end{array}\right.$$

(13)

$g_{block}^1\left(x,y\right)$ is the image after the first-time removal of bright and dark pixels. Continue analyzing $g_{block}^1\left(x,y\right)$ to calculate the mean $m_{block}^2$ and standard deviation ${\sigma }_{block}^2$ of the remaining points after the removal.

$$m_{block}^2=\frac{\sum _{(x,y)\in R_{block}\ne 0}{g}_{block}^1(x,y)}{N_{block}-N_0}$$

(14)

$${\sigma }_{block}^2=\sqrt{\frac{\sum _{(x,y)\in R_{block}\ne 0}{(g_{block}^1\left(x,y\right)-m_{block}^2)}^2}{N_{block}-N_0}}$$

(15)

where $N_0$ denotes the number of pixels where $g_{block}^1\left(x,y\right)=0$. Our method extracts the block background by removing pixels in each iteration. We update the pixel values iteratively until the image standard deviation satisfies the condition as follows:

$$\frac{{\sigma }_{block}^{n-1}-{\sigma }_{block}^n}{{\sigma }_{block}^{n-1}}<0.01$$

(16)

${\sigma }_{block}^{n-1}$ and ${\sigma }_{block}^n$ denote the standard deviation of the image after the n − 1st and nth removal. When the condition in Equation (16) is satisfied, the standard deviation of the image after the nth removal is close to the previous one, and the dark and bright spots in the image are considered to be removed cleanly. The mean value of $m_{block}^n$ at this point is the gray value of the background in the current region.

2.4. Full Image Background Estimation

There are stepped differences in gray values between the blocks of the background image acquired in Section 2.3. However, the gray value undulation of the background in the actual image is usually continuous. To recover the background of the image accurately, an interpolation filter is required. We choose bilinear interpolation as the method for full image recovery. Compared with nearest neighbor interpolation, bilinear interpolation can provide smoother results. It also performs significantly better than bicubic spline interpolation and meets the real-time requirement.

Similar to the up-sampling, our method takes the center pixel of each region as a sample point, then the estimated background can be obtained by interpolating the gray value over the pixels between the sample points. After, a corrected image can be obtained by differencing the original image from the estimated background.

3. Results and Discussion

In this section, we conduct experiments on the method to verify the effectiveness and practicality of the proposed method. We first introduce the evaluation indicators in Section 3.1. In Section 3.2, we evaluate the effectiveness using the simulated images and then analyze the experimental results. Then, we test our method on the image data captured by the experimental equipment in Section 2.1 and list the results of our method in Section 3.3. The experimental results of both simulated and captured images demonstrate the effectiveness of our method.

3.1. Evaluation Indicators

It should be noted that ground truth images are present in the simulated images but not in the captured images. Therefore, the evaluation indicators are applied to the simulated images. Structural similarity (SSIM) and mean square error (MSE) are used to evaluate the correction quality of the simulated inhomogeneous background. The MSE is defined as follows:

$$\mathrm{M}\mathrm{S}\mathrm{E}=\frac{1}{m\times n}\sum _{i=0}^{m-1}\sum _{j=0}^{n-1}{(I\left(i,j\right)-K(i,j\left)\right)}^2$$

(17)

where $m$ and $n$ are the resolutions of the image, respectively; both of them are set to 4096. $I$ and $K$ are the real images and the corrected image, respectively. The SSIM is defined as follows:

$$\mathrm{S}\mathrm{S}\mathrm{I}\mathrm{M}=\frac{(2{\mathsf{\mu }}_x{\mathsf{\mu }}_y+{\mathrm{c}}_1)\left(\right(2{\mathsf{\sigma }}_{xy}+{\mathrm{c}}_2\left)\right)}{({\mathsf{\mu }}_x^2+{\mathsf{\mu }}_y^2+{\mathrm{c}}_1)({\mathsf{\sigma }}_x^2+{\mathsf{\sigma }}_y^2+{\mathrm{c}}_2)}$$

(18)

where ${\mathsf{\mu }}_x$ and ${\mathsf{\mu }}_y$ represent the average of the input $x$ and the output $y$, ${\mathsf{\sigma }}_x$, and ${\mathsf{\sigma }}_y$ represent the standard deviation of $x$ and $y$, and ${\mathsf{\sigma }}_{xy}$ is the covariance between $x$ and $y$. ${\mathrm{c}}_1$ and ${\mathrm{c}}_2$ are the constant to avoid the error when the denominator is zero.

A smaller MSE indicates that the error between the ground truth image and the corrected image is smaller, which means the background noise is better removed. SSIM closer to 1 declares that the ground truth image and the corrected image are more similar, which denotes that the information of the original image is better preserved during the correction.

3.2. Experiment on Simulated Images

We generated a set of simulated images to evaluate the effectiveness of the proposed method. The simulated inhomogeneous backgrounds are generated according to the model [26]. We generated a background-free image and then composite the background onto it to obtain the test image. In this way, the ground truth image, the inhomogeneous background, and the simulated image are obtained. We apply the method in Section 2 to the simulated images to obtain the estimated background and corrected images; the results are shown in Figure 3.

Figure 4 compares the grey-level distribution in Figure 3a,d on the same row (scan line 350). Figure 4a shows the grey level distribution of the ground truth image on the scan line; Figure 4b shows the corrected image on the grey level distribution. We can see that the average grey level distribution of the corrected image remains the same as the ground truth image, and as for the superimposed noise, the method removes it well.

The effectiveness of the proposed method can be evaluated quantitatively on simulated images. We evaluate the effectiveness of the proposed method by comparing it with the method for a fixed empirical split region size. The results in

Table 1. Performance of different partition areas in the test set.

Split Area Size	MSE	SSIM
10 px × 10 px	33.62	0.9437
20 px × 20 px	23.02	0.9813
30 px × 30 px	8.65	0.9827
40 px × 40 px	9.55	0.9805
Ours	7.46	0.9837

show that the proposed method can find a reasonable split region size compared to the method based on empirical values. Additionally, it has a better correction result on the simulated degraded image, thus obtaining a noise-suppressed image.

3.3. Experiment on Captured Images

We have demonstrated the effectiveness of the method in Section 3.2. Different from the simulated images, it is difficult to obtain ground truth images of the captured images to evaluate our method quantitatively. However, the method processing effect can be identified using the improvement of the image quality before and after the method processing.

In Figure 5, the first column shows the captured images, and the second column lists the corrected images by our method. The third and fourth column represent the background estimated by the method and its gray-level distribution, respectively. The grey scale of the image in the figure has been stretched to present the detail. It is worth noting that the observations in rows 1–2 of the figure are set in cities, where the images contain background stray light scattered by the atmosphere due to the high level of city lighting. The observations in rows 3–4 are set in the countryside, where the stray light scattered by the atmosphere is weaker and the background is mainly caused by detector readout noise and diffraction from bright stars. The comparative results show that our method can correct the inhomogeneous background well in both strong and weak sky backgrounds, demonstrating the universality of the method of this paper.

We have also processed the case when a small amount of cloud cover is obscured. As shown in Figure 6, the method also works well, since the variations in the clouds are usually smooth over small areas.

Figure 7a,b displays the captured and corrected images, respectively, Figure 7c shows the background corrected by the method, and Figure 7d is the normalized background distribution. Due to the low elevation angle, more stray light is refracted and scattered by the atmosphere near the horizon, causing large undulations in the vertical direction of the shot. This variation is still relatively smooth in small areas. It can therefore be well removed by the method proposed in this paper.

We also compare the performance of the method in this paper with the single-image correction method of Zhang et al. [19]. It is worth noting that our method is run in MATLAB R2018a with a single thread. Therefore, only one core was used when testing our method. As shown in

Table 2. Comparison of average execution time.

Image-Size	CPU	Average Running Time	Method
2 kpx × 2 kpx	[email protected] GHz	9.29 s	Zhang
2 kpx × 2 kpx	[email protected] GHz	1.25 s	Ours
4 kpx × 4 kpx	[email protected] GHz	2.36 s	Ours

, in the sample average running time of a set of 40-star images, the average running time of our method is less than Zhang et al.’s method, under the same resolution. Meanwhile, for long-exposure images, the processing time requirement (2.5 s) can also be met under 4 Kpx × 4 Kpx resolution.

4. Conclusions

This paper presents a method for the inhomogeneous background correction of single-frame images from small telescopes. The method removes the skylight background and camera noise from the images. The method is a model-free method that can adapt to complex environments and extract background information from a single-frame image. The method has the advantages of a high performance, high correction accuracy, and the ability to deal with the heterogeneity of complex sky backgrounds, showing a strong practicality.

The experiments on the simulated and captured images prove the effectiveness of the method. The method can significantly improve image quality under the conditions of complex noise and complicated skylight backgrounds. It applies to astronomical observations and remote sensing images and is of great value for these applications. It can be applied to agricultural observations, ocean observations, and ground-based observations. This robust method can also be extended for medical detection and analysis and high dynamic range imaging. In future research, we will speed up our method with an optimized C++ implementation or GPU-accelerated computation.