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

Authors: Khalid M. Hosny; Mohamed A. Hafez
Journal: Mathematical Problems in Engineering
 (2012)

---

## 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.


---

## Body

## 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,Zn,ℓ,m, are orthogonal polynomials defined on a unit ball as follows:
(2.1)Zn,ℓ,m(ℜ)=Rn,ℓ(r)Yℓ,m(θ,ϕ),
where n∈[0,Max], ℓ∈[0,n], and m∈[-ℓ,ℓ]. The value (n-ℓ) must be even nonnegative integer number. Max is the maximum order considered in the computational process. Rn,ℓ(r) and Yℓ,m(θ,ϕ) are real-valued radial functions and spherical Harmonics, respectively. Any function, f(ℜ), defined within the unit ball could be expanded by using the 3D Zernike polynomials as follows:
(2.2)f(ℜ)=∑n=0∞∑ℓ=0n∑m=-ℓℓΩn,ℓ,mZn,ℓ,m(ℜ).
The expansion coefficients, Ωn,ℓ,m, represent the 3D Zernike moments and determined by using the complex conjugate of 3D Zernike polynomials as follows:
(2.3)Ωn,ℓ,m=∫01∫02π∫0πZn,ℓ,m(ℜ)¯f(ℜ)r2sinθdrdθdϕ.
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, [X]T=[ℜ]T, is defined in an explicit form as follows:
(2.4)[xyz]=[rsinθsinϕrsinθcosϕrcosϕ].
The 3D Zernike polynomials in the Cartesian coordinates are defined as follows:
(2.5)Zn,ℓ,m(X)=∑v=0kQk,ℓ,v|X|2veℓ,m(X),
where k is an integer equal to (n-ℓ)/2 and 0≤v≤k. The coefficients Qk,ℓ,v are defined as
(2.6)Qk,ℓ,v=(-1)k22k2ℓ+4k+33(2kk)(-1)v(kv)(2(k+ℓ+v)+12k)(k+ℓ+vk).
The harmonic polynomials, eℓ,m(X), are defined as
(2.7)eℓ,m(X)=Cℓ,mrℓ(i^x-y2)mzℓ-m∑μ=0⌊(ℓ-m)/2⌋(ℓμ)(ℓ-μm+μ)(-x2+y24z2)μ,
where i^=-1, z=x+i^y is the complex variable and Cℓ,m are the normalization factors and defined as
(2.8)Cℓ,m=(2ℓ+1)(ℓ+m)!(ℓ-m)!ℓ!.
Based on (2.8), the normalization factors Cℓ,m of positive and negative values of m are identical, Cℓ,-m=Cℓ,m while the harmonic polynomials with negative values of m are defined as
(2.9)eℓ,-m(X)=(-1)meℓ,m(X)¯.
The orthogonality relation of 3D Zernike polynomials is defined as follows:
(2.10)34π∫|X|≤1Zn,ℓ,m(X)Zn′,ℓ′,m′(X)¯dX=δn,n′δℓ,ℓ′δm,m′.
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:
(2.11)Ωn,ℓ,m=34π∫|X|≤1f(X)Zn,ℓ,m(X)¯dXf(X)=∑n=0∞∑ℓ=0n∑m=-ℓℓΩn,ℓ,mZn,ℓ,m(X),
where the 3D Zernike moments with negative values of m could be computed directly from their corresponding ones with positive values of m as follows:
(2.12)Ωn,ℓ,-m(X)=(-1)mΩn,ℓ,m(X)¯.
Canterkis [20] compactly formulate the 3D Zernike polynomials of order n as a linear combination of monomials of order up to n as follows:
(2.13)Zn,ℓ,m(X)=∑r+s+t≤nHn,ℓ,mrstxrysxt.
Consequently, the 3D Zernike moments Ωn,ℓ,m of order n could be expressed as a linear combination of the 3D geometric moments, Gr,s,t, by using the following relation:
(2.14)Ωn,ℓ,m=34π∑r+s+t≤nHn,ℓ,mrst¯Gr,s,t,
where the complex coefficients, Hn,ℓ,mrst, and the 3D geometric moments are defined as
(2.15)Hn,ℓ,mrst=Cℓ,m2-m·∑v=0kQk,ℓ,v·∑α=0v(vα)·∑β=0v-α(v-αβ)·∑u=0m(-1)m-u(mu)(i^)u·∑μ=0⌊(ℓ-m)/2⌋(-1)μ2-2μ(ℓμ)(ℓ-μm+μ)·∑η=0μ(μη),(2.16)Gr,s,t=∫|X|≤1f(X)xrysxtdX,
with r=2(η+α)+u, s=2(μ-η+β)+m-u, and t=2(ν-α-β-μ)+ℓ-m.The rotation invariance of 3D Zernike moments could be achieved where the moments are collected into(2ℓ+1) vectors Ωn,ℓ=(Ωn,ℓ,ℓ,Ωn,ℓ,ℓ-1,Ωn,ℓ,ℓ-2,…,Ωn,ℓ,-ℓ)T. The rotational invariant 3D Zernike descriptors, Fn,ℓ, 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:
(2.17)Total={(Max+22)2,Max  is  even(Max+1)(Max+3)4,Max  is  odd.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 axesx, y, and z divide the mentioned cube into eight equal small cubes.The 3D digital image/object of sizeN×N×N is a multidimensional array of voxels; centers of these voxels are the points (xi,yj,zk) where the intensity function is defined only for this discrete set of points  (xi,yj,zk)∈⌊-1/3,1/3⌋×⌊-1/3,1/3⌋×⌊-1/3,1/3⌋ with:
(3.1)xi=2i-N-1N3,yj=2j-N-1N3,zk=2k-N-1N3,
where the sampling intervals in the x-, y-, and z-directions are Δxi=xi+1-xi, Δyj=yj+1-yj, and Δzk=zk+1-zk with i,j,k=1,2,…,N.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 coefficientsHn,ℓ,mrst. 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 Hn,ℓ,mrst for each moment order required the computation of the normalization factors Cℓ,m, the coefficients Qk,ℓ,v, 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 factorsCℓ,m. 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
(3.2)C0,0=1,(3.3)Cℓ,0=2ℓ+1,(3.4)Cℓ,m=ℓ+mℓ-m+1Cℓ,m-1,
where ℓ=1,2,3,…Max and m=1,2,3,…ℓ.The second source of excessive computational complexity is the coefficientsQk,ℓ,v defined by (2.6). To overcome this problem, (2.6) will be rewritten as follows:
(3.5)Qk,ℓ,v=(-1)k+v22k2ℓ+4k+33Tk,ℓ,v,
with
(3.6)Tk,ℓ,v=(2kk)(kv)(2(k+ℓ+v)+12k)(k+ℓ+vk).
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 n∈[0,Max], ℓ∈[0,n], k is an integer equal to k=(n-ℓ)/2 and 0≤v≤k, the numerical values of the coefficients Tk,ℓ,v will be precalculated, stored, and recalled any time needed. Consequently, the coefficients Qk,ℓ,v 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, ℓ=n, the coefficients Qk,ℓ,v are easily computed using the following equation:
(3.7)Q0,ℓ,0=2ℓ+33.
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 n equals to ℓ.The third source of computational complexity is the direct computation of the complex coefficientsHn,ℓ,mrst 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 pointP1(xi,yj,zk) to the coordinate origin is xi2+yj2+zk2. 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.
Points
Coordinates
P
1
(
x
i
,
y
j
,
z
k
)
P
2
(
x
N
-
i
+
1
,
y
j
,
z
k
)
P
3
(
x
N
-
i
+
1
,
y
N
-
j
+
1
,
z
k
)
P
4
(
x
i
,
y
N
-
j
+
1
,
z
k
)
P
5
(
x
i
,
y
j
,
z
N
-
k
+
1
)
P
6
(
x
N
-
i
+
1
,
y
j
,
z
N
-
k
+
1
)
P
7
(
x
N
-
i
+
1
,
y
N
-
j
+
1
,
z
N
-
k
+
1
)
P
8
(
x
i
,
y
N
-
j
+
1
,
z
N
-
k
+
1
)Since the points{Pd,d=1,2,3,…,8} have the same radial distance, then the numerical values of xryszt will be dependent on whatever r, s, and t are even or odd. For quick proof of this assumption, we consider the following illustrative example. Let the first point P1 has the coordinates P1(x4,y3,z2)≡P1(-1/83,-3/83,-5/83). Based on coordinate relations defined in Table 1, the other seven points could be written as P2(1/83,-3/83,-5/83), P3(1/83,3/83,-5/83), P4(-1/83,3/83,-5/83), P5(-1/83,-3/83,5/83), P6(1/83,-3/83,5/83), P7(1/83,3/83,5/83), and P8(-1/83,3/83,5/83). Numerical values of xryszt for the points {Pd,d=1,2,3,…,8} with different possibilities of exponent indices, r, s, and t are listed in Table 2.Table 2
Numerical values ofxryszt are dependent on whatever r, s, and t are even or odd.
r
s
t
x
r
y
s
z
t
P
1
P
2
P
3
P
4
P
5
P
6
P
7
P
8
E
=
2
E
=
0
E
=
2
+
0.00156
+
0.00156
+
0.00156
+
0.00156
+
0.00156
+
0.00156
+
0.00156
+
0.00156
E
=
0
E
=
2
O
=
1
-
0.03107
−0.03107
−0.03107
−0.03107
+
0.03107
+
0.03107
+
0.03107
+
0.03107
E
=
2
O
=
1
E
=
0
-
0.002071
−0.002071
+
0.002071
+
0.002071
−0.002071
−0.002071
+
0.002071
+
0.002071
E
=
2
O
=
1
O
=
1
+
0.00092
+
0.00092
−0.00092
−0.00092
−0.00092
−0.00092
+0.00092
+
0.00092
O
=
1
E
=
2
E
=
0
-
0.00621
+
0.00621
+
0.00621
−0.00621
−0.00621
+
0.00621
+
0.00621
−0.00621
O
=
1
E
=
0
O
=
3
+
0.007629
−0.007629
−0.007629
+
0.007629
-
0.007629
+
0.007629
+
0.007629
−0.007629
O
=
1
O
=
1
E
=
2
+
0.004577
−0.004577
+
0.004577
-
0.004577
+
0.004577
-
0.004577
+
0.004577
-
0.004577
O
=
1
O
=
1
O
=
3
−0.002023
+
0.002023
−0.002023
+
0.002023
+
0.002023
-
0.002023
+
0.002023
-
0.002023Based on this symmetry property and the numerical results in Table2, we can define the different eight cases for what is called augmented intensity function as follows.Case 1.
r
=
E, s=E, and t=E;
(3.8a)fA(xi,yj,zk)=f1(xi,yj,zk)+f2(xi,yj,zk)+f3(xi,yj,zk)+f4(xi,yj,zk)+f5(xi,yj,zk)+f6(xi,yj,zk)+f7(xi,yj,zk)+f8(xi,yj,zk),Case 2.
r
=
E, s=E, and t=O;
(3.8b)fA(xi,yj,zk)=-f1(xi,yj,zk)-f2(xi,yj,zk)-f3(xi,yj,zk)-f4(xi,yj,zk)+f5(xi,yj,zk)+f6(xi,yj,zk)+f7(xi,yj,zk)+f8(xi,yj,zk),Case 3.
r
=
E, s=O, and t=E;
(3.8c)fA(xi,yj,zk)=-f1(xi,yj,zk)-f2(xi,yj,zk)+f3(xi,yj,zk)+f4(xi,yj,zk)-f5(xi,yj,zk)-f6(xi,yj,zk)+f7(xi,yj,zk)+f8(xi,yj,zk),Case 4.
r
=
E, s=O, and t=O;
(3.8d)fA(xi,yj,zk)=f1(xi,yj,zk)+f2(xi,yj,zk)-f3(xi,yj,zk)-f4(xi,yj,zk)-f5(xi,yj,zk)-f6(xi,yj,zk)+f7(xi,yj,zk)+f8(xi,yj,zk),Case 5.
r
=
O, s=E, and t=E;
(3.8e)fA(xi,yj,zk)=-f1(xi,yj,zk)+f2(xi,yj,zk)+f3(xi,yj,zk)-f4(xi,yj,zk)-f5(xi,yj,zk)+f6(xi,yj,zk)+f7(xi,yj,zk)-f8(xi,yj,zk),Case 6.
r
=
O, s=E, and t=O;
(3.8f)fA(xi,yj,zk)=f1(xi,yj,zk)-f2(xi,yj,zk)-f3(xi,yj,zk)+f4(xi,yj,zk)-f5(xi,yj,zk)+f6(xi,yj,zk)+f7(xi,yj,zk)-f8(xi,yj,zk),Case 7.
r
=
O, s=O, and t=E;
(3.8g)fA(xi,yj,zk)=f1(xi,yj,zk)-f2(xi,yj,zk)+f3(xi,yj,zk)-f4(xi,yj,zk)+f5(xi,yj,zk)-f6(xi,yj,zk)+f7(xi,yj,zk)-f8(xi,yj,zk),Case 8.
r
=
O, s=O, and t=O;
(3.8h)fA(xi,yj,zk)=-f1(xi,yj,zk)+f2(xi,yj,zk)-f3(xi,yj,zk)+f4(xi,yj,zk)+f5(xi,yj,zk)-f6(xi,yj,zk)+f7(xi,yj,zk)-f8(xi,yj,zk),
where f1,f2,f3,f4,f5,f6,f7, and f8 refer to the image intensity function defined in the cubes from 1 to 8, respectively. The letters “E” and “O” are the acronyms of even and odd, respectively.As shown in Figure1, 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(r+s+t) for the object intensity function f(x,y,z) are defined as the projection of the function f(x,y,z) onto the monomial xryszt as follows:
(3.9)Grst=∭-∞∞xrysztf(x,y,z)dxdydz.
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:
(3.10)Grst=∭-1/31/3xrysztf(x,y,z)dxdydz.
Equation (3.10) could be written as follows:
(3.11)Gpqr=∑i=1N∑j=1N∑k=1NSpqr(xi,yj,zk)f(xi,yj,zk),
where
(3.12)Spqr(xi,yj,zk)=∫xi-(Δxi/2)xi+(Δxi/2)∫yj-(Δyj/2)yj+(Δyj/2)∫zk-(Δzk/2)zk+(Δzk/2)xrysztdxdydz.
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:(3.13a)Ir(i)=∫xi-(Δxi/2)xi+(Δxi/2)xrdx=1r+1[(xi+Δxi2)r+1-(xi-Δxi2)r+1](3.13b)Is(j)=∫yj-(Δyj/2)yj+(Δyj/2)ysdy=1s+1[(yj+Δyj2)s+1-(yj-Δyj2)s+1](3.13c)It(k)=∫zt-(Δzt/2)zt+(Δzt/2)ztdz=1t+1[(zt+Δzt2)t+1-(zt-Δzt2)t+1].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
(3.14)Gpqr=∑i=1⌊N/2⌋∑j=1⌊N/2⌋∑k=1⌊N/2⌋Ir(i)Is(j)It(k)fA(xi,yj,zk),
where fA(xi,yj,zk) is the augmented intensity function defined by using (3.8a)–(3.8h). For even value of N, the operator ⌊N/2⌋ equals to N/2, while it equals to (N-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(r+s+t) defined by (3.14) are computed in three separate steps by successive computation of the 1D sth-order moment for each row, followed by the 2D (r+s)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:(3.15a)Grst=∑k=1⌊N/2⌋It(zk)Rrsk,where(3.16a)Rrsk=∑i=1⌊N/2⌋YiskIr(xi),(3.16b)Yisk=∑j=1⌊N/2⌋Is(yj)fA(xi,yj,zk).


## 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 Hn,ℓ,mrst for each moment order required the computation of the normalization factors Cℓ,m, the coefficients Qk,ℓ,v, 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 factorsCℓ,m. 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
(3.2)C0,0=1,(3.3)Cℓ,0=2ℓ+1,(3.4)Cℓ,m=ℓ+mℓ-m+1Cℓ,m-1,
where ℓ=1,2,3,…Max and m=1,2,3,…ℓ.The second source of excessive computational complexity is the coefficientsQk,ℓ,v defined by (2.6). To overcome this problem, (2.6) will be rewritten as follows:
(3.5)Qk,ℓ,v=(-1)k+v22k2ℓ+4k+33Tk,ℓ,v,
with
(3.6)Tk,ℓ,v=(2kk)(kv)(2(k+ℓ+v)+12k)(k+ℓ+vk).
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 n∈[0,Max], ℓ∈[0,n], k is an integer equal to k=(n-ℓ)/2 and 0≤v≤k, the numerical values of the coefficients Tk,ℓ,v will be precalculated, stored, and recalled any time needed. Consequently, the coefficients Qk,ℓ,v 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, ℓ=n, the coefficients Qk,ℓ,v are easily computed using the following equation:
(3.7)Q0,ℓ,0=2ℓ+33.
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 n equals to ℓ.The third source of computational complexity is the direct computation of the complex coefficientsHn,ℓ,mrst 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 pointP1(xi,yj,zk) to the coordinate origin is xi2+yj2+zk2. 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.
Points
Coordinates
P
1
(
x
i
,
y
j
,
z
k
)
P
2
(
x
N
-
i
+
1
,
y
j
,
z
k
)
P
3
(
x
N
-
i
+
1
,
y
N
-
j
+
1
,
z
k
)
P
4
(
x
i
,
y
N
-
j
+
1
,
z
k
)
P
5
(
x
i
,
y
j
,
z
N
-
k
+
1
)
P
6
(
x
N
-
i
+
1
,
y
j
,
z
N
-
k
+
1
)
P
7
(
x
N
-
i
+
1
,
y
N
-
j
+
1
,
z
N
-
k
+
1
)
P
8
(
x
i
,
y
N
-
j
+
1
,
z
N
-
k
+
1
)Since the points{Pd,d=1,2,3,…,8} have the same radial distance, then the numerical values of xryszt will be dependent on whatever r, s, and t are even or odd. For quick proof of this assumption, we consider the following illustrative example. Let the first point P1 has the coordinates P1(x4,y3,z2)≡P1(-1/83,-3/83,-5/83). Based on coordinate relations defined in Table 1, the other seven points could be written as P2(1/83,-3/83,-5/83), P3(1/83,3/83,-5/83), P4(-1/83,3/83,-5/83), P5(-1/83,-3/83,5/83), P6(1/83,-3/83,5/83), P7(1/83,3/83,5/83), and P8(-1/83,3/83,5/83). Numerical values of xryszt for the points {Pd,d=1,2,3,…,8} with different possibilities of exponent indices, r, s, and t are listed in Table 2.Table 2
Numerical values ofxryszt are dependent on whatever r, s, and t are even or odd.
r
s
t
x
r
y
s
z
t
P
1
P
2
P
3
P
4
P
5
P
6
P
7
P
8
E
=
2
E
=
0
E
=
2
+
0.00156
+
0.00156
+
0.00156
+
0.00156
+
0.00156
+
0.00156
+
0.00156
+
0.00156
E
=
0
E
=
2
O
=
1
-
0.03107
−0.03107
−0.03107
−0.03107
+
0.03107
+
0.03107
+
0.03107
+
0.03107
E
=
2
O
=
1
E
=
0
-
0.002071
−0.002071
+
0.002071
+
0.002071
−0.002071
−0.002071
+
0.002071
+
0.002071
E
=
2
O
=
1
O
=
1
+
0.00092
+
0.00092
−0.00092
−0.00092
−0.00092
−0.00092
+0.00092
+
0.00092
O
=
1
E
=
2
E
=
0
-
0.00621
+
0.00621
+
0.00621
−0.00621
−0.00621
+
0.00621
+
0.00621
−0.00621
O
=
1
E
=
0
O
=
3
+
0.007629
−0.007629
−0.007629
+
0.007629
-
0.007629
+
0.007629
+
0.007629
−0.007629
O
=
1
O
=
1
E
=
2
+
0.004577
−0.004577
+
0.004577
-
0.004577
+
0.004577
-
0.004577
+
0.004577
-
0.004577
O
=
1
O
=
1
O
=
3
−0.002023
+
0.002023
−0.002023
+
0.002023
+
0.002023
-
0.002023
+
0.002023
-
0.002023Based on this symmetry property and the numerical results in Table2, we can define the different eight cases for what is called augmented intensity function as follows.Case 1.
r
=
E, s=E, and t=E;
(3.8a)fA(xi,yj,zk)=f1(xi,yj,zk)+f2(xi,yj,zk)+f3(xi,yj,zk)+f4(xi,yj,zk)+f5(xi,yj,zk)+f6(xi,yj,zk)+f7(xi,yj,zk)+f8(xi,yj,zk),Case 2.
r
=
E, s=E, and t=O;
(3.8b)fA(xi,yj,zk)=-f1(xi,yj,zk)-f2(xi,yj,zk)-f3(xi,yj,zk)-f4(xi,yj,zk)+f5(xi,yj,zk)+f6(xi,yj,zk)+f7(xi,yj,zk)+f8(xi,yj,zk),Case 3.
r
=
E, s=O, and t=E;
(3.8c)fA(xi,yj,zk)=-f1(xi,yj,zk)-f2(xi,yj,zk)+f3(xi,yj,zk)+f4(xi,yj,zk)-f5(xi,yj,zk)-f6(xi,yj,zk)+f7(xi,yj,zk)+f8(xi,yj,zk),Case 4.
r
=
E, s=O, and t=O;
(3.8d)fA(xi,yj,zk)=f1(xi,yj,zk)+f2(xi,yj,zk)-f3(xi,yj,zk)-f4(xi,yj,zk)-f5(xi,yj,zk)-f6(xi,yj,zk)+f7(xi,yj,zk)+f8(xi,yj,zk),Case 5.
r
=
O, s=E, and t=E;
(3.8e)fA(xi,yj,zk)=-f1(xi,yj,zk)+f2(xi,yj,zk)+f3(xi,yj,zk)-f4(xi,yj,zk)-f5(xi,yj,zk)+f6(xi,yj,zk)+f7(xi,yj,zk)-f8(xi,yj,zk),Case 6.
r
=
O, s=E, and t=O;
(3.8f)fA(xi,yj,zk)=f1(xi,yj,zk)-f2(xi,yj,zk)-f3(xi,yj,zk)+f4(xi,yj,zk)-f5(xi,yj,zk)+f6(xi,yj,zk)+f7(xi,yj,zk)-f8(xi,yj,zk),Case 7.
r
=
O, s=O, and t=E;
(3.8g)fA(xi,yj,zk)=f1(xi,yj,zk)-f2(xi,yj,zk)+f3(xi,yj,zk)-f4(xi,yj,zk)+f5(xi,yj,zk)-f6(xi,yj,zk)+f7(xi,yj,zk)-f8(xi,yj,zk),Case 8.
r
=
O, s=O, and t=O;
(3.8h)fA(xi,yj,zk)=-f1(xi,yj,zk)+f2(xi,yj,zk)-f3(xi,yj,zk)+f4(xi,yj,zk)+f5(xi,yj,zk)-f6(xi,yj,zk)+f7(xi,yj,zk)-f8(xi,yj,zk),
where f1,f2,f3,f4,f5,f6,f7, and f8 refer to the image intensity function defined in the cubes from 1 to 8, respectively. The letters “E” and “O” are the acronyms of even and odd, respectively.As shown in Figure1, 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(r+s+t) for the object intensity function f(x,y,z) are defined as the projection of the function f(x,y,z) onto the monomial xryszt as follows:
(3.9)Grst=∭-∞∞xrysztf(x,y,z)dxdydz.
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:
(3.10)Grst=∭-1/31/3xrysztf(x,y,z)dxdydz.
Equation (3.10) could be written as follows:
(3.11)Gpqr=∑i=1N∑j=1N∑k=1NSpqr(xi,yj,zk)f(xi,yj,zk),
where
(3.12)Spqr(xi,yj,zk)=∫xi-(Δxi/2)xi+(Δxi/2)∫yj-(Δyj/2)yj+(Δyj/2)∫zk-(Δzk/2)zk+(Δzk/2)xrysztdxdydz.
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:(3.13a)Ir(i)=∫xi-(Δxi/2)xi+(Δxi/2)xrdx=1r+1[(xi+Δxi2)r+1-(xi-Δxi2)r+1](3.13b)Is(j)=∫yj-(Δyj/2)yj+(Δyj/2)ysdy=1s+1[(yj+Δyj2)s+1-(yj-Δyj2)s+1](3.13c)It(k)=∫zt-(Δzt/2)zt+(Δzt/2)ztdz=1t+1[(zt+Δzt2)t+1-(zt-Δzt2)t+1].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
(3.14)Gpqr=∑i=1⌊N/2⌋∑j=1⌊N/2⌋∑k=1⌊N/2⌋Ir(i)Is(j)It(k)fA(xi,yj,zk),
where fA(xi,yj,zk) is the augmented intensity function defined by using (3.8a)–(3.8h). For even value of N, the operator ⌊N/2⌋ equals to N/2, while it equals to (N-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(r+s+t) defined by (3.14) are computed in three separate steps by successive computation of the 1D sth-order moment for each row, followed by the 2D (r+s)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:(3.15a)Grst=∑k=1⌊N/2⌋It(zk)Rrsk,where(3.16a)Rrsk=∑i=1⌊N/2⌋YiskIr(xi),(3.16b)Yisk=∑j=1⌊N/2⌋Is(yj)fA(xi,yj,zk).


## 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 f(xi,yj,zk) is generated by using the Matlab8 statement f(xi,yj,zk)=r  and  (M,N,K) where 0≤f(xi,yj,zk)≤1 for all i,j,k with M=N=K=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).
Moment order
Conventional method [21]
Proposed method
Reduction percentage
1
0.629653
0.057412
90.8820%
3
6.742940
0.078032
98.8428%
5
32.336633
0.095085
99.7060%
10
212.330483
0.174056
99.9180%
12
386.394043
0.182601
99.9527%
15
847.327152
0.276992
99.9673%Figure 2
Elapsed CPU times for computing 3D Zernike moments (first numerical experiment).In the second experiment, a protein structure of dimensions70×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).
Moment order
Conventional method [21]
Proposed method
Reduction percentage
1
0.634301
0.068551
89.1927%
3
6.880023
0.080441
98.8308%
5
37.497548
0.132576
99.6464%
10
303.177357
0.198912
99.9344%
12
731.124886
0.282519
99.9614%
15
913.281611
0.342348
99.9625%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 coefficientsHn,ℓ,mrst. It is clear that the computation of 3D geometric moments is image-dependent, while the computation of Hn,ℓ,mrst is image-independent. Therefore, we concentrate on the complexity analysis of 3D geometric moments.For a 3D digital image of sizeN×N×N 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(4.1a)(Max+1)(Max+2)(Max+3)N36,additions(4.1b)(Max+1)(Max+2)(Max+3)N33,multiplications.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 Yiqk requires (N/2)3(Max+1) multiplications and (N/2)2(N/2-1)(Max+1) additions. The creation process of the matrix Rpqk using (3.16a) requires (N/2)((N/2)-1)(Max+1)(Max+2)/2 additions and (N/2)2(Max+1)(Max+2)/2 multiplications. Finally, the computational complexity of the 3D geometric moments using (3.15a) requires (N/2)(Max+1)2(Max+2)/2 multiplications and ((N/2)-1)(Max+1)2(Max+2)/2 additions. Therefore, computing the set of independent exact 3D geometric moments requires the following number of additions and multiplications:(4.2a)(Max+1)(Max+2)8[N2-2N+2Max(N-2)]+N2(N-2)(Max+1)8,additions(4.2b)(Max+1)(Max+2)8[N2+2N(Max+1)]+N3(Max+1)8,multiplications.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 N 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.
Image size and moment order
Conventional method [21]
Proposed method
No. of +
No. of *
No. of +
No. of *
Max = 5,N=64
14,680,064
29,360,128
27,063
222,144
Max = 10,N=64
74,973,184
149,946,368
91,388
451,264
Max = 15,N=64
213,909,504
427,819,008
206,088
733,184
Max = 20,N=64
464,257,024
928,514,048
382,788
1,079,904
Max = 5,N=128
1.0
e
+
009
*
0.1174
1.0
e
+
009
*
0.2349
103,383
1,666,944
Max = 10,N=128
1.0
e
+
009
*
0.5998
1.0
e
+
009
*
1.1996
329,868
3,200,384
Max = 15,N=128
1.0
e
+
009
*
1.7113
1.0
e
+
009
*
3.4226
709,128
4,890,624
Max = 20,N=128
1.0
e
+
009
*
3.7141
1.0
e
+
009
*
7.4281
1,264,788
6,761,664


## 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 coefficientsHn,ℓ,mrst. It is clear that the computation of 3D geometric moments is image-dependent, while the computation of Hn,ℓ,mrst is image-independent. Therefore, we concentrate on the complexity analysis of 3D geometric moments.For a 3D digital image of sizeN×N×N 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(4.1a)(Max+1)(Max+2)(Max+3)N36,additions(4.1b)(Max+1)(Max+2)(Max+3)N33,multiplications.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 Yiqk requires (N/2)3(Max+1) multiplications and (N/2)2(N/2-1)(Max+1) additions. The creation process of the matrix Rpqk using (3.16a) requires (N/2)((N/2)-1)(Max+1)(Max+2)/2 additions and (N/2)2(Max+1)(Max+2)/2 multiplications. Finally, the computational complexity of the 3D geometric moments using (3.15a) requires (N/2)(Max+1)2(Max+2)/2 multiplications and ((N/2)-1)(Max+1)2(Max+2)/2 additions. Therefore, computing the set of independent exact 3D geometric moments requires the following number of additions and multiplications:(4.2a)(Max+1)(Max+2)8[N2-2N+2Max(N-2)]+N2(N-2)(Max+1)8,additions(4.2b)(Max+1)(Max+2)8[N2+2N(Max+1)]+N3(Max+1)8,multiplications.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 N 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.
Image size and moment order
Conventional method [21]
Proposed method
No. of +
No. of *
No. of +
No. of *
Max = 5,N=64
14,680,064
29,360,128
27,063
222,144
Max = 10,N=64
74,973,184
149,946,368
91,388
451,264
Max = 15,N=64
213,909,504
427,819,008
206,088
733,184
Max = 20,N=64
464,257,024
928,514,048
382,788
1,079,904
Max = 5,N=128
1.0
e
+
009
*
0.1174
1.0
e
+
009
*
0.2349
103,383
1,666,944
Max = 10,N=128
1.0
e
+
009
*
0.5998
1.0
e
+
009
*
1.1996
329,868
3,200,384
Max = 15,N=128
1.0
e
+
009
*
1.7113
1.0
e
+
009
*
3.4226
709,128
4,890,624
Max = 20,N=128
1.0
e
+
009
*
3.7141
1.0
e
+
009
*
7.4281
1,264,788
6,761,664


## 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.



---

*Source: 353406-2012-10-15.xml*