Image geometrical moments are region-based integral features useful for describing the image. The moment invariants are combinations of the geometrical moments. They are invariant to scaling, rotation, and other distortions of the image. The orthogonal moments use orthogonal polynomials as the kernel functions in the moment integral. The image moments, moment variants, and orthogonal moments are basic image descriptors and extensively employed for image analysis, shape description, image indexing, annotation, retrieval, registration, matching, object detection, biometric recognition, high level vehicle navigation, image authentication, and image watermarking. Compared with other image descriptors, the moment-based descriptors have rigorous mathematical bases and show the best performance among the low dimensional descriptors for two-dimensional (2D) images^[1].

The first moment invariants have been introduced by Hu^[2] more than half century ago. After a long period of quiet development in the image moment theories, the research on the moment-based methods and their applications found a fast growing period since the 1990s with thousands of journal papers and book chapters published each year in the last ten years. A number of new invariant moments and orthogonal moments have been introduced, showing better performance for image description. In the same time, an important part of the research effort has been devoted to improving the speed and accuracy of the moment feature calculation.

In this Review, we give a short historic overview on the development of the image moment theories involving the moment invariants, the orthogonal moments, and the discrete orthogonal moments. It is impossible to cover the whole area of the image moments in a short overview. Moreover, several excellent textbooks on the image moments theory and applications have been recently published^[3–6]. We describe several selected representative moments, their performance, and the ideas behind their development. We analyze the information suppression drawback in the Zernike moments, which was not addressed in a recent and more complete overview^[7]. We present a moment family table. We also propose new criteria for assessing performance among the numerous orthogonal polynomials. This Review, however, does not cover issues of the computation speed and numerical stability of the image moments nor their applications.

The geometrical moments of a 2D image $f(x,y)$ in the rectangular coordinate system are defined as $$m_{p,q}=\iint x^p{y}^{q}f(x,y)\mathrm{d}x\mathrm{d}y,$$(1)where the moment orders $p$ and $q$ are positive integers. Thus, the zero order moment represents the image’s total energy. The first order moments specify the location of the image centroid. The moments computed in the coordinate system with the origin placed at the image centroid are the central moments, which are then translation invariant. The second order moments describe the rectangularity and ellipticity of the image and are used to determine the principle axes of the image^[2,8]. The low-order moment features are directly useful for the shape description in the object detection tasks^[9].

The moment invariants first introduced by Hu are mostly known as a set of seven parameters, which are invariant to rotation, scaling, and translation of the image^[2]. Hu’s moment invariants are algebraic combinations of the geometrical moments of up to second and third orders obtained by the algebraic invariants theory as $$\begin{gathered} h_1={\mu }_{20}+{\mu }_{02}, \\ {h}_2={({\mu }_{20}-{\mu }_{02})}^2+4{\mu }_{11}^2, \\ {h}_3={({\mu }_{30}-3{\mu }_{12})}^2+{(3{\mu }_{21}-{\mu }_{03})}^2, \\ {h}_4={({\mu }_{30}+{\mu }_{12})}^2+{({\mu }_{21}+{\mu }_{03})}^2, \\ {h}_5=({\mu }_{30}-3{\mu }_{12})({\mu }_{30}+{\mu }_{12})[{({\mu }_{30}+{\mu }_{12})}^2-3{({\mu }_{21}+{\mu }_{03})}^2] \\ +(3{\mu }_{21}-{\mu }_{03})({\mu }_{21}+{\mu }_{03})[3{({\mu }_{30}+{\mu }_{12})}^2-{({\mu }_{21}+{\mu }_{03})}^2], \\ {h}_6=({\mu }_{20}-{\mu }_{02})[{({\mu }_{30}+{\mu }_{12})}^2-{({\mu }_{21}+{\mu }_{03})}^2] \\ +4{\mu }_{11}({\mu }_{30}+{\mu }_{12})({\mu }_{21}+{\mu }_{03}), \\ {h}_7=(3{\mu }_{21}-{\mu }_{03})({\mu }_{30}+{\mu }_{12})[{({\mu }_{30}+{\mu }_{12})}^2-3{({\mu }_{21}+{\mu }_{03})}^2] \\ -({\mu }_{30}-3{\mu }_{12})({\mu }_{21}+{\mu }_{03})[3{({\mu }_{30}+{\mu }_{12})}^2-{({\mu }_{21}+{\mu }_{03})}^2], \end{gathered}$$(2)where ${\mu }_{p,q}$ is the central moments. The affine transform invariant moments were also later derived from the algebraic invariant theory^[10]. Hu’s moment invariants are highly concentrated image features and are extensively employed for a long period of time. However, Hu’s moment invariants composed from the low-order geometric moments do not describe the details of the image and suffer from information redundancy and information suppression drawbacks, as first pointed out in Ref. [11]. Hu’s moment invariants are information redundant because the kernels $x^p{y}^q$ in their component moments are not orthogonal. As a consequence, an assessment of the quality of the image description by Hu’s moment invariants is not straightforward. One does not know how many Hu’s moment invariants should be used and how well they describe the image for a given task. The information redundancy can be removed by using the orthogonal moments, as will be discussed in Section 3.

The information suppression and loss in Hu’s moment invariants may be shown in evidence by the mathematical equivalence between Hu’s moment invariants and the Fourier-Mellin moments. In the polar coordinate system, the Fourier-Mellin moments of an image $f(r,\theta )$ are defined as $$M_{n,m}={\int }_0^{2\pi }{\int }_0^{1}f(r,\theta )r^{n-1}{e}^{-jm\theta }r\mathrm{d}r\mathrm{d}\theta ,$$(3)where the radial coordinate $r$ is normalized for the scale invariance, the positive integer $n\ge 1$ is the radial moment order, and the integer $m$ is the circular harmonic order. A rotation of the image corresponds to a translation of the image along the angular coordinate and a phase shift of the circular Fourier transform, whose intensity is rotation invariant.

In fact, Hu’s moment invariants can be derived more easily in the polar coordinate system as the angular moments^[12] or the Fourier-Mellin moments^[13,14]. It can be shown^[13] that Hu’s moment invariants are equivalent to a special set of the Fourier-Mellin moments with the circular harmonic orders $m=0–3$ and the radial moment order $n=2$, 3. It is well known that the low circular harmonic orders $m=0–3$ do not describe angular variations in detail. In general, the circular harmonic orders $m=\mathrm{0,1},\dots ,M$, up to $M>10$ are required to reconstruct an image from its circular harmonic functions with acceptable detail^[15]. On the other hand, the high-order monomials $r^n$ suppresses contribution to the moments from the central part of the image and heavily weights the contribution of the background outside the object in the image. The information suppression drawback of the moment invariants can be overcame by using the Fourier-Mellin moments of low radial moment orders.

The complex moments are defined in the Cartesian coordinate system but with the complex coordinate kernels ($x+jy$) and ($x-jy$) instead of $x$ and $y$, as in Eq. (1)^[11,16]. The Fourier-Mellin moments in the polar coordinate system are equivalent to the complex moments in the rectangular coordinate system but with some fractional half-integer orders as $$\begin{gathered} M_{n,m}={\int }_0^{2\pi }{\int }_0^{\infty }f(r,\theta )r^{n-1}{e}^{-jm\theta }r\mathrm{d}r\mathrm{d}\theta \\ ={\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }f(x,y){(x+jy)}^{(n-m)/2}{(x-jy)}^{(n+m)/2}\mathrm{d}x\mathrm{d}y, \end{gathered}$$(4)so that the Fouirer-Mellin moments can be computed in the rectangular coordinate system via the complex moments without transforming the image from the Cartesian coordinate system to the polar coordinate system. The complex moments were also used to derive the affine invariant moments^[4].

The wavelet moments are defined in the polar coordinate system as^[17,18]$$F_{a,b,m}={\int }_0^{2\pi }{\int }_0^1{\psi }_{a,b}(r)e^{-jm\theta }f(r,\theta )r\mathrm{d}r\mathrm{d}\theta ,$$(5)where the circular Fourier transform gives the rotation invariance, and the wavelet transform is applied in the radial coordinate with the mother wavelets such as the Gabor wavelet, Mexican-hat wavelet, or the cubic B-spline function, the latter is shown as $$\psi (r)=\frac{4{\alpha }^{n+1}}{\sqrt{2\pi (n+1)}}\sigma \cos (2\pi f_0(2r-1))\exp (-\frac{{(2r-1)}^2}{2{\sigma }^2(n+1)}),$$(6)where $n=3$, $\alpha =0.697$, $f_0=0.409$, and ${\sigma }^2=0.561$. The wavelet transform is a convolution of the image’s circular harmonic function $$f_m(r)={\int }_0^{2\pi }{e}^{-jm\theta }f(r,\theta )\mathrm{d}\theta ,$$(7)with the wavelets, which are scaled, and translated mother wavelets that are shown as $${\psi }_{a,b}(r)=\frac{1}{\sqrt{a}}\psi \left(\frac{r-b}{a}\right),$$(8)where the discrete scaling factor is $a={(1/2)}^p$, the discrete translation factor is $b=(q/2){(1/2)}^p$, and $p$ and $q$ are positive integers. The continuous wavelets such as the cubic B-spline functions are bandpass filters^[19]. By adaptively choosing the appropriate scale $a$, the wavelet transform can highlight the local features such as edges in the image. However, the wavelet moments as defined in Eq. (5) are not orthogonal and, therefore, can be used as image features for the image description but not for the image reconstruction. Moreover, the wavelet moments are not scale invariant. The image must be resized to a fixed size before computing the wavelet moments^[17,18].

The orthogonal moments are an expansion of the image into the orthogonal polynomial bases. The orthogonal moments are information independent from each other, allowing the image reconstruction for assessing the quality of the image description.

The orthogonal polynomials result from the power series solutions of the ordinary differential equations with the orthogonality of the polynomials related to the Sturm-Liouville forms of the differential equations plus the boundary conditions. They are referred to as the special functions and are given special names because of their frequent uses. On the other hand, the orthogonal polynomials can be easily generated using the Gram-Schmidt orthonormalization, so that there can be an infinite number of orthogonal polynomials. The well-known orthogonal polynomials include the Jacobi, Hermite, Laguerre, and Bessel polynomials. The interrelations between the polynomials exist. For instance, the Gegenbauer, Legendre, Chebyshev, and Zernike polynomials are special cases of the Jacobi polynomials, which are the solutions of the hypergeometric differential equations.

Many orthogonal polynomials are proposed as the bases for the image analysis. For the 2D expansion of an image in the polar coordinate system, the radial orthogonal polynomials are the radial basis, and the Fourier kernels $\exp (jm\theta )$ are used for the circular harmonic analysis. As pointed out by Bhatia and Wolf^[20,21], for the rotation invariance in the form, the circular Fourier kernel is a unique solution to be included as the basis function. The scale invariance is achieved by using the normalized radial coordinate.

The Zernike moments are orthogonal on the unit circle. The Zernike radial polynomials play a dominant role for the aberration characterization in the optical design. The Zernike moments were proposed in 1980 for image analysis^[22] as $$A_{nm}=\frac{n+1}{\pi }{\int }_0^{2\pi }{\int }_0^1{R}_{nm}(r)\exp (-jm\theta )f(r,\theta )r\mathrm{d}r\mathrm{d}\theta ,$$(9)where the Zernike radial polynomials are $$R_{nm}(r)=\sum _{s=0}^{(n-|m|)/2}{(-1)}^s\frac{(n-s)!}{s!((n+|m|)/2-s)!((n-|m|)/2-s)!}{r}^{n-2s},$$(10)where $m$ is the circular harmonic order, and $n$ is the degree of the Zernike polynomial. The Zernike moments are designed to be presented in the Cartesian coordinate system, so that the Zernike radial polynomials must not contain powers of $r$ smaller than the circular harmonic order, $n\ge |m|$ and $n-|m|$ are even^[20]. In the Zernike moments, each term of the radial power in $R_{nm}(r)$ corresponds to a radial moment of a specific order, so that the Zernike moments can be computed as the algebraic combinations of the Fourier-Mellin moments.

Bhatia and Wolf derived^[20] another orthogonal set of polynomials in $x$ and $y$, which contains only the powers of $r$ higher than circular harmonic order $n\ge |m|$, but with the condition for even $n-|m|$ removed^[20]. This moment set is now referred to as the pseudo-Zernike moments^[23]. When high circular harmonics of the order $m$ are used to represent the image of the angular variation in the range $[\mathrm{0,2}\pi ]$^[15], the Zernike and pseudo-Zernike moments with the radial polynomial degrees $n\ge |m|$ must suffer from severe information suppression drawback by the high radial moment orders.

The orthogonal Fourier-Mellin moments^[24] are designed to avoid the information suppression issue in the Zernike moments. The orthogonal Fourier-Mellin moments are defined as $${\mathrm{\Phi }}_{nm}=\frac{1}{2\pi a_n}{\int }_0^{2\pi }{\int }_0^1{Q}_n(r)\exp (-jm\theta )f(r,\theta )r\mathrm{d}r\mathrm{d}\theta ,$$(11)where $a_n=½(n+1)$, and the radial polynomials $$Q_n(r)=\sum _{s=0}^n{\alpha }_{ns}{r}^s,\mathrm{with}{\alpha }_{ns}={(-1)}^{n+s}\frac{(n+s+1)!}{(n-s)!s!(s+1)!}.$$(12)are obtained by the Gram-Schmidt orthogonalization of the power series $r^n$ with $n=0,1,2,\dots$. The orthogonal Fourier-Mellin moments do not need to be expressed in the Cartesian coordinate system. Hence, the restriction of $n\ge |m|$ in the Zernike and pseudo-Zernike moments is removed. This removal of the restriction improves the performance in the image analysis.

The Fourier-Mellin polynomials $Q_n(r)$ have $n$ zeros in $0\le r\le 1$, while the Zernike polynomials have $(n-m)/2$ duplicated roots for $R_{nm}(r)=0$. To have the same number $n$ of zeros, the degree of the Zernike polynomial must be as high as $2n+|m|$. Figure 1 shows the Zernike polynomials $R_{nm}(r)$ for the circular harmonic order $m=10$ and the degree $n=10–28$ with the zeros distributed in the range of radial distance $r$$0.2\le r\le 1$. The closest endpoint to the origin of $R_{\mathrm{28,10}}(r)$ is at $r=0.2$. The image information in the region of the normalized radial distance $r\le 0.2$ is then completely suppressed in this set of Zernike moments. The set of orthogonal Fourier-Mellin radial polynomials $Q_n(r)$ with the degree $n=0–9$, independent of $m$, contain oscillations with the zeros distributed over the radial coordinate range of $0.04\le r\le 1$, as shown in Fig. 1(b). Thus, the set of orthogonal Fourier-Mellin moments favors the analysis of the image, especially for the central part of the image or the small scale object in the unit circle, which occurs in the context of the scale invariant pattern recognition. Figure 2 shows a capital letter E reconstructed from a maximum of 64 orthogonal Fourier-Mellin moments and those from the Zernike moments. The former in Fig. 2(c) shows much better detail than the latter in Fig. 2(b).

The key idea behind the design of the orthogonal Fourier-Mellin moments is to completely separate the radial polynomials from the circular Fourier kernel, and to avoid the constraint that the radial moment orders should be higher than the circular harmonic orders. Based on this idea, many new orthogonal radial polynomials can be proposed, such as Chebyshev-Fourier, Jacobi-Fourier, and Bessel-Fourier polynomials. All the moments of this type can be computed as algebraic combinations of the Fourier-Mellin moments. The latter can be computed from the equivalent complex moments as shown in Eq. (4) without the transformation of the image from the rectangular to the polar coordinate system.

The shifted Chebyshev polynomial of the second kind, which is orthogonal in the interval [0,1], is combined with the circular Fourier kernel to form the kernel of the Chebyshev-Fourier moments as^[25]$${\varphi }_{nm}=\frac{1}{2\pi }{\int }_0^{2\pi }{\int }_0^1{R}_n(r)\exp (-jm\theta )f(r,\theta )r\mathrm{d}r\mathrm{d}\theta ,$$(13)with $n=\mathrm{0,1},2\dots ,m=0,\pm 1,\pm 2,\dots$ and $$R_n(r)=\sqrt{\frac{8}{\pi }}{[w(r)]}^{1/2}{r}^{-1/2}{U}_n^{*}(r),$$(14)with shifted Chebyshev polynomials of the second kind $$U_n^{*}(r)=\sum _{k=0}^{[n/2]}{(-1)}^k\frac{(n-k)!}{k!(n-2k)!}{[2(2r-1)]}^{n-2k},$$(15)which are orthogonal over the range $0\le r\le 1$ as $${\int }_0^1{U}_n^{*}(r)U_k^{*}(r)w(r)\mathrm{d}r=\frac{\pi }{8}{\delta }_{nk},$$(16)where ${\delta }_{nk}$ is the Kronecker symbol, and $w(r)={(r-r^2)}^{1/2}$ is the weight function. As shown in Fig. 3, $R_n(r)$ possesses $n$ zero points evenly distributed over the interval [0,1]. The oscillation amplitude of the $R_n(r)$ of various degrees $n$ are close to constant in a large portion of the radial interval [0,1], except that when $r$ tends to be zero, $R_n(r)$ tends to be high absolute values, as shown in Fig. 3. This can be a potential drawback. Experiments showed^[25] that the performance for describing an image of the Chebyshev-Fourier moments is similar as that of the orthogonal Fourier-Mellin moments^[25].

The Jacobi-Fourier moments are referred to as the generic orthogonal moments^[26] with the kernel composed of the radial Jacobi polynomial and the angular Fourier complex exponential kernel, and are defined as $${\mathrm{\Phi }}_{nm}=\frac{1}{2\pi }{\int }_0^{2\pi }{\int }_0^1{J}_n(p,q,r)\exp (-jm\theta )f(r,\theta )r\mathrm{d}r\mathrm{d}\theta ,$$(17)where the radial Jacobi polynomials are generated by the orthogonalization of the natural power sequence of the radial coordinate with the weighting function^[20,27] as $$w(p,q,r)={(1-r)}^{p-q}{r}^{q-1},$$(18)where $p-q>-1,q>0$, such as $${\int }_0^1{G}_n(p,q,r)G_k(p,q,r)w(p,q,r)\mathrm{d}r=b_n(p,q){\delta }_{nk},$$(19)and the radial Jacobi polynomials are $$G_n(p,q,r)=\frac{n!\mathrm{\Gamma }(q)}{\mathrm{\Gamma }(n+p)}\sum _{s=0}^n{(-1)}^{n-s}\frac{\mathrm{\Gamma }(n+p+s)}{(n-s)!s!\mathrm{\Gamma }(q+s)}{r}^s,$$(20)which become, after the normalization, $$J_n(p,q,r)=G_n(p,q,r)\sqrt{w(p,q,r)/b_n(p,q,r)}.$$(21)where the normalization constant is $$b_n(p,q)=\frac{n!\mathrm{\Gamma }(n+p-q+1){[\mathrm{\Gamma }(q)]}^2}{(2n+p)\mathrm{\Gamma }(n+p)\mathrm{\Gamma }(n+q)}.$$(22)It can be shown that the Jacobi-Fourier moment corresponds to the Legendre-Fourier moments with $p=q=1$ and the unit weighting function, the Chebyshev-Fourier moments with $p=2$ and $q=3/2$, and the orthogonal Fourier-Mellin moments with $p=q=2$ and the weighting function equal to $r$^[26]. The Zernike radial polynomial corresponds to the Jacobi radial polynomial as $r^m{J}_k(p,q,r^2)$ with $p=q=m+1$, where $m$ is the circular harmonic order in the Zernike moments, and $k=(n-m)/2$. The pseudo-Zernike radial polynomial corresponds to the Jacobi radial polynomial $J_k(p,q,r)$ with $p=2m+2$, $q=m+2$ and $k=n-m$^[3]. Therefore, the Jacobi-Fourier moments provide a generic framework to express the rotation-invariant orthogonal moments. The common formulation is a benefit for the performance comparison and analysis of the orthogonal moments and for searching prime orthogonal moments.

The Bessel-Fourier moments were introduced recently, which are defined as^[28]$$B_{nm}=\frac{1}{2\pi {\alpha }_n}{\int }_0^{2\pi }{\int }_0^1{J}_v({\lambda }_{n}r)\exp (-jm\theta )f(r,\theta )r\mathrm{d}r\mathrm{d}\theta ,$$(23)where $J_v({\lambda }_{n}r)$ is the Bessel radial polynomials of order $v$ of the first kind, and ${\lambda }_n$ is the $n$th zero of $J_v({\lambda }_n)=0$. The Bessel polynomials of the same order $v$ are orthogonal over the radial interval [0,1], such as $${\int }_0^1{J}_v({\lambda }_{n}r)J_v({\lambda }_{k}r)\mathrm{d}r={\alpha }_n{\delta }_{nk},$$(24)where ${\alpha }_n={[J_{v+1}({\lambda }_n)]}^2/2$ is the normalization constant. Figure 4 is a plot of the Bessel function $v=1$. The Bessel polynomials of the first kind $J_1({\lambda }_{n}r)$ tend to be zero at both ends of the radial interval [0,1] and show a slow variation of the oscillation amplitude over the radial range [0,1], as shown in Fig. 4. Only in the range of $r$ close to $r=0$, the oscillation amplitudes are doubled.

Compared to the orthogonal Fourier-Mellin polynomials and the shifted Chebyshev polynomials, both of which tend to have high values when $r$ tends to be zero, as shown in Figs. 1 and 3, the Bessel-Fourier moments show a better performance for image analysis. They also are better than the Zernike moments, as the Zernike moments show a sharp increase to the unit when $r$ tends to be $r=1$, as shown in Fig. 1(a).

The radial-harmonic-Fourier moments^[29], or polar harmonic transform^[30,31], includes the Fourier harmonic analysis in both the radial and angular coordinates and are defined as $${\varphi }_{nm}=\frac{1}{2\pi }{\int }_0^{2\pi }{\int }_0^1{T}_n(r)\exp (-jm\theta )f(r,\theta )r\mathrm{d}r\mathrm{d}\theta .$$(25)where the radial kernels $T_n(r)$ are the orthogonal triangular functions, such as $T_n(r)=\sin [(n+1)\pi r$] for the odd $n$, and $T_n(r)=\cos (n\pi r)$ for the even $n$, or the complex exponential functions such as $T_n(r)=\exp (j2\pi n{r}^2)$ or $T_n(r)=\exp (jnr)$. In all cases, the $T_n(r)$ should be normalized, respectively. The corresponding transforms are referred to as the radial-harmonic-Fourier moments^[29], the polar harmonic transform^[30], and the exponent-Fourier moments^[31]. The angular radial transform, the polar complex exponential transform, the polar cosine transform, and the polar sine transform also belong to this group^[31–33]. As the Fourier kernels have uniform distributions of the zeros and a constant oscillation amplitude over the interval $0\le r\le 1$ and $0\le \theta \le 2\pi$, this type of transform shows a superior performance over the other rotation-invariant orthogonal moments with radial polynomial kernels^[34].

The radial-harmonic-Fourier moments are not scale invariant. An image must be rescaled to a given size before computing the moments. Moreover, the radial-harmonic-Fourier moments should preferably be referred to as the radial-harmonic-Fourier transforms, as the moments are essentially formulated with the monomial and polynomial kernels. Moreover, the moments with the radial polynomials can achieve the scale invariance by using the normalized radial coordinate.

There are an infinite number of complete sets of orthogonal polynomials. Many of them were proposed as the kernels of the orthogonal moments for 2D images in the Cartesian coordinate system. The orthogonal moments in this category include the Legendre moments^[22,35], the Gegenbauer moments^[36,37], the Gaussian-Hermite moments^[38–40], etc. Defined in the Cartesian coordinate system, these moments enjoy the orthogonality of the orthogonal polynomials in $x$ and $y$ in the interval $[-\mathrm{1,1}]$ or [0,1] for image analysis and reconstruction. They can also be scale invariant when using the scaling factor to normalize the coordinates.

A new type of the orthogonal moments uses two different polynomials as the basis for each dimension in the Cartesian coordinate system separately, such as the Chebyshev-Gegenbauer and the Gegenbauer-Legendre moments. These moments are called separable moments^[41,42].

The orthogonal moments in the Cartesian coordinate system are not rotational invariant by definition. However, the rotation invariance can still be introduced by the transformation between the Cartesian and polar coordinate systems.

Performances of the orthogonal moment image descriptors have been assessed in many Letters via the image reconstruction errors and the pattern recognition scores. However, these merits are dependent on the selected moment orders and the testing image, and therefore are the indirect measures. In fact, as the orthogonal moments are the expansions of the image into the orthogonal polynomial bases, the performance of the image description by the orthogonal moments is determined directly by the bases’ functions. In the case of rotation-invariant orthogonal moments, the Fourier harmonic analysis is applied along the angular coordinate. Hence, the performance of the image analysis along the radial coordinate depends uniquely on the set of the radial orthogonal polynomials. For the image analysis, the radial polynomial set should have a high number of real-valued zeros, a wide range of the distribution of zeros in the whole range $0\le r\le 1$, and a constant oscillation amplitude. From instance, the Zernike radial polynomials are restricted to having high degrees, and one set of Zernike moments, shown in Fig. 1(a), have the zeros distributed in the range of $0.2\le r\le 1$. The orthogonal Fourier-Mellin polynomials and shifted Chebyshev polynomials both tend to have high values when $r$ tends to be zero, as shown in Figs. 1(b) and 3. Those facts represent the weakness of these moments. The Bessel polynomials have $n+2$ zeros distributed in the range of $0\le r\le 1$, and a slow variation of the oscillation amplitudes, as shown in Fig. 4, so that they can have a better performance for the image analysis. In this sense, the radial-harmonic-Fourier transforms have the best performance, as this is simply the Fourier harmonic analysis in both the angular and radial dimensions. However, they are not image moments in the mathematical sense and are not scale invariant.

Furthermore, in terms of the space-frequency space analysis, the Fourier kernels, which are perfectly local in frequency but are global in space, are not efficient for describing local details in the image. The Fourier moments miss the location of an image detail, and very high frequency orders must be involved to reconstruct a sharp detail, such as an edge in the image. In this sense, the wavelet moment using the Gabor transform and the wavelet transform can be superior over the Fourier transform^[43]. However, the wavelet moments using the continuous mother wavelets are not orthogonal.

Table 1 shows a family of the moments with the classification. However, it is regrettable to not be able to list all image moments in Table 1 in the period of time when many novel moments are being proposed. For a given image description task, the choice among the available moments should depend on the specific tasks, the requested invariance, the computation performance, and many other considerations.

Most orthogonal moments discussed previously belong to continuous moments, as their bases are continuous functions with their orthogonality defined by the continuous integrals. Numerical computation of the continuous moments requires discretization of the continuous polynomials and discrete approximation of the continuous moment integrals, which causes the errors, affecting the accuracy of the calculation and the quality of image reconstruction with increased computational complexity.

The discrete orthogonal moments^[44–55] use discrete orthogonal polynomials as the basis set, which satisfies the orthogonality exactly without numerical approximation. For a smooth image intensity distribution, whose value is known only at a set of discrete coordinate points, the discrete polynomials $p_n(x)$ are also defined on these distinct real nodes, and satisfy the discrete orthogonality relation as $$\sum _{x=0}^{N-1}{p}_m(x)p_n(x)w(x)=\rho (n){\delta }_{mn},$$(26)where $0\le m,n\le N-1$, and $\rho (n)$ is the square norm of $p_n(x)$. The weight $w(x)$ in the orthogonality expression (26) of the discrete polynomials can be the Meixner, Charlier, Krawtchouk, and Hahn weight. Thus, the discrete orthogonal polynomials are uniquely determined by the given distinct nodes support, the weights, and the orthogonality condition, and can be built from monomials by the Gram-Schmidt process.

The discrete orthogonal moments of a 2D image $f(x,y)$ in the Cartesian coordinate system on a closed interval are defined as $$P_{mn}=\sum _{x=0}^{M-1}\sum _{y=0}^{N-1}f(x,y){\tilde{p}}_m(x){\tilde{p}}_n(y),$$(27)where the normalized discrete polynomial bases ${\tilde{p}}_n(x)$ are shown as $${\tilde{p}}_n(x)=p_n(x)\sqrt{w(x)/\rho (n)}.$$(28)As the discrete polynomials satisfy the orthogonal property precisely, no approximation errors are involved in the computation of moments. The image can be accurately reconstructed from the complete set of discrete moments. When two different discrete orthogonal polynomials serve as the bases in $x$ and $y$, separately, the moments are separable moments such as that using the Tchebichef-Krawtchouk, Hahn-Krawtchouk, and Charlier-Hahn polynomials combined with Tchebichef, Krawtchouk, Hahn, or Charlier polynomials^[41,42].

The first discrete orthogonal Tchebichef (Chebyshev) moments^[44] were proposed based on the discrete orthogonal Tchebichef polynomials. The discrete Tchebichef polynomials are shown as $$t_n(x)=n!\sum _{k=0}^n{(-1)}^{n-k}\left(\begin{array}{l}N-1-k\\ n-k\end{array}\right)\left(\begin{array}{l}n+k\\ n\end{array}\right)\left(\begin{array}{l}x\\ k\end{array}\right),$$(29)satisfying the orthogonality condition with the unit weight and the normalization constants $$\begin{gathered} \rho (n,N)=\frac{N(N^2-1)(N^2-2^2)\cdots (N^2-n^2)}{2n+1} \\ =(2n)!\left(\begin{array}{c}N+n\\ 2n+1\end{array}\right),n=\mathrm{0,1},\dots N-1. \end{gathered}$$(30)The normalized Tchebichef polynomials are $${\tilde{t}}_n(x)=t_n(x)/\beta (n,N),$$(31)with the normalized constant $$\tilde{\rho }(n,N)=\rho (n,N)/\beta {(n,N)}^2\mathrm{and}\beta (n,N)=N^n.$$(32)The discrete Tchebichef moments exhibit a superior performance for image reconstruction and feature representation compared to the continuous Legendre and Zernike moments^[44].

One of representative discrete orthogonal moments defined in the Cartesian coordinate system is the Krawtchouk^[46] moment, which are based on the Krawtchouk polynomials. The kernel functions are the separable discrete orthogonal polynomials in $x$ and $y$. The discrete Krawtchouk polynomials satisfy the orthogonality conditions with the Krawtchouk weight as $$w(x;p,N)=\left(\begin{array}{l}N\\ x\end{array}\right)p^x{(1-p)}^{N-x},$$(33)and the normalization constant $$\rho (n;p,N)={(-1)}^n{\left(\frac{1-p}{p}\right)}^n\frac{n!}{{(-N)}_n}.$$(34)The Krawtchouk polynomials are obtained as $$K_n(x;p,N)=\sum _{k=0}^N{a}_{k,n,p}{x}^k={}_{2}F_1(-n,-x;-N;\frac{1}{p}),$$(35)where $x,n=\mathrm{0,1},2\dots ,N$, $N>0$, and $p\in (\mathrm{0,1})$, the ${{}_{r}F}_s$ is the hypergeometric series, defined as $${{}_{r}F}_s(a_1,\cdots ,a_r;b_1,\cdots ,b_s;z)=\sum _{k=0}^{\infty }\frac{{(a_1)}_k{(a_2)}_k\cdots {(a_r)}_k}{{(b_1)}_k{(b_1)}_k\cdots {(b_s)}_k}\frac{z^k}{k!},$$(36)and ${(a)}_k$ is the Pochhammer symbol given by $${(a)}_k=a(a+1)\dots (a+k-1)=\frac{\mathrm{\Gamma }(a+k)}{\mathrm{\Gamma }(a)}.$$(37)The Krawtchouk moments perform better than the Tchebichef moments in terms of the reconstruction error. They were employed to extract local features of an image by varying the parameter of the binomial distribution associated with the Krawtchouk polynomials^[46].

Two types of Hahn polynomials with different standardizations were used to define two types of Hahn moments^[47,48], respectively. The Hahn polynomials are expressed as $$\begin{gathered} h_n^{(\mu ,\nu )}(x,N)={(N+\nu -1)}_n{(N-1)}_n \\ \times \sum _{k=0}^n{(-1)}^k\frac{{(-n)}_k{(-x)}_k{(2N+\mu +\nu -n-1)}_k}{{(N+\nu -1)}_k{(N-1)}_k}\frac{1}{k!}, \end{gathered}$$(38)where $\mu ,\nu (\mu >-1,\nu >-1)$ are the adjustable parameters controlling the shape of the polynomials. The Hahn polynomials satisfy the discrete orthogonal condition as $$\sum _{x=0}^{N-1}w(x)h_m^{(\mu ,\nu )}(x,N)h_n^{(\mu ,\nu )}(x,N)=d_n^2{\delta }_{mn},$$(39)with the weighting function as $$w(x)=\frac{1}{\mathrm{\Gamma }(x+1)\mathrm{\Gamma }(x+\mu +1)\mathrm{\Gamma }(N+\nu -x)\mathrm{\Gamma }(N-n-x)},$$(40)and the square norm $d_n^2$ has the following expression $$\begin{gathered} d_n^2=\frac{\mathrm{\Gamma }(2N+\mu +\nu -n)}{(2N+\mu +\nu -2n-1)\mathrm{\Gamma }(N+\mu +\nu -n)} \\ \times \frac{1}{\mathrm{\Gamma }(N+\mu -n)\mathrm{\Gamma }(N+\nu -n)\mathrm{\Gamma }(n+1)\mathrm{\Gamma }(N-n)}. \end{gathered}$$(41)Compared with the Tchebichef moments and the Krawtchouk moments, the first type of Hahn moments has better image reconstruction results^[47]. Another type of Hahn moments^[48] is based on another form of Hahn polynomials, shown as $$h_n(x;\alpha ,\beta ,N)={}_{3}F_2(-n,n+\alpha +\beta +1,-x;\alpha +1,-N;1),$$(42)where $x,n=\mathrm{0,1},2,\dots ,N-1$, $\alpha >-1,\beta >-1$, or $\alpha <-N$, $\beta <-N$, with the weighting function $$w(x;\alpha ,\beta ,N)=\left(\begin{array}{c}\alpha +x\\ x\end{array}\right)\left(\begin{array}{c}\beta +N-x\\ N-x\end{array}\right),$$(43)and the square norm $$d^2(n;\alpha ,\beta ,N)=\frac{{(-1)}^n{(n+\alpha +\beta +1)}_{N+1}{(\beta +1)}_{n}n!}{(2n+\alpha +\beta +1){(\alpha +1)}_n{(-N)}_{n}N!}.$$(44)This type of Hahn moments provides a unified understanding of the Tchebichef and Krawtchouk moments. The two latter moments can be obtained as particular cases of the Hahn moments with the appropriate parameter settings. This fact implies that the Hahn moments encompass all their properties and exhibit intermediate properties between the extremes set by the Tchebichef and Krawtchouk moments. Besides, the Hahn moments, as a generalization of the Tchebichef and Krawtchouk moments, can be used for global and local feature extraction^[48].

The discrete Tchebichef, Krawtchouk, and Hahn polynomials are orthogonal on a uniform lattice $\{x=\mathrm{0,1},2,\dots \}$. There is another category of discrete variable polynomials, which are orthogonal on a non-uniform lattice $\{x=x(s),s=\mathrm{0,1},2,\dots \}$. For different non-uniform lattice functions, there are different polynomials. The dual Hahn polynomials and the Racah polynomials are two types of discrete polynomials orthogonal on a non-uniform lattice $x(s)=s(s+1)$. Due to the finite domain of the definition, the dual Hahn polynomials and Racah polynomials were explored to construct the corresponding dual Hahn moments^[50] and Racah moments^[52]. They are defined in the Cartesian coordinate system, and their kernel functions are the product of two corresponding two corresponding 1D separable discrete orthogonal polynomials.

Given a uniform pixel lattice image $f(s,t)$ with a size $N\times N$, the dual Hahn moments^[50] are defined as $$W_{nm}=\sum _{s=a}^{b-1}\sum _{t=a}^{b-1}{\widehat{w}}_n^{(c)}(s,a,b){\widehat{w}}_m^{(c)}(t,a,b)f(s,t),$$(45)with $n,m=\mathrm{0,1},2,\dots ,N-1$, and the dual Hahn polynomials $w_n^{(c)}(s,a,b)$ defined as $$\begin{gathered} w_n^{(c)}(s,a,b)=\frac{{(a-b+1)}_n{(a+c+1)}_n}{n!} \\ \times {}_{3}F_2(-n,a-s,a+s+1;a-b+1,a+c+1;1), \end{gathered}$$(46)where parameters $a$, $b$, and $c$ are restricted to $-(1/2)<a<b$, $|c|<1+a,b=a+N$. The dual Hahn polynomials satisfy the following orthogonality property $$\sum _{s=a}^{b-1}\rho (s)w_n^{(c)}(s,a,b)w_m^{(c)}(s,a,b)\left[\mathrm{\Delta }x(s-\frac{1}{2})\right]=d_n^2{\delta }_{nm}.$$(47)where $\rho (s)$ is the weighting function $$\rho (s)=\frac{\mathrm{\Gamma }(a+s+1)\mathrm{\Gamma }(c+s+1)}{\mathrm{\Gamma }(s-a+1)\mathrm{\Gamma }(b-s)\mathrm{\Gamma }(b+s+1)\mathrm{\Gamma }(s-c+1)},$$(48)and $d_n^2$ is the square norm $$d_n^2=\frac{\mathrm{\Gamma }(a+c+n+1)}{n!(b-a-n-1)!\mathrm{\Gamma }(b-c-n)}.$$(49)The normalized dual Hahn polynomials are $${\widehat{w}}_n^{(c)}(s,a,b)=w_n^{(c)}(s,a,b)\sqrt{\frac{\rho (s)}{d_n^2}\left[\mathrm{\Delta }x(s-\frac{1}{2})\right]}.$$(50)The Tchebichef polynomials, Krawtchouk polynomials, and Hahn polynomials are special cases of the dual Hahn polynomials. Thus, the dual Hahn moments which contain more parameters have a generality and give more flexibility to promote the image describing ability. The dual Hahn moments perform better than the Legendre moments, Tchebichef moments, and Krawtchouk moments in terms of image reconstruction capability^[50].

Similar to the dual Hahn moments, for the image $f(s,t)$ with a size $N\times N$, the Racah moments^[52] are defined as $$U_{nm}=\sum _{s=a}^{b-1}\sum _{t=a}^{b-1}{\widehat{u}}_n^{(\alpha ,\beta )}(s,a,b){\widehat{u}}_m^{(\alpha ,\beta )}(t,a,b)f(s,t),$$(51)with $n,m=\mathrm{0,1},2,\dots ,N-1$. Racah polynomials $u_n^{(\alpha ,\beta )}(s,a,b)$ are defined as $$\begin{gathered} u_n^{(\alpha ,\beta )}(s,a,b)=\frac{1}{n!}{(a-b+1)}_n{(\beta +1)}_n{(a+b+\alpha +1)}_n \\ \times {}_{4}F_3(-n,\alpha +\beta +n+1,a-s,a+s+1;\beta +1,a-b+1,a+b+\alpha +1;1), \end{gathered}$$(52)where the parameters $a,b,\alpha$ and $\beta$ are restricted by $$-(1/2)<a<b,\alpha >-1,-1<\beta <2a+1,b=a+N.$$The Racah polynomials satisfy the following orthogonality property $$\sum _{s=a}^{b-1}\rho (s)u_n^{(\alpha ,\beta )}(s,a,b)u_m^{(\alpha ,\beta )}(s,a,b)\left[\mathrm{\Delta }x(s-\frac{1}{2})\right]=d_n^2{\delta }_{nm},$$(53)with the weight $$\rho (s)=\frac{\mathrm{\Gamma }(a+s+1)\mathrm{\Gamma }(s-a+\beta +1)\mathrm{\Gamma }(b+\alpha -s)\mathrm{\Gamma }(b+\alpha +s+1)}{\mathrm{\Gamma }(a-\beta +s+1)\mathrm{\Gamma }(s-a+1)\mathrm{\Gamma }(b-s)\mathrm{\Gamma }(b+s+1)},$$(54)and the square norm $$\begin{gathered} d_n^2=\frac{\mathrm{\Gamma }(\alpha +n+1)\mathrm{\Gamma }(\beta +n+1)\mathrm{\Gamma }(b-a+\alpha +\beta +n+1)}{(\alpha +\beta +2n+1)n!(b-a-n-1)!\mathrm{\Gamma }(\alpha +\beta +n+1)} \\ \times \frac{\mathrm{\Gamma }(a+b+\alpha +n+1)}{\mathrm{\Gamma }(a+b-\beta -n)}. \end{gathered}$$(55)The normalized Racah polynomials are $${\widehat{u}}_n^{(\alpha ,\beta )}(s,a,b)=u_n^{(\alpha ,\beta )}(s,a,b)\sqrt{\frac{\rho (s)}{d_n^2}\left[\mathrm{\Delta }x(s-\frac{1}{2})\right]}.$$(56)The Racah moments perform better than the Legendre moments, Tchebichef moments, and Krawtchouk moments for noise-free image reconstruction^[52].

Rotation-invariant discrete radial Tchebichef moments defined in polar coordinate systems were proposed based on the 1D radial discrete Tchebichef polynomials and the circular Fourier kernel^[53]. The primary advantage of the discrete radial moments is the rotational invariance. For an image $f(r,{\theta }_k)$, the explicit expressions of the radial Tchebichef moments^[53] are as follows: $$S_{pq}=\frac{1}{n\rho (p,m)}\sum _{r=0}^{m-1}\sum _{k=0}^{n-1}{t}_{p,m}(r)e^{-j\frac{2\pi qk}{n}}f(r,{\theta }_k),$$(57)where $r=\mathrm{0,1},\dots ,m-1$, $m=N/2,{\theta }_k=2\pi k/n$, and $k=\mathrm{0,1},\dots ,n-1$, with $n$ denoting the maximum number of pixels along the circumference of the circle of radius $m$, $$t_{p,m}(r)=\frac{p!}{m^p}\sum _{k=0}^p{(-1)}^{p-k}\left(\begin{array}{l}m-1-k\\ p-k\end{array}\right)\left(\begin{array}{l}p+k\\ p\end{array}\right)\left(\begin{array}{l}r\\ k\end{array}\right),$$(58)$$\rho (p,m)=\frac{m(1-\frac{1}{m^2})(1-\frac{2^2}{m^2})\cdots (1-\frac{p^2}{m^2})}{2p+1},$$(59)with $p=\mathrm{0,1},\dots m-1$. The rotation-invariant discrete radial Krawtchouk moments^[54] were introduced in a similar way to the radial Tchebichef moments. The kernel functions of the radial Krawtchouk moments are products of radial 1D Krawtchouk polynomials and angular circular Fourier functions.

We have reviewed some major steps in the development of the image moment theory, involving the geometrical moments, the moment invariants, and the orthogonal moments. It is impossible for this short Review to cover all the important moments in detail. Moreover, we have not included in this Review the works on the speed and accuracy of computation of the image moments and on their applications. A large body of excellent Letters and text books exist in the literature now. A more rigorous, accurate, and complete mathematical analysis on the image moments theory has been developed. However, some fundamental intuitive ideas, such as the use of low-order radial moments for avoiding information suppression, the definition of the moment invariants in the polar coordinate system and relation of them to the moments in the rectangular coordinate system via the complex moments, and the separation of the radial moment orders from the circular harmonic orders, have been widely accepted and continuously used in the development of the image moment-based methods. We proposed a moment family table and the performance measures for the orthogonal moments based on the number of zeros, the uniform distribution of zeros, and the oscillation amplitude variation of the orthogonal polynomials. The image reconstruction can depend on the choice of moment orders and on the testing image set so that they can be indirect measures for the performance.