Efficient Computation of Higher Order 2D Image Moments using the Discrete Radon Transform
EEfficient Computation of Higher Order 2D ImageMoments using the Discrete Radon Transform
William Diggin Michael DigginSeptember 2020
Abstract
Geometric moments and moment invariants of image artifacts havemany uses in computer vision applications, e.g. shape classification orobject position and orientation. Higher order moments are of interestto provide additional feature descriptors, to measure kurtosis or to re-solve n-fold symmetry. This paper provides the method and practicalapplication to extend an efficient algorithm, based on the Discrete RadonTransform, to generate moments greater than the 3 rd order. The mathe-matical fundamentals are presented, followed by relevant implementationdetails. Results of scaling the algorithm based on image area and its com-putational comparison with a standard method demonstrate the efficacyof the approach. In his seminal paper, Hu [6] defined 10 geometric moments up to the 3 rd orderwhich are used to derive 7 moment invariants, e.g. invariant to position, scale,orientation, etc. These invariants are useful as feature descriptors for objects inimages. The moment generating equation for a 2-D image I ( i, j ) of order ( p + q )is: M pq = (cid:88) i,j I ( i, j ) i p j q (1)Equivalent 4 th order moments have been utilised by Prokop et al [8] to mea-sure kurtosis. Higher order moments have been shown to be useful discrimina-tors for n-fold rotational symmetry [7]. Some previous work in image momentsfocuses on the algorithms and the computational cost of those algorithms togenerate moments up to the 3 rd order; the goal to optimise execution time withlittle or no loss to accuracy. An algorithm based on the Discrete Radon Trans-form [1] exhibits reduced complexity and faster execution time in comparisonto several other techniques, while retaining accuracy [4]. Other early work inusing the Discrete Radon Transform for image moments can be attributed toGindi et al [5] and Shen et al [9]. 1 a r X i v : . [ c s . C V ] S e p n section 2, we recall the Discrete Radon Transform and its application toimage moments calculation. In section 3, we outline the approach to extendthe DRT algorithm to the 4 th order. Then, some practical aspects are disclosedwhen implementing the algorithm, followed by a presentation of comparativeresults that demonstrate the efficacy of the algorithm relative to the OpenCVLibrary [2]. Finally, in section 6, we conduct some analysis of the results andcapture some of the learnings from the approach. The Radon Transform is the set of projections of a function f ( x, y ) onto lineintegrals at various angles over a range, 0 ≤ θ < π : R ( x (cid:48) , θ ) = (cid:90) f ( x, y ) dy (cid:48) where (cid:18) x (cid:48) y (cid:48) (cid:19) = (cid:18) cos θ sin θ − sin θ cos θ (cid:19) (cid:18) xy (cid:19) For the Discrete Radon Transform (DRT), let R x : y be a projection at angle θ = tan − yx of an image I ( i, j ), whose size is M · N . The quantities x : y arethe ratio of horizontal to vertical units, equivalent to slope. For an image, theseunits are pixels. Given this, we can efficiently compute the projections of theimage at angles 0 ° , 90 ° , 45 ° (diagonal) and 135 ° (anti-diagonal) as follows: R [ k ] = (cid:88) j I ( k, j ) where k = 0 ... ( M − R [ k ] = (cid:88) j I ( j, k ) where k = 0 ... ( N − R [ k ] = (cid:88) i + j = k I ( i, j ) where k = 0 ... ( M + N − R − [ k ] = (cid:88) j − i = k I ( i, j ) where k = ( − N +1) ... ... ( M − R − .Let M x : yr be the moment, of order r , for a 1-D projection of a DRT whosepixel ratio is x : y . Then the equivalent 1-D moments over the projection aregiven by: M x : yr = (cid:88) k R x : y [ k ] · k r (2)It was shown that various 1-D moments over the 4 projections R , R , R and R − , can be used to compute the 10 Hu 2-D moments up to 3 rd order[4]. The stated benefit is a reduction in computation time since the complexityof the DRT algorithm is O ( n ) versus O ( n ) for the 2-D computational approach.2 Extended DRT Algorithm
Extending the DRT algorithm to generate some th order 2-D moments doesnot require additional projections. For example, the moments, M and M ,are trivial to compute from (2) using the horizontal and vertical projections asfollows: M = M (3) M = M (4)The 4 th order M moment can be computed from the 4 th order diagonaland anti-diagonal 1-D moments as follows. Take the diagonal: M = (cid:88) k R [ k ] · k M = (cid:88) k (cid:88) i + j = k I ( i, j ) · k M = (cid:88) i,j I ( i, j )( i + j ) M = (cid:88) i,j I ( i, j )( i + j + 4 i j + 4 ij + 6 i j ) M = M + M + 4 M + 4 M + 6 M (5)Similarly, for the anti-diagonal: M − = (cid:88) k (cid:88) j − i = k I ( i, j ) · k M − = M + M − M − M + 6 M (6)Adding (5) and (6) and rearranging: M = ( M + M − − M − M ) /
12 (7)Computing M and M does require an additional projection. A conve-nient projection is: R [ k ] = (cid:88) i +2 j = k I ( i, j ) where k = 0 ... ( M +2 N −
1) (8)This is equivalent to a projection at a slope of 2 (angle ≈ ° ) onto a lineintegral, as depicted in figure 1. Then, to derive M and M , we use (2) and3igure 1: A Projection of the Radon Transform for Slope=2.(projection at angle ≈ ° )(8) and note: M = (cid:88) k R [ k ] · k M = (cid:88) k (cid:88) i +2 j = k I ( i, j ) · k M = (cid:88) i,j I ( i, j )( i + 2 j ) M = (cid:88) i,j I ( i, j )( i + 16 j + 8 i j + 32 ij + 24 i j ) M = M + 16 M + 8 M + 32 M + 24 M (9)Combining (5) with (9) and rearranging: M = ( M − M + M − M − M ) /
24 (10)By further substitution into (9): M = ( M − M − M − M − M ) / M · N .(for Slope=2, projection length = M +2 N ) In practical terms, the projection at slope 2 can be discrete (see figure 2) re-sulting in a linear projection. We benefit from not requiring interpolation andthe associated cost of that computation. This is at the expense of memory, asthe resulting projection is of length M +2 N .When executing the DRT algorithm, the underlying processor capabilitiesare a vital component in efficient computation. Over the years, the cost ofinteger multiplication has improved to equate to the cost of integer addition –in line with advancement in processor technologies and their architectures [3]. Sotoo has floating-point arithmetic. On the Intel processors, the x64 architecture,SIMD (Single Instruction, Multiple Data) and further instruction sets (SSE –Streaming SIMD Extensions or AVX – Advanced Vector Extensions) provideoptimisations where computations can be parallelised or pipelined, resulting ingreater algorithm performance. In the example presented herein, we encode thealgorithm to rely on and take advantage of the AVX optimisations providedthrough compiler settings (e.g. /arch:AVX2 in the Microsoft C++ compiler).5 Results
To demonstrate the improvement from the DRT algorithm, we compare it toan extended OpenCV implementation – where the base algorithm used in theOpenCV library is extended to compute the 4 th order moments. Using a largesample size and selecting the fastest exemplar, the comparison of executiontimes was gathered for differing image sizes. Both algorithms were executed onan Intel Core i5 7 th Gen CPU. Although we are concerned mainly with relativetimes, execution with and without AVX processor optimisations were also noted.Table 1 and figure 3 show the measurements and plot comparing the timeneeded to compute the raw image moments of a given image size. The smallestimage used was a 200x200 pixel image, with the largest at 4032x3024 pixels.The y-axis illustrates the execution time in µ -seconds (log base 10 for easierplotting). The x-axis is the square root of the number of pixels in the image (asthis is a better method to illustrate the performance in terms of image size).We also include the DRT method with and without SSE2/AVX optimisation(provided through the compiler settings). Image Size OCV4 DRT4 DRT4+AVX2 (pixels) µ sec µ sec µ sec4032 3024 20363 7426 44553000 3000 14644 5444 32602000 2000 6209 2194 12921500 1500 3500 1151 6651000 1000 1552 510 287750 750 878 296 175400 400 253 103 54200 200 64 27 16Table 1: Comparative Computational Time vs Image Size for OpenCV andDRT Algorithms. The na¨ıve implementation to compute all moments up to the 4 th order usingequation (1) requires (15 M · N ) multiplications and (15 M · N ) additions. Frominspection of the OpenCV implementation, the same computation requires (4 M · N +12 N ) multiplications with (5 M · N +15 N ) additions. The DRT algorithmpresented herein needs only (10 M +11 N ) multiplications and (5 M · N +11 M +11 N )additions – overall a significant improvement. A 4x–5.5x performance gainrelative to the OpenCV implementation is achieved across various image sizes.To extend the approach to higher order moments, e.g. 5 th or 6 th order,requires further projections. The standard horizontal and vertical projections6igure 3: Comparative Computational Cost Vs Image Area for OpenCV andDRT Algorithms.will continue to generate the desired M i , M i moments. The diagonal andanti-diagonal projections will contribute to computing the M ij moments, butthe number of unknowns will exceed the number of equations to resolve allmoments, as exemplified for the 4 th order moment computation. In that case,ratios of 2:1, 1:2, -1:2 and 1:-2 are convenient inputs to the transform. Whilethe 4 th order moments have 5 unknowns ( M , M , M , M , M ) to resolveand required 5 equations from 5 1-D moments ( M , M , M , M − and M ), the 5 th order moment set requires 6 equations for 6 unknowns ( M , M , M , M , M , M ); the 6 th order requiring 7 equations for ( M , M , M , M , M , M , M ) and so on.The selection of slope is not arbitrary. Indeed, it is selected to provide fora discrete summation into a linear vector. The desired slopes are derived fromthe integer ratio of vertical to horizontal pixels. For example, 3:1 or 3:2 andtheir corollary, 1:3 or 2:3, provide for convenient summations that do not requireinterpolation nor are so large as to make them unwieldy. Those ratio choicescan be further expanded upon via reflections about the X- and Y-axes.In execution, particularly on Intel CPUs, the cost of multiplication versusaddition, while not significant, does contribute additional execution time. Withoptimisation and careful algorithm construction, we take advantage of pipelin-ing and parallelism to counter that margin. Use of the AVX optimisationscontributes a ≈ . x gain for the DRT algorithm.Relative to the OpenCV algorithm, much of the computational savings arisefrom reducing the total number of operations, either multiplication or addition.The OpenCV algorithm has nearly twice as many operations as the DRT algo-rithm. The contributions from a reduction in the number of total operations,a reduction in the number of multiplications and careful algorithm/processoroptimisations leads to a significant computational advantage.7 Conclusion
The DRT algorithm provides increasing computational advantages as either ofimage size and/or moment order increases. Extending the algorithm to computethe full 4 th order moment set requires an additional projection in addition tothose required to compute all the 3 rd order moments. The resulting computa-tions are simple and elegant. The 4 th order moments of 5 projections generateall 15 2-D moments up to the 4 th order exactly. The DRT algorithm reduces thenet multiplicative complexity from O ( n ) to O ( n ) although the additive com-plexity remains O ( n ). The approach allows a practical means to extend thealgorithm to higher order moments, while retaining computational efficiency.Further work in this area may seek to determine optimal relationships be-tween higher order moments and the Inverse Radon Transform, e.g. an algo-rithm to reconstruct an image from its 1-D moments. References [1] G Beylkin. “Discrete radon transform”. In:
IEEE Transactions on Acous-tics, Speech, and Signal Processing
Dr. Dobb’s Journal of SoftwareTools (2000).[3] Wikipedia Contributors.
List of Intel microprocessors — Wikipedia, TheFree Encyclopedia . https://en.wikipedia.org/w/index.php?title=List_of_Intel_microprocessors&oldid=974087704 . [Online; accessed2-September-2020]. 2020.[4] W Diggin and M Diggin. Using the discrete radon transformation for grayscaleimage moments . 2020. arXiv: .[5] GR Gindi and AF Gmitro. “Optical feature extraction via the Radon trans-form”. In:
Optical Engineering
IRE Trans-actions on Information Theory
Moments and Moment Invariants in PatternRecognition . Wiley, 2009.[8] RJ Prokop and AP Reeves. “A survey of moment-based techniques forunoccluded object representation and recognition”. In:
CVGIP: GraphicalModels and Image Processing1996 IEEE International Conference onAcoustics, Speech, and Signal Processing Conference Proceedings