Thai V Hoang

Salvatore Tabbone

Fast Generic Polar Harmonic Transforms

Keywords: Polar harmonic transforms, orthogonal moments, recurrence relation, shortest addition chain, Chebyshev polynomials, fast computation

Generic polar harmonic transforms have recently been proposed to extract rotation-invariant features from images and their usefulness has been demonstrated in a number of pattern recognition problems. However, direct computation of these transforms from their definition is inefficient and is usually slower than some efficient computation strategies that have been proposed for other competing methods. This paper presents a number of novel strategies to compute these transforms rapidly. The proposed methods are based on the inherent recurrence relations among complex exponential and trigonometric functions used in the definition of the radial and angular kernels of these transforms. The employment of these relations leads to recursive and addition chain-based strategies for fast computation of harmonic function-based kernels. Experimental results show that the proposed method is about 10× faster than direct computation and 5× faster than fast computation of Zernike moments using the q-recursive strategy. Thus, among all existing rotation-invariant feature extraction methods, polar harmonic transforms are the fastest.

I. INTRODUCTION

Image moments extracted from a unit disk region are usually used in invariant pattern recognition problems as rotationinvariant features [START_REF] Teague | Image analysis via the general theory of moments[END_REF]. Let f be the image function, its moment H nm of order n and repetition m (n, m ∈ Z) over the unit disk region is computed as

H nm = ẍ2 +y 2 ≤1 f (x, y)V * nm (x, y) dxdy, (1)
where V nm is the decomposing kernel and the asterisk denotes the complex conjugate. It was shown in [START_REF] Bhatia | On the circle polynomials of Zernike and related orthogonal sets[END_REF] that in order to have the "invariance in form" property on H nm , the polar form of the kernel, V nm (r, θ), should take the form R n (r)A m (θ) with A m (θ) = e imθ . In addition, orthogonality is often preferred in many practical situations to reduce information redundancy in compression, to avoid difficulties in image reconstruction, and to increase accuracy in pattern recognition. Orthogonality between kernels means V nm , V n m = ´1 0 R n (r)R * n (r)r dr ´2π 0 A m (θ)A * m (θ) dθ = δ nn δ mm . From the orthogonality between angular kernels:

ˆ2π 0 A m (θ)A * m (θ) dθ = 2πδ mm , (2)
the remaining condition for radial kernels is

ˆ1 0 R n (r)R * n (r)r dr = 1 2π δ nn . (3)
T. V. Hoang is with Inria Nancy -Grand Est, 615 rue du Jardin Botanique, 54600 Villers-lès-Nancy, France S. Tabbone is with Université de Lorraine, LORIA, Campus scientifique, 54506 Vandoeuvre-lès-Nancy, France

The above equation represents the condition for the definition of a set of radial kernels in order to have orthogonality between kernels. There exists a number of methods that satisfies this condition. One direction is to employ polynomials of the variable r for R n , which turn out to be special cases of the Jacobi polynomials [START_REF] Ping | Generic orthogonal moments: Jacobi-Fourier moments for invariant image description[END_REF], [START_REF] Hoang | Errata and comments on "Generic orthogonal moments: Jacobi-Fourier moments for invariant image description[END_REF]. Popular methods are Zernike moments (ZM) [START_REF] Khotanzad | Invariant image recognition by Zernike moments[END_REF], pseudo-Zernike moments (PZM) [START_REF] Teh | On image analysis by the methods of moments[END_REF], orthogonal Fourier-Mellin moments (OFMM) [START_REF] Sheng | Orthogonal Fourier-Mellin moments for invariant pattern recognition[END_REF], Chebyshev-Fourier moments [START_REF] Ping | Image description with Chebyshev-Fourier moments[END_REF], and pseudo Jacobi-Fourier moments (PJFM) [START_REF] Amu | Image analysis by pseudo-Jacobi (p = 4, q = 3)-Fourier moments[END_REF] (see [START_REF] Flusser | Moments and Moment Invariants in Pattern Recognition[END_REF]Section 6.3] or [START_REF] Hoang | Image representations for pattern recognition[END_REF]Section 3.1] for a comprehensive survey). Despite its popularity, this class of orthogonal moments involves factorial computation, which results in high computational complexity and numerical instability, which often limit its practical usefulness.

The second direction defines V nm as the eigenfunctions of the Laplacian ∇ 2 on the unit disk, similar to the interpretation of Fourier basis as the set of eigenfunctions of the Laplacian on a rectangular domain. These eigenfunctions are obtained by solving the Helmholtz equation, ∇ 2 V + λ 2 V = 0, in polar coordinates to have the radial kernels defined based on Bessel functions of the first and second kinds [START_REF] Bowman | Introduction to Bessel Functions[END_REF]. By imposing the condition in Equation ( 3), a class of orthogonal moments are obtained [START_REF] Wang | Rotational invariance based on Fourier analysis in polar and spherical coordinates[END_REF] and different boundary conditions have been used for the proposal of different methods with distinct definition of λ: Fourier-Bessel modes (FBM) [START_REF] Guan | Fourier-Bessel analysis of patterns in a circular domain[END_REF], Bessel-Fourier moments (BFM) [START_REF] Xiao | Image analysis by Bessel-Fourier moments[END_REF], and disk-harmonic coefficients (DHC) [START_REF] Verrall | Disk-harmonic coefficients for invariant pattern recognition[END_REF]. However, the main disadvantage of these methods is the lack of an explicit definition of their radial kernels other than Bessel functions. And this leads to inefficiency in terms of computation complexity.

The third direction uses harmonic functions to define R n by taking advantage of their orthogonality. Harmonic functions are known to possess elegant orthogonality relations which are "similar in form" to the condition in Equation [START_REF] Ping | Generic orthogonal moments: Jacobi-Fourier moments for invariant image description[END_REF]. For example, the orthogonality relation between complex exponential functions in Equation (2) "differs in form" from that in Equation ( 3) by a single multiplicative term r. Adaptation to overcome r was attempted and led to a number of proposals. Complex exponential functions are used in polar complex exponential transform (PCET) [START_REF] Yap | Two-dimensional polar harmonic transforms for invariant image representation[END_REF] and trigonometric functions are used in radial harmonic Fourier moments (RFHM) [START_REF] Ren | Multidistortion-invariant image recognition with radial harmonic Fourier moments[END_REF], polar cosine transform (PCT), and polar sine transform (PST) [START_REF] Yap | Two-dimensional polar harmonic transforms for invariant image representation[END_REF]. A recent study [START_REF] Hoang | Generic polar harmonic transforms for invariant image representation[END_REF] gives a generic view on the use of harmonic functions to define the radial kernels, leading to the proposal of generic polar harmonic transforms (GPHT). The study also demonstrates some beneficial properties of GPHT on representation and recognition capabilities. Since harmonic functions are so much simpler than Jacobi polynomials and Bessel functions, computing GPHT thus promises to be faster than computing unit disk-based moments with radial kernels defined based on Jacobi polynomials or Bessel functions. It should be noted that the generalization of polar harmonic transforms by introducing a parameter is similar to the generalization of the R-transform published recently [START_REF]The generalization of the R-transform for invariant pattern representation[END_REF].

Direct computation of GPHT from their definition is inefficient and may be slower than some efficient computation strategies that have been proposed for methods based on Jacobi polynomials [START_REF] Chong | A comparative analysis of algorithms for fast computation of Zernike moments[END_REF], [START_REF] Singh | Fast and numerically stable methods for the computation of Zernike moments[END_REF], [START_REF] Gautschi | Orthogonal Polynomials: Computation and Approximation[END_REF]. This paper presents a number of novel strategies for fast computation of GPHT in order to further support them in terms of computational complexity. The proposed methods are based on the inherent recurrence relations among harmonic functions that are used in the definitions of the radial and angular kernels of GPHT. The employment of these relations leads to recursive and addition chain-based strategies for fast computation of GPHT kernels. The proposed strategies have been compared with direct computation and, in terms of kernel computation time, the best strategy is about 10× faster. When compared with recursive computation of Zernike kernels using the q-recursive strategy [START_REF] Chong | A comparative analysis of algorithms for fast computation of Zernike moments[END_REF], the proposed strategy is also 5× faster.

The content of this paper is a comprehensive extension of a previous research work [START_REF] Hoang | Fast computation of orthogonal polar harmonic transforms[END_REF], which only presents recursive computation of complex exponential functions. A similar work on fast computation of polar harmonic transforms has recently presented in [START_REF] Singh | Fast computation of polar harmonic transforms[END_REF]. This work, however, focuses only on recursive computation of trigonometric functions and does not cover complex exponential functions, efficient computation flows, and addition chain-based strategy as in this paper. The next section will give some background on GPHT, their discrete approximation, and revisit a fast computation strategy based on geometrical symmetry. Sections III and IV present the recurrence relations among harmonic functions and propose some strategies for their fast computation using recursion or the shortest addition chain. Experimental results to validate the proposed strategies are then given in Section V and conclusions are finally drawn in Section VI.

II. THE GENERIC POLAR HARMONIC TRANSFORMS

A generic view on the use of harmonic functions for the definition of radial kernels has recently led to the proposal of four generic classes of polar harmonic transforms (GPHT) [START_REF] Hoang | Generic polar harmonic transforms for invariant image representation[END_REF]: generic polar complex exponential transform (GPCET), generic radial harmonic Fourier moment (GRHFM), generic polar cosine transform (GPCT), and generic polar sine transform (GPST). While still owning "invariance in form" and orthogonality properties, these transforms have been shown to have some other beneficial properties in terms of representation and recognition capabilities. The radial kernels of GPCET, GRHFM, GPCT, and GPST are respectively defined as

R ns (r) = sr s-2 2π e i2πnr s R H ns (r) = sr s-2 2π      1, n = 0 √ 2 sin(π(n + 1)r s ), n odd > 0 √ 2 cos(πnr s ), n even > 0 R C ns (r) = sr s-2 2π 1, n = 0 √ 2 cos(πnr s ), n > 0 m = 0 m = 1 m = 2 n = 0 n = 1 n = 2 (a) s = 0.5 m = 0 m = 1 m = 2 n = 0 n = 1 n = 2 (b) s = 1 m = 0 m = 1 m = 2 n = 0 n = 1 n = 2 (c) s = 2 m = 0 m = 1 m = 2 n = 0 n = 1 n = 2 (d) s = 4
{(n, m) ∈ Z 2 : 0 ≤ n, m ≤ 2}. R S ns (r) = sr s-2 2π √ 2 sin(πnr s ), n > 0.
It can be seen that these definitions of radial kernels satisfy the condition in Equation (3). Also, GPCET is equivalent to GRHFM in terms of representation capability, similar to the equivalence between different forms of Fourier series (namely complex exponential and trigonometric functions).

In addition, the introduction of a parameter s ∈ Z [START_REF] Hoang | Generic polar harmonic transforms for invariant image description[END_REF] makes these transforms generic versions of existing polar harmonic transforms: PCET, RHFM, PCT, and PST are obtained from GPCET, RHFM, GPCT, and GPST when s = 2, 1, 2, 2, respectively. The two-dimensional view of the phase of some GPCET kernels is illustrated in Figure 1. Note that all radial and angular kernels of GPHT are defined based on complex exponential and trigonometric functions. Direct computation of these functions is time-consuming and constitutes a dominant part of the computation of GPHT due to their O(log 2 n) complexity, where n is the number of precision digits at which the functions are evaluated [START_REF] Borwein | On the complexity of familiar functions and numbers[END_REF]. The overall complexity may become excessively high when a large number of moments is needed, or large-sized images are used, or a highprecision computation is required. Since these requirements are common in practical applications, the existence of strategies for fast computation of GPHT kernels is vital for the applicability of GPHT. The remaining of this section will discuss the discrete approximation of GPHT and revisit a fast computation strategy based on geometrical symmetry. The next two sections (Sections III and IV) will propose some strategies to compute GPHT rapidly based on the inherent recurrence relations among harmonic functions.

Discrete approximation

The image processed in digital systems is not defined in a continuous domain, but a grid of pixels. Over each pixel region, a constant value is used to represent the intensity of the image function at that pixel location. Accordingly, the continuous integration in Equation (1) needs to be approximated in two separate ways: in the domain of integration (geometrically) and in the kernel values (numerically). To be precise, a lattice of pixels is selected to cover the unit disk for geometric approximation, and some values are assigned to each kernel at each pixel location for numerical approximation.

Geometric approximation is required to linearly map a rectangular grid of pixels onto a continuous and circular domain of the unit disk since an exact mapping is impossible due to their different geometric nature. In this work, the incircle mapping strategy [START_REF] Liao | On the accuracy of Zernike moments for image analysis[END_REF] which puts a disk inside the image's rectangular domain as illustrated in Figure 2 is used. In this way, most of the circular region is involved in the computation. However, some image pixels that lie outside the circular region (the white pixels) will not be used in the computation. Using incircle also leads to the consideration of image pixels that lie at the disk boundary since using these pixels causes an inherent error, called geometric error in the literature. In the figure, the blue and cyan pixels at the disk boundary have their centers lying outside and inside the circular region, respectively. While the error caused by the cyan pixels can be fixed, that by the blue pixels cannot be fixed and remains the main source of geometric error [START_REF] Pawlak | On the recovery of a function on a circular domain[END_REF]. The impact of this error type depends heavily on the behavior of the radial kernels around the point r = 1. Let f be the digital image of size M × N , the domain

[0, M ] × [0, N ] of f is mapped onto [-1, 1] × [-1, 1] by i -→ 2i -M M = i∆x -1, 0 ≤ i ≤ M j -→ 2j -N N = j∆y -1, 0 ≤ j ≤ N
where ∆x = 2 M and ∆y = 2 N . In order to avoid geometric error in this work, only the image pixels that lie entirely inside the unit disk (the yellow pixels) are used in the computation. The values of the image function over the intersection regions between the unit disk and the blue/cyan pixels are then assumed to take value 0. Note that an algorithm to compute ZM that eliminates the geometric error by performing interpolation and then doing computation in polar coordinates was proposed in [START_REF] Xin | Accurate computation of Zernike moments in polar coordinates[END_REF]. However, for simplicity and to avoid interpolation error, this algorithm is not used in this study.

Numerical approximation arises naturally when GPHT are computed from the approximated lattice of pixels over the unit disk. Let [i, j] denote one pixel region in the domain of f having a constant intensity value f [i, j], its mapped region in the unit disk is x i -∆x 2 , x i + ∆x 2 × y j -∆y 2 , y j + ∆y 2 , where (x i , y j ) are the coordinates of the center of the mapped region. The condition for this mapped region to be labeled yellow color is that its four corners lie inside the unit disk. Mathematically, the set C of pixels that satisfies this condition is defined as C = [i, j] :

x i ± ∆x 2 2 + y j ± ∆y 2 2 ≤ 1 .
The integration in Equation ( 1) can then be approximated as a discrete sum indexed by pixels in C:

Ĥnm = [i,j]∈C f [i, j] h nm [i, j],
where Ĥnm is an estimate of H nm and

h nm [i, j] = ˆxi+ ∆x 2 xi-∆x 2 ˆyj+ ∆y 2 yj -∆y 2 V * nm (x, y) dxdy (4)
represents the accumulation of V * nm over a region of size ∆x× ∆y in the unit disk that represents the pixel [i, j] (see Figure 3a). Since V nm is originally defined in polar coordinates while h nm [i, j] is evaluated in Cartesian ones, computing h nm [i, j] in practice usually relies on numerical integration techniques.

There exist many such techniques that can be employed to approximately evaluate h nm [i, j] [START_REF] Liao | On the accuracy of Zernike moments for image analysis[END_REF]. The simplest is the rectangle formula (scheme M 1 in Figure 3b) written as

h nm [i, j] V * nm (x i , y j )∆x∆y. (5)
This approximation assumes that the value of V nm over the mapped region of the pixel [i, j] is constant and equals its value at the central point (x i , y j ). The order of approximation error in this case is O (∆x∆y) 2 , thus depends on the area of the mapped region. Therefore, a more accurate approximation could be obtained by "pseudo" up-sampling the pixel [i, j] and then using Equation ( 5) on each up-sampled pixel (the L = 3 × 3 up-sampling is shown in scheme M 2 in Figure 3b). Other approaches use the L-dimensional cubature formulas [START_REF] Engels | Numerical Quadrature and Cubature[END_REF] defined as

h nm [i, j] L k=1 w k V * nm (u k , v k )∆x∆y,
where {(u k , v k ) : 1 ≤ k ≤ L} is a set of design points belonging to the mapped region that represents the pixel [i, j] and {w k : 1 ≤ k ≤ L} is a set of the corresponding weights. As an example, a five-dimensional cubature formula (scheme M 3 in Figure 3b) has the definition

h nm [i, j] 1 3 -V * nm (x i , y j ) + V * nm x i + ∆x 4 , y j +V * nm x i -∆x 4 , y j +V * nm x i , y j + ∆x 4 + V * nm x i , y j -∆x 4
∆x∆y. The order of approximation error now reduces to O (∆x∆y) L+1 . Thus, a higher accuracy is obtained when a larger value of L is used. Note that the gain in accuracy has to be paid by an increase in computational complexity, which is not always preferred in real applications. In this work, for simplicity and for the purpose of comparison, the scheme M 1 in Figure 3b is used to compute h nm [i, j] in GPHT and comparison methods.

In the literature, the aforementioned approximation error in computing h nm [i, j] is often called numerical error [START_REF] Liao | On the accuracy of Zernike moments for image analysis[END_REF]. The impact of this type of error depends heavily on the behavior of V nm over the mapped region that represents the pixel [i, j] and thus varies according to the actual moment and the orders n, m. In the case of GPHT, the radial and angular kernels that constitute V nm are defined based on harmonic functions and oscillate within their respective intervals, i.e., [0, 1] and [0, 2π). This means that V nm will oscillate at a higher frequency when n and m take higher values. In addition, since the aforementioned approximation schemes are based on the sampled values of V nm , the value of H nm is susceptible to information loss if V nm has been under-sampled. This is because the number of sampling points is fixed and determined by the size of the input image. Thus, when a moment of a too-high order n or m is considered, a fixed sampling of V nm becomes insufficient.

There exist some approaches that try to eliminate numerical error by avoiding direct approximation of h nm [i, j] over each mapped region of the pixel [i, j]. Some stick to the Cartesian space and perform the integration in Equation (4) analytically after converting V nm from being defined by polar coordinates to being defined by Cartesian ones [START_REF] Kotoulas | Accurate calculation of image moments[END_REF], [START_REF] Wee | On the computational aspects of Zernike moments[END_REF]. This is possible due to the relationship between Zernike and geometric moments [START_REF] Teh | On image analysis by the methods of moments[END_REF] or by piecewise polynomial interpolation of V nm . Another approach [START_REF] Xin | Accurate computation of Zernike moments in polar coordinates[END_REF] computes the moments directly in polar space by using bicubic interpolation to convert the digital image f from being defined as a piecewise constant function in the Cartesian space to being defined as a piecewise constant function in the polar space. Optimization techniques were recently used to improve the orthogonality of the approximated discrete kernels for more accurate image reconstruction [START_REF] Lin | Orthogonal rotation-invariant moments for digital image processing[END_REF]. However, these approaches to eliminate numerical error are much more computationally expensive than the aforementioned approximation schemes because interpolation makes any recurrence relation between radial and/or angular kernels invalid and also optimization requires iterative evaluations. In addition, these approaches also introduce new error into the value of H nm due to the interpolation or optimization process. For these reasons, in this work, these approaches are not used to compute GPHT and comparison methods.

Geometrical symmetry-based fast computation

It can be seen from Figure 4 that the seven symmetrical points of a point P 1 across the y-axis, the origin, the x-axis, x-axis

(x, -y) r -θ Q 1 y = x (y, x) r π 2 -θ Q 2 y = x, y-axis (-y, x) r π 2 + θ Q 3 y = x, origin (-y, -x) r -π 2 -θ Q 4 y = x, x-axis (y, -x) r -π 2 + θ
and the line y = x in the Cartesian coordinate system have the same radial and related angular coordinates as those of P 1 (see Table I). From the properties of complex exponential functions, the radial kernel takes the same value and the angular kernel has related values at these points. Thus, geometrical symmetry in the distribution of pixels inside the unit disk can be exploited to reduce the need of computing the radial and angular kernels for all pixels inside the unit disk to only pixels in one of the eight equal pies in Figure 4. This idea was initially proposed in [START_REF] Hwang | A novel approach to the fast computation of Zernike moments[END_REF] for Zernike moments and is applicable to all unit diskbased moments. In addition, due to its geometrical nature, this strategy is independent of the aforementioned pixel-based discrete approximation and thus can be used in combination with strategies for fast computation of radial/angular kernels proposed in the following two sections to further reduce the computational complexity. This observation will have experimental evidence in Section V.

III. FAST COMPUTATION OF EXPONENTIAL FUNCTIONS

Complex exponential functions are used in the definition of GPCET radial kernels and GPHT angular kernels. This section presents two different strategies to compute these functions efficiently using recursion and the shortest addition chain.

Recursive exponentiation

From the following recursive definition of exponentiation: base case: e 0α = 1, inductive clause: e kα = e (k-1)α e α , k ∈ Z, the complex exponential functions in the definition of GPCET radial kernels, R ns , and GPHT angular kernels, A m , can be computed recursively as

e i2πnr s = e i2π(n-1)r s e i2πr s , (6)
e imθ = e i(m-1)θ e iθ , using the base cases e i2π0r s = 1 and e i0θ = 1, respectively.

Assuming that

sr s-2
2π , e i2πr s , e iθ have been pre-computed and stored in memory for the polar coordinates (r, θ) of all pixel centers mapped inside the unit disk, the following recurrence relations of R ns and A m :

R ns (r) = sr s-2 2π e i2πnr s = R (n-1)s (r) e i2πr s , (7)
A m (θ) = e imθ = A m-1 (θ) e iθ , (8)
lead to their recursive computation with the base cases R 0s (r) =

sr s-2 2π
and A 0 (θ) = 1, respectively. Thus, computing R ns from R (n-1)s and A m from A m-1 each requires only one multiplication for each pixel location, which is very fast when compared to direct exponentiation. It should be noted that these forms of recurrence relations are simpler than those that were previously proposed for the radial kernels defined based on Jacobi polynomials (see [START_REF] Chong | A comparative analysis of algorithms for fast computation of Zernike moments[END_REF], [START_REF] Singh | Fast and numerically stable methods for the computation of Zernike moments[END_REF] and references therein). By using Equations ( 7) and ( 8), only two recursive computational threads are needed to reach every GPCET radial kernel and GPHT angular kernel, whereas many threads would be required to cover all radial kernels defined based on Jacobi polynomials. The corresponding computation flows for R ns and A m are illustrated in Figure 5. Note also that the method proposed here is much faster than the one mentioned in [START_REF] Yap | Two-dimensional polar harmonic transforms for invariant image representation[END_REF] where direct exponentiation is used to compute PCET radial and angular kernels.

The aforementioned recursive computation of R ns and A m could be employed for extremely fast computation of GPCET kernels when the order set S composes a square region in Z 2 and takes the origin as its center as

S = {(n, m) : n, m ∈ Z 2 , |n|, |m| ≤ K},
where K is a certain positive integer. By using the computational flow depicted in Figure 6a, computing V nms (r, θ) requires only two multiplications: one for updating R ns (r) or A m (θ) and one for computing V nms (r, θ) = R ns (r)A m (θ). As a result, computing a GPCET moment requires only three vector multiplications (two for computing V nms and one for multiplying V nms by f ), followed by a discrete sum over all the pixel locations. Vnms thus requires two multiplications, one for computing Rns or Am using Equations ( 7) or ( 8), depending on whether the value of n or m has just been changed, and the other for multiplying Rns by Am to get Vnms.

To further boost the computation speed, instead of letting the computational flow to visit all (n, m) ∈ S in the four quadrants in the Cartesian space, it is sufficient to visit only (n, m) ∈ S in one quadrant (n, m ≥ 0) as illustrated in Figure 6b. This is possible due to the following relations:

R -ns (r) = R * ns (r), A -m (θ) = A * m (θ).
Thus, whenever R ns and A m are available, computing the four related GPCET kernels requires only two multiplications and three conjugations:

V nms (r, θ) = R ns (r)A m (θ), V -nms (r, θ) = R * ns (r)A m (θ), V -n-ms (r, θ) = V * nms (r, θ), V n-ms (r, θ) = V * -nms (r, θ).
Since eight multiplications are needed to compute these four kernels if the computational flow in Figure 6a is used, the simple trick of staying in one quadrant leads to a 8 3 -time reduction in the number of multiplications.

When there is enough memory space to store all the radial kernels R ns (|n| ≤ K) and angular kernels A m (|m| ≤ K), then pre-computing and keeping their values in memory may lead to a further reduction in computational complexity by using the computational flow in Figure 7. It can be seen that, for each pixel location, the total number of multiplications required to update the values of R ns or A m using the computational flow in Figure 6a is (K + 1) 2 -1 whereas pre-computation using the computational flow in Figure 7 needs only 2K multiplications. However, besides the memory requirement to store pre-computed values, this strategy has a disadvantage concerning the data accessing time. The pre-computed values of all R ns (or A m ) are often stored in a matrix R s (or A) indexed by r and n (or θ and m). When the kernel V nms needs to be computed, the values in the nth column of R s and in the mth column of A are retrieved from memory. Since fetching data from memory is less efficient than directly using data in CPU cache as in the case of recursive computation, the pre- computing strategy may not provide any computational speedup as expected. This concern will be supported by experimental evidence in Section V.

In real situations, the speed-up in computing GPCET kernels leads to the possibility of increasing the number of GPCET moments without changing the system throughput. This increase is equivalent to an increase in the number of features or, more importantly, an increase in the number of pattern classes the system can discriminate. Another by-product of fast computation is the reduction in storage space requirement. In methods based on orthogonal polynomials the common practice is to pre-compute and store kernels in memory in order to avoid their repetitive expensive computation. For example, if L moments are needed then the L corresponding kernels need to be pre-computed and stored. In the case of GPCET, since the kernel of any order can be computed rapidly, there is merely a need to store the set

sr s-2
2π , e i2πr s , e iθ for the polar coordinates (r, θ) of all the pixel centers mapped into the unit disk, regardless the required number of moments. This practice also leads to a reduction in storage space by 2K 3 .

Addition-chain-based exponentiation

In some situations only some kernels V nms of certain orders n and m are needed, the recursive strategies proposed previously in this section are not optimal since all R ps and A q with 0 ≤ p ≤ n and 0 ≤ q ≤ m need to be computed (Figure 5), leading to n + m + 1 multiplications for each pixel location. There exists in the literature methods that allow computing R ns from R 0s and A m from A 0 quickly using addition-chain exponentiation [START_REF] Gordon | A survey of fast exponentiation methods[END_REF]. The basic idea is based on the minimal number of multiplications necessary to compute g n from g (n ∈ Z), given that the only permitted operation is multiplication of two already-computed powers of g. In fact, this minimal number is the length of the shortest addition chain for n, which is a list of positive integers a 1 , a 2 , . . . , a with a 1 = 1 and a = n such that for each 1 < i ≤ n, there are some j and k satisfying 1 ≤ j ≤ k < i and a i = a j + a k . It thus can be seen that a shorter addition chain for n results in computation of g n computing g a2 , g a3 , . . . g -1 , g n Algorithm Computation of g n by the "square-and-multiply" exponentiation method. 1: a ← 1 2: for i = → 0 do 3:

a ← a * a 4:

if c i = 1 then 5: a = a * g 6:
end if 7: end for in sequence. While the shortest addition chain for a relatively small number can be obtained from lookup table [START_REF] Knuth | The Art of Computer Programming[END_REF], finding this chain for a large number is known as an NP-complete problem [START_REF] Downey | Computing sequences with addition chains[END_REF]. For this reason, a near-optimal solution is often used as a substitution in practical situation.

A number of methods has been proposed to find nearoptimal addition chains for a number. One of the simplest and most frequently used is the binary method, also known as the "square-and-multiply" method, which is based on the following formula to compute g n (n > 0) recursively using squaring and multiplication as

g n =        1, if n = 0, g × g n-1 2 2 , if n is odd g n 2 2 , if n is even Let n = i=0 c i 2 i
be the binary expansion of a positive integer n, the pseudo-code for the computation of g n (n > 0) using the "square-and-multiply" method is given in Algorithm 1. At the beginning of each iteration i of the "for loop", a is equal to g t , where t has its binary expansion equal to the first -i prefix of the binary expansion of n. Squaring a has the effect of doubling t, and multiplying by g puts a 1 in the last digit if the corresponding bit c i is 1. It can be seen that this algorithm uses log n squaring operations and at most log n multiplications. In other words, the "square-and-multiply" method takes 2 log n multiplications in the worst case and 3 2 log n on average. Since the lower bound for the number of multiplications required to do a single exponentiation is log n, this method is usually considered as "good enough" and is computationally much more efficient than the recursive method of multiplying the base g with itself repeatedly to get g n , requiring n multiplications. The application of this "square-andmultiply" method for the computation of GPCET radial kernels and GPHT angular kernels when only some moments of certain orders need to be computed will speedup the computation. This observation will be supported by experimental evidence in Section V.

Comparison with the existing recursion-based method

In order to evaluate the improvement in computational complexity obtained by the proposed method over the method recently presented in [START_REF] Singh | Fast computation of polar harmonic transforms[END_REF], a theoretical comparison on the computational complexity of PCET (or GPCET at s = 2) has been carried out and the results are summarized in Table II. Note that [START_REF] Singh | Fast computation of polar harmonic transforms[END_REF] only presents recursive computation of complex exponential functions for the angular and radial kernels of PCET by employing the recurrence relations of trigonometric functions. It does not have efficient computation flows as in Figure 6 and the addition chain-based exponentiation presented in this paper. So a direct comparison with the method in [START_REF] Singh | Fast computation of polar harmonic transforms[END_REF] should only involve Equations ( 7), ( 8) and Figure 5.

Nevertheless, the comparison results demonstrate clearly the superiority of the proposed method over the method in [START_REF] Singh | Fast computation of polar harmonic transforms[END_REF].

IV. FAST COMPUTATION OF TRIGONOMETRIC FUNCTIONS

The trigonometric functions used in the definition of the radial kernels of GRHFM, GPCT, and GPST can also be computed rapidly. This is because cosine and sine functions are equivalent to complex exponential functions in terms of computational complexity due to cos(πnr s ) = Re{e iπnr s } = e iπnr s + e -iπnr s 2 , sin(πnr s ) = Im{e iπnr s } = e iπnr s -e -iπnr s 2i .

By using these identities, the recursive strategies to compute GPCET radial kernels and GPHT angular kernels presented in Section III can also be employed to compute the radial kernels GRHFM, GPCT, and GPST. In other words, the implementation of GRHFM, GPCT, and GPST can rely on that of GPCET for low computational complexity. In this way, all the computational gains claimed for GPCET in the previous section using and addition chain are still valid for GRHFM, GPCT, and GPST.

Recursive computation

Apart from relying on complex exponential functions, cosine and sine functions could also be fast computed by using the following recurrence relations:

T 0 (x) = 1 T 1 (x) = x T n+1 (x) = 2x T n (x) -T n-1 (x), for n ≥ 1, (9)
and

U 0 (x) = 1 U 1 (x) = 2x U n+1 (x) = 2x U n (x) -U n-1 (x), for n ≥ 1 (10) with
T n (cos(πr s )) = cos(πnr s ),

U n (cos(πr s )) = sin(π(n + 1)r s ) sin(πr s ) sin(πnr s ) = U n-1 (cos(πr s )) sin(πr s ), where T n and U n are the Chebyshev polynomials of the first and second kinds, respectively [START_REF] Koornwinder | Orthogonal polynomials[END_REF]. Accordingly, the computation of cos(πnr s ) and sin(πnr s ) can also be carried out recursively with the base cases {1, cos(πr s )} for Equation ( 9) and {1, 2 cos(πr s )} for Equation [START_REF] Flusser | Moments and Moment Invariants in Pattern Recognition[END_REF], provided that 2 cos(πr s ) and sin(πr s ) have been pre-computed and are available throughout the computation process. The computation flow for GRHFM, GPCT, and GPST radial kernels are illustrated in Figure 8. It can be seen that computing cos(πnr s ) and sin(πnr s ) each requires only one multiplication and one subtraction, which is almost equivalent to one multiplication required to compute e i2πnr s as in Equation ( 6).

Addition-chain-based computation

Fast computation of cosine and sine functions for certain n-multiple angles using the "square-and-multiply" method is also possible by using the recurrence relations of Chebyshev polynomials of the first and second kinds [40, page 782]

T n1+n2 (x) + T n1-n2 (x) = 2 T n1 (x) T n2 (x), U 2n-1 (x) = 2 T n (x) U n-1 (x), U n1+n2-1 (x) -U n1-n2-1 (x) = 2 T n1 (x) U n2-1 (x). By letting n 2 = n 1 = n and n 2 = n 1 -1 = n -1: T 2n (x) = 2 T 2 n (x) -1, T 2n-1 (x) = 2 T n (x) T n-1 (x) -x, U 2n-1 (x) = 2 T n (x) U n-1 (x), U 2n-2 (x) = 2 T n (x) U n-2 (x) + 1.
These identities are "similar in form" to g 2n = (g n ) 2 and g 2n-1 = g (g n-1 ) 2 , which are the bases for the "squareand-multiply" exponentiation. It can also be seen here that using addition chain to compute cos(πnr s ) and sin(πnr s ) is computationally more efficient than using recursion.

V. EXPERIMENTAL RESULTS

The computational complexity of existing unit disk-based orthogonal moments is evaluated in terms of the elapsed time required to compute their kernels/moments from an image of size 128 × 128. For this image size, the incircle illustrated in Figure 2 contains 12596 pixels that lie entirely inside the unit disk (the yellow pixels). Experiments are performed on a single core of a PC equipped with a 2.33 GHz CPU, 4 GB RAM, and running Linux kernel 2.6.38. MATLAB version 7.7 (R2008b) is used as the programming environment. Let K be a certain integer constant, all the kernels/moments of orders (n, m) of each method that satisfy the conditions listed in Table III are computed and the average elapsed time over all feasible orders is used as the kernel/moment computation time at that value of K. In order to study the trends in the dependence of kernel/moment computation time on the maximal kernel/moment order K, the value of K is varied in the range 0 ≤ K ≤ 20 in all experiments. In addition,

Method

Order range

ZM |m| ≤ n ≤ K, n -|m| = even PZM |m| ≤ n ≤ K OFMM/CHFM/PJFM 0 ≤ |m|, n ≤ K FBM/BFM/DHC 0 ≤ |m|, n ≤ K GPCET |m|, |n| ≤ K GRHFM |m| ≤ K, 0 ≤ n ≤ 2K GPCT 0 ≤ |m|, n ≤ K GPST |m| ≤ K, 1 ≤ n ≤ K
for more reliable results, all the computation times reported in this section are the average values over 100 trials. This section reports the computation time of kernel/moment by using definitions (Section V-A), recurrence relations (Section V-B), and the shortest addition chains (Section V-C). For fast computation GPHT, the experiments focus on complex exponential functions (i.e., GPCET) since harmonic functions are equivalent in terms of computational complexity. This is because the recurrence relations of complex exponential and trigonometric functions developed in Sections III and IV are similar. Thus, all conclusions drawn in this section for GPCET are also valid for GRHFM, GPCT, and GPST.

A. Direct computation

In this experiment, the kernels of all unit disk-based orthogonal moments are computed using their corresponding definitions, no recursive strategy is used. The computation times per kernel in milliseconds of comparison methods at different values of K are given in Figure 9. It can be seen from the figure that Jacobi polynomial-based methods (ZM, PZM, OFMM, CHFM, PJFM) have their kernel computation times increase almost linearly with the increase in K, meaning that a longer time is needed to compute a kernel of a higher order. This is because when K increases, factorials of larger integers need to be evaluated and more additive terms in the final summation need to be computed. Among Jacobi polynomialbased methods, OFMM and PJFM have the highest and similar complexity while ZM has the lowest. This complexity ranking of these methods is consistent with the ranking according to the number of multiplications required to compute their radial kernels by using definition.

Eigenfunction-based methods (FBM, BFM, DHC) require the longest times to compute their kernels over all comparison moments. Among these methods, FBM and DHC have almost the same complexity whereas that of BFM is slightly less. The high complexity of eigenfunction-based methods can be partly Kernel computation times of comparison methods by direct computation at different values of K. For each method and at a specific value of K, the kernel orders used in the computation should satisfy the conditions described in Table III and the average elapsed time over all feasible orders is used as the kernel computation time.

explained by the need of finding the zeros of (the derivative of) Bessel functions in order to define their radial kernels. And, there exists no formulation that allows fast computation of Bessel functions in the literature. In addition, the lower complexity of BFM is due to its use of a fixed-order Bessel function while FBM and DHC use Bessel functions of different orders. In this experiment, the evaluation of Bessel functions is facilitated by the MATLAB built-in function besselj.

Harmonic function-based methods (GPCET, GRHFM, GPCT, GPST) each requires an almost-constant time to compute its kernels of different orders. The constant kernel computation time is the direct result of the similarity in the form of their radial kernels at different orders. A change in the kernel order only leads to a change in the input to the complex exponential function (GPCET) or cosine/sine functions (GRHFM, GPCT, GPST) and, as a result, does not affect the kernel computation time. It can also be seen from Figure 9 that the kernel computation times of harmonic function-based methods are almost the same. This is because complex exponential, cosine, and sine functions are equivalent in terms of computational complexity [START_REF] Borwein | On the complexity of familiar functions and numbers[END_REF].

From the above observations on the direct kernel computation time of all unit disk-based orthogonal moments, it can be concluded here that the simple, resembling, and relating definitions of harmonic function-based kernels have resulted in an almostconstant kernel computation time, regardless of the maximal kernel order K. This makes a strong contrast with Jacobi polynomial-based and eigenfunction-based methods where a higher kernel order K means a longer kernel computation time.

B. Recursive computation

The proposed recursive strategies for fast computation of GPCET kernels are evaluated and compared with its direct computation and fast computation strategies for ZM kernels. The reason for comparing only to ZM lies in the lack of benchmarks on fast computation strategies for other methods and the popularity of recursive strategies for ZM computation in the literature. In this comparison, ZM kernels are computed using the q-recursive strategy [START_REF] Chong | A comparative analysis of algorithms for fast computation of Zernike moments[END_REF]. The comparison results are given in Figure 10 where the legends recursiveFour, recursiveOne, and recursivePre denote recursive computation of GPCET kernels using the computational flows in Figures 6a,6b, and 7, respectively. It can be seen from Figure 10 that ZM kernel computation time by using q-recursive is not constant. It increases in the range K = 0 → 5 and gradually decreases when K > 5. This is because the q-recursive strategy, which requires pre-computing two ZM radial kernels R nn and R n(n-2) for each order n, is applicable only when K ≥ 4 and is profitable when K ≥ 6. In addition, as K increases, more ZM radial kernels are computed by recursion and thus the proportion of pre-computed radial kernels decreases. This finally leads to a decrease in the average computation time of ZM radial kernels as K increases.

Using recursiveFour and recursiveOne to compute GPCET kernels leads to an almost-constant computation time, regardless of the maximal kernel order K. recursiveOne can be seen to be almost two-time faster than recursiveFour. The only exception is at K = 1 where the computation time suddenly drops below the average value. This might be because MATLAB simply copies the pre-computed value of e iθ into A 1 (θ) since A 0 (θ) = 1, instead of performing the more complex multiplication in order to optimize the calculation. A constant computation time of GPCET kernel is due to the fact that the recurrence relations in Equations ( 7) and ( 8) do not depend on the kernel orders and that there is no need to pre-compute GPCET radial kernels of any order as in the ZM q-recursive strategy. The "purely" recursive computation of radial and angular kernels is a distinct characteristic of harmonic functionbased methods. It can also be seen from the figure that the performance of recursivePre is not as good as expected; it is not faster than recursiveOne. The kernel computation time of recursivePre, unlike that of recursiveOne and recursiveFour, increases with the increase in K. The only possible explication for this phenomenon is the added burden of accessing columnwise data from the pre-computed matrices R s and A as described in Section III, instead of benefiting from CPU caching. On average, recursive computation of GPCET kernels using recursiveOne is approximately 10× faster than their direct computation and 5× faster than recursive computation of ZM kernels by using the q-recursive strategy.

In computing GPCET moments, that is a combination of computing GPCET kernels and performing the inner product between these kernels and image, recursiveOne could be combined with the strategy based on geometrical symmetry (Section II) to further reduce the moment computation time. Figure 11 provides the average elapsed times to compute one GPCET moment using recursiveOne without and with geometrical symmetry (symGeo). It can be seen from the figure that the combination of recursiveOne with symGeo is on average 7× faster than using recursiveOne alone, leading to an effective reduction in the computation time per moment from 0.4406 ms to 0.0634 ms at K = 20. These results clearly demonstrate that the strategies for fast computation of GPCET moments based on recursive computation of complex exponential functions and on geometrical symmetry are independent. Combining them definitely leads to a multiplication of the computational gains obtained individually by each of these two strategies.

C. Addition chain-based computation

Performance of fast computation of GPCET kernels using addition chain has also been evaluated at some selected values of order n. In this experiment, assuming that at each value of n, only the kernel V nns is computed, not all V n m s with |n |, |m | ≤ n as in the previous experiments. In the case of recursion, all the radial kernels R n s and angular kernels A m with n , m ≤ n need to be computed and the number of multiplications required to compute V nns at each pixel location is 2n + 1. On the contrary, the number of multiplications required to compute V nns at each pixel location in the case of addition chain is only 2 log n + 1.

The comparison results between recursive and addition chainbased computation for n = 1, 2, 4, 8, 16, 32, 64 are given in Figure 12. It can be seen from the figure that as n increases, the GPCET kernel computation time using recursion increases almost linearly with the increase in n. This linear relationship agrees with the recursion scheme illustrated in Figure 5 where a kernel of a higher order requires more recursion steps to compute. On the contrary, the GPCET kernel computation time using addition chain only increases logarithmically with the increase in n, agreeing with the theoretical estimation on the complexity of addition chain-based method. The benefit of using addition chain becomes more evident at a higher value of n. Thus, in situations where only harmonic kernels of certain orders need to be computed, the addition chain-based strategy should be used instead of the recursive strategy. VI. CONCLUSIONS Fast computation strategies for GPHT have been proposed in this paper by exploiting the recurrence relations between harmonic functions, leading to a method that is approximately 10× faster than direct computation. When compared with the q-recursive strategy for fast computation of ZM kernels, the proposed method is also approximately 5× faster. In addition, combination of the proposed method with the method based on geometrical symmetry leads to a multiplication of the computational gains obtained individually by the two methods. GPHT thus provide the fastest way to extract rotation-invariant features from images. Although other important issues such as robustness to discretization, noise, blurring, and symmetry preservation [START_REF] Pawlak | Image Analysis by Moments: Reconstruction and Computational Aspects[END_REF], [START_REF] Bissantz | Testing for image symmetries -with application to confocal microscopy[END_REF] need further studies, GPHT are the promising replacements of ZM in image analysis applications where low computational complexity is an important methodselection criterion.

It is worth noting here that the value of the parameter s only slightly affects the direct computation of GPHT and does not affect at all the proposed computation strategies. This is because s has no role in the recurrence relations in Equations ( 7) and ( 8), even though it appears in the definition of GPHT radial kernels. As a result, all the conclusions drawn in this paper hold for every s. In addition, the methods proposed in this paper can also be employed to compute the angular kernels of all existing unit disk-based moments and the radial kernels of ART [START_REF] Bober | Shape descriptors[END_REF] and GFD [START_REF] Zhang | Shape-based image retrieval using generic Fourier descriptor[END_REF] since they are defined based on harmonic functions albeit the radial kernels of ART and GFD do not satisfy the condition in Equation (3).

Figure 1 .

 1 Figure 1. 2D view of the phase of GPCET kernels using s = 0.5, 1, 2, 4 for {(n, m) ∈ Z 2 : 0 ≤ n, m ≤ 2}.

Figure 2 .Figure 3 .

 23 Figure 2. Lattice-point approximations of a circular region of an image of 32 × 32 pixels using incircle.

Figure 4 .

 4 Figure 4. Symmetrical points of a point P 1 inside the unit disk across the y-axis, the origin, and the x-axis.

Figure 5 .

 5 Figure 5. Computation of Rns and Am based on recursive computation of complex exponential functions.

Figure 6 .

 6 Figure 6. Computation flows for GPCET kernels starting from the order (0, 0). At each step, one of the orders n and m changes its value by adding or -1, depending on the current direction of the flow. Computing each kernel Vnms thus requires two multiplications, one for computing Rns or Am using Equations (7) or (8), depending on whether the value of n or m has just been changed, and the other for multiplying Rns by Am to get Vnms.

Figure 7 .

 7 Figure 7. Computation of GPCET kernels Vnms from the pre-computed and stored values of the radial kernels Rns and angular kernels Am (0 ≤ n, m ≤ K).

Figure 8 .

 8 Figure 8. Computation of GRHFM, GPCT, and GPST radial kernels based on recursive computation of cosine and sine functions.

  Figure 9.Kernel computation times of comparison methods by direct computation at different values of K. For each method and at a specific value of K, the kernel orders used in the computation should satisfy the conditions described in TableIIIand the average elapsed time over all feasible orders is used as the kernel computation time.

Figure 10 .

 10 Figure 10. ZM and GPCET kernel computation times at different values of K using the q-recursive and three recursive strategies proposed in this paper.

Figure 11 .

 11 Figure 11. GPCET moment computation times at different values of K without and with geometrical symmetry.

Figure 12 .

 12 Figure 12. GPCET kernel computation times at different values of n by using recursive and addition chain-based computation.

Table I CARTESIAN

 I AND POLAR COORDINATES OF THE SYMMETRICAL POINTS OF APOINT P 1 EXPRESSED VIA THOSE OF P 1 .

	Symmetrical	Symmetrical	Cartesian	Radial	Angular
	point	axis	coordinates	coordinate	coordinate
	P 1		(x, y)	r	θ
	P 2	y-axis	(-x, y)	r	π -θ
	P 3	origin	(-x, -y)	r	π + θ
	P 4

Table II COMPUTATIONAL

 II COMPLEXITY COMPARISON ON RECURSIVE COMPUTATION OF THE ANGULAR AND RADIAL KERNELS OF PCET.

		Addition Multiplication Table lookup
	[25]	3	4	2
	Proposed method	0	1	1

Table III THE

 III CONSTRAINTS ON THE MOMENT ORDERS (n, m) OF COMPARISON METHODS FOR A FIXED VALUE OF K IN THE EXPERIMENTS ON COMPUTATIONAL COMPLEXITY.

Salvatore Tabbone is a Professor of Computer Science at Université de Lorraine, Nancy, France. He received the Ph.D. degree in computer science from the Institut Polytechnique de Lorraine, Nancy, in 1994. Since 2007, he has been the Head of the QGAR team at the LORIA Laboratory, Université de Lorraine, Nancy, France. His research is related to image and graphics document processing and indexing, pattern recognition, image filtering and segmentation, and content-based image retrieval. He has been the leader of several national and international projects funded by the French and European institutes. He has authored and co-authored more than 100 articles in refereed journal and conferences. He serves as a PC Member for several international and national conferences.