A Cubic Convolution Interpolation-Based Chroma Subsampling Method for Bayer and RGBW CFA Raw Images

The Bayer and RGBW color filter array (CFA) raw images, denoted by <inline-formula> <tex-math notation="LaTeX">$I^{Bayer}$ </tex-math></inline-formula> and <inline-formula> <tex-math notation="LaTeX">$I^{RGB}$ </tex-math></inline-formula>, respectively, have been widely used in consumer markets. In the demosaicking-first compression scheme, chroma subsampling is necessary prior to compressing <inline-formula> <tex-math notation="LaTeX">$I^{RGB}$ </tex-math></inline-formula> and <inline-formula> <tex-math notation="LaTeX">$I^{RGBW}$ </tex-math></inline-formula>. Several linear interpolation-based chroma subsampling methods have been developed for <inline-formula> <tex-math notation="LaTeX">$I^{Bayer}$ </tex-math></inline-formula>, but no nonlinear interpolation-based chroma subsampling methods have targeted the above two CFA image types simultaneously. In this paper, we first propose a nonlinear interpolation-based, namely the cubic convolution interpolation-based (CCI-based), <inline-formula> <tex-math notation="LaTeX">$2\times 2$ </tex-math></inline-formula> block-distortion function for each <inline-formula> <tex-math notation="LaTeX">$2\times 2$ </tex-math></inline-formula> CFA block <inline-formula> <tex-math notation="LaTeX">$B^{CbCr}$ </tex-math></inline-formula>. Next, using the Cauchy-Schwarz inequality, we prove that the proposed block-distortion function is a convex function in the real domain, which further serves as the base of the initial subsampled chroma solution. Then, a CCI-based iterative method is proposed to improve the initial subsampled chroma solution. The results of comprehensive experimental tests using Bayer and RGBW CFA images created from the three RGB full-color datasets, namely IMAX, SCI (screen content images), and CI (classical images), demonstrate that on the Versatile Video Coding (VVC) platform VTM-11.0, the proposed method achieves substantial quality enhancement and quality-bitrate tradeoff merits of the reconstructed Bayer and RGBW CFA images compared with existing chroma subsampling methods.


I. INTRODUCTION
To save hardware cost, most color digital cameras employ a single-sensor technology with a Bayer color filter array (CFA) to capture real-world scenes. In addition, professional photographers and designers tend to prefer to work with raw Bayer CFA images directly because they can maximize control over the post-processing [22], [32]. At the same time, compressing raw Bayer CFA images also has the potential to improve compression performance relative to working with RGB full-color images.
Accordingly, several compression-first schemes were developed to compress raw Bayer CFA images directly, and these schemes include the geometric rotation approach [23], the integer-reversible spectral-spatial color domain The associate editor coordinating the review of this manuscript and approving it for publication was Wen Chen . approach [25], and the structure conversion approach [7], [10]. However, the intrinsic difficulty in converting CFA structures to the working domain limit the compression performance. Therefore, a demosaicking-first compression scheme [18] was proposed for compressing Bayer CFA images such that the quality of the reconstructed Bayer CFA images has better performance.
In this study, not only the Bayer CFA image I Bayer , we also consider the RGBW CFA image I RGBW . For the former, the four 2 × 2 Bayer CFA patterns [1], as depicted in Figs. 1(a)-(d), have been widely used in modern color digital cameras. Each pixel in I Bayer has only one R (red), G (green), or B (blue) color value. I Bayer contains 25% R, 50% G, and 25% B color values. However, in the low illumination condition, the thermal noise side effect of I Bayer results in low quality of the reconstructed RGB full-color image [16]. To receive more luminance, three RGBW CFA VOLUME 10, 2022 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ patterns, namely the RGBW-Sony CFA pattern [33], the RGBW-Kodak CFA pattern [5], and Yamagami et al.'s CFA pattern [37], were proposed and are depicted in Figs. 2(a)-(c), respectively. Each pixel in I RGBW has only one R, G, B, or W (white) color value. I RGBW contains 12.5% R, 25% G, 12.5% B, and 50% W color values.  The RGBW-Sony pattern [33]. (b) The RGBW-Kodak pattern [5]. (c) Yamagami et al.'s pattern [37].
In the demosaicking-first chroma subsampling scheme, as depicted in Fig. 3, both I Bayer and I RGBW should be first demosaicked to RGB full-color images at the server side, and then the demosaicked RGB full-color images are converted to YCbCr images which serve as the base of chroma subsampling. To demosaic I Bayer , the methods in [15], [26], [39], [41] can be used; to demosaic I RGBW , the methods in [6], [29] can be used. Here, we take Kiku et al.'s code [15] to demosaic I Bayer and take Condat's code [6] to demosaic I RGBW . After that, the demosaicked RGB full-color image I RGB is converted to a YCbCr image I YCbCr by using the BT.709-5 color conversion [13]: Then, for each 2 × 2 CbCr block B CbCr in I CbCr , one 4:2:0 subsampling method is used to obtain the subsampled (Cb, Cr)-pair of B CbCr . Accordingly, the subsampled CbCr image I sub,CbCr and the luma image I Y are fed into the encoder for compression. The lack of correlation between I Y and I sub,CbCr greatly contributes to the effectiveness of the image compression. After encoding I sub,CbCr and I Y , the encoded bitstream is transmitted to the decoder via the internet.
At the client side in Fig. 3, the decoded subsampled CbCr image I sub,CbCr is upsampled by the selected chroma upsampling process. Based on the 2 × 2 Bayer CFA pattern, the decompressed CbCr image and luma image can be directly transformed to the reconstructed Bayer CFA image I rec,Bayer using the following YCbCr-to-RGB conversion: Based on the 4 × 4 RGBW CFA pattern, the R, G, and B color pixels in the reconstructed RGBW CFA image I rec,RGBW can be obtained directly using Eq. (2), while the W color pixels in I rec,RGBW can be obtained by averaging the reconstructed R, G, and B color values at the same location.
Zhang et al. [43] proposed an IDID method that prefers the chroma upsampling process NEDI (new edge-directed interpolation) [17] used at the client side. To improve the IDID method, considering the palette mode used for screen content images [20], Wang et al. [36] proposed a JCDU method. The JCDU method prefers the chroma upsampling process BICU (bicubic interpolation) [27], [28] used at the client side.
2) FOUR BAYER CFA PATTERN-BASED CHROMA SUBSAMPLING METHODS Different from the eight chroma subsampling methods mentioned above Chen et al. [3] proposed a direct mapping (DM) method when the input image is I Bayer based on the following observation: in Eq. (2), the R value is dominated by the Y and Cr values, and the B value is dominated by the Y and Cb values. Therefore, the subsampled chroma pair of B CbCr equals (Cb 3 , Cr 2 ). Lin et al. [18] proposed a COPY-based 2 × 2 Bayer CFA block-distortion function in which the operator COPY is a special linear interpolation used to copy the subsampled chroma-pair parameter of B CbCr , denoted by (Cb s , Cr s ), as the four estimated chroma pairs of B CbCr . Then, based on the COPY-based block-distortion function, they applied the differentiation (DI) technique to determine the solution of (Cb s , Cr s ) as the subsampled chroma pair of B CbCr . The DI method prefers the chroma upsampling process 'COPY' used at the client side. In fact, the upsampling process 'COPY' equals the process 'NN (nearest neighbor)' supported by the VVC (Versatile Video Coding) platform VTM-11.0 [34].
To improve the DI method, Chung et al. [4] proposed a COPY-and gradient-descent (GD)-based chroma subsampling method. In the GD method, the block-distortion function is still based on the COPY-based approach, but they utilized the value verification technique to prove that their block-distortion function is a convex function, leading to the design of the GD-and BILI-based (bilinear interpolationbased) iteration method to improve the DI method. The GD method prefers the chroma upsampling process BILI [27], [28] used at the client side. Based on the DI method [18] but considering the demosaicked RGB full-color blockdistortion, Lin et al. [19] proposed a 'modified 4:2:0(A)' chroma subsampling method that selects the best case among the four subsampled chroma pairs of B CbCr by considering the ceiling operation-based 4:2:0(A) and the floor operationbased 4:2:0(A). The 'modified 4:2:0(A)' method prefers the chroma upsampling process TN (three neighboring reference pixel-based method) [40] used at the client side.

B. MOTIVATION
To the best of our knowledge, no chroma subsampling methods targeting I Bayer can also be used for I RGBW , and no chroma subsampling methods targeting only I RGBW have ever been proposed. In addition, for I Bayer , the block-distortion function used in the state-of-the-art chroma subsampling methods [4], [18], [19] is based on a special linear interpolation approach, i.e., the COPY-based approach, limiting the quality performance.
The above reasons motivated us to design a new and more effective nonlinear interpolation-based method, namely the cubic convolution interpolation-based (CCI-based) chroma subsampling method, to handle chroma subsampling for I Bayer and I RGBW simultaneously, achieving substantial quality enhancement and quality-bitrate tradeoff improvements of the reconstructed Bayer and RGBW CFA images.

C. CONTRIBUTIONS
The four major contributions of the proposed new CCI-based chroma subsampling method for I Bayer and I RGBW related to the previously published papers [3], [4], [18], [19] for only I Bayer are clarified as follows.
(1) In the previously published paper [3], based on the observation from the YCbCr-to-RGB color transformation, the subsampled chroma pair of B CbCr equals (Cb 3 , Cr 2 ). Instead of the COPY-based block-distortion functions used for only I Bayer [4], [18], [19], in this paper, a new and more effective nonlinear interpolation-based, namely the cubic convolution interpolation-based (CCI-based), blockdistortion function is proposed for I Bayer and I RGBW simultaneously.
(2) In [4], using the value verification approach for only the Bayer CFA pattern in Fig. 1(a), i.e., pat 1 , the determinant of the Hessian matrix of the COPY-based block-distortion function is 1.5577, proving the convex property of their blockdistortion function. In this paper, using the Cauchy-Schwarz inequality, we provide solid proof that our CCI-based blockdistortion function is a convex function for all ten considered Bayer and RGBW CFA patterns, i.e., pat i for 1 ≤ i ≤ 10. In particular, for the Bayer CFA pattern 'pat 1 ', the determinant of the Hessian matrix of our CCI-based block-distortion function is 22.0337 versus 1.5577 [4].
(3) Using the convex property of our CCI-based blockdistortion function, we derive a formula to determine the initial subsampled chroma solution of B CbCr . Furthermore, we propose an iterative method to refine the initial subsampled chroma solution for I Bayer and I RGBW simultaneously. In particular, for I Bayer , the experimental results have demonstrated that when compared with the two combinations, DI-COPY [18] and GD-BILI [4], the PSNR (peak signal-to-noise ratio) gains of our combination are 2.4841 dB and 2.0531 dB, respectively, when setting QP (quantization parameter) to zero.
(4) Based on the ground truth Bayer and RGBW CFA images which are created from the three RGB fullcolor datasets, namely IMAX [12], SCI (screen content images) [30], and CI (classical images) [44], the comprehensive experimental results show that on VTM-11.0 [34], our method achieves the best PSNR, SSIM (structural similarity index measure) [35], FSIM (feature similarity index measure) [42], visual effect, and BD-PSNR (Bjøntegaard delta PSNR difference) [2] performance of the reconstructed Bayer and RGBW CFA images when compared with the six traditional methods and the six state-of-the-art methods [3], [4], [18], [19], [36], [43]. VOLUME 10, 2022 The remainder of this paper is organized as follows. In Section II, the proposed CCI-based block-distortion function is presented. In Section III, we prove the convex property of the proposed block-distortion function, and then the proposed iterative chroma subsampling method is presented. In Section IV, the experimental results are provided to justify the quality and quality-bitrate tradeoff merits of our method. In Section V, some concluding remarks are offered.

II. THE PROPOSED CCI-BASED BLOCK-DISTORTION FUNCTION
Different from the COPY-based block-distortion functions used in the previous works [4], [18], [19] only for I Bayer , we propose a new and more effective CCI-based 2 × 2 blockdistortion function for I Bayer and I RGBW simultaneously.
Because the four Bayer CFA patterns in Fig. 1, denoted by pat i for 1 ≤ i ≤ 4, are all 2 × 2 blocks, to handle the 2 × 2 Bayer and 4×4 RGBW CFA blocks simultaneously, the three 4 × 4 RGBW CFA patterns in Fig. 2 are partitioned into six basic 2×2 RGBW CFA patterns, namely pat i for 5 ≤ i ≤ 10, in Figs. 4(a)-(f). According to the row-major scanning order, the three 4 × 4 RGBW CFA patterns in Fig. 2 are expressed as the following three ordered sets:

A. THE CCI-BASED METHOD TO ESTIMATE THE FOUR CHROMA PAIRS OF B CbCr
Before presenting the proposed CCI-based block-distortion function, we first propose a more effective CCI-based estimation method to estimate the four chroma pairs of B CbCr compared with the COPY-based estimation used in [4], [18], [19]. For simplicity, we only describe how to estimate the four Cb entries of B Cb . The estimation of the four Cr entries of B Cr can be easily followed. For easy exposition, only the top-left Cb component of B Cb , namely Cb est 1 in Fig. 5, is estimated, and the other three estimated Cb components, namely Cb est 2 , Cb est 3 , and Cb est 4 , can be obtained using the same method. Let each 2×2 Cb block in Fig. 5 be viewed as a 1×1 macropixel. Initially, we put the origin of the x-y coordinate system at the location of Cb est 1 . Therefore, Cb est 1 is located at (0, 0), and the locations of the subsampled Cb parameter Cb s and the neighboring reference subsampled Cb values are illustrated in Fig. 5, where the subsampled Cb parameter Cb s (= Cb 2,2 ) is located at (0.25, -0.25), Cb 1,2 is located at (0.25, 0.75), and so on. We deploy the cubic convolution interpolation, which was originally proposed by Keys [14], in the estimation of Cb est 1 . The value of the weight function assigned to each reference point Cb i,j , as depicted by one black bullet with yellow background in Fig. 5, is dependent on the vertical distance d v i,j and the horizontal distance d h i,j relative to the location of Cb est 1 . As mentioned before, Cb est 1 is located at (0, 0). For example, for the reference point Cb 2,2 (= Cb s ) in Fig. 5, The weight function W(d) is expressed as the following cubic function: Due to the constraint |d|<2, as depicted in Fig. 5, only 16 reference points in the 4 × 4 sliding window are needed. Among the 16 reference subsampled Cb values, the former 10 reference subsampled Cb values are obtained in advance by using the proposed chroma subsampling method, which will be formally presented in Section III. Except for the subsampled Cb parameter Cb s , the five future subsampled Cb values can be obtained by using one available subsampling method, e.g., 4:2:0(A). The value of the weight function W (i, j) assigned to Cb i,j is thus expressed aŝ Consequently, the estimation of Cb est 1 is expressed as In Eq. (6), the estimation of Cb est 1 is expressed as the sum of 16 products.
Furthermore, as depicted in Fig. 6(a), we first put the origin of the x-y coordinate system at the location of Cb est 2 , and then we move the 4 × 4 sliding window one pixel to the right to estimate Cb est 2 . By the same argument, we move the 4 × 4 sliding window one pixel downward to estimate Cb est 4 , and then we move the 4 × 4 sliding window one pixel to the left to estimate Cb est 3 . It is notable that the above cubic convolution interpolation-based computation method for estimating the four entries of B CbCr can explain why our method is called the cubic convolution interpolation-based, i.e., the CCI-based, method.

B. THE PROPOSED CCI-BASED CFA BLOCK-DISTORTION FUNCTION
After presenting the proposed CCI-based estimation method for estimating the four chroma pairs of B CbCr , the estimated 2 × 2 CbCr block is denoted by B est,CbCr . Based on B est,CbCr , the collocated 2×2 luma block B Y , and the 2×2 CFA pattern pat i , which corresponds to the input ground truth 2 × 2 CFA raw block B pat i , by Eq. (2), the estimated 2 × 2 CFA block B est,pat i can be obtained at the server side. We now present the proposed CCI-based 2 × 2 CFA block-distortion function to measure the sum of squared errors (SSE) between each 2×2 ground truth CFA block B pat i and the estimated 2 × 2 CFA block B est,pat i for 1 ≤ i ≤ 10. For easy exposition, we take the 2 × 2 Bayer CFA pattern pat 1 in Fig. 1(a) as the representative for I Bayer . Accordingly, our block-distortion between B pat 1 (= (G 1 , R 2 , B 3 , G 4 )) and B est,pat 1 where ||.|| denotes a 2-norm operation.
Using the YCbCr-to-RGB transformation in Eq. (2), our block-distortion function in Eq. (7) can be rewritten as Furthermore, we take the 2 × 2 RGBW CFA pattern pat 6 as the representative for I RGBW . Similarly, our block-distortion function to measure the SSE between B pat 6 (= (W 1 , G 2 , G 3 , W 4 )) and B est,pat 6 Using the YCbCr-to-RGB transformation in Eq. (2), Eq. (10) can be rewritten as From the explicit form of the block-distortion D pat 1 (Cb s , Cr s ) in Eqs. (8)- (9) for I Bayer , we deduce that the values of α k and β k , 1 ≤ k ≤ 4, are dependent on the color of the kth entry of B pat 1 . Similarly, from the explicit form of the block-distortion D pat 6 (Cb s , Cr s ) in Eqs. (11)-(12), we deduce that the values of α k and β k , 1 ≤ k ≤ 4, are also dependent on the color of the kth entry of B pat 6 . Let pat i,k denote the kth entry of pat i , 1 ≤ k ≤ 4 and 1 ≤ i ≤ 10. For example, for pat 1 = <pat 1,1 , pat 1,2 , pat 1,3 , pat 1,4 >, we have pat 1,1 = G, pat 1,2 = R, pat 1,3 = B, and pat 1,4 = G. Table 1 lists the values of α pat i,k and β pat i,k for 1 ≤ k ≤ 4 and 1 ≤ i ≤ 10.
Consequently, our general 2 × 2 CFA block-distortion function, which is used to measure the SSE between B pat i and B est,pat i , 1 ≤ i ≤ 10, is expressed as

III. THE PROPOSED CCI-BASED ITERATIVE CHROMA SUBSAMPLING METHOD
For the considered Bayer and RGBW CFA image types, we first provide solid proof to show that our 2×2 CFA blockdistortion function in Eq. (13) is a convex function in the real domain, which further serves as the base of the initial subsampled chroma solution of B CbCr . Next, a formula is derived to calculate the initial subsampled chroma solution. Furthermore, we propose an iterative method to refine the initial subsampled chroma solution.
Compared with the COPY-based block-distortion function [4], [18], [19] used only for the 2 × 2 Bayer CFA pattern pat 1 in Fig. 1(a), the proposed CCI-based block-distortion function for pat i , 1 ≤ i ≤ 10, is more effective and general. It is notable that for pat 1 in Fig. 1(a), the determinant of the Hessian matrix of our CCI-based block-distortion function is 22.0337, while the determinant of the Hessian matrix of the previous COPY-based block-distortion function is 1.5577.

B. DETERMINING THE INITIAL SUBSAMPLED CHROMA SOLUTION OF B CbCr
Proposition 1 implies that in the real domain, a critical point exists in the convex block-distortion function in Eq. (13) at which the convex block-distortion function has a minimal value. However, practically, we consider the nonnegative integer domain for chroma subsampling. Therefore, the critical point of Eq. (13) is used as the initial subsampled chroma solution in our iterative method.
To determine the critical point of Eq. (13), we take the first derivative on Eq. (13) with the subsampled chroma parameters Cb s and Cr s , respectively. Setting the two derivatives to zero yields Furthermore, the explicit form of Eq.  (21), as shown at the bottom of the next page, and is used as our initial subsampled chroma solution of B CbCr for pat i , 1 ≤ i ≤ 10. It is notable that for an RGBW CFA image I RGBW with the RGBW-Sony CFA pattern (= < pat 5 , pat 6 , pat 6 , pat 5 >), the two initial subsampled chroma solutions of B CbCr for pat 5 and pat 6 can be obtained by calling Eq. (21) with i = 5 and i = 6, respectively.

C. THE PROPOSED ITERATIVE CHROMA SUBSAMPLING METHOD
We take a 4 × 4 RGBW-Kodak CFA block example in Fig. 7(a) to sketch the idea of the proposed iterative chroma subsampling method. The four 2 × 2 partitioned CFA blocks of Fig. 7(a) are pat 7 , pat 6 , pat 6 , and pat 8 .
For simplicity, we only take the top-right 2×2 RGBW CFA block of Fig. 7(b) to sketch the idea of our iterative method. The demosaicked 2 × 2 RGB full-color block of the topright 2 × 2 RGBW CFA block of Fig. 7(b) is demonstrated in Fig. 7(c). By Eq. (1), the converted 2 × 2 YCbCr block of Fig. 7(c) is demonstrated in Fig. 7(d). According to the neighboring reference subsampled chroma pairs in Fig. 7(e), the grid plot of the block-distortion function D pat 6 (Cb s , Cr s ) is depicted in Fig. 7(f). By Eq. (21), the initial subsampled chroma solution of B CbCr is depicted by a black circle in Fig. 7(f). Clearly, there is room for improvement from the initial subsampled chroma solution to the desirable subsampled chroma solution, as depicted by the red circle.
To refine the initial subsampled chroma solution, we first consider the eight neighboring distance-1 reference points of the initial solution, namely (Cb  Fig. 8(a). Then, we select the best neighboring point with the smallest block-distortion.
If the block-distortion of one neighboring distance-1 reference point is smaller than that obtained in the previous iteration, it replaces the old subsampled chroma solution, and we go to the next solution-refinement process. Otherwise, to avoid getting stuck in a local minimum, we try an uphill-move step in the current iteration by considering the 16 neighboring distance-2 reference points, which are marked in blue in Fig. 8(b). We repeat the above solutionrefinement process until it reaches the termination condition. For this case in Fig. 8, the current subsampled (Cb, Cr)-pair solution (111, 181) with block-distortion 11 is replaced by a better subsampled (Cb, Cr)-pair solution (112, 183) with block-distortion 9.
In our experience, the uphill-move step in the proposed iterative method can achieve average PSNR gains of 0.04 dB, 0.4 dB, and 0.2 dB for one Bayer CFA image, RGBW-Kodak CFA image, and RGBW-Sony CFA image, respectively. Empirically, for each B CbCr , the average numbers of distance-1 iterations required in the proposed method for the Bayer, RGBW-Sony, and RGBW-Kodak CFA patterns are 1.57, 1.69, and 1.78, respectively. The average numbers of distance-2 iterations required in the proposed method for the three considered CFA patterns are 1.03, 1.1, and 1.15.
The proposed CCI-based iterative chroma subsampling method for each 2×2 CbCr block B CbCr with the CFA pattern pat i for 1 ≤ i ≤ 10, denoted by B CbCr,pat i , is listed in Algorithm 1.
The execution code of our combination 'Proposed-CCI' can be accessed from the website [9], where 'Proposed' denotes our iterative chroma subsampling method used at the server side and 'CCI' denotes the chroma upsampling process used at the client side. It is notable that after receiving the decompressed subsampled CbCr image, we first construct an initial upsampled CbCr image. Each time, we consider sixteen subsampled chroma pairs within a 4 × 4 window, which are marked by the 16 black bullets in Fig. 6, and then the chroma upsampling process 'CCI', which is quite similar to the CCI-based method to estimate the four chroma pairs of B CbCr described in Subsection II.A, fills the four missing chroma pixels one by one. The upsampling process 'CCI' continues the filling operations until all missing chroma pixels are constructed.

IV. EXPERIMENTAL RESULTS
Based on the ground truth Bayer and RGBW CFA images which are created from the three RGB full-color datasets, IMAX [12], SCI [30], and CI [44], the quality enhancement and quality-bitrate tradeoff improvement merits of our chroma subsampling method for I Bayer and I RGBW are demonstrated. The execution time comparison is also reported.
The proposed method and all comparative methods are implemented on a computer with an Intel Core i7-8700 CPU 3.2 GHz and 24 GB RAM. The operating system is the Microsoft Windows 10 64-bit operating system.

).
Step 3: s ), we perform m ← m + 1 and go to step 2; otherwise, we go to step 4.
Step 4: We calculate the 16 block-distortion values of the neighboring distance-2 reference points of (Cb Step 5: The program development environment is Visual C++ 2019. The VVC reference software platform used for compression is VTM-11.0.

A. QUALITY ENHANCEMENT MERIT OF OUR METHOD AND EXECUTION TIME COMPARISON
When setting QP to zero, we take the three quality metrics, namely PSNR, SSIM [35], and FSIM [42], to report the quality enhancement merit of our method for I Bayer and I RGBW relative to the existing methods. PSNR is used to evaluate the average quality of one reconstructed CFA image I rec,CFA and it is defined by 10 log 10 255 2 MSE (22) where N denotes the number of the CFA images in one dataset; MSE (mean square error) equals 1 (I CFA − I rec,CFA ) 2 where I CFA denotes a ground truth CFA image and XY denotes the image size. We first calculate the PSNR value of each dataset and then calculate the average PSNR value of the three PSNR values of the three considered datasets.
The quality metric SSIM is expressed as the product of the luminance mean similarity, the contrast similarity, and the structure similarity between the ground truth CFA image and the reconstructed CFA image. The quality metric FSIM utilizes the phase consistency and gradient magnitude to weight the local quality maps to obtain a quality score of the reconstructed image.

1) QUALITY ENHANCEMENT MERIT FOR I Bayer
For simplicity, we take the commercial Bayer CFA pattern in Fig. 1(a) as the representative to demonstrate the quality enhancement merit of the reconstructed Bayer CFA images using our iterative chroma subsampling method, namely 'Proposed'. For fairness, our two combinations, namely 'Proposed-BICU' and 'Proposed-CCI', are included in the comparative combinations. The two combinations 'Proposed-BICU' and 'Proposed-CCI' are implemented for the Bayer, RGBW-Sony, and RGBW-Kodak CFA images, respectively, and in total, there are six variants whose execution codes can be accessed from the website [9].

2) QUALITY ENHANCEMENT MERIT FOR I RGBW
We take the two commercial RGBW CFA patterns, namely the RGBW-Sony CFA pattern in Fig. 2(a) and the RGBW-Kodak CFA pattern in Fig. 2(b), as representatives to show the quality enhancement merit of our combination 'Proposed-CCI'. Table 4 indicates that our combination 'Proposed-CCI' has the highest PSNR, SSIM, and FSIM in boldface among the considered 11 combinations. Excluding the combination 'Proposed-BICU', the overall average PSNR gain of our combination 'Proposed-CCI' over the other nine comparative combinations is 4.033 dB, achieving a significant quality improvement.
For I RGBW with the RGBW-Sony CFA pattern in Fig. 2(a), as tabulated in the parentheses () in Table 4 Table 4, the PSNR gains of our combination 'Proposed-CCI' over the above four comparative combinations are 1.5503 (= 46.8499 -45.2996) dB, 3.6223 dB, 2.3390 dB, and 1.5129 dB, respectively, also indicating significant quality enhancement effects of our combination.

3) COMPUTATIONAL TIME ANALYSIS AND EXECUTION TIME COMPARISON
We first analyze the computational complexities required in the proposed iterative chroma subsampling method for the 2 × 2 Bayer, RGBW-Sony, and RGBW-Kodak CFA patterns. Next, we analyze the computational complexities required in the proposed method and the four state-of-the-art chroma subsampling methods, namely DM [3], DI [18], GD [4], and modified 4:2:0(A) [19], for the 2 × 2 Bayer CFA pattern. Finally, the execution times are compared. As mentioned before, for each 2 × 2 CbCr block B CbCr , the average numbers of distance-1 iterations required in the proposed method for the Bayer, RGBW-Sony, and RGBW-Kodak CFA patterns are 1.57, 1.69, and 1.78, respectively. The average numbers of distance-2 iterations required in the proposed method for the Bayer, RGBW-Sony, and RGBW-Kodak CFA patterns are 1.03, 1.  (13) in [18]) to determine the subsampled chroma pair of B CbCr .
The GD method [4] first takes 38 operations to determine the initial subsampled chroma solution using Eqs. (12)-(13) in [18]. Next, it calculates the distortion of the initial subsampled chroma solution. Furthermore, GD considers the distance-1 iterations. Because our CCI-based blockdistortion function in Eq. (13) is more complicated than that in GD and our method must consider distance-1 and distance-2 iterations, the computational complexity required in GD is less than ours.
For the 2 × 2 Bayer CFA pattern, let the computational complexity required in DM [3], DI [18], GD [4], and our method be denoted by T DM , T DI , T GD , and T Proposed , respectively. From the above analysis, we tabulate the time complexities of the four methods in Table 2, where #(memory accesses) denotes the number of memory accesses and #(BILI-based D-1 iterations), #(CCI-based D-1 iterations), and #(CCI-based D-2 iterations) denote the number of BILI-based distance-1 iterations, the number of CCI-based distance-1 iterations, and the number of CCI-based distance-2 iterations, respectively. From Table 2, we thus have T DM <T DI <T GD <T Proposed . The actual execution time cost required in the considered four methods, as tabulated in the last column of Table 3, justifies the computational complexity order of the four methods.

c: EXECUTION TIME COMPARISON
For I Bayer and I RGBW , the average execution time (in seconds) required in each considered combination is tabulated in the last columns of Tables 3 and 4, respectively. For each combination, the execution time requirement equals the sum of the execution time required in chroma subsampling and that required in chroma upsampling. In particular, when compared with the two combinations 'IDID-NEDI' [43] and 'JCDU-BICU' [36], our combination 'Proposed-CCI' takes much less time and has higher average PSNR, SSIM, and FSIM. Although our combination takes more time than the other combinations, our combination achieves much better quality enhancement effects.

B. BD-PSNR MERIT AND VISUAL EFFECT MERITS OF OUR METHOD
When setting four QP intervals, namely [4,16], [12,24], [20,32], and [28,40], we adopt the BD-PSNR metric [2] to demonstrate the quality-bitrate tradeoff merit of our combination 'Proposed-CCI' relative to the considered combinations. To calculate the BD-PSNR value of one combination, we take the combination '4:2:0(A)-BICU' as the basis of comparison. For one considered combination and the comparative combination, BD-PSNR indicates the average PSNR difference under the same bitrate. In addition, the visual effects of our combination are compared with those of some existing combinations.

1) VISUAL EFFECT MERIT
We take the 4th ground truth RGBW-Sony CFA image created from IMAX as the testing image, as depicted in Fig. 9(a), to demonstrate the visual effect merit of the reconstructed RGBW-Sony CFA image using our combination 'Proposed-CCI' relative to the nine comparative combinations. The magnified subimage in Fig. 9(b) is decoupled from the region containing the leaves in Fig. 9(a). For QP = 40, after performing the ten considered combinations on Fig. 9(a), the ten reconstructed magnified subimages for Fig. 9(b) are shown in Figs. 9(c)-(l), and we observe that our combination 'Proposed-CCI' has the best visual and texture-preserving effects of the ten combinations.

2) BD-PSNR MERIT
For I Bayer with the Bayer CFA pattern in Fig. 1(a), we demonstrate the BD-PSNR merit of the reconstructed Bayer CFA images using our combination 'Proposed-CCI' relative to  the existing combinations. In Table 5, we observe that our combination 'Proposed-CCI' has the best BD-PSNR performance in boldface among the 17 combinations. In particular, our combination has better BD-PSNR performance than the current work 'GD-BILI'.
Considering the RGBW-Sony CFA pattern in Fig. 2(a) and the RGBW-Kodak CFA pattern in Fig. 2(b), we demonstrate the BD-PSNR merit of the reconstructed RGBW CFA images using our combination 'Proposed-CCI' relative to the existing combinations. Here, the BD-PSNR results for I RGBW with the RGBW-Kodak CFA pattern are listed in the parentheses in Table 6. Table 6 indicates that for the three QP intervals, namely [4,16], [12,24], and [20,32], our combination 'Proposed-CCI' has the best BD-PSNR performance in boldface among the considered combinations. However, for the QP interval [28,40], our combination 'Proposed-BICU' has the best BD-PSNR performance, as shown in boldface, while our combination 'Proposed-CCI' is ranked second.

V. CONCLUSION
We have presented our new and effective CCI-based iterative chroma subsampling method for Bayer and RGBW CFA images. All three 4 × 4 RGBW CFA patterns in Fig. 2 are first partitioned into six 2 × 2 RGBW CFA patterns. Next, for I Bayer and I RGBW , a CCI-based block-distortion function is proposed. Then, we provide proof to show that the proposed CCI-based block-distortion function is a convex function that serves as the base of the initial subsampled chroma solution. Furthermore, a closed form is derived to obtain the initial subsampled chroma solution. To refine the initial subsampled chroma solution, an iterative chroma subsampling method is proposed. Finally, comprehensive experimental results demonstrate that on VTM-11.0, the proposed CCI-based chroma subsampling method achieves substantial quality enhancement and BD-PSNR improvement merits of the reconstructed Bayer and RGBW CFA images compared with existing methods. Future work will first deploy the abovementioned chroma subsampling works in the encrypted image for secret sharing [8], [38]. Second, the design of a deep learning-based chroma subsampling method will be investigated in future work. YA-CHI TSENG received the B.S. degree in computer science and engineering from Yuan Ze University, Taoyuan, Taiwan, in 2020. She is currently pursuing the M.S. degree in computer science and engineering with the National Taiwan University of Science and Technology, Taipei, Taiwan. Her research interests include image processing and video compression. VOLUME 10, 2022