An Unsupervised Method for Removing Shadows from Remote Sensing Images

Abstract

The present invention discloses an unsupervised remote sensing image shadow removal method, which includes the following steps: using the mean shift method to perform plot segmentation on remote sensing images respectively. A new image texture feature extraction method is proposed. Considering the distance constraint of the first law of geography, the first-order derivatives and second-order derivatives of the R, G, B channels and brightness in the x and y directions, a feature vector of pixel points is constructed, and then a feature matrix of the segmented plots is obtained by combining the covariance formula. Then, the singular value decomposition method is used to obtain the eigenvalues of the region, and the feature distances between plots are calculated. Using the matched illumination area, the method of illumination compensation is used to achieve the overall removal of the umbra. Combining the Manhattan distance and the position of the illumination area, a new dynamic weighted penumbra compensation method is proposed to achieve a natural transition from the illumination area to the shadow area, and finally a shadow-free remote sensing image is obtained. The present invention can effectively remove shadows for remote sensing regions with complex structures.

Description

Technical Field

[0001] The present invention relates to the field of image shadow removal, and particularly to an unsupervised remote sensing image shadow removal method. Background Art

[0002] Due to the existence of many shadow areas in remote sensing images, the ground object features are interfered by the shadows and become unclear. The existence of shadows seriously affects various applications of remote sensing images, such as target detection in the projected area and extraction of digital surface information. For high-spatial-resolution remote sensing images, shadow removal is crucial for target recognition and information restoration. Summary of the Invention

[0003] In order to remove the shadow areas in remote sensing images, the present invention provides an unsupervised remote sensing image shadow removal method, which includes the following steps:

[0004] S1. Obtain a remote sensing image. According to the difference in the ground object sizes between the shadow area and the illumination area in the remote sensing image, use the mean shift algorithm to adaptively segment the illumination area and the shadow area into segmented plots, where the segmented plots include shadow area plots and illumination area plots;

[0005] S2: Use the gradient reciprocal weighted method to extract the texture features of the segmented plots to obtain the texture features of the segmented plots;

[0006] S3: Add distance constraints, fuse the R, G, B color channels of the texture features of the segmented plots and their first-order and second-order derivatives in the x direction of the brightness, construct the feature vectors of the pixel points of the segmented plots, and use the feature vectors of the pixel points and the covariance formula to construct the feature matrix of the segmented plots;

[0007] S4: Use singular value decomposition to calculate the eigenvalues of the feature matrix of the segmented plots, and use the eigenvalues to achieve global search and matching between the shadow area plots and the illumination area plots to obtain the illumination area plots matched with the shadow area plots;

[0008] S5: Use the illumination compensation principle to calculate the illumination optimization factor, and use the illumination optimization factor to optimize the illumination area plots matched with the shadow area plots to remove the umbra of the illumination area plots and obtain the illumination area plots after removing the umbra;

[0009] S6: For the transition part between the illumination area and the shadow area, use the dynamic weighted penumbra compensation method combined with the Manhattan distance to achieve a natural transition from the illumination area to the shadow area, and finally output the shadow removal result.

[0010] The beneficial effects provided by the present invention are as follows: In remote sensing images with complex ground objects and diverse scenes, the original ground object structure within the area can be restored from the shadows. While eliminating the image shadows, the texture details of the shadow areas can be restored, enabling the restored shadow areas to be naturally integrated with the brightness and color of the surrounding non-shadow areas. BRIEF DESCRIPTION OF THE DRAWINGS

[0011] Figure 1 It is a schematic flow diagram of the method of the present invention. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0012] To make the objectives, technical solutions, and advantages of the present invention clearer, the embodiments of the present invention will be further described below in conjunction with the accompanying drawings.

[0013] Please refer to Figure 1 , Figure 1 which is a schematic flow diagram of the method of the present invention; An unsupervised remote sensing image shadow removal method includes the following steps:

[0014] S1. Obtain a remote sensing image. According to the difference in the ground object sizes between the shadow area and the illuminated area in the remote sensing image, use the mean shift algorithm to adaptively segment the land parcels in the illuminated area and the shadow area to obtain segmented land parcels, where the segmented land parcels include shadow area land parcels and illuminated area land parcels;

[0015] It should be noted that the mean shift method is an existing method, and only the corresponding formula principle is described here. The formulas involved in step S1 are as follows:

[0016]

[0017]

[0018]

[0019] h = {h s , h r} (4)

[0020] ms shift = y i,j ′ - y i.j (5)

[0021]

[0022] Equations (1)-(4) are the process of calculating the new clustering center. Among them, in Equation (1), y i.j is the initial clustering center, which is actually a 5-dimensional vector combining the spatial position and the LUV chromaticity space. y i.j = (i, j, L i.j , u i.j , v i.j), and then calculate the new cluster center y according to formula (1) i,j ′. Wherein formula (2) is the derivative of the Gaussian kernel function of formula (3). Formula (4) defines the spatial distance threshold h s and color distance threshold h r Formula (5) is used to calculate the new cluster center y i,j ′ to the initial cluster center y i,j The mean shift vector ms shift Formula (6) is the specific calculation method of the mean drift vector modulus. When the modulus of the drift vector does not meet the convergence condition, these calculated values ​​will be used as the new center point for the next iteration. The described process will be repeated until the convergence condition is met or the maximum number of iterations is reached. Setting appropriate parameters can achieve the segmentation of the illuminated area and the shadow area to obtain the segmented plot.

[0023] The specific algorithm steps are:

[0024] 1. First, a pixel point y of the RGB image i.j After mean shift filtering, the convergence point y is obtained. i.j * .

[0025] 2. Calculate the pixel y i.j * The color distance color_distance to the adjacent pixel, if it is less than the given color distance h r , then these pixels are grouped together.

[0026] 3. Create a linked list to save all blocks, calculate the average color distance between blocks, and use the color threshold h r Merge adjacent blocks of similar color to build a block set {C p}.

[0027] 4. Count the area of ​​the blocks and find the small blocks that do not meet the minimum area M. Use the color distance h r Merge to the adjacent block closest to the small block, and finally get the classification set L i ,i=1,…,p, thus realizing the image segmentation process.

[0028] Combined with the characteristics of remote sensing images, we set appropriate parameters h for the illuminated area and the shadow area respectively. s ,h r And the minimum area M. Thus, a more comprehensive and detailed remote sensing image plot segmentation effect is obtained. In the present invention, the above process is realized by the corresponding function of Matlab.

[0029] S2: Use the inverse gradient weighted method to extract texture features of the segmented plots to obtain the texture features of the segmented plots;

[0030] It should be noted that in a remote sensing image, the gradient change in adjacent regions is greater than that within the region. The gradient value is proportional to the difference in gray levels of neighboring pixels. That is, within the region of the image block, the gradient value changes little, while at the plot boundary region, the gradient change is large. Taking the reciprocal of the difference of neighboring pixels as the weight, the weight within the region is greater than that of the neighboring points near the edge or outside the region. Such processing can better preserve the details of the plot edge without being significantly damaged, and at the same time can achieve a smoothing effect and weaken the influence of image noise in the texture extraction process. The formulas involved in step S2 are as follows:

[0031]

[0032]

[0033]

[0034] In formulas (7)-(9), f(i, j) is a certain pixel point of the input image, and a 3*3 window is defined. G(i+m, j+n) is the reciprocal of the difference between the neighboring point (i+m, j+n) and the central point (i, j). It is stipulated that when |f(i+m, j+n)-f(i, j)| = 0, G(i+m, j+n) = 0. w(i+m, j+n) is the normalized weight, that is, the normalized value of the difference reciprocal G(i+m, j+n). It is stipulated that the normalized weight of the central point pixel is 0.5. f’(i, j) is the weighted sum of the normalized weight within the window and the pixel value at the current position, that is, the new pixel value under the reciprocal weighting idea. Further, in combination with the rotation angle, the pixel point f’(i, j) is re-binary encoded. The formulas involved are as follows:

[0035]

[0036]

[0037] F_final(i,j)=min(F’ θ (i,j)),θ∈[1,8] (12)

[0038] Equation (10) takes into account the eight-neighborhood of the pixel point and considers it in eight directions. In formula (10), 「」 is the floor function. Taking the pixel point corresponding to a certain direction as the initial pixel point, the other pixel points are compared with the central pixel point in the counterclockwise direction, so as to finally obtain a binary coding sequence in a certain direction. Using formula (11), the binary sequence is converted into a decimal value, denoted as the eigenvalue F' of f(i, j) in the θ direction θ(i, j). Through the above method, binary encodings in 8 directions are obtained, and then decimal values in 8 directions are obtained. Formula (12) selects the smallest value to represent the feature value at this position. Through Step 2, the texture details of the image can be effectively extracted.

[0039] S3: Add distance constraints, fuse the R, G, B color channels of the texture features of the segmented plot and the first and second derivatives of its brightness in the x direction, construct the feature vector of the pixel points of the segmented plot, and use the feature vector of the pixel points and the covariance formula to construct the feature matrix of the segmented plot;

[0040] It should be noted that considering the first principle of geography, distance constraints are added. Secondly, combined with the texture feature extraction method in Step S2, fuse the R, G, B color channels and the first and second derivatives of the brightness in the x and y directions to construct the feature vector of the pixel points. Use the texture feature extraction method in Step 2 to construct the 7-dimensional feature vector of the pixel points. Then use the covariance formula to construct the feature matrix of the plot. The formulas involved are:

[0041] z k =[R, G, B, g x , g y , g xx , g yy T (13)

[0042]

[0043] Among them, formula (13) represents that z k is the d-dimensional feature vector of the k-th pixel point. In formula (14), C r is the covariance matrix of region r, which describes the texture features of region r. n is the number of pixel points in region r, and μ is the average feature vector of region r. At this time, C r is a 7x7 covariance matrix. When calculating the covariance matrix, the feature vector of the point subtracts the average feature vector of the region. In addition, the covariance matrix is insensitive to the change of illumination, which can also reduce the influence caused by the different brightness between the shadow block and the non-shadow block. Through Step 2 and Step 3, a corresponding feature matrix can be constructed for each plot in the remote sensing image in Step 1.

[0044] S4: Use singular value decomposition to calculate the eigenvalues of the feature matrix of the segmented plot, and use the eigenvalues to achieve global search and matching between the plots in the shadow area and the plots in the illuminated area, and obtain the plots in the illuminated area that match the plots in the shadow area;

[0045] It should be noted that the formulas involved in Step S4 are as follows:

[0046] A = U∑V​T (15)

[0047]

[0048] Equation (15) represents the singular value decomposition process. Assuming that A is a 7x7 covariance matrix, then U is a 7x7 orthogonal matrix, and each row in the matrix is called a left singular vector. ∑ is a 7x7 diagonal matrix, and the elements on the diagonal are called singular values. V T is the transpose matrix of matrix V, which is a 7x7 orthogonal matrix, and each column in the matrix is called a right singular vector. Equation (16) represents the calculation of the characteristic distance. After performing singular value decomposition on the shaded plot, the eigenvalue matrix obtained is ∑ shadow = diag{α1, α2,..., α7}. Similarly, after performing singular value decomposition on a certain plot in the illuminated area, the eigenvalue matrix ∑ sunlit = diag{β1, β2,..., β7} is obtained, and then the characteristic distance between the illuminated block and the shaded block is calculated. By calculating the characteristic distance between each shaded block and the illuminated block, when the characteristic distance value is the smallest, it can be considered that the shaded block has found its similar illuminated block.

[0049] S5: Calculate the illumination optimization factor using the illumination compensation principle, and use the illumination optimization factor to optimize the illuminated area plot matched with the shaded area plot to achieve the removal of the umbra of the illuminated area plot, and obtain the illuminated area plot after removing the umbra;

[0050] The formula involved in the illumination compensation principle is:

[0051]

[0052]

[0053] In equations (17)-(18), where α x is the attenuation factor of the direct light source at point x, which can be understood as the transparency of the direct light source. shadow x is the RGB value of point x in the shadow area of the image, sunlit′ x is the RGB value in the corresponding illuminated area, removal x represents the RGB value after removing the shadow. Using the illumination compensation principle, and using the similar illuminated blocks matched with the shadow area, the overall illumination compensation for the umbra area is performed.

[0054] S6: For the transition part between the illuminated area and the shadow area, adopt a dynamic weighted penumbra compensation method combined with the Manhattan distance to achieve a natural transition from the illuminated area to the shadow area, and finally output the shadow removal result.

[0055] In step S6, the specific formula of the dynamic weighted penumbra compensation method combining Manhattan distance is as follows:

[0056]

[0057]

[0058]

[0059] Among them, the size of the rectangular window is M*N. Formula (9) combines the Manhattan distance to calculate the weights of adjacent points in the transition part between the illuminated area and the shadow area; formula (10) calculates the normalized weights of the rectangular window; in formula (11), removal is the plot of the illuminated area after removing the umbra. Combining the position of the illuminated area, the position of the pixel point (i, j) in the transition part in the rectangular window is adjusted, and weighted calculation is performed on the transition part in each channel, thereby realizing the illumination restoration of the transition part, and finally realizing the natural transition from the illuminated area to the shadow area to obtain a shadowless image.

[0060] The beneficial effects of the present invention are: in remote sensing images with complex ground objects and diverse scenes, the original ground object structure in the area can be restored from the shadow. While eliminating the image shadow, the texture details of the shadow area can be restored, so that the restored shadow area is naturally integrated with the brightness and color of the surrounding non-shadow areas.

[0061] The above are only the preferred embodiments of the present invention, and are not intended to limit the present invention. Any modifications, equivalent replacements, improvements, etc. made within the spirit and principle of the present invention shall be included in the protection scope of the present invention.

Claims

1. An unsupervised method for removing shadows from remote sensing images, characterized in that: It includes the following steps: S1. Obtain a remote sensing image. According to the difference in the ground object sizes between the shadow area and the illuminated area in the remote sensing image, use the mean shift algorithm to adaptively segment the plot in the illuminated area and the shadow area to obtain segmented plots, where the segmented plots include shadow area plots and illuminated area plots; S2: Use the gradient reciprocal weighted method to extract the texture features of the segmented plots to obtain the texture features of the segmented plots; S3: Add a distance constraint, fuse the R, G, B color channels of the texture features of the segmented plots and their first and second derivatives in the x direction of the brightness, construct the feature vector of the pixel points of the segmented plots, and use the feature vector of the pixel points and the covariance formula to construct the feature matrix of the segmented plots; S4: Use singular value decomposition to calculate the eigenvalues of the feature matrix of the segmented plots, use the eigenvalues to achieve global search and matching between the shadow area plots and the illuminated area plots, and obtain the illuminated area plots matched with the shadow area plots; S5: Calculate the illumination optimization factor according to the illumination compensation principle, use the illumination optimization factor to optimize the illuminated area plots matched with the shadow area plots, and achieve the removal of the umbra of the illuminated area plots to obtain the illuminated area plots after removing the umbra; S6: For the transition part between the illuminated area and the shadow area, use the dynamic weighted penumbra compensation method combined with the Manhattan distance to achieve a natural transition from the illuminated area to the shadow area, and finally output the shadow removal result.

2. The unsupervised remote sensing image shadow removal method according to claim 1, characterized in that: The calculation formula of step S2 based on the gradient reciprocal weighted method is as follows: m,n∈[-1,1] In formulas (1)-(3), f(i,j) is a certain pixel point of the segmented plot, and a 3*3 window is defined; G(i+m,j+n) is the reciprocal of the pixel difference between the pixel point (i+m,j+n) and the central point (i,j); when |f(i+m,j+n)-f(i,j)| = 0, G(i+m,j+n) = 0; w(i+m,j+n) is the normalized value of G(i+m,j+n); f’(i,j) is the weighted sum of the normalized weights within the window and the pixel value at the current position, that is, the new pixel value under the reciprocal weighted method, and all the new pixel values form the texture features of the segmented plot.

3. The unsupervised remote sensing image shadow removal method according to claim 2, characterized in that: After obtaining f’(i,j) in step S2, it is further updated in combination with the rotation angle, specifically as follows: F_final(i,j) = min(F’ θ (i,j)), θ ∈ [1, 8] (6) Formula (4) combines the eight neighborhoods of the pixel point f’(i,j) and is considered in eight directions; in formula (4), 「」 is the floor function; Take the pixel point corresponding to a certain direction as the initial pixel point, and compare the other pixel points with the central pixel point in the counterclockwise direction to finally obtain a binary coding sequence in a certain direction; Convert the binary sequence into a decimal value using formula (5), denoted as the eigenvalue \(F'\) of \(f'(i, j)\) in the \(\theta\) direction θ (i, j); Through the above method, the binary codes in eight directions are obtained, and then the decimal values corresponding to the eight directions are obtained. Finally, the smallest value F_final(i,j) is selected to represent the eigenvalue at this position, as shown in formula (6).

4. An unsupervised remote sensing image shadow removal method according to claim 1, characterized in that: The calculation formulas of the feature vector of the pixel points in step S3 and the feature matrix of the segmented plots are as follows: z k = [R, G, B, g x , g y , g xx , g yy T (7)​ Among them, formula (7) represents z k as the d-dimensional feature vector of the k-th pixel; C in formula (8) r is the feature matrix of region r, which describes the texture feature of region r. n is the number of pixels in region r, and μ is the average feature vector of region r.

5. The unsupervised remote sensing image shadow removal method according to claim 1, wherein: In step S6, the specific formula of the dynamic weighted penumbra compensation method combined with the Manhattan distance is as follows: Among them, the size of the rectangular window is M*N. Formula (9) is the weight of adjacent points in the transition part between the illuminated area and the shadow area calculated by combining the Manhattan distance; total_number = M*N; Formula (10) is to calculate the normalized weight of the rectangular window; in Formula (11), removal is the plot of the illuminated area after removing the umbra; combining the position of the illuminated area, adjust the position of the pixel point (i,j) in the transition part in the rectangular window, and perform weighted calculation on the transition part in each channel, so as to realize the illumination recovery of the transition part, and finally realize the natural transition from the illuminated area to the shadow area, and obtain a shadow-free image.

Images