Color Filter Arrays for Quanta Image Sensors

TL;DR

This paper introduces an optimization-based framework for designing color filter arrays tailored for quanta image sensors and small pixels, balancing light efficiency, anti-aliasing, and crosstalk.

Contribution

It unifies and generalizes existing design methodologies, enabling simultaneous optimization of multiple performance criteria for color filter arrays.

Findings

The framework outperforms existing methods in experimental comparisons.

It effectively balances light efficiency with anti-aliasing and crosstalk suppression.

The approach is flexible and applicable to various design constraints.

Abstract

Quanta image sensor (QIS) is envisioned to be the next generation image sensor after CCD and CMOS. In this paper, we discuss how to design color filter arrays for QIS and other small pixels. Designing color filter arrays for small pixels is challenging because maximizing the light efficiency while suppressing aliasing and crosstalk are conflicting tasks. We present an optimization-based framework which unifies several mainstream color filter array design methodologies. Our method offers greater generality and flexibility. Compared to existing methods, the new framework can simultaneously handle luminance sensitivity, chrominance sensitivity, cross-talk, anti-aliasing, manufacturability and orthogonality. Extensive experimental comparisons demonstrate the effectiveness of the framework.

Tables

Table 1. TABLE I: CFA Design Criteria

Criterion	Purpose	Regular Pixels				QIS
		[21]	[33]	[24]	[36]	[31]	Ours
Proposition 1	To minimize noise power after linear demosaicking	$\times$	$✓$	$\times$	$\times$	$\times$	$✓$
Proposition 2	To simplify denoising of luminance channel	$✓$	$✓$	$✓$	$✓$	$✓$	$✓$
Proposition 3	To maximize spatial resolution	$✓$	$✓$	$✓$	$\times$	$\times$	$✓$
Proposition 4	To mitigate cross-talk	$\times$	$\times$	$\times$	$✓$	$✓$	$✓$
Definition 4	Enforce total orthogonality	$✓$	$\times$	$✓$	$\times$	$\times$	$✓$
Definition 4	Enforce quadrature orthogonality	$\times$	$✓$	$\times$	$\times$	$\times$	$✓$

Table 2. TABLE II: CFA parameters and Reconstruction quality measured by YSNR and SMI metrics on Macbeth ColorChecker and average CPSNR on DIV2K evaluation dataset. An arrow is placed after each metric to show whether higher or lower is better

		Experiment 1					Experiment 2				Experiment 3
Size	CFA Pattern	CFA Parameters					YPSNR (dB) $\uparrow$		SMI (%) $\uparrow$		CPSNR (dB) $\uparrow$
		${\gamma }_l\uparrow$	${\gamma }_c\uparrow$	TV $\downarrow$	$J_l\downarrow$	$\kappa \left(𝑨\right)$ $\downarrow$	w/o Ctk	w/ Ctk	w/o Ctk	w/ Ctk	w/o Ctk	w/ Ctk
$4\times 4$	RGBCY [31]	0.679	0.125	0.156	6.53	2.376	23.24	22.15	93.65	92.90	30.65	30.63
	RGBCWY [31]	0.597	0.102	0.131	2.67	2.143	23.57	22.45	93.71	93.10	30.75	30.74
	Hao et al. [24]	0.586	0.094	0.206	1.00	2.305	24.13	18.50	94.34	91.31	31.34	30.83
	Ours	0.577	0.121	0.131	2.86	1.795	24.16	23.17	94.07	93.47	31.10	31.12
$3\times 3$	Biay-Cheng et al. [36]	0.612	0.167	0.175	1.65	1.892	23.92	22.04	94.44	93.88	31.07	31.14
	Ours	0.577	0.160	0.175	1.74	1.561	24.38	23.64	94.49	94.36	31.24	31.37
$3\times 2$	Condat [33]	0.866	0.250	0.317	1.16	1.915	25.33	21.68	94.71	93.87	31.30	31.17
	Ours	0.866	0.250	0.212	1.40	1.837	25.28	24.30	94.25	93.77	31.27	31.32
$4\times 2$	Hirakawa-Wolfe [21]	0.866	0.125	0.275	0.98	2.767	25.38	21.22	94.57	93.58	31.42	31.08
	Ours	0.866	0.187	0.256	1.27	2.166	24.89	22.34	94.36	94.03	31.40	31.31
$7\times 7$	Bai et al. [28]	0.577	0.091	0.204	1.87	2.242	24.26	19.68	94.56	92.94	31.41	31.09
	Ours	0.577	0.238	0.204	1.45	1.294	24.86	23.85	94.45	94.35	31.01	31.20

Equations

$${\mathbf{im}}_{rgb}(m,n)=\left[\begin{array}{c}{\mathrm{im}}_r(m,n)\\ {\mathrm{im}}_g(m,n)\\ {\mathrm{im}}_b(m,n)\end{array}\right],$$

$${\mathbf{c}}_{rgb}(m,n)=\left[\begin{array}{c}{c}_r(m,n)\\ c_g(m,n)\\ c_b(m,n)\end{array}\right],$$

$$\begin{array}{rl}\theta (m,n)& =\eta {\mathbf{c}}_{rgb}(m,n{)}^T{\mathbf{im}}_{rgb}(m,n)\\ & =\eta \sum _{i\in \{r,g,b\}}{c}_i(m,n){\text{im}}_i(m,n).\end{array}$$

$$\mathbb{P}(Y_j=y_j)=\frac{{\theta }_m^{y_j}{e}^{-{\theta }_j}}{y_j!}.$$

$$\mathbb{P}(B_j=b_j)={\Psi }_q({\theta }_j{)}^{1-b_j}{\left(1-{\Psi }_q({\theta }_j)\right)}^{b_j},$$

$${\mathbf{h}}_{rgb}(m,n)=\left[\begin{array}{c}{h}_r(m,n)\\ h_g(m,n)\\ h_b(m,n)\end{array}\right],$$

$$h_r=\left[\begin{array}{cc}0& 1\\ 0& 0\end{array}\right],h_g=\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right],h_b=\left[\begin{array}{cc}0& 0\\ 1& 0\end{array}\right],$$

$$\mathbf{T}=\left[\begin{array}{ccc}1/\sqrt{3}& 1/\sqrt{3}& 1/\sqrt{3}\\ -1/\sqrt{6}& 2/\sqrt{6}& -1/\sqrt{6}\\ 1/\sqrt{2}& 0& -1/\sqrt{2}\end{array}\right].$$

$${\mathbf{im}}_{l\alpha \beta }=\left[\begin{array}{c}{\text{im}}_l\\ {\text{im}}_{\alpha }\\ {\text{im}}_{\beta }\end{array}\right]=\mathbf{T}\left[\begin{array}{c}{\text{im}}_r\\ {\text{im}}_g\\ {\text{im}}_b\end{array}\right]=\left[\begin{array}{c}\left({\text{im}}_r+{\text{im}}_g+{\text{im}}_b\right)/\sqrt{3}\\ \left(-{\text{im}}_r+2{\text{im}}_g+{\text{im}}_b\right)/\sqrt{6}\\ \left({\text{im}}_r-{\text{im}}_b\right)/\sqrt{2}\end{array}\right],$$

$$\theta (m,n)$$

$$=\eta {\mathbf{c}}_{l\alpha \beta }(m,n{)}^T{\mathbf{im}}_{l\alpha \beta }(m,n)$$

$$=\eta \sum _{i\in \{l,\alpha ,\beta \}}{c}_i(m,n){\text{im}}_i(m,n),$$

$${\mathbf{c}}_{l\alpha \beta }(m,n)=\mathbf{T}{\mathbf{c}}_{rgb}(m,n).$$

$${\mathbf{h}}_{l\alpha \beta }(m,n)=\left[\begin{array}{c}{h}_l(m,n)\\ h_{\alpha }(m,n)\\ h_{\beta }(m,n)\end{array}\right],$$

$$h_l=\frac{1}{\sqrt{3}}\left[\begin{array}{cc}1& 1\\ 1& 1\end{array}\right],h_{\alpha }=\frac{1}{\sqrt{6}}\left[\begin{array}{cc}2& -1\\ -1& 2\end{array}\right],h_{\beta }=\frac{1}{\sqrt{2}}\left[\begin{array}{cc}0& 1\\ -1& 0\end{array}\right].$$

$${\tilde{h}}_i(u,v)=\sum _{m=0}^{M-1}\sum _{n=0}^{N-1}{h}_i(m,n)e^{-j2\pi \left(\frac{mu}{M}+\frac{nv}{N}\right)}$$

$${\tilde{h}}_l=\frac{1}{\sqrt{3}}\left[\begin{array}{cc}4& 0\\ 0& 0\end{array}\right],{\tilde{h}}_{\alpha }=\frac{1}{\sqrt{6}}\left[\begin{array}{cc}2& 0\\ 0& 6\end{array}\right],{\tilde{h}}_{\beta }=\frac{1}{\sqrt{2}}\left[\begin{array}{cc}0& -2\\ 2& 0\end{array}\right].$$

$${\tilde{c}}_i(\mathbf{\omega })=\sum _{u=0}^{M-1}\sum _{v=0}^{N-1}{\tilde{h}}_i(u,v)\delta \left(\mathbf{\omega }-\mathbf{\omega }(u,v)\right),$$

$$\boldsymbol{\omega}(u,v)=\left(\frac{2\pi u}{M},\frac{2\pi v}{N}\right)\;\forall\begin{array}[]{c}u\in\{0,\ldots,M-1\}\\ v\in\{0,\ldots,N-1\}\end{array}.$$

$$\tilde{\theta }(\mathbf{\omega })$$

$$=\sum _{i\in \{l,\alpha ,\beta \}}\sum _{u=0}^{M-1}\sum _{v=0}^{N-1}{\tilde{h}}_i(u,v){\widetilde{\text{im}}}_i(\mathbf{\omega }-\mathbf{\omega }(u,v)),$$

$$\mathbf{x}=\left[\begin{array}{c}{\mathbf{h}}_r\\ {\mathbf{h}}_g\\ {\mathbf{h}}_b\end{array}\right]\in {\mathbb{R}}^{3K\times 1}.$$

$$\left[\begin{array}{c}{\mathbf{h}}_r\\ {\mathbf{h}}_g\\ {\mathbf{h}}_b\end{array}\right]=\left[\begin{array}{c}{\mathbf{Z}}_r\\ {\mathbf{Z}}_g\\ {\mathbf{Z}}_b\end{array}\right]\mathbf{x}\text{and}\left[\begin{array}{c}{\mathbf{h}}_l\\ {\mathbf{h}}_{\alpha }\\ {\mathbf{h}}_{\beta }\end{array}\right]=\left[\begin{array}{c}{\mathbf{Z}}_l\\ {\mathbf{Z}}_{\alpha }\\ {\mathbf{Z}}_{\beta }\end{array}\right]\mathbf{x},$$

$$\displaystyle\begin{array}[]{lcl}\boldsymbol{Z}_{r}=[\boldsymbol{I},\boldsymbol{0},\boldsymbol{0}]&&\boldsymbol{Z}_{l}=[\boldsymbol{I},\boldsymbol{I},\boldsymbol{I}]/\sqrt{3},\\ \boldsymbol{Z}_{g}=[\boldsymbol{0},\boldsymbol{I},\boldsymbol{0}]&\mbox{ and }&\boldsymbol{Z}_{\alpha}=[-\boldsymbol{I},2\boldsymbol{I},-\boldsymbol{I}]/\sqrt{6},\\ \boldsymbol{Z}_{b}=[\boldsymbol{0},\boldsymbol{0},\boldsymbol{I}]&&\boldsymbol{Z}_{\beta}=[\boldsymbol{I},\boldsymbol{0},-\boldsymbol{I}]/\sqrt{2}.\end{array}$$

$${\tilde{\mathbf{h}}}_i=\mathbf{F}{\mathbf{h}}_i=\mathbf{F}{\mathbf{Z}}_i\mathbf{x},i\in \{l,\alpha ,\beta \}.$$

$${\tilde{h}}_i(u,v)=\text{mbox}{vec}^{-1}({\tilde{\mathbf{h}}}_i)$$

$${\gamma }_l\stackrel{\text{def}}{=}\frac{1}{K}||{\tilde{\mathbf{h}}}_l|{|}_2,\text{ and }{\gamma }_c\stackrel{\text{def}}{=}\frac{1}{K}\text{min}\left(||{\tilde{\mathbf{h}}}_{\alpha }|{|}_2,||{\tilde{\mathbf{h}}}_{\beta }|{|}_2\right).$$

$$\begin{array}{rl}{\gamma }_l(\mathbf{x})& =\frac{1}{K}{1}^T{\mathbf{Z}}_l\mathbf{x}={\mathbf{b}}^T\mathbf{x}\\ {\gamma }_c(\mathbf{x})& =\min \left(\sqrt{{\mathbf{x}}^T{\mathbf{Q}}_{\alpha }\mathbf{x}},\sqrt{{\mathbf{x}}^T{\mathbf{Q}}_{\beta }\mathbf{x}}\right),\end{array}$$

$${\tilde{\mathbf{h}}}_l-\text{mbox}diag({\mathbf{e}}_1){\tilde{\mathbf{h}}}_l=0,$$

$$(\mathbf{I}-\text{mbox}diag({\mathbf{e}}_1))\mathbf{F}{\mathbf{Z}}_l\mathbf{x}=0.$$

$$\left[\begin{array}{c}{\mathbf{W}}_{\alpha }\\ {\mathbf{W}}_{\beta }\end{array}\right]\mathbf{x}=\mathbf{W}\mathbf{x}=0$$

$$J_l\stackrel{\text{def}}{=}\frac{1}{HW}{\int }_{[-\pi ,\pi {)}^2}\frac{\left(S_l(\mathbf{\omega })S_{\alpha }(\mathbf{\omega })+S_l(\mathbf{\omega })S_{\beta }(\mathbf{\omega })\right)}{S_{\theta }(\mathbf{\omega })}d\mathbf{\omega },$$

Full text

I Introduction

I-A Quanta Image Sensor

Quanta Image Sensor (QIS) is a new type of image sensors proposed by E. Fossum in 2005 [1, 2, 3] as a candidate for the third generation digital image sensor after CCD and CMOS. The sensor comprises a massive array of sub-diffraction limit single-photon detectors, called “jots”, with a pixel pitch of $1.1\mu$m as of today. Having read-out noise of $0.21$e- and dark current of 0.16e-/sec, QIS can count incoming photons to produce a digital output of bit depth in the range of $1-5$ bits, assuming that the exposure does not saturate the sensor. As reported in [2] and [3], the latest QIS sensor can achieve a resolution of $1024\times 1024$ jots and frame rate up to $1000$ fps by using the commercial 45/65nm 3D-stacked backside illumination CMOS process. Unlike CMOS image sensors that accumulate photons, QIS oversamples the scene by producing binary measurements. Figure 1 depicts the image formation process. See also [4, 5, 6] for detailed discussions.

Despite the rapid advancement in QIS hardware [1, 2, 3] and algorithms [7, 4, 5, 6], all reported findings, to-date, are based on monochromatic data. The first color QIS imaging was recently presented by Gnanasambandam et al. [8], where they demonstrated how to reconstruct a color image from the sensor using a Bayer color filter array. In this paper, we discuss how to design color filter array for better image acquisition.

I-B Color Filter Arrays for QIS

A color filter array (CFA) is a mask placed on top of the sensor to select (filter) wavelengths. As light passes through the color filter array, the resulting image is a mosaic pattern of the three tri-stimulus RGB colors. Traditionally, CFA is organized as a periodic replica of a 2D kernel called the color atom. The de-facto color atom used in the industry is the Bayer pattern [9] because of its simplicity and the readily available demosaicking algorithms, see for instance [10, 11, 12, 13, 14, 15, 16, 17, 18, 19]. More sophisticated CFAs have been proposed [20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30] to improve the Bayer CFA.

When designing a CFA, there are three factors that should be taken into consideration:

• •

Aliasing: Since color filtering is a sampling process, aliasing happens when the sampling rate is less than Nyquist. Aliasing causes false color artifacts at color edges, called the Moiré artifacts [21]. Color filters that are susceptible to aliasing, such as the Bayer CFA, require sophisticated demosaicking algorithms to suppress the Moiré artifacts. In contrast, a robust CFA can use less complicated demosaicking algorithms.

• •

Sensitivity: Since CFA is a filter, it blocks part of the incoming light. This reduces the sensor sensitivity and makes the image more susceptible to noise. A good CFA design should maximize the sensitivity by allowing transparent or “panchromatic” color filters that block as few wavelengths as possible.

• •

Crosstalk: Crosstalk can be either optical or electrical [31]. If not treated, crosstalk will make colors look pale or de-saturated. Crosstalk desaturation is corrected by pixel-wise multiplication of the RGB color vector using a color correction matrix. However, color correction enhances residual noise in the image [32, 31]. The situation is worsen in QIS because of its small size.

The three factors above are conflicting: Optimizing one generally degrades the others. For conventional CMOS image sensors, crosstalk is not severe, and so most CFA designs in the literature consider aliasing and sensitivity only. The first QIS color filter array design is proposed by Anzagira and Fossum [31]. However, aliasing is not considered because QIS is sufficiently small. We generalize these prior work by considering all three factors so that our framework can be applied to both QIS and CMOS image sensors.

I-C Related Work

The design framework we propose in this paper is a unification of several mainstream CMOS-based color filter arrays. To put our paper in the proper context in the literature, we list a few of the known works.

Spatial CFA Design: By suppressing the Moiré artifacts and crosstalk while keeping the demosaicing algorithm simple, Lukac and Plataniotis [20] proposed a CFA and compared it with other CFAs using a universal demosaicking method. However, their work did not provide a mathematical framework to analyze the CFA optimality.

Spatio-Spectral CFA Design: Hirakawa and Wolfe [21] proposed a method through the spatial and spectral domain analysis. Their CFA reduces aliasing in the frequency domain, and possesses high sensitivity and numerical stability. Condat [33] extended the framework by optimizing luminance and chrominance sensitivity. He defined a new form of orthogonality between chrominance channels in the frequency domain. Hao et al. [24] and Wang et al. [25] proposed a framework based on symbolic discrete Fourier transform (DFT). Their CFA maximizes the numerical stability of linear demosaicking process under aliasing and physical constraints. Li et al. proposed spatio-spectral CFA design methods that are optimized for sensitivty [28, 29].

Learning-based CFA Design: By minimizing the average error on a color dataset, Lu and Vetterli [22] used an iterative algorithm to solve a least squares CFA design problem. Chakrabarti [34] and Henz et al. [35] proposed to learn the optimal CFA pattern by using a deep neural network.

Besides these mainstream CFA design frameworks, there are a number of other CFA designs such as [23, 27, 26, 30]. On the hardware side, Biay-Cheng et al. [36, 37] took into account that color filter fabrication technology lags the image sensor technology in terms of miniaturization. They proposed a hardware-friendly CFA assuming the color filter size is $1.5\times$ pixel size.

I-D Scope and Contributions

In this paper, we propose an optimization framework that encompasses aliasing, sensitivity and crosstalk in a unified model. This is the first time we can incorporates a quantitative crosstalk metric in an optimization framework for CFA design. Our framework is general as it works with any crosstalk kernel depending on the sensor. Moreover, it works with any color atom size by a mild change in the formulation.

The main contribution of this paper is a general and flexible framework for CFA design. Compared to the existing CFA design framework, the new framework is able to simultaneously (Section III)

• •

Improve CFA’s luminance and chrominance sensitivity,

• •

Reduce aliasing between luminance and chrominance channels,

• •

Suppress crosstalk between spectral sub-bands, and

• •

Enforce orthogonality between chrominance channels to permit simple linear demosaicking.

The design framework is presented in the form of optimization. We have two designs: A convex optimization and a non-convex optimization. In addition to the formulation, we also present an algorithm to solve the non-convex optimization. (Section IV)

For performance evaluation, we generalize the demosaicking algorithm of Condat [33] in Section V so that it now works with arbitrary CFAs. This pipeline comprises a demosaicking by frequency selection algorithm for removing the CFA masking effect followed by a color correction step for removing the desaturation effect of crosstalk. Experimental evaluation on the DIV2K evaluation dataset in Section VI shows the robustness of our proposed CFAs compared to other CFAs in literature.

II Background and Notations

To facilitate readers to understand the design framework, in this section we introduce a few notations and terminologies. We will start in Section II-A by describing the image formation using a CFA, then we discuss CFA in different domains in Sections II-B and II-C. Afterwards, in Section II-D, we define the optimization variables to simplify the design framework.

II-A Color Image Formation

Consider a color image $\mathbf{im}_{rgb}$ of size $H\times W$. We denote the normalized light intensities in the red, green and blue channels for the ($m$-th,$n$-th) pixel of the color image as

[TABLE]

where $m\in\{0,\ldots,H-1\}$, $n\in\{0,\ldots,W-1\}$.

Color Filtering: To obtain color, we place a color filter on top of each jot to collect light for one of the RGB colors. The CFA is a periodic pattern of the same resolution of $\mathbf{im}_{rgb}$, defined as

[TABLE]

where $c_{r}(m,n)$, $c_{g}(m,n)$, $c_{b}(m,n)\in[0,1]$ are the opacity rates for the red, green and blue pixels, respectively. For example, a red color filter is defined as $\boldsymbol{c}_{rgb}(m,n)=[1,0,0]^{T}$ as it only passes the red color. The light exposure on the QIS after passing through the CFA is denoted as $\theta(m,n)$, which is a linear combination of the tri-stimulus colors:

[TABLE]

Here, $\eta$ is a positive scalar representing the sensor gain factor.

Photon Arrival. The photon arrival is modeled as a Poisson process. Let $\boldsymbol{Y}\in\mathbb{N}^{HW}$ be a vector of non-negative random integers denoting the number of photons arriving at QIS jots according to the light exposure $\boldsymbol{\theta}$. Then, the probability of observing a photon count $Y_{j}=y_{j},j\in\{1,\ldots,HW\}$ is

[TABLE]

In this work, we assume single-bit QIS [38] that quantizes the photon counts by QIS jots to a binary values $\boldsymbol{B}\in\{0,1\}^{HW}$ with $B_{j}=1$ if $Y_{j}\geq q$ and $B_{j}=0$ if $Y_{j}<q$, where $q>0$ is a threshold. The probability of observing $B_{j}=b_{j}$ is

[TABLE]

where $\Psi_{q}(\cdot)$ is the incomplete Gamma function [6].

Temporal Oversampling. With frame rates that reach $1000$ fps, QIS is able to catch the scene movement by taking $T$ temporal samplings for the same scene. This allows us to utilize multiple independent measurements over time to improve the statistics and decrease noise. Hence, for every jot with light exposure $\theta_{j}$, we have a set of $T$ independent binary measurements $\mathcal{B}_{j}=\{b_{j,0},\ldots,b_{j,T-1}\}$.

II-B Color Filter Array Analysis in Different Color Spaces

Since the CFA $\boldsymbol{c}_{rgb}(m,n)$ is a periodic pattern, it is sufficient to use a color atom as the optimization variable when designing the CFA. The color atom takes the form

[TABLE]

where each of $h_{r}$, $h_{g}$ and $h_{b}$ is an $M\times N$ array. For instance, the GRBG Bayer pattern has the following color atom (when $M=N=2$):

[TABLE]

because the Bayer pattern has one red pixel and one blue pixel located at two opposite diagonals, and two green pixels located in the remaining two positions. Figure 2 illustrates the idea.

While the primal RGB color is common for making the CFA, it would be more convenient if the colors are decorrelated. To this end, we change the image representation from the canonical RGB basis to an orthornormal basis using a transformation matrix [39, 33]:

[TABLE]

This transformation maps an RGB image $\mathbf{im}_{rgb}$ to an image $\mathbf{im}_{l\alpha\beta}=[\textrm{im}_{l},\textrm{im}_{\alpha},\textrm{im}_{\beta}]^{T}$ as follows (we drop the spatial indices $(m,n)$ for simplicity)

[TABLE]

where $\textrm{im}_{l}$ is a luminance (luma) component that contains the high frequency components such as edges and textures, whereas $\textrm{im}_{\alpha}$ and $\textrm{im}_{\beta}$ are chrominance (chroma) components that carry the color information.

Since $\boldsymbol{T}$ is orthonormal (i.e., $\boldsymbol{T}^{T}\boldsymbol{T}=\boldsymbol{I}$), we can rewrite the sampling process in Equation (3) in the luma/chroma space:

[TABLE]

where $c_{l}(m,n)$, $c_{\alpha}(m,n)$ and $c_{\beta}(m,n)$ are the luma/chroma representation of the CFA, with

[TABLE]

The luma/chroma representation of the CFA has a corresponding color atom $h_{l}(m,n)$, $h_{\alpha}(m,n)$ and $h_{\beta}(m,n)$. For instance, the luma/chroma color atom of the GRBG Bayer pattern is

[TABLE]

where the individual components are

[TABLE]

Remark 1**.**

In principle, there are infinite choices for the luma/chroma basis $\boldsymbol{T}$. We choose the one in Equation (7) because it makes the components of natural images statistically independent in the first order approximation.

II-C Color Filter Array in Fourier Space

When analyzing the aliasing effects of the CFAs, we need to transform the color atom into the Fourier domain. For simplicity, we represent the Fourier transform of a signal by putting a tilde on top of the symbol, e.g., $h\overset{\mathcal{F}}{\rightarrow}\widetilde{h}$. The 2D discrete Fourier transform (DFT) of the $i$-th color atom is

[TABLE]

where $u\in\{0,\ldots,M-1\}$, $v\in\{0,\ldots,N-1\}$.

For example, the discrete Fourier transform of the luma/chroma color atoms in Equation (10) are

[TABLE]

Here, we observe that the Fourier transform of the color atom has the same size as the original color atom. The luminance channel has only one baseband components at $(0,0)$, whereas the $\alpha$ chrominance channel has one baseband component and a component at $(\pi,\pi)$. The $\beta$ chrominance channel has two components at $(0,\pi)$ and $(\pi,0)$. Figure 3 illustrates how these frequency locations are identified from a $3\times 3$ color atom.

While the Fourier transform of the color atom is useful, for demosaicing we also need to analyze the spectrum of the entire CFA. As shown by Hao et al. [24], the Fourier transform of the entire CFA can be written in terms of the Fourier transform of the color atoms:

[TABLE]

where $i\in\{l,\alpha,\beta\}$, $\boldsymbol{\omega}$ is the 2D angular frequency, and

[TABLE]

is the ($u$-th,$v$-th) 2D angular frequency. It is worth noting that the Fourier transform of the CFA comprises pure sinusoids of amplitudes $\widetilde{h}_{i}(u,v)$. These sinusoids are placed at $MN$ discrete 2D frequencies $\boldsymbol{\omega}(u,v)$ that divide the 2D frequency plane $[-\pi,\pi]\times[-\pi,\pi]$ into $MN$ equal regions. Therefore, the spectrum of the mosaicked image $\widetilde{\theta}(\boldsymbol{\omega})$ can be written as

[TABLE]

where $\circledast$ is the standard 2D convolution operator. Having the spectrum of the mosaicked image $\widetilde{\theta}(\boldsymbol{\omega})$, we can now discuss the optimization variables in our problem.

II-D Design Variables

To formulate the CFA design problem as an optimization problem, we define the following variables. We denote $\boldsymbol{h}_{r}$, $\boldsymbol{h}_{g}$ and $\boldsymbol{h}_{b}$ the vectorized representations of the red, green and blue color atoms, respectively. To ensure physical realizability, we require $\boldsymbol{h}_{r}$, $\boldsymbol{h}_{g}$, $\boldsymbol{h}_{b}\in\left[0,1\right]^{K\times 1}$, where $K\overset{\text{def}}{=}MN$, and we stack all design variables into one long vector

[TABLE]

The design variable $\boldsymbol{x}$ is related to the vectorized RGB and luma/chroma color atoms as

[TABLE]

where the $\boldsymbol{Z}$ matrices are defined by Equation (7) as

[TABLE]

Given the design variable $\boldsymbol{x}$, we also need to analyze its spectrum. We write the 2D Fourier transform equation (11) as a matrix-vector product by using the Fourier transform matrix $\boldsymbol{F}\in\mathcal{C}^{K\times K}$. Hence, the vectorized spectra of the luma/chroma color atoms can be written in terms of $\boldsymbol{x}$ as

[TABLE]

where $\widetilde{\boldsymbol{h}}_{i}\in\mathcal{C}^{K\times 1}$, for $i\in\{l,\alpha,\beta\}$. The relation between the matrix and the vector forms of the Fourier transform is:

[TABLE]

where $\widetilde{h}_{i}(u,v)$ is the Fourier coefficient.

III Design Criteria

We now present the design criteria. Our criteria unify the three major approaches in the literature: (i) Sensitivity of luma/chroma channels to noise by Condat [33]; (ii) Aliasing between different color components in the frequency domain by Hirakawa and Wolfe [21]; (iii) Crosstalk between neighboring pixels in the spatial domain by Anzagira and Fossum [31]. Note that the first two criteria were developed for CMOS, whereas the third criterion was developed for QIS. The proposed framework integrates all these criteria into a unified formulation. Table I summarizes the difference between this paper and the previous works. Figure 6 shows an example CFA, and its frequency representation.

In the following subsections, we present the design criteria and express them in terms of matrix-vectors for the optimization framework in Section IV.

III-A Luminance and Chrominance Sensitivity

Definition 1**.**

The luminance sensitivity $\gamma_{l}$ and the chrominance sensitivity $\gamma_{c}$ of a CFA with color atom $\{\boldsymbol{h}_{l},\boldsymbol{h}_{\alpha},\boldsymbol{h}_{\beta}\}$ of size $M\times N$ are defined as

[TABLE]

where $K=MN$ is a normalization factor.

Intuitively, the luminance and chrominance sensitivity are measures of the signal power that can be transmitted through the color filter. A more transparent color filter allows more light to pass through, and hence the signal power is higher. This is reflected by the magnitudes $\|\widetilde{\boldsymbol{h}}_{i}\|_{2}$ for $i\in\{l,\alpha,\beta\}$, which according to Parseval’s Theorem they are equivalent to $\|\boldsymbol{h}_{i}\|_{2}$.

The following proposition shows how can we compute $\gamma_{l}$ and $\gamma_{c}$ in terms of the optimization vector $\boldsymbol{x}$.

Proposition 1**.**

For a CFA with color atoms represented by the vector $\boldsymbol{x}$, the luminance and chrominance sensitivity can be calculated as

[TABLE]

where $\boldsymbol{b}=\frac{1}{K}\boldsymbol{Z}_{l}^{T}\boldsymbol{1}$, $\boldsymbol{Q}_{\alpha}=\boldsymbol{Z}_{\alpha}^{T}\boldsymbol{Z}_{\alpha}$ and $\boldsymbol{Q}_{\beta}=\boldsymbol{Z}_{\beta}^{T}\boldsymbol{Z}_{\beta}$.

Proof:

See Appendix A. ∎

The luminance sensitivity and the chrominance sensitivity cannot be arbitrarily chosen. One practical consideration is to ensure uniform noise power across the luma channel so that the denoising procedure can be simplified (because the noise will be i.i.d.). Thus, the luminance color atom $h_{l}$ should be constant, i.e., $h_{l}(m,n)=c,\forall m,n$, where $c$ is a positive constant. Taking Fourier transform, this means that $\widetilde{h}_{l}$ comprises only one impulse at baseband $\widetilde{h}_{l}(0,0)$, and no impulses at all other frequencies. In vector form, we need

[TABLE]

where $\boldsymbol{e}_{1}=[1,0,\ldots,0]^{T}$ is the standard basis. Putting in terms of the optimization variable $\boldsymbol{x}$, we have a constraint.

Proposition 2** (Uniform Luminance Constraint).**

If a CFA has a uniform luminance sensitivity, then $\boldsymbol{x}$ needs to satisfy

[TABLE]

Proof.

Using Equation (15), substitute $\widetilde{\boldsymbol{h}}_{l}=\boldsymbol{F}\boldsymbol{Z}_{l}\boldsymbol{x}$ into Equation (19). ∎

III-B Anti-Aliasing

In the frequency representation of a mosaicked image, the luminance controls the baseband whereas the chrominance controls the sideband of the spectrum. To minimize spectral interference, i.e., aliasing, CFA design methods modulate the chrominance as far as possible from the baseband. Typically, these methods do not put chrominance on the vertical and horizontal axes to prevent aliasing. However, since QIS is very small and it oversamples the scene, it is possible to relax the aliasing constraint. To this end, we allow chrominance channels to be placed on the vertical and horizontal axis as long as they are placed at frequencies higher than or equal $\pi/2$. This enlarges the feasible set of the solutions.

Mathematically, the anti-aliasing requirement is formulated by forcing the Fourier coefficients of the chrominance color atoms at $(\pm\pi,v)$ and $(u,\pm\pi)$ to zero for all $u$ and $v$ greater than $\pi/2$. In terms of our design variable $\boldsymbol{x}$, we require the following constraint.

Proposition 3** (Anti-aliasing Constraint).**

The chrominance placed on frequencies less than $\pi/2$ in the vertical and horizontal directions must be set to 0. Hence, $\boldsymbol{x}$ must satisfy

[TABLE]

where $\boldsymbol{W}_{\alpha}$ and $\boldsymbol{W}_{\beta}$ are the matrices formed by choosing the rows in $\boldsymbol{F}\boldsymbol{Z}_{\alpha}$ and $\boldsymbol{F}\boldsymbol{Z}_{\beta}$, respectively, that correspond to vertical and horizontal frequency components that are less than $\pi/2$.

To quantify the amount of aliasing for every CFA, we define the following aliasing criterion.

Definition 2**.**

For a CFA, aliasing between luminance and chrominance channels is measured by

[TABLE]

where $S_{l}$, $S_{\alpha}$, $S_{\beta}$ and $S_{\theta}$ denote the power spectral density of the luminance channel $\mathbf{im}_{l}$, the two chrominance channels $\mathbf{im}_{\alpha}$ and $\mathbf{im}_{\beta}$, and the mosaicked image $\boldsymbol{\theta}$, respectively.

III-C Crosstalk

Crosstalk is caused by the leakage of electrical and optical charge from a pixel to its adjacent pixels [40, 31]. Crosstalk leads to color de-saturation. To model crosstalk, we follow [31] by defining three scalars $\delta_{r}$, $\delta_{g}$, and $\delta_{b}$ representing the proportion of leaked charges to neighboring pixels. These three scalars then form a crosstalk kernel,

[TABLE]

which can be considered as a spatial lowpass filter of the mosaicked image. Applying the crosstalk kernel to the CFA is equivalent to a spatially invariant convolution

[TABLE]

where $h_{i}^{\mathrm{ctk}}$ denotes the effective CFA in the presence of crosstalk.

The effect of crosstalk is more severe when the adjacent colors are different. For example, in Figure 5, the red and blue pixels are surround by 8 neighbors of different colors and the green pixels are surrounded by 4 neighbors of different colors. This is equivalent to saying that there is a red pixel having a value $1$ and is surrounded by pixels having the value [math]. Using similar argument, we can see that if the color atoms have more rapid changes of colors, then the resulting CFA is more susceptible to crosstalk.

We propose to quantify the variation of the color atoms (and hence crosstalk) is by means of measuring the total variation of the color atom. The total variation is a proxy of the complexity of the color filter array. A color filter array with high total variation means a more complicated pattern and so it is more susceptible to crosstalk. Our total variation is defined as follows.

Definition 3** (Total Variation).**

For a CFA defined by the color atoms $\boldsymbol{h}_{r}$, $\boldsymbol{h}_{g}$ and $\boldsymbol{h}_{b}$, the weighted total variation is defined as

[TABLE]

where $\boldsymbol{D}\overset{\text{def}}{=}[\boldsymbol{D}_{x},\boldsymbol{D}_{y}]^{T}$ is an operator that computes the vertical and horizontal derivatives with circular boundary conditions, and $\delta_{i}$ is the leakage factor defined in the crosstalk kernel in Equation (23).

To control the amount of variations in the CFA (so that we can limit the amount of crosstalk), we upper bound the total variation by a scalar $\textrm{TV}_{\mathrm{max}}$. This leads to the following constraint.

Proposition 4** (Crosstalk Constraint).**

The crosstalk is limited by upper-bounding the total variation metric $\mathrm{TV}(\boldsymbol{x})$:

[TABLE]

Figure 6 shows two CFAs proposed in literature. The first one, proposed in [24] is more robust to aliasing than the second one proposed in [31]. This is because the chrominance channels are modulated at high frequencies which are far from baseband luminance. However, [31] is more robust to crosstalk than [24] because the color atom have less variation in colors. We can also see this in the total variation values: $0.206$ compared to $0.131$ at $(\delta_{r},\delta_{g},\delta_{b})=(0.23,0.15,0.10)$. This trade-off constitutes a gap in literature: Color filter designs can improve robustness of either aliasing or crosstalk, but not for both. Our proposed design framework allows us to optimize them simultaneously.

Since the best value for $TV_{max}$ changes for different atom sizes, we choose it according to the best value in the state-of-the-art CFA atoms. For example, we choose $TV_{max}=0.131$ for $4\times 4$ atoms since this is the lowest value achieved by Anzagira and Fossum $4\times 4$ RGBCWY CFA [31]. Other $4\times 4$ CFAs have higher $TV$ values.

III-D Orthogonality of Chrominance Channels

When designing a CFA, one should take into consideration of the complexity of the demosaicking process. Recall Equation (II-C) where we show that

[TABLE]

This is a modulation of the signal by a modulating frequency $\boldsymbol{\omega}(u,v)$. Therefore, to reconstruct the signal, one approach is to demodulate by shifting the channels to the baseband by multiplying pure sinusoids and then applying a lowpass filter [33]. Demodulation can be done efficiently if there is orthogonality between the channels. Following the literature, our optimization takes into account of two forms of orthogonality.

• •

Total Orthogonality [24] and [21]: The idea is to make one chroma channel zero and the other non-zero at any $(u,v)$. For example, $\widetilde{h}_{\alpha}(u,v)=0.9$ and $\widetilde{h}_{\beta}(u,v)=0$.

• •

Quadrature Orthogonality [33]: The idea is to make one chroma channel real and the other imaginary at any $(u,v)$, i.e., the two channels are modulated by a frequency $\boldsymbol{\omega}(u,v)$ but in quadrature phase. Translating the spatial domain, this means that

[TABLE]

where $m\in\{0,\ldots,M-1\}$, $n\in\{0,\ldots,N-1\}$, and $\phi$ is the phase angle. In this way, the two channels can be easily separated during the demosaicking process using the orthogonality of cosine and sine functions.

We formulate the orthogonality criteria as a penalty function that is a surrogate of both approaches.

Definition 4** (Orthogonality Penalty).**

For a CFA having chrominance channels with spectra $\widetilde{\boldsymbol{h}}_{\alpha}(\boldsymbol{\omega})$ and $\widetilde{\boldsymbol{h}}_{\beta}(\boldsymbol{\omega})$, the orthogonality penalty is defined as

[TABLE]

which can be written as a function in $\boldsymbol{x}$ as follows

[TABLE]

Looking at the first summation in Equation (4), we notice that for every 2D frequency $(u,v)$, the term $|\Re{\widetilde{h}_{\alpha}(u,v)}|+|\Re{\widetilde{h}_{\beta}(u,v)}|$ is the $\ell_{1}$-norm of a 2-dimensional vector $[\Re{\widetilde{h}_{\alpha}(u,v)},\Re{\widetilde{h}_{\beta}(u,v)}]^{T}$. Therefore, minimizing this $\ell_{1}$-norm promotes either one of the components to zero (or both). Similar argument applies for the imaginary components in the second summation. As a result, the total variation can be regarded as a proxy to the orthogonality condition.

III-E Condition Number

When designing a color filter array, one should also be aware of the simplicity of the demosaicking algorithm. Since the luminance/chrominance transformation, color filtering and crosstalk are all linear processes, we can represent them by an overall color acquisition matrix $\boldsymbol{A}$. To demosaic the image, in principle we need to invert the $\boldsymbol{A}$ matrix. To avoid the amplification of the estimation error of luminance and chrominance channels, the condition number of $\boldsymbol{A}$ should be minimized for numerical stability. This metric was discussed in [24], but the authors considered the condition number of the luminance/chrominance transformation matrix $\boldsymbol{T}$ only. In our case, we generalize this metric by taking the color filtering and crosstalk into account as well.

To represent the image acquisition in frequency domain as a linear process, we assume the crosstalk kernels for red, green and blue pixels are the same $g_{r}=g_{g}=g_{b}$. Define the following frequency domain variables:

[TABLE]

where $\widetilde{\boldsymbol{g}}\overset{\mathcal{F}}{\leftarrow}\boldsymbol{g}$ is the vectorized version of the $M\times N$ discrete Fourier transform of the crosstalk kernel $g$. Hence, the mosaicked image $\widetilde{\theta}$ can be written as

[TABLE]

where we define the color acquisition matrix as $\boldsymbol{A}\overset{\text{def}}{=}\widetilde{\boldsymbol{G}}\widetilde{\boldsymbol{H}}\boldsymbol{T}$. Denote by $\kappa(\boldsymbol{A})$ the condition number of $\boldsymbol{A}$, i.e.,

[TABLE]

Low values of $\kappa(\boldsymbol{A})$ imply stable demosaicking process that involves mild amplification of estimation errors in the luminance and chrominance components.

IV Formulation of Optimal CFA Design Problem

Using the variables and constraints defined in the previous section, we present two different optimization formulations of the CFA design problem in this section: (i) A non-convex formulation that integrates all the above information into a single optimization, and (ii) convex relaxation that makes the problem more tractable.

IV-A Non-Convex CFA Design

The non-convex CFA optimization puts the objectives and constraints defined in the previous section into a single optimization problem. This gives us

[TABLE]

where $\lambda_{l}$ and $\lambda_{\rho}$ are the regularization parameters controlling the relative weights of the luminance sensitivity and the orthogonality penalty. The penalty function $\rho(\boldsymbol{x})$ is added to the objective with a negative sign so that it is minimized. By lower bounding $\gamma_{c}(\boldsymbol{x})$ with a constant $\tau$, we can rewrite Equation (32) as

[TABLE]

In this optimization problem, the objective and constraints are convex except for (39)(f) and (39)(g). This is because these inequalities include convex quadratic form in the “$\geq$” side, where convexity comes from the fact that $\boldsymbol{Q}_{\alpha}$ and $\boldsymbol{Q}_{\beta}$ are positive semidefinite matrices. Hence, the optimization problem is non-convex.

IV-B Solving the Optimization

While problem (39) is non-convex, we can find a local minimum by successive convex approximations [41]. The idea of successive convex approximation is to replace the quadratic terms in the non-convex constrains (39)(f) and (39)(g) by first order approximations around the initial guess $\boldsymbol{x}^{(0)}$. Since the quadratic form is convex, its first order approximation constitutes a lower bound. Hence, we are replacing the non-convex constraints (39)(f) and (39)(g) with convex but tighter constraints that limit the feasible set of $\boldsymbol{x}$. The algorithm repeats until $\tau$ converges to a fixed-point, which is the final chrominance sensitivity.

The overall algorithm is summarized in Algorithm 1. Figure 7 shows the convergence of Algorithm 1 for designing a $4\times 4$ color atom. We notice the monotonic increase of $\tau$ until it converges to a fixed point. Since the original problem is non-convex, solution to the problem could be a local minimum depending on how the initialization is done. In practice, we solve the problem for multiple instances with different randomly generated initial guesses which approximately cover the design space (e.g., using the Latin hypercube sampling [42]), and pick the best solution among them.

IV-C Convex CFA Design

The relaxation from non-convex to convex can be done by explicitly forcing part of the chrominance components to zero. Specifically, we modulate the chrominance channels on the same frequency $\boldsymbol{\omega}(u,v)=(\frac{2\pi u}{M},\frac{2\pi v}{N})$ using the quadrature orthogonality mentioned in Equation (26). In terms of $\boldsymbol{x}$, these two equations can be written as:

[TABLE]

where $\boldsymbol{x}_{c}$ and $\boldsymbol{x}_{s}$ are constant vectors that represent the vectorized version of the cosine and sine signals on the right hand side of Equations (26a) and (26b), respectively, i.e.,

[TABLE]

Since we explicitly choose the modulation frequencies of chrominance channels manually, we can drop the aliasing constraint in Proposition 3. However, we still need the uniform luminance constraint in Propositions 2. Moreover, since the luminance and chrominance gains are adversarial, the objective of this formulation is to maximize their weighted sum. To this end, the problem is written as:

[TABLE]

In our terminology, the optimization problem of [33] is obtained from Equation (50) by removing the crosstalk constraint (50)(d). Hence, our optimization limits the search space of the optimization in [33] to get CFAs that have acceptable crosstalk performance.

Figure 8 shows two color atoms obtained using the convex and non-convex formalizations. In the convex formulation, we select the modulation frequency to be $\omega_{0}=[\pi,\pi]$ and the phase that maximizes $\gamma_{c}$ at this frequency is found to be $\phi=\pi/12$. Then, we solve the problem to get $(\gamma_{l},\gamma_{c},TV)=(0.573,0.08,0.263)$. As for the non-convex formulation, we let the optimization to choose modulation frequencies subject to crosstalk and aliasing constraints. Solving the non-convex formulation yields $(\gamma_{l},\gamma_{c},TV)=(0.573,0.09,0.263)$. We notice that the non-convex formulation achieves higher chrominance sensitivity because of its flexibility in choosing the modulation frequencies.

V Universal Demosaicking

In this section, we present a universal demosaicking algorithm which can be used to all CFAs presented in this paper. Our algorithm comprises two main parts: (i) a demosaicking step to remove the color filtering effect (Section V-B) and (ii) a color correction step to mitigate the crosstalk effect (Section V-C).

V-A Special Consideration for QIS.

Before we talk about the demosaicking algorithm, we should briefly discuss the photon statistics of QIS. In CMOS, the measured voltage can be modeled as a nominal value corrupted by i.i.d. Gaussian noise. For single-bit QIS, previous work showed that the measured photon counts follow a truncated Poisson process [5]. When averaging over a number of temporal frames, the truncated Poisson becomes a Binomial [6]. If the photon count is sufficiently high, this binomial will approximately approach a Gaussian. Applying the law of large numbers on the distribution of $\boldsymbol{B}$ in Equation (5), the average is

[TABLE]

and the maximum-likelihood estimate of the signal is

[TABLE]

As discussed in [6], we can regard this equation as a tone-mapping of the photon counts. We regard $\theta_{j}$ as the $j$-th pixel of the mosaicked image generated by the CFA. The goal of demosaicing is to reconstruct a color image from $\theta_{j}$.

V-B Demosaicking by Frequency Selection

Our demosaicking algorithm is based on frequency selection [11]. It generalizes [12] as it works for any CFA as long as it satisfies the orthogonality constraints in Section III-D

The key idea of the algorithm is to shift every chrominance channel to the baseband by multiplying with its carrier, then use a low-pass filter to reconstruct it. For chrominance components that are replicated over distinct carriers, we combine them by simple averaging. After obtaining the $\alpha$ and $\beta$ chrominance channels, they are re-modulated to their original positions and subtracted from the mosaicked image to obtain the luminance channel. This process is summarized in Algorithm 2 for a special case of a CFA that has strictly one replica of the $\alpha$ and $\beta$ chrominance channels. It is also illustrated by Figure 9. Extension of the algorithm to the general case is straightforward.

To apply Algorithm 2 on CFAs proposed in [24], [31] and [36], they must satisfy the orthogonality constraints in Section III-B. However, this is not satisfied with our choice of the luminance/chrominance basis defined by $\boldsymbol{T}$ in Equation (7). Hence, we use for every CFA the transformation matrix $\boldsymbol{T}$ that makes its luminance and chrominance channel orthogonal. To ensure fairness, we normalize the matrix rows to unity so that all luminance and chrominance have the same noise power. The transformation matrices are provided in the supplementary.

V-C Color Correction

The demosaicking algorithm in Algorithm 2 does not take into account of the crosstalk effect. Like most of the mainstream image and signal processing (ISP) pipelines, we reduce the cross-talk via a color correction step.

Given the demosaicked color pixel $\widehat{\textbf{im}}(m,n)$, the color correction multiplies $\widehat{\textbf{im}}(m,n)$ by a $3\times 3$ matrix $\boldsymbol{M}$ such that $\boldsymbol{M}\widehat{\textbf{im}}(m,n)$ is the color-corrected pixel. The matrix $\boldsymbol{M}$ is learned by comparing a measured color pixel to a known color chart value. Mathematically, suppose we have $K$ measured color pixels forming a $3\times K$ matrix $\boldsymbol{Q}_{\textrm{False}}$, and a corresponding true color values forming another $3\times K$ matrix $\boldsymbol{Q}_{\textrm{GT}}$, we can estimate $\boldsymbol{M}$ by solving

[TABLE]

where $\epsilon_{c}(\boldsymbol{M})=\mathrm{Tr}\left\{\left(\boldsymbol{M}\boldsymbol{Q}_{\textrm{False}}-\boldsymbol{Q}_{\textrm{GT}}\right)^{T}\left(\boldsymbol{M}\boldsymbol{Q}_{\textrm{False}}-\boldsymbol{Q}_{\textrm{GT}}\right)\right\}$ is the squared color error. To minimize the noise amplification, it is advised to add regularization when estimating $\boldsymbol{M}$ [43]. Since a standard color chart comprises $24$ color patches, we can estimate the noise by computing the norm of covariance matrix of every color patch, and get the average value over the $24$ color patches. Hence, the optimization problem is rewritten as

[TABLE]

where $\boldsymbol{Q}^{(i)}_{\textrm{False}}$ represents the pixels of the $i$th color patch, and $\mu$ is a positive scalar that controls the noise amplification effect. By varying $\mu$, we can draw a tradeoff curve between color reproduction accuracy and noise amplification.

Figure 10 shows reconstructed images before and after color correction. In this figure, the crosstalk parameters are $(\delta_{r},\delta_{r},\delta_{r})=(0.23,0.15,0.1)$. We can see the effect of color correction in the more saturated red and yellow feathers and in the green leaves in the background.

VI Experimental Evaluation

In this section, we present CFAs obtained using our optimization framework in Section VI-A. Afterwards, we evaluate the performance of different CFAs using Algorithm 2. First, we simulate in Section VI-B the reconstruction performance of a Macbeth ColorChecker under D65 illuminant. Second, we show in Section VI-C a quantitative and qualitative comparison of the reconstruction performance of all CFAs on DIV2K evaluation dataset [44]. This dataset comprises 100 color images of resolution $2040\times 1356$. We downsample them to $1020\times 678$ for our experiments. In all experiments, QIS parameters are chosen as $(q,\eta,T)=(1,2,1000)$, where the value of $T$ is set as a high value to minimize the noise power.

VI-A Experiment 1: Proposed Solutions of CFA Design Problem

We focus on the non-convex formulation (39) since it is more flexible than the convex formulation (50). We set the parameters of formulation (39) as $\lambda_{l}=0.1$ and $\lambda_{\rho}=0.02$. We run multiple instances of Algorithm 1 ($2000$ instances) using different random initializations for the color atoms $\boldsymbol{x}^{(0)}$. Then, we pick the solution with the highest chrominance sensitivity. To ensure that the initial guess spans the feasible set of $\boldsymbol{x}$, we use uniform Latin hypercube sampling of the domain $[0,1]^{3L}$.

Figure 11 shows our proposed CFAs and their accompanied spectra compared to other CFAs in the literature. For every CFA, we compute 1) the luminance and chrominance gains ($\gamma_{l}$ and $\gamma_{c}$) in Equation (17) to measure robustness to noise, 2) the total variation metric $TV(x)$ (Proposition 4) to measure robustness to crosstalk, and 3) the aliasing metric $J_{l}$ in Equation (22) to measure robustness to aliasing, and the condition number $\kappa(\boldsymbol{A})$ defined in Section III-E. To calculate the aliasing metric for [21], [24], [31] and [36], we use the transformation matrices that make the luminance and chrominance channels orthogonal as mentioned in Section V-B. Leakage factors $(\delta_{r},\delta_{g},\delta_{b})$ are chosen as $(0.23,0.15,0.10)$. Results are summarized in Table II.

• •

$4\times 4$: Among all $4\times 4$ CFAs in Table II, Hao et al. [24] is the most robust CFA to aliasing, but the least robust to crosstalk; whereas RGBCWY [31] is the most robust to crosstalk, but it is not as robust to aliasing. Our CFA achieves the best of both worlds by having the same total-variation like RGBCWY, and good aliasing metric. Moreover, it has the lowest condition number and it is more sensitive since it comprises 4 panchromatic or “white” pixels with $(r,g,b)=(1/3,1/3,1/3)$.

• •

$3\times 3$: Compared to Biay-Cheng et al. [36], our CFA has less aliasing. The high aliasing metric of [36] is attributed to its design which overlooks frequency domain aliasing.

• •

$3\times 2$: Compared to Condat [33], our CFA is more robust to crosstalk and aliasing. It is worth noting that we get exactly the same CFA as Condat’s automatically when we remove the crosstalk constraint.

• •

$4\times 2$: Compared to Hirakawa-Wolfe [21], our CFA has higher chrominance sensitivity and it is more robust to crosstalk.

• •

$7\times 7$: Compared to Bai et al. [28], our CFA is more robust to crosstalk and aliasing.

VI-B Experiment 2: Macbeth ColorChecker Reconstruction

In this experiment, we simulate the performance of different CFAs in reconstructing the Macbeth ColorChecker image. Pixel response is determined using the incident photon flux of D65 light and the spectral reflectance of Macbeth ColorChecker integrated over the visible light spectrum. QIS parameters and primary color filters are taken from [31]. For every CFA, we generate mosaicked images under two scenarios: 1) crosstalk kernels with leakage factors $(\delta_{r},\delta_{g},\delta_{b})=(0,0,0)$, i.e., no crosstalk, and 2) crosstalk kernels with leakage factors $(\delta_{r},\delta_{g},\delta_{b})=(0.45,0.30,0.20)$ as suggested by [31]. For fairness of comparison, we use Algorithm 2 for demosaicking all CFAs including RGBCY and RGBCWY CFAs proposed in [31]. Color correction with white balance is performed after color demosaicking for the crosstalk case for removing the crosstalk effect.

We use the following metrics [31] to evaluate the CFAs:

• •

Sensitivity metamerism index (SMI) which measures the drop in color reproduction accuracy due to crosstalk. It is obtained as a function of the CIEDE2000 color error metric which is obtained by calculating the mean square color difference in the CIELAB color space.

• •

Luminance SNR (YSNR) which measures the visual noise of luminance channel as defined in ISO 12232 [45].

Table II shows these metrics for different CFAs with and without crosstalk. Our CFAs achieve higher color reproduction accuracy compared to others since they are optimized for crosstalk. This gain in color accuracy happens by trading the noise performance as observed by the drop of YPSNR metric.

VI-C Experiment 3: Natural Color Image Reconstruction

In this experiment, we evaluate the performance of different CFAs for natural color image reconstruction. To this end, we use the DIV2K evaluation dataset to generate 100 mosaicked images according to QIS model. Two scenarios are assumed: 1) No crosstalk, and 2) Moderate crosstalk with leakage factors $(\delta_{r},\delta_{g},\delta_{b})=(0.23,0.15,0.10)$. The low pass filter in Algorithm 2 is chosen as $21\times 21$ Gaussian with standard deviation $\sigma=21/3$ and multiplied by a hamming window to mitigate the windowing effect. We simulate the diffraction of light occurring in QIS by blurring the ground truth image with a $5\times 5$ box kernel

Table II shows the reconstruction quality averaged on 100 images in the DIV2K evaluation dataset. We choose the color PSNR as a quality metric, which is calculated as $\mathrm{CPSNR}=10\log_{10}\left(\frac{1}{\mathrm{MSE}}\right)$ where $\mathrm{MSE}$ is the total mean squared error summed over the three color channels. Our CFAs achieve better quality for the crosstalk case. Visually, our CFAs obtain color images with less aliasing artifacts and better details as shown in Figure 14.

VII Conclusion

We presented a general and flexible optimization framework to design color filter arrays for QIS. Our framework unifies the CMOS-based color filter array designs and extends to QIS. The color filter arrays designed by our framework are robust to crosstalk between the primary color channels, robust to aliasing between the luminance and chrominance channels, and are robust to noise. We verified the designs both theoretically and numerically through extensive experiments. Our evaluation indicated the effectiveness of the framework in offering trade-offs between different design criteria.

Acknowledgment

The authors thank Professor Keigo Hirakawa (Univ. Dayton), Professor Eric Fossum (Dartmouth College), Jiaju Ma (GigaJot Tech Inc.) and Abhiram Gnanasambandam (Purdue Univ.) for many insightful discussions about CFA designs. This work is supported, in part, by the National Science Foundation under grants CCF-1763896 and CCF-1718007.

Appendix A

A-A Proof of Proposition 1

Since the luminance channel comprises only one baseband component in the frequency domain, the luminance sensitivity in the amplitude of this component, i.e.,

[TABLE]

Substituting in the DFT equation with $u=v=0$, we get

[TABLE]

where $\boldsymbol{b}\overset{\text{def}}{=}\frac{1}{K}\boldsymbol{1}^{T}\boldsymbol{Z}_{l}$. As for the chrominance sensitivity $\gamma_{c}$, by squaring the definition in Equation (17), we get

[TABLE]

where (a) follows from Parseval theorem, and $\boldsymbol{Q}_{\alpha}\overset{\text{def}}{=}\boldsymbol{Z}_{\alpha}^{T}\boldsymbol{Z}_{\alpha}$ and $\boldsymbol{Q}_{\beta}\overset{\text{def}}{=}\boldsymbol{Z}_{\beta}^{T}\boldsymbol{Z}_{\beta}$ are two positive semidefinite matrices.