Using a Hybrid of Interval Type-2 RFCMAC and Bilateral Filter for Satellite Image Dehazing

With increasing advancement of science and technology, remote sensing satellite imaging does not only monitor the Earth’s surface environment instantly and accurately but also helps to prevent destruction from inevitable disasters. The changing weather, e.g., cloudiness or haze formed from atmospheric suspended particles, results in low contrast satellite images, and partial information about Earth’s surface is lost. Therefore, this study proposes an effective dehazing method for one single satellite image, aiming to enhance the image contrast and filter out the region covered with haze, so as to reveal the lost information. First, the initial transmission map of the image is estimated using an Interval Type-2 Recurrent Fuzzy Cerebellar Model Articulation Controller (IT2RFCMAC) model. For the halo and color oversaturation resulted from the processing procedure, a bilateral filter and quadratic function nonlinear conversion are used in turn to refine the initial transmission map. For the estimation of atmospheric light, the first 1% brightest region is used as the color vector of atmospheric light. Finally, the refined transmission map and atmospheric light are used as the parameters for reconstructing the image. The experimental results show that the proposed satellite image dehazing method has good performance in the visibility detail and color contrast of the reconstructed image. In order to further validate the effectiveness of the proposed method, visual assessment and quantitative evaluation were implemented, respectively, and compared with the methods proposed by relevant scholars. The visual assessment and quantitative evaluation analysis demonstrated good results of the proposed approach.


Introduction
Remote sensing images play an important role in monitoring Earth's surface environment for weather prediction, forest fire management, and agricultural purposes. More recently, radar systems have been widely used in remote sensing, such as with synthetic aperture radar (SAR), since it produces high-resolution images which can be used in terrain segmentation [1][2][3][4]. On the other hand, optical or infrared instruments are unfortunately compromised by the effects of cloudiness and impacted by sunlight, especially haze in the atmosphere. To further alleviate these issues, improving the low contrast of image colors and increasing the visibility of observed objects is essential.
The dehazing method in unknown depth for one single image has been used by many scholars extensively. Fattal et al. [5] assumed that the transmission map was uncorrelated with the surface image dehazing, but it is difficult in parameter tuning. In practical applications, various types of noise exist in images; thus, the image denoising approaches were already studied in many fields [19][20][21]. Wang et al. [21] proposed three variational models combining the celebrated dark channel prior (DCP) and total variations (TV) models for simultaneously removing haze and noise. However, to denoise the image, the original image should be known in advance.
According to the aforesaid image dehazing methods, the reconstructed image after dehazing is clearer, but there is still a halo or color oversaturation in certain cases. To overcome these issues, this paper proposes an effective dehazing method for one satellite image. The image contrast of the haze region is enhanced, and the edge detail remains in the reconstructed image after dehazing. The major contributions of this study are as follows: (1) An Interval Type-2 Recurrent Fuzzy Cerebellar Model Articulation Controller (IT2RFCMAC) model is used for estimating the initial transmission map of the image. (2) The proposed dynamic grouping concept is added into a conventional artificial bee colony (ABC) algorithm for solving the falling into a local optimal problem. (3) For the halo and color oversaturation resulted from the processing procedure, the bilateral filter and quadratic function nonlinear conversion are used in turn to refine the initial transmission map. (4) In quantitative evaluation, the proposed method reaches better results in the visibility ratio of the reconstructed image edge and the geometric mean rations of the gradient norms of image's visible edge before and after dehazing.
The following text of this paper is described below. Section 2 introduces the proposed dehazing method, including training transmission map data and overcoming halo and color oversaturation. Section 3 shows the experimental results and analyzes the visual assessment and quantitative evaluation, which are compared with the methods proposed by other authors. Section 4 is the conclusion.

The Proposed Dehazing Method
In computer vision, in terms of enhancing the image contrast for image dehazing, the haze model of the general physical principle is used as the architecture of image dehazing, expressed as Equation (1): is the haze image, representing the pixel value of three color channels of each pixel in the image, which are red, green, and blue; variable is the haze-free image; variable is the color vector of atmospheric light; x is the pixel coordinate of the 2D image; ( ) represents the transmission map, describing the partial values free of scattering in the course of light projection through the medium to the lens, intensity ( ) ∈ [0,1].
This study proposes an effective dehazing method for satellite images. The processing procedure is shown in Figure 1. First, the initial transmission map is estimated by the Interval Type-2 Recurrent Fuzzy Cerebellar Model Articulation Controller (IT2RFCMAC) model, and then the transmission map is refined by using a bilateral filter and quadratic function in turn. Afterwards, the atmospheric light color vector is estimated. Finally, the refined transmission map and atmospheric light are integrated to the reconstructed image, so as to obtain a haze-free image.

Initial Transmission Map Estimation Using IT2RFCMAC Model
This section discusses the estimation of the initial transmission map in two parts. This study uses the IT2RFCMAC model as a learning architecture, which acts as trainer and estimator. First, the IT2RFCMAC model and dynamic group artificial bee colony are used for training data. The idea of training data is from the principle of dark channel prior mentioned in [8]. At least one color channel of the three color channels in the haze image has the minimum pixel value, which approaches 0. This pixel value is the estimated value of transmission map of the pixel. According to this principle, the haze image is observed, as shown in Figure 2. In Figure 2a, the region with a low minimum pixel value of three channels represents the color with high saturation, i.e., the region free of haze interference; hence, the estimated intensity value of transmission map in Figure 2b approaches 1. On the contrary, the region with a high minimum pixel value of three channels in Figure 2a represents the region with haze or cloud layer, meaning the color saturation is low. In order to avoid the low saturation of the haze region influencing the full image, the intensity value of the transmission map in Figure 2b will approach 0. According to the above description, the pixel value and transmission map value present inversely proportional linear variation. Therefore, the training data imported into the IT2RFCMAC model are the pixel values of three color channels (R, G, B), and the output is the transmission map value defined accordingly.   The neural network and fuzzy theory have become a popular research subject and have been proved applicable to many fields, such as image processing, control, and recognition. The cerebellar model articulation controller (CMAC) [22] is inspired by the conception of an artificial neural network (ANN) composed of a multilayer network architecture. The CMAC imitates the operation of the human cerebellum, which is characterized by fast learning, good local generalization capability, and easily implemented architecture. The CMAC has been used in many fields [23,24]; however, the CMAC needs large amounts of memory space for high dimension problems, and its ability for approximation of function is relatively poor. In order to enhance the search capability of CMAC, many scholars proposed some methods to remedy the deficiencies in the cerebellar model, such as improved controllers as the hierarchical fuzzy cerebellar model articulation controller (Hierarchical Fuzzy CMAC) [25], fuzzy cerebellar model articulation controller (FCMAC) [26], and recurrent cerebellar model articulation controller (Recurrent CMAC) [27][28][29]. Wang et al. [30] proposed a recurrent fuzzy cerebellar model articulation controller (RFCMAC) for solving single-image visibility in a degraded image.
This section introduces the IT2RFCMAC model proposed in this paper. The IT2RFCMAC model is combined with the interval type-2 fuzzy set, FCMAC model, and recurrent memory method to improve the conventional CMAC, as shown in Figure 4. The linear combination of input variables is used as the consequence of the fuzzy hypercube (Takagi-Sugeno-Kang function, TSK function), implementing the rule of fuzzy IF-THEN: Rule : IF is , and is , and … and is , where rule represents the th fuzzy hypercube (i.e., fuzzy rule), represents the th input, , represents the membership function of the th input resulted from the th fuzzy hypercube, y represents the consequent output resulted from TSK-function, represents the weight of the th fuzzy hypercube in the consequent part, , represents the weight in the th fuzzy hypercube of the th input in the consequent part.  Figure 5. Interval type-2 fuzzy sets.
and represent the center point and width of the th Gaussian membership function of the th input, respectively. The footprint of uncertainty is created by parallel translation of the grade of membership of the Gaussian membership function μ , expressed as upper membership function and lower membership function , defined as follows: and Therefore, the output of each block in layer 2 is an interval of μ , ( ) , μ , ( ) .

Layer 3 (Firing layer):
Each node of this layer gathers together to form the fuzzy hypercube according to the grade of membership of each rule of layer 2, mapping to each node of layer 3. The firing strength ( ) of each hypercube of this layer is obtained by algebraic product operation, where where Where where , , represent constants. Layer 6 (Defuzzification layer): This layer equates each node of layer 5 with each memory of CMAC, mapping the output fuzzy hypercube ( ) of layer 4 to the corresponding memory to form fuzzy output. Afterwards, type-1 fuzzy set is generated by type-reduction operation of interval type-2 fuzzy system using type-reduction [31]. Finally, the output of crisp value obtained by centroid defuzzification is [ ( ) , ( ) ].

The Proposed Dynamic Group Artificial Bee Colony
The conventional artificial bee colony (ABC) algorithm [32] has fast convergence behavior, but it is early to fall into the local optimal solution. Therefore, this study proposes the concept of dynamic group, considering the balance of global search and local search, and the algorithm efficiency is increased. First, all of food sources are ordered according to the fitness, and the best food source is used as the leader employed bees of the first group, the similarity and difference between the other ungrouped food sources and the leader are calculated. The food source more similar to the current leader will be distributed to the current group, playing the role of onlooker bees. When all food sources are grouped, the bees do their own jobs. The proposed dynamic group artificial bee colony (DGABC) algorithm is shown in Figure 6. The processing procedure is described below.

Step1 (initialization):
All parameters of the IT2RFCMAC model are coded as a bee, i.e., each bee represents an IT2RFCMAC model. Each fuzzy hypercube dimension contains Gaussian , , standard deviation , , displacement , of interval type-2, and weights and , of consequent linear equation, as well as the weight , in the recurrent layer (layer 4). The individual bee coding is shown in Figure  7: Step 2 (ordering): When the fitness of each bee is calculated, the fitness is ordered from the best to the worst, and all bee group numbers are initialized to 0, as shown in Figure 8: Step3 (grouping): When the fitness ordering is done, the bee with the best fitness for the moment is taken as the leader of the new group, and the group number is updated. Afterwards, the present leader is taken as group center to define the similarity threshold of the group according to the ungrouped bees, meaning the difference between the bees of the 0th group and the threshold distance and threshold fitness of the current leader is two thresholds. The threshold distance and the threshold fitness difference are expressed as follows: where is the total number of bees, is the bees of the 0th group among all bee groups, is the dimension, is the position of leader of the th group at present, and is the distance threshold of th group.
, is the position of ungrouped bee.
Best Worst where is the best fitness of leader of th group and is the fitness threshold of th group.
The ungrouped bee with the best fitness is set as leader and shown in Figure 9: Figure 9. The bee with the best fitness is set as leader.
The ungrouped bees are calculated again, i.e., the distance difference ( ) and fitness difference ( ) between the bees of the 0th group and current leader, and whether the ungrouped bees are of the group is judged according to the aforesaid thresholds. The distance difference is expressed as follows: If > and > , the bee is very similar to the current leader, and the bee is classified into the group of current leader, and the group number is updated as current group number. The bees mismatching the condition will not be classified into the group, as shown in Figure 10: If there are bees not yet grouped, the bee with the best fitness is selected from the ungrouped bees in step 3 as the leader of a new group. The grouping actions are repeated until all ungrouped bees are grouped. As shown in Figure 11: Step 4 (bee movement): When all the bees are grouped, the bee move search is implemented to find the best food source. In order to increase the efficiency of the algorithm, the traditional move mode of employed bees is improved in this study; hence, the information of optimal solution is added to the global search of employed bees. On the one hand, the employed bees advance towards the best food source. On the other hand, the random search mechanism of conventional ABC is maintained, so as to avoid too fast algorithm convergence which can easily to fall into the local solution. The improved employed bees search is expressed as follows: where ∅1 , and ∅2 , are random numbers within [0,1], , is the position of the employed bee with the best fitness among all employed bees groups, , is the position of the employed bee selected randomly from all employed bees groups.
The strategy of dynamic group is added to the local search stage of onlooker bees. Therefore, the group leader (employed bees) must lead the onlooker bees to implement a group search, the original exploration of onlooker bees is improved to refer to the direction of their group leader. The improved search is expressed as follows: where ∅ , is the random number within [0,1] and , is the position of the group leader.

Transmission Map Refinement Using Bilateral Filter and Quadratic Function
The aforesaid halo and color oversaturation resulted from using IT2RFCMAC model to estimate the transmission map are solved by using a bilateral filter and quadratic function.

Bilateral Filter
In order to inhibit the halo effect on the image, this study uses a smooth filter based on image edge preservation and noise abatement, known as a bilateral filter [33]. The bilateral filter considers the distance and intensity values of two pixels, known as spatial domain and intensity domain of image, respectively. Therefore, the ( ) transmission map estimated by the IT2RFCMAC model is used as the intensity value of pixel in the image. The initial transmission map is used as the target image for bilateral filter processing, and the weights of the center for peripheral pixels are calculated by the following equation: where is the weight of the spatial domain. The closer the adjacent pixel ( ) is to the center pixel ( ), the higher is the weight obtained. On the contrary, the farther the adjacent pixel ( ) is to the center pixel ( ), the lower is the weight obtained. represents the patch center coordinate. represents the coordinate of pixel around patch center pixel.
is the parameter for adjusting the Gaussian weight in the spatial domain; the larger the value is, the more blurred is the image. is set as 10 in this study. ‖ − ‖ represents the Euclidean distance between two pixels in the spatial domain, also known as Euclidean distance.
where is the weight of the intensity domain. The smaller the brightness difference between adjacent pixel ( ) and center pixel ( ), the higher is the weight obtained. On the contrary, the larger the brightness difference between adjacent pixel ( ) and center pixel ( ), the lower is the weight obtained.
( ) represents the intensity value of patch center pixel. ( ) represents the intensity value of pixel around patch center. is the parameter for adjusting the Gaussian weight in the intensity domain; is set as 5 in this study. ‖ ( ) − ( )‖ represents the Euclidean distance between two pixels in the intensity domain. Afterwards, the weight ( ) of the spatial domain and the weight ( ) of the intensity domain are integrated and multiplied by the transmission map value ( ) of the patch center ( ) adjacent pixel and normalized. The inhibition of halo of in the initial transmission map by bilateral filter is obtained, expressed as follows: where Ω represents the patch used. Figure 12a shows the halo phenomenon resulted from the estimated initial transmission map, mostly distributed over the vehicle and building edges. Figure 12b shows the halo smoothing effect of the bilateral filter in this study.

Quadratic Function
When the halo is inhibited by bilateral filter, there is still color over saturation in some reconstructed images. Therefore, this study proposes a quadratic function nonlinear conversion to refine the linear transmission map, expressed as follows: where ( ) is the transmission map processed by bilateral filter; the coefficients a and b adjust the curvature of the curve, set as 0.14 and 0.71, respectively. The value of the coefficient is compensated according to the gray level mean of the full image. is set as 0.14 in this study. According to the transmission map function curve in Figure 13, the horizontal axis is the referenced transmission map intensity [0,1], and the vertical axis is to refine the intensity value according to the corresponding horizontal axis intensity value. First, the horizontal axis transmission map intensity value approaches 0 in Figure 13, equivalent to thick haze, cloud layer, and white object regions in the image. This viewpoint can be obtained from Equation (1). The value in Equation (26) is used to increase the value of the curve approaching 0. The purpose is to keep some haze particles in the image, so that the image looks more natural. Secondly, the increase of the horizontal axis intensity value in Figure 13 controls a very important part from haze to haze-free intervals. It is observed that using a linear method to refine the transmission map cannot overcome the color oversaturation effectively. Therefore, the purpose of using the nonlinear quadratic function curve to refine the transmission map is to make the variation of transmission map more flexible; i.e., in Equation (26), adjustment of smoothness of the overall curve by coefficients a, b . Finally, the region where the horizontal axis transmission map intensity value approximates 1 is a region of the real image, so the transmission map is free of refined adjustment.

Estimation of Atmospheric Light
The atmospheric light is mainly from the light scattered by the atmospheric particles. It is observed in the haze model in Equation (1) that in the overall image reconstruction process, the effect of atmospheric light is a very important factor. According to Equation (1), if the value of the transmission map ( ) approaches 0, the haze image is equal to the color vector of the atmospheric light. According to this inference, the white region, cloud layer, and haze in the image can be the reference for atmospheric light estimation. However, it is found in the experiment that if the estimated atmospheric light value is too large or too small, the image after dehazing may be too dark or too bright; therefore, it is a big challenge to estimate the atmospheric light exactly in the haze image. The common method obtains atmospheric light referring to the method of He et al. [8]. The authors extract the first 0.1% of bright pixels from the image processed by the dark channel prior as the color vector of atmospheric light. Such a practice may obtain the color vector of atmospheric light accurately for satellite images, but it is found in the experiment that the first 0.1% of atmospheric light values of a haze satellite image are quite large. Thus, the image becomes quite dark.
After several times of experimental adjustment, the first 1% of atmospheric light is selected as the atmospheric light color vector, as shown in Figure 14. First, the minimum pixel value of three color channels is extracted from each pixel of image, and the first 1% bright pixels ( ) are extracted from the maximum to minimum pixel values, the pixels correspond to the pixel position in the original haze satellite image. The three-channel pixel values in the pixel position are added up and divided by 1% of pixels to obtain the final atmospheric light color vector. The virtual codes of the atmospheric light estimation method are shown in Algorithm 1.

Restored Image
The refined transmission map is obtained by using the aforesaid bilateral filter and quadratic function, so as to overcome the halo and color oversaturation problems, and the atmospheric light color vector is estimated. The refined transmission map and atmospheric light are used as parameters of the reconstructed image. The reconstructed image is expressed as follows: where ( ) is the refined transmission map, and is the estimated atmospheric light color vector.

Experimental Results
This section indicates that the parameters used in the IT2RFCMAC model for this study include the total number of bee groups, limited scout bees, generation number, and number of hypercubes, as shown in Table 1. The learning curve between the conventional ABC and improved ABC in training process is analyzed, as shown in Figure 15.  30 100 3000 6 It is difficult to adjust the number of hypercubes in the IT2RFCMAC model; therefore, this study uses 4, 6, and 8 hypercubes, respectively, and implements 10 cross validations, respectively. The results are shown in Table 2. The fitness is calculated as follows:

SN Scout bees limit Generation Number of hypercube
where RMSE represents the root mean square error value. According to the results in Table 2, using 4 hypercubes costs less memory space and shorter training time; however, the performance is not good. Using 8 hypercubes consumes more memory space and longer computing time, but the performance is better than using 4 hypercubes. However, according to the best solution, this study adopts 6 hypercubes in the IT2RFCMAC model.
It is observed in Figure 15 that the conventional ABC algorithm falls into the region solution due to the fast convergence of learning. The proposed DGABC algorithm in this study can find the optimal solution effectively by improving the employed bees' global searching ability and onlooker bees' local searching ability in the learning curve diagram. After the aforesaid module architecture and algorithm parameter analysis, the results of the proposed method will be presented by using the experimental images, including hazy satellite images from the Internet and NASA. People are subjective in image observation; hence, it is difficult to define better image quality. For this problem, this study uses two evaluation methods, which are visual assessment and quantitative evaluation, and compared them with the methods proposed by He et al. [8], Long et al. [16], Wang et al. [11], Pan et al. [17], and Wei et al. [12]. Figure 16 shows the visual assessment methods with respect to the original haze image, the stateof-the-art methods of [8,11,12,16,17], and the proposed method. For the images in Figure 16b,e, we used a guided filter to filter the halo of edges for dehazing. The effect is obvious; the region covered with haze can be removed, and the processed image is much brighter. However, a halo effect remains on the edges, e.g., the mountain edge. The image results by Long et al. [16], as shown in Figure 16c, show oversaturation, because the authors used a single patch size to estimate the transmission map. For the processed results in Figure 16d, the dehazing effect on the haze region is good due to fuzzy control for estimating the transmission map and weight strategy for overcoming the halo problem. However, the whole image has visually severe color oversaturation. Wei et al. [12] used linear conversion to adjust the estimated initial transmission map in the processing procedure; however, the light attenuation is a nonlinear problem. The adjustment of the transmission map is difficult to make be flexible which caused the poor performance of the brightness enhancement procedure after dehazing. The processed image, as shown in Figure 16f, was much brighter, the visual effect was unnatural, and the image lost many object colors. Last, in Figure 16g, the color contrast of the proposed method is much better than the other methods. Furthermore, the information covered with haze can be presented in the haze region, and the haze-free region preserves the original colors of the image. [ 16], (d) Wang et al. [11], (e) Pan et al. [17], (f) Wei et al. [12], (g) proposed method.

Quantitative Evaluation
This section uses ℯ and ̅ , two quantitative indexes proposed by Hautiere et al. [34], to evaluate the image visibility and detail after dehazing. The index ℯ evaluates the visibility ratio of the reconstructed image edge which is expressed as Equation (29): where is the total number of visible edges in the dehazing image. is the total number of visible edges in the original image.
Secondly, index ̅ represents the geometric mean rations of the gradient norms of the image's visible edge before and after dehazing. The gradient value of the reconstructed image is obtained by using the patch method, and the index ̅ is expressed as follows: where parameter is the th element of the corresponding set ℘ . ℘ represents the set of visible edges in the reconstructed image, is the th gradient ratio of the reconstructed image to the original haze image.
Based on these two indices, the comparison results of the aforementioned approaches are revealed in Tables 3 and 4. The aim is to inhibit the halo, increase the contrast without saturating so that the visual information will not be lost. The good results are described by high values of ℯ and ̅ . Long et al. [16] and Wang et al. [11] had color oversaturation in the visual assessment after image dehazing; hence, the value of index ℯ is much higher than those of other methods. He et al. [8] and Pan et al. [17] used a guided filter to inhibit the halo, and the image edge details remained. Therefore, the index ̅ had good results, but the visual assessment had the halo phenomenon. Comparatively, the proposed method successfully constrains the halo effect and retains the detail of the image edges. Moreover, the image color of the proposed method does not occur the oversaturation phenomenon.

Discussion
Owing to the purposes of the remote sensing satellite image, such as obtaining information related to the Earth's resources or monitoring the environment of the Earth's surface, the images should be processed accurately and transmitted instantly. The comparison results of the average image processing time using various methods are displayed in Table 5. In Table 5, a guide filter algorithm is able to reduce time complexity. However, both He et al. [8] and Pan et al. [17] reach longer processing times than other methods. Long et al. [16] and Wei et al. [12] use a single patch to estimate the transmission map which is the fastest in comparison. The proposed method uses the IT2RFCMAC model to estimate the transmission map, and then subsequently applies a bilateral filter to refine the transmission map. Although it costs more time, the reconstructed image shows superior results. The processing speed of the proposed method is not able to be implemented in a real-time system. Therefore, to achieve highspeed operation in real-time applications, the proposed method will be implemented on a field programmable gate array (FPGA), and several image transmission compression technologies [35,36] will be also considered. Table 5. Comparison results of the average image processing time using various methods.

Methods
Average image processing time(s) He et al. [8] 17.52 Long et al. [16] 2.58 Wang et al. [11] 8.01 Pan et al. [17] 17.51 Wei et al. [12] 2.55 Proposed method 5.25 After visual assessment and quantitative evaluation, the advantages and disadvantages using various methods are provided in Table 6. He et al. [8] and Pan et al. [17] maintain good image edge details, which means those methods can inhibit a halo successfully. However, because the image processing time is long in [8] and [17], both methods are hard to implement in real-time systems. Long et al. [16] and Wang et al. [11] have color oversaturation issues, which make an image look untrue. The proposed method is able to inhibit halo, and the color is not saturated; the image processing time is also acceptable.

Conclusions
This study proposes a dehazing method for satellite images. The approach estimates the initial transmission map by interval type-2 RFCMAC (IT2RFCMAC) model. The IT2RFCMAC model is treated as a trainer and estimator. The proposed dynamic group artificial bee colony is used to train the learning parameters of the IT2RFCMAC model. Subsequently, to inhibit the halo and color oversaturation after haze processing, the transmission map is refined by bilateral filter and quadratic function. Afterwards, the first 1% of bright pixels are used as the color vector of atmospheric light. Finally, the refined transmission map and atmospheric light are integrated into the reconstructed image. In the experimental results, the constructed image implemented by the proposed approach shows less haze, less halo, and better color balance than those of other methods.
In future works, the proposed method will be implemented on an FPGA for solving the timeconsuming task of image processing, and the image compression for reducing the transmission time will be also considered.