key: cord-0789112-9z2f16o5
authors: Daoui, Achraf; Karmouni, Hicham; Sayyouri, Mhamed; Qjidaa, Hassan
title: Robust 2D and 3D images zero - watermarking using dual Hahn moment invariants and Sine Cosine Algorithm
date: 2022-03-23
journal: Multimed Tools Appl
DOI: 10.1007/s11042-022-12298-0
sha: a48c2b0994e7fcd55a538d197d1d61443005f7b2
doc_id: 789112
cord_uid: 9z2f16o5

In this paper, we initially provide significant improvements on the computational aspects of dual Hahn moment invariants (DHMIs) in both 2D and 3D domains. These improvements ensure the numerical stability of DHMIs for large-size images. Then, we propose an efficient method for optimizing the local parameters of dual Hahn polynomials (DHPs) when computing DHMIs using the Sine-Cosine Algorithm (SCA). DHMIs optimized via SCA are used to propose new and robust zero-watermarking scheme applied to both 2D and 3D images. On one hand, the simulation results confirm the efficiency of the proposed computation of 2D and 3D DHMIs regarding the numerical stability and invariability. Indeed, the proposed computation method of 2D DHMIs allows to reach a relative error (RE) of the order ≈10(−10) for images of size 1024 × 1024 with an execution time improvement ratio exceeds 70% (ETIR ≥ 70%), which validates the efficiently of the proposed computation method. On the other hand, the simulation and comparison outcomes clearly demonstrate the robustness of the proposed zero-watermarking scheme against various geometric attacks (translation, rotation, scaling and combined transformations), as well as against other common 2D and 3D image processing attacks (compression, filtering, noise addition, etc.). GRAPHICAL ABSTRACT: [Image: see text]

substituting continuous orthogonal polynomials basis kernels by discrete ones such as Tchebichef [6, 9, 10, 26] , Krawtchouk [11, 26, 27] , Hahn [9, 38, 40] , Meixner [12, 22, 28, 29] , Charlier [13, 15, 19, 39] , Dual Hahn [25, 53] and Racah polynomials [4, 54] . The latter are defined on discrete intervals, and therefore no approximation or discretization is involved when using these polynomials in digital signal and image analysis.

Dual Hahn moment invariants (DHMIs) are defined from Dual Hahn polynomials (DHPs) and geometric moment invariants. DHMIs more general compared to the other types of moment invariants (with the exception of Racah moment invariants). Indeed, the discrete orthogonal Tchebichef, Meixner and Krawtchouk polynomials are derived from Dual Hahn polynomials (DHPs) by adapting the local parameters (a, b and c) of DHPs to specific cases [30] . For this reason, DHMIs are considered in the present work to propose a new and robust 2D and 3D image zero -watermarking scheme. The main advantages of the proposed scheme are: 1) the robustness against geometric and common image processing attacks, 2) numerically stable and suitable for large-size images, 3) it is a standard scheme valid for both 2D and 3D images, and 4) it offers a high level of security and robustness against various attacks. Moreover, to the best of our knowledge, this is the first time that the 3D moment invariants are used in the zero -watermarking application. The proposed scheme is mainly based on the computation of the 2D/3D DHMIs of the 2D/3D image. Next, the computed DHMIs with a binary logo (watermark) are used to generate the zero-watermark image. The latter is then used when verifying the copyrights of a transmitted image. This work also presents considerable improvements concerning the computation of both 2D and 3D DHMIs. Thanks to these improvements, the numerical stability of DHMIs is considerably improved through the elimination of all gamma and factorial functions that are required in the conventional computation of DHMIs. Moreover, the issue of the optimal selection of local DHPs parameter values when computing DHMIs is effectively addressed by using a recent population-based optimization algorithm, namely the Sine Cosine Algorithm (SCA).

The rest of the paper is organized as follows: the second section presents the theoretical background of DHPs used in DHMIs computation. In the third section, we present the computation method of both 2D and 3D DHMIs, as well as the proposed algorithms of DHMIs computation. The fourth section presents the proposed method for the optimal selection of the local DHPs parameter values. In the fifth section, we present the proposed zero -watermarking scheme based on DHMIs and SCA for 2D/3D image copyrights protection. Then, in the sixth section, the simulation and comparison results and are presented to prove the efficiency of the proposed computation of DHMIs, as well as to demonstrate the robustness of our zero-watermarking scheme compared to recent existing ones. The last section presents the conclusion and future work.

DHPs of order n and discrete variable s are defined by [53] : H n a;b;c ð Þ ðsÞ ¼ a À b þ 1 ð Þ n a þ c þ 1 ð Þ n n! 3 F 2 Àn; a À s; a þ s þ 1; a À b þ 1; a þ c þ 1; 1 ð Þ

where n = 0, 1, 2, …, N − 1, s = a, a + 1, …, b − 1, − 1/2 < a < b, |c| < a + 1, b = a + N, and 3 F 2 (.) is a generalized hypergeometric function given by:

The notation (a) k being the Pochhammer symbol that is defined from the gamma function Γ(.) as follows:

By using Eqs. (2)-(3) and the relation (−n) k = (−1) k (n) k , we can express Eq. (1) as follows: 

The symbol 〈a〉 k is the filling factorial defined as:

We notice that the expressions of the functions given by Eqs. ( 5 )- ( 7 ), depend on the factorial and gamma functions. The presence of the latter causes the problem of the numerical overflow of A n (a, b, c) , B nk (a, b, c) and F k (a) (s) values. Noting the following example of the overflow problem caused by the gamma function using Matlab:Γ(172) = Inf. To avoid this problem, A n (a, b, c) , B nk (a, b, c) and F k (a) (s) terms will be expressed independently of gamma and factorial functions. Indeed, to compute A n (a, b, c) , the following recursive formula is developed:

The recursive relation that allows to compute B nk (a, b, c) function is given by [25] :

To recursively compute the function F k (s, a), we use the following formula:

For k = 1, let's set x = s − a and a 1 = 2a + 1. Then the expression of F k (s) becomes as follows:

Using Eq. (11) for k = 2 with a 2 = 2a + 2, the following equation can be obtained:

In the general case for k ≥ 1 and a k = 2a + k, we show that F k (x, a) can be expressed as follows:

where C kr is defined as follows [25] :

where C 00 = 1, C 10 = 0, C 11 = a 1 , C 12 = 1 with a n = 2a + n. Using Eqs. (4), (9)- (10) and (14), DHPs can be expressed as follows [25] :

The DHPs satisfy an orthogonal relation of the following form [53] :

where Δx(s − 1/2) = 2s + 1 for s > − 1/2. ρ(s) and d 2 n are successively the weight and the square norm functions that are defined by Eqs.

Noting that ρ(s) and d 2 n functions depend on the gamma function, which provokes the numerical instability of the computed DHPs. We can avoid this problem by using the recursive computation of ρ(s) and d 2 n functions. Indeed, we show that ρ(s) can be expressed as:

We notice that the first term ρ(a) depends on Γ(.). Therefore, if the values of a, borc are large enough, then this term becomes numerically unstable. In order to overcome this problem, ρ(a) is computed based on logΓ(x) function as follows:

For d 2 n function is recursively computed with respect to the order n as follows:

The term d 2 0 is also computed based on logΓ(x) function as follows in order to avoid the numerical instability problem:

We notice that the expressions of A n (a, b, c) , B nk (a, b, c) , F k (a) (s), ρ(s) and d 2 n functions become independent of all gamma and factorial functions thanks to the proposed computation method. Consequently, the numerical stability of these functions is ensured when the values of n, s, a, b and c are high. As a result, the computation DHMIs of large-size images becomes possible.

The orthonormalized DHPs are computed by using the square norm and the weight functions as follows [53] :

For more details on the properties of the weighted DHPs and its recursive computation, the reader can consult the reference [53] . Figure 1 shows the orthonormalized DHPs computed up to the fifth order for different values of the parameters a, b and c.

From Fig. 1 we notice that the values of the local parameters a, b and c influence on the distribution of DHPs values, which means that these parameters can be used for the extraction of the characteristics of an image region of interest (ROI).

In this section, we briefly present the theoretical background of dual Hahn moments (DHMs) and the method of constructing dual Hahn moment invariants (DHMIs) in both 2D and 3D domains. In addition, the proposed improved algorithms for fast and stable computation of DHMIs are presented. Noting that in the rest of this article, we set x = s − a and y = t − a with s, t ∈ [a, b − 1] to compute DHMs and DHMIs of 2D image f(x, y). Furthermore, we set z = w − a and w ∈ [a, b − 1] to compute 3D DHMs and DHMIs of 3D image f(x, y, z).

For a 2D image f(x, y) of size N × M, the DHMs of order (n, m) are computed by [53] :

with n = 0, 1, 2, …, N max ≤ N − 1, m = 0, 1, 2, …, M max ≤ M − 1.

The computation time of DHMs is considerably reduced when using the following matrix formula [40] : .

The inverse transformation of the DHMs allows to obtain the reconstructed image through the following formula:

The matrix representation of Eq. (27) is given by:

wheref ¼f x; y ð Þ È É x¼N À1;y¼M À1

x;y¼0

being the reconstructed image.

a=0 , b=1024 and c=0 a=100 , b=2124 and c=50 a=200.5 , b=5225 and c=20.5 a=0, b=8000 and c=0 

This subsection presents the process used to derive DHMIs of a 2D image based on the geometric moment invariants (GMIs) and DHPs. The geometric moments (GMs) of order (n, m) are defined for a 2D image function f (x, y) of size N × M by [25] :

In order to reduce considerably the computation time of GMs, we use the following matrix formulation:

and

The central geometric moments denoted μ nm are derived from the GMs by the following formula [41] :

To improve the computation time of μ nm , we propose the use of the following matrix formulation:

The set of geometric invariant moments (GIMs) of order (n, m), denoted V nm , which are independent of translation, scaling and rotation are defined by the following relation [24] :

where ¼ nþm 2 þ 1 and θ ¼ 1 2 tan À1 2μ 11 μ 20 Àμ 02 .

In accordance with Eq. (33) , the values of the angle θ are limited to −π/4 ≤ θ ≤ π/4 . We refer the reader to see [41] for angle values in different ranges.

The computation time of V nm by Eq. (33) increases dramatically when the image size (N × M) becomes too large. In order to accelerate considerably the computation time of V nm , we propose the use of the following semi-matrix formulation: 10 GM 00 , y ¼ GM 01 GM 00 and ⊙ denotes the Hadamard product [5] .

The DHPs formula that is given by Eq. (4) allows to express DHMs (Eq. (27)) as follows [25] : (16) and (35), the DHMs are represented in terms of GMs as follows:

Finally, to compute 2D DHMIs up to the order (n, m), the expression of GMs in Eq. (36) is substituted by GMIs expression (V nm ) as follows [25] :

The conventional computation of 2D image DHMIs is summarized in Algorithm 1, and Algorithm 2 presents the proposed fast and stable computation of DHMIs based on the recursive formulas of ρ(s), d 2 n , A n (a, b, c) and B nk (a, b, c) functions, and on the matrix formulations.

For a given 3D image f(x, y, z) of size N × M × K, the DHMs of order (n, m, k) are computed as follows [33] :

The following relation is used for reconstruction a 3D image based on DHMs and DHPs:

To quantify the reconstruction errors between the reconstructed imagef x; y; z ð Þ and the original one f(x, y, z), the Mean Square Error (MSE) criterion can be used. The MSE is defined as:

The Peak Signal-to-Noise Ratio (PSNR) in decibels (dB) is another criterion that is widely used to measure the reconstruction error. PSNR is defined as:

with peak represents the maximum value of the image intensity matrix.

To derive 3D DHMIs, the 3D DHMs of the 3D image are expressed as a linear combination of 3D GMs. Then, the latter are replaced by the 3D GMIs. The 3D GMs of a 3D image f (x, y, z) of size N × M × K are defined by [18] :

To compute the central moments μ nmk , the following formula is used:

ð43Þ where x; y; z ð Þ are the coordinates of the image centroid. Noting that μ nmk are invariant with respect to translation.

The DHMIs that are invariant with respect to the rotation in the 3D domain are defined from the 3D rotation matrix associated with the Euler angle sequence. This matrix is defined by [18] :

where angle ϕ is associated with rotation on the x-axis, angle θ is associated with rotation on the y-axis, and angle ψ associated with rotation on the z-axis. It is possible to consider the 3D rotation matrix as a linear transformation of the 3D image coordinates, as follows:

where a 00 ,a 01 ,…,a 22 are coefficients of the 3D rotation matrix calculated from Eq. (45) as follows:

a 00 ¼ cos θ cos ; a 01 ¼ cos θ sin ; a 02 ¼ À sin θ a 10 ¼ sin sin θ cos À cos sin ; a 11 ¼ sin sin θ sin þ cos cos ; a 12 ¼ cos θ sin ; a 20 ¼ cos sin θ cos þ sin sin ;

By analogy with GMIs extraction method in the 2D domain, we can derive the GMIs up to (n, m, k) order in the 3D domain. These are invariants with respect to rotation, translation and scaling. The 3D GMIs are calculated by the following relation [24] :

For detailed information on the computation process of geometric moment invariants in the 3D domain, the reader can see [18] .

The use of Eqs. (16) and (38), allows to express 3D DHMs as follows [33] :

Using Eqs. (3), (33) and (37), the 3D DHMs are expressed in terms of 3D GMs as follows:

Finally, to compute DHMIs up to the order (n, m, k) in the 3D domain, we replace 3D GMs in Eq. (49) by 3D GMIs using the following formula [47] :

The conventional computation of 3D DHMIs is summarized in algorithm 3, and the proposed improved computation is presented in algorithm 4.

To quantify the invariance error between the feature vector values of the original image denoted V = DHMIs, and the feature vector values of the geometrically transformed image that is denoted V* = DHMIs*, the following Relative Error (RE) criterion can be used:

where ‖.‖ denotes the Euclidean norm. Noting that a low RE value indicates that the method used to extract the feature vector is accurate. It is worth noting that the computation of DHMIs (Algorithm 2 or 4) in both 2D and 3D domains depends on the local parameters of DHPs (a, b and c) , which influence the RE values. Therefore, the optimal selection of these parameters remains an open challenge.

The literature survey related to the image analysis by the invariant moments method shows that the choice of the local parameters is generally made by empirical methods for few cases of the local parameter values. In order to overcome this limitation, we propose to use a population-based optimization algorithm that is the Sine Cosine (SCA) for selecting the local parameters of the orthogonal polynomials during the computation of DHMIs.

The SCA is introduced by Seyedali Mirjalili [34] to solve optimization problems in various fields of engineering. This algorithm is considered as one of the most recent and efficient methods for global optimization. It is based on a mathematical model of the trigonometric functions of Sine and Cosine. SCA allows to generate and improve a set of solutions through Eqs. (52)-(53) [34] .

where V t i is a vector representing the current solution of the i-th solution in the t-th iteration, r 1 is a random vector defining the amplitude of the range of sin and cos, r 2 is a random variable used to define the range of sin or cos. The parameter r 3 represents a random value for the destination in the definition of the new position of the solution, P i represents the position of the destination point in the i-th dimension, and |.| indicates the absolute value.

SCA uses the Eqs.(52)-(53) with a 50% probability according to the following Eq. (48):

where r 4 is a random number ranging from 0 to 1. Figure 2 shows the effects of Sine and Cosine functions on Eqs. (52) and (53) . It shows how the equations describe a space between two solutions in the search area.

To stabilize exploitation and exploration, the Sine and Cosine dimensions are adaptively modified in Eqs. (52)-(54) using the following formula:

where t is the current iteration, T is the maximum number of iterations and a is a constant. In this work, we use SCA to optimize the values of the local parameters a, b and c of DHPs through the minimization of the relative error (RE), which is used as cost function of the SCA algorithm. The pseudo-code of the proposed method for optimal DHMIs computation is summarized in algorithm 5.

In this section, we present the proposed 2D/3D image zero-watermarking scheme based on 2D/3D DHMIs. This scheme is mainly designed to protect the intellectual property rights of both 2D and 3D images. The proposed scheme includes two fundamental phases, namely (i) Fig. 2 The effects of Sine and Cosine in Eqs. (52) and (53) on the following position [34] the phase of the zero-watermark image generation using a binary logo image (watermark) and the feature vector (DHMIs) that is extracted from the 2D/3D image to be protected. (ii) The second phase is the watermark recovering by using the zero-watermark image and the feature vector (DHMIs) that is computed from the attacked 2D/3D image. Therefore, the proposed method is blind, i.e. it does not require the original image. The graphical abstract shows the global scheme of the proposed zero-watermarking method, and the details are given below.

The generation of the zero-watermark image is based on the computation of DHMIs of the original image via algorithm 5. The steps of this phase are summarized in Fig. 4 , and further details are presented below.

For a 2D image, algorithm 5 is used to obtain the feature vector magnitude V = |HMI nm | of the original image up to the order (n, m). Then, we rearrange the elements of V in a square matrix of size L × L with L × L ≤ n × m. The latter will be referred as the feature matrix of the 2D image.

In the case of a 3D image, we calculate the feature vector V = |HMI nmk | of this image up to the order (n, m, k) using algorithm 5, then we rearrange V elements in a square matrix of size L × L (with L × L ≤ n × m × k). The latter referred as the feature matrix of the 3D image.

For security purpose, the values of the local parameters (a, b and c) are given as a security key denoted KEY 1. The sine and cosine functions having the dimension in [−2, 2] tolerate a solution to go beyond (outside the space between them) or around (inside the space between them) the destination [34] In this step, we rearrange the elements of the feature matrix of size L × L into a vector V c of size 1 × L 2 , then we compute the median T of this vector in order to binarize their elements as follows:

Next, the binary vector elements are rearranged in a matrix C of size L × L. The latter called the feature image.

To further increase the security of the feature image, we randomly generate a binary matrix B of the same size as C. Then, we apply the or-exclusive operation between C and B in order to obtain the permuted feature image denoted C p as follows:

where the operation "or-exclusive" satisfies the following rule:

Þwith A; B∈ 0; 1 f g ð58Þ

In this step, we use a binary logo image (watermark) of size L × L. Then, we improve the security of this image by using a pseudo-random 2D permutation based on the Arnold's cat map that is defined as [7] :

x' y'

where ε and η are positive integer parameters, (x, y) is the image pixel coordinates of the original logo and (x′, y′) is the image pixel coordinates of the permuted logo. The most important property of the chaotic Arnold's transformation is the ability to reorganize the 

To create the zero-watermark (ZW) image, we use the "or-exclusive" operation as follows:

where C p is the permuted feature image and L p is the permuted logo image.

The watermark image recovery phase allows to validate the intellectual property rights of the protected 2D/3D image. The summary of this phase is presented in Fig. 5 , and the details are given below.

For a protected 2D/3D image, KEY1 is used to compute the amplitude of DHMIs through algorithm 2 for a 2D image, or via algorithm 4 for a 3D image. Then, the feature matrix of the protected image is generated by the same method described in subsection 5.1 (a).

In this step, the feature matrix of the protected image is binarized by the same method described in subsection 5.1 (b) to obtain the feature image denoted C * .

The feature image C * with the binary matrix B that is generated in the feature image permutation step (section 5.1 (c)), are used to create the permuted feature image denoted C p * as follows: 

In this step, the permuted logo image is generated from the permuted feature image (C p * ) with the zero-watermark image (ZW) as follows:

Then we apply the inverse Arnold's transformation to L p * using KEY 2 to recover the original binary logo image (watermark). In order to measure the bit errors between the recovered watermark and the original one, we can use the Bit Error Rate (BER) criterion that is defined by:

where L × L is the binary logo size, B is the original binary logo and B * is the recovered binary logo. Noting that a low BER value indicates that the zero-watermarking method is efficient. The normalized cross correlation (NCC) criterion is also used to measure the similarity between B and B * watermarks. NCC is defined by the following relation [37] :

whereB andB * are the mean values of B and B * , respectively. Note that if the NCC value tends to the value 1, the similarity between B and B* is high. The structural similarity (SSIM) index can also be used to compute the difference between two image pixels. SSIM criterion is computed as follows [43] :

where μ, σandσ BB * represent successively the mean, deviation, and cross correlation between BandB * images. a 1 and a 2 are positive constants. If SSIM = 1, this means that there is no difference between B and B* images.

In this section, we present the results of simulations and comparisons that validate the efficiency of the present work. First, we perform a time analysis of the proposed algorithms. Then, we demonstrate the relevance of the local parameters of DHPs in the extraction of both global and local image characteristics. Next, we show the invariability property of DHMIs against geometric and noise attacks. In addition, we demonstrate the advantages of SCA in optimizing the local parameters of DHPs during the computation of DHMIs. Moreover, we show the robustness of the proposed zero-watermarking scheme against geometric and other common image processing attacks. Finally, we provide several between the proposed zerowatermarking scheme and recent existing ones to demonstrate the superiority of the proposed scheme.

Before applying 2D and 3D DHMIs in the zero-watermarking of 2D and 3D images, we show the superiority of the proposed algorithms (2 and 4) in terms of computation rapidity compared to the conventional algorithms (1 and 3) of DHMIs computation. For that, the 2D and 3D test images are resized to 64 × 64 pixels and 64 × 64 × 64 voxels, respectively (Fig. 6) to avoid the numerical instability of DHMIs due to the gamma and factorial functions involved in the conventional calculation methods (algorithms 1 and 3) . Table 1 shows a comparison in terms of execution time of 2D DHMIs using algorithms 1 and 3, and Table 3 illustrates the computation time of 3D DHMIs using algorithms 2 and 4. The improvement in the computation time is measured by the execution time improvement ratio criterion (ETIR (%)) that is defined as follows [20] :

where Time 1 represents the average of DHMIs computation time using the improved algorithms, and Time2 is the average computation time of DHMIs using the conventional algorithms. Noting that all the simulations in the present work are conducted using Matlab R2019a installed on a PC with a 2.4 GHz processor and RAM of 4 GB. On one hand, the results of the timing analysis presented in Table 1 show the proposed algorithm 3 is considerably rapid in comparison to algorithm 1. Indeed, an improvement of ETIR ≥ 70% is obtained for all the moment orders. This is explained by the use of matrix formulations to compute 2D DHMIs (algorithm 3) instead of the direct computation Fig. 6 a Standard grayscale test images resized to 64 × 64 pixels. b 3D test images resized to 64 × 64 × 64 voxels (algorithm 1). On the other hand, we observe in Table 2 a slight superiority of algorithm 4 over algorithm 2 in terms of 3D DHMIs computation time. This superiority is more noticeable in the lower orders of 3D DHMIs. We explain this result by the fact that the computation of 3D DHMs either by algorithm 2 or by algorithm 4 is performed through direct formulations, and the slight improvement of algorithm 4 is due to the recursive computation of A n (a, b, c) , B nk (a, b, c) and d 2 n instead of their direct formulations used in algorithm 2. However, the main advantage of algorithms 4 over algorithms 2 is the numerical stability. Therefore, both algorithms 2 and 4 are relevant for the computation of DHMIs of large-size images. This advantage will be exploited in the zero-watermarking application of large-size 2D and 3D images. 

In the first experiment of this subsection, we prove the invariance of 2D DHMIs against the rotation, scaling, translation, and noise adding. To do this, we select from the database [36] a real medical image named "Covid-19" of size 1024 × 1024 (Fig. 9 ). This image is then transformed by the following geometrical transformations: scale by a factor ranging from 0.5 (30, 30, 30) noise with different densities up to 5%. Figure 9 shows "Covid-19" test image with some applied attacks, and Fig. 10 shows the RE curves corresponding to the attacked test image by different transformations using various a, b and c parameter values.

The results presented in Fig. 10 show on one hand that the RE values are very low (≈10 −10 ), indicating that the DHMIs computed by the proposed algorithms are accurate and stable and "Gaussian" noise adding against geometric and noise attacks. On the other hand, we notice that the use of SCA allows obtaining the lowest RE values compared to the other cases when the local parameters of DHPs are empirically selected. Therefore, the SCA guarantees an optimal selection of the parameters a, b and c when computing DHMIs. After validating the efficiency of the proposed computation of DHMIs, we will exploit the latter in the copyright protection of large-size 2D and 3D images.

In the next test, we use a 2D color image named "Barbara" of size 1024 × 1024. Initially we compute DHMIs of this image up to the order (32, 32) using algorithm 5, and then we apply the proposed zero-watermarking scheme to generate the zero-watermark image based on the amplitude of the computed DHMIs and the binary logo image of size 32 × 32 (Fig. 11) .

Noting that the application of the proposed zero-watermarking scheme to a color image requires the decomposition of the latter into three channels R, G and B. Then, we apply the proposed scheme to a single channel that is considered as host image. In the present test, the Gchannel is used as host image. Then, "Barbara" image is attacked by 10 transformations illustrated in Fig. 11 . The recovered watermarks from the attacked images are shown in Fig. 13 with the corresponding BER, SSIM and NCC values. The results in Fig. 12 show that the BER values are very low (close to zero), and the SSIM and NCC values tend towards 1, which clearly indicates that the proposed method is of high efficiency for the zero-watermarking of large-size 2D images.

In the following test, we use the 3D "Vertebrae" test image of size 256 × 256 × 256. Next, algorithm 5 is used to compute DHMIs of this image up to the order (11, 11, 11) . Then, 1024 coefficients are randomly selected from the computed DHMIs, and the selected coefficients are rearranged in 2D matrix of size 32 × 32 (Fig. 13) for applying the proposed zerowatermarking scheme. The test image is attacked by 12 transformations including geometric attacks (rotation, translation, scaling and combined transformations), contamination by "Salt & Papper" and "Speckle" noise, cropping, compression, and filtering attacks. The attacked 3D image by various transformations is shown in Fig. 13 . Then, the proposed scheme is then used to recover the watermark from each attacked image, and the recovered watermarks with the corresponding BER, SSIM and NCC values are illustrated in Fig. 14 .

The achieved results presented in Fig. 14 

In this section, we present comparisons of the proposed scheme with recent zero-watermarking schemes of 2D and 3D images. In the first test, we compare the robustness of the proposed 2D zero-watermarking scheme with recent schemes presented in [24, [45] [46] [47] 51] . In the comparison experiments, we use four test images each of size 1024 × 1024 with a binary logo of size 32 × 32 (Fig. 15 ). It should be noted that the robustness of the test schemes is measured using Table 3 . The results illustrated in Table 3 indicate that the average BER values obtained when using the proposed scheme are lower than those obtained by the compared schemes, indicating that the proposed scheme presents excellent performance in comparison to the other schemes. This superiority can be justified by the strong robustness of DHMIs against both geometric and common image processing attacks. In addition, the use of SCA when computing DHMIs ensures the extraction of an optimal and stable feature vector, thus the SCA improves the performance of the proposed zero-watermarking scheme.

The 3D image zero-watermarking schemes remain relatively few. Therefore, the robustness of the proposed scheme is only compared to the latest three schemes presented in [31, 32, 42] . To conduct this comparison, we use six 3D images each of size 256 × 256 × 256 (Fig. 16) . These images are attacked by various transformations given in Table 4 . Then, the average BER is computed using on the recovered watermark from each attacked image.

The results of this comparison (Table 4) show that the proposed method performs better than the methods presented in [31, 32, 42] because the average BER values obtained by the proposed scheme are always lower than those obtained by the scheme used in this comparison. The superiority of the proposed method over the other compared schemes is justified by the fact that the 3D DHMIs are invariants with respect to the geometric transformations of the 3D image (translation, rotation and scale change). By contrast, the other compared 3D zerowatermarking schemes are influenced by the geometric attacks of the 3D image. In addition, the lower order 3D DHMIs are stable and robust against different types of noise and common 3D image processing attacks. In this paper, we have improved considerably the computation of 2D and 3D DHMIs. Then, SCA is efficiently used for the optimization of local DHPs parameters when computing DHMIs. These are used to propose a new and robust zero -watermarking scheme used in both 2D and 3D images copyrights protection. The simulation results and comparisons have shown that the numerical stability of DHMIs is significantly improved due to provided computational improvements. On the other hand, the high robustness of the proposed zerowatermarking scheme against different types of geometric and common attacks of 2D/3D images has been demonstrated. In the future, the proposed algorithms of 2D/3D DHMIs computation will be applied in other applications that are based on the invariance property, such as pattern recognition and classification of large-size images.

Chaos-based robust method of zerowatermarking for medical signals

New set of generalized legendre moment invariants for pattern recognition

Fractional-order orthogonal Chebyshev moments and moment invariants for image representation and pattern recognition

Fast and accurate computation of Racah moment invariants for image classification

Hadamard products of linear spaces

Some computational aspects of Tchebichef moments for higher orders

A symmetric image encryption scheme based on 3D chaotic cat maps

Biomedical signals reconstruction and zero-watermarking using separable fractional order Charlier-Krawtchouk transformation and sine cosine algorithm (SCA)

New algorithm for large-sized 2D and 3D image reconstruction using higher-order Hahn moments

Efficient Reconstruction and Compression of Large Size ECG Signal by Tchebichef Moments

Fast and stable bio-signals reconstruction using Krawtchouk moments

Efficient computation of high-order Meixner moments for large-size signals and images analysis

Stable computation of higher order Charlier moments for signal and image reconstruction

New robust method for image copyright protection using histogram features and sine cosine algorithm

Efficient methods for signal processing using Charlier moments and artificial bee Colony algorithm. Circuits, Systems, and Signal Processing

A zero watermarking algorithm for 3D mesh models based on shape diameter func-tion

New set of fractional-order generalized Laguerre moment invariants for pattern recognition

Image classification using separable invariant moments of Charlier-Meixner and support vector machine

Fast computation of accurate Zernike moments

New fractional-order Legendre-Fourier moments for pattern recognition applications

Image analysis by Meixner moments and a digital filter

A robust zero-watermarking algorithm for color image based on tensor mode expansion

Combining polar harmonic transforms and 2D compound chaotic map for distinguishable and robust color image zero-watermarking algorithm

Generalized dual Hahn moment invariants

Image analysis using separable Krawtchouk-Tchebichef's moments

Image reconstruction by Krawtchouk moments via digital filter

Fast computation of inverse Meixner moments transform using Clenshaw's formula

Fast computation of 3D Meixner's invariant moments using 3D image cuboid representation for 3D image classification

Hypergeometric orthogonal polynomials and their qanalogues

Robust 3D mesh zero-watermarking based on spherical coordinate and skewness measurement

Zero-watermarking method for resisting rotation attacks in 3D models

Robust reconstruction and generalized dual Hahn moments invariants extraction for 3D images

SCA: a sine cosine algorithm for solving optimization problems

Robot visual guide with Fourier-Mellin based visual tracking

Radiopaedia.org, the wiki-based collaborative Radiology resource

Simple low-dimensional features approximating NCC-based image matching

A fast computation of Hahn moments for binary and gray-scale images

A fast computation of charlier moments for binary and gray-scale images

Improving the performance of image classification by Hahn moment invariants

Image analysis via the general theory of moments

A zero-watermarking scheme for three-dimensional mesh models based on multifeatures

Image quality assessment: from error visibility to structural similarity

Based zero-watermark digital watermarking technology, in: the third China information hiding and multimedia security workshop (CIHW)

Local quaternion polar harmonic Fourier moments-based multiple zero-watermarking scheme for color medical images

A robust zero-watermarking algorithm for lossless copyright protection of medical images

Novel quaternion polar complex exponential transform and its application in color image zero-watermarking

Zero-watermarking scheme for 3D meshes based on geometric property

Fractional Charlier moments for image reconstruction and image watermarking

Robust zero-watermarking scheme based on novel quaternion radial fractional Charlier moments

Color image zero-watermarking based on fast quaternion generic polar complex exponential transform

Seismic signal matching and complex noise suppression by Zernike moments and trilateral weighted sparse coding

Image analysis by discrete orthogonal dual Hahn moments

Image analysis by discrete orthogonal Racah moments

Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.