A Thin-Cloud Mask Method for Remote Sensing Images Based on Sparse Dark Pixel Region Detection

Thin clouds in remote sensing images increase the radiometric distortion of land surfaces. The identification of pixels contaminated by thin clouds, known as the thin-cloud mask, is an important preprocessing procedure to guarantee the proper utilization of data. However, failure to effectively separate thin clouds and high-reflective land-cover features causes thin-cloud masks to remain a challenge. To overcome this problem, we developed a thin-cloud masking method for remote sensing images based on sparse dark pixel region detection. As a result of the effect of scattering, the path radiance is added to the radiance recorded by the sensor in the thin-cloud area, which causes the number of dark pixels in the thin-cloud area to be much less than that in the clear area. In this study, the area of a Thiessen polygon (a nonparametric measure) is used to evaluate the density of local dark pixels, and the region with the sparse dark pixel is selected as the thin-cloud candidate. Then, thin-cloud and clear areas are used as samples to train the background suppression haze thickness index (BSHTI) transform parameters, and convert the original multiband images into single-band images. Finally, an accurate thin-cloud mask is obtained for every buffered thin-cloud candidate, via the segmentation of the BSHTI band. Additionally, the multispectral images obtained by the Wide Field View (WFV), on board the Chinese GaoFen1, and the Operational Land Imager (OLI), on board the Landsat 8, are employed to evaluate the performance of the method. The results reveal that the proposed approach can obtain a thin-cloud mask with a high true-value ratio and detection ratio. Thin-cloud masks can satisfy various application demands.


Introduction
Thin-cloud areas in remote sensing images exhibit information regarding both thin clouds and the underlying land surface; the radiometric distortion of the land cover will result in image classification and target detection errors, which restricts land cover-based applications. The thin-cloud mask, which aims to delimit pixels contaminated by thin clouds, is a critical preprocessing step to ensure the accurate utilization of data. Jointly influenced by the thin-cloud thickness, diversity and complexity of obscured land covers, it is less likely to describe thin clouds via a uniform spectral characteristic. In addition, thin clouds are often confused with high-reflective land covers. As a result, it is difficult to utilize the traditional thick-cloud detection methods for thin-cloud masks. To overcome the aforementioned problems, several thin-cloud detection/mask approaches have been successfully developed according to the characteristics of thin clouds. Similar to those in Sun et al. [1], we can classify these approaches into three main categories: the transform method, decomposition method, and dark object method. then developed into a thin-cloud removal method called dark object subtraction (DOS), as presented in [2,3]. Meanwhile, features in dark pixels were also utilized for thin-cloud detection. For example, in the literature [15], images were classified into mutually disjoint sub-blocks of a definite size to search for dark pixels within them. Then, a haze thickness map (HTM), with a resolution identical to that of the pixel, was obtained by cubic interpolation of the dark pixel. Finally, the thin-cloud area could be identified via threshold segmentation of the HTM.
Instead of acquiring accurate thin-cloud masks directly, the methods described above mostly obtain haze/thin-cloud maps that indicate the thickness of the haze or thin clouds, which can be used for thin-cloud removal in later stages with the proper methods, such as the DOS method [2,3], air-light elimination [16], and the fusion-based method [17,18]. However, insufficient or excessive thin-cloud removal can easily result in new radiometric distortions in the thin-cloud removed image. In addition, the increased availability of remote sensing images with similar spatial resolutions, band settings, and short time intervals enables the distorted data incurred by thin clouds to be directly abandoned or dealt with later by users. Therefore, it is of great importance to develop accurate thin-cloud mask methods.
To address the problems mentioned above, a thin-cloud mask approach based on the sparse dark pixel region detection method is proposed in this study. The main contributions of our method are twofold. (1) We introduce a distinctive feature (the density of dark pixels) and a nonparametric measure (the area of the Thiessen polygon) for the thin-cloud mask. Extensive studies have confirmed that dark pixels are widely existent and scattered throughout the clear parts of images, whereas the density of dark pixels, which are sparse in thin-cloud areas, enables them to become a robust feature. This distinctive feature, combined with the nonparametric measure, can automate the thin-cloud mask process. (2) To obtain a pixel-based cloud mask, the BSHTI is used to transform the original image into the transformed space. The image is segmented locally for every thin-cloud candidate, which can suppress the salt-and-pepper noise and separate bright land covers.
The remaining parts of this paper are organized as follows. The principle of this approach will be introduced in Section 2. Then, we will present the detailed implementation of our method in Section 3. Section 4 will give the experimental settings and results. The applicability of our method is discussed in a wide range of settings, such as the different percentages of clouds and the underlying land cover, in Section 5. Finally, the entire study will be analyzed and discussed in Section 6.

Principle
The visual and near-infrared bands, which are commonly used in passive remote sensing, are highly sensitive to clouds. The effects of thin clouds on solar radiation include reflection, transmission, and absorption. Reflection includes two main effects-scattering and specular reflection. The intensity of these effects depends on the thickness of the clouds. In cases where clouds are rather thick, solar radiation is almost entirely reflected; the cloudy area in the remote sensing images is characterized by high reflectance. When the clouds are thin, the main effect is scattering reflection, which causes solar radiation to scatter in an approximately uniform manner in all directions in the visual and near-infrared bands. Consequently, the thin-cloud area pixels in remote sensing images are added with approximately equal scattered radiation values, which are recorded by a remote sensor. Visually, this process is equivalent to mixing a certain proportion of white components into the image. At the same time, only part of the solar radiation penetrates through the clouds (i.e., transmission) to reach the land surface, which results in the reduction of surface reflection received by the sensor. Visually, there is a corresponding contrast in thin-cloud area drops.
Consistent with the above analysis, under the circumstance of thin clouds, radiance received by the sensor can be simplified and expressed as follows: where R refers to the value of radiance (e.g., the DN value) measured by the sensor, A represents the value of the solar radiation received by the Earth's surface, α is the surface reflectance ratio, B and C represents the scattered radiance and the reflected radiance caused by thin-cloud, respectively. The reflection and scattering phenomena of clouds may lead to changes in the transmission direction of solar radiation; this portion of radiance is recorded by the sensor before it mutually interacts with the land cover. This part of the scattering radiation is generally referred to as the path radiance. Since the path radiance is received by the corresponding sensor before it arrives at the Earth's surface, it contains no surface information and serves as a significant source for information loss and passive remote sensing noise. Common dark pixels include mountain shadows, soil (black soil), pure water bodies, dense vegetation, and brightly colored objects. In other words, dark pixels exist widely in richly textured clear images [14,19] and are scattered throughout the whole image. The effects of thin clouds on atmospheric radiance are mainly characterized by scattering. Therefore, all pixels in the thin-cloud area contain approximately equivalent scattered radiations in the visual and near-infrared channels. Since the radiometric quantities of dark pixels also increase, no dark pixels with low gray levels exist in thin-cloud areas. For this reason, knowing whether or not dark pixels exist in local areas can be used to evaluate the existence of thin clouds. In this paper, a thin-cloud mask method based on the sparse dark pixel region detection method is proposed. Novelties for this proposed method are mainly reflected by (1) the proposal of a multiresolution dark pixel detection method, which is suitable for thin-cloud masks; and (2) the introduction of the Thiessen polygon area as a dark pixel density measure, where the candidate for thin-cloud areas can be identified, and the thin clouds are masked locally for each thin-cloud candidate. In this manner, the accuracy of the cloud mask can be improved.

Methods
The thin-cloud mask method presented in this paper included three major steps: dark pixel extraction, sparse dark pixel region detection, and the BSHTI transform and local thin-cloud segmentation (as shown in Figure 1). Below, the implementation details are illustrated. represents the scattered radiance and the reflected radiance caused by thin-cloud, respectively. The reflection and scattering phenomena of clouds may lead to changes in the transmission direction of solar radiation; this portion of radiance is recorded by the sensor before it mutually interacts with the land cover. This part of the scattering radiation is generally referred to as the path radiance. Since the path radiance is received by the corresponding sensor before it arrives at the Earth's surface, it contains no surface information and serves as a significant source for information loss and passive remote sensing noise. Common dark pixels include mountain shadows, soil (black soil), pure water bodies, dense vegetation, and brightly colored objects. In other words, dark pixels exist widely in richly textured clear images [14,19] and are scattered throughout the whole image. The effects of thin clouds on atmospheric radiance are mainly characterized by scattering. Therefore, all pixels in the thin-cloud area contain approximately equivalent scattered radiations in the visual and near-infrared channels. Since the radiometric quantities of dark pixels also increase, no dark pixels with low gray levels exist in thin-cloud areas. For this reason, knowing whether or not dark pixels exist in local areas can be used to evaluate the existence of thin clouds. In this paper, a thin-cloud mask method based on the sparse dark pixel region detection method is proposed. Novelties for this proposed method are mainly reflected by (1) the proposal of a multiresolution dark pixel detection method, which is suitable for thin-cloud masks; and (2) the introduction of the Thiessen polygon area as a dark pixel density measure, where the candidate for thin-cloud areas can be identified, and the thin clouds are masked locally for each thin-cloud candidate. In this manner, the accuracy of the cloud mask can be improved.

Methods
The thin-cloud mask method presented in this paper included three major steps: dark pixel extraction, sparse dark pixel region detection, and the BSHTI transform and local thin-cloud segmentation (as shown in Figure 1). Below, the implementation details are illustrated.

Dark Pixel Extraction
For a thin-cloud region, the electromagnetic radiation intensity (e.g., the DN value, reflectance, and the gray level, which we denoted hereafter as the gray level for simplicity) received by the sensor was contaminated by the path radiance. The path radiance increased the observation value for all of the pixels in the thin-cloud area of the remote sensing image, indicating that no dark pixels with extremely low gray levels existed in the visible and near-infrared bands of that area. However, dark

Dark Pixel Extraction
For a thin-cloud region, the electromagnetic radiation intensity (e.g., the DN value, reflectance, and the gray level, which we denoted hereafter as the gray level for simplicity) received by the sensor was contaminated by the path radiance. The path radiance increased the observation value for all of the pixels in the thin-cloud area of the remote sensing image, indicating that no dark pixels with Remote Sens. 2018, 10, 617 5 of 18 extremely low gray levels existed in the visible and near-infrared bands of that area. However, dark pixels in the clear areas of remote sensing images were widely existent and scattered throughout the image. Based on this distinct feature, the thin-cloud area detection problem could be resolved by detecting areas with sparse dark pixels. Consequently, the definition and extraction of dark pixels were crucial to the successful implementation of this algorithm.
In traditional methods, a pixel whose gray value was smaller than the user defined threshold, such as the smallest percent (P) of pixels, was recognized as a dark pixel. However, large areas of shadows, pure water bodies, and dark colored vegetation were often clumped together and formed dark pixel patches, which had a negative effect when measuring the density of dark pixels. To overcome this problem, a dark pixel was defined in a local window (i.e., an area of w × w pixels) to assure the relatively uniform distribution. Accordingly, the dark pixel extraction method was designed as follows: (1) A cloudy image with n band(s), denoted as I = {B 1 , B 2 , . . . , B n }, needed to carry out the thin-cloud mask, where B j represents the ith band of the image. The minimum values of all of the bands were extracted and a new band was formed, which is represented by B min . For a pixel in the ith row and the jth column, B (i,j) min could be expressed by the following equation: where min represents the minimization operation.
(2) For the band used for B min , a pixel whose gray level was smaller than the user defined threshold, T Dark , was selected as the dark pixel candidate.
T Dark is defined using the following formula: where H(B min ) represents the cumulative histogram of band B min , and j represents the smallest gray level, whose cumulative percentage in H(B min ) is larger than P. The problem for the selection of T Dark is transformed to the determination of P, where P indicates a pixel has the potential to become a dark pixel. As mentioned above, the image content varied; some images contained large portions of dark pixels clumped together; therefore, it was unreasonable to use an identical P for all of the images. Considering this problem, a multiresolution method was utilized in this paper to define the value of P. The image was divided into a group of non-overlapped patches, which were expressed as R = {R 1 , R 2 , . . . , R m }, where R i represents a patch with a size of i. To eliminate the negative effects of large areas of dark pixels clumped together, only one pixel in the center of the patches was used to calculate the cumulative histogram.
For the resolution of R i , the threshold T Dark (R i ) of the dark pixels was determined via Equation (3). Lastly, the final threshold of T Dark was obtained by selecting the maximum threshold of all the resolutions: On this basis for threshold determination, the pixels of B min ≤ T Dark were selected as dark pixels candidates.
(3) For a dark pixel candidate along the ith row and jth column, we defined a window surrounding the test pixel with a size of w (i.e., a w × w size). When this pixel was the smallest in the window, it was identified as a dark pixel and flagged with 1; otherwise, it was flagged with 0. The dark pixel bands that were acquired are denoted as B Dark : The dark pixels obtained by this method were defined as feature points, assuming that there were m feature points extracted, which is expressed as follows: Therefore, only one dark pixel was selected in the local widow to ensure that all these feature points were uniformly distributed with similar densities in the clear areas. However, in the thin-cloud areas, there were either none or only a few dark pixels, which resulted in a low density. Figure 2 shows the dark pixels extracted by the method that was proposed in this study (w = 7, P = 30%). It could be observed that the feature points of the dark pixels had a large quantity and were more uniformly distributed in urban areas and in areas of dense vegetation. By contrast, none or few dark pixels had been found in thin-cloud areas. In other words, dark pixels in these areas were sparse. Depending on the above analysis, thin-cloud detection became a problem when extracting regions where dark pixels were sparse.
It should be noted that w and P were two essential parameters that decided whether or not a pixel could become a dark pixel, which determined the quantity and distribution of the dark pixels. The parameter set method and its influence on the accuracy of thin-cloud masks will be discussed in detail in Section 5.1.
Remote Sens. 2018, 10, x FOR PEER REVIEW 6 of 18 Therefore, only one dark pixel was selected in the local widow to ensure that all these feature points were uniformly distributed with similar densities in the clear areas. However, in the thin-cloud areas, there were either none or only a few dark pixels, which resulted in a low density. Figure 2 shows the dark pixels extracted by the method that was proposed in this study (w = 7, P = 30%). It could be observed that the feature points of the dark pixels had a large quantity and were more uniformly distributed in urban areas and in areas of dense vegetation. By contrast, none or few dark pixels had been found in thin-cloud areas. In other words, dark pixels in these areas were sparse. Depending on the above analysis, thin-cloud detection became a problem when extracting regions where dark pixels were sparse.
It should be noted that w and P were two essential parameters that decided whether or not a pixel could become a dark pixel, which determined the quantity and distribution of the dark pixels. The parameter set method and its influence on the accuracy of thin-cloud masks will be discussed in detail in Section 5.1.

Sparse Dark Pixel Region Extraction
Traditionally, point density was mostly measured parametrically based on counting the number of points in a local window. The accuracy obtained by these methods relied on the window shape design and parameters (such as window size), which should be selected carefully. When the feature points were distributed uniformly in a regular shape, all of the approaches could yield favorable results. However, sparse feature point areas formed by thin clouds had irregular shapes, different sizes, and notable differences in their densities. Therefore, it was a challenge to adopt a windowbased feature point density measurement, as failure in the effective estimation of sparse feature point densities might have led to the missing detection of thin clouds or a low accuracy of the thin-cloud boundary. To address the problem presented above, a nonparametric point density measurement method was utilized in this study.
The Thiessen polygon was a commonly used zoning method. In the Thiessen polygon result, the distance from any position in the polygon to its corresponding point was shorter than that to any of the other sides of the polygon. Therefore, only one feature point existed within the range of each Thiessen polygon. In other words, a one-to-one corresponding relation existed between the feature points and their Thiessen polygons. Figure 3 shows the Thiessen polygon corresponding to the

Sparse Dark Pixel Region Extraction
Traditionally, point density was mostly measured parametrically based on counting the number of points in a local window. The accuracy obtained by these methods relied on the window shape design and parameters (such as window size), which should be selected carefully. When the feature points were distributed uniformly in a regular shape, all of the approaches could yield favorable results. However, sparse feature point areas formed by thin clouds had irregular shapes, different sizes, and notable differences in their densities. Therefore, it was a challenge to adopt a window-based feature point density measurement, as failure in the effective estimation of sparse feature point densities might have led to the missing detection of thin clouds or a low accuracy of the thin-cloud boundary. To address the problem presented above, a nonparametric point density measurement method was utilized in this study. The Thiessen polygon was a commonly used zoning method. In the Thiessen polygon result, the distance from any position in the polygon to its corresponding point was shorter than that to any of the other sides of the polygon. Therefore, only one feature point existed within the range of each Thiessen polygon. In other words, a one-to-one corresponding relation existed between the feature points and their Thiessen polygons. Figure 3 shows the Thiessen polygon corresponding to the feature points in Figure 2. It was clear that the area of the polygon was large when the feature points were sparse. According to this fact, a sparse dark pixel region extraction method was designed as follows: Remote Sens. 2018, 10, x FOR PEER REVIEW 7 of 18 feature points in Figure 2. It was clear that the area of the polygon was large when the feature points were sparse. According to this fact, a sparse dark pixel region extraction method was designed as follows: (a) (b) (1) Thiessen polygons, which divided image regions into disjoint subsets, were established according to the feature point set Voronoi diagrams and Thiessen polygons were dual diagrams that interacted with each other, which meant that the Thiessen polygon could be established through the Voronoi diagram. We constructed the Voronoi diagram using a maximum least-angle rule, which could be implemented by multiple classic methods and an improved method.
(2) The Thiessen polygon was divided into two mutually disjoint subsets, which were referred to as the sparse set ( were calculated using the following formulas: where a i ∈ A represents the area of the Thiessen polygon corresponding to the feature point, d i ∈ D, indicates the density of the feature points around point d i .
Voronoi diagrams and Thiessen polygons were dual diagrams that interacted with each other, which meant that the Thiessen polygon could be established through the Voronoi diagram. We constructed the Voronoi diagram using a maximum least-angle rule, which could be implemented by multiple classic methods and an improved method.
(2) The Thiessen polygon was divided into two mutually disjoint subsets, which were referred to as the sparse set (A SP ) and the dense set (A DS ). The former represented the sparse areas of feature points, which indicated thin-cloud candidates, and the latter represented their dense areas, which were also known as clear areas. The initial values of A SP and A DS were set as A 0 SP = ∅ and A 0 DS = A, respectively. The superscript 0 denoted the number of iterations.
(3) The average µ(A 0 DS ) and the standard deviation σ(A 0 DS ) of the polygon for a dense set were calculated using the following formulas: where A 0 DS represents the number of elements in set A 0 DS . (4) The large area of the Thiessen polygon indicated the sparse feature points near the location of the feature point. Considering this, polygons satisfying Equation (10) were treated as abnormal values of a i ∈ A 0 DS and were eliminated from the set A 0 DS .
where λ represents a multiple of the standard deviation, which we set to 3 in this paper. As we selected only one dark pixel in the local window w, whose size was a w × w pixel. Ideally, the area of the Thiessen polygon in the clear part of the image obeyed a normal distribution, whose mean µ(A DS ) approached the number of w × w pixels. By repeating the experiment with different parameters, λ = 3 was deemed as appropriate to separate the dense, dark pixel region from the sparse region. The eliminated polygons were added into set A SP . In this case, two new sets were acquired and denoted as A 1 SP and A 1 DS . (5) A 1 SP and A 1 DS were utilized as initial values to repeat steps (3) and (4); the iteration process was repeated i times until no more polygons met the condition of Equation (10) in set A i DS . (6) The region enclosed by a polygon in the sparse set was recognized as a thin-cloud area candidate. The connected pixels formed a thin-cloud patch and were treated as thin-cloud candidates, which could be expressed as follows: As certain dark pixels existed near the thin clouds, the candidate thin-cloud area was larger than the real thin-cloud area in most cases, which incurred false detection. Meanwhile, the corresponding detection might have missed some thin clouds with shadows, as an improper area threshold was adopted. Therefore, the boundary of the thin cloud should have been further refined by additional methods.

BSHTI Transform and Local Segmentation
As a result of the diversity of the land cover and the various thicknesses of thin clouds, it was less likely to represent thin clouds and backgrounds via a unique feature with wide applicability. In addition, single-band segmentation had drawbacks regarding fully exploiting the information contained in the multiband remote sensing image. According to [8], the BSHTI transform, which synthesized the multiband image into a single band, had the capacity to suppress the background and highlight thin clouds. The new band, B BSHTI , obtained by the BSHTI transform, had a linear weighted sum that was calculated by the multiband images. The expression was given as follows: where K = (k 0 , k 1 , . . . , k n ) represents a series of conversion coefficients for multiple bands. The objective of this transform was to enable the mean value of the thin-cloud area (σ CL (B BSHTI )) in the BSHTI band to be as high as possible, whereas the values in the clear area should have been distributed in a concentrated manner to the greatest extent This meant that the standard variation of the clear area for σ CL (B BSHTI ) should have been as small as possible. This problem could be solved by searching for an optimal solution to Q 1 with a constraint condition of Q 2 : Remote Sens. 2018, 10, 617 9 of 18 The problem could be optimized by the following formula: where K = (k 1 , k 2 , . . . , k n ) represents the series of conversion coefficients for all of the bands, which was the parameter to be determined. S is the correlation coefficient for all of the bands denoted by an n × n matrix. Both µ CL = (µ CL (B 1 ), µ CL (B 2 ), . . . , µ CL (B n )) and µ TC = (µ TC (B 1 ), µ TC (B 2 ), . . . , µ TC (B n )) are mean values for all of the bands in the clear and thin-cloud areas, respectively. On the basis of obtaining the value K, µ TC (B BSHTI ) could be set to 0 for the convenience of solving k 0 . As a supervised method, thin-cloud and clear areas should have been manually interpreted to calculate the statistics. In this study, candidates of the thin-cloud areas and clear area obtained in Section 3.2 were adopted to train the B BSHTI parameters directly. The conversion results of the image in Figure 2 are presented in Figure 4, which improved the capacity to highlight thin-cloud areas and concentrate various background areas.  The false and missing detection of thin-cloud areas might have appeared in regions adjacent to those for candidates of thin-cloud areas. Therefore, thin-cloud areas were iteratively buffered with 1 pixel until the following condition was satisfied: refers to the number of pixels additionally found due to buffering, and represents the number of a candidate cloudy area in the i T C pixels. Equation (16) indicated that the proportion of thin clouds to background pixels approached a 1:1 ratio. Then, the Otsu algorithm was employed to define a segmentation threshold for the area of  The false and missing detection of thin-cloud areas might have appeared in regions adjacent to those for candidates of thin-cloud areas. Therefore, thin-cloud areas were iteratively buffered with 1 pixel until the following condition was satisfied: where N(Bu f (TC i )) refers to the number of pixels additionally found due to buffering, and N(TC i ) represents the number of a candidate cloudy area in the TC i pixels. Equation (16) indicated that the proportion of thin clouds to background pixels approached a 1:1 ratio. Then, the Otsu algorithm was employed to define a segmentation threshold for the area of TC i ∪ Bu f (TC i ) in order to locally separate the thin-cloud area from the clear area. Finally, the thin-cloud mask results were obtained by repeating the morphological operations with a set number of times (in this paper, we set it to four) to remove the rough boundaries and interpolate the missing data.

Data and Experimental Method
Multi-spectral images, acquired by sensors of the Wide Field View (WFV) onboard the Chinese Gao Fen 1 (GF1) satellite and the Operational Land Imager (OLI) onboard the Landsat 8 satellite, were utilized to evaluate the performance of the proposed method. Despite minor differences that existed in the spectral channel settings of these sensors, they all contained near-infrared, red, green, and blue bands. The following processing procedures were performed via these four bands. Five images that were acquired by each sensor were selected individually. For simplicity, a subset of size 2000 × 2000 with typical thin-cloud areas and underlying land surfaces was clipped. The meta-information of the test data is shown in Table 1. As the images obtained by GF1 WFV and Landsat 8 OLI were characterized with 10-bit and 16-bit resolutions, respectively, all of the bands underwent a linear stretch by means of a percent clip (where the parameter was set to 1%) to force the gray level to remain between 1 and 255. This process might have compressed the gray level of the images and removed the information that was beneficial for separating land covers from thin clouds; better transformation may exist, but for simplicity, this paper utilized a simple, yet acceptable method to normalize the different images. w and P were two important parameters utilized in this experiment. In this paper, w = 7 and P = 30% were used for all of the experiments.
To evaluate the method that was proposed in this study, the results were compared with those via the method that was developed by Makarau et al. [15], which was denoted with the haze thickness map (HTM). In that study, dark pixels were adopted to estimate the haze thickness map, which was similar in principle to that presented in our paper, but with a different implementation method. It was implemented by dividing the images into mutually disjointed windows of definite size to search for the minimum values in all of the windows and use them as dark pixels. By utilizing cubic interpolation, a haze thickness identical to the image resolution could be estimated. Then, the haze was detected using a 21 × 21 window, and the mean value of the haze thickness map was utilized to perform the threshold segmentation for the HTM and define the haze area.
For the ground truth cloud-mask acquisition, the images that were acquired by Landsat 8 OLI were distributed with a Quality Assessment (QA) band, which indicated the pixelwise radiometric quality and included information on cloud confidence, cloud shadow confidence, and cirrus confidence. Since the OLI sensor had a cirrus band, supplemental information was provided to separate thin-cirrus clouds from confusing land covers [20]. In the QA band, two bits were used to indicate the radiometric quality; for example, 00 indicated a non-determined situation, 01 referred to low confidence (0% to 33%), 10 represented medium confidence (34% to 66%), and 11 represented high confidence (67% to 100%). We created a thin-cloud mask with the following rule: cloud or cirrus confidences equal to 11 were assigned as clouds; pixels equal to 10 were assigned as suspected clouds, the remaining conditions were treated as clear pixels.
For the GF1 WFV images, we obtained a ground truth thin-cloud mask via visual interpretation. In terms of substantially thin-cloud images, different interpreters (or the same interpreter at different times) might have generated contradictory thin-cloud interpretation outcomes. Considering this, three interpreters were used in this paper to carry out thin-cloud interpretations and labeled thin clouds independently. The shared areas of these results were selected as the ground truth values of thin clouds. In other words, only thin clouds that were identified by all three interpreters could be seen as thin clouds. Thin clouds that were interpreted by one or two interpreter(s) were recognized as suspected clouds. This interpretation adopted a strategy similar to that of people observe objects from a whole scale to a local scale. First, on a small scale, the entire thin-cloud area was obtained. Then, the images were scaled up to 1:1. In this case, the boundaries of thin-cloud areas were modified to closely fit the boundaries of the thin clouds.
The thin cloud results extracted by our method and the contrasting method were compared with the corresponding ground truth values acquired by the interpreters. In addition, three indexes were acquired: C D , which represents thin clouds that were correctly detected, C F , which represents thin clouds that were identified by the algorithm, but with relevant ground truth values that revealed that they were not clouds or thin clouds, and C O , which refers to thin clouds that were not identified by the algorithm. The unit for these three indexes was a pixel. It is worth noting that the suspected thin clouds could have been classified into thin clouds or clear conditions, without the calculation of false or missing detections. On this basis, the three indexes of true ratio (TR), detection ratio (DR), and overall accuracy (OA) were statistically summarized to evaluate the proposed method in this paper and its contrasting methods. The computational methods are as follows:

Experimental Results and Performance Evaluation
The thin-cloud mask results via the proposed method and its contrasting methods are illustrated in Figures 5 and 6. In general, our method was able to detect thin clouds with high accuracy and fewer missing data errors. By carefully observing the boundaries between thin clouds and land surface features, we found that the extracted edge adhered to the true edge between the clouds and the clear part of the image. Because our method detected thin clouds within a large range via the Thiessen polygon, and the extraction of thin-cloud areas was based on the buffering of candidate thin-cloud areas, this method was able to effectively avoid confusion regarding high-reflectance surface features, such as small patches of bare soil and buildings. In this way, no salt-and-pepper noise could be generated. It should be noted that sea surfaces or broad water surfaces, which had large amounts of low reflectance dark pixels, combined with an undetermined threshold, detected a certain amount of false dark pixels, which led to the missing detection of thin clouds.
The HTM used a window to extract dark pixels and estimate the thickness of haze; the large-size thin-cloud detection results were similar to our method. As a result of the influence of window size, however, some missing detection problems could be found in small pieces of the clouds or in thin and long clouds. Similarly, thin-cloud boundaries extracted by the HTM were roughly affected by the window size, and there were mosaic characteristics in the results. Notably, the HTM should have excluded high reflectance objects through the segmentation of buildings and bare soil. Since there were a certain amount dark pixels in urban and bare soil regions, the method proposed in this paper could be adopted to effectively solve confusion regarding high-reflectance surface features, including urban areas.  The HTM used a window to extract dark pixels and estimate the thickness of haze; the large-size thin-cloud detection results were similar to our method. As a result of the influence of window size, however, some missing detection problems could be found in small pieces of the clouds or in thin and long clouds. Similarly, thin-cloud boundaries extracted by the HTM were roughly affected by the window size, and there were mosaic characteristics in the results. Notably, the HTM should have In Figure 7, we present the statistical metrics of the proposed method and its contrasting methods (TR, DR, and OA). We could observe that the TR average in the proposed method reached 93.22%, and its DR reached 88.7%. Generally, the method proposed in this paper was superior to the HTM method for OA and DR with similar TR values. In addition, the results confirmed that the method could be utilized to perform valid thin-cloud masking for images with medium-spatial resolutions and quantization types. On the basis of removing the pixels influenced by thin clouds and preserving clear pixels to the greatest extent, accurate thin-cloud masks that were acquired by this method could be applied to further thin-cloud removal and land cover-based applications.
could be adopted to effectively solve confusion regarding high-reflectance surface features, including urban areas.
In Figure 7, we present the statistical metrics of the proposed method and its contrasting methods (TR, DR, and OA). We could observe that the TR average in the proposed method reached 93.22%, and its DR reached 88.7%. Generally, the method proposed in this paper was superior to the HTM method for OA and DR with similar TR values. In addition, the results confirmed that the method could be utilized to perform valid thin-cloud masking for images with medium-spatial resolutions and quantization types. On the basis of removing the pixels influenced by thin clouds and preserving clear pixels to the greatest extent, accurate thin-cloud masks that were acquired by this method could be applied to further thin-cloud removal and land cover-based applications.

Parameter Setting
The number and distribution of the dark pixels imposed a direct influence on the separation ability of sparse dark pixel areas from the whole image, using threshold segmentation. By determining whether or not a pixel could become a dark pixel, the window size w and the percent threshold P were two key parameters, which should be set carefully by the user. Of these, w represents the size of the local window where a dark pixel exists; P defines the maximum acceptable gray level required by a pixel to become a dark pixel. Excessively low w and large P values could have easily resulted in too many dark pixels; if they were distributed in a thin-cloud area, the corresponding thin area might have been identified as a clear area. On the other hand, when w was too large and P was too small, the number of dark pixels might decrease significantly, resulting in dark pixels that would not be sufficient for separating high-reflectance land surfaces. This caused high reflective land surfaces to be identified as thin clouds.
To evaluate the influence of the parameter setting on the accuracy of the thin-cloud mask, the overall accuracy (OA) of the thin-cloud detection was calculated with a different combination of w and P, where w changed between 7 and 19 with step of 2, and P varied from 5% to 40% with step of 5. Figure 8 shows the OA error bar with changing w and P values; it could be observed that parameter variations within the above range had a minor influence on the OA of the algorithm, as it was simply a tradeoff between the missing detection and the false detection. The results demonstrated that the algorithm was robust over a wide range of parameter settings; therefore, we suggest that the values of w and P could be set to 7 and 30% to obtain an average accuracy with an acceptable true ratio and detection ratio.
The number and distribution of the dark pixels imposed a direct influence on the separation ability of sparse dark pixel areas from the whole image, using threshold segmentation. By determining whether or not a pixel could become a dark pixel, the window size w and the percent threshold P were two key parameters, which should be set carefully by the user. Of these, w represents the size of the local window where a dark pixel exists; P defines the maximum acceptable gray level required by a pixel to become a dark pixel. Excessively low w and large P values could have easily resulted in too many dark pixels; if they were distributed in a thin-cloud area, the corresponding thin area might have been identified as a clear area. On the other hand, when w was too large and P was too small, the number of dark pixels might decrease significantly, resulting in dark pixels that would not be sufficient for separating high-reflectance land surfaces. This caused high reflective land surfaces to be identified as thin clouds.
To evaluate the influence of the parameter setting on the accuracy of the thin-cloud mask, the overall accuracy (OA) of the thin-cloud detection was calculated with a different combination of w and P, where w changed between 7 and 19 with step of 2, and P varied from 5% to 40% with step of 5. Figure 8 shows the OA error bar with changing w and P values; it could be observed that parameter variations within the above range had a minor influence on the OA of the algorithm, as it was simply a tradeoff between the missing detection and the false detection. The results demonstrated that the algorithm was robust over a wide range of parameter settings; therefore, we suggest that the values of w and P could be set to 7 and 30% to obtain an average accuracy with an acceptable true ratio and detection ratio.

Separation of Thin Clouds from Confusing Land Covers
It was well recognized that thin clouds were often confused with high reflective landscapes, such as urban areas or bright bare soils. As a mixed landscape, urban areas were composed of impervious surfaces (e.g., buildings and roads), vegetation, soil, and waterbodies. Of these, dense vegetation, waterbodies, and shadows of tall buildings and trees could become dark pixels, which were scattered in the urban area around the high reflective pixels in a mosaic pattern. Similarly, a certain number of dark pixels existed in the high reflective bare soil area, such as shadows and low-reflection pixels. The proposed thin-cloud detection method captured the unique features when the portion of thinclouds in the remote sensing image were mixed with a white component in the visual and nearinfrared bands; the result was that no pixels with extremely low values existed in these areas. Since this method detected thin-cloud areas with dark pixels instead of using high reflective values, the problem regarding confusion between high reflective landscapes and thin clouds was resolved.

Separation of Thin Clouds from Confusing Land Covers
It was well recognized that thin clouds were often confused with high reflective landscapes, such as urban areas or bright bare soils. As a mixed landscape, urban areas were composed of impervious surfaces (e.g., buildings and roads), vegetation, soil, and waterbodies. Of these, dense vegetation, waterbodies, and shadows of tall buildings and trees could become dark pixels, which were scattered in the urban area around the high reflective pixels in a mosaic pattern. Similarly, a certain number of dark pixels existed in the high reflective bare soil area, such as shadows and low-reflection pixels. The proposed thin-cloud detection method captured the unique features when the portion of thin-clouds in the remote sensing image were mixed with a white component in the visual and near-infrared bands; the result was that no pixels with extremely low values existed in these areas. Since this method detected thin-cloud areas with dark pixels instead of using high reflective values, the problem regarding confusion between high reflective landscapes and thin clouds was resolved. Additionally, the experimental results confirmed that the proposed method possessed the capacity to separate confusing land surfaces from thin clouds. In addition, the method identified thin clouds via local segmentation of the thin-cloud candidate, which suppressed the salt-and-pepper noise that was induced by high reflective landscapes.
As a common approach for thin-cloud detection, this method was robust for both clear images and images composed of clouds. For a clear texture-rich image, where dark pixels were scattered throughout the whole image, the areas of the Thiessen polygons obeyed the normal distribution and were concentrated in the range of the mean, plus or minus three times the standard deviation. Through threshold segmentation of the Thiessen polygon area, few or no sparse dark pixel areas were identified, which caused few or no false detections in the clear texture-rich images. For an image composed of clouds, the dark pixels were sparsely distributed, which exhibited similar features with those in the clear image. However, the minimum gray value or the absolute dark pixel density was larger than that of the clear image, so an additional threshold was required to separate images full of clouds from clear images.
Similar to thin clouds, thick clouds provided missing or distorted landscape information in the remote sensing image. In this respect, there was no need to separate thin clouds from thick clouds. However, semitransparent thin clouds in remote sensing images contained a certain proportion of information about the landscape; then, the thin clouds could be unveiled from the landscape information to obtain a synthetically clear image via an appropriate method (e.g., DOS [2,3], air-light elimination [16], and the fusion based method [17,18]). Since the main objective of our method focused on the identification of the distorted area, there was no clear separation between thin clouds and thick clouds. This method was robust for large, thick clouds. However, this method might have omitted small, thick clouds with corresponding shadows that were distributed near thick clouds.
However, when an image was covered with a large area of waterbody, this method might have failed to mask thin clouds accurately. The spectral signature of clear water was characterized with low reflectance (especially in the near-infrared band) when thin clouds obscure this type of water, and an inappropriate global gray level threshold P. This might have led to the identification of several dark pixels in the thin-cloud area, which resulted in the missing detection problem at the water surface. However, because of the effects of specular reflection and high sediment concentration, several waterbodies were characterized with high-or middle-reflective values, which resulted in none or fewer dark pixels on the smooth and bright water surface. Consequently, false detection could incur.

Conclusions
Thin cloud, which is a common radiometric degradation factor in passive remote sensing, reveals land cover information that is lost or distorted in remote sensing images. Accurate thin-cloud masks, which aim to delineate abnormal areas, are critical to guarantee the proper utilization of partially contaminated remote sensing data, which also serve as a precondition for the reconstruction of thin-cloud areas. Conventional thin-cloud mask approaches, which are combined with thin-cloud removal, require a relatively low thin-cloud mask accuracy. To address the problems described above, a thin-cloud mask method based on the sparse dark pixel area extraction has been developed in this paper. This method possesses the following characteristics: 1.
The scattering effects of thin clouds on solar radiation cause the path radiance to increase in thin-cloud areas, which causes the number of pixels with relatively low values (i.e., dark pixels) in a cloudy area to be much less than those in clear areas. Accordingly, a method has been proposed to extract the candidate of thin-cloud area. Images with diverse backgrounds are used to test this method, and the experimental results indicate that this method is robust under different underlying land covers.

2.
As a result of the window size setting and the shape selection, traditional template-based approaches may fail to extract thin, long, and small clouds effectively. However, the area of the Thiessen polygon, which is a nonparametric parameter, has been selected to measure the density of dark pixels. Because the Thiessen polygon is constructed with the spatial distribution dark pixels, the area of the Thiessen polygon can measure the density of dark pixels in a nonparametric manner.

3.
Since dark pixels and brightly colored objects are distributed in a mosaic pattern, conventionally confusing land surfaces, such as high-reflective urban areas and bare soil, can be effectively separated from thin clouds. Additionally, when the final thin-cloud mask is obtained via the local segmentation of the candidate thin cloud, salt-and-pepper noises can also be suppressed.
Although the proposed method is effective for thin-cloud masks compared with the state of art method, it is far from perfect. The false or missing detection of thin clouds also exists, which further restricts the application of data. Improvement of this method is needed to overcome the problems mentioned above.