Multi-method fusion SAR (Synthetic Aperture Radar) image shadow detection method and system

Abstract

The invention discloses a multi-method fusion SAR image shadow detection method and system, and the method comprises the steps: processing an original image set according to a wavelet algorithm, a Markov random field and a region growing algorithm, and obtaining a first-scene cut image set and a second-scene cut image set; processing the first-scene cut image set and the second-scene cut image set according to an amplitude detection method and a coherence coefficient detection method to obtain a first shadow region set, a first shadow undetermined region set, a second shadow region set and a second shadow undetermined region set; according to the four region sets, performing set operation to obtain a first partial shadow region set and a third shadow undetermined region set; processing the third shadow undetermined region set according to a radar sight line and a digital elevation model analysis algorithm to obtain a second partial shadow region set; and according to the first partial shadow region set and the second partial shadow region set, a final detection result is obtained through union. According to the method, the accuracy and the reliability of shadow detection are improved by fusing multiple algorithms.

Description

Technical Field

[0001] The present invention belongs to the technical field of computer vision, and in particular relates to a multi-method fusion SAR image shadow detection method and system. Background Art

[0002] Spaceborne synthetic aperture radar interferometry (InSAR) is an advanced space geodetic method that obtains information such as small deformation and elevation of the surface by calculating the phase information between multi-view coherent radar images. However, since synthetic aperture radar (SAR) is an active radiation side-view imaging technology, it is easy to produce shadows and overlap in areas with large terrain undulations. Shadows can lead to the loss of information in SAR images, and also affect the recognition of ground objects and phase unwrapping in InSAR data processing. Therefore, accurately detecting and marking shadow areas in SAR images is crucial to ensure data quality and the accuracy of subsequent processing.

[0003] However, existing shadow detection is often performed together with overlapping detection. The complete single shadow detection technology has not received enough attention in previous studies, and each shadow detection algorithm has its own shortcomings. Summary of the invention

[0004] In view of the shortcomings of the existing technology, the present application provides a multi-method fusion SAR image shadow detection method and system, which integrates multiple algorithms, solves the problem of accurate detection of shadow areas in SAR images, and improves the accuracy and reliability of shadow detection.

[0005] To achieve the above objectives, in a first aspect, the present invention provides a multi-method fusion SAR image shadow detection method, comprising:

[0006] According to the preset wavelet multi-scale bilateral filtering algorithm, the original two-view image set is processed to obtain the two-view image set;

[0007] According to a preset Markov random field and region growing algorithm, the two-view image sets are processed to obtain a one-view cut image set and a two-view cut image set;

[0008] According to a preset amplitude adaptive threshold detection method, the one-scene cut image set is processed to obtain a first shadow area set and a first shadow undetermined area set;

[0009] According to a preset coherence coefficient detection method, the one-scene cut image set and the two-scene cut image set are processed to obtain a second shadow area set and a second shadow undetermined area set;

[0010] According to the first shadow area set, the first shadow undetermined area set, the second shadow area set and the second shadow undetermined area set, a set operation is performed to obtain a first partial shadow area set and a third shadow undetermined area set;

[0011] According to a preset radar line of sight and digital elevation model analysis algorithm, the third shadow undetermined area set is processed to obtain a second partial shadow area set;

[0012] The final shadow area detection result is obtained by calculating the union of the first partial shadow area set and the second partial shadow area set.

[0013] Compared with the prior art, the embodiments of the present application have the following beneficial effects: by processing the original image through the wavelet transform decomposition algorithm, the noise can be effectively removed while retaining important features such as edges, and the bilateral filtering process can smooth the image while maintaining the edges, thereby improving the quality of the image; a Markov random field model is constructed to make full use of the spatial dependence between pixels to improve the accuracy of segmentation, and a region growing algorithm is used to more accurately segment the target area by selecting seed points and gradually expanding the area; an amplitude adaptive threshold detection method is adopted to adapt to the image characteristics in different scenes and improve the robustness of detection; coherence calculation can effectively identify shadow areas caused by physical characteristics, especially those caused by terrain or object occlusion; the results obtained by multiple methods are set and operated, which can combine the advantages of various methods, reduce false detection and missed detection, and improve the overall detection effect; through radar line of sight and digital elevation model analysis, the shadow area caused by insufficient radar depression angle is further detected, and the detection accuracy in areas with large terrain undulations is improved; finally, the SAR image shadow detection result fused by multiple methods is output, which significantly improves the accuracy and comprehensiveness of shadow detection.

[0014] In some embodiments of the first aspect of the present application, the processing of the original two-view image set according to a preset wavelet multi-scale bilateral filtering algorithm to obtain the two-view image set includes:

[0015] According to a preset wavelet transform decomposition algorithm, the original two-view image in the original two-view image set is processed respectively to obtain each first high-frequency sub-band and first low-frequency sub-band;

[0016] According to a preset bilateral filtering algorithm, the first high frequency sub-bands are processed respectively to obtain second high frequency sub-bands;

[0017] According to a preset smoothing algorithm, the first low-frequency sub-bands are processed respectively to obtain second low-frequency sub-bands;

[0018] According to a preset wavelet inverse transform reconstruction algorithm, each of the second high-frequency sub-bands and the second low-frequency sub-bands is processed to obtain each of the two-scene images, that is, to obtain a two-scene image set.

[0019] Compared with the prior art, the above embodiment has the following beneficial effects: by processing the original image through the wavelet transform decomposition algorithm, noise can be effectively removed while retaining important features such as edges; bilateral filtering can smooth the image while maintaining the edges, thereby improving the image quality.

[0020] In some embodiments of the first aspect of the present application, the processing of the two-view image sets according to a preset Markov random field and region growing algorithm to obtain a one-view cut image set and a two-view cut image set includes:

[0021] According to the normalization algorithm, the one-view image set and the two-view image set are processed respectively to obtain the one-view normalized image set and the two-view normalized image set;

[0022] According to the one-scene normalized image set and the two-scene normalized image set, a Markov random field model is constructed respectively; wherein the Markov random field model includes a model of a label field;

[0023] According to a preset optimization algorithm, the Markov random field models are processed respectively to obtain a one-scene preliminary segmentation image set and a two-scene preliminary segmentation image set;

[0024] According to a preset automatic region growing algorithm, the first scene preliminary segmentation image set and the second scene preliminary segmentation image set are processed respectively to obtain the first scene secondary segmentation image set and the second scene secondary segmentation image set;

[0025] According to a preset morphological algorithm, the one-scene secondary segmentation image set and the two-scene secondary segmentation image set are processed respectively to obtain a one-scene cutting image set and a two-scene cutting image set.

[0026] Compared with the prior art, the above embodiments have the following beneficial effects: the Markov random field model can well capture the spatial correlation between pixels in the image, perform global optimization, effectively reduce noise interference, and provide fine preliminary segmentation results; region growing can further refine the segmented area, especially in the edge or transition area, to ensure accurate division of shadow and non-shadow areas; morphological processing can eliminate isolated noise points in the segmentation results, smooth the boundaries, and improve the segmentation accuracy and integrity.

[0027] In some embodiments of the first aspect of the present application, a Markov random field model is constructed according to the one-scene normalized image set and the two-scene normalized image set, respectively; wherein the Markov random field model includes a model of a label field, including:

[0028] A Markov random field is defined, each pixel of each image in the one-view normalized image set and the two-view normalized image set is regarded as a random variable, the grayscale distribution in the image is modeled according to a finite Gaussian mixture model, and an original Markov random field model is constructed;

[0029] Optimizing the parameters of the original Markov random field model using an expectation-maximization algorithm;

[0030] The Markov random field model is obtained by using Gibbs distribution as the prior probability distribution of the marking field and using the Porter model as the model of the marking field.

[0031] Compared with the prior art, the above embodiment has the following beneficial effects: defining a Markov random field model, representing the spatial relationship between pixels by establishing a neighborhood system, and constructing a probability distribution using the Gibbs distribution, which is beneficial for capturing the spatial correlation between pixels in the image and facilitating subsequent global optimization.

[0032] In some embodiments of the first aspect of the present application, the Markov random field models are processed separately according to a preset optimization algorithm to obtain a first-scene preliminary segmented image set and a second-scene preliminary segmented image set, including:

[0033] Adaptive weights are introduced, where the weight estimation model is as follows:

[0034]

[0036] Where α(t) is the energy weight of the feature field, C is a constant, t is the number of iterations, L is the number of categories, and x i and i represents the i-th pixel of the first scene image and the i-th pixel of the second scene image respectively, the subscript i represents the pixel position of the image, U(x i ) represents pixel x i The corresponding label field energy function, P represents the probability function, N i Represents the set of neighboring pixels of pixel i;

[0037] Introducing penalty function The definition is as follows:

[0038]

[0039] The label field energy function U(x i )for: Where β is a custom weighting coefficient ranging from 0 to 1, V2(x i ,x l ) represents the weight of the edge corresponding to the pixels, and exp(*) represents the exponential operation;

[0040] The final energy function expression of the posterior probability solved pixel by pixel is:

[0041]

[0042] According to the preset number of classifications and the K-means algorithm, the images in the one-scene normalized image set and the two-scene normalized image set are respectively initially segmented to obtain an initial label for each pixel;

[0043] According to the preset number of iterations, each of the Markov random field models is iteratively optimized to obtain a set of preliminary segmented images of one scene and a set of preliminary segmented images of two scenes, and in each iteration, according to the preset expectation algorithm, the model parameters of the Markov random field model are estimated and calculated, and according to the preset iterative algorithm, the maximum likelihood parameter estimate of each category model is calculated until the preset convergence standard is met, and the characteristic field energy weight under the current number of iterations and the model parameters of each category calculated by the expectation maximization algorithm are calculated, and substituted into the final energy function expression to update the segmentation result of each pixel.

[0044] Compared with the prior art, the above embodiment has the following beneficial effects: the Bode model, as a model of the marking field, achieves global optimal segmentation through iterative optimization, effectively reduces noise interference, and provides a refined preliminary segmentation result.

[0045] In some embodiments of the first aspect of the present application, the first scene preliminary segmentation image set and the second scene preliminary segmentation image set are processed respectively according to a preset automatic region growing algorithm to obtain the first scene secondary segmentation image set and the second scene secondary segmentation image set, including:

[0046] According to a preset rule, a number of seed points are selected from the images in each of the one-scene preliminary segmentation image set and the two-scene preliminary segmentation image set;

[0047] According to the preset initial window size, the window size is gradually increased, and the texture difference of each pixel under each window size is calculated in turn, and the window size with the smallest texture difference is taken as the optimal window size;

[0048] According to a preset initial window, the probability distribution of each pixel belonging to each texture class is calculated, and according to the probability distribution, the entropy of the initial window is calculated;

[0049] According to the optimal threshold and the optimal window size, the texture difference between each pixel point in the optimal window where various sub-points are located and the corresponding seed point is calculated, and the pixel points whose texture difference meets the optimal threshold are screened and merged into the initial area of ​​the corresponding seed point until all pixel points are screened, thereby obtaining a one-scene secondary segmentation image set and a two-scene secondary segmentation image set respectively; wherein the texture difference is obtained by calculating the difference in grayscale co-occurrence values ​​of each pixel point in the window, and the optimal threshold is calculated by the standard deviation of each initial window class.

[0050] Compared with the prior art, the above embodiment has the following beneficial effects: using the automatic region growing method, by selecting seed points and gradually expanding the region, the target area can be segmented more accurately, especially in the edge or transition area, ensuring accurate division of shadow and non-shadow areas.

[0051] In some embodiments of the first aspect of the present application, the step of processing the scene cut image set according to a preset amplitude adaptive threshold detection method to obtain a first shadow area set and a first shadow undetermined area set includes:

[0052] Obtaining the original amplitude image corresponding to each image in the scene cut image set;

[0053] Compressing the original amplitude image according to a grayscale compression algorithm to obtain an amplitude image;

[0054] According to a preset initial threshold, the shadow area and the non-shadow area of ​​the amplitude image are repeatedly divided, and after each division, the mean and standard deviation of the shadow area and the non-shadow area are calculated, and the initial threshold is updated according to the mean and standard deviation until the updated change of the initial threshold meets the preset difference range, thereby obtaining an optimal threshold;

[0055] The amplitude image is processed according to the optimal threshold to obtain a first shadow area set and a first shadow undetermined area set.

[0056] Compared with the prior art, the above embodiment has the following beneficial effects: by calculating the mean and variance of the shadow area and the background area, using an iterative algorithm to optimize and obtain the optimal threshold, and using the optimal threshold to separate the shadow area and the non-shadow area, when the difference in grayscale values ​​between the shadow and non-shadow areas is small, the false detection rate can be effectively reduced and the accuracy of shadow detection can be improved.

[0057] In some embodiments of the first aspect of the present application, the step of processing the first scene cut image set and the second scene cut image set according to a preset coherence coefficient detection method to obtain a second shadow area set and a second shadow undetermined area set includes:

[0058] According to the one-scene cut image set and the two-scene cut image set, the coherence value of each two-scene image after the phase correction factor is introduced is calculated, and the calculation method is as follows:

[0059] Where z1, z2 represent the complex values ​​of the two images; * represents the conjugate operator; i represents the pixel number; represents the phase correction factor, Indicates the phase caused by subtracting the terrain factor from the original phase δ flattenedIt represents the coherence value after introducing the phase correction factor, LL is the number of pixels in the window, and j is the imaginary unit;

[0060] According to the preset spatial averaging algorithm, process the coherence value after introducing the phase correction factor to obtain the final coherence value;

[0061] According to the final coherence value, process the first scene of segmented image set and the second scene of segmented image set to obtain the second shadow area set and the second shadow pending area set.

[0062] Compared with the prior art, the above embodiments have the following beneficial effects: calculating the coherence coefficient of two scenes of images, using the phase correction factor and the spatial neighborhood average value to improve the unbiasedness of coherence estimation; setting a coherence threshold to extract the shadow area, which can more accurately identify the shadow area with low coherence and reduce the deviation caused by terrain or noise.

[0063] In some embodiments of the first aspect of the present application, the processing of the third shadow pending area set according to the preset radar line of sight and digital elevation model analysis algorithm to obtain the second partial shadow area set includes:

[0064] The preset radar line of sight and digital elevation model analysis algorithm includes a geometric relationship formula:

[0065] H / M≥ksinβ1; 0<k≤1; where M, H, and β1 are the slant range, elevation difference, and radar depression angle corresponding to two pixel points respectively;

[0066] According to the geometric relationship, process the images of the third shadow pending area set, and classify the pixel points that meet the geometric relationship formula as the shadow area to obtain the second partial shadow area set.

[0067] Compared with the prior art, the above embodiments have the following beneficial effects: combining the radar line of sight and digital elevation model data, analyzing the geometric relationship of the ground objects in the image, judging the shadow area by calculating the relationship between the ground height difference and the slant range, further detecting the shadow area caused by insufficient radar depression angle, and improving the detection accuracy in areas with large terrain undulations.

[0068] In the second aspect, the embodiments of the present application also provide a SAR image shadow detection system integrating multiple methods, including: a preprocessing module, a two-step segmentation module, an amplitude detection module, a coherence detection module, a set operation module, a radar elevation detection module, and a result output module;

[0069] Among them, the preprocessing module is used to process the original two scenes of image sets according to the preset wavelet multi-scale bilateral filtering algorithm to obtain two scenes of image sets;

[0070] The two-step segmentation module is used to process the two-view image set according to a preset Markov random field and region growing algorithm to obtain a one-view cut image set and a two-view cut image set;

[0071] The amplitude detection module is used to process the scene cut image set according to a preset amplitude adaptive threshold detection method to obtain a first shadow area set and a first shadow undetermined area set;

[0072] The optical coherence detection module is used to process the one-scene cut image set and the two-scene cut image set according to a preset coherence coefficient detection method to obtain a second shadow area set and a second shadow undetermined area set;

[0073] The set operation module is used for performing set operation to obtain a first partial shadow area set and a third shadow area set according to the first shadow area set, the first shadow area set to be determined, the second shadow area set and the second shadow area set to be determined;

[0074] The radar elevation detection module is used to process the third shadow undetermined area set according to a preset radar line of sight and digital elevation model analysis algorithm to obtain a second partial shadow area set;

[0075] The result output module is used to calculate the final shadow area detection result based on the union of the first partial shadow area set and the second partial shadow area set.

[0076] Compared with the prior art, the embodiments of the present application have the following beneficial effects: by processing the original image through the wavelet transform decomposition algorithm, the noise can be effectively removed while retaining important features such as edges, and the bilateral filtering process can smooth the image while maintaining the edges, thereby improving the quality of the image; a Markov random field model is constructed to make full use of the spatial dependence between pixels to improve the accuracy of segmentation, and a region growing algorithm is used to more accurately segment the target area by selecting seed points and gradually expanding the area; an amplitude adaptive threshold detection method is adopted to adapt to the image characteristics in different scenes and improve the robustness of detection; coherence calculation can effectively identify shadow areas caused by physical characteristics, especially those caused by terrain or object occlusion; the results obtained by multiple methods are set and operated, which can combine the advantages of various methods, reduce false detection and missed detection, and improve the overall detection effect; through radar line of sight and digital elevation model analysis, the shadow area caused by insufficient radar depression angle is further detected, and the detection accuracy in areas with large terrain undulations is improved; finally, the SAR image shadow detection result fused by multiple methods is output, which significantly improves the accuracy and comprehensiveness of shadow detection. BRIEF DESCRIPTION OF THE DRAWINGS

[0077] Figure 1: A flowchart of a multi-method fusion SAR image shadow detection method provided in some embodiments of the present invention.

[0078] Figure 2 : A schematic diagram of the structure of a multi-method fusion SAR image shadow detection system provided in some embodiments of the present invention.

[0079] Figure 3 : A schematic flow chart of a wavelet multi-scale bilateral filtering algorithm provided in some embodiments of the present invention.

[0080] Figure 4 : A flow chart of a Markov random field and region growing algorithm provided in some embodiments of the present invention.

[0081] Figure 5 : A schematic diagram of the geometric relationship between a radar line of sight and a digital elevation model analysis algorithm provided in some embodiments of the present invention. DETAILED DESCRIPTION

[0082] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0083] Embodiment 1:

[0084] Please refer to Figure 1 , a multi-method fusion SAR image shadow detection method provided by an embodiment of the present invention, comprising steps S1 to S7:

[0085] Step S1: According to a preset wavelet multi-scale bilateral filtering algorithm, the original two-view image set is processed to obtain a two-view image set.

[0086] In some embodiments of the present application, Figure 3 The flowchart of a wavelet multi-scale bilateral filtering algorithm is shown in FIG. 1 , step S1 can be implemented by the following preferred implementation:

[0087] Step S11: According to a preset wavelet transform decomposition algorithm, the original two-view image in the original two-view image set is processed respectively to obtain each first high-frequency sub-band and first low-frequency sub-band.

[0088] In this step, the wavelet function of the wavelet transform decomposition algorithm is The following conditions should be met:

[0089] in, is the basis function or mother wavelet, which can form a cluster of function systems through scaling and translation on the time axis: in, It is called a sub-wavelet, a is the scale factor, which reflects the period length of the wavelet, b is the translation factor, which reflects the time translation, R represents a real number, L 2 (R) represents real two-dimensional space.

[0090] In this step, the wavelet transform function may be selected from Haar wavelet, Dobecy wavelet, Mexican hat wavelet, Mollet wavelet and Mayer wavelet, or other similar wavelet functions may be selected for implementation, and no further limitation is made here.

[0091] Preferably, the calculation formula of wavelet transform can be as follows:

[0092] Among them, w f (a,b) represents continuous wavelet transform, f(t) represents the function to be transformed by wavelet, a represents the scale factor, b represents the time translation factor, Represents a wavelet transform function that has undergone time translation and scale expansion, also known as a sub-wavelet.

[0093] Step S12: processing the first high frequency sub-bands respectively according to a preset bilateral filtering algorithm to obtain the second high frequency sub-bands.

[0094] Bilateral filtering can preserve the edge of the image and smooth the noise at the same time. The calculation formula of bilateral filtering can be as follows:

[0095] Where p and q represent the coordinates of two pixels, namely p = (x0, y0), q = (x1, y1), BF[I] P represents the bilateral filtering centered at pixel p, I p ,I q Respectively represent the pixel values ​​at positions p and q, that is, the intensity values ​​of the pixels, represents Gaussian weighting in the spatial domain, Represents Gaussian weighting within the pixel value domain, W p is a normalization factor used to ensure that the weighted sum of the pixels is 1, W p The expression is as follows: Q represents the set of all pixels in a local neighborhood centered on pixel p.

[0096] Step S13: processing the first low-frequency sub-bands respectively according to a preset smoothing algorithm to obtain second low-frequency sub-bands.

[0097] Among them, the smoothing process can remove noise while retaining the main structure of the image. The smoothing algorithm can be mean filtering: Among them, M represents a certain pixel, Num represents the total number of pixels, and f(i,j) represents the pixel value.

[0098] Preferably, when designing a linear smoothing filter, it should be ensured that the filter has a single peak, ie, a main lobe, and is symmetrical in both vertical and horizontal directions. For example, a typical smoothing filter has a window of 3×3.

[0099] Preferably, the smoothing algorithm may also use a Gaussian filter to remove noise that obeys a normal distribution; wherein a first-order zero-mean Gaussian function may be expressed as: The two-dimensional discrete zero-mean Gaussian function can be expressed as: The distribution exponent σ represents the width of the Gaussian distribution.

[0100] Step S14: processing the second high-frequency sub-band and the second low-frequency sub-band according to a preset wavelet inverse transform reconstruction algorithm to obtain two images, that is, to obtain a two-image set.

[0101] In step S1 of this embodiment, the advantages of wavelet transform and bilateral filtering are combined, so that image details can be retained while removing noise, and it is suitable for application scenarios requiring high-quality SAR images.

[0102] Step S2: processing the two-view image sets according to a preset Markov random field and region growing algorithm to obtain a one-view cut image set and a two-view cut image set.

[0103] In some embodiments of the present application, Figure 4 The schematic diagram of a flow chart of a Markov random field and region growing algorithm shown in FIG. 1 , step S2 can be implemented by the following preferred implementation manner:

[0104] Step S21: According to a normalization algorithm, the one-view image set and the two-view image set are processed respectively to obtain a one-view normalized image set and a two-view normalized image set.

[0105] Preferably, the normalization algorithm may be as follows:

[0106] Among them, I represents the pixel value in the original image; Min and Max represent the maximum and minimum values ​​of the pixels in the original image, and I N Indicates the normalized value.

[0107] Step S22: constructing Markov random field models respectively according to the one-view normalized image set and the two-view normalized image set; wherein the Markov random field model includes a model of a label field.

[0108] Furthermore, the construction of the Markov random field model can be achieved by the following preferred implementations:

[0109] Define a Markov random field, regard each pixel of each image in the one-view normalized image set and the two-view normalized image set as a random variable, model the grayscale distribution in the image according to the finite Gaussian mixture model, and construct an original Markov random field model, wherein the probability distribution of each pixel in the image can be expressed as:

[0110]

[0111] Where S represents the set of random fields, i represents the subset of random fields; P(Y|X=x) represents the probability distribution function or likelihood function of the feature field when the label field is known; y i represents the features in the random field subset i; k = 1, 2, 4, ... L represents the classification label; w represents the dimension of the observed feature y; μ k and∑ k is the expectation and variance of the k-th feature, π k It represents the proportion of pixels with segmentation mark k, and T represents vector transpose.

[0112] Preferably, the expectation maximization algorithm can be used to optimize the parameters of the original Markov random field model, which is achieved by the following preferred implementation:

[0113] The energy function corresponding to the likelihood function of a pixel in the image is expressed as follows:

[0114]

[0115] The parameters θ of the likelihood function are estimated according to the maximum pseudo-likelihood method based on expected maximization, iterating in the following way:

[0116]

[0117]

[0118] in, Indicates the proportion of pixels with k-type features in this iteration; θ k represents the parameters of the likelihood function, θ k ={μ k ,Σ k}; P(k|y i ,θ k ) represents the conditional probability of each pixel; Represents the expectation of k-class features in this iteration; Represents the expectation of k-class features in this iteration; Represents the variance of the feature marked as k in this iteration.

[0119] Preferably, the Gibbs distribution can be used as the prior probability distribution of the labeling field, which is expressed as follows:

[0120] P(x)=Z -1 e -U(x) ; Z = ∑ x∈Ω e -U(x) ; U(x)=∑ c∈C V c (x); where Z is the normalized value of the partitioning function, and the energy function U(x) is a series of potential functions V defined in the group c. c (x); the complement C represents the total set of groups; P(x) is the probability of a certain configuration x appearing in the random field X, and Ω is the set of all possible state vectors.

[0121] Among them, the energy function U(x) can be calculated by the following sum:

[0122]

[0123] Among them, V1, V2, V3 represent potential functions of different orders, and the order represents the complexity of the Gibbs distribution model; C1, C2, C3 represent the group sets corresponding to different orders; x i ,x j ,x k Represents the random variable in the pixel or group at position i, j, k.

[0124] Preferably, when only the second-order neighborhood system is considered, the energy function U(x) is expressed as:

[0125] Preferably, the Bode model can be used as the model of the marking field, and finally the Markov random field model is obtained, wherein the potential function is defined as follows:

[0126]

[0127] The local probability of the Porter model is:

[0128] Where n i (x i ) is the number of cells in the neighborhood of position i that is equal to x i The number of pixel neighborhood positions, β ′ is the coupling coefficient.

[0129] Step S23: According to a preset optimization algorithm, the Markov random field models are processed respectively to obtain a one-scene preliminary segmented image set and a two-scene preliminary segmented image set.

[0130] The preferred algorithm can be implemented as follows:

[0131] Adaptive weights are introduced, where the weight estimation model can be as follows:

[0132]

[0133] Where α(t) is the energy weight of the feature field, C is a constant, t is the number of iterations, L is the number of categories, and x i and i represents the i-th pixel of the first scene image and the i-th pixel of the second scene image respectively, the subscript i represents the pixel position of the image, U(x i ) represents pixel x i The corresponding label field energy function, P represents the probability function, N i Represents the set of neighboring pixels of pixel i.

[0134] Introducing penalty function The definition is as follows:

[0135]

[0136] The label field energy function U(x i )for: Where β is a custom weighting coefficient ranging from 0 to 1, V2(x i ,x l ) represents the weight of the edge corresponding to the pixels, and exp(*) represents the exponential operation;

[0137] The final energy function expression of the posterior probability solved pixel by pixel is:

[0138]

[0139] According to the preset number of classifications and the K-means algorithm, the images in the one-scene normalized image set and the two-scene normalized image set are respectively initially segmented to obtain an initial label for each pixel;

[0140] According to a preset number of iterations, each of the Markov random field models is iteratively optimized to obtain a set of preliminary segmented images of one scene and a set of preliminary segmented images of two scenes, and in each iteration, according to a preset expectation algorithm, each model parameter of the Markov random field model is estimated and calculated, and according to a preset iterative algorithm, the maximum likelihood parameter estimate of each category model is calculated until a preset convergence criterion is met, and the characteristic field energy weight under the current number of iterations and the model parameters of each category calculated by the expectation maximization algorithm are calculated, and substituted into the final energy function expression to update the segmentation result of each pixel; wherein the convergence criterion may be that the absolute value of the parameter difference before and after the update is sufficiently small.

[0141] Step S24: According to a preset automatic region growing algorithm, the first scene preliminary segmentation image set and the second scene preliminary segmentation image set are processed respectively to obtain the first scene secondary segmentation image set and the second scene secondary segmentation image set.

[0142] The automatic region growing algorithm can be implemented in the following preferred manner:

[0143] According to a preset rule, a number of seed points are selected from the images in each of the one-scene preliminary segmentation image set and the two-scene preliminary segmentation image set;

[0144] According to the preset initial window size, the window size is gradually increased, and the texture difference of each pixel under each window size is calculated in turn, and the window size with the smallest texture difference is taken as the optimal window size; wherein, the initial window size may be 3×3, and the texture difference may be obtained by calculating the joint entropy, and the calculation method is as follows:

[0145]

[0146] in, represents the joint co-occurrence matrix of adjacent pixels in the same window, Represents the window texture corresponding to adjacent pixels p and q; Represents the gray level co-occurrence value; N c is the number of texture types, E n represents the joint entropy.

[0147] If the textures of the two windows are equally rough or smooth, the joint matrix will have only a single pixel, that is, the value of the joint entropy is zero.

[0148] According to the preset initial window, the probability distribution of each pixel belonging to each texture class is calculated, and the entropy of the initial window is calculated according to the probability distribution, wherein the calculation method is as follows:

[0149] Among them, α represents the set of texture classes, n represents the size of the window group centered on each pixel, and ε represents the entropy value.

[0150] Furthermore, we define the distance between two entropies as

[0151] in, represents the standard deviation of each initial window class; this entropy distance expression measures the average similarity of two texture objects and quantifies the effectiveness of the segmentation. Once the best window is established, the average value of each initial window class is calculated. and standard deviation Therefore, if the entropy value of a pixel (εp ) is less different from its average value, then the pixel is aggregated into one region.

[0152] According to the optimal threshold and the optimal window size, the texture difference between each pixel point of the optimal window where various sub-points are located and the corresponding seed point is calculated, and the pixel points whose texture difference meets the optimal threshold are screened and merged into the initial area of ​​the corresponding seed point until all the pixel points are screened, and a one-scene secondary segmentation image set and a two-scene secondary segmentation image set are obtained respectively; wherein, the texture difference is obtained by calculating the difference in the grayscale co-occurrence value of each pixel point in the window; the optimal threshold is calculated by the standard deviation of each initial window class, which can be obtained according to the following preferred implementation mode:

[0153] The optimal threshold is calculated as follows: Among them, ∈ α represents the optimal threshold of texture type a, ε′ p represents the joint entropy value of pixel p, μ′ n represents the entropy mean of the initial window of this class, p→Classτ means aggregating pixels p into τ classes, Represents the standard deviation of the texture of type A for each initial window.

[0154] Step S25: According to a preset morphological algorithm, the one-scene secondary segmentation image set and the two-scene secondary segmentation image set are processed respectively to obtain a one-scene cutting image set and a two-scene cutting image set.

[0155] Morphological operations are defined by sets. Generally, two types of pixel sets are used in image processing: target elements and structural elements. The target element is defined as a set of foreground pixels, and the structural element can be determined by foreground pixels and background pixels. Basic morphological operations include erosion, dilation, opening, closing, etc. In this solution, the opening and closing operations are used to process the results, which can remove speckle noise and fill holes.

[0156] For example, the morphological algorithm can be completed by opening and closing operations:

[0157] Among them, the opening operation of the structure element B on the set A (expressed as ) is defined as: It means that B first corrodes A, and then B expands the corrosion result;

[0158] Geometrically, the opening operation of B on A is the union of all translations of B. When the structural element B cannot fit the narrow line segments or figures in A, these small line segments or figures will be deleted. Therefore, the characteristic of the opening operation is to delete areas narrower than the structural element, which can usually smooth the contours of the object, disconnect narrow gaps, and eliminate slender protrusions.

[0159] The closing operation of the structure element B on the set A (expressed as A·B) is defined as: It means that B first expands A, and then B corrodes the expansion result. represents the dilation operation, Represents the erosion operation.

[0160] The closing operation is the complement of the union of all translations of B that do not overlap with A. This operation can bridge narrow breaks and elongated gullies, fill gaps in contours, smooth object edges, and fill small holes inside objects.

[0161] In step S2 of this embodiment, by constructing a Markov random field model, the spatial correlation between pixels in the image can be well captured, global optimization can be performed, noise interference can be effectively reduced, and a fine preliminary segmentation result can be provided; region growing can further refine the segmented area, especially in the edge or transition area, to ensure accurate division of shadow and non-shadow areas; morphological processing can eliminate isolated noise points in the segmentation results, smooth the boundaries, and improve segmentation accuracy and integrity.

[0162] Step S3: According to a preset amplitude adaptive threshold detection method, the scene cut image set is processed to obtain a first shadow area set and a first shadow undetermined area set.

[0163] In some embodiments of the present application, step S3 can be implemented by the following preferred implementations:

[0164] S31: Acquire the original amplitude image corresponding to each image in the scene cut image set;

[0165] S32: compressing the original amplitude image according to a grayscale compression algorithm to obtain an amplitude image;

[0166] S33: According to a preset initial threshold, the shadow area and the non-shadow area of ​​the amplitude image are repeatedly divided, and after each division, the mean and standard deviation of the shadow area and the non-shadow area are calculated, and the initial threshold is updated according to the mean and standard deviation, until the updated change of the initial threshold meets the preset difference range, and the optimal threshold is obtained; wherein the preset difference range is, for example, 10 -6 .

[0167] In this preferred implementation step, an improved threshold segmentation method is used. It is assumed that the amplitude image contains two grayscale areas, z represents the grayscale value in the image, p(z) represents the probability density function corresponding to the histogram, and the mixed probability density function can be expressed as: p(z)=P1p1(z)+P2p2(z); wherein P1 represents the probability that the pixel belongs to the shadow area, P2 represents the probability that the pixel belongs to the non-shadow area, and P1+P2=1.

[0168] Furthermore, assuming that there is a threshold Γ, which minimizes the average error in determining pixels as shadow areas or non-shadow areas, then the probability of mistakenly determining a pixel point belonging to a non-shadow area as a pixel point belonging to a shadow area is: The probability of judging a pixel belonging to the shadow area as a pixel in the non-shadow area is: Then the overall error probability is: E(Γ)=P2E1(Γ)+P1E2(Γ).

[0169] Furthermore, in order to find the threshold that minimizes E(Γ), it is necessary to find the first-order derivative of E(Γ) with respect to Γ and make the first-order derivative equal to 0, and we can obtain: P1p1(Γ)=P2p2(Γ).

[0170] Substituting the Gaussian probability density function into the mixed probability density function, we can obtain:

[0171] Among them, μ1 and represents the mean and variance of the Gaussian probability density function corresponding to the pixels in the shadow area; μ2 and Represents the mean and variance of the Gaussian probability density function corresponding to the pixels in the non-shadow area.

[0172] Furthermore, the mean and standard deviation of the shaded and non-shaded areas can be calculated as follows:

[0173]

[0174] Among them, μ S , represents the mean and standard deviation of the shaded area, ∑ (i，j) f S (i, j) represents the sum of the pixel values ​​in the shadow area, Num S Represents the total number of pixels in the shadow area, μ B , represents the mean and standard deviation of the non-shaded area, ∑ (i，j) f B (i, j) represents the sum of the pixel values ​​in the non-shadow area, Num B Represents the total number of pixels in the non-shadow area, ∑ (i，j) (f S (i, j)-μ S ) 2 represents the sum of squares of the difference between the pixels in the shadow area and the mean, ∑ (i，j) (f B (i, j)-μ B ) 2 Represents the sum of squares of the differences between the pixels in the non-shadow area and the mean.

[0175] According to P1p1(Γ)=P2p2(Γ), we can get: AΓ 2 +BΓ+C=0; where

[0176] Furthermore, when the variances are the same, We can get: When P1=P2, the threshold Γ is the average of the means of the two distributions.

[0177] When the above improved algorithm is applied to a small area in an image, the threshold Γ corresponding to a pixel can be obtained by calculating the grayscale average of the pixels in a window centered on the pixel, as shown below:

[0178] Among them, z(i, j) represents the pixel value at (i, j), and m and n represent the size of the window.

[0179] S34: Processing the amplitude image according to the optimal threshold to obtain a first shadow area set and a first shadow undetermined area set.

[0180] In step S3 of this embodiment, by calculating the mean and variance of the shadow area and the background area, an iterative algorithm is used to optimize the optimal threshold, and the optimal threshold is used to separate the shadow area and the non-shadow area. When the difference in grayscale values ​​between the shadow and non-shadow areas is small, the false detection rate can be effectively reduced and the accuracy of shadow detection can be improved.

[0181] Step S4: According to a preset coherence coefficient detection method, the one-scene cut image set and the two-scene cut image set are processed to obtain a second shadow area set and a second shadow undetermined area set.

[0182] In some embodiments of the present application, step S4 can be implemented by the following preferred implementations:

[0183] S41: According to the one-scene cut image set and the two-scene cut image set, the coherence value of each two-scene image after the phase correction factor is introduced is calculated, and the calculation method is as follows:

[0184] Where z1, z2 represent the complex values ​​of the two images; * represents the conjugate operator; i represents the pixel number; represents the phase correction factor, Indicates the phase caused by subtracting the terrain factor from the original phase δ flattened represents the coherence value after the phase correction factor is introduced, LL is the number of pixels in the window, and j is the imaginary unit.

[0185] In step S41, if the phase correction factor is not introduced, the conventional coherence value Δ is directly calculated or the maximum likelihood coherence value δ of the pixels in the fixed window is calculated. L , as shown below:

[0186] The calculation of these two coherence values ​​does not fully take into account the image of terrain factors, and the deviation will be larger. Therefore, this application introduces a phase correction factor for calculation.

[0187] S42: processing the coherence value after the phase correction factor is introduced according to a preset spatial averaging algorithm to obtain a final coherence value;

[0188] S43: Processing the first scene cut image set and the second scene cut image set according to the final coherence value to obtain a second shadow area set and a second shadow undetermined area set.

[0189] In step S4 of the present embodiment, the coherence coefficient of the two images is calculated, and the phase correction factor and the spatial neighborhood average are used to improve the unbiasedness of the coherence estimation; the coherence threshold is set and the shadow area is extracted, so that the low-coherence shadow area can be identified more accurately and the deviation caused by terrain or noise can be reduced.

[0190] Step S5: Based on the first shadow area set, the first shadow area set to be determined, the second shadow area set and the second shadow area set to be determined, a set operation is performed to obtain a first partial shadow area set and a third shadow area set to be determined.

[0191] In some embodiments of the present application, step S5 can be implemented by the following preferred implementations:

[0192] S51: according to the first shadow area set and the second shadow area set, obtaining an intersection to obtain a first partial shadow area set;

[0193] S52: after finding the intersection of the first shadow area set to be determined and the second shadow area set to be determined, subtract the first part of the shadow area set from the intersection to obtain a third part of the shadow area set to be determined;

[0194] S53: according to the third partial shadow undetermined region set, the first shadow undetermined region set and the second shadow undetermined region set, obtaining a third shadow undetermined region set by finding the intersection.

[0195] In step S5 of this embodiment, the advantages of different methods are integrated through set operations to improve the accuracy and robustness of shadow detection.

[0196] Step S6: processing the third shadow undetermined area set according to a preset radar line of sight and digital elevation model analysis algorithm to obtain a second partial shadow area set.

[0197] The pixel value in the SAR image reflects the echo energy value received by the receiver after the electromagnetic wave interacts with the target point. The incident angle has a greater impact on the scattering coefficient. According to the radar line of sight and combined with the digital elevation model data, the shadow area concentrated in the third shadow undetermined area can be effectively detected.

[0198] In some embodiments of the present application, Figure 5 A schematic diagram of the geometric relationship between radar sight line and digital elevation model analysis algorithm is shown, and step S6 can be implemented by the following preferred implementation:

[0199] S61: Along the slant range direction, at two points A and B, assuming that the slant range between the two pixel points is M, the elevation difference between the two points is H, and the radar depression angle is set to β1, when these three variables satisfy: H / M=sinβ1, then it can be determined that the point is a shadow area.

[0200] Furthermore, due to the error in the measurement data, the H / M value of the shadow is relatively large, so the threshold condition of H / M should be changed to: H / M ≥ ksinβ1; where 0 <k≤1。

[0201] S62: According to the above threshold condition, the third shadow undetermined area set is processed to obtain a second partial shadow area set.

[0202] In step S6 of the present embodiment, the geometric relationship of the objects in the image is analyzed by combining the radar line of sight and the digital elevation model data, and the shadow area is determined by calculating the relationship between the ground height difference and the slant range. The shadow area caused by the insufficient radar depression angle is further detected, thereby improving the detection accuracy in areas with large terrain undulations.

[0203] Step S7: obtaining a final shadow area detection result by performing a union calculation based on the first partial shadow area set and the second partial shadow area set.

[0204] In step S7 of this embodiment, the SAR image shadow detection result fused by multiple methods is output, which significantly improves the accuracy and comprehensiveness of shadow detection.

[0205] In summary, compared with the prior art, the embodiments of the present application have the following beneficial effects: by processing the original image through the wavelet transform decomposition algorithm, the noise can be effectively removed while retaining important features such as edges, and the bilateral filtering process can smooth the image while maintaining the edges, thereby improving the quality of the image; a Markov random field model is constructed to make full use of the spatial dependence between pixels to improve the accuracy of segmentation, and a region growing algorithm is used to more accurately segment the target area by selecting seed points and gradually expanding the area; an amplitude adaptive threshold detection method is adopted to adapt to the image characteristics in different scenes and improve the robustness of detection; coherence calculation can effectively identify shadow areas caused by physical characteristics, especially those caused by terrain or object occlusion; the results obtained by multiple methods are set and operated, which can combine the advantages of various methods, reduce false detection and missed detection, and improve the overall detection effect; through radar line of sight and digital elevation model analysis, the shadow area caused by insufficient radar depression angle is further detected, which improves the detection accuracy in areas with large terrain undulations; finally, the SAR image shadow detection result fused by multiple methods is output, which significantly improves the accuracy and comprehensiveness of shadow detection.

[0206] Embodiment 2:

[0207] Please refer to Figure 2 Based on the same inventive concept, a multi-method fusion SAR image shadow detection system disclosed in an embodiment of the present invention includes: a preprocessing module M1, a two-step segmentation module M2, an amplitude detection module M3, a coherence detection module M4, a set operation module M5, a radar elevation detection module M6 and a result output module M7;

[0208] The preprocessing module M1 is used to process the original two-view image set according to a preset wavelet multi-scale bilateral filtering algorithm to obtain a two-view image set.

[0209] The preprocessing module M1 includes: a wavelet decomposition unit, a bilateral filtering unit, a smoothing processing unit and an inverse transformation unit; the wavelet decomposition unit processes each original two-scene image in the original two-scene image set according to a preset wavelet transform decomposition algorithm to obtain each first high-frequency sub-band and a first low-frequency sub-band; the bilateral filtering unit is used to process each first high-frequency sub-band according to a preset bilateral filtering algorithm to obtain each second high-frequency sub-band; the smoothing processing unit is used to process each first low-frequency sub-band according to a preset smoothing processing algorithm to obtain each second low-frequency sub-band; the inverse transformation unit is used to process each second high-frequency sub-band and a second low-frequency sub-band according to a preset wavelet inverse transform reconstruction algorithm to obtain each two-scene image, that is, to obtain a two-scene image set.

[0210] Preferably, the processing of the original two-view image in the original two-view image set according to a preset wavelet transform decomposition algorithm to obtain the first high-frequency sub-band and the first low-frequency sub-band can be achieved by the following preferred implementations:

[0211] Wavelet functions of the wavelet transform decomposition algorithm Should meet:

[0212] in, is the basis function or mother wavelet, which can form a cluster of function systems through scaling and translation on the time axis: in, It is called a sub-wavelet, a is the scale factor, reflecting the period length of the wavelet, b is the translation factor, reflecting the time translation, R represents a real number, and L2(R) represents a real two-dimensional space.

[0213] Preferably, the calculation formula of wavelet transform can be as follows:

[0214] Among them, w f (a, b) represents continuous wavelet transform, f(t) represents the function to be transformed by wavelet, a represents the scale factor, b represents the time shift factor, Represents a wavelet transform function that has undergone time translation and scale expansion, also known as a sub-wavelet.

[0215] Preferably, the processing of each of the first high frequency sub-bands according to a preset bilateral filtering algorithm to obtain each of the second high frequency sub-bands can be achieved by the following preferred implementations:

[0216] Among them, the calculation formula of bilateral filtering can be as follows:

[0217] Where p and q represent the coordinates of two pixels, namely p = (x0, y0), q = (x1, y1), BF[I] P represents the bilateral filtering centered at pixel p, I p ,I q Respectively represent the pixel values ​​at positions p and q, that is, the intensity values ​​of the pixels, represents Gaussian weighting in the spatial domain, Represents Gaussian weighting within the pixel value domain, W p is a normalization factor used to ensure that the weighted sum of the pixels is 1, W p The expression is as follows: Q represents the set of all pixels in a local neighborhood centered on pixel p.

[0218] Preferably, the processing of each of the first low-frequency sub-bands according to a preset smoothing algorithm to obtain each of the second low-frequency sub-bands can be achieved by the following preferred implementations:

[0219] The smoothing algorithm may be a mean filter: Among them, M represents a certain pixel, Num represents the total number of pixels, and f(i,j) represents the pixel value.

[0220] Preferably, when designing a linear smoothing filter, it should be ensured that the filter has a single peak, ie, a main lobe, and is symmetrical in both vertical and horizontal directions. For example, a typical smoothing filter has a window of 3×3.

[0221] Preferably, the smoothing algorithm may also use a Gaussian filter, wherein a first-order zero-mean Gaussian function may be expressed as: The two-dimensional discrete zero-mean Gaussian function can be expressed as: The distribution exponent σ represents the width of the Gaussian distribution.

[0222] The preprocessing module M1 of this embodiment processes the original image through the wavelet transform decomposition algorithm, which can effectively remove noise while retaining important features such as edges. The bilateral filtering process can smooth the image while maintaining the edges, thereby improving the image quality.

[0223] The two-step segmentation module M2 is used to process the two-scene image sets according to a preset Markov random field and region growing algorithm to obtain a one-scene cut image set and a two-scene cut image set.

[0224] The two-step segmentation module includes: a normalization unit, a random field construction unit, a first optimization unit, a region growing unit and a morphological processing unit; the normalization unit is used to process the one-scene image set and the two-scene image set respectively according to the normalization algorithm to obtain the one-scene normalized image set and the two-scene normalized image set; the random field construction unit is used to construct Markov random field models according to the one-scene normalized image set and the two-scene normalized image set respectively; wherein the Markov random field model includes a model of a labeled field; the first optimization unit is used to process each of the Markov random field models according to a preset optimization algorithm to obtain a one-scene preliminary segmentation image set and a two-scene preliminary segmentation image set; the region growing unit is used to process the one-scene preliminary segmentation image set and the two-scene preliminary segmentation image set respectively according to a preset automatic region growing algorithm to obtain a one-scene secondary segmentation image set and a two-scene secondary segmentation image set; the morphological processing unit is used to process the one-scene secondary segmentation image set and the two-scene secondary segmentation image set respectively according to a preset morphological algorithm to obtain a one-scene cut image set and a two-scene cut image set.

[0225] Preferably, the normalization algorithm may be as follows:

[0226] Among them, I represents the pixel value in the original image; Min and Max represent the maximum and minimum values ​​of the pixels in the original image, and I N Indicates the normalized value.

[0227] Preferably, the construction of the Markov random field model can be achieved by the following preferred implementations:

[0228] Define a Markov random field, regard each pixel of each image in the one-view normalized image set and the two-view normalized image set as a random variable, model the grayscale distribution in the image according to the finite Gaussian mixture model, and construct an original Markov random field model; wherein the probability distribution of each pixel in the image can be expressed as:

[0229]

[0230] Where S represents the set of random fields, i represents the subset of random fields; P(Y|X=x) represents the probability distribution function or likelihood function of the feature field when the label field is known; y i represents the features in the random field subset i; k = 1, 2, 4, ... L represents the classification label; w represents the dimension of the observed feature y; μ k and∑ k is the expectation and variance of the k-th feature, π k It represents the proportion of pixels with segmentation mark k, and T represents vector transpose.

[0231] Preferably, the expectation maximization algorithm can be used to optimize the parameters of the original Markov random field model, which is achieved by the following preferred implementation:

[0232] The energy function corresponding to the likelihood function of a pixel in the image is expressed as follows:

[0233]

[0234] The parameters θ of the likelihood function are estimated according to the maximum pseudo-likelihood method based on expected maximization, iterating in the following way:

[0235]

[0236] in, Indicates the proportion of pixels with k-type features in this iteration; θ k represents the parameters of the likelihood function, θ k ={μ k ,∑ k}; P(k|y i ,θ k) represents the conditional probability of each pixel; Represents the expectation of k-class features in this iteration; Represents the expectation of k-class features in this iteration; Represents the variance of the feature marked as k in this iteration.

[0237] Preferably, the Gibbs distribution can be used as the prior probability distribution of the labeling field, which is expressed as follows:

[0238] P(x)=Z -1 e -U (x); Z = ∑ x∈Ω e -U (x); U(x) = ∑ c∈C V c (x); where Z is the normalized value of the partitioning function, and the energy function U(x) is a series of potential functions V defined in the group c. c (x); the complement C represents the total set of groups; P(x) is the probability of a certain configuration x appearing in the random field X, and Ω is the set of all possible state vectors.

[0239] Among them, the energy function U(x) can be calculated by the following sum:

[0240]

[0241] Among them, V1, V2, V3 represent potential functions of different orders, and the order represents the complexity of the Gibbs distribution model; C1, C2, C3 represent the group sets corresponding to different orders; x i , x j , x k Represents the random variable in the pixel or group at position i, j, k.

[0242] Preferably, when only the second-order neighborhood system is considered, the energy function U(x) is expressed as:

[0243] Preferably, the Bode model can be used as the model of the marking field, and finally the Markov random field model is obtained, wherein the potential function is defined as follows:

[0244]

[0245] The local probability of the Porter model is:

[0246] Where n i (x i ) is the number of cells in the neighborhood of position i that is equal to x i is the number of pixel neighborhood positions, and β′ is the coupling coefficient.

[0247] Preferably, the preferred algorithm can be implemented by the following preferred implementations:

[0248] Adaptive weights are introduced, where the weight estimation model can be as follows:

[0249]

[0250] Where α(t) is the energy weight of the feature field, C is a constant, t is the number of iterations, L is the number of categories, and x i and i represents the i-th pixel of the first scene image and the i-th pixel of the second scene image respectively, the subscript i represents the pixel position of the image, U(x i ) represents pixel x i The corresponding label field energy function, P represents the probability function, and Ni represents the neighborhood pixel set of pixel i.

[0251] The penalty function is introduced and defined as follows:

[0252]

[0253] The label field energy function U(x i )for:

[0254] The final energy function expression of the posterior probability solved pixel by pixel is:

[0255]

[0256] According to the preset number of classifications and the K-means algorithm, the images in the one-scene normalized image set and the two-scene normalized image set are respectively initially segmented to obtain an initial label for each pixel;

[0257] According to a preset number of iterations, each of the Markov random field models is iteratively optimized to obtain a set of preliminary segmented images of one scene and a set of preliminary segmented images of two scenes, and in each iteration, according to a preset expectation algorithm, each model parameter of the Markov random field model is estimated and calculated, and according to a preset iterative algorithm, the maximum likelihood parameter estimate of each category model is calculated until a preset convergence criterion is met, and the characteristic field energy weight under the current number of iterations and the model parameters of each category calculated by the expectation maximization algorithm are calculated, and substituted into the final energy function expression to update the segmentation result of each pixel; wherein the convergence criterion may be that the absolute value of the parameter difference before and after the update is sufficiently small.

[0258] Preferably, the automatic region growing algorithm can be implemented in the following preferred manner:

[0259] According to a preset rule, a number of seed points are selected from the images in each of the one-scene preliminary segmentation image set and the two-scene preliminary segmentation image set;

[0260] According to the preset initial window size, the window size is gradually increased, and the texture difference of each pixel under each window size is calculated in turn, and the window size with the smallest texture difference is taken as the optimal window size; wherein the texture difference can be obtained by calculating the joint entropy, and the calculation method is as follows:

[0261]

[0262] in, represents the joint co-occurrence matrix of adjacent pixels in the same window, Represents the window texture corresponding to adjacent pixels p and q; Represents the gray level co-occurrence value; N c is the number of texture types, E n represents the joint entropy.

[0263] According to the preset initial window, the probability distribution of each pixel belonging to each texture class is calculated, and the entropy of the initial window is calculated according to the probability distribution, wherein the calculation method is as follows:

[0264] Among them, α represents the set of texture classes, n represents the size of the window group centered on each pixel, and ε represents the entropy value.

[0265] Furthermore, we define the distance between two entropies as

[0266] in, represents the standard deviation of each initial window class; this entropy distance expression measures the average similarity of two texture objects and quantifies the effectiveness of the segmentation. Once the best window is established, the average value of each initial window class is calculated. and standard deviation Therefore, if the entropy value of a pixel (ε p ) is less different from its average value, then the pixel is aggregated into one region.

[0267] According to the optimal threshold and the optimal window size, the texture difference between each pixel point in the optimal window where various sub-points are located and the corresponding seed point is calculated, and the pixel points whose texture difference meets the optimal threshold are screened and merged into the initial area of ​​the corresponding seed point until all pixel points are screened, thereby obtaining a one-scene secondary segmentation image set and a two-scene secondary segmentation image set respectively; wherein the texture difference is obtained by calculating the difference in grayscale co-occurrence values ​​of each pixel point in the window, and the optimal threshold is calculated by the standard deviation of each initial window class.

[0268] Furthermore, the optimal threshold can be calculated in the following way:

[0269] The optimal threshold is calculated as follows: Among them, ∈ α represents the optimal threshold of texture type a, ε′ p represents the joint entropy value of pixel p, μ′ n represents the entropy mean of the initial window of this class, p→Classτ means aggregating pixels p into τ classes, Represents the standard deviation of the texture of type A for each initial window.

[0270] Preferably, the morphological algorithm can be completed by opening and closing operations:

[0271] Among them, the opening operation of the structure element B on the set A (expressed as ) is defined as: It means that B first corrodes A, and then B expands the corrosion result;

[0272] The closing operation of the structure element B on the set A (expressed as A·B) is defined as: It means that B first expands A, and then B corrodes the expansion result. represents the dilation operation, Represents the erosion operation.

[0273] The two-step segmentation module M2 of this embodiment can well capture the spatial correlation between pixels in the image by constructing a Markov random field model, perform global optimization, effectively reduce noise interference, and provide fine preliminary segmentation results; region growing can further refine the segmented area, especially in the edge or transition area, to ensure accurate division of shadow and non-shadow areas; morphological processing can eliminate isolated noise points in the segmentation results, smooth the boundaries, and improve segmentation accuracy and completeness.

[0274] The amplitude detection module M3 is used to process the scene cut image set according to a preset amplitude adaptive threshold detection method to obtain a first shadow area set and a first shadow undetermined area set.

[0275] The amplitude detection module M3 includes: an original acquisition unit, a compression unit, an optimal threshold calculation unit and an amplitude segmentation unit. The original acquisition unit is used to obtain the original amplitude image corresponding to each image in the one-scene cut image set; the compression unit is used to compress the original amplitude image according to the grayscale compression algorithm to obtain the amplitude image; the optimal threshold calculation unit is used to repeatedly divide the shadow area and the non-shadow area of ​​the amplitude image according to the preset initial threshold, and after each division, calculate the mean and standard deviation of the shadow area and the non-shadow area, and update the initial threshold according to the mean and standard deviation until the update change of the initial threshold meets the preset difference range to obtain the optimal threshold; the amplitude segmentation unit is used to process the amplitude image according to the optimal threshold to obtain the first shadow area set and the first shadow undetermined area set.

[0276] Preferably, the optimal threshold calculation unit can be implemented by the following preferred implementations:

[0277] Assuming that the amplitude image contains two grayscale areas, z represents the grayscale value in the image, and p(z) represents the probability density function corresponding to the histogram, the mixed probability density function can be expressed as: p(z) = P1p1(z) + P2p2(z); where P1 represents the probability that the pixel belongs to the shadow area, P2 represents the probability that the pixel belongs to the non-shadow area, and P1+P2=1;

[0278] Furthermore, assuming that there is a threshold Γ, which minimizes the average error in determining pixels as shadow areas or non-shadow areas, then the probability of mistakenly determining a pixel point belonging to a non-shadow area as a pixel point belonging to a shadow area is: The probability of judging a pixel belonging to the shadow area as a pixel in the non-shadow area is: Then the overall error probability is: E(Γ)=P2E1(Γ)+P1E2(Γ).

[0279] Furthermore, in order to find the threshold that minimizes E(Γ), it is necessary to find the first-order derivative of E(Γ) with respect to Γ and make the first-order derivative equal to 0, and we can obtain: P1p1(Γ)=P2p2(Γ).

[0280] Substituting the Gaussian probability density function into the mixed probability density function, we can obtain:

[0281] Among them, μ1 and represents the mean and variance of the Gaussian probability density function corresponding to the pixels in the shadow area; μ2 and Represents the mean and variance of the Gaussian probability density function corresponding to the pixels in the non-shadow area.

[0282] Furthermore, the mean and standard deviation of the shaded and non-shaded areas can be calculated as follows:

[0283]

[0284] Among them, μ S , represents the mean and standard deviation of the shaded area, ∑ (i，j) f S (i, j) represents the sum of the pixel values ​​in the shadow area, Num S Represents the total number of pixels in the shadow area, μ B , represents the mean and standard deviation of the non-shaded area, ∑ (i，j) f B (i, j) represents the sum of the pixel values ​​in the non-shadow area, Num B Represents the total number of pixels in the non-shadow area, ∑ (i，j) (f S (i, j)-μ S ) 2 represents the sum of squares of the difference between the pixels in the shadow area and the mean, ∑ (i，j) (f B (i, j)-μ B ) 2 Represents the sum of squares of the differences between pixels in the non-shadow area and the mean;

[0285] According to P1p1(Γ)=P2p2(Γ), we can get: AΓ 2 +BΓ+C=0; where

[0286] Furthermore, when the variances are the same, We can get:

[0287] When P1=P2, the threshold Γ is the average of the means of the two distributions.

[0288] When the above improved algorithm is applied to a small area in an image, the threshold Γ corresponding to a pixel can be obtained by calculating the grayscale average of the pixels in a window centered on the pixel, as shown below:

[0289] Among them, z(i, j) represents the pixel value at (i, j), and m and n represent the size of the window.

[0290] The amplitude detection module M3 of this embodiment calculates the mean and variance of the shadow area and the background area, uses an iterative algorithm to optimize and obtain the optimal threshold, and uses the optimal threshold to separate the shadow area and the non-shadow area. When the difference in grayscale values ​​between the shadow and non-shadow areas is small, it can effectively reduce the false detection rate and improve the accuracy of shadow detection.

[0291] The coherence detection module M4 is used to process the one-scene cut image set and the two-scene cut image set according to a preset coherence coefficient detection method to obtain a second shadow area set and a second shadow undetermined area set.

[0292] The coherence detection module M4 includes: a coherence value calculation unit, a spatial average algorithm processing unit and a coherence cutting unit; the coherence value calculation unit is used to calculate the coherence value of each two-scene image after introducing the phase correction factor according to the one-scene cutting image set and the two-scene cutting image set; the spatial average algorithm processing unit is used to process the coherence value after introducing the phase correction factor according to a preset spatial average algorithm to obtain a final coherence value; the coherence cutting unit is used to process the one-scene cutting image set and the two-scene cutting image set according to the final coherence value to obtain a second shadow area set and a second shadow undetermined area set.

[0293] Preferably, the coherence value of each of the two images after the phase correction factor is introduced is calculated as follows:

[0294] Where z1, z2 represent the complex values ​​of the two images; * represents the conjugate operator; i represents the pixel number; represents the phase correction factor, Indicates the phase caused by subtracting the terrain factor from the original phase δ flattened represents the coherence value after the phase correction factor is introduced, LL is the number of pixels in the window, and j is the imaginary unit.

[0295] The coherence detection module M4 of this embodiment calculates the coherence coefficient of the two images and uses the phase correction factor and the spatial neighborhood average to improve the unbiasedness of the coherence estimation; sets the coherence threshold and extracts the shadow area, which can more accurately identify the low-coherence shadow area and reduce the deviation caused by terrain or noise.

[0296] The set operation module M5 is used for performing set operation to obtain a first partial shadow area set and a third shadow area set according to the first shadow area set, the first shadow area set to be determined, the second shadow area set and the second shadow area set to be determined.

[0297] The set operation module M5 includes: a first operation unit, a second operation unit and a third operation unit; the first operation unit is used to obtain a first partial shadow area set by finding the intersection of the first shadow area set and the second shadow area set; the second operation unit is used to obtain a third partial shadow area set by subtracting the first partial shadow area set from the first shadow area set to be determined and the second shadow area set to be determined; the third operation unit is used to obtain a third partial shadow area set by finding the intersection of the third partial shadow area set, the first shadow area set to be determined and the second shadow area set to be determined.

[0298] The set operation module M5 of this embodiment integrates the advantages of different methods through set operations to improve the accuracy and robustness of shadow detection.

[0299] The radar elevation detection module M6 is used to process the third shadow undetermined area set according to a preset radar line of sight and digital elevation model analysis algorithm to obtain a second partial shadow area set.

[0300] Preferably, the processing of the third shadow undetermined area set according to the preset radar line of sight and digital elevation model analysis algorithm to obtain the second partial shadow area set can be achieved by the following preferred implementations:

[0301] Along the slant range direction, at points A and B, assuming that the slant range between the two pixel points is M, the elevation difference between the two points is H, and the radar depression angle is set to β1, when these three variables satisfy: H / M = sinβ1, then the point can be determined to be a shadow area;

[0302] Furthermore, due to the error in the measurement data, the H / M value of the shadow is relatively large, so the threshold condition of H / M should be changed to: H / M ≥ ksinβ1; where 0 <k≤1;

[0303] According to the above threshold condition, the third shadow undetermined area set is processed to obtain a second partial shadow area set.

[0304] The radar elevation detection module M6 of this embodiment combines the radar line of sight and the digital elevation model data to analyze the geometric relationship of the objects in the image, determines the shadow area by calculating the ground height difference and the slant range relationship, and further detects the shadow area caused by insufficient radar depression angle, thereby improving the detection accuracy in areas with large terrain undulations.

[0305] The result output module M7 is used for calculating the final shadow area detection result based on the union of the first partial shadow area set and the second partial shadow area set.

[0306] The result output module M7 of this embodiment significantly improves the accuracy and comprehensiveness of shadow detection by outputting the SAR image shadow detection result fused by multiple methods.

[0307] In summary, compared with the prior art, the embodiments of the present application have the following beneficial effects: by processing the original image through the wavelet transform decomposition algorithm, the noise can be effectively removed while retaining important features such as edges, and the bilateral filtering process can smooth the image while maintaining the edges, thereby improving the quality of the image; a Markov random field model is constructed to make full use of the spatial dependence between pixels to improve the accuracy of segmentation, and a region growing algorithm is used to more accurately segment the target area by selecting seed points and gradually expanding the area; an amplitude adaptive threshold detection method is adopted to adapt to the image characteristics in different scenes and improve the robustness of detection; coherence calculation can effectively identify shadow areas caused by physical characteristics, especially those caused by terrain or object occlusion; the results obtained by multiple methods are set and operated, which can combine the advantages of various methods, reduce false detection and missed detection, and improve the overall detection effect; through radar line of sight and digital elevation model analysis, the shadow area caused by insufficient radar depression angle is further detected, which improves the detection accuracy in areas with large terrain undulations; finally, the SAR image shadow detection result fused by multiple methods is output, which significantly improves the accuracy and comprehensiveness of shadow detection.

[0308] The specific working process of each module described above can refer to the corresponding process in the aforementioned method embodiment, which will not be repeated here. The module division is only a logical function division, and there may be other division methods in actual implementation, such as multiple modules can be combined or integrated into another system.

[0309] The specific embodiments described above further illustrate the purpose, technical solutions and beneficial effects of the present invention. It should be understood that the above description is only a specific embodiment of the present invention and is not intended to limit the scope of protection of the present invention. It is particularly pointed out that for those skilled in the art, any modification, equivalent replacement, improvement, etc. made within the spirit and principle of the present invention should be included in the scope of protection of the present invention.

Claims

1. A multi-method fusion SAR image shadow detection method, characterized in that: include: According to the preset wavelet multi-scale bilateral filtering algorithm, the original two-view image set is processed to obtain the two-view image set; According to a preset Markov random field and region growing algorithm, the two-view image sets are processed to obtain a one-view cut image set and a two-view cut image set; According to a preset amplitude adaptive threshold detection method, the one-scene cut image set is processed to obtain a first shadow area set and a first shadow undetermined area set; According to a preset coherence coefficient detection method, the one-scene cut image set and the two-scene cut image set are processed to obtain a second shadow area set and a second shadow undetermined area set; According to the first shadow area set, the first shadow undetermined area set, the second shadow area set and the second shadow undetermined area set, a set operation is performed to obtain a first partial shadow area set and a third shadow undetermined area set; According to a preset radar line of sight and digital elevation model analysis algorithm, the third shadow undetermined area set is processed to obtain a second partial shadow area set; The final shadow area detection result is obtained by calculating the union of the first partial shadow area set and the second partial shadow area set.

2. The multi-method fusion SAR image shadow detection method according to claim 1, characterized in that: The method of processing the original two-view image set according to the preset wavelet multi-scale bilateral filtering algorithm to obtain the two-view image set includes: According to a preset wavelet transform decomposition algorithm, the original two-view image in the original two-view image set is processed respectively to obtain each first high-frequency sub-band and first low-frequency sub-band; According to a preset bilateral filtering algorithm, the first high frequency sub-bands are processed respectively to obtain second high frequency sub-bands; According to a preset smoothing algorithm, the first low-frequency sub-bands are processed respectively to obtain second low-frequency sub-bands; According to a preset wavelet inverse transform reconstruction algorithm, each of the second high-frequency sub-bands and the second low-frequency sub-bands is processed to obtain each of the two-scene images, that is, to obtain a two-scene image set.

3. The multi-method fusion SAR image shadow detection method according to claim 2, characterized in that: The method of processing the two-view image set according to a preset Markov random field and a region growing algorithm to obtain a one-view cut image set and a two-view cut image set includes: According to the normalization algorithm, the one-view image set and the two-view image set are processed respectively to obtain the one-view normalized image set and the two-view normalized image set; According to the one-scene normalized image set and the two-scene normalized image set, a Markov random field model is constructed respectively; wherein the Markov random field model includes a model of a label field; According to a preset optimization algorithm, the Markov random field models are processed respectively to obtain a one-scene preliminary segmentation image set and a two-scene preliminary segmentation image set; According to a preset automatic region growing algorithm, the first scene preliminary segmentation image set and the second scene preliminary segmentation image set are processed respectively to obtain the first scene secondary segmentation image set and the second scene secondary segmentation image set; According to a preset morphological algorithm, the one-scene secondary segmentation image set and the two-scene secondary segmentation image set are processed respectively to obtain a one-scene cutting image set and a two-scene cutting image set.

4. The multi-method fusion SAR image shadow detection method according to claim 3, characterized in that: The Markov random field model is constructed according to the one-scene normalized image set and the two-scene normalized image set respectively; wherein the Markov random field model includes a model of a label field, including: A Markov random field is defined, each pixel of each image in the one-view normalized image set and the two-view normalized image set is regarded as a random variable, the grayscale distribution in the image is modeled according to a finite Gaussian mixture model, and an original Markov random field model is constructed; Optimizing the parameters of the original Markov random field model using an expectation-maximization algorithm; The Markov random field model is obtained by using Gibbs distribution as the prior probability distribution of the marking field and using the Porter model as the model of the marking field.

5. The multi-method fusion SAR image shadow detection method according to claim 4, characterized in that: The Markov random field models are processed according to a preset optimization algorithm to obtain a first scene preliminary segmentation image set and a second scene preliminary segmentation image set, including: Adaptive weights are introduced, where the weight estimation model is as follows: Where α(t) is the energy weight of the feature field, C is a constant, t is the number of iterations, L is the number of categories, and x i and i represents the i-th pixel of the first scene image and the i-th pixel of the second scene image respectively, the subscript i represents the pixel position of the image, U(x i ) represents pixel x i The corresponding label field energy function, P represents the probability function, N i Represents the set of neighboring pixels of pixel i; Introducing penalty function The definition is as follows: The label field energy function U(x i )for: Where β is a custom weighting coefficient ranging from 0 to 1, V2(x i ,x l ) represents the weight of the edge corresponding to the pixels, and exp(*) represents the exponential operation; The final energy function expression of the posterior probability solved pixel by pixel is: According to the preset number of classifications and the K-means algorithm, the images in the one-scene normalized image set and the two-scene normalized image set are respectively initially segmented to obtain an initial label for each pixel; According to the preset number of iterations, each of the Markov random field models is iteratively optimized to obtain a set of preliminary segmented images of one scene and a set of preliminary segmented images of two scenes, and in each iteration, according to the preset expectation algorithm, the model parameters of the Markov random field model are estimated and calculated, and according to the preset iterative algorithm, the maximum likelihood parameter estimate of each category model is calculated until the preset convergence standard is met, and the characteristic field energy weight under the current number of iterations and the model parameters of each category calculated by the expectation maximization algorithm are calculated, and substituted into the final energy function expression to update the segmentation result of each pixel.

6. The multi-method fusion SAR image shadow detection method according to claim 5, characterized in that: The method of processing the first scene preliminary segmentation image set and the second scene preliminary segmentation image set respectively according to a preset automatic region growing algorithm to obtain the first scene secondary segmentation image set and the second scene secondary segmentation image set comprises: According to a preset rule, a number of seed points are selected from the images in each of the one-scene preliminary segmentation image set and the two-scene preliminary segmentation image set; According to the preset initial window size, the window size is gradually increased, and the texture difference of each pixel under each window size is calculated in turn, and the window size with the smallest texture difference is taken as the optimal window size; According to a preset initial window, the probability distribution of each pixel belonging to each texture class is calculated, and according to the probability distribution, the entropy of the initial window is calculated; According to the optimal threshold and the optimal window size, the texture difference between each pixel point in the optimal window where various sub-points are located and the corresponding seed point is calculated, and the pixel points whose texture difference meets the optimal threshold are screened and merged into the initial area of ​​the corresponding seed point until all pixel points are screened, thereby obtaining a one-scene secondary segmentation image set and a two-scene secondary segmentation image set respectively; wherein the texture difference is obtained by calculating the difference in grayscale co-occurrence values ​​of each pixel point in the window, and the optimal threshold is calculated by the standard deviation of each initial window class.

7. The multi-method fusion SAR image shadow detection method according to claim 6, characterized in that: The step of processing the scene cut image set according to a preset amplitude adaptive threshold detection method to obtain a first shadow area set and a first shadow undetermined area set includes: Obtaining the original amplitude image corresponding to each image in the scene cut image set; Compress the original amplitude image according to the grayscale compression algorithm to obtain an amplitude image; According to a preset initial threshold, repeatedly divide the shadow area and non-shadow area of the amplitude image. After each division, calculate the mean and standard deviation of the shadow area and non-shadow area, and update the initial threshold according to the mean and standard deviation until the update change of the initial threshold meets the preset difference range to obtain an optimal threshold; Process the amplitude image according to the optimal threshold to obtain a first set of shadow areas and a first set of shadow-undetermined areas.

8. The multi-method fusion SAR image shadow detection method according to claim 7, characterized in that: The processing of the one-scene cut image set and the two-scene cut image set according to the preset coherence coefficient detection method to obtain a second set of shadow areas and a second set of shadow-undetermined areas includes: According to the one-scene cut image set and the two-scene cut image set, calculate the coherence values of each pair of two-scene images after introducing a phase correction factor, and the calculation method is as follows: Where z1, z2 represent the complex values ​​of the two images; * represents the conjugate operator; i represents the pixel number; represents the phase correction factor, Indicates the phase caused by subtracting the terrain factor from the original phase δ flattened represents the coherence value after the phase correction factor is introduced, LL is the number of pixels in the window, and j is the imaginary unit; Process the coherence value after introducing the phase correction factor according to the preset spatial averaging algorithm to obtain a final coherence value; Process the one-scene cut image set and the two-scene cut image set according to the final coherence value to obtain a second set of shadow areas and a second set of shadow-undetermined areas.

9. The multi-method fusion SAR image shadow detection method according to claim 8, characterized in that: The processing of the third set of shadow-undetermined areas according to the preset radar line-of-sight and digital elevation model analysis algorithm to obtain a second set of partial shadow areas includes: Among them, the preset radar line-of-sight and digital elevation model analysis algorithm includes a geometric relationship formula: H / M≥ksinβ1; 0<k≤1; where M, H, and β1 are the slant range, elevation difference, and radar depression angle corresponding to two pixel points respectively; Process the images in the third set of shadow-undetermined areas according to the geometric relationship, and classify the pixel points that meet the geometric relationship formula as shadow areas to obtain a second set of partial shadow areas.

10. A multi-method fusion SAR image shadow detection system, characterized in that: Includes: A preprocessing module, a two-step segmentation module, an amplitude detection module, a coherence detection module, a set operation module, a radar elevation detection module, and a result output module; Among them, the preprocessing module is used to process the original two-scene image set according to the preset wavelet multi-scale bilateral filtering algorithm to obtain a two-scene image set; The two-step segmentation module is used to process the two-scene image set according to the preset Markov random field and region growing algorithm to obtain a one-scene cut image set and a two-scene cut image set; The amplitude detection module is used to process the one-scene cut image set according to the preset amplitude adaptive threshold detection method to obtain a first set of shadow areas and a first set of shadow-undetermined areas; The optical coherence detection module is used to process the one-scene cut image set and the two-scene cut image set according to the preset coherence coefficient detection method to obtain a second set of shadow areas and a second set of shadow-undetermined areas; The set operation module is used to perform set operations on the first set of shadow areas, the first set of shadow-undetermined areas, the second set of shadow areas, and the second set of shadow-undetermined areas to obtain a first set of partial shadow areas and a third set of shadow-undetermined areas; The radar elevation detection module is used to process the third set of shadow-undetermined areas according to the preset radar line-of-sight and digital elevation model analysis algorithm to obtain a second set of partial shadow areas; The result output module is used to calculate the final shadow area detection result based on the union of the first partial shadow area set and the second partial shadow area set.

Images