Digital watermarking algorithm based on generalized multivariate gamma distribution

Abstract

The invention discloses a digital watermarking method based on generalized multivariate gamma distribution, and the method comprises the steps: 1, carrying out the three-stage non-sampling dual-tree complex wavelet transformation of a host image, carrying out the subtraction of a real part coefficient and an imaginary part coefficient in the same direction at a third scale, calculating a difference value, and selecting a sub-band with the maximum energy as a target difference value sub-band, thirdly, carrying out third-order moment transformation on the target difference value sub-band by using an accurate pseudo Jacobi-Fourier moment to obtain a UDTCWT Difference + APJFMs amplitude domain as a watermark carrier, and carrying out third-order moment transformation on the target difference value sub-band by using the UDTCWT Difference + APJFMs amplitude domain as a watermark carrier; secondly, embedding the watermark into a selected area by utilizing a multiplicative embedding rule; thirdly, selecting a three-scale target difference value sub-band and a two-scale corresponding difference value sub-band for modeling, deriving corresponding parameters, and performing accurate estimation on shape and position parameters by adopting a low-complexity expectation maximization algorithm; and 4, designing a global optimal decoder according to an ML criterion, and extracting a watermark bit according to a decision threshold. Results of a large number of repeated experiments prove that the method has good decoding performance.

Description

Technical Field

[0001] The invention belongs to the technical field of image copyright protection and relates to a digital watermark algorithm based on generalized multivariate gamma distribution. Background Art

[0002] The advent of the third technological revolution has enabled network communications to penetrate into every aspect of life. The widespread dissemination of text, images, music, and pictures on the Internet has facilitated users' sharing and storage, but it has also given rise to illegal image editing and copyright infringement. Digital watermarking technology can protect the copyright of multimedia works on a large scale.

[0003] The performance evaluation indicators of digital watermarking algorithms include imperceptibility, robustness, and watermark capacity. Robustness can be defined as the ability of the watermarking scheme to extract and verify the watermark after the image is attacked; imperceptibility refers to the ability to maintain the perceptual characteristics of the host image after it is distorted by the addition of the watermark; and watermark capacity can be defined as the amount of secret information embedded in the image. All three indicators are important, but there are inherent contradictions between them that are not easy to resolve.

[0004] Compared with detectors, decoders are more widely used in watermark algorithms. The role of decoders is to accurately extract the watermark information hidden in the image. Although watermark schemes based on statistical models can effectively ensure the balance between imperceptibility, robustness and watermark capacity, the following problems still exist through the analysis of existing watermark algorithms based on statistical modeling: first, the robustness of the watermark carrier is insufficient; second, the traditional distribution model has weak description ability and the correlation characteristics between signals are not fully utilized; third, the parameter estimation method is inefficient, which affects the decoding accuracy of the watermark. Summary of the invention

[0005] The present invention aims to solve the above technical problems existing in the prior art and provides a digital watermarking algorithm based on generalized multivariate gamma distribution.

[0006] The technical solution of the present invention is: a digital watermark algorithm based on generalized multivariate gamma distribution, which is performed according to the following steps:

[0007] Convention: L represents the watermark length; L r1 and L r2 Represents the two low-frequency subbands of the real part obtained by UDTCWT decomposition, L i1 and L i2 represents the two low-frequency sub-bands of the imaginary part obtained by UDTCWT decomposition; H r1 ~H r18 and H i1 ~H i18They represent the high-frequency subbands of the real part and the imaginary part decomposed by UDTCWT respectively; M represents the host image; M′ represents the image after embedding the watermark; F d Represents the real high-frequency subband coefficient of each scale, S d represents the imaginary high frequency subband coefficient of each scale, D d Indicates the difference high-frequency subband of each scale calculated, d represents the direction, d = 1, 2...6; N represents the size of the image block; r, j represent the horizontal and vertical coordinates respectively; y i Represents the amplitude coefficient containing the watermark; x i represents the amplitude coefficient without watermark; w l is the information sequence that needs to hide the watermark; λ refers to the weighting factor; and Represent the watermark sequence variance and amplitude coefficient variance respectively; δ represents a custom minimum number; α is a shape parameter, β is a location parameter, Σ is a 4×4 covariance scale matrix; Γ(.) is a gamma function; tr is the rank of the matrix; represents the lth high entropy block l = 1, 2...L; U represents the binary quantity composed of the i-th coefficients of the sub-blocks with and without watermarks at the same position; H 0 , H 1 Refers to the watermark position that needs to be hidden; Z l represents the response value; τ l represents the decision threshold;

[0008] a. Initial Setup

[0009] Get the original image I and initialize the variables;

[0010] b. Watermark Embedding

[0011] b.1 Perform M three-level UDTCWT decomposition on the host image to obtain two low-frequency subbands L in the real part r1 , L r2 And the six high-frequency sub-bands H in different directions of scale 1, scale 2, and scale 3 r1 ~H r18 , imaginary part: 2 low frequency sub-bands L i1 , L i2 And the six high-frequency sub-bands H in different directions of scale 1, scale 2, and scale 3 i1 ~H i18 , where the resulting high-frequency subband size is the same as the host image size;

[0012] b.2 Calculate the three-scale one to six-direction UDTCWT difference subbands using the following formula:

[0013] D d (r,j)=F d (r,j)-Sd (r, j);

[0014] b.3 Select the subband with the maximum energy difference in scale 3. The specific calculation method is as follows:

[0015]

[0016] b.4 Divide the scale three energy maximum difference subband into N×N blocks, calculate the entropy value of each block, sort the entropy values ​​from large to small according to the embedded watermark size L, select the first L blocks and perform APJFMs transformation (n=3, |m|=4) to obtain UDTCWT-Difference+APJFMs amplitude blocks, and select the precise amplitude domain in each amplitude block as the watermark carrier;

[0017] b.5 Embed watermarks. Perform watermark information embedding operations according to the multiplicative embedding rule of watermarks:

[0018] y i =(1+λw l )x i

[0019]

[0020] b.6 Perform APJFMs reconstruction, modify the coefficients, put them back to their original positions, and perform inverse UDTCWT transformation to obtain the watermarked image M′;

[0021] c. Generalized Multivariate Gamma Distribution Modeling

[0022] c.1 Perform three-level UDTCWT decomposition on the watermarked image M′ to obtain two low-frequency subbands L in the real part r1 ',L r2 ' and 6 high-frequency subbands H in different directions at each scale r1 '~H r18 ', imaginary part 2 low frequency subbands L i1 ',L i2 ' and 6 high-frequency subbands H in different directions at each scale i1 '~H i18 ';

[0023] c.2 Select the subband with the maximum energy difference of the three scales, and select the corresponding two-scale difference subband. If the subband with the maximum energy difference of the three scales is in the first direction, then the two-scale difference subband should also select the first direction, divide the target subband into equal-sized non-overlapping blocks, perform APJFMs moment transformation on each frequency domain block to obtain the UDTCWT Difference-APJFMs amplitude, bind the UDTCWT Difference-APJFMs amplitude of each frequency domain block to form a vector, and then use the generalized multivariate gamma model to perform statistical modeling on these vectors; the probability density function f based on the generalized multivariate gamma model is expressed as follows:

[0024]

[0025] d. Low complexity expectation maximization algorithm

[0026] d.1 Select noisy and noise-free samples for parameter estimation, derive the pseudo-log likelihood function and select the initial value: give the initial estimate α of the shape parameter α and the location parameter β (0) , β (0) , where the parameter Σ is the covariance matrix of the selected training samples;

[0027] d.2 Loop through the E and M steps, given any small positive number δ, |α m+1 -α m |<δ,|β m+1 -β m |<δ, when both equations hold true, the unknown parameters α, β can be obtained;

[0028] e. Construct ML decoder to extract watermark

[0029] e.1 The extraction of digital watermark is regarded as a binary hypothesis testing process. Assume that H 1 With H 0 Respectively, they indicate that the information embedded in the current coefficient to be tested is "+1" and "-1", then the hypothesis test is defined as follows:

[0030] H 1 :y i =(1+λ)x i ,ω l =+1

[0031] H 0 :y i =(1-λ)x i ,ω l =-1

[0032] e.2 According to the Maximum Likelihood decision criterion, the designed likelihood ratio expression is as follows:

[0033]

[0034] e.3 The expression of conditional probability under the two assumptions Substituting into the above formula, the log-likelihood ratio is further calculated, f X (x) is the PDF of the selected subband coefficients, then the decoder based on generalized multivariate gamma is expressed as follows:

[0035]

[0036] e.4 According to the following decision expression, extract the lth watermark bit in the watermark image:

[0037]

[0038] e.5 According to the following expression, if Z l (y)≥τ l , then the watermark information is regarded as "+1"; otherwise, the extracted watermark is "-1";

[0039]

[0040] The present invention discloses a digital watermarking method based on generalized multivariate gamma distribution. First, a three-level unsampled dual-tree complex wavelet transform is performed on the host image, and the real and imaginary coefficients in the same direction under the third scale are subtracted. The subband with the largest variance is selected as the target difference subband after calculating the difference. Then, the target difference subband is transformed by the third-order moment using the accurate pseudo-Jacobi-Fourier moment to obtain the UDTCWT Difference+APJFMs amplitude domain as the watermark carrier; second, the watermark is embedded into the selected area using the multiplicative embedding rule; third, the three-scale target difference subband and the corresponding difference subband of the two scales are selected for modeling, and the corresponding parameters are derived, and the low-complexity expectation maximization algorithm is used for accurate estimation of shape and position parameters; fourth, according to the ML criterion, a global optimal decoder is designed, and the watermark position is extracted according to the decision threshold. The present invention has been repeatedly tested in large quantities, and the results show that it has good decoding performance.

[0041] Compared with the existing technology, the present invention has the following advantages:

[0042] First, a unique watermark carrier. The UDTCWT decomposition difference is combined with the APJFMs moment transform, and the former's multi-directional and displacement-invariant characteristics are combined with the latter's strong reconstruction ability and good description ability to design an optimal watermark carrier that is highly resistant to attack and imperceptible. At the same time, this carrier also improves the performance of the decoder;

[0043] Second, a generalized multivariate gamma distribution model is constructed. This model fully considers the edge distribution characteristics of the modeled object and the multiple correlations within the subband and between scales, and can accurately fit the amplitude distribution such as non-Gaussian and asymmetric constant positive data, thus improving the accuracy of watermark information extraction;

[0044] Third, accurate and efficient parameter estimation method: The present invention uses a low-complexity expectation maximization algorithm to estimate model parameters, which greatly reduces the time complexity and greatly improves the accuracy. BRIEF DESCRIPTION OF THE DRAWINGS

[0045] Figure 1 This is a non-Gaussian characteristic histogram of the watermark carrier verified in the embodiment of the present invention.

[0046] Figure 2 This is a graph of the PDF fitting results of the generalized multivariate gamma distribution according to an embodiment of the present invention.

[0047] Figure 3 This is a watermarked image in which a 256-bit watermark is embedded in the grayscale image according to an embodiment of the present invention.

[0048] Figure 4 The grayscale image of the embodiment of the present invention contains a 256-bit watermark image and a 20-fold relief effect image of the original image.

[0049] Figure 5 This is a line chart comparing the 1024-bit watermark of an embodiment of the present invention with other algorithms under geometric attacks.

[0050] Figure 6 The following is a flow chart of watermark embedding according to an embodiment of the present invention.

[0051] Figure 7 This is a flow chart of watermark extraction according to an embodiment of the present invention. DETAILED DESCRIPTION

[0052] The present invention provides a digital watermarking algorithm based on generalized multivariate gamma distribution, such as Figure 6 , 7 As shown, follow the steps below:

[0053] Convention: L represents the watermark length; L r1 and L r2 Represents the two low-frequency subbands of the real part obtained by UDTCWT decomposition, L i1 and L i2 represents the two low-frequency sub-bands of the imaginary part obtained by UDTCWT decomposition; H r1 ~H r18 and H i1 ~H i18 They represent the high-frequency subbands of the real part and the imaginary part decomposed by UDTCWT respectively; M represents the host image; M′ represents the image after embedding the watermark; Fd Represents the real high-frequency subband coefficient of each scale, S d represents the imaginary high frequency subband coefficient of each scale, D d Indicates the difference high-frequency subband of each scale calculated, d represents the direction, d = 1, 2...6; N represents the size of the image block; r, j represent the horizontal and vertical coordinates respectively; y i Represents the amplitude coefficient containing the watermark; x i represents the amplitude coefficient without watermark; w l is the information sequence that needs to hide the watermark; λ refers to the weighting factor; and Represent the watermark sequence variance and amplitude coefficient variance respectively; δ represents a custom minimum number; α is a shape parameter, β is a location parameter, Σ is a 4×4 covariance scale matrix; Γ(.) is a gamma function; tr is the rank of the matrix; represents the lth high entropy block l = 1, 2...L; U represents the binary quantity composed of the i-th coefficients of the sub-blocks with and without watermarks at the same position; H 0 , H 1 Refers to the watermark position that needs to be hidden; Z l represents the response value; τ l represents the decision threshold;

[0054] a. Initial Setup

[0055] Get the original image I and initialize the variables;

[0056] b. Watermark Embedding

[0057] b.1 Perform M three-level UDTCWT decomposition on the host image to obtain two low-frequency subbands L in the real part r1 , L r2 And the six high-frequency sub-bands H in different directions of scale 1, scale 2, and scale 3 r1 ~H r18 , imaginary part: 2 low frequency sub-bands L i1 , L i2 And the six high-frequency sub-bands H in different directions of scale 1, scale 2, and scale 3 i1 ~H i18 , where the resulting high-frequency subband size is the same as the host image size;

[0058] b.2 Calculate the three-scale one to six-direction UDTCWT difference subbands using the following formula:

[0059] D d (r,j)=F d (r,j)-S d (r, j);

[0060] b.3 Select the subband with the maximum energy difference in scale 3. The specific calculation method is as follows:

[0061]

[0062] b.4 Divide the scale three energy maximum difference subband into N×N blocks, calculate the entropy value of each block, sort the entropy values ​​from large to small according to the embedded watermark size L, select the first L blocks and perform APJFMs transformation (n=3, |m|=4) to obtain UDTCWT-Difference+APJFMs amplitude blocks, and select the precise amplitude domain in each amplitude block as the watermark carrier;

[0063] b.5 Embed watermarks. Perform watermark information embedding operations according to the multiplicative embedding rule of watermarks:

[0064] y i =(1+λw l )x i

[0065]

[0066] b.6 Perform APJFMs reconstruction, modify the coefficients, put them back to their original positions, and perform inverse UDTCWT transformation to obtain the watermarked image M′;

[0067] c. Generalized Multivariate Gamma Distribution Modeling

[0068] c.1 Perform three-level UDTCWT decomposition on the watermarked image M′ to obtain two low-frequency subbands L in the real part r1 ',L r2 ' and 6 high-frequency subbands H in different directions at each scale r1 '~H r18 ', imaginary part 2 low frequency subbands L i1 ',L i2 ' and 6 high-frequency subbands H in different directions at each scale i1 '~H i18 ';

[0069] c.2 Select the subband with the maximum energy difference of the three scales, and select the corresponding two-scale difference subband. If the subband with the maximum energy difference of the three scales is in the first direction, then the two-scale difference subband should also select the first direction, divide the target subband into equal-sized non-overlapping blocks, perform APJFMs moment transformation on each frequency domain block to obtain the UDTCWT Difference-APJFMs amplitude, bind the UDTCWT Difference-APJFMs amplitude of each frequency domain block to form a vector, and then use the generalized multivariate gamma model to perform statistical modeling on these vectors; the probability density function f based on the generalized multivariate gamma model is expressed as follows:

[0070]

[0071] d. Low complexity expectation maximization algorithm

[0072] d.1 Select noisy and noise-free samples for parameter estimation, derive the pseudo-log likelihood function and select the initial value: give the initial estimate α of the shape parameter α and the location parameter β (0) , β (0) , where the parameter Σ is the covariance matrix of the selected training samples;

[0073] d.2 Loop through the E and M steps, given any small positive number δ, |α m+1 -α m |<δ,|β m+1 -β m |<δ, when both equations hold true, the unknown parameters α, β can be obtained;

[0074] e. Construct ML decoder to extract watermark

[0075] e.1 The extraction of digital watermark is regarded as a binary hypothesis testing process. Assume that H 1 With H 0 Respectively, they indicate that the information embedded in the current coefficient to be tested is "+1" and "-1", then the hypothesis test is defined as follows:

[0076] H 1 :y i =(1+λ)x i ,ω l =+1

[0077] H 0 :y i =(1-λ)x i ,ω l =-1

[0078] e.2 According to the Maximum Likelihood decision criterion, the designed likelihood ratio expression is as follows:

[0079]

[0080] e.3 The expression of conditional probability under the two assumptions Substituting into the above formula, the log-likelihood ratio is further calculated, f X (x) is the PDF of the selected subband coefficients, then the decoder based on generalized multivariate gamma is expressed as follows:

[0081]

[0082] e.4 According to the following decision expression, extract the lth watermark bit in the watermark image:

[0083]

[0084] e.5 According to the following expression, if Z l (y)≥τ l , then the watermark information is regarded as "+1"; otherwise, the extracted watermark is "-1";

[0085]

[0086] Experimental test and parameter setting:

[0087] The experimental environment is MATLAB R2018a, the grayscale image size is 512×512, and the image database is selected: http: / / decsai.ugr.es / cvg / dbimagenes / index.php.

[0088] Figure 1 This is a non-Gaussian characteristic histogram for verifying the watermark carrier in the embodiment of the present invention.

[0089] Figure 2 This is a graph of the PDF fitting results of the generalized multivariate gamma distribution according to an embodiment of the present invention.

[0090] Figure 3 This is a watermarked image in which a 256-bit watermark is embedded in the grayscale image according to an embodiment of the present invention.

[0091] Figure 3 (a) Original Lena image; (b) Original Mandrill image; (c) Original Boat image; (d) Original Peppers image; (e) Lena image with 256 bits watermark; (f) Mandrilln image with 256 bits watermark; (g) Boat image with 256 bits watermark; (h) Peppers image with 256 bits watermark

[0092] Figure 4 This is a relief effect image (magnified 20 times) of the grayscale image in the embodiment of the present invention after embedding the 256-bit watermark and comparing it with the original image.

[0093] Figure 4 (a) Lena relief image; (b) Mandrill relief image; (c) Boat relief image; (d) Peppers relief image.

[0094] Figure 5 This is a line chart comparing the 1024-bit watermark under geometric attack in an embodiment of the present invention with other excellent algorithms.

[0095] Figure 5Among them, (a) Scaling; (b) Translation; (c) Cropping; (d) Rotation;

[0096] Figure 5 The comparative documents used: M Amini, M O Ahmad, M N S Swamy. A robust multibit multiplicative watermark decoder using vector-based hidden Markov model in wavelet domain. IEEE Transactions on Circuits & Systems for Video Technology, 2018, 28(2): 402-413

[0097] Xiang-yang WANG, Xin SHEN, Jia-lin TIAN, Pan-pan NIU, Hong-ying YANG. Statistical image watermark decoder using high-order difference coefficients and bounded generalized Gaussian mixtures-based HMT. Signal Processing, 2022, 192: 108371.

Claims

1. A digital watermarking algorithm based on generalized multivariate gamma distribution, characterized in that Follow these steps: Convention: L represents the watermark length; L r1 and L r2 Represents the two low-frequency subbands of the real part obtained by UDTCWT decomposition, L i1 and L i2 represents the two low-frequency sub-bands of the imaginary part obtained by UDTCWT decomposition; H r1 ~H r18 and H i1 ~H i18 They represent the high-frequency subbands of the real part and the imaginary part decomposed by UDTCWT respectively; M represents the host image; M′ represents the image after embedding the watermark; F d Represents the real high-frequency subband coefficient of each scale, S d represents the imaginary high frequency subband coefficient of each scale, D d Indicates the difference high-frequency subband of each scale calculated, d represents the direction, d = 1, 2...6; N represents the size of the image block; r, j represent the horizontal and vertical coordinates respectively; y i Represents the amplitude coefficient containing the watermark; x i represents the amplitude coefficient without watermark; w l is the information sequence that needs to hide the watermark; λ refers to the weighting factor; and Represent the watermark sequence variance and amplitude coefficient variance respectively; δ represents a custom minimum number; α is a shape parameter, β is a location parameter, Σ is a 4×4 covariance scale matrix; Γ(.) is a gamma function; tr is the rank of the matrix; represents the lth high entropy block l = 1, 2...L; U represents the binary quantity composed of the i-th coefficient of the sub-block with and without watermark in the same position; H0, H1 refer to the watermark position to be hidden; Z l represents the response value; τ l represents the decision threshold; a. Initial Setup Get the original image I and initialize the variables; b. Watermark Embedding b.1 Perform M three-level UDTCWT decomposition on the host image to obtain two low-frequency subbands L in the real part r1 , L r2 And the six high-frequency sub-bands H in different directions of scale 1, scale 2, and scale 3 r1 ~H r18 , imaginary part: 2 low frequency sub-bands L i1 , L i2 And the six high-frequency sub-bands H in different directions of scale 1, scale 2, and scale 3 i1 ~H i18 , where the resulting high-frequency subband size is the same as the host image size; b.2 Calculate the three-scale one to six-direction UDTCWT difference subbands using the following formula: D d (r,j)=F d (r,j)-S d (r,j); b.3 Select the subband with the maximum energy difference in scale 3. The specific calculation method is as follows: b.4 Divide the scale three energy maximum difference subband into N×N blocks, calculate the entropy value of each block, sort the entropy values ​​from large to small according to the embedded watermark size L, select the first L blocks and perform APJFMs transformation (n=3, |m|=4) to obtain UDTCWT-Difference+APJFMs amplitude blocks, and select the precise amplitude domain in each amplitude block as the watermark carrier; b.5 Embed watermarks. Perform watermark information embedding operations according to the multiplicative embedding rule of watermarks: y i =(1+λw l )x i b.6 Perform APJFMs reconstruction, modify the coefficients, put them back to their original positions, and perform inverse UDTCWT transformation to obtain the watermarked image M′; c. Generalized Multivariate Gamma Distribution Modeling c.1 Perform three-level UDTCWT decomposition on the watermarked image M′ to obtain two low-frequency subbands L in the real part r1 ',L r2 ' and 6 high-frequency subbands H in different directions at each scale r1 '~H r18 ', imaginary part 2 low frequency subbands L i1 ',L i2 ' and 6 high-frequency subbands H in different directions at each scale i1 '~H i18 '; c.2 Select the subband with the maximum energy difference of the three scales, and select the corresponding two-scale difference subband. If the subband with the maximum energy difference of the three scales is in the first direction, then the two-scale difference subband should also select the first direction, divide the target subband into equal-sized non-overlapping blocks, perform APJFMs moment transformation on each frequency domain block to obtain the UDTCWTDifference-APJFMs amplitude, bind the UDTCWT Difference-APJFMs amplitude of each frequency domain block to form a vector, and then use the generalized multivariate gamma model to perform statistical modeling on these vectors; the probability density function f based on the generalized multivariate gamma model is expressed as follows: d. Low complexity expectation maximization algorithm d.1 Select noisy and noise-free samples for parameter estimation, derive the pseudo-log likelihood function and select the initial value: give the initial estimate α of the shape parameter α and the location parameter β (0) , β (0) , where the parameter Σ is the covariance matrix of the selected training samples; d.2 Loop through the E and M steps, given any small positive number δ, |α m+1 -α m |<δ,|β m+1 -β m |<δ, when both equations hold true, the unknown parameters α, β can be obtained; e. Construct ML decoder to extract watermark e.1 The extraction of digital watermark is regarded as a binary hypothesis test process. Assuming that H1 and H0 represent the embedded information of the current coefficient to be tested as "+1" and "-1" respectively, the hypothesis test is defined as follows: H1:y i =(1+λ)x i ,oh l =+1 H0:y i =(1-λ)x i ,oh l =-1 e.2 According to the Maximum Likelihood decision criterion, the designed likelihood ratio expression is as follows: e.3 The expression of conditional probability under the two assumptions Substituting into the above formula, the log-likelihood ratio is further calculated, f X (x) is the PDF of the selected subband coefficients, then the decoder based on generalized multivariate gamma is expressed as follows: e.4 According to the following decision expression, extract the lth watermark bit in the watermark image: e.5 According to the following expression, if Z l (y)≥τ l , then the extracted watermark information is regarded as "+1"; otherwise, the extracted watermark is "-1";

Images