A) Lossless coding system

In this paper, we consider the lossless coding system illustrated in Fig. 1, and discuss on a RCT for decorrelation of four color components, G ₀, G ₁, R, and B, acquired through Bayer pattern color filter array illustrated in Fig. 2. The reversible wavelet transform (RWT) is based on the 5/3 tap integer-type filter bank. The embedded block coding with optimal truncation points (EBCOT) is an entropy encoder [Reference Taubman18]. Note that both of them are defined by the widely used international standard JPEG 2000 for lossless image coding [8].

Fig. 1. Lossless coding system for Bayer pattern images.

Fig. 2. Bayer pattern image and its color components.

B) Non-adaptive RCT (existing method I)

Figure 3 illustrates the RCT with fixed coefficients (existing method I) [9]. It converts four color components by

(1)$$\left[\matrix{ Y_0\cr Y_1\cr Y_2\cr Y_3 }\right]= \left[\matrix{ G_1 -G_0\cr F\lsqb \lpar R+B+F\lsqb W\rsqb \cdot 2\rpar /4\rsqb \cr R-F\lsqb W\rsqb \cr B-F\lsqb W\rsqb }\right]\comma \;$$

where

$$\left\{\matrix{ W = F\lsqb \lpar G_1 +G_0\rpar /2\rsqb \comma \; \cr F\lsqb x\rsqb = x-\lpar x \hbox{mod} \ 1\rpar }\right.$$

and F[x] denotes a floor function that converts x into integer. When the errors generated by this function is negligible, the RCT is expressed as

(2)$${\bf Y} = {\bf K}_4^T \cdot {\bf X}\comma \;$$

where

(3)$${\bf K}_4^T = \left[\matrix{ {-1} & 1 & 0 & 0\cr {1/4} & {1/4} & {1/4} & {1/4}\cr {-1/2} & {-1/2} & 1 & 0\cr {-1/2} & {-1/2} & 0 & 1 }\right]$$

and

(4)$$\eqalign{{\bf Y}& = \lsqb Y_0 \quad Y_1 \quad Y_2 \quad Y_3\rsqb ^{T}\comma \; \cr {\bf X}& = \lsqb G_0 \quad G_1 \quad R \quad B\rsqb ^{T}.}$$

Fig. 3. Non-adaptive RCT (existing method I).

Since the multiplier coefficients are fixed as expressed in (3), this RCT cannot be adaptive to each of input images.

C) Adaptive RCT (existing method II)

Figure 4 illustrates the RCT based on KLT (existing method II) [Reference Iwahashi and Kiya16]. The matrix K₄ in (2) is determined as eigenvectors of the covariance matrix R_X defined as

(5)$${\bf R}_{{\bf X}} = E\lsqb {\bf X}\cdot {\bf X}^T\rsqb \comma \;$$

where E[ ] denotes the ensemble average of each row and column. The KLT converts R_X to

(6)$$\eqalign{{\bf R}_{{\bf Y}} &= E\lsqb {\bf Y}\cdot {\bf Y}^T\rsqb \cr &= {\bf K}_4^T \cdot E\lsqb {\bf X}\cdot {\bf X}^T\rsqb \cdot {\bf K}_4\cr &= {\bf K}_4^T \cdot {\bf R}_{{\bf X}} \cdot {\bf K}_4 \cr &= diag \lsqb \lambda_{0} \, \, \lambda_{1} \, \, \lambda_{2} \, \, \lambda_{3}\rsqb \comma \; }$$

where λ_i, i ∈ {0, 1, 2, 3} denote eigenvalues. As a result, Y is decorrelated adaptively to each of input color images.

Fig. 4. Adaptive RCT (existing method II).

The 4 × 4 matrix K₄ is factorized into a product of 2 × 2 matrices

(7)$${\bf F}\lpar \theta_i\rpar = {\bf P}_t \cdot {\bf G}\lpar \varphi_{i\comma t}\rpar$$

for i ∈ {0, 1, 2, 3, 4, 5} as illustrated in Fig. 4. Each of the matrices is composed of Givens rotation

(8)$${\bf G}\lpar \varphi_{i\comma t}\rpar = \left[\matrix{ {\cos \varphi_{i\comma t}} & {-\sin \varphi_{i\comma t}}\cr {\sin \varphi_{i\comma t}} & {\cos \varphi_{i\comma t}} }\right]$$

and one of the permutations

(9)$$\eqalign{{\bf P}_1 & = \left[\matrix{ {+1} & 0\cr 0 & {+1} }\right]= {\bf G}\lpar 0\rpar \comma \; \cr {\bf P}_2 & = \left[\matrix{ 0 & {+1}\cr {-1} & 0 }\right]= {\bf G}\lpar -\pi/2\rpar \comma \; \cr {\bf P}_3 & = \left[\matrix{ {-1} & 0\cr 0 & {-1} }\right]= {\bf G}\lpar \pi\rpar \comma \; \cr {\bf P}_4 &= \left[\matrix{ 0 & {-1}\cr {+1} & 0 }\right]= {\bf G}\lpar +\pi /2\rpar }$$

to avoid the singular point problem [Reference Iwahashi and Kiya16]. Since all the rotation angles θ_i in (7) are determined according to each of input images, this RCT is adaptive.

A) Semi-adaptive RCT and its advantages

In this paper, we propose an RCT illustrated in Fig. 5. A redundant part of the existing method II in Fig. 4 is eliminated. The rotation angle θ₀ is fixed, and the total number of rotations F is reduced from six to four. In this sense, computational load is reduced. Advantages of the proposed method are summarized as below.

Fig. 5. Semi-adaptive RCT (proposed method).

(i) The optimum angles to be included into the bit-stream are reduced from six to three (50.0%). (ii) Size of the eigenvalue problem to be solved is reduced from 4 × 4 for the six angles to 3 × 3 for the three angles (56.3%). (iii) The total number of lifting steps, which means each equation in (10), is reduced from 18 to 12 (66.7%), since one rotation contains three lifting steps. (iv) Total delay due to the lifting steps is reduced, since a lifting step waits for a calculation result of the previous lifting step. In RCT, it accumulates and causes delay from input to output. (v) Total amount of rounding error in Y _i, i ∈ {0, 1, 2, 3} is expected to be reduced, since the rounding functions which generate the error in a lifting step are reduced. We confirm it in Section IV-A. Rational of our simplification for reduction of computational load is described in Section III-C.

B) Implementation issue

To make the transform reversible in Figs 4 and 5, the rotation G(φ_i,t) in (7) is implemented in the lifting structure. In case of t = 1, output signals are calculated as

(10)$$\left\{\matrix{ x'_1 =x_1 + O\lsqb f_0 \cdot x_0 \rsqb \comma \; \cr y_0 = x_0 +O\lsqb f_1 \cdot x'_1 \rsqb \comma \; \cr y_1 =x'_1 +O\lsqb f_2 \cdot y_0 \rsqb \comma \; }\right.$$

with multiplier coefficients f ₀, f ₁, and f ₂ as illustrated in “type 1” in Fig. 6. Neglecting rounding errors generated by the rounding function O[ ], Right-hand side of (7) in the lifting structure is described as

(11)$$\lsqb y_0 \quad y_1\rsqb ^{T} = {\bf P}_t \cdot {\bf G}\lpar \phi_{i\comma t}\rpar \cdot \lsqb x_0 \quad x_1\rsqb ^{T}\comma \;$$

where

(12)$${\bf G}\lpar \varphi_{i\comma t}\rpar = \left[\matrix{ 1 & 0\cr {f_2 } & 1 }\right]\left[\matrix{ 1 & {f_1}\cr 0 & 1 } \right]\left[\matrix{ 1 & 0\cr f_0 & 1 }\right]$$

and

(13)$$\lsqb f_0 \, \, f_1 \, \, f_2\rsqb = \left[\tan {\phi_{i\comma t} \over 2}\comma \; \, \, {-\sin \phi_{i\comma t}}\comma \; \, \, \tan {\phi_{i\comma t}\over 2}\right].$$

Fig. 6. Implementation of the rotation F.

It indicates that the absolute value of f ₀ and f ₂ becomes extremely huge when the rotation angle φ_{i, t} is close to π [rad]. It magnifies rounding error, and therefore it degrades coding efficiency. This singular point problem is avoided by selecting a proper t ∈ {1, 2, 3, 4}. From (7) and (9),

(14)$$\varphi_{i\comma t} =\theta_i -\Theta_t$$

for

(15)$$\lsqb \Theta_1\ \Theta_2\ \Theta_3\ \Theta_4\rsqb = \lsqb 0 \, \, -\pi /2 \, \, \pi \, \, \pi/2\rsqb$$

is derived. Therefore, a proper t is the one which maximizes the distance between the singular point φ_i,t=π rad and the compensated angle θ_i–Θ_t. In this paper, we follow this procedure reported in [Reference Iwahashi and Kiya16].

C) Rational of the simplification

Firstly, the angle θ₀ in the existing method II is fixed to π/4 rad in our approach. This is because variance of G ₀ and that of G ₁ are almost the same. For example, the optimum angle θ_opt of a rotation F(θ) is given as

(16)$$\theta_{opt} = {{1} \over {2}}\left({{\pi} \over {2} }+ \arctan{ {r-1} \over {2\rho}}\right)\comma \;$$

where

(17)$$\left\{\matrix{ r = {{\sigma_{x_1}^2} \over {\sigma_{x_0}^2}}\comma \; \hfill & \rho = {{E\lsqb \lpar x_0 -\bar{x}_0\rpar \lpar x_1 -\bar{x}_1\rpar \rsqb } \over {\sigma_{x_0}^2}}\comma \; \hfill \cr \bar{x}_i =E\lsqb x_i \rsqb \comma \; \hfill & \sigma_{x_i}^2 =E\lsqb \lpar x_i -\bar{x}_i\rpar ^2\rsqb \comma \; \quad i\in \lcub 0\comma \; 1\rcub }\right.$$

for input signals x = [x ₀x ₁].

Equation (16) is derived as follows. Covariance matrix R_X of X = [x ₀x ₁] is converted to R_Y by the rotation F(θ) as

(18)$${\bf R}_{{\bf Y}} = {\bf F}^{\comma T}\lpar \theta\rpar \cdot {\bf R}_{{\bf X}} \cdot {\bf F}\lpar \theta\rpar.$$

Substituting

(19)$${\bf F}\lpar \theta\rpar = \left[\matrix{ {\cos \theta} & {-\sin \theta}\cr {\sin \theta} & {\cos \theta} }\right]\comma \; \quad {\bf R}_{{\bf X}} = \left[\matrix{ 1 & \rho\cr \rho & r }\right]\sigma_{x_0}^2\comma \;$$

we have

(20)$${{\bf R}_{{\bf Y}} \over {\sigma_{x_0}^2}} = {\bf I}_2 + \left[\matrix{ {1-c} & s\cr s & {1+c} }\right]\cdot {{r-1} \over {2}} + \left[\matrix{ s & c\cr c & {-s} }\right]\cdot \rho$$

for [c s] = [cos 2θ sin 2θ] and I₂ = [1 Reference Matsuda, Kaneko, Minezawa and Itoh1]. Therefore, when it becomes uncorrelated,

(21)$$s\cdot {{r-1} \over {2}}+c\cdot \rho =0$$

holds. Since it is equivalent to

(22)$$\sqrt{\left({{r-1} \over {2}}\right)^2+\rho ^2} \cdot \cos \left(2\theta -\hbox{atctan} {{r-1} \over {2\rho}}\right)= 0\comma \;$$

we have the equation equivalent to (16) as

(23)$$2\theta_{opt} - \hbox{atctan} {{r-1} \over {2\rho}} = {{\pi}\over {2}}.$$

According to (16) and (17), θ_opt=π/4 [rad] for r = 1 in which variance of x ₀ and that of x ₁ is the same. This is experimentally endorsed for several images as summarized in Fig. 7. Note that the tested images are subsampled to generate Bayer-type four color components and used for evaluation. The angle varies depending on input images. Standard deviation was observed to be 6.58° for [x ₀x ₁] = [R B] at maximum. It is 0.20° for [x ₀x ₁] = [G ₀G ₁], namely the angle can be fixed.

Fig. 7. Optimum angle of the rotation F(θ) for images.

Secondly, both of the angles θ₁ and θ₂ in the existing method II are set to zero. This is because variance of Y ₀ in Fig. 5 is close to zero. It comes from a fact that correlation between G ₀ and G ₁ is close to one. For example, variance of Y ₀ is given as

(24)$$\sigma_{Y_0}^2 = {{\sigma_{x_0}^2} \over {2}} \lpar 1+r-2\rho\rpar$$

for x = [x ₀x ₁]. This equation is derived as follows. Since Y ₀ is the first row of the matrix

(25)$${\bf Y} = {\bf F}^{T}\lpar \pi /4\rpar \cdot {\bf X}\comma \;$$

its variance is given as

(26)$$\sigma_{Y_0}^2 = E\left[\left({{x_0 -x_1} \over {\sqrt{2}}}\right)^2\right]= {{E\lsqb x_0^2\rsqb +E\lsqb x_1^2\rsqb -2E\lsqb x_0 x_1 \rsqb } \over {2}}\comma \;$$

which is equivalent to (24).

As described above, it becomes almost zero for r = 1 and ρ = 1. It is no use to rotate this component with other component. Therefore, rotations F(θ₁) and F(θ₂) can be eliminated. Namely, we simplified the existing method II utilizing a fact that statistical relation between the two green components G ₀ and G ₁ is fixed for various input images. The procedure itself is very simple. In the next section, we experimentally investigate its coding performance based on (1) rounding error, (2) adaptability, and (3) compaction rate.

A) Rounding error

Figure 8 illustrates evaluation results of rounding errors included in Y ₀, Y ₁, Y ₂, and Y ₃. The error is defined as difference between output of an RCT with rounding function and that of the same RCT without rounding. Total amount of variance of the error is measured in peak signal to noise ratio (PSNR) defined as

(27)$${PSNR} = 10 \hbox{log}_{10} 255^2-10 \hbox{log}_{10} \sum_{i=0}^3 {\sigma_i^2}\quad \lsqb \hbox{dB}\rsqb \comma \;$$

where σ_i denotes the standard deviation of the error in Y _i. PSNR of the proposed method was observed to be 54.1 dB in average for images. It is greater than the existing method II by 2.2 dB. Reduction of the rounding functions is considered to contribute to this improvement.

Fig. 8. Evaluation of rounding error.

B) Adaptability and compaction rate

Table 1 summarizes evaluation results of data compaction rate for KODAK test images. Since this data set was reported to have considerably high correlations [Reference Zhang, Wu, Yang, Zhang and Zhang19], we also used SIDBA data set and RAW images taken with Nikon D80 camera. For these images, results are summarized in Tables 2 and 3, respectively.

Table 1. The average code length of images in [bpp] for KODAK image set.

Table 2. The average code length of images in [bpp] for SIDBA image set.

Table 3. The average code length of images in [bpp] for RAW images taken with Nikon D80 camera.

The compaction rate was measured with the average code length defined as the ratio of total amount of compressed binary data and the total number of pixels in the image. Note that it indicates 8 bit per pixel per component [bpp], when the data volume is not compressed at all. “Tint” indicates parameter of color balancing in Paint Shop Pro X3. Figures 9–11 illustrate examples of sample images at different value of “Tint”.

Fig. 9. Sample image “Kod05” at different value of “Tint”.

Fig. 10. Sample image “Airplane” at different “Tint”.

Fig. 11. Sample image “Scenery” at different “Tint”.

Table 4 indicates the average code length for all the images in Tables 1–3. In case of “Tint 0”, the existing method II and the proposed method are the same. Those are superior to the existing method I by 0.05 bpp. However, there is no significant difference. On the contrary, in case of “Tint +100”, the proposed method is the best, and the existing method I is the worst. The difference was observed to be 0.20 bpp in average.

Table 4. The average code length of all the images in [bpp].

Figures 12–14 illustrate the average code length of the existing method I and the proposed method for a sample image in Figs 9–11, respectively. These are subtracted by that of the existing method II. There is no significant difference between the existing method II and the proposed method. It means that the proposed method keeps adaptability of the existing method II, even though signal processing is simplified. It also indicates that the existing method I is not adaptive. This is because its multiplier coefficients are fixed as indicated in (3).

Fig. 12. The average code length subtracted by the existing method II for a sample image “Kod05”.

Fig. 13. The average code length subtracted by the existing method II for a sample image “airplane”.

Fig. 14. The average code length subtracted by the existing method II for a sample image “scenery”.

C) Discussions on improvement of compaction rate

As a result of our experiments, it was observed that the proposed method attains almost the same coding performance as the existing method II, even though its computational load is reduced. Remarkably, the proposed method is slightly better than the existing method II for greenish images with “Tint +100”. In the rest of this paper, we endorse it with theoretical analysis on entropy rate and probability density function (PDF) of signals inside the RCT.

Figure 15 illustrates the first-order entropy H, which is defined as

(28)$$H={{1} \over {4}}\sum_{i=0}^3 \sum_j W_j \cdot H_{i\comma j}\comma \;$$

where

(29)$$H_{i\comma j} =-\sum_{x_{i\comma j} } P\lpar x_{i\comma j} \rpar \hbox{log}_2 P\lpar x_{i\comma j}\rpar$$

and x _{i, j} denotes the pixel value in jth frequency band of the RWT in ith component. W _j denotes the ratio of the total number of pixels in jth band to all the number of pixels in each component. The entropy H is determined by the PDF of x _{i, j}, and it estimates the average code length illustrated in Fig. 12.

Fig. 15. The first-order entropy H subtracted by that of the existing method II for “for a sample image “Kod05”.

Figures 16 and 17 illustrate the sigma term H _σ and the epsilon term H_ε of the first-order entropy H, respectively. Those are defined as

(30)$$H = H_\sigma +H_\varepsilon$$

where

(31)$$\left\{\matrix{ H_\sigma =-\displaystyle{{1}\over{2}}\hbox{log}_2 \prod \limits_{i=0}^3 \prod_j \lpar \sigma_{i\comma j}^2\rpar ^{W_j /4} \cr H_\varepsilon =-\displaystyle{{1} \over {2}}\log_2 \prod \limits_{i=0}^3 \prod_j \lpar \varepsilon_{i\comma j}^2\rpar ^{W_j /4} }\right.$$

for standard deviation σ_{i, j} of x _{i, j}. A parameter ɛ_i,j is determined by shape of the PDF [Reference Jayant and Noll20]. Figure 16 indicates that there is no difference between the existing method II and the proposed method in respect of the coding gain, which is often used for theoretical evaluation of lossless coding assuming that shape of all the PDF are the same. On the contrary, Fig. 17 indicates superiority of the proposed method to the existing method II. As a result, it was found that superiority of the proposed method comes from suitability of shape of PDF.

Fig. 16. The sigma term H _σ in the first-order entropy H. It is subtracted by that of the existing method II.

Fig. 17. The epsilon term H _ɛ in the first-order entropy H, which is subtracted by that of the existing method II. (a) Existing method II (KLD = 0.09) (b) Proposed method (KLD = 0.33).

Finally, we compare the existing method II and the proposed method in respect of shape of PDF. In Fig. 18(a), “Existing II” indicates histogram of signal values of Y ₀ in Fig. 4, and “Gaussian” indicates the Gaussian function:

(32)$$G\lpar x\rpar = {{1} \over {\sqrt{2\pi} \sigma}} \exp \left\{-\left({{\lpar x-\mu\rpar ^2} \over {2\sigma^2}}\right)\right\}$$

with zero mean μ = 0 and the same variance σ as the signal values of Y ₀. Similarly, “Proposed” in Fig. 18(b) indicates histogram of signal values of Y ₀ in Fig. 5. Comparing these figures, it is observed that “existing II” is closer to “Gaussian” than “Proposed”. This is numerically endorsed with the Kullback–Leibler divergence defined by

(33)$${KLD}\lpar P\Vert Q\rpar = \sum_k P\lpar k\rpar \log {{P\lpar k\rpar } \over {Q\lpar k\rpar }}\comma \;$$

where k denotes the pixel value. In this case, Q(k) is set to Gaussian function with the same variance as that of Y ₀. It is observed that KLD of the existing method II is 0.09, which is smaller than that of the proposed method by 0.24.

Fig. 18. PDF and KLD.

It can be considered that PDF tends to become Gaussian-like shape according to the central limit theorem in the existing method II, due to more lifting steps than the proposed method. In conclusion, it was indicated that simplification of the proposal in this paper (reduction of the number of lifting steps) keeps “good” shape of PDF, and contributes to improving coding performance. However, there is some room for further investigation.