Efficient and Stable Voxel-Based Algorithm for Computing 3D Zernike Moments and Shape Reconstruction

3D Zernike moments based on 3D Zernike polynomials have been successfully applied to the field of voxelized 3D shape retrieval and have attracted more attention in biomedical image processing. As the order of 3D Zernike moments increases, both computational efficiency and numerical accuracy decrease. Due to this phenomenon, a more efficient and stable method for computing high-order 3D Zernike moments was proposed in this study. The proposed recursive formula for computing 3D Zernike radial polynomials combines the recursive calculation of spherical harmonics to develop a voxel-based algorithm for the calculation of 3D Zernike moments. The algorithm was applied to the 3D shape Michelangelo's David with a size of 150×150×150 voxels. As compared to the method without additional acceleration, the proposed method uses a group action of order sixteen orthogonal group and saving unnecessary iterations, the factor of speed-up is 56.783±3.999 when the order of Zernike moments is between 10 and 450. The proposed method also obtained an accurate reconstructed shape with the error rate (normalized mean square error) of 0.00 (4.17×10 -3 ) when the reconstruction was computed for all moments up to order 450.


Introduction
The 3D Zernike moments as shape descriptor based on 3D Zernike polynomials have been widely applied in the many fields of 3D image processing [1,2], such as terrain matching [3,4], shape retrievals [1,2,5], human motion analysis [6,7], and proteinprotein interface prediction [8], etc. Due to the orthogonal property of 3D Zernike polynomials, the resultant Zernike moments do not contain redundancy or overlap of information between the moments, and the magnitudes of these moments are invariant to the rotation of the represented 3D shape.
The 3D Zernike moments can be traced back to the work of Canterakis [9,10].
The 3D moments are expressed by a complex number linear combination of geometric moments, which are widely used in image processing and pattern recognition. In [11][12][13], the algorithms for calculating geometric moments was proposed by using exact integration in each voxel. Based on the theory in [1,2], Novotni and Klein developed the algorithm for calculating Zernike-Canterakis moments. The algorithm has been commonly adopted to calculate 3D Zernike moments. A reformulation of Zernike-Canterakis moments and its implementation were provided in the article by Callahan and De Graef [11]. The coefficients of geometric moments range from 10 0 to10 9 for the maximal order 20 and the geometric moments also have a great dynamic range [12].
This method of computing 3D Zernike moments is subject to severe instability of the geometric moment in the case of higher orders [13]. On the other hand, [14,15] provided the algorithms for calculating geometric moments based on the exact volumetric integrations for surface representations. The reconstructed images with high resolution from the moments up to order 40 were shown in [14]. These surface-based methods calculate 3D Zernike moments more effectively and accurately than Novotni and Klein's method.
Geometric moments using exact integration in each pixel are also used for computing 2D Zernike moments. However, when the moment order is greater than 45, this method of calculating Zernike moments cannot provide numerical stability [16].
Another approach is to recursively calculate 2D Zernike radial polynomials. Even if the zeroth-order approximation is used for the integration, these methods are stable and efficient [16,17]. Based on the recursive approach of 3D Zernike radial polynomials, we proposed an algorithm for calculating 3D Zernike moments. 3D Zernike polynomials are the product of 3D Zernike radial polynomials and the spherical harmonics. Spherical harmonics are generally suitable for surface modeling [18], which are the base functions of the 2 functions on the unit sphere S 2 [19]. If the surface of the 3D shape is compact and can be one-to-one projected onto a sphere (e.g., a star shaped surface), it is sufficient to describe it using the Fourier coefficients of spherical harmonics. However, the surface of such a 3D object is topologically identical to the unit sphere S 2 [20]. In the real world, not all surfaces of 3D objects are homeomorphic to S 2 , for example, the torus. As described in [1], multiple concentric shells were used to define spherical functions, and spherical harmonics are not suitable for the torus.
In this study, we introduced a series of stable algorithms to compute high-order 3D Zernike moments, which is an extension of our work in [21]. A recursive formula analogue to Kintner's p-method [22] was proposed for calculating 3D Zernike radial polynomials. Based on the recursion, in combination with the calculation of the spherical harmonics, a voxel-based algorithms for computing 3D Zernike moments were developed. The use of the symmetries generated by group actions effectively reduces the calculations in the proposed algorithm. The group theory was comprehensively introduced in [23], where groups of order 4 or 8 were used to calculate 2D Zernike moments [16,17], and an abelian group was suggested in the computation of 3D Zernike moments [24]. Orthogonal groups were also used to further reduce the amount of calculation. The concept of group actions was used to introduce the definition of shape symmetry coefficient for a 3D shape, which resulted in better acceleration effects for larger coefficients. The proposed methods not only provide an accurate reconstructed voxelized 3D shape but is also very effective and stable when calculating using group action under a matrix Lie group of order sixteen.
The 3D Zernike polynomials form an orthogonal basis for the Hilbert space 2 ( ) of square-integrable complex functions over the unit ball . The projection coefficients of the 3D image function onto the 3D Zernike polynomials are called 3D Zernike moments. The low-order Zernike moments describe the global contour of a shape and the higher-order moments describe the local detail. Low-order shape descriptor provides sufficient 3D shape information for the search and classification of objects; however, the more detailed shape features are required to differentiate the similar shape objects. The proposed method was applied to voxelized 3D shapes having a size of 150×150×150 voxels with the maximal moment order up to 450. The simulation results have shown that the proposed methods are fast and stable as compared to the methods provided in the literature.
The rest of this article is organized as follows. Section 2, a brief introduction to 2D Zernike and 3D Zernike moments were given. Some theories for calculating 3D Zernike moments have also been proposed with proofs, including a summary of spherical harmonic calculations. In Section 3, two kinds of error measurements for the reconstructed voxelized 3D shape were presented. The proposed algorithms for computing 3D Zernike moments were listed in Section 4. Based on the group action of some specific orthogonal groups, the group theory was applied to accelerate the calculation of 3D Zernike moments in Section 5. A finite group of order sixteen was used to speed up the algorithm by a factor of sixteen. To understand the relationship more clearly between acceleration effects and voxelized 3D shapes, the definition of acceleration factor and the shape symmetry coefficient of 3D shape were given in Section 6. In Section 7, we summarized the experimental simulations with numerical results. In section 8, a brief conclusion was given.

2D Zernike and 3D Zernike Moments
3D Shape analysis is getting more and more attention in the medical community because it can potentially evaluate morphological changes between healthy and pathological structures. Before introducing the 3D Zernike transformation, let us briefly review the case of 2D. The 2D version of Zernike polynomials were introduced by the Nobel prize winner F. Zernike, which form an orthogonal basis for the Hilbert space 2 ( ) of square integrable complex functions over the unit disc D [25]. By the definition of Teague [26], the coefficients of the image function projected onto the Zernike polynomials are called Zernike moments. The unit disc D is the topological product of the closed unit interval I = [0, 1] and the unit circle 1 . For a non-negative integer (called order) and an integer (called repetition) for which | | ≤ and ≡ (mod 2), the 2D Zernike polynomials ( ) is the product of 2D Zernike radial polynomial 2 ( ) defined over I and the complex exponential function e = cos + sin defined over the unit circle 1 , where = ∈ stands for a complex number, is its distance to the origin (0, 0) and is the angle from the x-axis. The computations of 2D Zernike polynomials, Zernike moments and their variants were comprehensively studied in [17,[27][28][29][30].
The complex exponentials functions e ( = ⋯ , −3, −2, −1, 0, 1, 2, 3, … ) build an orthonormal basis for the Hilbert space 2 ( 1 ). In the case of 2 , the spherical harmonics ℓ (•) form an orthonormal basis for the Hilbert space 2 ( 2 ) of the square-integrable functions over the unit sphere 2 [31,32]. Any function of 2 ( 2 ) can be expressed as follows: where ℓ are the coefficients. For a non-negative integer ℓ and an arbitrary integer m with | | ≤ ℓ, the spherical harmonic ℓ (•) of degree ℓ with order is given by where ℓ (•) is the associated Legendre polynomial of degree ℓ, ℓ is the normalizing factor given by The orthogonal properties hold in Eq. (5) and (6), thus the spherical harmonics ℓ (•) are mutually orthogonal.
where ′ is the Kronecker delta.
Theorem B. The notations as stated above. It holds where ′ is the Kronecker delta.
It is worth noting that due to the orthogonal properties stated in Eq. (6) and Eq. (13), a classic argument for the orthogonal polynomials (Cf. Theorem 3.2.1 in [31]) can be used to derive the three-term recurrence formulas of the associated Legendre polynomials and the one for 3D Zernike radial polynomials described in Theorem A and Theorem D, respectively. The definition of the radial polynomial given in (12) looks different from that of Canterakis [9,10], but these two definitions are equivalent [34].
After expanding the factors in Eq. (12) and performing some algebraic calculations, a representation of the 3D Zernike polynomial is obtained as shown in Eq. (16).
Due to the numerical accuracy consideration, it is recommended to use a recursive approach instead of computing the Zernike polynomial stated in Eq. (16) directly. Based on the recurrence of the Jacobi polynomials [31], the recursive approach is provided.
This 3D recursive formula can be regarded as similar to the Kintner's p-method used in the case of 2D Zernike polynomials [22,30]. For the completeness, we provide a proof.
To apply this recursive formula, the following initial relations are also required: and , −2 ( ) = ( + Let ( , , ) be a 3D image function within the unit ball B whose values lie over the closed interval [0, 1]. The 3D Zernike moment ℓ can be regarded as the inner product of the image function ( , , ) with the basis function ℓ ( , , ), which is given by Eq. (24).
where ̅ denotes the complex conjugate of . Since the skew symmetry property of 3D Zernike polynomials ℓ ( , , ) is inherited from spherical harmonics as Eq. (9), so is the 3D Zernike moments All moments with the same order n and same degree ℓ build (2ℓ 11)-dimensional vector ⃑ ℓ given by The l 2 -norm of ⃑ ℓ , denoted by is rotationally invariant and represents the shape descriptor of a voxelized 3D image.

Error Measurements for Reconstruction
Let the function be defined as Since the 3D image function is piecewise constant, and hence squared integrable.
The function is convergent to in the sense of L 2 -norm as tends to infinity.
For large is considered to be the reconstructed function ̂ of the image function . The normalized mean square error (NMSE) between the original image and its reconstructed image, shown in Eq. (29), is used to measure the computational stability and accuracy of the proposed method in calculation 3D Zernike moments.
This measurement is well suited for the image function carrying different intensity information. However, the voxelized image is the characteristic function for a 3D object ⊂ the value of ( , , ) ∈ is , otherwise . For a given moment order , the reconstructed function ̅ is defined as follows: So, the error measurement called the relative recovery error rate ̅ is defined by

Algorithms for computing 3D Zernike Moments
The algorithms here we proposed are all voxel-based. For simplicity, we assume that is even, if is odd, the following discussion needs a slight modification with replaced by + 1 in (32). Suppose that the entire 3D image is inside the cube with edge length = 2/(√3 ) , and ( , , ) denotes its corresponding spherical coordinates. Of note that for different ≠ ′, on the plane = (or = ′ ), the angles from the x-axis and ′ are the same which can be obtained by the following equation, Therefore, the spherical coordinate system can be rewritten as (  ,  ,  ). Let the image function be piecewise constant for which has the same value on the same voxel. Hence, the 3D Zernike moment ℓ can be approximated by the following discrete form: where 0 ≤ , , ≤ − 1. The algorithms for computing all 3D Zernike moments The calculation complexity of (3) in Algorithm A is ( 3 3 ) and takes the most computational time, which executes sequentially. In the above algorithms, the look-up tables for coefficients, trigonometric functions and exponentials can be built in advance, and the algorithm can be revised slightly to execute in parallel.

Computation Speed-up via Group Action
In this section, we introduced the method to speed up the algorithms described previously by using group action. In [23], some concepts of group theory for this To explicit the use of group action by the matrix group G 16 , we summarize some useful properties in the following theorem; most of them can be verified directly by matrix multiplication calculation.
Theorem E. The notations as stated above. The following properties hold: 1. G 16 ( ) ⊂ C, i.e. for any ∈ G 16 and for any voxel in C, it holds ( ) ∈ C. Proof. A proof sketch for (3) of Theorem E was given. By direct matrix multiplication, it follows that 2 = 2 = 2 = 2 = Let G ′ be the group generated by those four involutions. By (2) of Theorem E, it concludes that these four involutions are elements of G 16 , so G ′ ⊆ G 16 . Since both groups have the same group order, G ′ = G 16 .

Q.E.D.
The cube C can be equally divided into sixteen parts. In Figure 1, the H is defined as follows: The visual representation of the G16-orbit of a voxel inside the H is shown in Figure 2. Each voxel in the set has its associated stabilizer 〈 〉 that is a group of order two, and its G16-orbit is composed of 8 distinct points: The computation of the exponential sum has only four cases to be considered with ≡ 0, 1, 2, 3 (mod 4). By the use of Table 1, the calculations are listed in Table 2.
It is larger than the set . To measure the degree of shape symmetry under the group action by G, we define the G-symmetry coefficient as follows: There is one-to-one correspondence between the set J (≡ J H 0 ∪ J D ) and the set of G16orbits in G 16 ( ). Thus, the acceleration factor under G16 is given by The ACC( ) value depends on the shape of 3D image and is slightly smaller than the 16 × SC(G 16 , ) value that is the theoretical value of the acceleration factor of the G16-conjugation applied to the Algorithm A'.

Experimental Setup
These experiments were performed on a personal computer with 32 GB of RAM, 3  Table 3 with parameters whose definition were described in the previous sections. In the implementation of Algorithm A' with the acceleration of G16conjugation, the index set J, or equivalently the set of G16-orbits in G 16 ( ) , is constructed before calculating Eq. (59) to run the algorithm smoothly.

Numerical Accuracy
There are four methods for calculating 3D Zernike moments with different settings for the performance comparison are listed in Table 4. With the proposed method, when the maximum moment order becomes larger (up to 450), the NMSE (recovery error rate) of the shape reconstruction becomes smaller rapidly for the 3D test images as show in Figure 4. We also used algorithm provided by Novotni and Klein in [1,2] and applied it to the image Michelangelo's David. As reported in [12], when the maximum order is 20, the absolute value of the geometric moment coefficient ranges from 10 0 to 10 9 , and the geometric moments also have a great dynamic range. The calculation error accumulates very quickly and is directly reflected in the shape reconstruction at = 20, as shown in Figure 3. The result of the shape reconstruction has many connected components, most of which are caused by computational errors except for the middle one. When the order is greater than 20, the algorithm reconstructs the shape is very time-consuming. To evaluate the performance of numerical precision, we first calculated the 3D Zernike moments for the 3D voxelized images in the experiments and then used those moments to reconstruct the 3D shapes. The two types of differences between the original shape and the reconstructed, ε 2 (Eq. (29)) and ̅ (Eq. (30)), are calculated to assess the numerical accuracy of the reconstruction. Figure 4 shows that the NMSE ε 2 for the shape reconstructions of Angel, Buddha, and Michelangelo's David 3D test images have similar error behaviors. The values of ε 2 of the shape reconstructions keep decreasing up to the maximum order = 450, which are 0.003831291, 0.003683659, and 0.0041664, respectively.    Table 6 shows three different angle views of reconstructed shapes using different Zernike moment orders. As mentioned before, the low-order Zernike moments describe the global contour of the shape; the higher-order moments describe the local detail.

Computation Speed
The methods A0 and A16 are performed for every voxel ( , , ) with center ( , , ) inside C ( H ). Therefore, the elapsed times for all three test 3D images differ only slightly. Table 7 only shows the elapsed times of applying the methods A0 and A16 to the voxelized 3D Angel image. The value of speed-up entitled A0/A16 and displayed in the fourth column in Table 7 is between 14.125 and 16.544, which is very close to the theoretical value 16. Unit of elapsed times: second.
The elapsed time of method Anz and method Anz,16 is proportional to the numbers of non-zero value voxels and the value of J, respectively. The experimental results for the methods Anz and Anz,16 applied to the three images are listed in Table 8 and Table 9 respectively. For different Zernike moment orders, the ratios of the elapsed time over the non-zero value voxels (elapsed time/|J|) shown in the same row in Table 8 (Table 9) are only slightly different. Since the values of NZV and |J| are equivalent to the total number of iterations processed in the two algorithms respectively, the elapsed times of these two methods are highly related to the total number of iterations in the algorithm.  From the experimental results in Table 9, as long as the |J| value of both the reference image and the query image is known in advance, we can estimate the elapsed time of the Zernike moments calculation for the query 3D image based on the elapsed time of calculating Zernike moments taken by the reference image for a given maximal order. The elapsed time of the method Anz,16 applied to the image Angel at the maximal order 450 takes 1638.303 seconds. The calculation of the proposed method outperforms other methods in terms of speed and accuracy in the literature.
Comparing the speed-up of method Anz,16 relative to method Anz with the acceleration factor calculated in Table 3, we list the ratio of the elapsed time of method Anz over the one of Anz,16 in Table 10. The experiment shows that the speed-up results are not far from the expected values.  result shows excellent improvement in the calculation. In particular, 3D images have a small number of non-zero voxels and have a large ACC value. This type of image will make the |J| value smaller, which will speed up the 3D Zernike moments calculation.

Conclusions
The proposed method can achieve accurate numerical results, especially for high-order moments. Applying the group actions of the finite orthogonal group G16 can significantly reduce the computational cost and accelerate the calculation. The experimental results showed that the proposed method Anz,16 (Algorithm A' and using G16-conjugation) took 896.056 seconds to compute the top 450-order 3D Zernike moments of Michelangelo's David with a size of 150×150×150 voxels. If the shape is reconstructed using the Zernike moments up to order = 450, the normalized mean square error (NMSE) ε 2 is 0.0041664, and the recovery error rate ̅ is 0. When compared to the proposed method A0 (Algorithm A without any acceleration), the acceleration factor of applying method ANZ16 to Michelangelo's David is between 50.219 and 61.822 when calculating moments at the maximal order between 10 and 450.
This proposed method is superior to other compared voxel-based methods in terms of speed and accuracy when calculating high-order 3D Zernike moments. The algorithm is also suitable for applying 3D Zernike transform to 3D objects containing voxel intensities. For the future, the application of this analysis would be explored, for example, potentially assessing the morphological changes of biological structures between healthy and pathology, including not only the shape but also texture.