AUTOMATIC MOSAICKING OF SATELLITE IMAGERY CONSIDERING THE CLOUDS

With the rapid development of high resolution remote sensing for earth observation technology, satellite imagery is widely used in the fields of resource investigation, environment protection, and agricultural research. Image mosaicking is an important part of satellite imagery production. However, the existence of clouds leads to lots of disadvantages for automatic image mosaicking, mainly in two aspects: 1) Image blurring may be caused during the process of image dodging, 2) Cloudy areas may be passed through by automatically generated seamlines. To address these problems, an automatic mosaicking method is proposed for cloudy satellite imagery in this paper. Firstly, modified Otsu thresholding and morphological processing are employed to extract cloudy areas and obtain the percentage of cloud cover. Then, cloud detection results are used to optimize the process of dodging and mosaicking. Thus, the mosaic image can be combined with more clear-sky areas instead of cloudy areas. Besides, clear-sky areas will be clear and distortionless. The Chinese GF-1 wide-field-of-view orthoimages are employed as experimental data. The performance of the proposed approach is evaluated in four aspects: the effect of cloud detection, the sharpness of clear-sky areas, the rationality of seamlines and efficiency. The evaluation results demonstrated that the mosaic image obtained by our method has fewer clouds, better internal color consistency and better visual clarity compared with that obtained by traditional method. The time consumed by the proposed method for 17 scenes of GF-1 orthoimages is within 4 hours on a desktop computer. The efficiency can meet the general production requirements for massive satellite imagery.


INTRODUCTION
Image mosaicking is the construction of a large orthoimage by splicing multiple orthoimages which have overlaps between each other.It is an important part of remote sensing imagery production (Pan, 2015;Cresson, 2015;Chen, 2014a).However, the existence of clouds is a problem that cannot be ignored.Clouds not only cover the ground information, but also change the spectral characteristics of the image.Existing mosaic methods for satellite imagery are often proposed based on the assumption that the orthoimages are cloud-free.When dealing with cloudy images, the assumption may cause a lot of problems mainly in two aspects: 1) Image blurring may be caused during the process of image dodging.Both visual effects and the capability of quantitative assessment were greatly affected.2) Cloudy areas may be passed through by automatically generated seamlines.Many ground objects were obscured by clouds.
Generally, images with high cloud cover will be removed artificially or automatically in the preprocessing stage.However, this treatment may be unreasonable considering that the clearsky areas of these images still contain lots of ground information.Multi-spectral synthesis (Amato, 2008 andIris, 2006), texture analysis (Yuan, 2015;Hégarat-Mascle, 2009) and Homomorphism-Filtering (Wu, 2013) were widely used for automatic cloud detection.However, land resources satellites imagery, such as WorldView-2, Pléiades-1, and Chinese GF-1 imagery, have limitations in wave band and spectral range (Li, 2012).Existing methods are hardly able to meet the production requirements simultaneously on applicability, accuracy and efficiency.On the other side, it means huge costs of manpower and time if we trace clear-sky areas of orthoimages or modify seamlines of the mosaic image artificially for optimal results.Therefore, it is practical and important to explore an automatic mosaicking method considering the clouds for land resources satellite imagery.
In this paper, modified Otsu thresholding and morphological processing were employed to extract cloudy areas and obtain the percentage of cloud cover, then cloud detection results were used to optimize the process of dodging and mosaicking.Details of this method are discussed in Section 2. The Chinese GF-1 wide-field-of-view (WFV) orthoimages are employed as experimental data.The performance of the proposed approach was evaluated in four aspects: the effect of cloud detection, the sharpness of clear-sky areas, the rationality of seamlines and efficiency.Details of this experiment are discussed in Section 3.

Modified Otsu thresholding
Otsu is also known as maximum variance between clusters.Otsu thresholding divides the image into two parts with the principle of maximizing the variance between target and background (Otsu, 1979).Otsu thresholding can be used for cloud detection and extraction under the precondition that the image is cloudy.But in a general workflow of satellite imagery production, quality control of the images was executed before geometric correction.Therefore, orthoimages are generally with no or few clouds.In this condition, miscalculations will be caused in clear-sky pixels if Otsu thresholding was directly used.To solve this problem, qualifications are set for Otsu thresholding.Steps are as follows:

1) Get previous knowledge of cloud-free images:
For images received by a specific satellite sensor, such as Chinese GF-1 WFV images, adequate (more than 200) cloudfree images were selected artificially as a sample set.Images of different areas and different seasons should be collected to ensure its representativeness.Furthermore, homogeneous images (such as ocean images) should be avoided.Then, a Gaussian Mixture Model was employed for fitting the grey histogram of each cloud-free image: (1) where GMM(x) = the Gaussian Mixture Model M = the number of Gaussian-components. ( | , ) i i G x µ σ = the i-th Gaussian-component i τ = weight of the i-th Gaussian-component i µ = mean value of the i-th Gaussian-component i σ = standard deviation of the i-th Gaussian-component Generally, the histogram can be precisely fitted by setting M to 5, and EM (Expectation-Maximum) algorithm was employed for iterative computation until parameters convergence (Permuter, 2006;Bilmes, 1998).Then, [ 1.3 , 1.3 ] was defined as the "Intensive-section" of ( | , ) According to the probability density function of Gaussian distribution, the probability within "Intensive-section" is about 80%.And then, [Gmin, Gmax] was defined as the "Intensivesection" of the histogram by merging "Intensive-section" of all the Gaussian components.In figure 1, Gaussian Mixed Models are represented as red curves; Gaussian components are represented as green curves; "Intensive-section" of Gaussian components are represented as orange curves; Gmin and Gmax are represented as blue lines and blue numerals.

2) Get qualifications of Otsu:
After calculating the "Intensive-section" of all sample images, the minimum value of all Gmax is used as the qualification for Otsu thresholding, recorded as Gini: where G = digital number For a multispectral image, the qualifications of blue, green and red band should be recorded respectively as Gini-B, Gini-G and Gini-R.Most of the clear-sky pixels were ruled out and miscalculations were significantly reduced through this qualification.
3) Otsu thresholding: Otsu threshold was calculated under the limit of {G>Gini}, recorded as O.For a multispectral image, the Otsu thresholds of blue, green and red band should be recorded respectively as OB, OG and OR.It is worth noting that for images received by a specific satellite sensor, section 1) and 2) are only needed once.

Morphology-based integration of cloudy areas
Clouds could be differentiated from most of the ground objects (such as water, soil and vegetation) through the modified Otsu thresholding.But miscalculations may be caused in man-made objects with high brightness, such as buildings and roads.Therefore, images of which percentage of high brightness pixels were less than 1% could be considered as cloud-free images, and the rest of the images needed a morphology-based integration of cloudy areas.Steps are as follows: 1) Erosion: eliminate man-made features Buildings and roads have a size much less than clouds.Morphology erosion can be used to eliminate these smaller areas (Benediktsson, 2003).The size of structural element adopted by the morphology erosion (SE-1) was determined as follows: where Gsd = image ground resolution Man-made objects with a size less than 200 metres would be eroded according to formula (3).For GF-1 WFV imagery (Gsd = 16m), 13*13 pixels sized structure element was adopted to operate the initial segmentation result (figure 3-c).Processing result is represented in figure 3-d.

2) Dilation: fill gaps of cloud
Pixels in the gaps of clouds have little value since the gaps are small and usually covered by mist or shadow.Morphology dilation can be used to fill gaps and integrate cloudy areas.The size of structural elements adopted by the morphology dilation (SD) was determined as follows: Gaps of clouds with a size less than 2000 metres would be filled according to formula (4).For GF-1 WFV imagery, 125*125 pixels sized structure element was adopted.Processing result is represented in figure 3-e.
3) Erosion once more: shrink cloudy areas Morphology dilation filled gaps of clouds, but it also led to the amplification of cloudy areas.Erosion was used once more to shrink cloudy areas in order to keep accuracy.The size of structural elements (SE-2) was determined as follows: For GF-1 WFV imagery, 51*51 pixels sized structure element was adopted.Processing result is represented in figure 3-f It can be seen that misjudgements of man-made objects were eliminated, gaps of clouds were filled, and cloudy areas were integrated through morphological processing.

Image dodging considering the clouds
Image dodging is a vital step to achieve radiometric consistency among images for mosaic (Chen, 2014b).Wallis filter based dodging method is often used due to its superiority and efficiency (Li, 2006).Wallis filter adjusts the grey mean and variance to a given mean and variance (Sun, 2008).In a Wallis process, a desirable image is selected as the standard image while its mean value and variance were set as the standard values.Mean values and variances of other images were adjusted to standard values.The formula is as follows: For cloud-free images, desirable effect can be obtained through Wallis transform.But for images with clouds, mean and variance were affected by clouds and image blurring may be caused during the process of Wallis transform.To address this problem, only clear-sky pixels were used in the establishment of the Wallis filter.In formula (6), mean and variance of clear-sky areas, instead of the whole image, will be adopted.

Seamline determination considering the clouds
Seamline determination is one of the key steps in orthoimage mosaicking and has an important impact on the quality of the final mosaic (Pan, 2015).Area Voronoi Diagrams is often used to generate a seamline network (Pan, 2009), but the seamlines are unable to bypass the clouds.Cloud detection results are used to optimize these seamlines: 1) Clear-sky pixels were preferred in the process of mosaicking, 2) A high priority was set to the images whose percentage of cloud cover is lower, because the images with high cloud cover are always along with mist and cloud shadows that are hard to be detected.In summary, if the position (i, j) is covered by several orthoimages, the pixel P(i, j) will be determined as follows:

EXPERIMENT
We use Chinese GF-1 wide-field-of-view (WFV) orthoimages as experimental data.GF-1 is the first satellite of the high resolution land-observation system of China launched on April 26th, 2013.WFV images achieved by GF-1 satellite are used for large-scale observation with the ground resolution of 16 metres and the swath width of 200 kilometres.WFV images are multispectral images with 4 bands (blue, green, red and near-infrared) and a spectral range of 0.45~0.90micron.With the support of ground control points and the Shuttle Radar Topography Mission, the geometric distortion in WFV orthoimages is less than 2 pixels (Zhang, 2015).
In this paper, 17 orthoimages located in Henan Province, China are employed.It can be seen from figure 6 It can be seen that two mosaic images covering whole province were generated, both of which have a better radiometric consistency compared with the original images.But the traditional mosaic result is covered by clouds in some areas and is blurred in clear-sky areas, while the improved mosaic result is cloud-free and clear.The following paragraphs in chapter 3 will evaluate the performance of the proposed approach from four respects: the effect of cloud detection, the sharpness of clear-sky areas, the rationality of seamlines and efficiency.

Effect of cloud detection
For 17 orthoimages, the percentage of cloud cover after threshold segmentation (intermediate result) and morphological processing (final result) are recorded in table 1

Sharpness of clear-sky areas
According to the results of cloud detection, image-14 (upper image in figure 8-a-0) was selected as the standard image for Wallis processing.Image-09 (under image in figure 8-a-0) was taken as an example and traditional and improved Wallis transform were performed respectively.The influence of image blurring is evaluated qualitatively (figure 8) and quantitatively (table 2).Two parameters were employed: Average Grads (marked as g ∇ below) can be seen as a reflection of image sharpness and Shannon Entropy (marked as H(X) below) can be seen as a reflection of the richness of image details (Wu, 2013).Generally, a blurred image has a smaller g ∇ and H(X) compared with a clear one.Formulas are as follows: where R = image rows C = image columns 2 ( , ) i f i j ∇ = horizontal gradient of pixel (i j) 2 ( , ) j f i j ∇ = vertical gradient of pixel (i j)   8 that the result obtained by traditional Wallis transform has a good radiometric consistency on the whole, but the seamlines in clear-sky areas are easy to be observed.What makes the matter worse was that the image got blurred and the details were difficult to distinguish.Comparatively, the result obtained by improved Wallis transform has a better performance in clear-sky areas while some overexposure in cloudy areas.It also can be seen from table 2 that the result obtained by improved Wallis transform has a higher Average Grads and Shannon Entropy in clear-sky areas.Given that the clear-sky areas are more valuable than cloudy areas, the result obtained by improved method is preferable.

Rationality of seamlines
Seamlines based on Area Voronoi Diagrams are represented as yellow lines in figure 9-a-0 and three fragments are displayed behind.Seamlines based on this method are represented as purple lines in figure 9

Efficiency
The traditional method can be divided into two steps: dodging and mosaicking.The improved method can be divided into three steps: cloud detection, dodging and mosaicking.Some manual work was needed to get the qualification of Otsu, but these work was excluded from efficiency statistics considering they are once and for all.That means in table 3, cloud detection included Otsu thresholding and morphological processing only.For 17 GF-1 WFV Orthoimages (about 37 GB), all the processing work was finished on a desktop computer made up by Intel Core-i5 It can be seen that there is no significant difference in time consumption between two methods, though the improved method has an additional step of cloud detection.The reason is that the use of cloud detection results increased the efficiency of subsequent processing.Thus, the improved method can meet the production requirements for massive satellite imagery.

CONCLUSION
An automatic mosaicking method considering the clouds was proposed for satellite imagery.The feasibility and superiority of this method were verified by 17 scenes of Chinese GF-1 WFV orthoimages.The evaluation results demonstrated that the mosaic image obtained by our method has fewer clouds, better internal color consistency and better visual clarity compared with that obtained by traditional method.The time consumed by the proposed method for 17 scenes of GF-1 orthoimages is within 4 hours on a desktop computer.The efficiency can meet the general production requirements for massive satellite imagery.
At the present stage, the improvements against the clouds were just focused on dodging method based on Wallis filter and seamline determination method based on Area Voronoi Diagrams.Further researches are required in optimizing various dodging algorithms and seamline determination algorithms.

Figure 2 .
Figure 2. Otsu thresholds and qualifications (from top to bottom: blue, green and red band) In figure 2, Otsu thresholds O were represented as red points.Qualifications Gini were represented as blue lines.
Morphology-based integration of cloudy areas

(
= the grey value of original image ( , ) f x y = the grey value after Wallis transform C m = the mean of original image C v = the variance of original image S m = the mean value of standard image S v = the variance of standard image

Figure 4 .
Figure 4. Flow chart of modified Wallis transform

Figure 5 .
Figure 5. Flow Chart of Seamline Determination Thus, the seamlines could bypass the clouds and the mosaic image could be combined with more clear-sky areas instead of cloudy areas.
-a that the majority of the original orthoimages are cloudy and the radiometric consistency among them is poor.Traditional mosaic method based on Wallis transform and Area Voronoi Diagrams was carried out as controlled experiments.Traditional mosaic result is represented in figure 6-b and improved mosaic result considering the clouds is represented in figure 6-c. a Original Orthoimages b Traditional Mosaic Result c Improved Mosaic Result Figure 6.Henan's GF-1 wide-field-of-view images Besides, image-02/04/09/11 were selected as models.Original images are shown in figure7-a (from top to bottom).Cloudy areas after threshold segmentation are shown in figure 7-b and cloudy areas after morphological processing are shown in figure 7-c.It can be seen that the cloudy areas after threshold segmentation and morphological processing are accurate.a b c Figure 7. Cloud detection of Orthoimages the maximum grey-scale of image p(k) = proportion of pixels with the grey-scale of k -b-0 and three fragments of the same areas as 9-a-0 are displayed behind.It can be seen that the mosaic image obtained by traditional method has a higher proportion of cloud cover and the seamlines are obvious in cloudy areas.Comparatively, the mosaic image by improved method proposed in this paper is almost cloud-free and the seamlines are unobvious.a-0 Seamlines obtained by Area Voronoi Diagrams a

Table . 1
. Percentage of Cloud Cover

Table . 2
Quantitative Parameters of Clear-sky Areas It can be indicated from figure CPU and Solid-State-Disk.Efficiency statistics is recorded in table 3: