Watermarking method based on multi-scale sparse interest points in the fresnel domain

Abstract

This invention relates to the fields of digital image processing and information security technology, and particularly to a watermarking method based on Fresnel domain multi-scale sparse interest points, comprising the following steps: First, the host image is preprocessed using the Fresnel domain to extract the host subband singular value matrix; then, the watermark image is decomposed by wavelet decomposition and gradient calculation to extract multi-scale interest point features, which are then embedded into the host subband singular value matrix after Arnold scrambling and singular value decomposition to reconstruct the watermarked host image; watermark extraction and decryption are performed through preprocessing, difference operations, and inverse scrambling, and multi-scale interest point reconstruction and wavelet correction are achieved through row and column direction interpolation and cropping correction. This reduces watermark data redundancy while effectively completing the continuous structure under sparse anchor point constraints, and improves the stability and reliability of watermark recovery under complex attack conditions such as geometric distortion.

Description

Technical Field

[0001] This invention relates to the fields of digital image processing and information security technology, specifically to digital watermarking and copyright protection technology, and particularly to a watermarking method based on Fresnel domain multi-scale sparse interest points. Background Technology

[0002] With the rapid development of the Internet and social platforms, digital images have become one of the main carriers of information dissemination. However, images are easily tampered with, forged, or illegally disseminated, posing serious challenges to their ownership authentication, copyright protection, and content traceability. Digital watermarking technology, as an effective means of information hiding and authentication, has been widely researched and applied in related fields.

[0003] Currently, mainstream watermarking technologies are mainly divided into spatial domain, transform domain, and deep learning methods. Spatial domain watermarking embeds information by directly modifying pixel values, which has the advantages of simple implementation and large capacity, but its robustness to common image processing attacks (such as compression and noise) is poor. Transform domain watermarking (such as DCT and DWT domains) modulates coefficients in the frequency domain, resulting in better visual concealment and greater robustness to common signal processing attacks, but it is more sensitive to geometric distortions (such as rotation and scaling). Methods represented by deep learning achieve watermark embedding and extraction through end-to-end training, showing good performance under specific attacks, but it relies on a large amount of data and complex models, and its generalization ability is insufficient under strong geometric attacks and unknown scenarios. In addition, although derivative schemes such as optical watermarking improve the resistance to geometric distortion, the systems are complex and difficult to popularize.

[0004] To reduce data redundancy and the impact of embedding, embedding methods based on sparse feature points or key regions have been proposed. These methods embed only a small number of structural "anchor points" rather than complete watermarks, effectively reducing storage and transmission overhead. However, in the recovery stage, relying solely on discrete anchor points makes it difficult to reconstruct continuous structural details. Especially when subjected to geometric attacks, local occlusion, or noise interference, anchor points are easily lost or mismatched, leading to structural breaks and missing details in the recovery results. Existing recovery methods mostly rely on feature matching and interpolation, lacking effective modeling and robust completion mechanisms for the continuous structure between anchor points, resulting in a significant decrease in recovery quality in complex textures or non-uniform deformation regions.

[0005] In summary, existing technologies still face significant contradictions and technical bottlenecks in achieving stable completion and geometric distortion-resistant restoration of continuous structures under sparse anchor point constraints. On the one hand, full-image embedding introduces data redundancy and perceptible perturbations. On the other hand, while sparse feature point embedding can reduce redundancy, it is insufficient in covering continuous structures and limits the stability and accuracy of restoration under complex distortion conditions. Summary of the Invention

[0006] To address the aforementioned technical issues, this invention employs a watermarking method based on Fresnel domain multi-scale sparse interest points. This method reduces watermark data redundancy while effectively completing continuous structures under sparse anchor point constraints, and improves the stability and reliability of watermark recovery under complex attack conditions such as geometric distortion. Watermarking methods for Fresnel domain multi-scale sparse interest points include: S1, Fresnel domain preprocessing: Host image A is preprocessed to obtain the host subband singular value matrix. and ; S2. Multi-scale interest point extraction: The watermark image is decomposed by wavelet and the gradient magnitude and orientation angle of the image are calculated. After obtaining the spatial orientation code of each pixel through spatial orientation coding, the gradient magnitude is subjected to multi-scale local maximum detection with consistent direction and effective features are extracted. S3, Watermark Embedding and Reconstruction: Select the watermark features to be embedded at the target scale, perform Arnold scrambling and singular value decomposition to obtain the singular value diagonal matrix of the watermark features, and then sort them according to the watermark strength parameters. Embedded in the singular value matrix and In the process, the final watermarked host image is obtained through inverse transformation reconstruction. ; S4. Watermark Extraction and Decryption: Watermarked host image The singular value matrix of watermark recovery is obtained through preprocessing. The singular value matrix is ​​used to obtain the singular value estimation matrix corresponding to the watermark through difference operation. Inverse singular value decomposition is performed to reconstruct the watermark gradient magnitude at the target scale. The watermark gradient magnitude at the target scale is then subjected to Arnold inverse scrambling to obtain the decrypted gradient magnitude. S5. Multi-scale interest point reconstruction: For the The results of the scale maxima detection are used to construct an anchor point set. The anchor points are then segmented in both the row and column directions to form an interpolation interval. One-dimensional interpolation is performed on each pair of adjacent anchor points in the interpolation interval. Based on the orientation angle selection criterion, the gradient magnitude is pruned and corrected in both the row and column directions. While keeping the orientation angle unchanged, the gradient magnitude is restored to the horizontal and vertical components, thus obtaining a gradient component field with consistent direction and continuous details, i.e., subband coefficients. S6, Wavelet Reconstruction Correction of Multi-Scale Gradients: Wavelet reconstruction correction is performed on the gradient component field, i.e., the subband coefficients j, to obtain a new set of subband coefficients, i.e., multi-scale wavelet coefficients. These multi-scale wavelet coefficients are then combined to obtain the image. .

[0007] Furthermore, the preprocessing operations S1 and S4 include: at a given wavelength and transmission distance Under the given conditions, the complex-valued coefficient matrix is ​​obtained by mapping to the Fresnel domain through Fresnel transform, and then converted to real-valued form by real-preservation transform to obtain a real-valued coefficient matrix. The real-valued coefficient matrix is ​​then normalized and quantized to obtain a grayscale image, i.e., the quantized Fresnel domain image. The Fresnel domain image is then subjected to discrete wavelet transform to obtain four sub-bands: the low-frequency sub-band of the Fresnel domain image, the high-frequency sub-band in the horizontal direction of the Fresnel domain image, the high-frequency sub-band in the vertical direction of the Fresnel domain image, and the high-frequency sub-band in the diagonal direction of the Fresnel domain image. Singular value decomposition is then performed on the low-frequency sub-band and the high-frequency sub-band in the diagonal direction of the Fresnel domain image to obtain a singular value matrix.

[0008] Furthermore, the wavelet decomposition of the watermarked image and the calculation of the gradient magnitude and orientation angle of the image include: performing discrete wavelet multi-scale decomposition on the watermarked image to obtain the sub-band coefficients at each scale. ,in, This represents the horizontal information of high-frequency components. Represents the vertical information of high-frequency components. Indicates low-frequency components; in the first At each decomposition scale, the corresponding gradient magnitude is calculated using the horizontal and vertical gradients of the image. With direction angle The calculation formula is as follows: ; ; in, Represents the gradient in the horizontal direction. Represents the gradient in the vertical direction. Indicates the first Location at each decomposition scale The gradient magnitude of a pixel. Indicates the first Location at each decomposition scale The orientation angle of a pixel is the first Location at each decomposition scale The phase of a pixel; The spatial orientation encoding operation obtains the spatial orientation code for each pixel, including: the orientation angle. The image is mapped to a preset angle range, and the preset angle range is divided into multiple sectors. Based on the multiple sectors, a spatial orientation code is assigned to each pixel. For gradient magnitude Execution of multi-scale local maxima detection and extraction of effective features with consistent direction includes: for the first Location at each decomposition scale The spatial encoding of pixels involves selecting neighboring pixels along the corresponding direction and comparing them with the first... Location at each decomposition scale The gradient magnitude of a pixel and its neighboring pixels are used to determine whether they are local maxima. Pixels that are greater than their neighbors on both sides are marked as local maxima with the same direction, while the rest are suppressed. After suppressing non-maxima at all scales, the points that remain in each decomposition layer that are local maxima at multiple scales are the effective features.

[0009] Furthermore, the inverse transform reconstruction process yields the final watermarked host image. Including: by watermark strength parameters Embedded in the singular value matrix and Afterwards, the low-frequency subband with embedded watermark is obtained through the Inverse Singular Value Decomposition (ISVD) operation. and the high-frequency subband in the horizontal direction after embedding watermark The low-frequency subband and high-frequency subband Together with the remaining subbands not involved in embedding, an inverse wavelet transform is performed to obtain the Fresnel domain real-valued coefficient matrix containing the watermark. This Fresnel domain real-valued coefficient matrix undergoes inverse quantization and is then remapped back to the spatial domain via an inverse Fresnel transform to obtain the final watermarked host image. .

[0010] Furthermore, obtaining the singular value estimation matrix corresponding to the watermark through difference operations on the singular value matrix includes: using the singular value matrix... and and watermark strength parameters The estimated value of the watermark component is recovered in the singular value domain; the calculation formula is as follows: ; ; in, This indicates that after preprocessing and Discrete Wavelet Transform (DWT) decomposition, the low-frequency subband of the watermarked image is... The singular value matrix obtained by performing singular value decomposition (SVD); This indicates that after preprocessing and Discrete Wavelet Transform (DWT) decomposition, the high-frequency subband of the watermarked image is... The singular value matrix obtained by performing singular value decomposition (SVD); This represents the singular value estimation matrix corresponding to the horizontal features at the second scale of the watermarked image; This represents the singular value estimation matrix corresponding to the vertical features at the second scale of the watermarked image.

[0011] Furthermore, regarding the first... The results of the scale maxima detection are used to construct an anchor point set, and the anchor points are segmented in both the row and column directions, including: Row direction: For fixed rows Collect the fixed rows The column indexes of all anchor points are used, and sorted from left to right to obtain the anchor point sequence; ; in, For the first Line up Anchor point column index, This represents the total number of anchor points in that row. Indicates the first The set of anchor indexes for rows. , , ..., These represent the column indexes of the first anchor point, the second anchor point, and the third anchor point, respectively. Column indexes for anchor points; Column direction: For fixed columns Collect the fixed column The row indices of all anchor points are sorted from top to bottom to obtain: ; in, Indicates the first List of the first Row index of anchor points This represents the total number of anchor points in this column. Indicates the first Anchor index set of columns, , , ..., These represent the row indexes of the first anchor point, the second anchor point, and the third anchor point, respectively. Row index of anchor points.

[0012] Furthermore, the interpolation interval performs one-dimensional interpolation for each pair of adjacent anchor points, including: Row-wise interpolation: S5.1: Each pair of adjacent anchor points naturally forms an interpolation interval, for the first... Each scale, based on anchor point segmentation, is used for fixed rows. Adjacent anchor points An interpolation interval is defined. Based on the current benchmark estimate The residuals at both ends are calculated using the following formula: ; ; in, Indicates the first In terms of scale, Gradient magnitude at position, Indicates in The baseline estimate of the location, Indicates in Positional residuals Indicates the first In terms of scale, Gradient magnitude at position, express The baseline estimate of the location, express Positional residuals; S5.2: Introduce a scale-dependent unit step size decay factor Defined as: ; in, Indicates scale. This represents the baseline attenuation coefficient, i.e., the value selected in the experiment; S5.3: Utilize the residuals described at both ends With residual and attenuation factor Constructing a two-ended exponential interpolation, we obtain the intra-segment correction term. Its analytical form can be written as a weighted combination of the exponential terms at both ends, as shown in the following expression: ; in, Indicates the current position. This represents the scale-dependent unit step size decay factor. S5.4: Update the current position The gradient strength is calculated using the following formula: ; in, This represents the magnitude of the row gradient after interpolation update. Indicates the current benchmark estimate; Column interpolation: S5.5: For fixed columns Adjacent anchor point pairs An interpolation interval is defined. Based on the current benchmark estimate The residuals at both ends are calculated using the following formula: ; ; in, Indicates the first In terms of scale, Gradient magnitude at position, Indicates in The baseline estimate of the location, Indicates in Positional residuals Indicates the first In terms of scale, Gradient magnitude at position, express The baseline estimate of the location, express Positional residuals; S5.6: Introduce a scale-dependent unit step size decay factor Defined as: ; in, Indicates scale. This represents the baseline attenuation coefficient, the value selected in the experiment; S5.7: Utilize the residuals described at both ends With residual and unit step size decay factor Constructing a two-ended exponential interpolation, we obtain the intra-segment correction term. Its analytical form can be written as a weighted combination of the exponential terms at both ends, as shown in the following expression: ; S5.8: Update the current position The gradient strength is calculated using the following formula: ; in, This represents the magnitude of the column gradient after interpolation.

[0013] Furthermore, for boundary intervals with anchor points on only one side, a single-ended exponential extension is used: the residuals on the side closest to the anchor point are attenuated along the segment by a unit step size decay factor. Attenuation propagation: When the number of anchor points in a row or column is less than 2, interpolation updates are not performed on the corresponding row or column.

[0014] Furthermore, the gradient magnitude is pruned and corrected in both row and column directions, including: S5.9: Select two constant thresholds, calculate the degree of deviation between the direction angle and the vertical reference direction, determine whether the corresponding direction is horizontal or vertical, and achieve precise orientation correction; The criterion for determining whether the corresponding direction is horizontal or vertical is: when When the corresponding pixel is considered to be moving away from the vertical direction and closer to the horizontal direction, line clipping is triggered; when When the direction of the corresponding pixel is considered to be close to vertical, column clipping is triggered; where Indicates the first Location at each decomposition scale The direction angle, , This represents the two selected constant thresholds; S5.10: Row-wise clipping: For fixed rows A segment, i.e., the interpolation interval The gradient magnitude sequence within a segment is denoted as: in, This represents the sequence of row gradient magnitudes after interpolation update; The minimum gradient magnitudes and positions of the gradient magnitude sequence within the segment are as follows: ; in, Indicates the first Line number The minimum gradient magnitude within the column interpolation interval, This indicates all pixels within the interval. gradient magnitude sequence Take the minimum value. This represents the minimum value. Indicates the phase at the minimum value; When the "horizontal" direction threshold is met, perform two monotonic scans along the row direction: Left sweep: from the left end position Towards the minimum position Scan, location ,like ; in, Let represent the gradient magnitude at the left end of the row direction, then let ; This represents the gradient magnitude in the direction following the left sweep. Right sweep: from the right end Towards Scan: Location ,like ; Let represent the gradient magnitude at the right end of the row direction, then let ; in, This represents the gradient magnitude in the direction following the right sweep. S5.11: Column-wise clipping: For fixed rows A segment, i.e., the interpolation interval The gradient magnitude sequence within a segment is denoted as: ; in, This represents the sequence of column gradient magnitudes after interpolation update; The minimum gradient magnitudes and positions of the gradient magnitude sequence within the segment are as follows: ; in, Indicates the first Line number The minimum gradient magnitude within the column interpolation interval, This indicates all pixels within the interval. gradient magnitude Take the minimum value. This represents the minimum value. Indicates the phase at the minimum value; When the "vertical" direction threshold is met, perform two monotonic scans along the column direction: Left sweep: from the left end position Towards the minimum position Scan, location ,like ; in, Let represent the gradient magnitude at the left end of the column direction, then let ; in, This represents the gradient magnitude in the column direction after the left sweep. Right sweep: from the right end Towards Scan: Location ,like ; in, Let the gradient magnitude at the right end of the column direction be represented; then let ; in, This represents the gradient magnitude in the column direction after the right sweep.

[0015] Furthermore, wavelet reconstruction correction is performed on the gradient component field, i.e., the subband coefficients, to obtain a new set of subband coefficients, i.e., multi-scale wavelet coefficients. These multi-scale wavelet coefficients are then combined to obtain the reconstructed image. Includes: gradient component fields, i.e., sub-band coefficients, based on S5. Perform inverse discrete wavelet transform to reconstruct the intermediate image. ;in, This represents the high-frequency horizontal components that maintain a consistent direction and continuous detail after cropping and trimming. This represents the high-frequency components in the vertical direction that maintain a consistent direction and continuous details after cropping and trimming. This represents low-frequency component information that maintains consistent direction and continuous detail after cropping and correction; the intermediate image. The grayscale range is linearly normalized, and a new set of multi-scale wavelet coefficients is obtained through discrete wavelet transform. ,in, This represents the updated high-frequency components in the horizontal direction. This represents the high-frequency components in the vertical direction after the update. This indicates the updated low-frequency component; Updated low-frequency components Replace with the original low-frequency component obtained in S2 This allows for updating only the high-frequency gradient subbands at each scale while maintaining consistency between the low-frequency structure and the initial image. and Finally, the original low-frequency components... With the updated high-frequency components , By combining the results, the final reconstructed image is obtained. .

[0016] Compared with the prior art, the advantages of the present invention are as follows: 1. Reduce watermark representation redundancy and system overhead. This invention no longer directly embeds the complete watermark image. Instead, based on discrete wavelet multi-scale decomposition and directional nonmaximum suppression, it extracts only structurally stable multi-scale sparse interest points and their gradient features as watermark information, and performs modulation embedding through singular value domain. This significantly reduces the amount of embedded data while retaining key structural information, thereby reducing storage, transmission and computational overhead. 2. Improve adaptability and robustness to geometric distortion. This invention unifies watermark embedding and extraction within the Fresnel domain. The Fresnel domain representation can more stably characterize the structural features after propagation / transformation. Combined with wavelet domain and singular value embedding mechanism, it helps to improve the stability and extractability of the method under geometric distortion conditions such as rotation and scaling, as well as conventional signal processing attacks. 3. Achieving effective completion of continuous structures under sparse anchor point constraints. Addressing the problem that existing sparse embedding methods rely solely on discrete feature points and struggle to recover continuous structures between anchor points, this invention proposes a scale-adaptive exponential directional interpolation strategy. This strategy propagates bidirectionally along the main direction of the interest point and smoothly merges within segments, effectively completing continuous details between anchor points without destroying the main direction information, thus improving the completeness of the recovery results. 4. Suppress interpolation artifacts and enhance reconstruction stability. To address overshoot, oscillation, and secondary peak artifacts that may be introduced by interpolation completion, this invention proposes an amplitude clipping and correction algorithm with directional consistency. While maintaining the direction angle unchanged, the amplitude within the segment is subject to monotonic constraints and neighborhood correction, thereby suppressing local anomalous abrupt changes and improving the stability and reliability of the reconstruction process. 5. Ensuring structural consistency and improving reconstruction quality. This invention further introduces a consistency correction mechanism based on wavelet reconstruction, which updates the high-frequency gradient sub-band while maintaining the stability of the low-frequency structure. This ensures that the completed multi-scale gradient meets the constraints that can be generated from the real image, reducing structural breaks and artifacts caused by inconsistencies and improving the quality of watermark restoration and reconstruction. Attached Figure Description

[0017] To more clearly illustrate the technical solutions of the embodiments of the present invention, the drawings used in the following description of the embodiments will be briefly introduced. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort. Figure 1 This is a flowchart of the watermark embedding and extraction process of the present invention, wherein... Figure 1 (a) represents the watermark embedding process. Figure 1 (b) indicates the watermark extraction process; Figure 2 This is a flowchart of the multi-scale interest point extraction and reconstruction process of the present invention; Figure 3 This is a diagram illustrating the spatial orientation encoding and directional gradient maxima analysis of the present invention. Figure 4 This is a diagram illustrating the interpolation technology architecture of the present invention. Figure 5 This is a diagram illustrating the cutting technology architecture of the present invention; Figure 6 This is a multi-scale point of interest distribution map of the present invention, wherein (a1) represents the original image; (a2) represents the low-frequency approximation component; (b1)–(c1) represent the horizontal high-frequency details at different scales; and (b2)–(c2) represent the vertical high-frequency details at different scales. Figure 7 This is a multi-scale interest point extraction and reconstruction result image of a 256×256 image of the present invention, where (a1)-(e1) represent the original image, (a2)-(e2) represent the interest point extraction result, the green area represents the location of the extracted interest points, and (a3)-(e3) represent the reconstruction result using the extracted interest points; Figure 8This is a multi-scale interest point extraction and reconstruction result of a 512×512 image of the present invention, where (a1)-(e1) represent the original image, (a2)-(e2) represent the interest point extraction result, the green area represents the location of the extracted interest points, and (a3)-(e3) represent the reconstruction result using the extracted interest points; Figure 9 The images shown are reconstruction results of the images under different interest point retention ratios according to the present invention. (a1)-(f1) represent the interest point distribution when 10%, 30%, 50%, 70%, 90%, and 100% of interest points are retained, respectively; (a2)-(f2) represent the corresponding reconstruction watermark results, where green represents all detected interest points and red represents interest points randomly selected under the retention ratio. Figure 10 The host diagram used in this invention is shown in which (a) represents a woman, (b) represents a baboon, (c) represents a house, (d) represents a sailboat, (e) represents a man, and (f) represents an airplane. Figure 11 The watermark image used in this invention is shown in the image below, where (a) represents a rose, (b) represents text, (c) represents a teacher, and (d) represents a cat. Figure 12 This is a comparison of the grayscale histograms of six host images before and after watermark embedding according to the present invention; Figure 13 The average normalized correlation coefficient value under signal processing attack of the present invention is given by (a) representing salt-and-pepper noise; (b) representing speckle noise; (c) representing Gaussian noise; (d) representing median filtering; (e) representing Gaussian low-pass filtering; (f) representing mean filtering; (g) representing JPEG compression; and (h) representing JPEG2000 compression. Figure 14 Here is the normalized correlation coefficient value of the present invention under geometric attack, where (a) is the rotation angle; and (b) represents the clipping ratio. Figure 15 This is the watermark image extracted under the combined attack of the present invention. Detailed Implementation

[0018] The present invention will be further described below with reference to the accompanying drawings and specific embodiments.

[0019] Example 1: like Figures 1-5 As shown, this embodiment provides a watermarking method based on Fresnel domain multi-scale sparse interest points, specifically including the following steps: S1, Fresnel domain preprocessing: Host image A is preprocessed to obtain the host subband singular value matrix. and .

[0020] As a specific implementation of the present invention: host image A at a given wavelength and transmission distance Under the given conditions, the complex-valued coefficient matrix B of the host image A is obtained by mapping to the Fresnel domain through Fresnel transform. Since the complex-valued coefficient matrix B is in complex form, to facilitate subsequent wavelet and singular value analysis, this invention uses a real-valued transform to convert it to real form, obtaining the real-valued coefficient matrix C of the host image A. The real-valued coefficient matrix C is discretized to the grayscale range of [0, 255] through normalization, and the grayscale image, i.e., the quantized Fresnel domain image H, is obtained through a nearest-neighbor quantization operation. The Fresnel domain image H is then transformed into four sub-bands through a single-layer two-dimensional discrete wavelet transform. These are the low-frequency subband of the Fresnel domain image H, the high-frequency subband in the horizontal direction of the Fresnel domain image H, the high-frequency subband in the vertical direction of the Fresnel domain image H, and the high-frequency subband in the diagonal direction of the Fresnel domain image H, respectively. To improve the stability of watermark embedding, the low-frequency subband of the Fresnel domain image H is... High-frequency subbands diagonally opposite to the Fresnel domain image H Singular value decomposition yields the singular value matrix of the Fresnel domain image H. and As the main carrier for subsequent watermark embedding; the expressions for two-dimensional discrete wavelet transform and singular value decomposition are as follows: ; ; ; in, Represents Fresnel domain image The low-frequency subband, Represents Fresnel domain image Horizontal high-frequency sub-band, Represents Fresnel domain image Vertical high-frequency sub-bands Represents Fresnel domain image High-frequency sub-bands in the diagonal direction; Indicates the low-frequency subband The left singular vector matrix obtained by performing singular value decomposition; Indicates the low-frequency subband The singular value matrix obtained by performing singular value decomposition; Indicates the low-frequency subband The right singular vector matrix obtained by performing singular value decomposition; Indicates the high-frequency subband The left singular vector matrix obtained by performing singular value decomposition; Indicates the high-frequency subband The singular value matrix obtained by performing singular value decomposition; Indicates the high-frequency subband The right singular vector matrix obtained by performing singular value decomposition.

[0021] In scalar diffraction theory, Fresnel diffraction is used to describe the incident light field, i.e., the host image. Propagation distance in free space Then, observe the plane coordinates. Complex amplitude distribution formed at the location The propagated light field can be expressed by the Fresnel diffraction integral as follows: ; in, It can be regarded as the input field A first integral transformation, representing plane coordinates The complex amplitude distribution formed at that point is the complex-valued coefficient matrix. This represents the two-dimensional coordinate variables on the input plane. Represents two-dimensional coordinate variables on the observation plane; This represents the amplitude decay and phase change factor of the entire integral. This represents the quadratic phase factor on the observation plane. This represents the quadratic phase factor on the input plane. The cross-phase factor represents the coordinate relationship between the input plane and the view plane; The Niel diffraction integral used in this invention has a kernel function that includes the propagation distance. With wavelength Determined quadratic phase term Therefore, the free-space propagation process is equivalent to the propagation of the input field. Applying a Fresnel transform yields a Fresnel domain representation.

[0022] S2. Multi-scale interest point extraction: The watermark image is decomposed by wavelet and the gradient magnitude and orientation angle of the image are calculated. After obtaining the spatial orientation code of each pixel through spatial orientation coding, multi-scale local maxima detection with consistent orientation is performed on the gradient magnitude and effective features are extracted.

[0023] The process of performing wavelet decomposition on a watermarked image and calculating its gradient magnitude and orientation angle involves: performing discrete wavelet multi-scale decomposition on the watermarked image to obtain the sub-band coefficients at each scale. ,in, Indicates the first High-frequency components in the horizontal direction at each decomposition scale Indicates the first High-frequency components in the vertical direction at each decomposition scale Indicates the first Low-frequency components at each decomposition scale.

[0024] In this embodiment, the input for discrete wavelet multi-scale decomposition of the watermark image is a two-dimensional watermark image. ,in, , , These are the height and width of the image; the discrete wavelet multi-scale decomposition process is represented as follows: ; Among them, the filter The first Low-pass decomposition filter and high-pass decomposition filter at each decomposition scale; Indicates the image at the 1st Low-frequency subbands in the horizontal direction of each decomposition scale. Indicates the first High-frequency components in the horizontal direction at each decomposition scale Indicates the first High-frequency components in the vertical direction at each decomposition scale Indicates the first Low-frequency components at each decomposition scale.

[0025] In the At each decomposition scale, the corresponding gradient magnitude is calculated using the horizontal and vertical gradients of the image. With direction angle The calculation formula is as follows: ; ; in, This represents the horizontal gradient of the watermarked image. This represents the vertical gradient of the watermark image. Indicates the first Location at each decomposition scale The gradient magnitude of a pixel. Indicates the first Location at each decomposition scale The orientation angle of the pixel, i.e., the first... Location at each decomposition scale The phase of a pixel.

[0026] The spatial orientation encoding operation obtains the spatial orientation code for each pixel as follows: Since the orientation angle is directly calculated... The range of values ​​is Therefore, the negative angle needs to be added Mapping to a preset angle range in a certain way and preset angle range Divided into 8 equal sectors, each sector having a width of [missing information]. Based on the eight sectors to which the orientation angle of a pixel belongs, a spatial orientation code of 1-8 is assigned to each pixel; the continuous gradient orientation angle is transformed into a finite number of discrete directions, providing a directional prior for subsequent local maximum detection along a given direction.

[0027] For gradient magnitude Execution of multi-scale local maxima detection and extraction of effective features with consistent direction specifically includes: for the first Location at each decomposition scale Spatial encoding of pixels involves selecting neighboring pixels along the corresponding direction and comparing the first pixel with the second pixel. Location at each decomposition scale The gradient magnitude of a pixel and its neighboring pixels is used to determine whether it is a local maximum. Pixels that are greater than their two neighboring pixels are marked as local maxima with the same direction, while the rest are suppressed. After suppressing non-maximum values ​​at all scales, the points that remain as local maxima at multiple scales in each decomposition layer are the effective features.

[0028] The specific rules for determining local maxima are as follows: Horizontal direction (code 1 or 5): If the gradient magnitude of the current pixel If the gradient magnitude of a pixel is greater than that of its left and right neighbors, then the pixel is considered a local maximum. ; in, This represents the gradient magnitude of the left neighbor in the horizontal direction. This represents the gradient magnitude of the right neighbor in the horizontal direction; Diagonal direction (code 2 or 6): If the gradient magnitude of the current pixel If the gradient magnitude is greater than both the top-left and bottom-right neighbors, it is determined to be a local maximum, i.e.: ; in, This represents the gradient magnitude of the top-left neighbor in the diagonal direction. This represents the gradient magnitude of the lower right neighbor in the diagonal direction; Vertical direction (code 3 or 7): If the gradient magnitude of the current pixel If the gradient magnitude is greater than that of its upper and lower neighbors, it is determined to be a local maximum, that is: ; in, This represents the gradient magnitude of the neighboring units in the vertical direction. This represents the gradient magnitude of the lower neighbor in the diagonal direction; Diagonal direction (code 4 or 8): If the gradient magnitude of the current pixel If the gradient magnitude is greater than both the lower-left and upper-right neighbors, it is determined to be a local maximum, i.e.: ; in, This represents the gradient magnitude of the lower left neighbor in the vertical direction. This represents the gradient magnitude of the top-right neighbor in the diagonal direction.

[0029] S3, Watermark Embedding and Reconstruction: Select the watermark features to be embedded at the target scale, perform Arnold scrambling and singular value decomposition to obtain the singular value diagonal matrix of the watermark features, and then sort them according to the watermark strength parameters. Embedded in the singular value matrix With singular value matrix In the process, the final watermarked host image is obtained through inverse transformation reconstruction. ; The watermark embedding and reconstruction process is as follows: S3.1: Select the horizontal details of the high-frequency components at the second scale. Vertical details with the second-scale high-frequency components As a watermark feature to be embedded.

[0030] S3.2: To enhance the security of the watermark, Arnold scrambling is performed, and the spatial distribution of the pixels is shuffled through iterative mapping of the pixel positions to obtain the horizontal details of the high-frequency components at the second scale. Rear matrix Vertical details of two-scale high-frequency components Scrambled matrix Meanwhile, to improve the numerical stability of the embedding process and its robustness to attacks, the matrix... With matrix Perform singular value decomposition to obtain the singular value diagonal matrix. and singular value diagonal matrix The expression for singular value decomposition is as follows: ; ; in, This represents the second-scale horizontal (row direction) watermark feature map after Arnold scrambling; This represents the second-scale vertical (column) watermark feature map after Arnold scrambling. , , These represent the watermark feature maps respectively. The transpose of the left singular vector matrix, the singular value diagonal matrix, and the right singular vector matrix obtained by singular value decomposition; , , These represent the watermark feature maps respectively. The transpose of the left singular vector matrix, the singular value diagonal matrix, and the right singular vector matrix obtained by singular value decomposition.

[0031] S3.3: Based on the preset watermark strength parameters Singular value diagonal matrix and singular value diagonal matrix Embedded into singular value matrix With singular value matrix In Chinese, the expression is as follows: ; ; in, Represents the singular value matrix in the low-frequency subband of the host. Embedded singular value diagonal matrix The resulting updated singular value matrix; Represents the singular value matrix in the host high-frequency subband. Embedded singular value diagonal matrix The resulting updated singular value matrix.

[0032] The inverse transform reconstruction process yields the final watermarked host image. Specifically, based on watermark strength parameters: Embedded in the singular value matrix and The low-frequency subband with embedded watermark is obtained by inverse singular value decomposition (ISVD) in the middle and later stages. and the high-frequency subband in the horizontal direction after embedding watermark The low-frequency subband and horizontal high-frequency subband Together with the remaining subbands not involved in embedding, an inverse wavelet transform (IDWT) is performed to obtain the watermarked Fresnel domain real-valued coefficient matrix. This Fresnel domain real-valued coefficient matrix is ​​then inversely quantized and remapped back to the spatial domain using an inverse Fresnel transform to obtain the final watermarked host image. .

[0033] S4. Watermark Extraction and Decryption: Watermarked host image The singular value matrix for watermark recovery is obtained through preprocessing. The singular value matrix is ​​then used to obtain the singular value estimation matrix corresponding to the watermark through difference operations. Inverse singular value decomposition is then performed to reconstruct the watermark gradient magnitude at the target scale. The watermark gradient magnitude at the corresponding scale is then subjected to Arnold inverse scrambling to obtain the decrypted gradient magnitude.

[0034] As a specific implementation of the present invention: S4.1: Watermarked Host Image At a given wavelength and transmission distance Under the given conditions, the host image containing the watermark is obtained by mapping to the Fresnel domain through Fresnel transformation. Complex coefficient matrix Due to the complex-valued coefficient matrix Since the value is in complex form, a real-valued transformation is used to convert it to real form, resulting in the watermarked host image. Real-valued coefficient matrix Real-valued coefficient matrix The values ​​are discretized to the grayscale range of [0, 255] by normalization, and then the grayscale image, i.e., the quantized watermarked host image, is obtained by rounding down the nearest integer. Fresnel domain image Fresnel domain image Four sub-bands are obtained through a single-level two-dimensional discrete wavelet transform. These are Fresnel domain images. Low-frequency subband, Fresnel domain image Horizontal high-frequency subband and Fresnel domain images Vertical high-frequency subband and Fresnel domain images High-frequency subbands in the diagonal direction; to improve the stability of watermark embedding, Fresnel domain images are processed. low-frequency subband With Fresnel domain image High-frequency subbands in diagonal directions Singular value decomposition is performed to obtain the singular value matrix for subsequent watermark recovery. and singular value matrix The expressions for the two-dimensional discrete wavelet transform and singular value decomposition are as follows: ; ; ; in, Represents Fresnel domain image The low-frequency subband, Represents Fresnel domain image Horizontal high-frequency sub-band, Represents Fresnel domain image Vertical high-frequency sub-bands Represents Fresnel domain image High-frequency sub-bands in the diagonal direction; Indicates the low-frequency subband The left singular vector matrix obtained by performing singular value decomposition; Indicates the low-frequency subband The singular value matrix obtained by performing singular value decomposition; Indicates the low-frequency subband The right singular vector matrix obtained by performing singular value decomposition; Indicates the high-frequency subband The left singular vector matrix obtained by performing singular value decomposition; Indicates the high-frequency subband The singular value matrix obtained by performing singular value decomposition; Indicates the high-frequency subband The right singular vector matrix obtained by performing singular value decomposition.

[0035] S4.2: Singular Value Matrix and singular value matrix The singular value estimation matrix corresponding to the watermark is obtained through difference operations. and and embedding strength parameters The estimated value of the watermark component is obtained by recovering it in the singular value domain; the calculation formula is as follows: ; ; in, This indicates that after preprocessing and Discrete Wavelet Transform (DWT) decomposition, the low-frequency subband of the watermarked image is... The singular value matrix obtained by performing singular value decomposition (SVD); This indicates that after preprocessing and Discrete Wavelet Transform (DWT) decomposition, the high-frequency subband of the watermarked image is... The singular value matrix obtained by performing singular value decomposition (SVD); This represents the singular value estimation matrix corresponding to the horizontal features at the second scale of the watermarked image; This represents the singular value estimation matrix corresponding to the vertical features at the second scale of the watermarked image.

[0036] S4.3: Obtaining the singular value estimation matrix and singular value estimation matrix Then, it is combined with the left and right singular vector matrices saved during the embedding stage. , , ,and The watermark gradient magnitude at the second scale is reconstructed using inverse singular value decomposition (ISVD): ; ; in, This represents the high-frequency horizontal information of the reconstructed second scale, i.e., the horizontal gradient magnitude of the watermark. This represents the high-frequency vertical information of the reconstructed second scale, i.e., the vertical gradient magnitude of the watermark.

[0037] S4.4: Adjusting the horizontal gradient magnitude of the watermark Vertical gradient magnitude of watermark Perform Arnold inverse scrambling separately to restore the original spatial distribution from the scrambled coordinate system, and obtain the decrypted gradient magnitude.

[0038] S5. Multi-scale interest point reconstruction: For the The results of the scale maxima detection are used to construct an anchor point set. The anchor points are then segmented in both the row and column directions to form an interpolation interval. One-dimensional interpolation is performed on each pair of adjacent anchor points in the interpolation interval. Based on the orientation angle selection criterion, the gradient magnitude is clipped and corrected in both the row and column directions. While keeping the orientation angle unchanged, the gradient magnitude is restored to the horizontal and vertical components, thus obtaining a multi-scale gradient field with consistent orientation and continuous details.

[0039] Since maximum detection only retains non-zero detail coefficients at local maxima and sets them to zero at other locations, the gradient field between anchor points is relatively sparse. To restore continuous details without destroying the principal direction of the edges, this invention... The results of the scale maxima detection are used to construct an anchor point set, and the anchor points are segmented in both the row and column directions as follows: Row direction: For fixed rows Collect the fixed rows The column indexes of all anchor points are used, and sorted from left to right to obtain the anchor point sequence; ; in, For the first Line up Anchor column index, This represents the total number of anchor points in that row. Indicates the first The set of anchor indexes for rows; , , ..., These represent the column indexes of the first anchor point, the second anchor point, and the third anchor point, respectively. Column indexes for anchor points; Column direction: For fixed columns Collect the fixed column The row indices of all anchor points are sorted from top to bottom to obtain: ; in, Indicates the first List of the first Row index of anchor points This represents the total number of anchor points in this column. Indicates the first The set of anchor indexes for rows; , , ..., These represent the row indexes of the first anchor point, the second anchor point, and the third anchor point, respectively. Row index of anchor points.

[0040] After local maxima detection, non-zero gradient intensities are retained only at local maxima locations, while those at other locations are set to zero, thus creating numerous "empty segments" between anchor points. To restore continuous details without disrupting the main edge direction, this invention performs one-dimensional interpolation on each pair of adjacent anchor points in both row and column directions, specifically as follows: Row-wise interpolation: S5.1: Each pair of adjacent anchor points naturally forms an interpolation interval, for the first... Each scale, based on anchor point segmentation, is used for fixed rows. Adjacent anchor points An interpolation interval is defined. Based on the current benchmark estimate The residuals at both ends are calculated using the following formula: ; ; in, Indicates the first In terms of scale, Gradient magnitude at position, Indicates in The baseline estimate of the location, Indicates in Positional residuals Indicates the first In terms of scale, Gradient magnitude at position, express The baseline estimate of the location, express Positional residuals; S5.2: To reflect scale differences, a scale-related unit step size decay factor is introduced. Defined as: ; in, Indicates scale. This represents the baseline attenuation coefficient, i.e., the value chosen in the experiment; scale. The larger, The closer to 1, the slower the residual decays within the segment, resulting in a smoother transition with a wider effective range; scale The smaller, Smaller size results in faster decay and a more localized transition; selection can be achieved through experiments. The value of achieves a balance between smoothness and locality.

[0041] S5.3: Utilize the residuals described at both ends With residual and attenuation factor Constructing a two-ended exponential interpolation, we obtain the intra-segment correction term. Its analytical form can be written as a weighted combination of the exponential terms at both ends, as shown in the following expression: ; in, Indicates the current position. This represents the scale-dependent unit step size decay factor. S5.4: Update the current position The gradient strength is calculated using the following formula: ; in, This represents the magnitude of the row gradient after interpolation update. Indicates the current benchmark estimate; Column interpolation: S5.5: For fixed columns Adjacent anchor point pairs An interpolation interval is defined. Based on the current benchmark estimate The residuals at both ends are calculated using the following formula: ; ; in, Indicates the first In terms of scale, Gradient magnitude at position, Indicates in The baseline estimate of the location, Indicates in Positional residuals Indicates the first In terms of scale, Gradient magnitude at position, express The baseline estimate of the location, express Positional residuals; S5.6: Introduce a scale-dependent unit step size decay factor Defined as: ; in, Indicates scale. This represents the baseline attenuation coefficient, the value selected in the experiment; S5.7: Utilize the residuals described at both ends With residual and unit step size decay factor Constructing a two-ended exponential interpolation, we obtain the intra-segment correction term. Its analytical form can be written as a weighted combination of the exponential terms at both ends, as shown in the following expression: ; S5.8: Update the current position The gradient strength is calculated using the following formula: ; in, This represents the magnitude of the column gradient after interpolation update; For boundary intervals with anchor points on only one side, a single-ended exponential extension is used: the residual on the side closest to the anchor point is attenuated along the outer edge of the segment by a unit step size decay factor. The propagation is attenuated, and the residual contribution is kept at a distance that is far from the anchor point and can be ignored. When the number of anchor points in a row or column is less than 2, the corresponding row or column is not updated by interpolation.

[0042] Row and column clipping correction: After the interpolation process is completed, the row and column gradient fields at each scale have been filled with dense distribution instead of sparse anchor points. However, interpolation inevitably introduces secondary peaks, steps, and local oscillations between anchor points. Therefore, this invention maintains the orientation angle... Assuming the gradient magnitude remains unchanged, only the gradient magnitude is considered. Perform uniform monotonicity and de-oscillation correction along the main direction to make the profile between anchor point pairs show a shape that "monotonically converges from both ends to the intra-segment minimum value", thereby suppressing the pseudo texture introduced by interpolation. S5.9 To avoid over-correction in directions inconsistent with the current edge orientation, this invention designs a direction selection criterion based on the orientation angle and selects two constant thresholds. , It calculates the degree of deviation between the direction angle and the vertical reference direction, and determines whether the corresponding direction is horizontal or vertical, thereby achieving precise orientation correction; The orientation angle design orientation selection criterion is: when When the direction of the corresponding pixel is considered to be far from the vertical and close to the horizontal, row clipping is triggered. when When the direction of the corresponding pixel is considered to be close to vertical, column clipping is triggered; where Indicates the first Each decomposition scale at position The direction angle, , This represents the two selected constant thresholds; S5.10: Row-wise clipping: For fixed rows A segment, i.e., the interpolation interval The gradient magnitude sequence within a segment is denoted as: in, This represents the sequence of row gradient magnitudes after interpolation update; The minimum gradient magnitudes and positions of the gradient magnitude sequence within the segment are as follows: ; in, Indicates the first Line number The minimum gradient magnitude within the column interpolation interval, This indicates all pixels within the interval. gradient magnitude sequence Take the minimum value. This represents the minimum value. Indicates the phase at the minimum value; When the "horizontal" direction threshold is met, perform two monotonic scans along the row direction: Left sweep: from the left end position Towards the minimum position Scan, location ,…, ,like ; in, Let represent the gradient magnitude at the left end of the row direction, then let ; in, This represents the gradient magnitude in the direction following the left sweep. Right sweep: from the right end Towards Scan: Location ,like ; in, Let represent the gradient magnitude at the right end of the row direction, then let ; in, This represents the gradient magnitude in the direction following the right sweep. Through two monotonic scans in the row direction, the gradient intensity curve within the segment monotonically moves from both ends towards the segment's minimum along the row direction. Convergence effectively suppresses secondary peaks and jagged edges caused by interpolation, while maintaining the direction angle unchanged to ensure a smooth transition consistent with the edge direction.

[0043] S5.11: Column-wise clipping: For fixed rows A segment, i.e., the interpolation interval The gradient magnitude sequence within a segment is denoted as: ; in, This represents the sequence of column gradient magnitudes after interpolation update; The minimum gradient magnitudes and positions of the gradient magnitude sequence within the segment are as follows: ; in, Indicates the first Line number The minimum gradient magnitude within the column interpolation interval, This indicates all pixels within the interval. gradient magnitude Take the minimum value. This represents the minimum value. Indicates the phase at the minimum value; When the "vertical" direction threshold is met, perform two monotonic scans along the column direction: Left sweep: from the left end position Towards the minimum position Scan, location ,…, ,like ; in, Let represent the gradient magnitude at the left end of the column direction, then let ; in, This represents the gradient magnitude in the column direction after the left sweep. Right sweep: from the right end Towards Scan: Location ,like ; in, This indicates the gradient magnitude at the right end of the column direction; Then let ; in, This indicates the gradient magnitude in the column direction after the right sweep; By performing two monotonic scans along the column direction, the gradient magnitudes that do not satisfy monotonicity are projected onto nearby smaller values, causing the column profile to also monotonically converge from both ends to [the desired value]. Remove secondary peaks and jagged edges introduced by interpolation.

[0044] After row and column cropping corrections, the uneven transitions in the image are effectively suppressed, and the pseudo-textures and noise introduced by interpolation are significantly reduced; the corrected anchor points are smoother and more consistent in the main direction; finally, while maintaining the orientation angle... Without changing the gradient magnitude, restore it to its horizontal and vertical components: ; ; in, The cosine of the direction angle is used to measure the gradient magnitude. Project back to the horizontal component. The sine of the direction angle is used to measure the gradient magnitude. Project back to the vertical component. Represents the restored horizontal component This represents the restored vertical component; This results in a multi-scale gradient field with consistent direction and continuous details, providing a stable coefficient basis for the final watermark image reconstruction.

[0045] S6, Wavelet Reconstruction Correction of Multi-Scale Gradients: Wavelet reconstruction correction is performed on the gradient component field, i.e., the subband coefficients, to obtain a new set of subband coefficients, i.e., multi-scale wavelet coefficients. These multi-scale wavelet coefficients are then combined to obtain the image. .

[0046] Specifically, it is based on the gradient component field, i.e., the sub-band coefficient, in S5. Perform inverse discrete wavelet transform to reconstruct the intermediate image. ;in, This represents the high-frequency horizontal components that maintain a consistent direction and continuous detail after cropping and trimming. This represents the high-frequency components in the vertical direction that maintain a consistent direction and continuous details after cropping and trimming. This represents low-frequency component information that maintains consistent direction and continuous detail after cropping and correction; the intermediate image. The grayscale range is linearly normalized, and a new set of multi-scale wavelet coefficients is obtained through discrete wavelet transform. ,in, This represents the updated high-frequency components in the horizontal direction. This represents the high-frequency components in the vertical direction after the update. This indicates the updated low-frequency component; Updated low-frequency components Replace with the original low-frequency component obtained in S2 This allows for updating only the high-frequency gradient subbands at each scale while maintaining consistency between the low-frequency structure and the initial image. and Finally, the original low-frequency components... With the updated high-frequency components , By combining the results, the final reconstructed image is obtained. The updated multi-scale gradient should satisfy the integral consistency that can be generated from a real image, while preserving as much of the overall contour and brightness distribution of the original image as possible.

[0047] Among them, the gradient component field based on S5 is the sub-band coefficient. Perform inverse discrete wavelet transform to reconstruct the intermediate image. The reconstruction process is as follows: ; in, The first Low-pass and high-pass reconstruction filters in the row and column directions at each decomposition scale; This indicates the reconstruction of the low-frequency components in the row and column directions; This indicates the reconstruction of the high-frequency components in the row direction; This indicates the reconstruction of the high-frequency components in the column direction; A new set of multiscale wavelet coefficients is obtained through discrete wavelet transform. The decomposition process is as follows: ; Among them, the filter The first Low-pass decomposition filter and high-pass decomposition filter at each decomposition scale; Indicates the image at the 1st Low-frequency subbands in the horizontal direction of each decomposition scale. Indicates the first The updated high-frequency components in the horizontal direction at each decomposition scale Indicates the first The high-frequency components in the vertical direction after updating at each decomposition scale Indicates the first The updated low-frequency components at each decomposition scale.

[0048] Experimental evaluation: This invention focuses on the experimental evaluation of a watermarking method based on Fresnel domain multi-scale sparse interest points. First, the effectiveness of the multi-scale interest point detection and reconstruction process is analyzed from the perspective of sparse representation and reconstruction. Then, under the complete watermark embedding and extraction process, the comprehensive performance of the invention is evaluated from the aspects of invisibility, robustness, and method comparison. All experiments were completed in the MATLAB R2016b environment and were uniformly conducted in grayscale image scenarios.

[0049] Multi-scale interest point extraction analysis: To verify the effectiveness of multi-scale interest point extraction in the feature extraction stage, this invention experimentally evaluated the multi-scale interest point extraction process; such as... Figure 6 As shown, this is a multi-scale interest point distribution map of the present invention, where (a1) represents the original image; (a2) represents the low-frequency approximation component; (b1)–(c1) represent horizontal high-frequency details at different scales; and (b2)–(c2) represent vertical high-frequency details at different scales. The present invention demonstrates the extraction results at the first and second scales. Stable interest points can be extracted from the image through spatial orientation encoding and directional gradient maximum detection. Interest points are mainly distributed in structurally significant regions such as edges and textures. At the same time, multi-scale decomposition makes high-frequency details mainly concentrated in contours and local change locations, while low-frequency components retain a smoother overall contour. Through this multi-scale representation, the key structural information of the image is preserved, and the representation of watermark features is sparser, thereby significantly reducing the computation and storage overhead of subsequent embedding and restoration stages.

[0050] Multi-scale interest point reconstruction quality analysis of images of different sizes: To verify the reconstruction capability on images of different sizes, this invention conducts multi-scale interest point extraction and reconstruction experiments on images of various sizes, and analyzes the impact of interest point sparsity on reconstruction quality. It should be noted that all portraits in the accompanying drawings are drawn images, not real people. In images of different sizes, although interest points only cover partial areas, through interpolation and cropping correction, the details, textures, and edge structures of the reconstructed images can be recovered relatively fully, with minimal subjective visual differences. Figure 7 and Figure 8 The reconstruction results of images at different sizes are shown; Figure 7 This is a multi-scale interest point extraction and reconstruction result image of a 256×256 image of the present invention, where (a1)-(e1) represent the original image, (a2)-(e2) represent the interest point extraction result, the green area represents the location of the extracted interest points, and (a3)-(e3) represent the reconstruction result using the extracted interest points; Figure 8This is a multi-scale interest point extraction and reconstruction result diagram of a 512×512 image of the present invention, where (a1)-(e1) represent the original image, (a2)-(e2) represent the interest point extraction result, the green area represents the location of the extracted interest points, and (a3)-(e3) represent the reconstruction result using the extracted interest points.

[0051] Table 1 shows the number of interest points (OPPs), the proportion of OOPs, and the corresponding peak signal-to-noise ratio (PSNR) and structural similarity (SSIM) for each image. The results show that when the number of OOPs accounts for about 20% of the total number of pixels, the present invention can still achieve high reconstruction quality, with the PSNR typically exceeding 26 dB and reaching a maximum of 37.23 dB. This result indicates that even under conditions of sparse OOPs, the present invention can still effectively recover the main structure and details of the image, thus verifying the effectiveness of the multi-scale OOP extraction and reconstruction strategy in sparse reconstruction tasks. Furthermore, under the experimental settings, as the image size increases, the OOP coverage ratio and reconstruction index generally show an upward trend, indicating that the present invention has good adaptability and stability to images of different sizes.

[0052] Table 1. Reconstruction quality of images of different sizes under multi-scale interest point reconstruction.

[0053] Analysis of the impact of multi-scale interest point sparsity on reconstruction quality: In the above experiments, for larger images, there are more interest points; for example, in a 256×256 image, more than 10,000 interest points were detected. Although they only account for about 20% of the total pixels, they have covered the key information of the image. Therefore, the reconstruction quality based on these interest points is still relatively high. In order to verify the impact of a smaller number of interest points on the image reconstruction quality, a 64×64 image was selected, and 5%–100% of the detected interest points were randomly retained in proportion to analyze the impact of interest point density on reconstruction quality and its stability.

[0054] Figure 9 This is an image reconstruction result diagram of the image under different interest point retention ratios according to the present invention, where (a1)-(f1) represent the interest point distribution when 10%, 30%, 50%, 70%, 90%, and 100% of interest points are retained, respectively; (a2)-(f2) represent the corresponding reconstruction watermark results, where green represents all detected interest points, and red represents randomly selected interest points under the retention ratio; from Figure 9It can be observed that when the interest point retention ratio is 10%, 30%, 50%, 70%, 90%, and 100%, the spatial distribution of interest points gradually changes from sparse to dense, and the corresponding reconstructed image gradually transitions from only being able to outline the general contour to having basically complete details. Table 2 shows the objective evaluation results under different retention ratios. When the retention ratio is 5%, the number of interest points is 38, accounting for 0.93% of the total pixels, the peak signal-to-noise ratio (PSNR) is 19.056 dB, and the structural similarity (SSIM) is 0.9429. Less than 1% of the pixels can already recover a relatively complete contour structure. Overall, a small number of interest points show good stability in terms of image quality restoration, and as the retention ratio increases, the reconstruction quality is continuously and significantly improved.

[0055] Table 2 Peak signal-to-noise ratio and structural similarity of reconstructed images under different interest point retention ratios.

[0056] Imperceptibility: Imperceptibility is one of the core metrics for evaluating the performance of digital watermarking algorithms. It is used to assess the impact of watermark embedding on the visual quality of the host image. An ideal watermarking scheme should maintain perceptual consistency with the original image as much as possible while ensuring robustness and a certain embedding capacity, making it difficult for observers to detect the presence of the watermark.

[0057] Figure 10 The host images used in this invention are as follows: (a) represents a woman, (b) represents a baboon, (c) represents a house, (d) represents a sailboat, (e) represents a man, and (f) represents an airplane. In this invention, the host image size is 256×256, representing a woman, a baboon, a house, a sailboat, a man, and an airplane, respectively. Four 64×64 images were selected for the watermark: a rose, text, a teacher, and a cat, as shown below. Figure 11 The watermark images used in this invention are shown in the following diagrams: (a) represents a rose, (b) represents text, (c) represents a teacher, and (d) represents a cat.

[0058] While balancing robustness and imperceptibility, this application selects an embedding strength with superior overall performance. This method is compared to the robust image watermarking method based on chaotic sequences disclosed in JAP Artiles, DPB Chaves, C. Pimentel, Robust imagewatermarking algorithm using chaotic sequences, Journal of Information Security and Applications 68 (2022) (hereinafter referred to as Method 1), and NA Azam, T. Haider, U. Hayat, An optimized watermarking scheme based on genetic algorithm and elliptic curve, Swarm and Evolutionary Computation 91. The method for optimizing watermarking based on genetic algorithm and elliptic curve (hereinafter referred to as Method 2) disclosed in (2024) was compared; Table 3 summarizes the peak signal-to-noise ratio (PSNR) and structural similarity (SSIM) results under different host image and watermark combinations; it can be seen that the method of the present invention has an average peak signal-to-noise ratio (PSNR) of about 53-61dB on each test image, and the structural similarity (SSIM) is close to 1.000. Compared with Method 1 and Method 2, it has a higher peak signal-to-noise ratio (PSNR) and better structural similarity (SSIM), indicating that it can maintain high visual quality in both complex and smooth scenes.

[0059] Table 3 Comparison of Peak Signal-to-Noise Ratio and Structural Similarity Results of Watermarked Host Images under the Invention and Related Algorithms

[0060] To further examine the impact of watermark embedding on pixel statistical characteristics, this application uses a teacher's watermark as an example and compares the grayscale histograms of the host image before and after embedding. Figure 12 As shown in the watermark embedding of this invention, the grayscale distribution curves of the two methods basically overlap, with only slight shifts in local grayscale ranges, indicating that the watermark embedding has minimal disturbance to the overall grayscale structure. Table 3 compares the peak signal-to-noise ratio and structural similarity results of the watermarked host image under the proposed method and related algorithms. Figure 13 The results of the average normalized correlation coefficient under signal processing attacks of the present invention show that the present invention can maintain a high peak signal-to-noise ratio (PSNR) and structural similarity (SSIM) under various host images and different types of watermarks, while taking into account the fidelity of visual and statistical characteristics, and exhibiting superior performance in terms of imperceptibility.

[0061] Robustness analysis: To verify the stability of this invention under various distortion conditions, system experiments were conducted from three aspects: signal processing attack, geometric attack, and combined attack. Finally, the results were compared with several representative watermarking algorithms. The normalized correlation coefficient was used as the main evaluation index for all experiments. The closer the value is to 1, the more reliable the watermark recovery is under the corresponding attack.

[0062] Signal processing attacks: First, we examine the impact of common signal processing operations on the proposed algorithm, including scenarios such as noise perturbation, smoothing filtering, and lossy compression. The normalized correlation coefficient variation curves are shown below. Figure 13 The average normalized correlation coefficient value of the present invention under signal processing attacks is shown.

[0063] like Figure 13 The average normalized correlation coefficient values ​​of this invention under signal processing attacks are shown in (a) salt-and-pepper noise, (b) speckle noise, and (c) Gaussian noise. Under salt-and-pepper noise, speckle noise, and Gaussian noise interference, the average normalized correlation coefficient shows a monotonically decreasing trend with increasing noise intensity, but it still remains at a relatively high level overall. When the noise density or variance does not exceed 0.01, the average normalized correlation coefficient is higher than 0.99, indicating that the watermark can be recovered almost without distortion in a low-noise environment. Even when the noise intensity increases to 0.05, the average normalized correlation coefficient is still around 0.97, and the difference between different host images is very small, indicating that the algorithm has relatively consistent noise resistance performance under different texture conditions.

[0064] like Figure 13 The average normalized correlation coefficient values ​​of this invention under signal processing attacks are shown in (d) median filtering; (e) Gaussian low-pass filtering; and (f) mean filtering. Under the effects of median filtering, low-pass Gaussian filtering, and mean filtering, the average normalized correlation coefficient value gradually decreases as the filter kernel size or filtering intensity increases. For medium-intensity filters such as 3×3 or 5×5, the average normalized correlation coefficient mostly remains in the range of 0.97-0.98, still effectively recovering the watermark; when the filter kernel is expanded to 7×7 or 9×9, the average normalized correlation coefficient drops to approximately 0.88-0.93. Overall, the algorithm can still maintain a high watermark similarity under light to moderate filtering, but its performance slightly decreases under strong smoothing conditions.

[0065] like Figure 13The average normalized correlation coefficient values ​​of this invention under signal processing attacks are shown in (g) JPEG compression and (h) JPEG2000 compression. The proposed algorithm also exhibits good stability in both JPEG and JPEG2000 compression scenarios. In JPEG compression, when the quality factor increases from 10 to 90, the average normalized correlation coefficient increases from approximately 0.90 to 0.99. When the quality factor is greater than 70, the average normalized correlation coefficient remains relatively stable above 0.98. In JPEG2000, when the compression ratio increases from 2:1 to 15:1, the average normalized correlation coefficient remains consistently above 0.985. These results demonstrate that the proposed method has strong robustness and stability against mainstream lossy compression.

[0066] like Figure 15 The watermark image extracted under the combined attack of this invention shows that the proposed algorithm also exhibits good stability in JPEG and JPEG2000 compression scenarios. In JPEG compression, when the quality factor increases from 10 to 90, the average normalized correlation coefficient increases from approximately 0.90 to 0.99, and when the quality factor is greater than 70, the average normalized correlation coefficient remains relatively stable above 0.98. In JPEG2000, when the compression ratio increases from 2:1 to 15:1, the average normalized correlation coefficient remains consistently above 0.985. These results demonstrate that the proposed method has strong robustness and stability against mainstream lossy compression.

[0067] Geometric attack analysis: Geometric distortion mainly includes two typical operations: global rotation and local cropping. Their impact on the algorithm is as follows: Figure 14 The average normalized correlation coefficient value of this invention under geometric attacks is shown. As... Figure 14 The average normalized correlation coefficient values ​​of this invention under geometric attacks are shown in (a) by rotation angle. Under different rotation angles, this invention maintains a consistently high correlation. At a rotation of 30°, the average normalized correlation coefficient exceeds 0.945, with the average normalized correlation coefficients for the Barbara and man images reaching 0.9773 and 0.9723 respectively, indicating almost no significant distortion of the watermark. When the rotation angle increases to 70° and 90°, the average normalized correlation coefficient decreases slightly but remains above 0.91. Even when further increased to 110° and 150°, the average normalized correlation coefficients for some images still reach 0.8980 and 0.9625, with an overall decrease of no more than 0.07, demonstrating that the algorithm maintains relatively stable extraction performance under large-angle rotations.

[0068] like Figure 14The average normalized correlation coefficient (ARC) under geometric attacks of this invention is shown in (b), which represents the cropping ratio. Under different cropping ratios, the ARC gradually decreases with increasing cropping area, but the overall level remains acceptable. With a slight cropping of 1 / 64, the ARC is approximately 0.98, hardly affecting watermark recovery. When the cropping ratio increases to 1 / 16 and 1 / 8, the ARC is approximately 0.93 and 0.91, respectively. Under a more severe 1 / 2 cropping condition, the ARC remains around 0.86, indicating that the watermark information is not completely destroyed. The differences between different host images are small. House and Barbara perform relatively best, with an ARC exceeding 0.89 at a cropping ratio of 1 / 4. Sailboat, with its more complex texture, is slightly more affected by cropping, but the overall difference does not exceed 0.04. These results demonstrate that even with a large area of ​​missing local information, the proposed scheme can still achieve high watermark recoverability by relying on embedding redundancy.

[0069] Combination attacks: To further verify the robustness of the proposed algorithm under multiple distortion conditions, this invention conducts experimental evaluations on several typical combined attacks. The combined methods mainly cover three scenarios: (1) superposition of noise and filtering, such as joint median or mean filtering of salt-pepper noise / speckle noise; (2) superposition of compression and spatial operations, such as JPEG / JPEG2000 compression combined with cropping; (3) joint of geometric transformation and smoothing filtering, such as the combination of rotation, cropping and low-pass Gaussian filtering. The normalized correlation coefficient results under different attack configurations are shown in Table 4. The normalized correlation coefficient values ​​of this invention under combined attacks are shown in Table 4. The extracted watermark images are as follows. Figure 15 The watermark image extracted under the combined attack of this invention is shown in the figure.

[0070] As can be seen from the table, the algorithm performs relatively stably under combined noise and filtering attacks. Taking "Salt&pepper(0.01)+Median(5×5)" as an example, the normalized correlation coefficients of the four watermarks are all above 0.92, with an average normalized correlation coefficient of 0.9414. Under the combined attacks of "Speckle(0.01)+Average(5×5)" and "JPEG(50)+Average(5×5)", the average normalized correlation coefficient remains around 0.93, indicating that the watermark can still be reliably recovered even when point noise, multiplicative noise, and smoothing filtering are superimposed. In the combined attacks of compression and geometry, the performance decreases to some extent but remains within an acceptable range. For example, when the JPEG2000 compression ratio is 4:1 and the cropping ratio is 6.25%, the average normalized correlation coefficient is 0.9067. When the image undergoes a 45° rotation accompanied by a 25% cropping, the average normalized correlation coefficient drops to 0.8769. If a 25% cropping is performed first, followed by a 5° rotation, the average normalized correlation coefficient further decreases to 0.8569, indicating that the superposition of geometric operations has a more significant impact on the accuracy of watermark restoration. Under the triple attack of cropping, rotation, and Gaussian low-pass filtering, "Cropping(25%)+Rotation(15°)+Low-pass Gaussian(5×5)", the average normalized correlation coefficient is 0.8641, which is lower than that of a single attack, but still significantly higher than the level of complete failure.

[0071] Overall, the proposed algorithm has an average normalized correlation coefficient higher than 0.90 under most combined attack conditions. Its performance only shows a significant decline in strong geometric composite attacks involving large-scale pruning or large-angle rotation, demonstrating good robustness in complex distortion environments.

[0072] Table 4 Normalized correlation coefficient values ​​of the present invention under combined attacks

[0073] Comparative analysis: To verify the comprehensive robustness performance of this invention, it is compared with three representative watermarking algorithms: the optimized watermarking method based on genetic algorithm and elliptic curve disclosed in JAP Artiles, DPB Chaves, C. Pimentel, Robust image watermarking algorithm using chaotic sequences, Journal of Information Security and Applications 68 (2022) (hereinafter referred to as Method 1); the optimized watermarking method based on genetic algorithm and elliptic curve disclosed in NA Azam, T. Haider, U. Hayat, An optimized watermarking scheme based on genetic algorithm and elliptic curve, Swarm and Evolutionary Computation 91 (2024) (hereinafter referred to as Method 2); and R. Thanki, A. Kothari, D. Trivedi, Hybrid and blind watermarking scheme in DCuT-RDWT domain, Journal of Information Security and Applications 46 (2019). The hybrid blind watermarking method based on the DCuT-RDWT domain disclosed in 231-249 (hereinafter referred to as Method 3) was compared under various typical attacks; Table 5 shows the normalized correlation coefficients of each algorithm under different attack scenarios.

[0074] Under attacks from salt-pepper noise and Gaussian noise, the average normalized correlation coefficient of the four watermarks in this invention is around 0.985. Among them, the average normalized correlation coefficient is close to 0.997 when the noise intensity is 0.001, which is higher than or close to the existing methods, indicating that the watermark can be stably recovered under mild to moderate random noise interference. In JPEG compression attacks, the present invention maintains an average normalized correlation coefficient of approximately 0.94–0.97 under various quality factor settings of 70–40, which is roughly equivalent to the performance of the comparison algorithm, and has a slight advantage under some strong compression conditions, indicating that the present invention has good adaptability to lossy compression. In terms of filtering attacks, the average normalized correlation coefficient of the present invention under 3×3 and 5×5 median filtering is approximately 0.98 and 0.96, respectively, which is significantly better than the results of method 3, indicating that the embedding structure obtained by using multi-scale interest points and direction interpolation has stronger resistance to smoothing filtering. In geometric attacks, the average normalized correlation coefficient of the present invention under 20° and 45° rotation and 0.25 scaling conditions is close to or higher than 0.98, which is significantly better than method 3 and method 2, and has a significant improvement over some algorithms under scaling attacks, demonstrating robustness advantages in spatial transformation scenarios.

[0075] In summary, the present invention maintains a high average normalized correlation coefficient under various typical attacks such as noise, compression, filtering, and geometric transformation. In most cases, it is superior to or no worse than the comparison method, indicating that the present invention has good comprehensive performance in terms of robustness and stability.

[0076] Table 5 compares the normalized correlation coefficients of the present invention and existing methods under typical attacks.

[0077] Of course, the above description is not intended to limit the present invention, and the present invention is not limited to the examples given above. Any changes, modifications, additions or substitutions made by those skilled in the art within the scope of the present invention should be protected by the present invention.

Claims

1. A watermarking method based on Fresnel domain multi-scale sparse interest points, characterized in that, Includes the following steps: S1, Fresnel domain preprocessing: Host image A is preprocessed to obtain the host subband singular value matrix. and ; S2. Multi-scale interest point extraction: The watermark image is decomposed by wavelet and the gradient magnitude and orientation angle of the image are calculated. After obtaining the spatial orientation code of each pixel through spatial orientation coding, the gradient magnitude is subjected to multi-scale local maximum detection with consistent direction and effective features are extracted. S3, Watermark Embedding and Reconstruction: Select the watermark features to be embedded at the target scale, perform Arnold scrambling and singular value decomposition to obtain the singular value diagonal matrix of the watermark features, and then sort them according to the watermark strength parameters. Embedded in the singular value matrix and In the process, the final watermarked host image is obtained through inverse transformation reconstruction. ; S4. Watermark Extraction and Decryption: Watermarked host image The singular value matrix of watermark recovery is obtained through preprocessing. The singular value matrix is ​​used to obtain the singular value estimation matrix corresponding to the watermark through difference operation. Inverse singular value decomposition is performed to reconstruct the watermark gradient magnitude at the target scale. The watermark gradient magnitude at the target scale is then subjected to Arnold inverse scrambling to obtain the decrypted gradient magnitude. S5. Multi-scale interest point reconstruction: For the The results of the scale maxima detection are used to construct an anchor point set. The anchor points are then segmented in both the row and column directions to form an interpolation interval. One-dimensional interpolation is performed on each pair of adjacent anchor points in the interpolation interval. Based on the orientation angle selection criterion, the gradient magnitude is pruned and corrected in both the row and column directions. While keeping the orientation angle unchanged, the gradient magnitude is restored to the horizontal and vertical components, thus obtaining a gradient component field with consistent direction and continuous details, i.e., subband coefficients. S6, Wavelet Reconstruction Correction of Multi-Scale Gradients: Wavelet reconstruction correction is performed on the gradient component field, i.e., the subband coefficients, to obtain a new set of subband coefficients, i.e., multi-scale wavelet coefficients. These multi-scale wavelet coefficients are then combined to obtain the image. .

2. The method according to claim 1, characterized in that, The preprocessing operations S1 and S4 include: at a given wavelength and transmission distance Under the given conditions, the complex-valued coefficient matrix is ​​obtained by mapping to the Fresnel domain through Fresnel transform, and then converted to real-valued form by real-preservation transform to obtain a real-valued coefficient matrix. The real-valued coefficient matrix is ​​then normalized and quantized to obtain a grayscale image, i.e., the quantized Fresnel domain image. The Fresnel domain image is then subjected to discrete wavelet transform to obtain four sub-bands: the low-frequency sub-band of the Fresnel domain image, the high-frequency sub-band in the horizontal direction of the Fresnel domain image, the high-frequency sub-band in the vertical direction of the Fresnel domain image, and the high-frequency sub-band in the diagonal direction of the Fresnel domain image. Singular value decomposition is then performed on the low-frequency sub-band and the high-frequency sub-band in the diagonal direction of the Fresnel domain image to obtain a singular value matrix.

3. The method according to claim 1, characterized in that, The process of wavelet decomposing the watermarked image and calculating its gradient magnitude and orientation angle includes: performing discrete wavelet multi-scale decomposition on the watermarked image to obtain the sub-band coefficients at each scale. ,in, This represents the horizontal information of high-frequency components. Represents the vertical information of high-frequency components. Indicates low-frequency components; in the first At each decomposition scale, the corresponding gradient magnitude is calculated using the horizontal and vertical gradients of the image. With direction angle The calculation formula is as follows: ； ； in, Represents the gradient in the horizontal direction. Represents the gradient in the vertical direction. Indicates the first Location at each decomposition scale The gradient magnitude of a pixel. Indicates the first Location at each decomposition scale The orientation angle of a pixel is the first Location at each decomposition scale The phase of a pixel; The spatial orientation encoding operation obtains the spatial orientation code for each pixel, including: the orientation angle. The image is mapped to a preset angle range, and the preset angle range is divided into multiple sectors. Based on the multiple sectors, a spatial orientation code is assigned to each pixel. For gradient magnitude Execution of multi-scale local maxima detection and extraction of effective features with consistent direction includes: for the first Location at each decomposition scale The spatial encoding of pixels involves selecting neighboring pixels along the corresponding direction and comparing them with the first... Location at each decomposition scale The gradient magnitude of a pixel and its neighboring pixels are used to determine whether they are local maxima. Pixels that are greater than their neighbors on both sides are marked as local maxima with the same direction, while the rest are suppressed. After suppressing non-maxima at all scales, the points that remain in each decomposition layer that are local maxima at multiple scales are the effective features.

4. The method according to claim 1, characterized in that, The inverse transform reconstruction process yields the final watermarked host image. Including: by watermark strength parameters Embedded in the singular value matrix and Afterwards, the low-frequency subband with embedded watermark is obtained through the Inverse Singular Value Decomposition (ISVD) operation. and the high-frequency subband in the horizontal direction after embedding watermark The low-frequency subband and high-frequency subband Together with the remaining subbands not involved in embedding, an inverse wavelet transform is performed to obtain the Fresnel domain real-valued coefficient matrix containing the watermark. This Fresnel domain real-valued coefficient matrix undergoes inverse quantization and is then remapped back to the spatial domain via an inverse Fresnel transform to obtain the final watermarked host image. .

5. The method according to claim 1, characterized in that, The singular value matrix is ​​used to obtain the singular value estimation matrix corresponding to the watermark through difference operations, including: using the singular value matrix... and and watermark strength parameters The estimated value of the watermark component is recovered in the singular value domain; the calculation formula is as follows: ； ； in, This indicates that after preprocessing and Discrete Wavelet Transform (DWT) decomposition, the low-frequency subband of the watermarked image is... The singular value matrix obtained by performing singular value decomposition (SVD); This indicates that after preprocessing and Discrete Wavelet Transform (DWT) decomposition, the high-frequency subband of the watermarked image is... The singular value matrix obtained by performing singular value decomposition (SVD); This represents the singular value estimation matrix corresponding to the horizontal features at the second scale of the watermarked image; This represents the singular value estimation matrix corresponding to the vertical features at the second scale of the watermarked image.

6. The method according to claim 1, characterized in that, For the first The results of the scale maxima detection are used to construct an anchor point set, and the anchor points are segmented in both the row and column directions, including: Row direction: For fixed rows Collect the fixed rows The column indexes of all anchor points are used, and sorted from left to right to obtain the anchor point sequence; ； in, For the first Line up Anchor column index, This represents the total number of anchor points in that row. Indicates the first The set of anchor indexes for rows. , , ..., These represent the column indexes of the first anchor point, the second anchor point, and the third anchor point, respectively. Column indexes for anchor points; Column direction: For fixed columns Collect the fixed column The row indices of all anchor points are sorted from top to bottom to obtain: ； in, Indicates the first List of the first Row index of anchor points This represents the total number of anchor points in this column. Indicates the first Anchor index set of columns, , , ..., These represent the row indexes of the first anchor point, the second anchor point, and the third anchor point, respectively. Row index of anchor points.

7. The method according to claim 3, characterized in that, The interpolation interval performs one-dimensional interpolation for each pair of adjacent anchor points, including: Row-wise interpolation: S5.1: Each pair of adjacent anchor points naturally forms an interpolation interval, for the first... Each scale, based on anchor point segmentation, is used for fixed rows. Adjacent anchor points An interpolation interval is defined. Based on the current benchmark estimate The residuals at both ends are calculated using the following formula: ； ； in, Indicates the first In terms of scale, Gradient magnitude at position, Indicates in The baseline estimate of the location, Indicates in Positional residuals Indicates the first In terms of scale, Gradient magnitude at position, express The baseline estimate of the location, express Positional residuals; S5.2: Introduce a scale-dependent unit step size decay factor Defined as: ； in, Indicates scale. This represents the baseline attenuation coefficient, i.e., the value selected in the experiment; S5.3: Utilize the residuals described at both ends With residual and attenuation factor Constructing a two-ended exponential interpolation, we obtain the intra-segment correction term. Its analytical form can be written as a weighted combination of the exponential terms at both ends, as shown in the following expression: ； in, Indicates the current position. This represents the scale-dependent unit step size decay factor. S5.4: Update the current position The gradient strength is calculated using the following formula: ； in, This represents the magnitude of the row gradient after interpolation update. Indicates the current benchmark estimate; Column interpolation: S5.5: For fixed columns Adjacent anchor point pairs An interpolation interval is defined. Based on the current benchmark estimate The residuals at both ends are calculated using the following formula: ； ； in, Indicates the first In terms of scale, Gradient magnitude at position, Indicates in The baseline estimate of the location, Indicates in Positional residuals Indicates the first In terms of scale, Gradient magnitude at position, express The baseline estimate of the location, express Positional residuals; S5.6: Introduce a scale-dependent unit step size decay factor Defined as: ； in, Indicates scale. This represents the baseline attenuation coefficient, the value selected in the experiment; S5.7: Utilize the residuals described at both ends With residual and unit step size decay factor Constructing a two-ended exponential interpolation, we obtain the intra-segment correction term. Its analytical form can be written as a weighted combination of the exponential terms at both ends, as shown in the following expression: ； S5.8: Update the current position The gradient strength is calculated using the following formula: ； in, This represents the magnitude of the column gradient after interpolation.

8. The method according to claim 7, characterized in that, For boundary intervals with anchor points on only one side, a single-ended exponential extension is used: the residual on the side closest to the anchor point is attenuated along the outer edge of the segment by a unit step size decay factor. Attenuation propagation: When the number of anchor points in a row or column is less than 2, no interpolation update is performed on the corresponding row or column.

9. The method according to claim 7, characterized in that, The gradient magnitude is pruned and corrected in both row and column directions, including: S5.9: Select two constant thresholds, calculate the degree of deviation between the direction angle and the vertical reference direction, determine whether the corresponding direction is horizontal or vertical, and achieve precise orientation correction; The criterion for determining whether the corresponding direction is horizontal or vertical is: when When the corresponding pixel is considered to be moving away from the vertical direction and closer to the horizontal direction, line clipping is triggered; when When the direction of the corresponding pixel is considered to be close to vertical, column clipping is triggered; where Indicates the first Location at each decomposition scale The direction angle, , This represents the two selected constant thresholds; S5.10: Row-wise clipping: For fixed rows A segment, i.e., the interpolation interval The gradient magnitude sequence within a segment is denoted as: in, This represents the sequence of row gradient magnitudes after interpolation update; The minimum gradient magnitudes and positions of the gradient magnitude sequence within the segment are as follows: ； in, Indicates the first Line number The minimum gradient magnitude within the column interpolation interval, This indicates all pixels within the interval. gradient magnitude sequence Take the minimum value. This represents the minimum value. Indicates the phase at the minimum value; When the "horizontal" direction threshold is met, perform two monotonic scans along the row direction: Left sweep: from the left end position Towards the minimum position Scan, location ,…, ,like ； in, Let represent the gradient magnitude at the left end of the row direction, then let ； This represents the gradient magnitude in the direction following the left sweep. Right sweep: from the right end Towards Scan: Location ,like ； in, Let represent the gradient magnitude at the right end of the row direction, then let ； in, This represents the gradient magnitude in the direction following the right sweep. S5.11: Column-wise clipping: For fixed rows A segment, i.e., the interpolation interval The gradient magnitude sequence within a segment is denoted as: ； in, This represents the sequence of column gradient magnitudes after interpolation update; The minimum gradient magnitudes and positions of the gradient magnitude sequence within the segment are as follows: ； in, Indicates the first Line number The minimum gradient magnitude within the column interpolation interval, This indicates all pixels within the interval. gradient magnitude Take the minimum value. This represents the minimum value. Indicates the phase at the minimum value; When the "vertical" direction threshold is met, perform two monotonic scans along the column direction: Left sweep: from the left end position Towards the minimum position Scan, location ,…, ,like ； in, This represents the gradient magnitude at the left end of the column direction. Then let ； in, This represents the gradient magnitude in the column direction after the left sweep. Right sweep: from the right end Towards Scan: Location ,like ； in, This indicates the gradient magnitude at the right end of the column direction; Then let ； in, This represents the gradient magnitude in the column direction after the right sweep.

10. The method according to claim 1, characterized in that, Wavelet reconstruction correction is performed on the gradient component field, i.e., the subband coefficients, to obtain a new set of subband coefficients, i.e., multi-scale wavelet coefficients. The multi-scale wavelet coefficients are then combined to obtain the reconstructed image. Includes: gradient component fields, i.e., sub-band coefficients, based on S5. Perform inverse discrete wavelet transform to reconstruct the intermediate image. ;in, This represents the high-frequency horizontal components that maintain a consistent direction and continuous detail after cropping and trimming. This represents the high-frequency components in the vertical direction that maintain a consistent direction and continuous details after cropping and trimming. This represents low-frequency component information that maintains consistent direction and continuous detail after cropping and correction; the intermediate image. The grayscale range is linearly normalized, and a new set of multi-scale wavelet coefficients is obtained through discrete wavelet transform. ,in, This represents the updated high-frequency components in the horizontal direction. This represents the high-frequency components in the vertical direction after the update. This indicates the updated low-frequency component; Updated low-frequency components Replace with the original low-frequency component obtained in S2 This allows for updating only the high-frequency gradient subbands at each scale while maintaining consistency between the low-frequency structure and the initial image. and Finally, the original low-frequency components... With the updated high-frequency components , By combining the results, the final reconstructed image is obtained. .

Images