Water ripple pre-elimination method based on fusion of ripple and environmental characteristics and iterative restoration

Abstract

The invention discloses a method for pre-eliminating water ripples based on the fusion and iterative restoration of ripples and environmental characteristics. The original image is obtained by underwater shooting, and the image is subjected to secondary wavelet decomposition based on the sym2 wavelet to separate the image frequency band based on the ripple and environmental characteristics; the image frequency band based on the ripple characteristics is detected based on the Markov random field, the water ripple is removed based on the expectation maximization method, and the water ripple area is detected and judged based on the Otsu method. If it does not meet the requirements, it is repeated until the area range is met; the image frequency band based on the environmental characteristics is enhanced based on the grayscale linear transformation, and the color balance is integrated based on the background light to obtain the underwater enhanced frequency band after color correction. Then, the underwater enhanced frequency band after color correction is subjected to bilateral filtering denoising, and the image is subjected to wavelet inverse transformation. The method can more accurately identify the water ripple area, solve the influence of water ripples on underwater images, and ensure the quality of underwater image shooting.

Description

Technical Field

[0001] The invention belongs to the technical field of underwater image processing, and in particular relates to a water ripple pre-elimination method based on fusion and iterative restoration of ripples and environmental characteristics. Background Art

[0002] Water ripples are often caused by the movement of underwater objects. Water ripples cause obvious ghost areas to appear in the captured calibration plate image. The existing method requires waiting for the water surface to calm down before shooting. This method is time-consuming and has very demanding requirements on the shooting environment.

[0003] Existing watermark image removal methods are mainly divided into two types: hardware removal and software removal.

[0004] The first hardware method is to remove it. Patent CN 104200447 B discloses a real-time low-light color image enhancement device. This method can effectively improve the visual effect of video images, but it does not take the underwater environment into consideration when performing nonlinear transformation of color components. Patent CN 109253863 B discloses a water ripple test system and method. This method can automatically identify water ripples on the device screen, but it does not solve the visual quality impact of water ripples on the image. Patent CN 105430290 A discloses a water ripple pre-elimination method and camera device. This method eliminates water ripples by adjusting the exposure time of the camera device to a positive integer multiple of the photosensitive period, but the camera device of this method is greatly disturbed by the indoor environment and it is difficult to stably eliminate water ripples.

[0005] The second software method for removing water ripples can be divided into two categories. The first category is a method designed based on the characteristic that water ripples will produce shadows: Patent CN107742171A designs a photovoltaic power station power generation prediction method based on mobile shadow image recognition, using irradiance sensors and shadow recognition methods to obtain irradiance data and shadow images on the photovoltaic array, but this method is only applicable to the field of photovoltaic power generation technology; Gong et al. proposed a user-assisted shadow removal method (Gong H, Darren A, et al. Interactive removal and ground truth for difficult shadow scenes [J]. Journal of the Optical Society of America A Optics Image Science & Vision, 2016), using sampling-based shadow scale to estimate the shadow area, but this method requires a lot of manual assistance and cannot meet the requirements of automatic removal of shadow areas; Gryka et al. proposed a data-driven method to model shadow areas (Parisotto S et al. Anisotropic osmosis filtering for shadow removal in images [J]. Inverse Problems, 2019: 36-39), but this method requires shadow annotations as input, and the conditional constraints are relatively harsh. The second type of software method is designed based on the characteristic that water ripples will reduce the clarity of the image.Traditional image enhancement methods include: histogram enhancement method (KANSAL S et al, Image contrast enhancement using unsharp masking and histogram equalization [J], 2018, 77 (20): 26919-26938) which will increase image detail information and water ripple effect at the same time, resulting in image quality degradation; Retinex-based image contrast enhancement method (LIU C et al, Enhancement of low visibility aerial images using histogram truncation and an explict Retinex representation for balancing contrast and color consistency [J], 2017, 128: 16-26) which can eliminate image distortion under non-uniform lighting conditions, but this method will produce image over-enhancement or halo artifacts during the brightness change process; for this reason, Ancuti et al. proposed a fusion technology for underwater images and videos (Ancuti C, et al. Fusion-based restoration of the underwater images [C]. IEEE Intermational Conference on Image Processing.IEEE, 2011:1557-1560), which has a good enhancement effect on underwater images; He Huan et al. proposed an image defogging method based on the MALLAT method (He Huan. An image defogging method based on the MALLAT method [J]. Journal of Xinjiang Normal University (Natural Science Edition), 2020, 39, 110(01):28-32), which achieves the defogging effect by decomposing and reconstructing the image using two-dimensional discrete wavelets, but the above two methods cannot completely eliminate the influence of water ripples on the image.

[0006] To sum up, the existing hardware methods are greatly constrained by the surrounding environment; the existing software methods have the problem of limited applicable scenarios; the methods designed based on the characteristic that water ripples will produce shadows cannot adaptively and completely eliminate the impact of water ripples on images; the methods designed based on the characteristic that water ripples will reduce the clarity of the picture do not solve the impact of water ripples on underwater images.

[0007] Therefore, the above methods cannot completely remove the influence of water ripples on the image. Summary of the invention

[0008] The purpose of the present invention is to provide a water ripple pre-elimination method based on the fusion and iterative restoration of ripples and environmental characteristics, extract the high and low frequency components of the image based on the sym2 wavelet, eliminate the water ripples using the expectation maximization method, and consider the energy influence of the water ripples to perform image fusion enhancement. The method can well eliminate the influence of water ripples and effectively improve the shooting quality of underwater calibration plates.

[0009] The technical solution to achieve the purpose of the present invention is: a water ripple pre-elimination method based on ripple and environmental characteristics fusion and iterative restoration, which is used to process low-light images taken underwater with water ripples, including the following steps:

[0010] Step (1): Underwater shooting to obtain a water ripple image I0;

[0011] Step (2): Data processing of the water ripple pre-elimination method based on the fusion of ripples and environmental characteristics and iterative restoration, specifically including:

[0012] Step (21): performing secondary wavelet decomposition on the image I0 obtained in step (1) based on the sym2 wavelet to separate an image frequency band based on ripple characteristics and an image frequency band based on environmental characteristics;

[0013] Step (22): using a Markov random field-based method to detect water ripples on the image frequency band based on ripple characteristics separated in step (21);

[0014] Step (23): removing water ripples from the image frequency band obtained in step (22) based on the expectation maximization method;

[0015] Step (24): Detect and judge the water ripple area based on the Otsu method for the image frequency band obtained in step (23). If the water ripple area range is not met, return to step (22) and re-detect the water ripple based on the Markov random field, and iterate until the area range is met;

[0016] Step (25): performing frequency band enhancement on the image frequency band separated based on the environmental characteristics in step (21) based on grayscale linear transformation;

[0017] Step (26): performing color balancing based on the background light fusion on the image frequency band obtained in step (25) to obtain a color-corrected underwater enhanced frequency band;

[0018] Step (27): performing bilateral filtering based denoising on the image frequency band obtained in step (26);

[0019] Step (3): Perform inverse wavelet transform on the image.

[0020] Furthermore, the image I0 in step (1) is expressed as:

[0021] I0(x,y)=D(x,y)+J(x,y)+O(x,y) (1)

[0022] Among them, x and y are the independent variable and dependent variable of the two-dimensional image signal respectively, D(x,y) is the image frequency band based on ripple characteristics, J(x,y) is the image frequency band based on environmental characteristics, and O(x,y) is the noise influence, which includes the scattering and absorption of light by suspended matter in water, interference from artificial ambient light source, and interference from aquatic organisms.

[0023] Furthermore, step (21) "performing secondary wavelet decomposition on the image I0 acquired in step (1) based on the sym2 wavelet to separate the image frequency band based on the ripple characteristics and the image frequency band based on the environmental characteristics" is specifically as follows:

[0024] Assume a separable two-dimensional scaling function And wavelet function ψ(x), construct three two-dimensional basic wavelet functions: in is the level ψ of the wavelet function 1 , vertical ψ 2 , the diagonal ψ 3 The corresponding scaling function satisfies

[0025] The two-dimensional wavelet is scaled and translated by the following formula:

[0026]

[0027] c, d are the stretching and translation of the wavelet;

[0028] The two-dimensional image signal f(x,y) is at scale 2 i The first-level decomposition of the image band D(x,y) based on the ripple characteristics is the low-frequency component L i and high frequency components They are:

[0029]

[0030] Where (a, b) is the pixel point in the image, the two-dimensional sequence is the low-frequency component of the i-th order image of the first-level decomposition, Respectively represent the high-frequency components in the horizontal, vertical, and diagonal directions of the first-order decomposition of the i-th order image;

[0031] Used in scale 2 i The second-order decomposition of the image is obtained by the second-order decomposition low-frequency component and the high-frequency components in the horizontal direction of the second-level decomposition High frequency components in the vertical direction High frequency components in the diagonal direction Low frequency components after secondary decomposition Corresponding to the image frequency band based on environmental characteristics, the high-frequency component in the horizontal direction High frequency components in the vertical direction High frequency components in the diagonal direction Corresponding to the image frequency bands based on the ripple characteristics.

[0032] Furthermore, the step (22) of "using the Markov random field to detect water ripples on the image frequency band based on ripple characteristics separated in step (21)" is specifically:

[0033] Assume that the high-frequency component in the horizontal direction of the image frequency band based on the ripple characteristic The gradient field is Then the gradient field distribution image I0 is: ζ0 is the Markov random variable at the pixel (a, b), and the posterior energy of ζ0 is It is expressed as:

[0034]

[0035] Among them, F0 represents the probability of the gradient value at the pixel point (a, b) in the image, represents the 8-neighborhood pixel value of pixel point (a, b), (1-ζ0) is the prior energy of the possibility that pixel point (a, b) is an edge pixel, ζ1 is the Markov random variable of the 8-neighborhood pixel value around pixel point (a, b), C0 is the posterior energy compensation, κ(ζ1,ζ0) is the likelihood energy of pixel point (a, b), and the specific expression is:

[0036]

[0037] I0 is the gradient field distribution image of the pixel point (a, b), I1 is the 8-neighborhood gradient field distribution image of the pixel point (a, b), C1 and C2 are the local energy thresholds of the pixel point (a, b) marking the boundary pixels, and C1<C2;

[0038] Based on the minimum a posteriori energy Obtain ζ0 and obtain the high-frequency component in the horizontal direction of the image frequency band based on the ripple characteristics Water ripple area

[0039] Formula (4) is used to obtain the high-frequency component in the vertical direction of the secondary decomposition: High frequency components in the diagonal direction The minimum value of the posterior energy is obtained to obtain the high-frequency component in the vertical direction of the image frequency band based on the ripple characteristics after the second-level decomposition. Water ripple area High frequency components in the diagonal direction Water ripple area

[0040] Furthermore, the step (23) of "removing water ripples from the image frequency band obtained in step (22) based on the expectation maximization method" is specifically as follows:

[0041] Step (231): setting the initial value θ0 of the probability distribution parameter θ;

[0042] Step (232): Desired step: Assume that the high-frequency component in the horizontal direction of the image frequency band based on the ripple characteristic Water ripple area There are N pixels in the n}, F represents the enhancement type corresponding to pixel N, recorded as {F1, F2, ..., F n}; Substitute the initial value θ0 in step (231) into the likelihood formula R(θ0) to express it as:

[0043]

[0044] Among them, w(N,F|θ) is the joint distribution probability model of variables N and F, which is specifically the high-frequency component in the horizontal direction Water ripple area The probability of a pixel belonging to class enhancement type F is

[0045] The above formula is simplified to Calculate the current estimate of the enhancement type F;

[0046] Step (233): Maximization step: Substitute the current estimate of F calculated in step (232) into the maximum likelihood estimation formula:

[0047]

[0048] Repeat the iteration steps (232) and (233) until the two steps converge at the same time, and maximize the likelihood function of formula (6) through iteration to obtain the maximum value of parameter θ1; then calculate the accurate value of the class enhancement type F for N pixels based on the probability model method of w(N|θ1);

[0049] Step (234): The vertical high-frequency component water ripple region of the secondary decomposition and diagonal direction The enhancement type F is obtained and substituted into formula (6) to calculate the vertical high-frequency component of the secondary decomposition in the image frequency band based on the ripple characteristics: Diagonal high frequency components The accurate value of the regional enhancement type; get the estimated horizontal high-frequency component of the secondary decomposition without the influence of water ripples High frequency components in the vertical direction Diagonal high frequency components

[0050] Furthermore, step (24) specifically includes the following steps:

[0051] Step (241): Calculate the global threshold T0 based on the equation maximization

[0052] Based on the first-order Otsu method, the global function V(x) is maximized. The global function V(T0) is expressed as:

[0053]

[0054] in u ε It is the high-frequency component in the horizontal direction of the image frequency band based on the ripple characteristics Water ripple area The probability of a pixel with gray level ε in ε P ε is a probability density function based on grayscale, M(T0) is the probability that the pixel is within the threshold T0, λ(T0) is the probability that the pixel value is higher than the threshold T0, Represents the average value of pixels in the image; this is used to separate the high-frequency components in the horizontal direction of the image frequency band based on the ripple characteristics The shadow and non-shadow areas to be determined

[0055] The vertical high-frequency components of the image frequency band that meet the ripple characteristics are decomposed by the second level High frequency components in the diagonal direction Water ripple area Substitute the probability of pixels with gray level ε in the image into formula (7) to obtain the high-frequency component in the vertical direction of the image frequency band based on the ripple characteristic after the second-level decomposition: High frequency components in the diagonal direction The shadow and non-shadow areas to be determined

[0056] Step (242): determining the optimal threshold T1 based on minimization of regional variance;

[0057] The high-frequency component in the horizontal direction of the image frequency band obtained in step (241) based on the ripple characteristic The shadow and non-shadow areas to be determined The optimal threshold is determined by maximizing the variance; the high-frequency component in the horizontal direction of the image frequency band based on the ripple characteristics The shadow and non-shadow areas to be determined The average intensity u w and regional variance It is expressed as:

[0058]

[0059] L represents the high-frequency component in the horizontal direction of the image frequency band based on the ripple characteristics The shadow and non-shadow areas to be determined The gray level, the global threshold T0 divides the gray level into pixels {1,2,...,l} and pixels {l+1,...,L}, γ k ,u k Represents the high-frequency component in the horizontal direction of the image frequency band based on the ripple characteristics The shadow and non-shadow areas to be determined The cumulative probability and average grayscale of Minimize the high-frequency components in the horizontal direction in this smaller image band based on the ripple characteristics The shadow and non-shadow areas to be determined

[0060] The variance maximization method is used to obtain the high-frequency component in the vertical direction of the secondary decomposition in the image frequency band based on the ripple characteristics after the secondary decomposition. Diagonal high frequency components The shadow and non-shadow areas to be determined

[0061] Step (243): iteratively eliminating the water ripple area;

[0062] According to step (242), the high-frequency components in the horizontal direction of the image frequency band based on the ripple characteristic are respectively Vertical Diagonal direction The shadow and non-shadow areas to be determined The size of steps (22) and (23) is continuously iterated. When the shadow and non-shadow areas are less than the preset threshold condition ν, it indicates that the water ripple area has been eliminated; the absolute value of the preset threshold condition ν satisfies that the grayscale pixel value inside the shadow and non-shadow areas to be determined is not greater than 20% of the external pixel value; it indicates that the water ripple area has been eliminated, and the high-frequency component of the horizontal direction of the secondary decomposition without the influence of water ripples is obtained. High frequency components in the vertical direction High frequency components in the diagonal direction

[0063] Furthermore, step (3) performs inverse wavelet transform on the image as follows:

[0064] The Mallat method using two-dimensional wavelet decomposition is:

[0065]

[0066] Among them, m0, n0 represent the coefficient matrix The rows and columns of

[0067] For low frequency components and the high-frequency components in the horizontal direction High frequency components in the vertical direction High frequency components in the diagonal direction Wavelet decomposition to the highest decomposition order gives Multi-level wavelet transform; when it is known After these four decomposition coefficient matrices, they are reconstructed by the following formula

[0068]

[0069] Fusion and Get the image without water ripple enhancement

[0070]

[0071] Among them, T i c is the fusion coefficient, T i c ∈(0,1), and

[0072] Assume that there is no water ripple enhanced image The fused background light is G1, and the image after water ripple enhancement is processed using a guided filter and substituted into the restoration formula:

[0073]

[0074] Among them, t c (a, b) are RGB channel transmission maps, is the background light of the RGB channel, This is an enhanced image after color balancing combined with background light.

[0075] Compared with the prior art, the present invention has the following significant advantages:

[0076] (1) The water ripple pre-elimination method based on the fusion and iterative restoration of ripple and environmental characteristics processes the ripple characteristics and environmental characteristics separately, which better eliminates the influence of water ripples on the image. At the same time, it can ensure the restoration of image colors underwater, which has the advantage of improving the quality of underwater images.

[0077] (2) The water ripple pre-elimination method based on the fusion and iterative restoration of ripples and environmental characteristics can adaptively identify the water ripple area. It does not need to rely on too much human-computer interaction, nor does it need to rely too much on the shooting operation of multi-angle images. It has the advantages of low complexity and strong applicability. BRIEF DESCRIPTION OF THE DRAWINGS

[0078] Figure 1 This is a flow chart of the water ripple pre-elimination method based on the fusion and iterative restoration of ripples and environmental characteristics of the present invention.

[0079] Figure 2 This is a diagram of the image frequency band processing method based on ripple characteristics of the present invention.

[0080] Figure 3 This is a diagram of the image frequency band processing method based on environmental characteristics of the present invention.

[0081] Figure 4 Schematic diagram of the two-layer wavelet decomposition of an image based on the sym2 wavelet.

[0082] Figure 5 Schematic diagram of wavelet reconstruction based on Mallat method. DETAILED DESCRIPTION

[0083] The following is combined with Figure 1-5 The present invention is described in further detail below.

[0084] like Figure 1 As shown, the present invention provides a water ripple pre-elimination method based on ripple and environmental characteristics fusion and iterative restoration, comprising the following steps:

[0085] Step 1: Get the water ripple image;

[0086] Let an image be I0, which can be expressed as:

[0087] I0(x,y)=D(x,y)+J(x,y)+O(x,y) (1)

[0088] Among them, x and y are the independent variable and dependent variable of the two-dimensional image signal respectively. D(x,y) is the image frequency band based on the ripple characteristics, J(x,y) is the image frequency band based on the environmental characteristics, and O(x,y) is the noise influence. O(x,y) is generally a low-frequency component in the environment.

[0089] Step 2: Data processing of the water ripple pre-elimination method based on the fusion of ripples and environmental characteristics and iterative restoration;

[0090] Step 2.1: Secondary wavelet decomposition of the image based on sym2 wavelet

[0091] Assume a separable two-dimensional scaling function And wavelet function ψ(x), the three two-dimensional basic wavelet functions constructed are: in is the level ψ of the wavelet function 1 , vertical ψ 2 , the diagonal ψ 3 The corresponding scaling function satisfies

[0092] The two-dimensional wavelet is scaled and translated by the following formula:

[0093]

[0094] c, d are the stretching and translation of the wavelet.

[0095] like Figure 4 As shown, the two-dimensional image signal f(x,y) is at scale 2 i The first-level decomposition of the image band D(x,y) based on the ripple characteristics is the low-frequency component L i and high frequency components They are:

[0096]

[0097] Where (a, b) is the pixel point in the image, the two-dimensional sequence is the low-frequency approximation of the i-th order image, Represent the high-frequency details of the image in the horizontal, vertical, and diagonal directions of the i-th order. i The second-order decomposition of the image can be obtained by the second-order decomposition low-frequency component and horizontal direction Vertical Diagonal direction high frequency components.

[0098] Low frequency components after secondary decomposition Corresponding to the image frequency band based on environmental characteristics, the high-frequency component in the horizontal direction Vertical Diagonal direction Corresponding to the image frequency band based on the ripple characteristics. After the secondary decomposition, the following is obtained: the high-frequency component without noise influence of the secondary decomposition and the low-frequency component with noise influence of the secondary decomposition. Finally, in step 2.7, bilateral filtering is used to eliminate the noise influence O(x, y) in the low-frequency component with noise influence of the secondary decomposition. In the following content, it is referred to as the high-frequency component of the secondary decomposition and the low-frequency component of the secondary decomposition.

[0099] Step 2.2: Detect water ripples based on the image frequency band based on ripple characteristics using the Markov random field

[0100] like Figure 2 As shown, let the high-frequency component in the horizontal direction of the image frequency band based on the ripple characteristics in step 2.1 be The gradient field Then the gradient field distribution image I0 is: ζ0 is the Markov random variable at the pixel (a, b), and the posterior energy of ζ0 is It is expressed as:

[0101]

[0102] Among them, F0 represents the probability of the gradient value at the pixel point (a, b) in the image. Represents the 8-neighborhood pixel value of pixel point (a, b). (1-ζ0) is the prior energy of the possibility that pixel point (a, b) is an edge pixel, ζ1 is the Markov random variable of the 8-neighborhood pixel value around pixel point (a, b), C0 is the posterior energy compensation, κ(ζ1,ζ0) is the likelihood energy of pixel point (a, b), and the specific expression is:

[0103]

[0104] I0 is the gradient field distribution image of the pixel point (a, b), I1 is the 8-neighborhood gradient field distribution image of the pixel point (a, b), C1 and C2 are the local energy thresholds of the pixel point (a, b) marking the boundary pixels, and C1<C2.

[0105] Based on the minimum a posteriori energy Obtain ζ0 and obtain the high-frequency component in the horizontal direction of the image frequency band based on the ripple characteristics Water ripple area Use formula (4) to obtain the vertical direction Diagonal direction The minimum value of the posterior energy can be obtained. The vertical direction of the high-frequency component in the image frequency band based on the ripple characteristics after the secondary decomposition in step 2.1 Water ripple area Diagonal direction Water ripple area

[0106] Step 2.3: Remove water ripples based on expectation maximization

[0107] 2.3.1 Setting the initial value θ0 of the probability distribution parameter θ

[0108] 2.3.2 Expected Steps

[0109] Assume that the high-frequency component in the horizontal direction of the image frequency band based on the ripple characteristics obtained in step 2.1 is Water ripple area There are N pixels in the n}, F represents the enhancement type corresponding to pixel N, recorded as {F1, F2, ..., F n}. Substituting the initial value θ0 in step 2.3.1 into the likelihood formula R(θ0) can be expressed as:

[0110]

[0111] Among them, w(N,F|θ) is the probability model of variables N and F.

[0112] The above formula can be simplified to Calculate the current estimate of the enhancement type F.

[0113] 2.3.3 Maximization Step

[0114] Substitute the current estimate of F calculated in step 2.3.2 into the maximum likelihood estimation formula:

[0115]

[0116] Repeat steps 2.3.2 and 2.3.3 until both steps converge, maximizing the likelihood function of formula (6) and obtaining the maximum value of parameter θ1. Then, calculate the exact value of each enhancement type F for N pixels based on the probability model of w(N|θ1). and diagonal direction Obtain the enhancement type F and substitute it into formula (6) to calculate the high-frequency component in the vertical direction of the image frequency band based on the ripple characteristics after the secondary decomposition in step 2.1 High frequency components in the diagonal direction The accurate value of the regional enhancement type. Get the estimated horizontal high-frequency component without water ripples Vertical Diagonal direction

[0117] Step 2.4: Detect water ripples based on Otsu's method

[0118] 2.4.1 Based on the equation maximization, the global threshold T0 can be calculated

[0119] Based on the first-order Otsu method, the global function V(x) is maximized. The global function V(T0) can be expressed as:

[0120]

[0121] in u ε It is the high-frequency component in the horizontal direction of the image frequency band based on the ripple characteristics Water ripple area The probability of a pixel with gray level ε in . ε P ε is a probability density function based on grayscale, M(T0) is the probability that the pixel is within the threshold T0, λ(T0) is the probability that the pixel value is higher than the threshold T0, Represents the average value of pixels in the image. This is used to separate the high-frequency components in the horizontal direction of the image frequency band based on the ripple characteristics. The shadow and non-shadow areas to be determined The high-frequency components in the image frequency band that meet the ripple characteristics are vertically Diagonal direction Water ripple area Substituting the probability of pixels with gray level ε into formula (7) can obtain the vertical direction of the high-frequency component in the image frequency band based on the ripple characteristics after the secondary decomposition Diagonal direction The shadow and non-shadow areas to be determined

[0122] 2.4.2 Determining the optimal threshold T1 based on regional variance minimization

[0123] The high-frequency component in the horizontal direction of the image frequency band based on the ripple characteristics obtained in step 2.4.1 The shadow and non-shadow areas to be determined The optimal threshold is determined by maximizing the variance. The high-frequency component in the horizontal direction of the image frequency band based on the ripple characteristics The shadow and non-shadow areas to be determined The average intensity u w and regional variance It is expressed as:

[0124]

[0125] L represents the high-frequency component in the horizontal direction of the image frequency band based on the ripple characteristics The shadow and non-shadow areas to be determined The gray level, the global threshold T0 divides the gray level into pixels {1,2,...,l} and pixels {l+1,...,L}, γ k ,u k Represents the high-frequency component in the horizontal direction of the image frequency band based on the ripple characteristics The shadow and non-shadow areas to be determined The cumulative probability and average grayscale of Minimize the high-frequency components in the horizontal direction in this smaller image band based on the ripple characteristics The shadow and non-shadow areas to be determined The variance maximization method can be used to obtain the high-frequency components in the vertical direction of the image frequency band based on the ripple characteristics after the second-level decomposition. High frequency components in the diagonal direction The shadow and non-shadow areas to be determined

[0126] 2.4.3 Iterative elimination of water ripple area

[0127] According to step 2.4.2, the high-frequency components in the horizontal direction of the image frequency band based on the ripple characteristics are respectively High frequency components in the vertical direction High frequency components in the diagonal direction The shadow and non-shadow areas to be determined The size of steps 2.2 and 2.3 is continuously iterated. When the shadow and non-shadow areas are less than a preset threshold condition ν, it means that the water ripple area has been eliminated. The absolute value of the threshold condition ν satisfies that the grayscale pixel value inside the shadow and non-shadow areas to be determined is not greater than 20% of the external pixel value. It means that the water ripple area has been eliminated. The high-frequency component of the high-frequency component in the horizontal direction without the influence of water ripples is obtained. High frequency components in the vertical direction High frequency components in the diagonal direction

[0128] Step 2.5: Enhance the image band J(x,y) based on environmental characteristics based on grayscale linear transformation

[0129] like Figure 3 As shown, the low-frequency component of the image based on environmental characteristics is obtained in step 2.1 The grayscale linear transformation function is:

[0130]

[0131] Where [g,h] is the low-frequency component of the image based on environmental characteristics Gray value range, [s, t] is the low-frequency component of the image based on environmental characteristics after stretching The grayscale value range is obtained by using the image frequency band based on the environmental characteristics and the low-frequency component enhanced based on the grayscale linear transformation.

[0132] Step 2.6: Blending color balance based on background light

[0133] Assume that the low-frequency component of the image frequency band based on environmental characteristics in step 2.5 is enhanced based on grayscale linear transformation The global background color is G0, and the low-frequency components of the image band based on the environmental characteristics are enhanced based on the grayscale linear transformation The average intensity value of the background color after background color extraction is recorded as g0. Then the water ripple ratio value q is: The hue tolerance threshold can be expressed as:

[0134]

[0135] Where (a, b) is the pixel value in the water ripple area S, grand(·) represents the gradient operation, which can be expressed as:

[0136]

[0137] Adjust the background color average intensity value of the low-frequency component image. Get the low-frequency component based on the image band of the environmental characteristics and the background light fusion color balance.

[0138] Step 2.7: Denoise the image O(x,y) based on bilateral filtering

[0139] For image bands based on environmental characteristics, low-frequency components based on background light fusion color balance Bilateral filtering is used to remove O(x,y) in the image. Bilateral filtering is based on Gaussian filtering, and adds the influence of image pixel values ​​on the central pixel to jointly affect the value of the central pixel. The two-dimensional Gaussian kernel function G d It is expressed as:

[0140]

[0141] Among them, the pixel (a0, b0) is the low-frequency component of the image frequency band based on the environmental characteristics and the background light fusion color balance The position of the center element. α1, α2 are the parameter factors in the x and y directions.

[0142] By adjusting the parameters α1, α2, the grayscale value of the central pixel and the edge area generate a larger weight. The low-frequency component of the image frequency band based on environmental characteristics is obtained based on bilateral filtering denoising

[0143] Step 3: Perform inverse wavelet transform on the image

[0144] like Figure 5 As shown, the Mallat method using two-dimensional wavelet decomposition is:

[0145]

[0146] Among them, m0, n0 represent the coefficient matrix rows and columns.

[0147] For low frequency components and the high-frequency components in the horizontal direction High frequency components in the vertical direction High frequency components in the diagonal direction When wavelet decomposition is performed to the highest order, we can obtain Multi-level wavelet transform of. When it is known After these four decomposition coefficient matrices, we can reconstruct

[0148]

[0149] Fusion and You can get an image without water ripple enhancement

[0150]

[0151] Among them, T i c is the fusion coefficient, T i c ∈(0,1), and

[0152] Assume that the image without water ripple enhancement obtained in step 2.3 The fused background light is G1, and the guided filter is used to process the image I2 after water ripple enhancement. Substitute into the restoration formula:

[0153]

[0154] Among them, t c (a, b) are RGB channel transmission maps, is the background light of the RGB channel, This is an enhanced image after color balancing combined with background light.

[0155] The above shows and describes the basic principles, main features and advantages of the present invention. Those skilled in the art should understand that the present invention is not limited to the above embodiments. The above embodiments and descriptions are only for explaining the principles of the present invention. Without departing from the spirit and scope of the present invention, the present invention may have various changes and improvements, which fall within the scope of the present invention. The scope of protection of the present invention is defined by the attached claims and their equivalents.

Claims

1. A water ripple pre-elimination method based on ripple and environmental characteristics fusion and iterative restoration, used to process low-light images with water ripples taken underwater, characterized in that: The steps include: Step (1): Underwater shooting to obtain a water ripple image I0; Step (2): Data processing of the water ripple pre-elimination method based on the fusion of ripples and environmental characteristics and iterative restoration, specifically including: Step (21): performing secondary wavelet decomposition on the image I0 obtained in step (1) based on the sym2 wavelet to separate an image frequency band based on ripple characteristics and an image frequency band based on environmental characteristics; Step (22): using a Markov random field-based method to detect water ripples on the image frequency band based on ripple characteristics separated in step (21); Step (23): removing water ripples from the image frequency band obtained in step (22) based on the expectation maximization method; Step (24): Detect and judge the water ripple area based on the Otsu method for the image frequency band obtained in step (23). If the water ripple area range is not met, return to step (22) and re-detect the water ripple based on the Markov random field, and iterate until the area range is met; Step (25): performing frequency band enhancement on the image frequency band separated based on the environmental characteristics in step (21) based on grayscale linear transformation; Step (26): performing color balancing based on the background light fusion on the image frequency band obtained in step (25) to obtain a color-corrected underwater enhanced frequency band; Step (27): performing bilateral filtering based denoising on the image frequency band obtained in step (26); Step (3): Perform inverse wavelet transform on the image.

2. The method according to claim 1, characterized in that The image I0 in step (1) is represented as: I0(x,y)=D(x,y)+J(x,y)+O(x,y) (1) Among them, x and y are the independent variable and dependent variable of the two-dimensional image signal respectively, D(x,y) is the image frequency band based on ripple characteristics, J(x,y) is the image frequency band based on environmental characteristics, and O(x,y) is the noise influence, which includes the scattering and absorption of light by suspended matter in water, interference from artificial environmental light sources, and interference from aquatic organisms.

3. The method according to claim 2, characterized in that Step (21) "performing secondary wavelet decomposition on the image I0 obtained in step (1) based on the sym2 wavelet to separate the image frequency band based on the ripple characteristics and the image frequency band based on the environmental characteristics" is specifically: Assume a separable two-dimensional scaling function And wavelet function ψ(x), construct three two-dimensional basic wavelet functions: ψ 1 (x,y)=ψ(x)ψ(y), in is the level ψ of the wavelet function 1 , vertical ψ 2 , the diagonal ψ 3 The corresponding scaling function satisfies The two-dimensional wavelet is scaled and translated by the following formula: c, d are the stretching and translation of the wavelet; The two-dimensional image signal f(x,y) is at scale 2 i The first-level decomposition of the image band D(x,y) based on the ripple characteristics is the low-frequency component L i and high frequency components They are: Where (a, b) is the pixel point in the image, and the two-dimensional sequence L 1 i is the low-frequency component of the i-th order image of the first-level decomposition, Respectively represent the high-frequency components of the horizontal, vertical, and diagonal directions of the first-order decomposition of the i-th order image; Used in scale 2 i The second-order decomposition of the image is obtained by the second-order decomposition low-frequency component and the high-frequency components in the horizontal direction of the second-level decomposition High frequency components in the vertical direction High frequency components in the diagonal direction Low frequency components after secondary decomposition Corresponding to the image frequency band based on environmental characteristics, the high-frequency component in the horizontal direction High frequency components in the vertical direction High frequency components in the diagonal direction Corresponding to the image frequency bands based on the ripple characteristics.

4. The method according to claim 3, characterized in that The specific step (22) of "using the Markov random field to detect water ripples on the image frequency band based on ripple characteristics separated in step (21)" is: Assume that the high-frequency component in the horizontal direction of the image frequency band based on the ripple characteristic The gradient field is Then the gradient field distribution image I0 is: ζ0 is a Markov random variable at the pixel (a, b), and the posterior energy P about ζ0 is ζ0 It is expressed as: Among them, F0 represents the probability of the gradient value at the pixel point (a, b) in the image, represents the 8-neighborhood pixel value of pixel point (a, b), (1-ζ0) is the prior energy of the possibility that pixel point (a, b) is an edge pixel, ζ1 is the Markov random variable of the 8-neighborhood pixel value around pixel point (a, b), C0 is the posterior energy compensation, κ(ζ1,ζ0) is the likelihood energy of pixel point (a, b), and the specific expression is: I0 is the gradient field distribution image of the pixel point (a, b), I1 is the 8-neighborhood gradient field distribution image of the pixel point (a, b), C1 and C2 are the local energy thresholds of the pixel point (a, b) marking the boundary pixels, and C1<C2; Based on the minimum a posteriori energy Obtain ζ0 and obtain the high-frequency component in the horizontal direction of the image frequency band based on the ripple characteristics Water ripple area Formula (4) is used to obtain the high-frequency component in the vertical direction of the secondary decomposition: High frequency components in the diagonal direction The minimum value of the posterior energy is obtained to obtain the high-frequency component in the vertical direction of the image frequency band based on the ripple characteristics after the second-level decomposition. Water ripple area High frequency components in the diagonal direction Water ripple area 5. The method according to claim 4, characterized in that The specific steps of step (23) of "removing water ripples from the image frequency band obtained in step (22) based on the expectation maximization method" are as follows: Step (231): setting the initial value θ0 of the probability distribution parameter θ; Step (232): Desired step: Assume that the high-frequency component in the horizontal direction of the image frequency band based on the ripple characteristic Water ripple area There are N pixels in the n }, F represents the enhancement type corresponding to pixel N, recorded as {F1, F2, ..., F n }; Substitute the initial value θ0 in step (231) into the likelihood formula R(θ0) to express it as: Among them, w(N,F|θ) is the joint distribution probability model of variables N and F, which is specifically the high-frequency component in the horizontal direction Water ripple area The probability of a pixel belonging to class enhancement type F, The above formula is simplified to Calculate the current estimate of the enhancement type F; Step (233): Maximization step: Substitute the current estimate of F calculated in step (232) into the maximum likelihood estimation formula: Repeat the iteration steps (232) and (233) until the two steps converge at the same time, and maximize the likelihood function of formula (6) through iteration to obtain the maximum value of parameter θ1; then calculate the accurate value of the class enhancement type F for N pixels based on the probability model method of w(N|θ1); Step (234): The vertical high-frequency component water ripple region of the secondary decomposition and diagonal direction The enhancement type F is obtained and substituted into formula (6) to calculate the vertical high-frequency component of the secondary decomposition in the image frequency band based on the ripple characteristics: Diagonal high frequency components The accurate value of the regional enhancement type; get the estimated horizontal high-frequency component of the secondary decomposition without the influence of water ripples High frequency components in the vertical direction Diagonal high frequency components 6. The method according to claim 5, characterized in that Step (24) specifically comprises the following steps: Step (241): Calculate the global threshold T0 based on the equation maximization Based on the first-order Otsu method, the global function V(x) is maximized. The global function V(T0) is expressed as: in u ε It is the high-frequency component in the horizontal direction of the image frequency band based on the ripple characteristics Water ripple area The probability of a pixel with gray level ε in ; P ε is a probability density function based on grayscale, M(T0) is the probability that the pixel is within the threshold T0, λ(T0) is the probability that the pixel value is higher than the threshold T0, Represents the average value of pixels in the image; this is used to separate the high-frequency components in the horizontal direction of the image frequency band based on the ripple characteristics The shadow and non-shadow areas to be determined The vertical high-frequency components of the image frequency band that meet the ripple characteristics are decomposed by the second level High frequency components in the diagonal direction Water ripple area Substitute the probability of pixels with gray level ε in the image into formula (7) to obtain the high-frequency component in the vertical direction of the image frequency band based on the ripple characteristic after the second-level decomposition: High frequency components in the diagonal direction The shadow and non-shadow areas to be determined Step (242): determining the optimal threshold T1 based on minimization of regional variance; The high-frequency component L in the horizontal direction in the image frequency band based on the ripple characteristic obtained in step (241) 1 i H i 1 The shadow and non-shadow areas to be determined The optimal threshold is determined by maximizing the variance; the high-frequency component in the horizontal direction of the image frequency band based on the ripple characteristics The shadow and non-shadow areas to be determined The average intensity u w and regional variance It is expressed as: L represents the high-frequency component in the horizontal direction of the image frequency band based on the ripple characteristics The shadow and non-shadow areas to be determined The gray level, the global threshold T0 divides the gray level into pixels {1,2,...,l} and pixels {l+1,...,L}, γ k ,u k Represents the high-frequency component in the horizontal direction of the image frequency band based on the ripple characteristics The shadow and non-shadow areas to be determined The cumulative probability and average grayscale of Minimize the high-frequency components in the horizontal direction in this smaller image band based on the ripple characteristics The shadow and non-shadow areas to be determined The variance maximization method is used to obtain the high-frequency component in the vertical direction of the secondary decomposition in the image frequency band based on the ripple characteristics after the secondary decomposition. Diagonal high frequency components The shadow and non-shadow areas to be determined Step (243): iteratively eliminating the water ripple area; According to step (242), the high-frequency components in the horizontal direction of the image frequency band based on the ripple characteristic are respectively Vertical Diagonal direction The shadow and non-shadow areas to be determined The size of steps (22) and (23) is continuously iterated. When the shadow and non-shadow areas are smaller than a preset threshold condition ν, it indicates that the water ripple area has been eliminated. The absolute value of the preset threshold condition ν satisfies that the grayscale pixel value inside the shadow and non-shadow areas to be determined is not greater than 20% of the pixel value outside; It means that the water ripple area has been eliminated, and the high-frequency component in the horizontal direction of the second-level decomposition without the influence of water ripples is obtained. High frequency components in the vertical direction High frequency components in the diagonal direction 7. The method according to claim 6, characterized in that Step (3) performs inverse wavelet transform on the image as follows: The Mallat method using two-dimensional wavelet decomposition is: Among them, m0, n0 represent the coefficient matrix The rows and columns of For low frequency components and the high-frequency components in the horizontal direction High frequency components in the vertical direction High frequency components in the diagonal direction Wavelet decomposition to the highest decomposition order gives Multi-level wavelet transform; when it is known After these four decomposition coefficient matrices, they are reconstructed by the following formula Fusion and Get the image without water ripple enhancement Among them, T i c is the fusion coefficient, T i c ∈(0,1), and Assume that there is no water ripple enhanced image The fused background light is G1, and the image after water ripple enhancement is processed using a guided filter and substituted into the restoration formula: Among them, t c (a, b) are RGB channel transmission maps, is the background light of the RGB channel, This is an enhanced image after color balancing combined with background light.

Images