Modified Hilbert Curve for Rectangles and Cuboids and Its Application in Entropy Coding for Image and Video Compression

In our previous work, by combining the Hilbert scan with the symbol grouping method, efficient run-length-based entropy coding was developed, and high-efficiency image compression algorithms based on the entropy coding were obtained. However, the 2-D Hilbert curves, which are a critical part of the above-mentioned entropy coding, are defined on squares with the side length being the powers of 2, i.e., 2n, while a subband is normally a rectangle of arbitrary sizes. It is not straightforward to modify the Hilbert curve from squares of side lengths of 2n to an arbitrary rectangle. In this short article, we provide the details of constructing the modified 2-D Hilbert curve of arbitrary rectangle sizes. Furthermore, we extend the method from a 2-D rectangle to a 3-D cuboid. The 3-D modified Hilbert curves are used in a novel 3-D transform video compression algorithm that employs the run-length-based entropy coding. Additionally, the modified 2-D and 3-D Hilbert curves introduced in this short article could be useful for some unknown applications in the future.


Introduction
Entropy coding plays a critical role in data compression, such as image and video compression, etc. Two commonly used algorithms for entropy coding are Huffman coding [1] and arithmetic coding [2]. In terms of compression efficiency, arithmetic coding is preferred. However, arithmetic coding has a higher computational complexity because it requires multiplication and division during the coding process. To resolve the complexity issue, approximations are used in binary arithmetic coding algorithms, such as [3][4][5], etc. These binary arithmetic coding algorithms are practical algorithms because the coding of a multiple symbol source can always be converted to coding of a sequence of binary symbol sources. For example, in image compression, the bit-plane coding method [6] and the symbol grouping coding method [7,8] eventually convert the quantized coefficients to binaries to code.
Although the existing binary arithmetic coding algorithms solved the computational complexity issue, it is still not computationally efficient in extremely low entropy conditions because arithmetic coding algorithms encode symbols one by one. For example, to code a binary source with the probabilities p = 0.999 for the symbol "0" and 1 − p = 0.001 for the symbol "1", arithmetic coding needs to code 999 "0"s on average before it codes a "1". On the other hand, run-length-based entropy coding [7,8] is much more computationally efficient for low-entropy coding situations, as it does not need to code the "0"s one by one. Note, low-entropy sources are very common in compression. For example, in subband image compression, most of the quantized coefficients in a subband are zeros. The positions of the non-zero coefficients in a subband normally form an extremely low-entropy source.
For non-stationary binary sources, the binary arithmetic coding algorithms use probability estimators to adapt to the probability variations; whereas the run-length-based entropy coding uses the symbol grouping method to handle non-stationary binary sources [7]. For coding 2-dimensional (and higher-dimensional) subband coefficient arrays, binary arithmetic coding can estimate the probabilities from the coded adjacent coefficients in different spatial directions (context modeling). However, for the run-length-based binary entropy coding, the 2-D coefficient array needs to be scanned into a 1-D array before the run-length coding can be performed. As a result, for run-length-based binary entropy coding, exploiting probability estimations in different spatial directions on the 2-D array before scanning is very difficult. Thus, to achieving coding efficiency, variations in the original 2-D signal need to be maximumly kept into the scanned 1-D array, which requires that nearby elements in the 2-D array are still nearby in the 1-D scanned array. Ideally, adjacent elements in the 2-D array are required to be adjacent elements in the 1-D scanned array, which is impossible, as can be easily shown. However, different scan routes lead to different scatterings of the 2-D nearby elements. Thus, scan routes with small scattering are desired. The Hilbert curve [9,10] is such a route. Figure 1 shows a 2-D Hilbert curve. As can be seen, the Hilbert curve tries to keep the 2 k × 2 k (k = 1, 2, . . .) elements in the 2-D array together in the scanned 1-D array. In fact, the locality-preserving feature of the Hilbert curve has been extensively studied [10][11][12][13]. Thus, using the Hilbert curve scan route, the variations within a 2-D subband coefficient array are greatly kept into the 1-D scanned coefficient array. Indeed, combining the Hilbert scan with the symbol grouping method, an efficient entropy coding was achieved, and high-efficiency image compression algorithms were obtained [7,8]. However, 2-D Hilbert curves are defined on a square of sizes 2 i × 2 i (i = 1, 2, . . .) [9,14]. In other words, not only the array shape is a square, but also the side length of the square can only be the powers of 2, i.e., 2 i . Yet, a subband is normally a rectangle of arbitrary sizes. It is not straightforward to modify the Hilbert curve from the 2 i × 2 i squares to an arbitrary rectangle. In [7,8], details of this modification are not provided.
In this short article, we provide the details of constructing the modified 2-D Hilbert curve of arbitrary rectangle sizes. Furthermore, we extend the method from a 2-D rectangle to a 3-D cuboid. The entropy coding in the 3-D transform video compression algorithm introduced in [15] uses the 3-D Hilbert curve. Test results show that the algorithm is promising. However, the original 3-D Hilbert curve is defined on a cube of side length of power of 2, i.e., size 2 i × 2 i × 2 i (i = 1, 2, . . .) [16]. Because the 3-D modified Hilbert curves for cuboids were not available at the time, videos were cropped to the size of 1024 × 1024 for testing the algorithm prototype in [15]. The extension from a 2-D rectangle to a 3-D cuboid makes the prototype proposed in [15] a practical video compression algorithm that accommodates arbitrary rectangle video sizes. Further, Hilbert curves have been widely used in many applications, such as image data encryption, query, and retrieval [17,18], etc. Extending the original Hilbert curves to arbitrary size rectangles and cuboids could be useful for some unknown applications in the future.

The 2-D Modified Hilbert Curve
The original 2-D Hilbert curve connects a square array of the size of 2 i × 2 i . We denote it as the i th order Hilbert curve. Hilbert curves of orders i = 1, i = 2, and i = 3 are respectively shown in Figure 1a-c. There is an important property of the Hilbert curve. As can be seen from Figure 1a-c, the starting and the ending points (green and red points in the graphs) are always on one side of the 2 i × 2 i square. Since the starting and the ending points are at the ends of one side of the Hilbert curve square, a Hilbert curve of any order i can easily be represented by a square with labeled starting and ending points like Figure  1d, when the internal structure does not need to be shown. Now, we can show that an (i + 1) th order Hilbert curve can be easily constructed from four i th order Hilbert curves. First, replace the 4 points in the 1 st order Hilbert curve in Figure 1a with four i th order Hilbert curves whose starting and ending points are arranged as indicated in Figure 1e; then, connect the ending points and starting points of the i th Hilbert curves orderly as indicated in Figure 1e, and the (i + 1) th order Hilbert curve is constructed. With the 1 st order Hilbert curve provided by Figure 1a and the method of constructing the (i + 1) th order Hilbert curve from the i th Hilbert curve, we can construct the Hilbert curve of any order by mathematical induction. Now, our task is to construct a scanning route that is close to the Hilbert curve for a 2-D array of size W × H, where W and H are the element numbers along the vertical and horizontal direction, respectively.
The basic idea is to divide the W × H rectangle array into N small square arrays of the size 2 i n × 2 i n , n = 1, 2, . . . , N, and i 1 ≥ i 2 ≥ . . . ≥ i N . For example, a 12 × 8 rectangle array can be divided into three 2 i n × 2 i n square arrays with i 1 = 3 and i 2 = i 3 = 2. Apparently, within each square, an (i n ) th order Hilbert curve can be easily constructed as shown in Figure 2a. By appropriately arranging the directions of each Hilbert curve, the 3 Hilbert curves are connected to form the desired route, as shown in Figure 2b.
To form a route close to the Hilbert curve, there are two requirements. First, the largest i n needs to be selected in order, i.e., select the largest i 1 first, then the largest i 2 , . . . , etc. Without this restriction, the constructed curves may deviate from the Hilbert curve significantly. As an extreme example, for the 12 × 8 = 96 points in Figure 2, one can simply choose i 1 = i 2 = . . . = i 96 = 0, i.e., use the 96 points as 96 small square arrays. In this case, there are too many possible routes. Most of them are not close to the Hilbert curve at all, for example, the raster scan. Secondly, the ending point of the Hilbert curve in the n th square must be adjacent to the starting point of the Hilbert curve of the (n + 1) th square, like the example shown in Figure 2b. For design convenience, the first requirement may not be satisfied strictly sometimes. However, the second requirement must be satisfied.  Based on the above observations, the procedures of our design method are provided as follows (note, the method is not unique). Without loss of generality, we consider rectangle arrays with W ≥ H: 1.

2.
Construct a modified Hilbert curve within the sub-rectangle S 1 of size 2 m 1 × H. S 1 is a special rectangle with the width being 2 m 1 . The starting and the ending points of S 1 need to be at the ends of the S 1 's top width, as indicated in Figure 3. The construction details for this step is provided shortly. (Note, we use "width" and "height" to represent the number of elements in the two orthogonal directions of a rectangle array throughout the paper. They are not the geometric lengths. Do not get confused with the illustrating diagrams used in the paper.) 3.
Once the modified Hilbert curve for the rectangle array S 1 is constructed, the construction of the remaining rectangle array S 2 goes back to step (1) with the starting and the ending points indicated in Figure 3. However, the array size of S 2 is less than half of the original W × H rectangle.

4.
By iterating steps 1 to 3, the subsequent remaining S 2 's become smaller and smaller very quickly, with speed faster than 0.5 k , where k is the iteration number. The iteration stops when the remaining S 2 is of size 2 l × H , and the construction is complete. Note, a size of 2 l × H can always be achieved because the smallest H can be 1.
Note, when iterating the three steps 1 to 3, if for S 2 , the height H is larger than W − 2 m 1 like the situation shown in Figure 3, then H plays the role of the width W for the new iteration on S 2 because we assumed the initial condition of W ≥ H for each iteration. In this case, the design route changes direction because it is always along the width direction, see Figure 3. If for S 2 the height H is still smaller than W − 2 m 1 , the iteration continues along the same design route direction. In the following example, we show a specific design to provide a more intuitive understanding of the procedures.
Suppose we want to design a modified 2-D Hilbert curve for a practical size of W = 1920/8 = 240 and H = 1080/8 = 135, which is the subband size from the 8 × 8 subband decomposition on a 1920 × 1080 image, the standard HDTV size. Following the design procedures: Iteration 1: m 1 = 7, and sub-rectangle arrays S 1 and S 2 for the 1st iteration are: We start from the (m 1 ) th order Hilbert curve, whose height is 2 m 1 . Thus, we need to reduce the height by ∆H = 2 m 1 − H. Recall that an integer B can be converted into its binary format b s b s−1 . . . b 1 b 0 , i.e., B can be decomposed as where the b i 's are either 0 or 1. Decomposing ∆H using (1), we can perform a reduction of ∆H step by step, with each step achieving a reduction of 2 i points. For example, suppose ∆H = 13. Then, from (1), we have ∆H = 8 + 4 + 1, i.e., b 3 = 1, b 2 = 1, b 1 = 0, and b 0 = 1. Thus, we need to reduce 8 points, 4 points, and 1 point on H to achieve the total reduction of ∆H = 13 points.
To perform a reduction of 2 i points, first, a reduction of 2 i points on the (i + 1) th order Hilbert curve is straightforward. By inspection, it can be seen that the top two sub-squares in Figure 1e can be removed, leading to Figure 1f. Then, a reduction of 2 i points on the (i + 1) th order Hilbert curve is achieved. Next, we observe the following property: For an opening-toward-top Hilbert curve of any order, the bottom sub-Hilbert curves always have the same opening-toward-top orientation. As an example, in Figure 5a, a 4 th order opening-toward-top Hilbert curve is plotted. The 4 th order Hilbert curve can be represented using the structure of Figure 5b, i.e., the main structure is an opening-toward-top 1 st order Hilbert curve with 4 sub-curves, which are four 3 rd order Hilbert curves denoted by four small shaded squares. As just described above, the removal of the top two sub-squares, or, equivalently, a reduction of 8 points in H in this case, can be easily achieved. Now, it is important to observe from Figure 5b that the bottom two shaded squares, i.e., the bottom two 3 rd order sub-Hilbert curves, are also opening-toward-top Hilbert curves. When the original 4 th order Hilbert curve is represented using the structure of Figure 5c, for each of the bottom two opening-toward-up 3rd order sub-Hilbert curves, the removal of the top two sub-squares, or equivalently a reduction of 4 points on H in this case, can be achieved. Similarly, it can be seen from Figure 5a,d that reductions of 2 points and 1 point on H can be achieved. Figure 6 shows a specific example of how the 4 th order Hilbert curve is reduced by ∆H = 13 = 8 + 4 + 1. In Figure 5a, the original 4 th order H = 16 Hilbert curve is shown. Reductions on H by 8, 4, and 1 in each step are respectively shown in Figure 6a-c. After all the sub-reductions are completed, the total reduction of ∆H = 13 is achieved in Figure 6d.

Condition (B) H > 2 m 1 : We need to increase the height by
Thus, similar to condition (A), ∆H is decomposed by (1), and we can increase H by 2 i points in each step because the increase in height by 2 i points on the (i + 1) th order Hilbert curve can be achieved using the modification from Figure 1e-g.
With the height of the (m 1 ) th order Hilbert curve reduced or increased to H, procedure 2 of the iterative design method described earlier is performed. The modified 2-D Hilbert curve on a W × H rectangle array using the iterative design procedures is completed.
To provide more intuitions, the 2-D modified Hilbert curves of sizes 27 × 17, 27 × 18, and 27 × 19, constructed using the proposed method, are shown in Figure 7. In addition, the MATLAB codes implementing the proposed method are available in [19]. As can be easily checked, the algorithm runs reasonably fast, and the results can be obtained instantly.

The 3-D Modified Hilbert Curve
In [15], the binary run-length-based symbol grouping entropy coding method is used in video compression. For the 3-D transform video compression algorithm introduced in [15], conventional motion compensation is not used in order to improve the computational complexity. Instead, a 4-band SCWP transform [20] is performed along the time dimension. In other words, the first step of the video compression algorithm is a 3-D transform. Thus, the transformed coefficients are 3-D subband arrays.
In entropy coding of the quantized 3-D subband coefficient arrays, the 3-D Hilbert curve scan was used to maximally keep the correlations in the 3-D subband into the 1-D scanned array. Because the original 3-D Hilbert curves are for cubes of side length 2 i , in [15], the 1920 × 1080 test videos were cropped to a size of 1024 × 1024 for testing. Apparently, to accommodate an arbitrary rectangle video size, the original 3-D Hilbert curve needs to be modified. Below, we extend the modification method introduced in Section 2.1 for 2-D arrays to 3-D conditions. The 3-D arrays are of size W × H × D, where W and H are the width and height of the cuboid array, D is the third dimension denoting the depth here, which corresponds to the time dimension of the input video.
The depth D of the 3-D decomposed subband is determined by the parameter "group of pictures" (GOP), which is normally selected to be the powers of 2. As a result, the depths D are also the powers of 2, i.e., D = 2 d , d = 0, 1, 2, . . . For example, the GOP used in [15] was 32, which leads to the depths D of the 3-D subbands being 8, 4, 2, or 1 (for details, please refer to [15]). Furthermore, compared with W and H, D is normally much smaller in our case.
We begin from the original 3-D Hilbert curve. Again, we denote a 3-D Hilbert curve of size 2 i × 2 i × 2 i , the i th order 3-D Hilbert curve. Figure 8a,b, respectively, show the 1 st order and the 2 nd order 3-D Hilbert curve. Similar to the 2-D situation, the starting and the ending points of any order 3-D Hilbert curves are at the two ends of one side of the cube. Therefore, as shown in Figure 8c, when the internal structure is not needed, a 3-D Hilbert curve of any order i can be represented by a cube with the starting and ending points labeled. With this simple and intuitive representation, the construction of the (i + 1) th order 3-D Hilbert curve from the i th order 3-D Hilbert curve can be easily demonstrated by Figure 8d. By mathematical induction, given (1) the 1 st order 3-D Hilbert curve and (2) the method of constructing the (i + 1) th order 3-D Hilbert curve from the i th order 3-D Hilbert curve, the 3-D Hilbert curve of any order can be constructed.
Without loss of generality, assume W ≥ H. As mentioned in our application, D is the powers of 2 and is normally much smaller than W and H. In other words, the modified 3-D Hilbert curve is of size W × H × D = W × H × 2 d , and D is much smaller than W and H. Next, we consider the situation where W and H are not multiples of D, but W ≥ D. In this case, we can still borrow the 2-D construction structure, i.e., use the 2-D iterative route similar to Figures 3 and 4c for the height-width (W-H) surface of the 3-D cuboid. This is similar to what we performed in Figure 9c,d, where W and H are multiples of D. The difference is that, in this case, ∆H is not a multiple of D. In this case, we need to consider the non-zero terms in ∆H = ∑ b i 2 i that are smaller than D, i.e., the 2 i < D terms. The 2 i ≥ D terms are multiples of D, which is the situation previously considered. The 2 i < D terms in ∆H, however, need to be handled on the 3-D cubes at the bottom. For example, if in Figure 9d we want to construct a modified 3-D Hilbert curve of size 2D × 1.5D × D instead of 2D × 2D × D, then the bottom cubes in Figure 9d Figure 4c, which was used for constructing the 240 × 135 2-D modified Hilbert curve. Nevertheless, a construction may end up with a residue S 2 on the W-H surface, whose height and width are both smaller than D. For example, using the above method to construct the modified 3-D Hilbert curve of the size 243 × 135 × 8, we end up with a residue cuboid of size 3 × 7 × 8. Then, the 2-D Figure 3 iterative procedure cannot proceed for the W-H surface anymore. We have to consider the situation of constructing a modified 3-D Hilbert curve of the size W × H × D with D = 2 d > W ≥ H.
However, in the 3-D transform video application, D is very small. As mentioned above, in [15], the maximum D = 8 even if the GOP used is 32. When D = 8, the maximum residue cuboid is only 7 × 7 × 8, which is very small. For such tiny residue cuboids, using some other routes, such as the raster scan, would not lead to any noticeable effect on the final video compression results. On the other hand, the design for the situation of D = 2 d > W ≥ H is complex, and thus, for the application of coding the 3-D transformed coefficients in video compression, we can just use a simple scan route for the residue cuboids with D > W ≥ H. We implemented in MATLAB such 3-D extension with small D > W ≥ H residue cuboid connected using the raster scan, which is available at [19]. Figure 11 shows a modified 3-D Hilbert curve of size 9 × 6 × 4 (i.e., D = 2 2 ) produced by MATLAB codes.
We will not lengthily go into the design on the condition D = 2 d > W ≥ H. For completeness, we only briefly describe that the design is possible using similar ideas we have used up to now. Note, there can be some other methods to handle the D = 2 d > W ≥ H situation because the design method is not unique.
For the D = 2 d > W ≥ H situation, first, consider the situation where either W or H is a power of 2. Without loss of generality, assume H = 2 h . Observe that the sizes W × H × D, D × H × W, . . . , etc., i.e., all the 6 permutations, are the same for our curve construction task. In order to exploit our previously developed construction techniques, we need to change the roles of W, H, and D. Because D is the longest side, we need to use D = 2 d as the width. Since H = 2 h < D, use H as the depth. Then, the construction finishes nicely in one step.
For the more difficult situation, where both W and H are not powers of 2, decompose the shortest side H into a sum of 2 i using equation (1): H = 2 h 0 + 2 h 1 + . . . + 2 h n , where h 0 > h 1 . . . > h n (h 0 corresponds to the most significant bit, i.e., 2 h 0 > 1 2 H). Then, the construction on the D × W × 2 h 0 = 2 d × W × 2 h 0 cuboid is immediately achieved as described above. To increase the thickness from 2 h 0 to H, the point-increasing operation needs to be along the direction as illustrated in Figure 10d. For the 4 length-increased sub-cubes at the back in Figure 10d, the bottom 2 sub-cubes can use the point-increasing operation we already used, i.e., the one from Figure 10a to Figure 10c, but the top 2 sub-cubes need to use a different point-increasing structure, which is skipped here. We may also need to perform multiple point-increasing operations and then perform a point-decreasing operation to achieve the desired value H, and the operations need to be performed individually for the sub-cubes at the back of the D × W × 2 h 0 cuboid. The sizes of the sub-cubes can be different depending upon the W value. As a result, the implementation is complex. We will not go into the details further since currently, there is no immediate application.

Conclusions and the Near Future Work
We have shown the method of modifying the 2-D Hilbert curve to fit an arbitrary W × H rectangle array and the method of modifying the 3-D Hilbert curve to fit a cuboid array of size W × H × 2 d . These modified Hilbert curves can be used in entropy coding for image and video compression. Furthermore, since the construction of the modified 2-D and 3-D Hilbert curves is not straightforward, the methods presented in this short article could be useful for some unknown applications in the future.
The 2-D modified Hilbert curve has already been used in the run-length-based symbol grouping entropy coding method for lossy and lossless image compression. High compression efficiency is achieved, as shown in [7,8].
Because of using the 3-D Hilbert curve, the video compression algorithm prototype introduced in [15] only tested videos with the cropped size of 1024 × 1024, although some promising results were shown. On applying the 3-D modified Hilbert curve for coding to the 3-D subband coefficients so that the algorithm can handle arbitrary video sizes, together with some other fine tunings, we are completing the video compression algorithm very soon. We will systematically compare the performances of the new video compression algorithm with state-of-the-art video compression algorithms, such as HEVC, etc., in terms of compromise between complexity and compression efficiency. From the preliminary test results shown in [15], we expect that the final completed video compression algorithm using the 3-D modified Hilbert curve developed in this paper will be competitive to state-of-the-art video compression algorithms in certain important situations, such as the compression at the high video quality.