Double-Opponent Vectorial Total Variation

Abstract

We present a new vectorial total variation (VTV) method that addresses the problem of color consistent image filtering. Our approach combines insights based on the double-opponent cell representation in the visual cortex with state-of-the-art variational modelling using VTV regularization. Existing methods of vectorial total variation regularizers have insufficient (even no) coupling between the color channels and thus may introduce color artifacts. We address this problem by introducing a novel color channel coupling inspired from a pullback-metric from an opponent space to the observation space. We show existence and uniqueness of a solution in the space of vectorial functions of bounded variation. In experiments, we demonstrate that our novel approach compares favorably to state-of-the-art methods w.r.t. to structure coherence and color consistency.

These keywords were added by machine and not by the authors. This process is experimental and the keywords may be updated as the learning algorithm improves.

1 Introduction

Color image processing poses several challenges, since the notion of a “color edge” or a “color boundary” has no unique natural characterization. In this work we address the problem of adaptive color image filtering exploring biological findings in the visual cortex. We use these findings for image enhancement and design our new approach named double-opponent vectorial total variation.

Fig. 1.

Example of denoising with the proposed double-opponent vectorial total variation (right) of the noisy image (left).

The connection between observed visual stimuli and color space models is naturally described using tools from differential geometry. The model used in this work is inspired from recent findings in color experience models and the psychophysics of human color perception. Accordingly, we adopt a geometric viewpoint to explore the relation between color edges and a regularizer based on the color space geometry. Our presentation and treatment of the “geometry of color” leads us to a perceptually plausible model for color image enhancement via discontinuity preserving filtering. Figure 1 illustrates that our approach provides excellent recovery of color image data, retaining sharp edges, without introducing color artifacts.

Contributions. We propose a VTV-based regularizer, derived from a double-opponent space. We prove that the variational problem is convex and that its solution is unique and exists in the space of vectorial functions of bounded variation. Our experiments take into account multiple noise-levels and competing state-of-the-art denoising methods. We demonstrate improved structural coherence and improved color consistency.

2 Related Works

Color Space Representation. The choice of a color space representation is application dependent, yet without reaching any general consensus, so far. At the same time, modeling psychophysical effects of color in scientific applications is a highly non-trivial problem and many spaces have been proposed, e.g., RGB, sRGB, HSV, YPbPr and the YCbCr, CIELAB to mention few. It is by now widely accepted that image processing application in the RGB (Red, Green, Blue) space is suboptimal due to the high color channel correlation. The YCbCr and the YPbPr color spaces were introduced for analog and digital television transmission, respectively [1]. The HSV (Hue, Saturation, Value) color space was developed in the 1970’s for applications related to color display systems, and largely influenced by that time’s computer display systems [2]. The sRGB system was proposed for consistent image rendering over a wide range of imaging devices [3]. The CIELAB color space was proposed to yield a color space which is perceptually uniform. However, as shown by [4] (and references therein) this is not the case. This motivates the use of non-Euclidean metrics, even in supposedly perceptually uniform spaces. In this work we consider the double-opponent space.

Double-Opponent Color Representation. The double-opponent color space is thought to describe the representation of color in the human visual cortex, see [5–8]. Therefore, and due to its geometric structure, it is of great interest to investigate this color space for image enhancement applications. Previous works using this color space include, e.g., [9, 10]. Let $\varvec{u} = (r,g,b)^\top$ denote the red, green and blue color components of the RGB (observation) space, then the mapping from the observation space to the double-opponent space is given by the linear mapping $\varvec{Ou} : \mathbb {R}^{3} \rightarrow \mathbb {R}^{3}$ where

$$\begin{aligned} \varvec{O} = \begin{pmatrix} 1/\sqrt{3} &{} 0 &{} 0 \\ 0 &{} 1/\sqrt{6} &{} 0 \\ 0 &{} 0 &{} 1/\sqrt{2} \end{pmatrix} \begin{pmatrix} 1 &{} 1 &{} 1 \\ 1 &{} 1 &{} -2 \\ 1 &{} -1 &{} 0 \end{pmatrix} . \end{aligned}$$

(1)

The matrix $\varvec{O}\varvec{u} = (o_1, o_2, o_3)^\top$ produces a rotation and scaling of the RGB coordinate system. The opponent component $o_1$ is nothing else than the gray-scale value, $o_2$ is the subtraction of blue from yellow (mixing red and green equals yellow), and the last component $o_3$ is the subtraction of green from red. The components $o_2$ and $o_3$ consists of the image chroma and further decomposition of $o_2,o_3$ yields the corresponding hue (the color, red, green, yellow, etc. expressed as an angle) and saturation (colorfulness).

Physiological studies have shown the existence of double-opponent cell structures which are orientation-selective w.r.t. color discrimination and color boundaries [11]. Therefore, by preserving color discontinuities in this color space, we hypothesize that color borders trigger the activation of these double-opponent cells and thus yields the perception of crisp color edges in the image. This motivates the use of discontinuity preserving filtering introduced next.

Discontinuity Preserving Filtering. Let $\varvec{u} : \varOmega \rightarrow \mathbb {R}^M$ and $\varvec{u} = (u_1, ..., u_M)^\top$ where M is the image dimension and $\varOmega \subset \mathbb {R}^{2}$ is the image domain. The total variation (TV) for scalar-valued functions ($M=1$) [12] is defined as

$$\begin{aligned} {{\mathrm{TV}}}(u) := |\nabla u|(\varOmega ) = \int _\varOmega | \nabla u |, \end{aligned}$$

(2)

where $\nabla u$ denotes the generalized gradient which is a $\mathbb {R}^{2}$-valued Borel measure, $|\nabla u|$ the total variation of this measure, and ${{\mathrm{TV}}}(u)$ the total mass. We refer to [13, Ch. 10] for a detailed account. The function u belongs to the space of bounded variation $\mathrm {BV}(\varOmega )$ if $u \in L^1(\varOmega )$ and ${{\mathrm{TV}}}(u) < +\infty$. In this case, the following representation is valid, with $\mathcal {C}^1_c(\varOmega ,\mathbb {R}^2)$ denoting the space of compactly supported, continuously differentiable vector fields:

$$\begin{aligned} {{\mathrm{TV}}}(u) = \sup _{\Vert \varvec{\varphi }\Vert _{\infty } \le 1} \left\{ \int _\varOmega u {\text {div}} \left( \varvec{\varphi } \right) d\varvec{x} :\varvec{\varphi } \in \mathcal {C}^1_c(\varOmega ,\mathbb {R}^2) \right\} . \end{aligned}$$

(3)

Equation (3) shows that ${{\mathrm{TV}}}(u)$ is a support function in the sense of convex analysis. In particular, it is convex and lower semicontinuous.

Extending the scalar TV to color images is a non-trivial problem. If a color edge is insufficiently preserved in the smoothing process, artificial colors may emerge at the smooth transition between these colors as demonstrated by Fig. 2. The same figure also illustrates the problem of color shimmering, i.e., insufficient filtering of homogeneous regions.

Fig. 2.

Typical problems in color image denoising. (a) Introduction of artificial colors at edges due to insufficient color channel coupling. (b) Color shimmering due to insufficient smoothing. Our approach (right images) exhibits neither of these drawbacks. (Color figure online)

Vectorial TV. One of the first generalizations of TV to the vector-valued case was done by Blomgren and Chan [14]. They proposed the VTV measure

$$\begin{aligned} {{\mathrm{TV}}}_{BC}(\varvec{u}) = \left( \sum _{i=1}^M {{\mathrm{TV}}}(u_i)^2 \right) ^{1/2}, \end{aligned}$$

(4)

which sums up the squared scalar TV (2) over the image channels. However, due to no handling of the RGB-space intra-channel correlation, this model produces significant color smearing artifacts and therefore does not give color consistent filtering results [15].

Riemannian Geometry. Sapiro and Ringach [16] observed that the first fundamental form in the given (Riemannian) metric signals the presence of color edges. They proposed a regularizer ${{\mathrm{TV}}}_{SR}$ with an integrand $\sqrt{\lambda _+ + \lambda _-}$ and $\lambda _+ > \lambda _- \ge 0$ are the eigenvalues of the metric tensor. Building on this framework, Goldluecke and Cremers [17] introduced a vectorial total variation regularizer ${{\mathrm{TV}}}_J$ based “on the largest singular value of the derivative matrix” and show the existence of a solution. In an altogether geometric setting Sochen et al. [18] introduced the Beltrami framework. In this approach, they consider images as embedding maps between Riemannian manifolds characterized by the Polyakov action [18–20]. Although there is a coupling between the color channels, these approaches produce undesired artifacts such as color shimmering in homogeneous regions.

Dual VTV. Bresson and Chan [21] proposed a color total variation formulation that naturally extends the dual TV formulation (3) to the vectorial case. Based on the work of Chambolle [22] and Fornasier and March [23], Bresson and Chan presented a coherent framework for vectorial total variation together with a study of well-posedness. The dual VTV is defined as follows: Let $\varvec{u} \in L^1(\varOmega ;\mathbb {R}^M)$ and $\varvec{\xi } = (\varvec{\xi }_{1},\cdots ,\varvec{\xi }_{M})$ with ${\varvec{\xi }}_{i} \in \mathcal {C}^1_c(\varOmega ,\mathbb {R}^2),\, \forall i$. Then

$$\begin{aligned} \mathrm {VTV}(\varvec{u}) = \sup _{\Vert \varvec{\xi }\Vert _{F} \le 1} \left\{ \int _\varOmega \langle \varvec{u}, {\text {Div}} \left( \varvec{\xi } \right) \rangle \; d\varvec{x} \right\} , \end{aligned}$$

(5)

where $\Vert \cdot \Vert _{F}$ denotes the Frobenius norm and with the extension ${\text {Div}} \left( \varvec{\xi } \right) :=({\text {div}} \left( \varvec{\xi }_{1} \right) ,\cdots ,{\text {div}} \left( \varvec{\xi }_{M} \right) )^{\top } \in \mathbb {R}^{M}$ of the divergence operator. The VTV approach, however, still exhibits some color smearing as it does not properly take into account the color channel coupling.

Weighted VTV. One of the most recent approaches to incorporate color into a TV formulation was presented by Ono and Yamada [9]. They define a mixed $\ell _{1,2}$-norm where the intensity and chroma, obtained via the transformation (1), are weighted independently. However, it should be noted that the subspace defined by the chroma is not decorrelated but actually consists of the components hue and saturation. As a consequence, their framework does not treat the non-uniformity of the opponent space. The implication is that direct regularizing on the chroma via an Euclidean distance metric violates the non-Euclidean structure of this opponent space. Furthermore, it is easy to construct scenarios where the image saturation changes independently, and thus further motivates that chroma should be decomposed into hue and saturation [24].

3 Geometry of the Double-Opponent Space

In this section we consider the double-opponent space geometry and show how to extract the color information.

Connecting Observation and Double-Opponent Space. Recall that we denote by the linear mapping $\varvec{O} :\mathbb {R}^{3} \rightarrow \mathbb {R}^{3},\quad \varvec{u} = (r,g,b)^{\top } \mapsto \varvec{O u} = \varvec{o} = (o_{1}, o_{2}, o_{3})^{\top }$, with $\varvec{O}$ defined by (1), the transformation from the observation (RGB) color space to the double-opponent space. This linear mapping has full rank. The non-linear mapping to the hue (h), saturation (s) and lightness (L) representation of the opponent space is given by

$$\begin{aligned} \varvec{\psi } :\mathbb {R}^{3} \rightarrow \mathbb {R}^{3},\quad \varvec{o} \mapsto \varvec{c} =(L,h,s)^{\top }, \end{aligned}$$

(6)

where $L = o_{1}, h = \arctan (o_{2}/o_{3} ), s = \left\| ( o_{2}, o_{3}) \right\|$. We further set $\varvec{\varphi } :\varvec{u} \rightarrow \varvec{\varphi }(\varvec{u}) := \varvec{\psi }(\varvec{O u}) = (L,h,s)^{\top }$. Then, the Euclidean inner product $\langle \cdot ,\cdot \rangle$ on the Lhs-space induces via $\varvec{\varphi }$ the pullback metric on the RGB-space by

$$\begin{aligned} \langle \varvec{u}_{1}, \varvec{u}_{2} \rangle _{\varvec{u}}&:= \langle D\varvec{\varphi }(\varvec{u}) \varvec{u}_{1}, D\varvec{\varphi }(\varvec{u}) \varvec{u}_{2} \rangle = \langle \varvec{u}_{1}, \varvec{G}(\varvec{u}) \varvec{u}_{2} \rangle , \nonumber \\ \varvec{G}(\varvec{u})&:= \big (D\varvec{\varphi }(\varvec{u})\big )^{\top } D\varvec{\varphi }(\varvec{u}), \end{aligned}$$

(7)

where $\varvec{G}(\varvec{u})$ is the double-opponent metric tensor and $D\varvec{\varphi }(\varvec{u})$ is the Jacobian. Strictly speaking, we regard the RGB space as a linear Riemannian manifold $\mathcal {M}$ equipped with the above metric, which is an inner product on the tangent space $T_{\varvec{u}}\mathcal {M}$ that smoothly varies with $\varvec{u} \in \mathcal {M}$. Since every tangent space $T_{\varvec{u}}\mathcal {M}$ can be identified with $\mathcal {M}$, however, it makes sense to regard the Riemannian metric as inner product defined on the space itself. We refer, e.g., to [25] for background and further details.

Extracting Color Information. Next we compute the metric tensor $\varvec{G}(\varvec{u})$ and its eigendecomposition later used to define our novel color channel coupling. Setting

$$\begin{aligned} \varvec{\alpha }&= \varvec{\alpha }(\varvec{u}) = (b-g, r-b, g-r)^{\top }, \end{aligned}$$

(8)

$$\begin{aligned} \varvec{\beta }&= \varvec{\beta }(\varvec{u}) = (b+g-2 r,b+r-2g, r+g-2b)^{\top }, \end{aligned}$$

(9)

one easily verifies the relations $\varvec{\alpha } \perp \varvec{\beta }, \ \varvec{u}\perp \varvec{\alpha }, \ \langle \varvec{u}, \varvec{\beta } \rangle = \Vert \varvec{\alpha }\Vert ^{2}, \ \Vert \varvec{\beta }\Vert ^{2} = 3 \Vert \varvec{\alpha }\Vert ^{2}.$ In the subsequent analysis we will return to the following decomposition

$$\begin{aligned} \Vert \varvec{\alpha }\Vert ^{2} = \varvec{u}^\top \varvec{Pu}, \; \varvec{P} = \left( \begin{array}{ccc} 2 &{} -1 &{} -1 \\ -1 &{} 2 &{} -1 \\ -1 &{} -1 &{} 2 \\ \end{array} \right) , \end{aligned}$$

(10)

where $\varvec{P}$ is a symmetric and positive semi-definite matrix.

The Jacobian of the mapping $\varvec{\varphi }$ reads

$$\begin{aligned} D\varvec{\varphi }(\varvec{u}) = \dfrac{1}{\sqrt{3}} \bigg (\mathbbm {1}, 3 \frac{\varvec{\alpha }}{\Vert \varvec{\alpha }\Vert ^2}, \dfrac{\varvec{\beta }}{\Vert \varvec{\alpha }\Vert } \bigg )^{\top } , \end{aligned}$$

(11)

and the corresponding metric tensor (7) is

$$\begin{aligned} \varvec{G}(\varvec{u}) = \frac{1}{3} \left( \varvec{I} + \frac{9}{\Vert \varvec{\alpha }\Vert ^{4}} \varvec{\alpha } \varvec{\alpha }^{\top } + \frac{1}{\Vert \varvec{\alpha }\Vert ^{2}} \varvec{\beta } \varvec{\beta }^{\top } \right) , \end{aligned}$$

(12)

where $\varvec{G}(\varvec{u})$ has non-normalized eigenvectors $\mathbbm {1}, \varvec{\alpha }, \varvec{\beta }$ and eigenvaluematrix

$$\begin{aligned} \varvec{\varLambda } = \frac{1}{3} \varvec{I} + {{\mathrm{diag}}}\left( 0, \frac{3}{\Vert \varvec{\alpha }\Vert ^{2}}, 1 \right) . \end{aligned}$$

(13)

Saprio and Ringach [16] also derived the first fundamental form (i.e., $D\varvec{\varphi }(\varvec{u})$), yet in the Euclidean space, and they concluded that the tensor’s eigenvalues capture the color edge information. By contrast, we adopt the inverse principal directional change obtained from the eigendecomposition of the double-opponent metric tensor. Just as in the case of Saprio and Ringach, the interpretation of the decomposition is that a large eigenvalue of the tensor indicates the presence of image color change. Next we confirm this statement while investigating the information encoded in $\Vert \varvec{\alpha } \Vert ^2$ which will constitute the basis for our filtering scheme presented in Sect. 4.

Fig. 3.

Detected color structure in real images extracted by $\Vert \varvec{\alpha }\Vert ^2$. In these examples, primary colors such as red, green and blue and the opponent color yellow are well characterized. (Color figure online)

Encoded Information. The function $\Vert \varvec{\alpha }\Vert ^2$, given in (10), represents the principal change of color. To illustrate this, Fig. 3 shows the corresponding response for few natural images. To further understand $\Vert \varvec{\alpha }\Vert ^2$ we exploit the decomposition

$$\begin{aligned} \gamma (\varvec{u})&:=\Vert \varvec{\alpha }\Vert ^2 = (b - r)^2 + (r-g)^2 + (g - b)^2 \nonumber \\&= \Vert \varvec{Q}\varvec{u} \Vert ^2, \qquad \varvec{Q} = \begin{pmatrix} 1 &{} -1 &{} 0 \\ 0 &{} 1 &{} -1 \\ -1 &{} 0 &{} 1\\ \end{pmatrix}. \end{aligned}$$

(14)

Note that $\varvec{P} = \varvec{QQ}^\top$ (cmp. (10)). The coefficients of $\varvec{Qu}$ have previously appeared in an early work by Chambolle [26]. Chambolle defined a PDE with a directional diffusivity orthogonal to $(g - b)\nabla r + (b - r)\nabla g + (r - g)\nabla b$. However, as noted by Sapiro and Ringach [16], if two channels are equiluminant and if the third channel has an edge, this edge will remain unaffected by the filter. We remedy this drawback by considering a two component regularizer introduced in the next section.

Proposition 1

The function, $\gamma$ in (14), has the properties $\mathrm {(P1)} \; \gamma (\varvec{u} + c \mathbbm {1} ) = \gamma (\varvec{u})$ and $\mathrm {(P2)} \; \gamma (c \varvec{u}) = c^2 \gamma (\varvec{u})$, where c is a constant.

Proof

The result follows immediately from (14).

The above result yields the following interpretation of the $\gamma$-function: (a) (P1) shows that $\gamma$ is invariant to intensity shifts. (b) (P2) shows that $\gamma$ has a quadratic dependency on intensity changes. (c) it follows from (a) that $\gamma$ depends on color changes and (d) it follows from (b) that $\gamma$ depends on color shifts. Under constant intensity, $\gamma$ captures change of color as illustrated in Fig. 4 (a). In this figure we show equiluminant discs at constant intensity along with the corresponding response of $\gamma$. It is clearly visible that $\gamma$ describes the structure of the color change as there is a stronger response for highly saturated colors. In the lower half of the intensity range we predominantly detect the primary colors red, green and blue. As the intensity increases $\gamma$ shows primary responses from yellow, cyan and magenta. The intensity axis is located in the center of these discs and, as expected, we do not obtain a value of colorfulness.

Fig. 4.

(a) Color discs with corresponding response of $\gamma$, (14), in (b). The largest magnitude (red color) is obtained at the primary colors (red, green and blue) and the opponent colors (yellow, cyan and magenta). As expected, the response on the intensity axis (center of discs) is 0 (black). (c) Interpretation of the vector $r-g$ as an orthogonal component to yellow. (Color figure online)

The geometric interpretation of $\gamma$ is illustrated as an example via the $r-g$-component. The other two color difference terms follow with similar reasoning. We know that the color yellow, y, is composed as a sum of red and green, i.e., $y=r+g$, and written in vector form we have $r-g = (1, -1, 0)^\top$ and $y = r + g = (1, 1, 0)^\top$. We see that yellow is perpendicular to the difference $r-g$, i.e., $y \perp (r-g) = 0$. This is illustrated in Fig. 4 (b). Analogous argument hold for the other terms of $\gamma$, i.e., $b-r$ is orthogonal to magenta, and $g-b$ orthogonal to cyan. In this way $\gamma$ covers the RGB space. Moreover, as $\gamma$ describes the color structure, preserving its edge information prevents color distortion in the filtering process. Based on this analysis, we are now prepared to introduce the double-opponent VTV regularization term.

4 Double-Opponent Vectorial Total Variation

The Energy. Channel-by-channel filtering of the RGB space is prone to introduce color artifacts [14, 21]. On the other hand, purely decorrelating the color channels without considering the geometry is also sub-optimal, see e.g., [9, 16]. We propose a two-component regularizer: one component performs channel-by-channel filtering penalizing all intra-channel content and one component which explicitly targets the color information. The color specific prior, $J_{OPP}$, defines a natural inter-channel coupling from the geometry of the double-opponent space. Let $\varvec{g}$ be the observed (noisy) image data. The energy we introduce is

$$\begin{aligned}&\min _{\varvec{u}} \left\{ E(\varvec{u}) = \frac{\mu }{2}\Vert \varvec{u}-\varvec{g}\Vert ^2_{L^2(\varOmega )} + R(\varvec{u}) \right\} \end{aligned}$$

(15a)

$$\begin{aligned}&R(\varvec{u}) = \alpha \sum _{i=1}^3 {{\mathrm{TV}}}(u_i) + \beta J_{OPP}(\varvec{u}), \end{aligned}$$

(15b)

where $\mu ,\alpha ,\beta > 0$.

Furthermore, let $\nabla \varvec{u}= (\nabla u_1, \nabla u_2, \nabla u_3) : \varOmega \rightarrow \mathbb {R}^{2\times 3}$ be the vectorial gradient of an color image ($M=3$) in the generalized sense as discussed in connection with Eq. (2), and $\varvec{p} = (\varvec{p}_1, \varvec{p}_2, \varvec{p}_2) \in C_{c}^{1}(\varOmega ; \mathbb {R}^{2 \times 3})$ with ${\text {Div}} \left( \varvec{p} \right) = ({\text {div}} \left( \varvec{p}_1 \right) , {\text {div}} \left( \varvec{p}_2 \right) , {\text {div}} \left( \varvec{p}_3 \right) )^{\top }$.

Definition 1

(Double-Opponent VTV). The double opponent regularizer is defined as

$$\begin{aligned} J_{OPP}(\varvec{u})&:= \int _\varOmega \Vert \nabla \varvec{Qu} \Vert \end{aligned}$$

(16a)

$$\begin{aligned}&:= \sup _{\Vert \varvec{p}\Vert _{\infty }\le 1} \left\{ \int _\varOmega \langle \varvec{Qu}, {\text {Div}} \left( \varvec{p} \right) \rangle \; d\varvec{x} \right\} \end{aligned}$$

(16b)

where $\Vert \varvec{p}\Vert _{\infty } = \max \{|\varvec{p}_{1}|,|\varvec{p}_{2}|,|\varvec{p}_{3}|\}$.

Remark 1

Note that although this looks like an anisotropic TV formulation, it does incorporate a proper coupling of the channels through the matrix $\varvec{Q}$.

Theorem 1

(Invariance and Convexity). $J_{OPP}$ is rotationally and intensity invariant, 1-homogeneous and convex.

Proof

Rotational invariance follows from the isotropy of the feasible set of the dual variable $\varvec{p}$, that is $\Vert \varvec{p}\Vert _{\infty }=\Vert (\varvec{p}_{1},\varvec{p}_{2},\varvec{p}_{3})\Vert _{\infty } \le 1 \implies \Vert (\varvec{R} \varvec{p}_{1},\varvec{R} \varvec{p}_{2},\varvec{R} \varvec{p}_{3})\Vert _{\infty } \le 1$, for any orthogonal matrix $\varvec{R}$. As a consequence of property (P1) and (P2) of Proposition 1, $J_{OPP}$ is invariant to intensity shifts, and the relation $J_{OPP}(c \varvec{u}) = c J_{OPP}(\varvec{u})$ is immediate, for any positive constant $c > 0$. Finally, convexity follows from the definition of $J_{OPP}$ as pointwise supremum of affine functions.

Existence of Solution. Next we show that the variational approach (15a) is well posed.

Lemma 1

(Bounded Variation). Let $\varvec{u} \in BV(\varOmega ;\mathbb {R}^3)$ then $\varvec{Qu} \in BV(\varOmega ;\mathbb {R}^3)$.

Proof

If $\varvec{u}$ is in $L^{1}(\varOmega ;\mathbb {R}^{3})$, then so is $\varvec{Q u}$, because $\varvec{Q}$ is a constant matrix. Furthermore, expanding the bilinear form under the integral of (16b) results in a linear combination of terms of the form (3), which are finite by the assumption $\varvec{u} \in BV(\varOmega ;\mathbb {R}^3)$.

As a consequence, the objective function $E(\varvec{u})$ (15a) is well defined. We next show that there is a unique color image $\varvec{u}$ minimizing $E(\varvec{u})$.

Theorem 2

(Uniqueness and Existence of Solution). Let $\varvec{g} \in L^\infty (\varOmega ,\mathbb {R}^3)$ and $\varvec{u} \in BV(\varOmega ,\mathbb {R}^3)$. Then there exists a unique minimizer $\varvec{u}^{*}$ of $E(\varvec{u})$ in (15a).

Proof

We adapt and sketch a standard proof pattern from [13]. Due to $\varvec{g} \in L^\infty (\varOmega ,\mathbb {R}^M)$, we may assume that all admissible $\varvec{u}$ are uniformly bounded in the sense that $|u_{i}(x)| \le \Vert g_{i}\Vert _{L^{\infty }(\varOmega )},\, i=1,2,3,\, \forall x \in \varOmega$. Let $(\varvec{u}_{n})_{n \in \mathbb {N}}$ be a minimizing sequence with respect to $E(\varvec{u})$. Then, after passing to a subsequence $(\varvec{u}_{n_{k}})_{k \in \mathbb {N}}$, there exists a $\varvec{u}^{*} \in BV(\varOmega ;\mathbb {R}^3)$ with $\varvec{u}_{n_{k}} \rightarrow \varvec{u}^{*}$ strongly in $L^{1}_{loc}(\varOmega ;\mathbb {R}^{3})$, $\nabla (u_{i})_{n_{k}} \rightarrow \nabla u^{*}_{i}$ in an appropriate weak sense, and $J_{opp}(\varvec{u}_{n_{k}}) \rightarrow J_{opp}(\varvec{u}^{*})$ in view of Lemma 1. It follows from Fatou’s lemma and the lower-semicontinuity of $E(\varvec{u})$ that $\varvec{u}^{*}$ minimizes $E(\varvec{u})$, whereas uniqueness of $\varvec{u}^{*}$ is a consequence of the strict convexity of $E(\varvec{u})$ due to the data term of (15a).

Next, we derive an efficient numerical scheme which optimize our novel energy.

Discretization. With slight abuse of notation, we denote again by $\varvec{u}, \varvec{g} \in \mathbb {R}^{3 N}$ the discretized representations of $\varvec{u}$ and $\varvec{g}$ as column vectors where the image channels are stacked in the order $u_1, u_2, u_3$. N is the number of pixels in one channel. We define the discrete image gradient for one channel as $\varvec{D}_1 = \begin{bmatrix} \varvec{D}_x \\ \varvec{D}_y \end{bmatrix} \in \mathbb {R}^{2N\times N}$, $\varvec{D}_x, \varvec{D}_y \in \mathbb {R}^{N\times N}$ and subscript denotes the forward finite difference operator in x and y directions, respectively. Furthermore, we let $q \in \mathbb {N}^+$ s.t. $\varvec{I}_q \in \mathbb {R}^{q\times q}$ denotes the identity matrix. In this notation, the three channel derivative matrix for a color image is $\varvec{D} = \varvec{D}_1 \otimes \varvec{I}_3 \in \mathbb {R}^{6N\times 3N}$ where $\otimes$ is the Kronecker product. Then $\varvec{D}\varvec{u} : \mathbb {R}^{3N} \rightarrow \mathbb {R}^{6N}$ is the derivative for the three channels. The discrete representation of the Laplacian is denoted $\varvec{L} = (\varvec{D}_1^\top \varvec{D}_1) \otimes \varvec{I}_3 \in \mathrm {R}^{3N\times 3N}$ such that $\varvec{L} \varvec{u} : \mathbb {R}^{3N} \rightarrow \mathbb {R}^{3N}$. The channel coupling matrix for the discretized image is denoted as $\varvec{C} = \varvec{Q} \otimes \varvec{I}_N \in \mathbb {R}^{3N \times 3N}$.

Optimization. With the notation introduced above, we write the corresponding discretized form of (15a) as

$$\begin{aligned} \min _{\varvec{u},\varvec{d},\varvec{e}}\; \frac{\mu }{2}\Vert \varvec{u}-\varvec{g}\Vert ^{2}_2 + \alpha \Vert \varvec{d}\Vert _{1} + \beta \Vert \varvec{e}\Vert _{1} \quad \text {s.t.}\quad \varvec{d} = \varvec{D u},\quad \varvec{e}=\varvec{D C u}. \end{aligned}$$

(17)

To minimize this objective function we may use any standard optimization technique, e.g., [27–30]. Here we adopt the Split Bregman approach [31] as it yields a simple numerical scheme. In the following, we let $\Vert \varvec{v} \Vert _{\varvec{W}}^2 := \langle \varvec{v}, \varvec{W}\varvec{v} \rangle$ be a weighted Euclidean norm and set

$$\begin{aligned} \varvec{B}(\varvec{u},\varvec{b},\varvec{d},\varvec{e}) := \begin{pmatrix}\varvec{d} \\ \varvec{e} \end{pmatrix}-\begin{pmatrix}\varvec{D} \\ \varvec{D C} \end{pmatrix}\varvec{u} - \varvec{b}. \end{aligned}$$

(18)

Applying the Split Bregman approach yields the iteration

$$\begin{aligned} \nonumber&(\varvec{u}^{k+1},\varvec{d}^{k+1},\varvec{e}^{k+1}) \end{aligned}$$

(19a)

$$\begin{aligned}&= \min _{\varvec{u},\varvec{d},\varvec{e}}\; \frac{\mu }{2}\Vert \varvec{u}-\varvec{g}\Vert ^{2}_2 + \Vert \varvec{d}\Vert _{1} + \Vert \varvec{e}\Vert _{1} + \frac{1}{2} \left\| \varvec{B}(\varvec{u},\varvec{b},\varvec{d},\varvec{e}) \right\| ^{2}_{\varvec{W}}, \end{aligned}$$

(19b)

$$\begin{aligned} \varvec{W}&= \begin{pmatrix} \alpha \varvec{I}_{6N} &{} 0\\ 0 &{} \beta \varvec{I}_{6N} \end{pmatrix}, \quad \alpha ,\beta > 0 \\ \varvec{b}^{k+1}&= \varvec{b}^{k} + \begin{pmatrix}\varvec{D} \\ \varvec{D C} \end{pmatrix}\varvec{u}^{k+1} - \begin{pmatrix}\varvec{d}^{k+1} \\ \varvec{e}^{k+1} \end{pmatrix}. \end{aligned}$$

(19c)

The former problem is solved iteratively by

$$\begin{aligned} \varvec{u}^{k+1}&= \min _{\varvec{u}}\;\frac{\mu }{2}\Vert \varvec{u}-\varvec{g}\Vert ^{2}_2 + \frac{1}{2} \left\| \varvec{B}(\varvec{u},\varvec{b}^k,\varvec{d}^k,\varvec{e}^k) \right\| ^{2}_{\varvec{W}} \end{aligned}$$

(20)

followed by two shrinkage updates for $\varvec{d}^{k+1},\varvec{e}^{k+1}$. Regarding, subproblem (20), we set $\varvec{b} = [\varvec{b}_1^\top , \; \varvec{b}_2^\top ]^\top$, $\varvec{b}_1,\varvec{b}_2 \in \mathbb {R}^{6N}$ for notational convenience and obtain the update step

$$\begin{aligned} ((\mu \varvec{I}&+ \alpha \varvec{I})\varvec{L} + \beta \varvec{C}^\top \varvec{L}\varvec{C} )\varvec{u}^{k+1} \nonumber \\&= \mu \varvec{g} + \alpha \varvec{D}^\top ( \varvec{d}^k - \varvec{b}^k_1) + \beta (\varvec{D} \varvec{C})^\top (\varvec{e}^k-\varvec{b}^k_2). \end{aligned}$$

(21)

Our experiments confirm the observation of [31] that only computing an approximate solution accelerates the overall iterative scheme without compromising convergence. Consequently, we merely apply few conjugate gradient iterative steps to compute $\varvec{u}^{k+1}$. This is computationally cheap since all matrices involved are sparse.

Finally, we update $\varvec{b}^{k+1}$ according to (19c) and iterate all steps until $\Vert \varvec{u}^{k}-\varvec{u}^{k+1}\Vert ^2_2/\Vert \varvec{u}^{k+1}\Vert ^2_2 < 0.9\sqrt{3N\sigma ^2}/255^2$ and $\sigma$ is the noise level standard deviation.

5 Experiments

Setup. The experimental evaluation use all 100 images of the Berkeley validation dataset [34]. The image data is normalized to the range [0, 1] from an 8-bit representation. In addition to a qualitative evaluation we include the peak signal-to-noise ratio (PSNR), the structural similarity index (SSIM) [35] and the CIEDE 2000, a measure of color consistency [36]. We optimize the parameter settings in a given feasible range for each method, image and noise level with respect to the best obtained SSIM value. The following methods and parameter ranges are included in the evaluation and we refer to the respective works for further details:

• Decorrelated VTV [9] (DVTV): Search space for optimal parameter configuration is $\tau \in \{0.95, 1, 1.05\}$, $w \in \{0.3, 0.4, 0.5, 0.6, 0.7\}$.

• Primal-dual VTV [21] (PDVTV): The regularization parameter was optimized for 5 uniformly sampled values in the range $10^{-3}$ to 0.2.

• Double Opponent VTV (ours) (OVTV): Optimized parameter space of $\mu$ are 5 uniformly sampled values from 1 / 255 to 30 / 255, $\alpha = 1$ and $\beta$ was uniformly sampled from 5 values in the range 1 / 255 to 5 / 255.

• Total generalized variation [33] (TGV): Applied componentwise and only included for comparison. The regularization parameter was uniformly sampled with 5 values in the range $10^{-3}$ and 0.25.

• Color BM3D [32] (BM3D): Standard deviation of the additive Gaussian noise was given as input.

Fig. 5.

Visual comparison of the double-opponent regularization and influence of $\mu$ and $\beta$ for fixed $\alpha =1$. Error measures are SSIM/CIEDE. As $\mu$ increases details are oversmoothed. $\beta$ influence the double-opponent term and we see (in this example) the best performance at $\beta =10/255$. (Color figure online)

Impact of Parameters. We check the sensitivity of different parameter settings for our double-opponent regularizer in relation to the dataterm in Fig. 5. The image was corrupted with standard deviation 20 of additive Gaussian noise and in this example we consider the visual quality and the color consistency measured with the CIEDE measure (lower value is better) and SSIM value (higher is better). The first row illustrates that the noise is not accurately removed as there is considerable amount of speckle-effect, i.e., the regularizer treats noise as structure. Increasing the influence of the data term improves the visual quality significantly and we do not observe any color shimmer. In the last row $\mu$ is clearly too large as image details are oversmoothed.

Color Denoising. Common artifacts in color image denoising is color shimmer in homogeneous regions and the introduction of artificial colors. Examples of these artifacts were given in Fig. 2. In this work, we introduce a challenging image recovering problem where current state of art denoising algorithms consistently show poor performance. Rather than corrupting all image data with additive noise, we corrupt the color components ($\varvec{o}_2$ and $\varvec{o}_3$) with 20, 50 or 80 standard deviations of Gaussian noise and ignore the intensity channel. Remarkably, after transforming from the (now noisy) opponent representation to the RGB space, one can show that the r, g, b components retain a Gaussian noise distribution of zero mean but with $\sqrt{2/3}$ scaled standard deviation. For this reason we also evaluate BM3D (our main competing method) which only requires an accurate estimation of the image noise, with this scaled noise variance, we denote this with

• Color BM3D with scaled noise distribution (BM3DS)

Fig. 6.

Visual comparison of the compared methods and corresponding error values. The result of our OVTV produces the most accurate result, only marginally beaten by BM3D in terms of color accuracy. Yet, the visual quality of OVTV is more clear and does not suffer from desaturated colors as in DVTV. (Color figure online)

Table 1. Error measures for the evaluated methods. The best values are marked in bold. For lower noise levels, BM3D is the only competitor for our OVTV. As the noise level increases OVTV consistently shows the best performance. See text for details.

Fig. 7.

In this image, accurate restoration of the surfer is highly non-trivial due to the large uniform and non-texturized regions surrounding it. Our OVTV approach produces sharp color borders which are similar to the original noise free image. Although the state of the art VTV method DVTV does restore the uniform background accurately the red color appears desaturated. PDVTV produces color shimmering. TGV and BM3D(S) both oversmooth the image. (Color figure online)

Fig. 8.

(left) Inpainting of 85 % missing data and (right) deconvolution examples. These illustrations show that the OVTV approach can accurately restore edges without introducing color artifacts also in related variational problems. (Color figure online)

Fig. 9.

Additional denoising results (ssim/psnr/ciede). The first and second row show recovery of data corrupted with standard deviation 50 and 80 of additive Gaussian noise. Our OVTV performs well for a wide range of images and noise levels.

Results. Table 1 shows the average error and standard deviation values for each method and noise level for the 100 Berkeley images. Our double-opponent approach (OVTV) compares well to DVTV and BM3D. Close-ups from few result images are given in Figs. 6 and 7. It is clearly visible that the our OVTV produces clean images, does not introduce artificial color artifacts, does not oversmooth details and does not suffer from color shimmering.

Figure 9 shows additional denoising results for standard deviation 50 (first row) and 80 (second row) of Gaussian noise. Although BM3D and DVTV show better color consistency in these examples, the structural coherence measured in PSNR is similar and SSIM is significantly improved for OVTV.

Inpainting and Deconvolution. We adopt our scheme to include inpainting and deconvolution in Fig. 8. As seen, homogeneous regions and edges are accurately restored without introducing color shimmering and artificial colors.

6 Conclusion

We have shown that the double-opponent theory can greatly improve the performance of VTV-based methods. Motivated by recent and classical results in color theory we let the mapping from the opponent-space to the observation space serve as a basis of our vectorial formulation. If the aim is consistent image filtering, then clearly the double-opponent VTV is preferable for image restoration.