AU2012100257A4 - Method for Radiometric Information Restoration of Mountainous Shadows in Remotely Sensed Images - Google Patents

Abstract

The present invention disclosed a method for radiometric information restoration of mountainous shadows in remotely sensed images, includes shadow region detection, non-shaded similar pixel searching, and shadow information retrieving. This method makes full use of the spectral information derived from both the shaded pixel and its neighboring non-shaded pixels without the aid of DEM, it is an alternative for topographic correction when adequate DEM data are unavailable or spatial resolution and quality of DEM cannot satisfy the requirements of traditional topographic correction algorithms.

Description

1 DESCRIPTION Title Method for Radiometric Information Restoration of Mountainous Shadows in Remotely Sensed Images FIELD OF THE INVENTION The present invention relates to an image processing technique for imagery, especially to a method for the mountainous shadows information restoration. BACKGROUND OF THE INVENTION Shadows in remotely sensed imagery occur when objects totally or partially occlude direct light from a source of illumination, which include cast shadows (shadows cast on the ground by high-rise objects), and self shadows (part of the object not illuminated). Shadows are sometimes used to aid in the reconstruction of three-dimensional geometry, such as measurement of the shape and height of buildings. However, they generate great difficulty in most applications of land cover interpretation and classification because of the reduction or total loss of spectral information of the shaded objects. Particularly, the problem of cast shadows is very severe in mountainous environment with rough terrain, occurring everywhere due to the discrepancy of direct illumination between sunny and shady slopes. Topographic correction is generally used to tackle this problem. Therefore, various algorithms for topographic correction have been proposed and validated in past decades. Teillet proposed the cosine correction, C-correction, and statistical-empirical method for the Lambertian condition; Minnaert correction is one of the most extended methods for non-Lambertian condition, and it was further improved by Richter into Modified Minnaert algorithm with factoring in the influence of various land cover types. In addition to these geometrical optics models, the Image Processing Workbench (IPW) method developed by Dozier is another effective algorithm based on both atmosphere and land surface radioactive transfer models. Various advanced algorithms have also been proposed in the recent years. However, there is a common problem when applying these topographic correction algorithms, that is, DEM data with adequate spatial resolution and quality to satisfy the requirements of topographic correction algorithms are difficult to obtain. For example, DEM data with a 30-meter spatial resolution is essential for topographic correction of Landsat TM/ETM+ data, however DEM data in global scope are currently unavailable. Moreover, the spatial resolution or geometric accuracy of DEM does not precisely match the image under correction even if 2 DEM data are available in most cases. Therefore, errors may be introduced during resampling and geometric registration, which usually result in generating a fragmented and discontinuous image after topographic correction. SUMMERY OF THE INVENTION The objects of present invention are to provide a method for radiometric information restoration of mountainous shadows in remotely sensed images to avoid the above problem. In order to solve the above problem, the present invention disclosed a new method to restore the radiometric information of mountainous shadows in images, includes following steps: A) Shadow region detection: Distinguish shadows and non-shadows according to a fact that shaded areas usually look darker than non-shaded areas if a mountain obstructs the direct light illuminating the surface, First, the Tasseled Cap (TC) transformation was employed to calculate brightness for each pixel, whose brightness was defined as the first band of the TC bands. Then image segmentation based on this calculated brightness band was performed to split an image into separated regions or objects according to spectral and spatial property similarities, After image segmentation, the average brightness value for each object was calculated, and it was taken as the brightness value of the corresponding object. Finally, an automatic thresholding method was employed to find the optimal threshold of brightness where shadow objects can be detected with an acceptable accuracy; B) Non-shaded similar pixel searching: For each shaded pixel, based on the result of image segmentation, a certain amount of non-shaded similar pixels, which could be used to perform the shadow information restoration, were found within a proximate spatial distance to the target shaded pixel. Continuum removal is a brightness normalization technique for spectral absorption feature analysis, where a spectrum is separated into two parts: continuum information (CI) represents brightness, and continuum removed (CR) for spectrum shape information. Shadows were recognized at region scale after image segmentation, for each object, an outward buffer having a width of two pixels was generated around a shadow object, and within the buffer area, N spectrally nearest pixels to a target shaded pixel, where the spectral root mean square deviation (RMSD) satisfies the conditions for spectral similarity, were selected as non-shaded similar pixels. Here, RMSD is calculated as: 3 RMSDI (CR(xi,yi,b) - CR(x,y,b)) 2 n where CR (xi, yi, b) is the CR value of ith non-shaded pixel located in (xi, yi) in band b, CR (x, y, b) shares the same definition but for a target pixel in the shadow region, and n is the number of spectral bands. A large RMSD denotes a large difference in the shape of the spectrum, whereas a small RMSD indicates high similarity between the two pixels. The above mentioned condition for the spectral similarity between target shaded pixel and alternative non-shaded pixel is defined as: RMSDi '(b)X 2 I / n whereo-(b) is the standard deviation of CR value for the whole input image in band b, and m is the number of land cover classes. Moreover, the estimated number of classes (m) needs to be predefined. This value is an empirical parameter and varies with the complexity of the landscape, which can be estimated by visual interpretation of the input images, or using an existing land cover map. In the current study the value of four for m was used. When the number of similar pixels within a width buffer of two pixels does not up to 20, the buffer can be expanded by two pixels repeatedly to search for similar pixels in a larger area until the required number of similar pixels 20 was met. C, Shadow information retrieving: Restoring the DN value of a target shaded pixel using its own CR and CI provided by the non-shaded similar pixels. Hence, the CI of the target shaded pixel can be replaced by the weighted average CI derived from the non-shaded similar pixels to uplift its brightness to a comparable level with the surrounding similar pixels. The following equation represents this process: OIre(b) = Cwavg(b)X CR(h) where OIre(b) is the restored DN value of the target shaded pixel at a specific wavelength (2), or band (b), CR(b) is the continuum removed value of the target shaded pixel at corresponding wavelength( ), or band (b). CIwavg(b) is the weighted average continuum values of non-shaded similar pixels as calculated from: N CLavg(b) = Wi x CI(b) J=1 where CIj(b) is the continuum value of similar pixel j; the weight Wj denotes contribution of similar pixelj to CIwavg(b).

4 Preference in step A, image segmentation can be performed by eCognition Developer 8.0*. In this procedure, shape factor Ssh can be set as 0.1 to emphasize spectral importance, and compactness factor Scm can be set as 0.5 to give equal importance on compactness and smoothness of the region. The most crucial parameter is the scale factor, which controls the region size, can be empirically decided according to the heterogeneous degree of specific areas. Preference in step A, the automatic thresholding method can be described as: the 2 separability of the two classes aB can be maximized when applying the threshold on two separate classes, (B2= 0 0 TP 2 +01( 1PT 2 where po and pi are the mean of the two separated classes, respectively, PT is the total mean of the original data; and wo and wi are the proportion of the total number of pixels with the two classes, respectively. Therefore, the optimal threshold is the one which can induce the 2 maximum B Preference in step B, the maximum buffer width of 10 was set to save on searching time. If the widest buffer still could not provide at least 20 similar pixels, all selected similar pixels were used, regardless of their actual number. The present invention provided a new method to restore the radiometric information of mountainous shadows in remotely sensed images. This method makes full use of the spectral information derived from both the shaded pixel and its neighboring non-shaded pixels without the aid of DEM, it is an alternative for topographic correction when adequate DEM data are unavailable or spatial resolution and quality of DEM cannot satisfy the requirements of traditional topographic correction algorithms. DESCRIPTION OF PREFERRED EMBODIMENT In order that the present invention can be more readily understood, the embodiments of the present invention will be explained in details. The present invention disclosed a new method to restore the radiometric information of mountainous shadows in remotely sensed images, includes following steps: Step A: Shadow region detection, Step B: Non-shaded similar pixel searching, Step C: Shadow information retrieving based on spectral similarity between shaded pixel, i.e., weak information, and its neighboring pixels with similar spectral features but unaffected 5 by shadows, i.e., non-shaded similar pixels. Considering that atmospheric correction probably changes the spectrum of shaded pixels by total or partial loss of their weak spectral information, this new method is recommended for application on raw Digital Number (DN) image before atmospheric correction. Step A: Shadow region detection: Accurate identification of shadows is a prerequisite for shadow information restoration. Many methods are available for shadow detection, such as: Q9Salvador, E., Cavallaro, A., & Ebrahimi, T. (2001). Shadow identification and classification using invariant color models. In IEEE International Conference on Acoustics, Speech and Signal Processing, 3, 1545-1548; @Giles, P. T. (2001). Remote sensing and cast shadows in mountainous terrain. Photogrammetric Engineering and Remote Sensing, 67, 833-839; @Tsai, V. J. D. (2006). A comparative study on shadow compensation of color aerial images in invariant color models. IEEE Transactions on Geoscience and Remote Sensing, 44, 1661-1671; @Chen, Y., Wen, D., Jing, L., & Shi, P. (2007). Shadow information recovery in urban areas from very high resolution satellite imagery. International Journal ofRemote Sensing, 28, 3249-3254. However, these methods are too complex and they need auxiliary data, such as DEM. A relatively simple method was developed to perform shadow detection in the current invention, according to a fact that shaded areas usually look darker than non-shaded areas if a mountain obstructs the direct light illuminating the surface. First, the Tasseled Cap (TC) transformation was employed to calculate brightness for each pixel, whose brightness was defined as the first band of the TC bands. Considering that water body also remains lower brightness which is probably confused with shadows, water body masks should be used to exclude water regions from shadows. A shadow in mountain terrains is often depicted by a certain shape and size, thus, it is reasonable to detect shadows at region scale rather than at pixel scale. After TC transformation, image segmentation can be performed based on calculated brightness band to split an image into separated regions or objects according to spectral and spatial property similarities. Image segmentation can be performed by eCognition Developer 8.0*. In this procedure, shape factor Ssh can be set as 0.1 to emphasize spectral importance, and compactness factor Scm can be set as 0.5 to give equal importance on compactness and smoothness of the region.

6 The most crucial parameter is the scale factor, which controls the region size and should match the required level of detail of a specific study. The scale factor can be empirically decided according to the heterogeneous degree of the study area. The more heterogeneous the study area is, the smaller the scale factor should be set. In this invention, a scale factor optimization method, developed by Chabrier et al. (2006) based on segmentation quality criteria, was used in the current method. The main idea of this method is to select an appropriate scale factor to maximize intra-segment homogeneity and inter-segment heterogeneity. The intra-segment homogeneity can be calculated using the global weighted variance (wVar), whereas intra-segment heterogeneity can be calculated using a spatial autocorrelation measure, Global Moran's I (MI). Weighted variance (wVar) can be defined as: n $ a, * vi wVar = i=1n Y a,

Z=

1 where n is the total number of objects, vi is the variance of spectral value and ai is the area of segment i. Global Moran's I (MI) can be defined as: n$$W~ w,(y - Ayi - Y) MI - =1 j=1 M= i, j where n is the total number of regions, yi is the mean spectral value of region Ri, and Yis the mean spectral value of the image. Each weight wij is a measure of the spatial adjacency of regions Ri and Rj, if regions Ri and Rj are neighbors, wij = 1, otherwise, wij = 0. Combining both the effects of intra-segment homogeneity and inter-segment heterogeneity, a synthetic index "Global Score" (GS) is designed to identify the optimal scale parameter whose optimal scale factor can be achieved as the one with the lowest GS. Global Score (GS) is defined as: GS = V,,,, + MI where Vnorm is the normalized weighted variance and MInorm is the normalized Moran's I. The average of brightness value for each object was calculated after image segmentation, and was taken as the brightness value of the present object. Then, shadow areas can be extracted by an optimal threshold of brightness value generated from an automatic 7 thresholding method. The automatic thresholding method used in the current method is the one developed by Otsu (1979), the related reference literature is: Otsu, N. (1979). A threshold selection method from grey-level histograms. IEEE Transactions on Systems, Man, and Cybernetics, 9, 62-66. The idea of the method is to find an optimal threshold where the separability of the two classes can be maximized when applying the threshold on two separate classes. In the current method, the separability is measured by: ~B 2 0 OG-0 _PT )2 1 1~ (P T )2 Where po and pi are the mean of the two separated classes, respectively, P 1 is the total mean of the original data; and wo and wi are the proportion of the total number of pixels with the two classes, respectively. The optimal threshold is the one which can induce the 2 maximum B Step B: Non-shaded similar pixel searching Shaded pixel in remotely sensed image manifests lower brightness (lower DN-value in each band) than non-shaded pixel because of the lack of illumination as earlier mentioned. However, the spectral shape of shaded pixel can be kept similar with non-shaded pixel belonging to the same class land cover type due to the presence of diffused light scattering. A reasonable method is to separate a spectrum into two parts, one represents brightness information and another one holds spectral shape information, the latter will be helpful to search for non-shaded pixels that have similar land cover type with the shaded one. Accordingly, "Continuum Removal" (CR) method can be used for this purpose. CR is a brightness normalization technique for spectral absorption feature analysis. The related reference literatures are: CDClark, R. N., & Roush, T. L. (1984). Reflectance spectroscopy: Quantitative analysis techniques for remote sensing applications. Journal of Geophysical Research, 89, 6329-6340.; @Youngentob, K. N., Roberts, D. A., Held, A. A., Dennison, P. E., Jia, X., & Lindenmayer, D. B. (2011). Mapping two Eucalyptus subgenera using multiple endmember spectral mixture analysis and continuum-removed imaging spectrometry data. Remote Sensing ofEnvironment, 115, 1115-1128. CR method can separate a spectrum into two parts: continuum information (C]) represents brightness, and continuum removed (CR) for spectrum shape information. In detail, a continuum is fitted over a spectrum to connect the points with local maximum values using 8 a straight line. CR information can be achieved by dividing the actual spectrum value (Original information, 01) at a specific wavelength (2), or band (b), by the value of the continuum line CI at a corresponding wavelength or band. As a result, all points on a spectrum can be normalized at the range of [0, 1], forming a new curve representing the spectral shape of the original. Accordingly, the CR of a shaded pixel should be identical or close to the non-shaded pixels belonging to the same land cover type in spite of the variations in their CIs by applying continuum removal. The separated CR information provides theoretical foundation for similar pixels searching. Particularly in the current method, for each shaded pixel, based on the result of image segmentation, a certain amount of non-shaded similar pixels, which could be used to perform the shadow information restoration, were found within a proximate spatial distance to the target shaded pixel. Shadows were recognized at region scale after image segmentation, for each object, an outward buffer having a width of two pixels was generated around a shadow object, and within the buffer area, N spectrally nearest pixels to a target shaded pixel, where the spectral root mean square deviation (RMSD) satisfies the conditions for spectral similarity, were selected as non-shaded similar pixels. Here, RMSD is calculated as: $ (CR (xi, yi, b) - CR (x, y, b))2 RMSDi= = n where CR (xi, yi, b) is the CR value of ith non-shaded pixel located in (xi, yi) in band b, CR (x, y, b) shares the same definition but for a target pixel in the shadow region, and n is the number of spectral bands. A large RMSD denotes a large difference in the shape of the spectrum, whereas a small RMSD indicates high similarity between the two pixels. The above mentioned condition for the spectral similarity between target shaded pixel and alternative non-shaded pixel is defined as: RMSDi [ [ '(b)x 2/m 7 I In where o(b) is the standard deviation of CR value for the whole input image in band b, and m is the number of land cover classes. Moreover, the estimated number of classes (m) needs to be predefined. This value is an empirical parameter and varies with the complexity of the landscape, which can be estimated by visual interpretation of the input images, or using an existing land cover map. In the current study the value of four for m was used.

9 When the number of similar pixels within a width buffer of two pixels does not up to 20, the buffer can be expanded by two pixels repeatedly to search for similar pixels in a larger area until the required number of similar pixels 20 was met. Meanwhile, the maximum buffer width of 10 was set to save on searching time. If the widest buffer still could not provide at least 20 similar pixels, all selected similar pixels were used, regardless of their actual number. Step C: Shadow information retrieving. According to the spectral relationship between a shaded pixel and its non-shaded similar pixels sharing a similar spectral shape but obtaining different brightness, it is reasonable to restore the DN value of a target shaded pixel using its own CR and CI provided by the non-shaded similar pixels. Hence, the CI of the target shaded pixel can be replaced by the weighted average CI derived from the non-shaded similar pixels to uplift its brightness to a comparable level with the surrounding similar pixels. The DN value of the target shaded pixel was then re-calculated as: 0 r(b) = Cwag() x CR(h) where OIre(b) is the restored DN value of the target shaded pixel at a specific wavelength (2), or band (b), CR(b) is the continuum removed value of the target shaded pixel at corresponding wavelength(t), or band (b). CIwavg(b) is the weighted average continuum values of non-shaded similar pixels as calculated from: N CI,,ag4) = W§. x CI(h) j=1 where the weight Wj denotes contribution of similar pixel j. Each of the non-shaded similar pixels has a different contribution in proportion to both the geographic distance between the similar pixel and target pixel and the spectral similarity between them. Accordingly, the weight can be calculated from the product of the spectral distance (RMSDj) and geographic distance (Dj), where Dj is the Euclidean distance between the jth similar pixel at location (xj, yj) and the target pixel at location (x, y), which can be calculated as: Di= x xf+y -y Considering that the geographic distance is incomparable with the spectral distance due to their different ranges, thus, both spatial distances (Dj) and spectral distance (RMSDj) were normalized and rescaled as the following two equations: D*= D.-Dmi +1 'Dma -Dmi 10 RMSD. - RMSD,, ±1 RMSD RMSD ax- RMSDn where the subscripts "min" and "max" are associated with the minimum and maximum values, respectively; A multiplication operator was then used to combine spatial distance and spectral distance to calculate weight (Wj) as: 1/(D x RMSDI*) W, = N 1(D * x RMSDI*) j=1 The invention provided a new method to restore the radiometric information of mountainous shadows in remotely sensed images. This method makes full use of the spectral information derived from both the shaded pixel and its neighboring non-shaded pixels without the aid of DEM, it is an alternative for topographic correction when adequate DEM data are unavailable or spatial resolution and quality of DEM cannot satisfy the requirements of traditional topographic correction algorithms. Whilst the above has been given by way of illustrative examples of the present invention, many variations and modifications thereto will be apparent to those skilled in the art without departing from the broad ambit and scope of the invention as herein set forth in the following claims.

Claims (4)

1. A method for radiometric information restoration of mountainous shadows in remotely sensed images, includes following steps: A) shadow region detection: distinguish shadows and non-shadows according to a fact that shaded areas usually darker than non-shaded areas if a mountain obstructs the direct light illuminating the surface; first, the Tasseled Cap (TC) transformation was employed to calculate brightness for each pixel, whose brightness was defined as the first band of the TC bands; then image segmentation based on this calculated brightness band was performed to split an image into separated regions or objects according to spectral and spatial property similarities; after image segmentation, the average brightness value for each object was calculated, and it was taken as the brightness value of the corresponding object; finally, an automatic thresholding method was employed to find the optimal threshold of brightness where shadow objects can be detected with an acceptable accuracy; B) non-shaded similar pixel searching: for each shaded pixel, based on the result of image segmentation, a certain amount of non-shaded similar pixels, which could be used to perform the shadow information restoration, were found within a proximate spatial distance to the target shaded pixel; continuum removal is a brightness normalization technique for spectral absorption feature analysis, where a spectrum is separated into two parts: continuum information (CI) represents brightness level, and continuum removed (CR) for spectrum shape information; shadows were recognized at region scale after image segmentation, for each object, an outward buffer having a width of two pixels was generated around a shadow object, and within the buffer area, N spectrally nearest pixels to a target shaded pixel, where the spectral root mean square deviation (RMSD) satisfies the conditions for spectral similarity, were selected as non-shaded similar pixels. Here, RMSD is calculated as: n RMYi (CR (xi, yi, b) -CR (x, y, b))2 RMSDi r n where CR (xi, yi, b) is the CR value of ith non-shaded pixel located in (xi, yr) in band b, CR (x, y, b) shares the same definition but for a target pixel in the shadow region, and n is the number of spectral bands; a large RMSD denotes a large difference in the shape of the spectrum, whereas a small RMSD indicates high similarity between the two pixels; the above mentioned condition for the spectral similarity between target shaded pixel and 12 alternative non-shaded pixel is defined as: RMSD 0[u(b)X 2/m /n where,7(b) is the standard deviation of CR value for the whole input image in band b, and m is the number of land cover classes; moreover, the estimated number of classes (m) needs to be predefined; this value is an empirical parameter and varies with the complexity of the landscape, which can be estimated by visual interpretation of the input images, or using an existing land cover map; in the current study the value of four for m was used; when the number of similar pixels within a width buffer of two pixels does not up to 20, the buffer can be expanded by two pixels repeatedly to search for similar pixels in a larger area until the required number of similar pixels 20 was met; C) shadow information retrieving: restoring the DN value of a target shaded pixel using its own CR and CI provided by the non-shaded similar pixels; hence, the CI of the target shaded pixel can be replaced by the weighted average CI derived from the non-shaded similar pixels to uplift its brightness to a comparable level with the surrounding similar pixels; the following equation represents this process: OIres() = Cwavg(b)X CR(b) where OIre(b) is the restored DN value of the target shaded pixel at a specific wavelength (X), or band (b), CR(b) is the continuum removed value of the target shaded pixel at corresponding wavelength (X), or band (b); CIwavg(b) is the weighted average continuum values of non-shaded similar pixels as calculated from: N Cwavg()= Y W X CI(h) J=1 where CIj(b) is the continuum value of similar pixelj; the weight Wj denotes contribution of similar pixelj to CIwavg(b).

2. The method according to claim 1, wherein in step A, the image segmentation can be performed by eCognition Developer 8.0; shape factor Ssh can be set as 0.1, and compactness factor Scm can be set as 0.5; the scale factor Sc can be decided according to the heterogeneous degree of specific areas.

3. The method according to claim 1, wherein in step A, the automatic thresholding method can be described as: the separability of the two classes aB 2 can be maximized when applying the threshold on two separate classes, 13 (TB 0. 0 T2 1T) a 1 C -PT )2 where po and pi are the mean of the two separated classes, respectively, pT is the total mean of the original data; and wo and wi are the proportion of the total number of pixels with the two classes.

4. The method according to claim 1, wherein in step B, the maximum buffer width was set to 10; if the widest buffer still could not provide at least 20 similar pixels, all selected similar pixels were used, regardless of their actual number.