An Algorithm for Fast Computation of 3D Zernike Moments for Volumetric Images

Abstract

An algorithm was proposed for very fast and low-complexity computation of three-dimensional Zernike moments. The 3D Zernike moments were expressed in terms of exact 3D geometric moments where the later are computed exactly through the mathematical integration of the monomial terms over the digital image/object voxels. A new symmetry-based method was proposed to compute 3D Zernike moments with 87% reduction in the computational complexity. A fast 1D cascade algorithm was also employed to add more complexity reduction. The comparison with existing methods was performed, where the numerical experiments and the complexity analysis ensured the efficiency of the proposed method especially with image and objects of large sizes.

1. Introduction

Moments of images are generally defined as projections of the image function onto a set of basis functions. In his pioneer work, Hu [1] popularized the usage of image moments in 2D pattern recognition. The set of Hu's moments gain the interest of scientists and widely applied during the last five decades. Teague [2] suggested the usage of orthogonal basis functions such as Legendre and Zernike to construct moments. Orthogonal moments are used to represent images with a minimum amount of information redundancy. Teh and Chin [3] evaluated different orthogonal and nonorthogonal moments. They found 2D Zernike moments superior over others moments in terms of noise sensitivity and discrimination power. 2D Zernike moments are used in a wide range of applications such as pattern recognition applications [4–6], content-based image retrieval [7–9], image watermarking [10–12], biometrics [13, 14], analysis and recognition of medical images [15–17], and edge detection [18].

Reconstruction, recognition, and discrimination of 3D objects gained more research interest during the last decade. 3D moment invariants are used as shape descriptors [19]. The superiority of orthogonal 2D Zernike moments over the nonorthogonal moments motivate Canterakis [20] to generalize the classical 2D Zernike polynomials to 3D. Canterakis paid his attention to theoretical aspects of deriving 3D Zernike polynomials and moments. Novotni and Klein [21] addressed 3D Zernike descriptors for content-based shape retrieval where the information content of the recovered 3D shape has no redundancy because of the orthonormality. In addition to this property, a group of 3D Zernike moments are rotation-invariant and the shape reconstruction from 3D Zernike is a very simple process. Also, 3D Zernike moments have the advantage of capturing global information about the 3D shape without requiring closed boundaries as in boundary-based methods.

Due to the attractive properties, applications of 3D Zernike descriptor gain more interest. Millán et al. [22] used 3D Zernike moment invariants in morphological characterization of intracranial aneurysms. Qiuting and Bing [23] used the 3D Zernike moments to derive a 3D terrain matching algorithm. Since, shape plays a crucial role in molecular recognition and function. The development of shape analysis techniques is important for understanding protein structure-function relationships. Recently, 3D Zernike moments gain more interest from peoples working in the area of bioinformatics and molecular biology. Sael et al. [24, 25] applied 3D Zernike moments for protein tertiary structure retrieval and comparison of properties on protein surface. 3D Zernike moments are promising descriptors in the field of biological imaging. X-ray diffraction and electronic microscopic imaging are examples where these 3D descriptors could be used to build up 3D views of biological entities such as proteins, nucleic acids, cells, tumors, tissues and whole organs or organisms.

Unfortunately, high computational demands of 3D Zernike moments hindered the wide applications of these applications. The need for computational approaches to efficiently store, display, and compare this data is the motivation of the present work. In the literature, there are many methods and approaches for fast computation of 2D Zrnike moments. The extendibility of these methods to compute 3D Zernike moments will be discussed.

In the first approach, Papakostas et al. [26] proposed a modified direct method for the computation of the Zernike moments. They computed factorial terms by using Sirling's formula. This method cannot be extended to compute 3D Zernike moments. Using the recurrence relations is another approach that should be consider. Al-Rawi [27] proposed an algorithm for fast computation of 2D Zernike moments. Unfortunately, this approach cannot be extended to compute 3D Zernike moments.

Exact computation of 2D Zernike moments via exact 2D geometric moments is the third approach. Wee and Paramesran [28] proposed a novel method in which they expressed 2D Zernike moments as expansion of exact 2D geometric moments of the same order. In fact, the method of Chong is highly accurate and very time consuming. Hosny [29] modified the method of Chong to be a very fast method. Recently, Hosny [30] proposed a novel method for fast and accurate computation of the full and subsets 2D Zernike moments. Fortunately, this approach is extendable to compute 3D Zernike moments.

In this work, a fast, low-complexity method is proposed for efficient computation of 3D Zernike moments. The entire set of 3D Zernike moments and the selected set of rotationally invariant 3D Zernike moments are computed as a combination of exact 3D geometric moments. A 3D symmetry-based method is applied where 87% of the computational complexity is reduced. A fast algorithm is applied to accelerate the computational process. The heavy computational expensive binomial coefficients terms are avoided by using a very simple computational method. The proposed method significantly reduces the whole computational complexity and seems to be more suitable for large objects and databases.

2. Three-Dimensional Zernike Moments

Three-dimensional Zernike polynomials, 𝑍_{𝑛,ℓ,𝑚}, are orthogonal polynomials defined on a unit ball as follows: 𝑍_{𝑛,ℓ,𝑚}(ℜ)=𝑅_𝑛,ℓ(𝑟)𝑌_ℓ,𝑚𝜃,𝜙, (2.1) where 𝑛∈[0,Max], ℓ∈[0,𝑛], and 𝑚∈[−ℓ,ℓ]. The value (𝑛−ℓ) must be even nonnegative integer number. Max is the maximum order considered in the computational process. 𝑅_𝑛,ℓ(𝑟) and 𝑌_ℓ,𝑚(𝜃,𝜙) are real-valued radial functions and spherical Harmonics, respectively. Any function, 𝑓(ℜ), defined within the unit ball could be expanded by using the 3D Zernike polynomials as follows: 𝑓(ℜ)=^∞𝑛_𝑛=0^ℓ_{ℓ=0𝑚=−ℓ}Ω_{𝑛,ℓ,𝑚}𝑍_{𝑛,ℓ,𝑚}(ℜ). (2.2) The expansion coefficients, Ω_{𝑛,ℓ,𝑚}, represent the 3D Zernike moments and determined by using the complex conjugate of 3D Zernike polynomials as follows: Ω_{𝑛,ℓ,𝑚}=¹₀₀^2𝜋𝜋₀𝑍_{𝑛,ℓ,𝑚}(ℜ)𝑓(ℜ)𝑟²sin𝜃𝑑𝑟𝑑𝜃𝑑𝜙. (2.3) According to the extreme complexity of computation in the 3D spherical coordinates, Canterakis [20] formulated the 3D Zernike polynomials in the Cartesian coordinates where the conversion between spherical and Cartesian coordinates, [𝑋]^𝑇=[ℜ]^𝑇, is defined in an explicit form as follows: 𝑥𝑦𝑧=𝑟sin𝜃sin𝜙𝑟sin𝜃cos𝜙𝑟cos𝜙. (2.4) The 3D Zernike polynomials in the Cartesian coordinates are defined as follows: 𝑍_{𝑛,ℓ,𝑚}(𝑋)=^𝑘_𝑣=0𝑄_{𝑘,ℓ,𝑣}|𝑋|^2𝑣𝑒_ℓ,𝑚(𝑋), (2.5) where 𝑘 is an integer equal to (𝑛−ℓ)/2 and 0≤𝑣≤𝑘. The coefficients 𝑄_{𝑘,ℓ,𝑣} are defined as 𝑄_{𝑘,ℓ,𝑣}=(−1)^𝑘2^2𝑘2ℓ+4𝑘+33𝑘()2𝑘−1^{𝑣𝑘𝑣2()𝑘+ℓ+𝑣+12𝑘}_{𝑘𝑘+ℓ+𝑣}. (2.6) The harmonic polynomials, 𝑒_ℓ,𝑚(𝑋), are defined as 𝑒_ℓ,𝑚(𝑋)=𝐶_ℓ,𝑚𝑟^ℓ𝑖𝑥̂−𝑦2^𝑚𝑧^{ℓ−𝑚⌊(ℓ−𝑚)/2⌋}_𝜇=0ℓ𝜇−ℓ−𝜇𝑥𝑚+𝜇²+𝑦²4𝑧^2𝜇, (2.7) where 𝑖̂=√−1, 𝑧=𝑥+𝑖𝑦̂ is the complex variable and 𝐶_ℓ,𝑚 are the normalization factors and defined as 𝐶_ℓ,𝑚=(2ℓ+1)(ℓ+𝑚)!(ℓ−𝑚)!.ℓ! (2.8) Based on (2.8), the normalization factors 𝐶_ℓ,𝑚 of positive and negative values of 𝑚 are identical, 𝐶_ℓ,−𝑚=𝐶_ℓ,𝑚 while the harmonic polynomials with negative values of 𝑚 are defined as 𝑒_ℓ,−𝑚(𝑋)=(−1)^𝑚𝑒_ℓ,𝑚(𝑋). (2.9) The orthogonality relation of 3D Zernike polynomials is defined as follows: 34𝜋_|𝑋|≤1𝑍_{𝑛,ℓ,𝑚}(𝑋)𝑍_{𝑛′,ℓ′,𝑚′}(𝑋)𝑑𝑋=𝛿_𝑛,𝑛′𝛿_ℓ,ℓ′𝛿_𝑚,𝑚′. (2.10) The 3D Zernike moments and the image/object reconstruction in spherical coordinates as defined by (2.3) and (2.2) are converted into the Cartesian coordinates and rewritten as follows: Ω_{𝑛,ℓ,𝑚}=34𝜋_|𝑋|≤1𝑓(𝑋)𝑍_{𝑛,ℓ,𝑚}(𝑋)𝑓(𝑋)=𝑑𝑋_{∞𝑛𝑛=0ℓℓ=0𝑚=−ℓ}Ω_{𝑛,ℓ,𝑚}𝑍_{𝑛,ℓ,𝑚}(𝑋), (2.11) where the 3D Zernike moments with negative values of 𝑚 could be computed directly from their corresponding ones with positive values of 𝑚 as follows: Ω_{𝑛,ℓ,−𝑚}(𝑋)=(−1)^𝑚Ω_{𝑛,ℓ,𝑚}(𝑋). (2.12) Canterkis [20] compactly formulate the 3D Zernike polynomials of order 𝑛 as a linear combination of monomials of order up to 𝑛 as follows: 𝑍_{𝑛,ℓ,𝑚}(𝑋)=_{𝑟+𝑠+𝑡≤𝑛}𝐻^𝑟𝑠𝑡_{𝑛,ℓ,𝑚}𝑥^𝑟𝑦^𝑠𝑥^𝑡. (2.13) Consequently, the 3D Zernike moments Ω_{𝑛,ℓ,𝑚} of order 𝑛 could be expressed as a linear combination of the 3D geometric moments, 𝐺_{𝑟,𝑠,𝑡}, by using the following relation: Ω_{𝑛,ℓ,𝑚}=34𝜋_{𝑟+𝑠+𝑡≤𝑛}𝐻^𝑟𝑠𝑡_{𝑛,ℓ,𝑚}𝐺_{𝑟,𝑠,𝑡}, (2.14) where the complex coefficients, 𝐻^𝑟𝑠𝑡_{𝑛,ℓ,𝑚}, and the 3D geometric moments are defined as 𝐻^𝑟𝑠𝑡_{𝑛,ℓ,𝑚}=𝐶_ℓ,𝑚2^−𝑚⋅^𝑘_𝑣=0𝑄_{𝑘,ℓ,𝑣}⋅^𝑣_𝛼=0𝑣𝛼⋅^𝑣−𝛼_𝛽=0𝛽⋅𝑣−𝛼^𝑚_𝑢=0(−1)^𝑚−𝑢𝑚𝑢𝑖̂^𝑢⋅^{⌊(ℓ−𝑚)/2⌋}_𝜇=0(−1)^𝜇2^−2𝜇ℓ𝜇⋅ℓ−𝜇𝑚+𝜇^𝜇𝜂=0𝜇𝜂,𝐺_{𝑟,𝑠,𝑡}=_|𝑋|≤1𝑓(𝑋)𝑥^𝑟𝑦^𝑠𝑥^𝑡𝑑𝑋, (2.15)(2.16) with 𝑟=2(𝜂+𝛼)+𝑢, 𝑠=2(𝜇−𝜂+𝛽)+𝑚−𝑢, and 𝑡=2(𝜈−𝛼−𝛽−𝜇)+ℓ−𝑚.

The rotation invariance of 3D Zernike moments could be achieved where the moments are collected into (2ℓ+1) vectors Ω_𝑛,ℓ=(Ω_𝑛,ℓ,ℓ,Ω_{𝑛,ℓ,ℓ−1},Ω_{𝑛,ℓ,ℓ−2},…,Ω_{𝑛,ℓ,−ℓ})^𝑇. The rotational invariant 3D Zernike descriptors, 𝐹_𝑛,ℓ, are defined as the norms of these vectors. To avoid any misunderstanding, an example is illustrated. For maximum moment order, Max=10, the total number of 3D Zernike moments is equal to 286. The number of independent 3D Zernike moments is 161, while the number of 3D Zernike descriptors for this moment order is 36. The number of 3D Zernike descriptors could be easily determined using the following form: Total=Max+22²(),MaxisevenMax+1)(Max+34,Maxisodd. (2.17)

The independent 3D Zernike moments are used in 3D image reconstruction, while the 3D Zernike descriptors are used in comparing similar structure by simply comparing the vectors of these descriptors. Theoretically, if a 3D object is rotated with any angle, the vectors of 3D Zernike descriptors must be the same for both original and rotated objects.

3. The Proposed Algorithm

3D Zernike moments are defined as the projection of the digital image/object onto the 3D Zernike polynomials. These polynomials are defined in spherical coordinates within a unit ball. By converting to Cartesian coordinates, the 3D image/object is defined within a cube; this cube is completely surrounded by this unit ball where its centre coincides with the centre of the unit ball. The coordinate axes 𝑥, 𝑦, and 𝑧 divide the mentioned cube into eight equal small cubes.

The 3D digital image/object of size 𝑁×𝑁×𝑁 is a multidimensional array of voxels; centers of these voxels are the points (𝑥_𝑖,𝑦_𝑗,𝑧_𝑘) where the intensity function is defined only for this discrete set of points  (𝑥_𝑖,𝑦_𝑗,𝑧_𝑘)∈⌊−1/√3,1/√3⌋×⌊−1/√3,1/√3⌋×⌊−1/√3,1/√3⌋ with: 𝑥_𝑖=2𝑖−𝑁−1𝑁√3,𝑦_𝑗=2𝑗−𝑁−1𝑁√3,𝑧_𝑘=2𝑘−𝑁−1𝑁√3, (3.1) where the sampling intervals in the 𝑥-, 𝑦-, and 𝑧-directions are Δ𝑥_𝑖=𝑥_𝑖+1−𝑥_𝑖, Δ𝑦_𝑗=𝑦_𝑗+1−𝑦_𝑗, and Δ𝑧_𝑘=𝑧_𝑘+1−𝑧_𝑘 with 𝑖,𝑗,𝑘=1,2,…,𝑁.

Computation of 3D Zernike moments in Cartesian coordinated is completely dependent on two computational modules. The first one consists of the computational processes required to compute the complex coefficients 𝐻^𝑟𝑠𝑡_{𝑛,ℓ,𝑚}. The second is the process of computing the 3D geometric moments. In order to reduce the overall computational complexity, both modules must be designed and executed by using efficient methodology.

The extremely time-consuming computational process of the first module could be significantly reduced through the recurrence relations and avoiding the repeated evaluation of factorial terms. The computational complexity of the second module is significantly reduced through the implementation of a symmetry property and the successive computation of 1D cascade for each moment order. Through the next subsection, a detailed description of the proposed efficient method is presented.

3.1. Computational Aspects

Fast algorithms are generally desired in the computational processes of 3D moments. In this subsection, the attention is paid to reduce the computational complexity of the first module. Analysis of (2.15) shows that the computation of the complex coefficients 𝐻^𝑟𝑠𝑡_{𝑛,ℓ,𝑚} for each moment order required the computation of the normalization factors 𝐶_ℓ,𝑚, the coefficients 𝑄_{𝑘,ℓ,𝑣}, and a huge number of factorial terms. Computation of each of these coefficients is a source of excessive computational complexity. These time-consuming processes are repeated with each moment order. In order to reduce the overall computational complexity, an efficient method must be applied to overcome these aforementioned sources of computational complexity.

A recurrence relation is derived and employed to compute the normalization factors 𝐶_ℓ,𝑚. The derived relation completely ignored the heavy computational costs of factorial terms. It is clear that these normalization factors are image-independent; so, their values could be precalculated, stored, and recalled whenever they are needed. This is another attractive point. According to (2.8), the derived recurrence relations are 𝐶_0,0𝐶=1,_ℓ,0=√𝐶2ℓ+1,_ℓ,𝑚=ℓ+𝑚𝐶ℓ−𝑚+1^ℓ,𝑚−1, (3.2)(3.3)(3.4) where ℓ=1,2,3,…Max and 𝑚=1,2,3,…ℓ.

The second source of excessive computational complexity is the coefficients 𝑄_{𝑘,ℓ,𝑣} defined by (2.6). To overcome this problem, (2.6) will be rewritten as follows: 𝑄_{𝑘,ℓ,𝑣}=(−1)^𝑘+𝑣2^2𝑘2ℓ+4𝑘+33𝑇_{𝑘,ℓ,𝑣}, (3.5) with 𝑇_{𝑘,ℓ,𝑣}=^{𝑘2𝑘𝑘𝑣2()𝑘+ℓ+𝑣+12𝑘}_{𝑘𝑘+ℓ+𝑣}. (3.6) The time-consuming direct computations of combinational terms are avoided by using the recurrence relations [29] where an image-independent matrix could be precomputed, stored, and recalled whenever it is needed. The combinational terms of (3.6) could be easily computed using the stored values of this matrix. For moment order 𝑛∈[0,Max], ℓ∈[0,𝑛], 𝑘 is an integer equal to 𝑘=(𝑛−ℓ)/2 and 0≤𝑣≤𝑘, the numerical values of the coefficients 𝑇_{𝑘,ℓ,𝑣} will be precalculated, stored, and recalled any time needed. Consequently, the coefficients 𝑄_{𝑘,ℓ,𝑣} could be efficiently computed without any combination or factorial terms. So, we could overcome the second source of excessive computational complexity. For the special case, ℓ=𝑛, the coefficients 𝑄_{𝑘,ℓ,𝑣} are easily computed using the following equation: 𝑄_0,ℓ,0=2ℓ+33. (3.7) To see the importance of (3.7), the 3D image reconstruction using 3D Zernike moments of maximum order 10 required the computation of 161 independent moments. Equation (3.7) computes 66 of these moments where the highest number of moments is found when the order 𝑛 equals to ℓ.

The third source of computational complexity is the direct computation of the complex coefficients 𝐻^𝑟𝑠𝑡_{𝑛,ℓ,𝑚} using (2.15). Factorial terms represent the main source of the complexity in (2.15) where ignoring these factorial terms results in a great simplicity. Replacing factorial terms in (2.15) by their numerical values will achieve this goal.

3.2. Symmetry Property

According to the mapping of the 3D object inside the unit ball, the radial distance from any point 𝑃₁(𝑥_𝑖,𝑦_𝑗,𝑧_𝑘) to the coordinate origin is 𝑥_𝑖²+𝑦_𝑗²+𝑧_𝑘². Based on the definition of radial distance, there are eight points in the different eight small cubes having the same radial distance to the coordinate origin. These points and their coordinates are shown in Table 1.

Table 1 

3D symmetry points and their coordinates.

Since the points {𝑃_𝑑,𝑑=1,2,3,…,8} have the same radial distance, then the numerical values of 𝑥^𝑟𝑦^𝑠𝑧^𝑡 will be dependent on whatever 𝑟, 𝑠, and 𝑡 are even or odd. For quick proof of this assumption, we consider the following illustrative example. Let the first point 𝑃₁ has the coordinates 𝑃₁(𝑥₄,𝑦₃,𝑧₂)≡𝑃₁(−1/8√3,−3/8√3,−5/8√3). Based on coordinate relations defined in Table 1, the other seven points could be written as 𝑃₂(1/8√3,−3/8√3,−5/8√3), 𝑃₃(1/8√3,3/8√3,−5/8√3), 𝑃₄(−1/8√3,3/8√3,−5/8√3), 𝑃₅(−1/8√3,−3/8√3,5/8√3), 𝑃₆(1/8√3,−3/8√3,5/8√3), 𝑃₇(1/8√3,3/8√3,5/8√3), and 𝑃₈(−1/8√3,3/8√3,5/8√3). Numerical values of 𝑥^𝑟𝑦^𝑠𝑧^𝑡 for the points {𝑃_𝑑,𝑑=1,2,3,…,8} with different possibilities of exponent indices, 𝑟, 𝑠, and 𝑡 are listed in Table 2.

Table 2 

Numerical values of 𝑥^𝑟𝑦^𝑠𝑧^𝑡 are dependent on whatever 𝑟, 𝑠, and 𝑡 are even or odd.

Based on this symmetry property and the numerical results in Table 2, we can define the different eight cases for what is called augmented intensity function as follows.

Case 1. 𝑟=𝐸, 𝑠=𝐸, and 𝑡=𝐸; 𝑓_𝐴𝑥_𝑖,𝑦_𝑗,𝑧_𝑘=𝑓₁𝑥_𝑖,𝑦_𝑗,𝑧_𝑘+𝑓₂𝑥_𝑖,𝑦_𝑗,𝑧_𝑘+𝑓₃𝑥_𝑖,𝑦_𝑗,𝑧_𝑘+𝑓₄𝑥_𝑖,𝑦_𝑗,𝑧_𝑘+𝑓₅𝑥_𝑖,𝑦_𝑗,𝑧_𝑘+𝑓₆𝑥_𝑖,𝑦_𝑗,𝑧_𝑘+𝑓₇𝑥_𝑖,𝑦_𝑗,𝑧_𝑘+𝑓₈𝑥_𝑖,𝑦_𝑗,𝑧_𝑘, (3.8a)

Case 2. 𝑟=𝐸, 𝑠=𝐸, and 𝑡=𝑂; 𝑓_𝐴𝑥_𝑖,𝑦_𝑗,𝑧_𝑘=−𝑓₁𝑥_𝑖,𝑦_𝑗,𝑧_𝑘−𝑓₂𝑥_𝑖,𝑦_𝑗,𝑧_𝑘−𝑓₃𝑥_𝑖,𝑦_𝑗,𝑧_𝑘−𝑓₄𝑥_𝑖,𝑦_𝑗,𝑧_𝑘+𝑓₅𝑥_𝑖,𝑦_𝑗,𝑧_𝑘+𝑓₆𝑥_𝑖,𝑦_𝑗,𝑧_𝑘+𝑓₇𝑥_𝑖,𝑦_𝑗,𝑧_𝑘+𝑓₈𝑥_𝑖,𝑦_𝑗,𝑧_𝑘, (3.8b)

Case 3. 𝑟=𝐸, 𝑠=𝑂, and 𝑡=𝐸; 𝑓_𝐴𝑥_𝑖,𝑦_𝑗,𝑧_𝑘=−𝑓₁𝑥_𝑖,𝑦_𝑗,𝑧_𝑘−𝑓₂𝑥_𝑖,𝑦_𝑗,𝑧_𝑘+𝑓₃𝑥_𝑖,𝑦_𝑗,𝑧_𝑘+𝑓₄𝑥_𝑖,𝑦_𝑗,𝑧_𝑘−𝑓₅𝑥_𝑖,𝑦_𝑗,𝑧_𝑘−𝑓₆𝑥_𝑖,𝑦_𝑗,𝑧_𝑘+𝑓₇𝑥_𝑖,𝑦_𝑗,𝑧_𝑘+𝑓₈𝑥_𝑖,𝑦_𝑗,𝑧_𝑘, (3.8c)

Case 4. 𝑟=𝐸, 𝑠=𝑂, and 𝑡=𝑂; 𝑓_𝐴𝑥_𝑖,𝑦_𝑗,𝑧_𝑘=𝑓₁𝑥_𝑖,𝑦_𝑗,𝑧_𝑘+𝑓₂𝑥_𝑖,𝑦_𝑗,𝑧_𝑘−𝑓₃𝑥_𝑖,𝑦_𝑗,𝑧_𝑘−𝑓₄𝑥_𝑖,𝑦_𝑗,𝑧_𝑘−𝑓₅𝑥_𝑖,𝑦_𝑗,𝑧_𝑘−𝑓₆𝑥_𝑖,𝑦_𝑗,𝑧_𝑘+𝑓₇𝑥_𝑖,𝑦_𝑗,𝑧_𝑘+𝑓₈𝑥_𝑖,𝑦_𝑗,𝑧_𝑘, (3.8d)

Case 5. 𝑟=𝑂, 𝑠=𝐸, and 𝑡=𝐸; 𝑓_𝐴𝑥_𝑖,𝑦_𝑗,𝑧_𝑘=−𝑓₁𝑥_𝑖,𝑦_𝑗,𝑧_𝑘+𝑓₂𝑥_𝑖,𝑦_𝑗,𝑧_𝑘+𝑓₃𝑥_𝑖,𝑦_𝑗,𝑧_𝑘−𝑓₄𝑥_𝑖,𝑦_𝑗,𝑧_𝑘−𝑓₅𝑥_𝑖,𝑦_𝑗,𝑧_𝑘+𝑓₆𝑥_𝑖,𝑦_𝑗,𝑧_𝑘+𝑓₇𝑥_𝑖,𝑦_𝑗,𝑧_𝑘−𝑓₈𝑥_𝑖,𝑦_𝑗,𝑧_𝑘, (3.8e)

Case 6. 𝑟=𝑂, 𝑠=𝐸, and 𝑡=𝑂; 𝑓_𝐴𝑥_𝑖,𝑦_𝑗,𝑧_𝑘=𝑓₁𝑥_𝑖,𝑦_𝑗,𝑧_𝑘−𝑓₂𝑥_𝑖,𝑦_𝑗,𝑧_𝑘−𝑓₃𝑥_𝑖,𝑦_𝑗,𝑧_𝑘+𝑓₄𝑥_𝑖,𝑦_𝑗,𝑧_𝑘−𝑓₅𝑥_𝑖,𝑦_𝑗,𝑧_𝑘+𝑓₆𝑥_𝑖,𝑦_𝑗,𝑧_𝑘+𝑓₇𝑥_𝑖,𝑦_𝑗,𝑧_𝑘−𝑓₈𝑥_𝑖,𝑦_𝑗,𝑧_𝑘, (3.8f)

Case 7. 𝑟=𝑂, 𝑠=𝑂, and 𝑡=𝐸; 𝑓_𝐴𝑥_𝑖,𝑦_𝑗,𝑧_𝑘=𝑓₁𝑥_𝑖,𝑦_𝑗,𝑧_𝑘−𝑓₂𝑥_𝑖,𝑦_𝑗,𝑧_𝑘+𝑓₃𝑥_𝑖,𝑦_𝑗,𝑧_𝑘−𝑓₄𝑥_𝑖,𝑦_𝑗,𝑧_𝑘+𝑓₅𝑥_𝑖,𝑦_𝑗,𝑧_𝑘−𝑓₆𝑥_𝑖,𝑦_𝑗,𝑧_𝑘+𝑓₇𝑥_𝑖,𝑦_𝑗,𝑧_𝑘−𝑓₈𝑥_𝑖,𝑦_𝑗,𝑧_𝑘, (3.8g)

Case 8. 𝑟=𝑂, 𝑠=𝑂, and 𝑡=𝑂; 𝑓_𝐴𝑥_𝑖,𝑦_𝑗,𝑧_𝑘=−𝑓₁𝑥_𝑖,𝑦_𝑗,𝑧_𝑘+𝑓₂𝑥_𝑖,𝑦_𝑗,𝑧_𝑘−𝑓₃𝑥_𝑖,𝑦_𝑗,𝑧_𝑘+𝑓₄𝑥_𝑖,𝑦_𝑗,𝑧_𝑘+𝑓₅𝑥_𝑖,𝑦_𝑗,𝑧_𝑘−𝑓₆𝑥_𝑖,𝑦_𝑗,𝑧_𝑘+𝑓₇𝑥_𝑖,𝑦_𝑗,𝑧_𝑘−𝑓₈𝑥_𝑖,𝑦_𝑗,𝑧_𝑘, (3.8h) where 𝑓₁,𝑓₂,𝑓₃,𝑓₄,𝑓₅,𝑓₆,𝑓₇, and 𝑓₈ refer to the image intensity function defined in the cubes from 1 to 8, respectively. The letters “𝐸” and “𝑂” are the acronyms of even and odd, respectively.

As shown in Figure 1, the removed small cube refers to the first subcube. Therefore, only one-eighth of the whole object space is required to compute the entire set of 3D Zernike moments. The implementation of this property results in 87% reduction in the computational cost. A detailed discussion of this will be found through the following subsections.

Figure 1 

Input 3D image is mapped inside the unit ball.

3.3. Fast Computation of Exact 3D Geometric Moments

General 3D geometric moments of order (𝑟+𝑠+𝑡) for the object intensity function 𝑓(𝑥,𝑦,𝑧) are defined as the projection of the function 𝑓(𝑥,𝑦,𝑧) onto the monomial 𝑥^𝑟𝑦^𝑠𝑧^𝑡 as follows: 𝐺_𝑟𝑠𝑡=^∞_−∞𝑥^𝑟𝑦^𝑠𝑧^𝑡𝑓𝑥,𝑦,𝑧𝑑𝑥𝑑𝑦𝑑𝑧. (3.9) Based on the circumstances of the present problem, where the input 3D image is mapped inside the unit ball, the upper and lower limits of triple integrals in (3.9) must be rewritten as follows: 𝐺_𝑟𝑠𝑡=^√1/3_√−1/3𝑥^𝑟𝑦^𝑠𝑧^𝑡𝑓𝑥,𝑦,𝑧𝑑𝑥𝑑𝑦𝑑𝑧. (3.10) Equation (3.10) could be written as follows: 𝐺_𝑝𝑞𝑟=^𝑁𝑁_𝑖=1^𝑁_𝑗=1𝑘=1𝑆_𝑝𝑞𝑟𝑥_𝑖,𝑦_𝑗,𝑧_𝑘𝑓𝑥_𝑖,𝑦_𝑗,𝑧_𝑘, (3.11) where 𝑆_𝑝𝑞𝑟𝑥_𝑖,𝑦_𝑗,𝑧_𝑘=^{𝑥𝑖+(Δ𝑥𝑖}_𝑥^/2)_{𝑖−Δ𝑥𝑖/2}^{𝑦𝑗+(Δ𝑦𝑗}_𝑦^/2)_{𝑗−Δ𝑦𝑗/2}^{𝑧𝑘+(Δ𝑧𝑘}_𝑧^/2)_{𝑘−Δ𝑧𝑘/2}𝑥^𝑟𝑦^𝑠𝑧^𝑡𝑑𝑥𝑑𝑦𝑑𝑧. (3.12) The triple integral defined by (3.12) is the source of approximation error. For exact computation of 3D geometric moments, this triple integral could be divided into three separate individual integrals as follows:𝐼_𝑟(𝑖)=^{𝑥𝑖+(Δ𝑥𝑖}_𝑥^/2)_{𝑖−Δ𝑥𝑖/2}𝑥^𝑟1𝑑𝑥=𝑥𝑟+1^𝑖+Δ𝑥_𝑖2^𝑟+1−𝑥_𝑖−Δ𝑥_𝑖2^𝑟+1𝐼_𝑠𝑗=^{𝑦𝑗+(Δ𝑦𝑗}_𝑦^/2)_{𝑗−Δ𝑦𝑗/2}𝑦^𝑠1𝑑𝑦=𝑦𝑠+1^𝑗+Δ𝑦_𝑗2^𝑠+1−𝑦_𝑗−Δ𝑦_𝑗2^𝑠+1𝐼_𝑡(𝑘)=^{𝑧𝑡+(Δ𝑧𝑡}_𝑧^/2)_{𝑡−Δ𝑧𝑡/2}𝑧^𝑡1𝑑𝑧=𝑧𝑡+1^𝑡+Δ𝑧_𝑡2^𝑡+1−𝑧_𝑡−Δ𝑧_𝑡2^𝑡+1. (3.13a)(3.13b)(3.13c)The upper and lower limits of these integrals are created and stored in a vector form. By applying the symmetry property, the set of geometric moments can thus be computed exactly by 𝐺_𝑝𝑞𝑟=^⌊𝑁/2⌋_𝑖=1^⌊𝑁/2⌋_𝑗=1^⌊𝑁/2⌋_𝑘=1𝐼_𝑟(𝑖)𝐼_𝑠𝑗𝐼_𝑡(𝑘)𝑓_𝐴𝑥_𝑖,𝑦_𝑗,𝑧_𝑘, (3.14) where 𝑓_𝐴(𝑥_𝑖,𝑦_𝑗,𝑧_𝑘) is the augmented intensity function defined by using (3.8a)–(3.8h). For even value of 𝑁, the operator ⌊𝑁/2⌋ equals to 𝑁/2, while it equals to (𝑁−1)/2 for odd value. The kernel of the exact 3D geometric moments is defined by (3.13a)–(3.13c). This kernel is image-independent. Therefore, this kernel could be precomputed, stored, and recalled whenever it is needed to avoid repetitive computation.

The computational complexity of exact 3D geometric moments could be greatly reduced by applying the successive computation process. The 3D geometric moments of order (𝑟+𝑠+𝑡) defined by (3.14) are computed in three separate steps by successive computation of the 1D 𝑠th-order moment for each row, followed by the 2D (𝑟+𝑠)th-order moment. Then, the required 3D moment is calculated as the sum of the different 2D moments. This approach was successfully applied in recent works of Hosny [31, 32] where it significantly reduces the total number of required addition and multiplication processes. Equation (3.14) is rewritten in separable forms as follows:𝐺_𝑟𝑠𝑡=^⌊𝑁/2⌋_𝑘=1𝐼_𝑡𝑧_𝑘𝑅_𝑟𝑠𝑘, (3.15a)where𝑅_𝑟𝑠𝑘=^⌊𝑁/2⌋_𝑖=1𝑌_𝑖𝑠𝑘𝐼_𝑟𝑥_𝑖,𝑌_𝑖𝑠𝑘=^⌊𝑁/2⌋_𝑗=1𝐼_𝑠𝑦_𝑗𝑓_𝐴𝑥_𝑖,𝑦_𝑗,𝑧_𝑘. (3.16a)(3.16b)

4. Experimental Results

The computation of 3D Zernike moments as a liner combination of exact 3D geometric moments ensures the accuracy. This approach was proved through our previously published works of 2D Zernike moments [29, 30]. Based on this fact, this work concentrates on the issue of efficiency in computational time and complexity. The conducted numerical experiments concentrate on the efficiency of the proposed method against the existing methods. The full set of 3D Zernike moments is computed by using the proposed method and the conventional method [21]. Elapsed CPU times for both methods are used to judge the efficiency. Two numerical experiments are conducted. In the first experiment, a randomly generated 3D image with intensity function 𝑓(𝑥_𝑖,𝑦_𝑗,𝑧_𝑘) is generated by using the Matlab8 statement 𝑓(𝑥_𝑖,𝑦_𝑗,𝑧_𝑘)=𝑟  and  (𝑀,𝑁,𝐾) where 0≤𝑓(𝑥_𝑖,𝑦_𝑗,𝑧_𝑘)≤1 for all 𝑖,𝑗,𝑘 with 𝑀=𝑁=𝐾=64. Selected orders of 3D Zernike moments and the corresponding elapsed CPU times are for both methods shown in Table 3. A graphical representation of these elapsed CPU times is plotted in Figure 2.

Table 3 

Elapsed CPU times and the reduction percentage for selected moment orders (first numerical experiment).

Figure 2 

Elapsed CPU times for computing 3D Zernike moments (first numerical experiment).

In the second experiment, a protein structure of dimensions 70×70×70 is used. All computational processes are performed by using a code designed with Matlab8 and operated on a Lenovo R400 laptop. Similar to the first numerical experiment, selected orders of 3D Zernike moments and their elapsed CPU times for both methods are shown in Table 4. The graphical representation of the elapsed CPU times for moment orders ranging from 1 to 15 is plotted in Figure 3. The logarithmic scale is more suitable to clearly show the big differences in the elapsed CPU times.

Table 4 

Elapsed CPU times and the reduction percentage for selected moment orders (second numerical experiment).

Figure 3 

Elapsed CPU times for computing 3D Zernike moments (second numerical experiment).

It is clear that the proposed method outperformed the conventional method where the reduction in the elapsed CPU times exceeds 95%. In addition to this big reduction, the implementation of symmetry property achieved 87% of memory saving. On the other side, the conventional method does not have any kind of memory saving. Generally, the comparison ensures the superiority of the proposed method.

4.1. Computational Complexity

Complexity analysis of any numerical method presents a simple and clear way to judge the efficiency of this method. Complexity analysis mainly concentrates on the number of operations required by such a method to achieve its goal. The total number of multiplication and addition operations is the core of the complexity analysis. Evaluation of additional operations such as factorial terms, exponential and power functions is also considered.

The complexity analysis of the proposed and the conventional methods for computing 3D Zernike moments is performed. The computational process of 3D Zernike moments required the computation of 3D geometric moments and the complex coefficients 𝐻^𝑟𝑠𝑡_{𝑛,ℓ,𝑚}. It is clear that the computation of 3D geometric moments is image-dependent, while the computation of 𝐻^𝑟𝑠𝑡_{𝑛,ℓ,𝑚} is image-independent. Therefore, we concentrate on the complexity analysis of 3D geometric moments.

For a 3D digital image of size 𝑁×𝑁×𝑁 and a maximum moment order equal to Max, the complexity analysis of the conventional method represented by (2.16) is discussed first. The total number of arithmetic operations (addition and multiplication) required by the conventional method for computing 3D geometric moments is(Max+1)(Max+2)(Max+3)𝑁³6()𝑁,additionsMax+1)(Max+2)(Max+3³3,multiplications. (4.1a)(4.1b)The corresponding total number of arithmetic operations required by the proposed method will be evaluated. The computational process of the 3D geometric moments by using the proposed method consists of three steps represented by (3.15a), (3.16a), and (3.16b). The total number of arithmetic operations for each step will be discussed individually; then the whole computational complexity will be evaluated. Starting with (3.16b), the creation of the matrix 𝑌_𝑖𝑞𝑘 requires (𝑁/2)³(Max+1) multiplications and (𝑁/2)²(𝑁/2−1)(Max+1) additions. The creation process of the matrix 𝑅_𝑝𝑞𝑘 using (3.16a) requires (𝑁/2)((𝑁/2)−1)(Max+1)(Max+2)/2 additions and (𝑁/2)²(Max+1)(Max+2)/2 multiplications. Finally, the computational complexity of the 3D geometric moments using (3.15a) requires (𝑁/2)(Max+1)²(Max+2)/2 multiplications and ((𝑁/2)−1)(Max+1)²(Max+2)/2 additions. Therefore, computing the set of independent exact 3D geometric moments requires the following number of additions and multiplications:(Max+1)(Max+2)8𝑁²(+𝑁−2𝑁+2Max𝑁−2)²(𝑁−2)(Max+1)8(),additionsMax+1)(Max+28𝑁²(+𝑁+2𝑁Max+1)³(Max+1)8,multiplications. (4.2a)(4.2b)A quick and clear comparison of the complexity of both methods is shown in Table 5 where the total number of arithmetic operations is required by both methods for computing 3D geometric moments with the dimension 𝑁 and maximum moment order Max. It is clear that the proposed method tremendously reduced the total number of arithmetic operations.

Table 5 

Number of addition and multiplication processes required by both methods.

5. Conclusion

This paper proposes very fast and computationally efficient method for computing 3D Zernike moments. In the proposed method, 3D Zernike moments are expressed as a linear combination of 3D geometric moments. Numerical experiments show the efficiency of the proposed method, where it achieves more than 95% saving in elapsed CPU times. The implementation of the symmetry property achieves 87% memory saving, which is a very attractive property especially in the processing of 3D images.

Acknowledgments

k. M. Hosny and M. A. Hafez are gratefully acknowledging the financial support of Najran University (Grant no. NU43/2010).