Abstract

Performing optimal bit-allocation with 3D wavelet coding methods is difficult, because energy is not conserved after applying the motion-compensated temporal filtering (MCTF) process and the spatial wavelet transform. The problem has been solved by extending the 2D wavelet coefficients weighting method with the consideration of the complicated pixel connectivity resulting from the lifting-based MCTF process. However, the computation of weighting factors demands much computing time and vast amount of memory. In this paper, we propose a method, which facilitates the property of sparse matrices, to compute the subband weighing factors efficiently. The proposed method is employed for the 5–3 temporal filter and the 9–7 temporal filter in the lifting-based MCTF process. Experiments on various video sequences show that the computing time and memory needed for the computation of weighting factors are dramatically reduced by applying the proposed approach.

1. Introduction

Scalable video coding, which is developed to support multiple resolutions in a single stream, is an excellent solution to the compression task of the video broadcasting over the heterogeneous network [1, 2]. Compared to H.264/SVC [3], which is recently developed based on the prevailing conventional close-loop codec H.264/AVC, the multiresolution property of 3D wavelet representation based on motion-compensated temporal filtering (MCTF) is a more natural solution to the scalability issue in video coding [4–7]. However, to compete with the great success of scalable coding methods based on H.264, the MCTF-based 3D wavelet video codec must be constantly improved.

For both of H.264/AVC and the wavelet-based encoder, the coding performance is greatly dependent on the technique used to solve the optimal bit-allocation problem. The discrete cosine transform (DCT), which is an orthogonal transform used in H.264/AVC, benefits the procedure of searching the optimal bit-allocation because the variances of the distortions caused by the quantization are equal in the pixel and the frequency domains. However, the wavelet transforms used in the video/image coding (ex: 5–3 and 9–7 filters) are usually not orthogonal so the variances of the distortions in the pixel domain and the frequency domain are not the same. For image coding, the problem has been elegantly solved [8, 9] by assigning different weights to the subbands, resulting in equivalent reconstruction error variance and quantization error variance in the frequency domain. The method was directly extended for 3D wavelet coding previously [10], but the results are not satisfied because the energy difference between the pixel and frequency domains results not only from using a biorthogonal wavelet but also from the MCTF process, which imposes a different connectivity status (single-connected, multiple-connected, or unconnected) on each pixel during motion prediction. Thus, the direct extension of the methods in [8, 9] from 2D wavelet coefficients to 3D wavelet coefficients cannot solve the 3D wavelet optimal bit-allocation problem because they do not take account of the status of the pixel connectivity in the MCTF process.

To fix this problem, a method, which takes account of both the biorthogonal property of the wavelet and the pixel connectivity resulted from the motion compensation in the MCTF process, is proposed to equal the variances of the distortion in the pixel and frequency domain by assigning weighting to subbands and, therefore, solves the bit-allocation problem in 3D wavelet coding [11]. Although the results show the method is effective, lots of computing time and memory are needed during the computation of the weighting factors and it causes that the method becomes a heavy load of the encoding process. To make the computation of weighting factors more efficient, we must develop an approach to decrease the time and memory used by the process. For this purpose, a special data structure and corresponding operations are proposed to save the time and memory during the computation of the subband weighting factors in this paper.

The remainder of this paper is organized as follows. In Section 2, we briefly review the MCTF-EZBC (embedded zero block coding) [12] wavelet coding scheme. In Section 3, we give a summary of how to derive the spatial and temporal weighting factors. In Section 4, we explain the proposed method used to save time and memory during the computation of the subband weighting factors. We also report the experiment results that compare the time and memory used by the traditional method and our method in Section 5. And finally the conclusion is given in Section 6.

2. A MCTF-EZBC Wavelet Video Coding Scheme

In a MCTF-EZBC wavelet video coding scheme, the video frames are first decomposed into multiple wavelet subbands by spatial and temporal wavelet transforms, then the quantization and entropy coding are sequentially applied to the wavelet subbands. According to the order of the spatial and the temporal decompositions, the wavelet coding schemes can be categorized into two categories: 2D+T (spatial filtering first) and T+2D (temporal filtering first). However, no matter 2D+T or T+2D scheme is applied, the spatial and the temporal filterings can be described independently.

The purpose of spatial filtering is separating low and high frequency coefficients from a video frame. The spatial filtering usually consists of multiple sequential 2D wavelet decompositions. In a 2D wavelet decomposition, the input signal, which is represented by a by two-dimensional matrix, is decomposed into four by two-dimensional matrices, which are denoted by , , , and . For these subbands, the previous letter means that the subband contains the low () or high () frequency coefficients after the horizontal 1D wavelet transform, and the following letter means the subband contains the low () or high () frequency coefficients after the vertical 1D wavelet transform. After the decomposition, subband is sequentially taken as the input of the next level spatial decomposition.

If we let and denote the analyzing matrices in the -th spatial decomposition, the corresponding wavelet subbands can be computed as where , , , and correspond to the , , , and subbands, respectively. If there is another spatial decomposition applied to the frame, the subband is further decomposed, by the -th analyzing matrices and . According to this definition, the subband is the frame before any spatial decomposition, and any wavelet subband after the spatial filtering can be indexed by , which represent the subband is the -th subband after the -th spatial decomposition. In Figure 1(a), the example shows the wavelet subbands with the indices after performing the spatial filtering consisting of two spatial decompositions.

(a) (x,y)

(b) (m,n)

(a) (x,y)
(b) (m,n)

Figure 1 

An example of indexing the spatial and temporal decomposed subbands. The above figure shows the indices of the spatially decomposed subbands, in which the subband indexed by represents the subband is the -th after the -th spatial decomposition. The below figure shows the indices of the temporally decomposed subbands, in which the subband indexed by represents the subband is the -th after the -th spatial decomposition. Accordingly, any subband obtained by the MCTF can be indexed by .

To reconstruct from the wavelet subbands , , , and , the synthetic matrices and are used in the synthetic procedure, which can be represented by several matrix operations as

Compared to the spatial filtering, which takes a frame as its input, the temporal filtering takes multiple frames (T+2D) or the subands with the same spatial indices (2D+T) as its input. The temporal filtering consists of multiple temporal decompositions. However, unlike the spatial filtering, which only needs several 2D wavelet transforms during the procedure, the temporal filtering needs multiple 1D wavelet transforms and motion estimation/compensation to specify the input coefficients of each wavelet transform. So, the temporal filtering is more complicated compared to the spatial filtering.

Usually, a temporal decomposition is implemented by the lifting structure, which consists of the predicting stage and the updating stage. As depicted in in Figure 2, the H-frames (denoted by ) are firstly obtained after the predicting stage, then the L-frames (denoted by ) are obtained after the updating stage. The lifting structure can be described by the matrix operations, in which a frame is represented by a one-dimensional matrix . If each input frame of the temporal filtering has the size by and can be represented by a by two-dimensional matrix , the one dimensional matrix has the size and can be mapped from by letting .

Figure 2 

An example of the lifting structure used in the temporal decomposition. The subbands , which contain high frequency coefficients, are first obtained after the prediction stage. Then, the subbands , which contain the low frequency coefficients, are obtained after the update stage.

The motion vectors are computed before the predicting and updating stages. We use the , which is a two-dimensional by matrix, to denote the motion vectors obtained by predicting the -th frame from the -th frame in the -th temporal decomposition. Then, the predicting stage of the -th temporal decomposition can be written as where is the corresponding coefficient of the wavelet filter. Since the index represents the -th lifting stage, the term represents the -th frame which is obtained after the -th lifting stage. If the 5–3 wavelet transform is used, the value of in the predicting stage is −0.5.

In the updating stage, the inverse motion vector matrix , which has the same size as and can be calculated from the matrix [13], is used to compute the decomposed L-frame as where is the corresponding coefficient of the wavelet filter. If the 5–3 wavelet transform is used, the value of in the updating stage is 1.

To understand how to represent the motion vectors with a matrix, an example is given in Figure 3. And these two matrices are constructed as follows:

Figure 3 

An example of MCTF motion estimation. The types of connectivity pixels in the example are single-connected, multiple-connected, and unconnected pixels. The corresponding prediction and update matrices are given in (5), where indicates that the -th pixel in frame is predicted by the -th pixel in frame . Note that means the -th pixel in frame is updated by the -th pixel in frame .

After the -th temporal decomposition, if another temporal decomposition is needed for the temporal filtering, the L-frames generated by the current temporal decomposition are taken as its input. So, we let , which is the -th frame of the -th temporal decomposition, as , which is the output -th frame of the -th temporal decomposition. Although the H-frames do not participate in the -th temporal decomposition, we still arrange the indices to them by letting , in which is the number of the input frames in the -th temporal decomposition. So, any frames obtained by the temporal filtering can be indexed by , which represents it is the -th subband after the -th temporal decomposition.

Usually, the temporal decompositions are sequentially performed until the output has only one L-frame. Since the decomposed frame can be synthesized, only the frames which cannot be synthesized are necessary to reconstruct all the frames. In Figure 1(b), an example, in which the size of the group of pictures (GOP) is 8, shows the indices of the frames that cannot be synthesized from the decomposed frames.

To recover the original frames, the synthesis is applied to the decomposed frames. The frame which is decomposed last is synthesized first, and vice versa. In the procedure of the synthesis, the inverse updating is firstly performed as where is the coefficient of the wavelet transform used in the inverse updating stage. If the temporal 5–3 filter is used, the value of is −1 in the inverse updating stage. After the inverse updating stage, the inverse predicting stage is performed as If the temporal 5–3 filter is used, the value of is 0.5 in the inverse updating stage. For some wavelet filters, such as 9–7 filter, a temporal decomposition needs multiple level lifting structure and can be easily extended by cascading multiple one level lifting structures.

After the spatial and temporal filterings, the quantization and entropy coding are applied to these wavelet subbands. The coefficients of these subbands may be quantized by scalar or vector quantization, then the quantized coefficients are coded without loss by the entropy coder. The method is common in the still image compression standard, such as JPEG 2000 [14]. In the decoding process, the quantized coefficients are obtained by decoding the received bitstream, then the subbands are rescaled according to the quantization step used in the encoder. The advantage of separating the quantization and the entropy coding is that the quality of the reconstructed video can be predicted according to the quantization step.

The quantization and the entropy coding can also be combined with the bitplane coding method, such as the EZBC entropy coder [12]. In these methods, the rates allocated to the subbands are calculated first, then the entropy coder encodes the subbands with these rates. The advantage of the scheme is that the rates of the subbands can be any nonnegative integers and the performance of the bit-allocation can be improved accordingly.

No matter the quantization and entropy coding are combined or separated, the bit-allocation greatly affects the coding efficiency. However, because the energy in the pixel domain can be altered after application of the spatial wavelet transform, temporal wavelet transform, and motion estimation in the MCTF process, the bit-allocation is prevented from achieving optimum. To preserve the energy between the pixel and the wavelet domain, we derive the weighting factors of the decomposed subbands. We explain how to derive these weighting factors in the next section.

3. Subband Weighting

The weighting factor indicates how much a unit quantization power in the subband contributes to the overall distortion in the reconstructed GOP, that is, for the subband indexed by , which means the subband is the spatial -th subband after the -th spatial decomposition, and the temporal -th subband after the -th temporal decomposition (an example of the index method is depicted in Figure 1), a weighting is given to satisfy where is the variance of the distortion in the pixel domain, is the variance of subband's distortion in the wavelet domain, and is the set including the indices of the subbands which cannot be synthesized during the reconstruction. The weighting factor can be computed by where is the spatial weighting factor and is the temporal weighting factor with respect to the subband indexed by .

The spatial weighting factors can be directly computed according to the error propagation model mentioned in [15]. However, to compute the temporal weighting factors, the effect of the motion compensation must also be considered, and the approach is mentioned in [11]. In the following paragraph, we explain how to derive the spatial and temporal weighting factors, respectively.

3.1. Spatial Subband Weighting

In this section, we review the method used to derive the weighting factors of spatial wavelet transforms in [15]. According to (2), the reconstructed subband can be written as Let denote the error matrix resulting from the quantization of an image ; according to (10), we have Using the property of Kronecker products, where , are column vectors constructed row-wise from the matrices , , respectively. Applying this identity to (11), we now have where is the column vector constructed row-wise from the matrices . The reconstruction mean square error (MSE) of an image is The equation can be solved under the high bit rate assumption. At a high bit rate, it is assumed that the quantization errors of wavelet subbands are white and mutually uncorrelated [16]. So, the following identities for the vector representation of errors in subbands and are obtained: where is the MSE of any element in the vector . By substituting (13) into (14) and applying the identities in (15), we obtain

The above derivation shows that the MSE measured before the wavelet transform is the weighted sum of subbands' MSE after the wavelet transform. Note that the weighting factors are determined only by the filters. According to the discussion in the previous section, the subband indexed by is further decomposed in the -th spatial decomposition. So, the spatial weighting factors of a multilevel decomposition case can be solved by the recursive substitution. To summarize, the MSE of the original frame can be computed as in which is the total number of the spatial wavelet transforms. So, the multilevel spatial weighting factor can be computed from one-level spatial weighting factor as

3.2. Temporal Subband Weighting

We follow the scheme proposed in [11] to explain how to derive and compute the temporal subband weighting factors. The temporal 5–3 filter, which is usually used in the wavelet video coding, is used to explain the procedure. And the derivation of the weighting factors with respect to the temporal 9–3 filter, which is the other filter usually used in the wavelet video coding, is detailed in [11].

We first explain how to compute the weighting factor , and the weighting factor is expected to satisfy where is the total number of the input subbands in the -th temporal decomposition.

If the temporal 5–3 filter is used in the MCTF process, the filter coefficients used in the predicting and the updating stages are −0.5 and 1, respectively. To simplifying the notation, is denoted by , and is denoted by , then (3) and (4) can be written as

Let represents the quantization error resulting from lossy source coding, and denotes the reconstructed frame. From (21), we can obtain the reconstruction error of the -th frame as follows:

Substituting (22) into (20) for and , we obtain the following reconstruction error of the -th frame:

Equations (22) and (24) represent a motion-dependent error propagation model for a one-level MCTF process using the 5–3 temporal wavelet filter. The reconstructed MSE of the -th frame can be derived as follows: by (14), By defining and using the high bit rate assumption that derives the identities in (15), the last three cross-terms are zero and (24) becomes Following the same derivation, the reconstruction error of the -th frame is

By applying derivations similar to those in (25) and (26) to each input frame, we can relate the reconstruction errors before and after the temporal decomposition by a linear relation:

where The matrix is used to adjust the positions of the subbands. To measure the consequence of quantizing a subband after the temporal decomposition, we should aggregate the errors induced in all the subbands before the temporal decomposition. This is exactly the summation of the corresponding column of the matrix for the subband. For example, the weighting factor of the reconstruction error resulting from subband is the summation of the values in the column of and denoted by . The following equation shows how to compute the weighting factor :

If the temporal filtering consists of total temporal decompositions, the variances of the distortions of the frames in the pixel domain can be computed as Please note the index of subband can be changed from to , which represents the subband is obtained as the -th subband after the -th temporal decomposition (see Figure 1(b)). So, the temporal weighting factor used in (9) can be computed as

Since the spatial and temporal weighting factors are both derived, the performance of the bit-allocation method can be improved accordingly. In our previous work [11], we compared three bit-allocation methods, which consider spatial weighting factors, spatial-temporal weighting factors without motion vectors, and spatial-temporal weighting factors with motion vectors, respectively. And the the results showed the bit-allocation method considering spatial-temporal weighting factors with motion vectors greatly outperforms the other two methods. However, the computation of the best weighting method costs lots of computing resources. So, the method, which reduces the needed resources, is proposed in the next section.

4. Fast Computation of the Weighting Factors

As our discussion in the previous section, (16) is used to compute the spatial weighting factors, and (24) and (25) are used to compute the temporal weighting factors. However, a lot of computing resources are needed to obtain these weighting factors. Let the size of an input frame is by . In (16), (24), and (25), if we use matrices to represent , , and , the size of each matrix is by . If bytes are used to store a coefficient in a matrix, a matrix would take bytes in storage. For a CIF frame (of size ), the matrix representation needs at least 39204 megabytes for only one matrix. Thus, using the matrix representation to compute the weighting factors is extremely inefficient since it needs to store many by matrices.

Another issue that handicap using the matrix representation for weighting factor is the computing time to perform matrix addition and multiplication. If we use the conventional matrix operations, which are defined as follows: the time complexity of the matrix multiplication in computing the weighting factors is bounded by . For a CIF frame, it costs approximately real number multiplications for performing one matrix multiplication. Fortunately, these matrices are usually very sparse, that is, only few coefficients in these matrices are not zero. By taking advantage of the property, we use an efficient data structure to store the required data in computing the weighting factors.

Thus, we propose to use arraylist representation. An arraylist consists of a one-dimensional array comprised of pointers to linklists. Figure 4 is an example to demonstrate the arraylist structure. In the proposed data structure, two arraylists, which, respectively, record the positions and values of the nonzero coefficients, are used to represent a sparse matrix of large size. In Figure 4, the matrix is an abbreviation of the matrix . The arraylists and , respectively, record the positions and values of nonzero coefficients in the matrix . Because the matrix has 6 rows, both arraylists have 6 elements. For the arraylist , the -th linklist links to the locations of nonzero elements at the -th row of the matrix. The arraylist has a similar structure at that of , instead of linking to locations of nonzero elements, links to the values of the nonzero elements. That if the position of the nonzero coefficient is recorded in , we can deduce that the -row and -column element of the matrix is a nonzero coefficient, and whose value is stored in .

Figure 4 

An example of the arraylist representation. The matrix is in (5). According to whether row (above) or column (below) is considered first, the arraylist representation can be categorized into row representation and column representation. The arraylists in the left are used to store the locations of nonzero coefficients, and the arraylists in the right are used to store the values of nonzero coefficients.

If the matrix is sparse enough, the proposed arraylist data structure can greatly reduce the memory because only the nonzero locations and values are stored in the arraylists. If we assume that there are nonzero coefficients in each row of an by matrix. To construct an arraylist, we need a one-dimensional array of elements, and linklists, one for each element in the arraylist. Each linklist has elements. Because our structure uses two arraylists, it needs bytes, where is the number of bytes used to store a coefficient. Thus, our data structure has reduced to size to represent an by matrix from bytes to bytes. The proposed data structure in many examples can reduce the memory to store a of size matrix from 39204 megabytes to about 1 megabytes.

After solving the matrix storage, we use the matrices in (5) to explain how to perform matrix addition and multiplication with the proposed data structure. To simplify the notation, and are used to represent and in (5), respectively. In matrix addition, we first create two new arraylists and to store the result of the addition with For all , we move the linklists and to the tails of the linklists and , respectively. Then, all linklists in are examined to search if any two elements have the same value. If the two elements, denoted by and with , have the same value, then the value of is assigned to the element . Finally, the element is removed for the consistency of the data structure. If no elements of a linklist in have the same value for all the rows, the matrix addition with respect to the proposed data structure is completed. The results of the matrix are stored in and . Figure 5 demonstrates the matrix addition by using the proposed data structure.

Figure 5 

An example of adding two matrices with the arraylist representation. First, two arraylists and are created by cloning and . For all, we move the linklists and to the tails of the linklists and , respectively. For any linklist in , if any two elements, which are denoted by and , have the same value, the elements and are removed (if ), then value of the element is assigned as . This procedure continues until no elements in a linklist have the same value. Then, the arraylists and are the arraylist representation of .

To obtain the matrix multiplication with the proposed data structure, the coefficients are needed for all . We create two empty arraylists and , each of which arraylists has linklists to store the result of matrix multiplication. To compute , we examine the linklists and to search for the elements of the same value, then we multiply their corresponding values in and . The value of is the summation of the products obtained by the process. To store the value of , the values and are appended to the tail of the -th link lists of and , respectively. The result of performing matrix multiplication with the specified data structure is depicted in Figure 6.