Joint Reconstruction of Multi-channel, Spectral CT Data via Constrained Total Nuclear Variation Minimization
David S. Rigie and Patrick J. La Rivi\`{e
We explore the use of the recently proposed "total nuclear variation" (TNV) \cite{Rigie2014,Holt2014 as a regularizer for reconstructing multi-channel, spectral CT images. This convex penalty is a natural extension of the total variation (TV) to vector-valued images and has the advantage of encouraging common edge locations and a shared gradient direction among image channels. We show how it can be incorporated into a general, data-constrained reconstruction framework and derive update equations based on the first-order, primal-dual algorithm of Chambolle and Pock. Early simulation studies based on the numerical XCAT phantom indicate that the inter-channel coupling introduced by the TNV leads to better preservation of image features at high levels of regularization, compared to independent, channel-by-channel TV reconstructions. %In our first experiment, we directly reconstruct five energy-bin images from an ideal photon-counting system, using both our TNV joint reconstruction model and independent channel-by-channel TV reconstructions. In the second experiment, we consider an ideal, dual-kVp system and jointly reconstruct the low-energy, high-energy, bone, and soft-tissue images. This "hybrid" reconstruction technique represents a novel means for controlling noise in basis-material reconstructions and demonstrates that the TNV works equally well in the energy and material domains.
paper_draft4_revised.tex
\documentclass{article}
\usepackage{amsmath}
\usepackage{bm}
\usepackage[left=1in,top=1in,right=1in,bottom=1in]{geometry}
%\usepackage{showframe}
\usepackage{subcaption,graphicx}
\usepackage{titlesec}
\usepackage{comment}
\usepackage[section]{placeins}
\usepackage{amsfonts}
\usepackage{algorithm}% http://ctan.org/pkg/algorithms
\usepackage{algpseudocode}% http://ctan.org/pkg/algorithmicx
\usepackage{tikz}
\usepackage{csquotes}
\usepackage{caption}% http://ctan.org/pkg/caption
\usepackage{appendix}
\MakeOuterQuote{"}
\usepackage{todonotes}
\usepackage{ulem}
\usepackage{setspace}
\onehalfspacing


\title{\textbf{Joint Reconstruction of Multi-channel, Spectral CT Data via Constrained Total Nuclear Variation Minimization}}
\author{ David S. Rigie and Patrick J. La Rivi\`{e}re}

\DeclareMathOperator*{\argmin}{argmin}
\DeclareMathOperator*{\argmax}{argmax}

\begin{document}
\bibliographystyle{ieeetr}

\newcommand{\utrue}{\bm{u}_\mathrm{true}}
\newcommand{\epsref}{\epsilon_{\mathrm{ref}}}
\newcommand{\Grad}[1]{\mathrm{Grad}\,#1}
\newcommand{\Div}[1]{\mathrm{Div}\,#1}
\newcommand{\TV}{\mathrm{TV}}
\newcommand{\VTV}{\mathrm{VTV}}
\newcommand{\TNV}{\mathrm{TNV}}
\newcommand{\inner}[1]{\left<#1\right>}
\renewcommand{\div}[1]{\mathrm{div}\, #1}
\newcommand{\unoisy}{\bm{u}_{\mathrm{noisy}}}
\newcommand{\uu}{{\bm{u}}}
\newcommand{\uubar}{{\bar{\bm{u}}}}
\newcommand{\zz}{{\bm{z}}}
\newcommand{\qq}{{\bm{q}}}
\newcommand{\deltak}{\delta_{\mathcal{N}}}
\newcommand{\jth}{j^{\mathrm{th}}}
\newcommand{\ff}{\bm{f}}
\newcommand{\xx}{\bm{x}}
\newcommand{\yy}{\bm{y}}
\newcommand{\prox}[2]{\mathbf{prox}_{#1}\hspace{-2pt}\left(#2\right)}
\renewcommand{\vec}[1]{\mathrm{vec}\,#1}


\maketitle

\begin{abstract}
\textit{
We explore the use of the recently proposed "total nuclear variation" (TNV) \cite{Rigie2014,Holt2014} as a regularizer for reconstructing multi-channel, spectral CT images. This convex penalty is a natural extension of the total variation (TV) to vector-valued images and has the advantage of encouraging common edge locations and a shared gradient direction among image channels. We show how it can be incorporated into a general, data-constrained reconstruction framework and derive update equations based on the first-order, primal-dual algorithm of Chambolle and Pock. Early simulation studies based on the numerical XCAT phantom indicate that the inter-channel coupling introduced by the TNV leads to better preservation of image features at high levels of regularization, compared to independent, channel-by-channel TV reconstructions. %In our first experiment, we directly reconstruct five energy-bin images from an ideal photon-counting system, using both our TNV joint reconstruction model and independent channel-by-channel TV reconstructions. In the second experiment, we consider an ideal, dual-kVp system and jointly reconstruct the low-energy, high-energy, bone, and soft-tissue images. This "hybrid" reconstruction technique represents a novel means for controlling noise in basis-material reconstructions and demonstrates that the TNV works equally well in the energy and material domains.
}
\end{abstract}

\section{Introduction}
Since the influential work of Rudin, Osher, and Fatemi \cite{Rudin1992}, the total-variation (TV) has been used for a variety of image processing tasks, such as denoising \cite{Rudin1992}, deconvolution \cite{Chan1998}, and demosaicking\cite{Komatsu2005}. Because the TV is a nonsmooth penalty that disfavors oscillations while allowing for sharp discontinuities, it has been found to effectively suppress noise while preserving edges. The TV has more recently been rediscovered as a prior for tomographic inverse problems suffering from incomplete projection data \cite{Candes2006,Sidky2006}. Consequently, many works have investigated TV-regularized reconstruction algorithms for sparse-view, limited angle, and low-dose CT \cite{Bian2013,Song2007,Tian2011}. 

In last couple of years, several papers have proposed generalizations of the TV for multi-channel or vector-valued images, in order to extend popular image-processing techniques to the RGB color space \cite{Sapiro1997,Blomgren1998,Goldluecke2012,Lefkimmiatis2013,Holt2014}. Coincidentally, there is also a fast-growing interest in multi-spectral CT imaging, due to rapid developments in energy-resolving, photon-counting detectors. Since spectral CT is essentially just CT with multiple "color" channels, many of these theoretical advances can be directly applied to spectral CT reconstruction.

The primary consideration when formulating a vectorial TV (VTV) is how the color channels should be coupled. One could na\"{\i}vely apply the conventional TV to each channel independently, but this would fail to leverage the strong inter-channel correlations inherent to multi-spectral images. In this work, we advocate for a VTV that penalizes the nuclear norm of the Jacobian derivative because it has several desirable properties. This convex prior encourages all image channels to have common edge locations as well as for their gradient vectors to point in the same direction. The notion that gradient vectors of different color channels should be approximately aligned has been suggested previously \cite{keren1998,Holt2011,Rigie2014}. Additionally, we will demonstrate that the inter-channel coupling helps prevent edge smearing in the noisiest color channels, allowing for stronger regularization without significant feature loss. This specific VTV has been recently proposed by Lefkimmiatis \cite{Lefkimmiatis2013}, Rigie \cite{Rigie2014}, and Holt \cite{Holt2014}. We borrow from Holt in referring to it as the "Total Nuclear Variation" (TNV). 

The main objective of this paper is to demonstrate the usefulness of the TNV as a regularizer for multi-energy CT reconstruction. Our study contains several important contributions. To our knowledge, this is the first study with empirical results demonstrating the use of the TNV for spectral CT reconstruction, and our results indicate a clear advantage over independent, channel-by-channel reconstructions. Secondly, we propose a novel, noise-balancing transform that improves noise-suppression for multi-channel data for which the color channels have mismatched noise levels. This is a situation that is likely to occur in photon-counting CT, where the lowest and highest energy windows tend to suffer from relatively low count rates. Finally, we demonstrate that the TNV can be applied not only to the energy-binned data, but also to synthesized basis-material data. In most color image processing tasks, the measured RGB basis is the same as the observed basis. The analog in spectral CT would be to directly reconstruct the energy-bin data, so we will refer to this as the "color" basis. However, often the goal is, instead, to reconstruct density maps of different basis materials \cite{Fessler2002,Long2014}, which can be obtained through a non-linear transformation of the measured data \cite{Alvarez1976}. In this work, we will refer to this space as the "material" basis. We show empirically that the TNV works well in either basis and that it may actually be advantageous to co-reconstruct the color basis and the material basis. We refer to this as the "hybrid" basis, which is formed by concatenating the color and material basis channels. By simultaneously reconstructing the low-noise, color basis and the high-noise, material basis, the TNV is able to better suppress much of the noise that is amplified during the decomposition. 

First, we will review some theoretical concepts and define the total nuclear variation (TNV). Next, we will describe our data-constrained, TNV reconstruction model and show how the first-order, primal-dual algorithm of Chambolle and Pock \cite{Chambolle2011} can be used to derive efficient update equations. We compare our TNV reconstruction model to performing separate channel-by-channel reconstructions using simulated, dual-energy and photon-counting data using the anthropomorphic XCAT phantom \cite{Segars2010}. In the photon-counting study, we present reconstruction results from the color basis, where the color channels correspond to five equispaced energy windows. In the dual kVp study, we present reconstructions of the hybrid basis, where 80 kVp, 140 kVp, bone, and soft-tissue images are simultaneously reconstructed. Overall, we find that the TNV regularization results in less edge smearing at the same noise levels, while the computation time is very similar to the na\"{\i}ve channel-by-channel reconstruction. Additionally, the hybrid TNV reconstruction technique seems to be a promising method for controlling noise in the basis-material images.

\section{Theory}

\subsection{Notation and definitions}
In this work, we deal only with 2D images, but the generalization to 3D is straightforward. Consider a single-energy CT image vector $u \in\mathbb{R}^{\mathcal{L}}$, where $\mathcal{L}$ is the total number of image pixels.  We will also refer to its discrete derivative $Du\in \mathbb{R}^{2\cdot\mathcal{L}}$, which consists of a 2D derivative vector $\left[Du\right]_\ell$ at every pixel location. where $[.]_\ell$ denotes the a particular pixel-index. This quantity is also directly related to the discrete gradient by the following relationship: 
%
\begin{equation}
\left[Du \right]_\ell^T = \left[\nabla u \right]_\ell.
\end{equation}
%
 The quantities $u$ and $Du$ are both 1D vectors, so we will use $\|.\|$ and $\inner{.,.}$ to refer to the usual vector norm and inner product. We will also refer to a multi-channel image, $\bm{u} \in \mathbb{R}^{\mathcal{M}\cdot\mathcal{L}}$, with $\mathcal{M}$ spectral channels and its discrete derivative $D\bm{u} \in \mathbb{R}^{(\mathcal{M}\times 2)\cdot\mathcal{L} }$. The quantity $D\bm{u}$ is a discrete analog to the Jacobian derivative of a vector-function, and at every pixel location the first-order derivatives are fully characterized by the $\mathcal{M}\times 2$ rectangular matrix $\left[D\bm{u} \right]_\ell$. For two tensors $\bm{x},\bm{y} \in \mathbb{R}^{(\mathcal{M}\times 2)\cdot\mathcal{L}}$, we define the inner product $\inner{\xx,\yy} = \inner{\vec{\xx},\vec{\yy}}$, where $\mathrm{vec}(.)$ unfolds a tensor into a 1D vector. 

\begin{comment}
\subsection{Regularization in inverse problems}
\paragraph{Regularized form}
Typically, when solving inverse problems for an unknown image, the objective function will contain both a "data-fidelity" term and an adjustable regularization term. 
%
\begin{equation}\label{eq:regularized}
\min_u \quad F(u) + \lambda R(u)
\end{equation}
%
The data-fidelity $F(u)$ term encompasses the measurement model, capturing any relevant physical, geometric, or statistical considerations. The regularization term $R(u)$ may be used to incorporate prior knowledge we have about the object (e.g. positivity), encourage desired image properties (e.g. smoothness), improve the conditioning of the problem, and prevent overfitting to noisy data. 
\paragraph{Data-constrained form} 
A useful result from convex analysis demonstrates the equivalence of the aforementioned regularized form and a "data-constrained" form, which is given by 
%
\begin{equation}
\begin{split}
\min_u &\quad R(u)\\
\mathrm{s.t.} &\quad D(u)\leq \epsilon\,.
\end{split}
\label{eq:dataconstrained}
\end{equation}
%
The parameter $\epsilon$ should be adjusted based on the magnitude of noise and model error. Both Eqns. (\ref{eq:regularized}) and (\ref{eq:dataconstrained}) define the same set of solutions, but it is usually difficult to determine the mapping between $\lambda$ and $\epsilon$. Though the regularized form is often easier to solve, we will use the data-constrained form in this work because it provides an efficient and fair way to compare different regularization schemes. Specifically, we can fix the data-constraint parameter $\epsilon$, and compare the solutions for different regularization functions.

\end{comment}

\subsection{The scalar total variation}
Total variation regularized CT reconstruction algorithms have been studied extensively due to the approximate gradient sparsity of CT images. Many works have demonstrated that such schemes can often yield high quality images from severely undersampled projection data \cite{Candes2006,Sidky2006,Han2012,Bian2013,Xu2014}. 

\paragraph{The anisotropic TV}
The "anisotropic" TV is simply defined as the $\ell_1$-norm of the derivative of the image $u$,
%
\begin{equation}
\TV_a(u) = \|Du\|_1,
\end{equation}
%
where it is useful to think of the $\|.\|_1$ as a "sparsity-inducing" norm because it is a convex relaxation of $\|.\|_0$.
\paragraph{The isotropic TV}
Somewhat more common is the "isotropic" TV, which is defined in terms of the gradient-magnitude image,
%
\begin{equation}
\TV(u) = \sum_{\ell=1}^\mathcal{L} \|\left[D u\right]_\ell\|_2
\end{equation}
%
This penalty function has similar properties to the anisotropic TV, but is rotationally invariant, whereas the anisotropic TV tends to under-penalize horizontal and vertical edges. 

\subsection{The vectorial total variation}
Now, we consider generalizing the TV to the $\mathcal{M}$-channel image $\bm{u}$.
\paragraph{Na\"{\i}ve generalization to multi-channel images}
The simplest way to extend the definition of the TV to this vector-valued image, $\bm{u}$, would be to take the summation of the total-variation of each channel. This is given by
%
\begin{equation}\label{eq:VTVS}
\TV_S(\bm{u}) = \sum_{m=1}^{\mathcal{M}} \TV(\bm{u}_m).
\end{equation}
%
Since this penalty is additively separable, this would be akin to just reconstructing each channel separately, so we will refer to this approach as the "channel-by-channel" TV. 

\paragraph{A class of vectorial TV regularizers.}

Now, we will describe a more principled way to generalize the TV to multi-channel images. We mentioned earlier that the first-order derivatives of $\bm{u}$ can be characterized at every pixel-location by the Jacobian matrix $\left[D\bm{u} \right]_{\ell}$. An interesting class of regularizers is described by
%
\begin{equation}\label{eq:VTV}
\VTV(\bm{u}) = \sum_{\ell=1}^{\mathcal{L}} \|\left[D\bm{u}\right]_\ell\|_{S_p},
\end{equation}
%
where the matrix norm $\|.\|_{S_p}$ is known as the "Schatten p-norm," which is equivalent to the $\ell_p-$norm of the singular values, $\|X\|_{S_p} = \|\sigma(X)\|_{p}$, where $\sigma(X)$ is a vector comprised of the singular values of matrix $X$. Some special cases include $S_\infty$ (spectral norm), $S_2$ (Frobenius norm), and $S_1$ (nuclear norm). All three of these norms have been suggested as potential VTV penalties \cite{Blomgren1998,Goldluecke2012,Lefkimmiatis2013,Holt2014,Rigie2014}, and they share a few important properties. Because they are norms defined by the singular values of the Jacobian, they are rotationally invariant and convex. Additionally, for a single-channel image, they all reduce exactly to the conventional, isotropic TV. 

\paragraph{The Total Nuclear Variation (TNV)}
The total nuclear variation (TNV), defined as the nuclear norm of the Jacobian, is the only one of these VTV's that encourages the gradient vectors of all of the image channels to point in a common direction. We have found empirically that this is a very useful prior for spectral CT images, and that it generally outperforms the other VTV penalties. Therefore, in this work we will focus solely on the TNV, which is defined by 
%
\begin{equation}\label{eq:TNV}
\TNV(\bm{u}) = \sum_{\ell=1}^{\mathcal{L}} \|\left[D\bm{u}\right]_\ell\|_{S_1}.
\end{equation}
%
To see how this TNV couples the different energy channels of $\bm{u}$, note that the Jacobian at pixel $\ell$ can also be expressed in terms of the gradient vectors of the various image channels:
%
\begin{equation}\label{eq:Jacobian}
\left[D\bm{u}\right]_\ell = \begin{bmatrix}
\longleftarrow & \left[\nabla \uu_1\right]_\ell^T & \longrightarrow \\
\longleftarrow & \vdots & \longrightarrow \\
\longleftarrow & \left[\nabla \uu_{\mathcal{M}}\right]_\ell^T & \longrightarrow
\end{bmatrix}_{\mathcal{M}\times\mathcal{N}}.
\end{equation}
%
\begin{figure}[htbp]
\centering
\newcommand{\pwidth}{0.24\linewidth}
\subcaptionbox{All singular values are zero}{\includegraphics[width=\pwidth]{Tikz/sv_edges_a.pdf}} \hspace{0.5cm}
\subcaptionbox{1 non-zero singular value}{\includegraphics[width=\pwidth]{Tikz/sv_edges_b.pdf}} \hspace{0.5cm}
\subcaptionbox{All non-zero singular values}{\includegraphics[width=\pwidth]{Tikz/sv_edges_c.pdf}}
\caption{This illustrates the relationship between the directionality of the gradient vectors of a 3-channel image and the singular values of its Jacobian. When all image channels are constant valued, all singular values of the Jacobian are zero (a). If the gradient vectors are parallel or anti-parallel, one singular value is non-zero (b). All singular values will be nonzero if the gradients have unique directions (c).}
\label{fig:gradientdirection}
\end{figure}
%
\noindent The nuclear norm encourages sparsity in the singular values of $\left[D\bm{u}\right]_\ell$ \cite{Cai2010,Candes2011}. If pixel $\ell$ lies within a constant-valued region of all of the image channels then this Jacobian matrix will be entirely null-valued, and all of its singular values will be zero. Further, if all of the gradient vectors of the various image channels are parallel or anti-parallel, such as at a common edge, then the Jacobian will be rank one and thus have only one non-zero singular value. This is because parallel or anti-parallel vectors are not linearly independent. The fact that the TNV treats parallel and anti-parallel vectors the same means that it is robust to contrast inversions. We also point out that if one of the gradient vectors has zero magnitude, the rank is still one, so the TNV also appears to be robust in handling edges that may not exist in every image channel. In sum, encouraging singular-value sparsity in the Jacobian is akin to encouraging common edge locations and gradient directions. This notion is illustrated graphically in Figure \ref{fig:gradientdirection}.

\paragraph{Converting to a "noise-balanced" image space}

In spectral CT imaging it is often the case that certain color channels are significantly noisier than others, and we have found the TNV is much more effective if a noise-balancing transform is applied prior to reconstruction or denoising. For example, suppose we had a two-channel image, $\bm{u}$, with zero-mean, Gaussian noise with variances given by $\sigma_1^2$ and $\sigma_2^2$ for channels 1 and 2, respectively. If we wanted to use the TNV to denoise $\bm{u}$, we should first define a noise-balanced image $\bm{u'}$ as in equation \ref{eq:noisebalancing}.
%
\begin{equation}\label{eq:noisebalancing}
\bm{u'} = \begin{bmatrix} u_1/\sigma_1 \\ u_2/\sigma_2 \end{bmatrix}
\end{equation}
%
Each channel is globally scaled so that the per-channel noise is the same. After denoising, the inverse scale factors are applied in order to return to the original image space. For image reconstruction, the same type of scaling would be applied to the sinogram data. In general, the projection data will not have a uniform noise level within a single color channel, so we use some average measure of channel noise to determine the global scale factor. We find that this noise-balancing procedure is a very important step, as it greatly improves the noise suppression in multi-spectral images with unequal noise levels. In many cases, even images that are corrupted by extremely high noise levels may be well recovered using this technique as long as at least one of the color channels is reasonably clean. 

\begin{figure}[htbp]
\centering
\includegraphics[width=2in]{Figures/noise_balancing_noisy.png}
\includegraphics[width=2in]{Figures/noise_balancing_denoised.png}
\caption{A noisy two-channel image (left) with grossly mismatched noise levels. The image is successfully denoised using the TNV in conjunction with the noise-balancing procedure (right). Even the extremely noisy channel is recovered nearly perfectly due to the tight coupling to the clean image channel.}
\end{figure}

\section{Materials and methods}

\subsection{General reconstruction model}

Assume we have a set of $\mathcal{M}$ sinograms, corresponding either to the directly measured energy channels of a spectral CT system (color basis) or to the material basis projections. We jointly reconstruct these channels by solving the following data-constrained optimization problem:
%
\begin{equation}\label{eq:constrainedVTV}
\begin{split}
\min_{\bm{u}} &\quad \VTV(\bm{u})\\
\mathrm{s.t.} &\quad \|A\bm{u}-\bm{f}\|_W \leq \epsilon,
\end{split}
\end{equation}
%
where $\ff$ is the multi-channel projection data related to $\uu$ by the discrete fan-beam or cone-beam projection operator $A$. The projector is defined on a per channel basis,
%
\begin{equation}
\left(A\uu\right)_m = A_m\uu_m,
\end{equation}
%     
where $A_m$ characterizes the measurement geometry associated with channel $m$. The norm on the data constraint $\|.\|_W$ is a weighted $\ell_2$ norm, defined by 
%
\begin{equation}\label{eq:wl2}
\|x\|^2_W = \inner{x,Wx}.
\end{equation}
%
This penalized weighted least squares (PWLS) data model is often used in CT reconstruction because when the weights are selected such that $W = K^{-1}_{\bm{f}}$, the data fidelity term is a maximum likelihood estimator for Gaussian data with covariance $K_{\bm{f}}$. The single adjustable parameter $\epsilon$ controls the balance between data-fidelity and regularity. By fixing $\epsilon$ we can directly compare reconstructions using the na\"{i}ve channel-by-channel $\TV_S$ to our proposed $\TNV$ subject to the same data fidelity constraint. 

\subsection{The Chambolle Pock Primal Dual Algorithm}
In this section, we will provide an overview of the first-order, primal-dual algorithm of Chambolle and Pock \cite{Chambolle2011} and demonstrate how to apply it to our reconstruction model. For a tutorial on how to apply the CPPD algorithm to various CT reconstruction schemes, we refer the reader to Sidky et al \cite{Sidky2012}. 

\paragraph{The proximal mapping}
In order to describe the CPPD algorithm, we must first introduce the concept of a proximal mapping. For a function $f(x)$, the proximal mapping is defined by
%
\begin{equation}\label{eq:prox}
\prox{f}{x} = \argmin_u \,\, f(u) + \frac{1}{2}\|u-x\|_2^2.
\end{equation}
%
It will also be useful to consider the proximal mapping of the scaled function $\lambda f$, which can be written as
%
\begin{equation}\label{eq:proxscaled}
\prox{\lambda f}{x} = \argmin_u \,\, f(u) + \frac{1}{2\lambda}\|u-x\|_2^2.
\end{equation}
%
The proximal mapping can be interpreted as a generalized projection operator because for the special case that $f(x)$ is an indicator function, then it is precisely a euclidean projection. For an extensive overview of the prox operator and applications, refer to Parikh and Boyd \cite{Parikh2013}. 

\paragraph{A general saddle-point problem}

Let $X$ and $Y$ be two finite-dimensional, real vector spaces. The CPPD algorithm is  designed to solve the saddle-point problem described by
%
\begin{equation}\label{eq:saddle_xy}
\min_{x\in X}\,\max_{y\in Y} \,\, \inner{Kx,y} - F(y) + G(x),
\end{equation}
%
where $F(y)$ and $G(x)$ are convex functions with very simple proximal mappings, and $K$ is a general linear operator. In particular, the proximal mappings associated with $F$ and $G$ should have a closed form or be easily solvable to high precision. It turns out that many interesting image processing and reconstruction problems fit this description, including the VTV reconstruction model already described in the previous section. The CPPD update equations are summarized by Algorithm \ref{alg:CPPD}.
%
\begin{algorithm}
\caption{CPPD Algorithm}\label{alg:CPPD}
\begin{algorithmic}[1]
\State Initialize: $\theta\in \left[-1,1\right]$, $\,\, \sigma\tau\|K\|_2^2 < 1$, $\,\, x^{(0)},\bar{x}^{(0)} \in X, \, y^{(0)} \in Y$ \Comment{$\|K\|_2 \equiv \sigma_\mathrm{max}(K)$}
\Repeat
\State $y^{(k+1)} = \prox{\sigma F}{y^{(k)} + \sigma K \bar{x}^{(k)}}$
\State $x^{(k+1)} = \prox{\tau G}{x^{(k)} - \tau K^Ty^{(k+1)}}$
\State $\bar{x}^{(k+1)} = x^{(k+1)} + \theta\left(x^{(k+1)}-x^{(k)} \right)$
\Until $x^{(k+1)},y^{(k+1)} = x^{(k)},y^{(k)}$ \Comment{stop when convergence criteria met}
\end{algorithmic}
\end{algorithm}
%
The parameters $\theta,\sigma,\tau$ can be thought of as step-size parameters that affect convergence speed but not the final solution, and $\|K\|_2$ is the "spectral norm" of the operator $K$, which is equivalent to its largest singular value. The power method, outlined in Sidky et al \cite{Sidky2012}, provides a very efficient iterative procedure for computing $\|K\|_2$ for very large-scale problems, involving only matrix-vector multiplications. This algorithm has close ties with several other popular algorithms, including ADMM, split-bregman, and proximal-point. This is further discussed in Chambolle and Pock \cite{Chambolle2011}.

\paragraph{Applying the CPPD algorithm to VTV reconstruction}

Now we will outline how the CPPD algorithm can be applied to the data-constrained optimization problem of (\ref{eq:constrainedVTV}). First we rewrite (\ref{eq:constrainedVTV}) as 
%
\begin{equation}
\min_{\uu} \quad \VTV(\uu) + \delta_{\epsilon}\left(A\uu-\ff\right),
\end{equation}
%
where the indicator function $\delta_\epsilon(\bm{x})$ is defined by
%
\begin{equation}\label{eq:dataconstraint}
\delta_\epsilon(\xx) = \left\{\begin{array}{lr} 0 & \|\xx\|_W \leq \epsilon\\
\infty &  \|\xx\|_W > \epsilon
\end{array}\right..
\end{equation}
%
Now, we introduce an auxiliary variable $\qq$ and rewrite the data constraint as an optimization problem:
%
\begin{equation}\label{eq:dataconjugate}
\delta_\epsilon(A\uu-\ff) = \max_\qq \,\, \inner{\qq,A\uu-\ff} - \epsilon\|W^{-1/2}\qq\|_2.
\end{equation}
%
Similarly, we can rewrite the VTV penalty function as an optimization over another auxiliary variable, $\zz $. For the channel-by-channel TV, we have
%
\begin{equation}\label{eq:vtvsdual}
\TV_S(\uu) = \max_\zz \,\, \inner{D \uu,\zz} - \delta_\mathcal{S}(\zz),
\end{equation}
%
where the definition of the set $\mathcal{S}$ is given by
%
\begin{equation}\label{eq:setS}
\mathcal{S} = \left\{\zz \in \mathbb{R}^{(\mathcal{M}\times 2)\cdot\mathcal{L}} \,\, : \,\, \mathcal{}\|\left[\zz_m\right]_\ell\|_2 \leq 1 \,\,\, \forall\,\,\, m\in 1,\dotsc,\mathcal{M},\,\,;\ell\in 1,\dotsc,\mathcal{L}   \right\}.
\end{equation}
%
The quantity $\left[\zz_m \right]_\ell$ refers to image channel $m$ at pixel $\ell$. Similarly, we can write the proposed $\TNV$ penalty function in terms of an optimization over auxiliary variable $\zz$: 
%
\begin{equation}\label{eq:vtvndual}
\TNV(\uu) = \max_\zz \,\, \inner{D\uu,\zz} - \delta_{N}(\zz),
\end{equation}
%
where the set $N$ is defined by 
%
\begin{equation}\label{eq:setN}
N = \left\{\zz \in \mathbb{R}^{(\mathcal{M}\times 2)\cdot\mathcal{L}} \,\, : \,\, \sigma_\mathrm{max}(\left[\zz\right]_\ell) \leq 1 \,\,\, \forall\,\,\, \ell\in 1,\dotsc,\mathcal{L}  \right\}.
\end{equation}
%
The quantity $[\zz]_\ell$ is an $\mathcal{M}\times2$ matrix, with the same dimensions as the Jacobian derivative at pixel $\ell$, and $\sigma_\mathrm{max}$ is its largest singular value. These transformations, which are detailed in Appendix A, allow us to write our original VTV reconstruction model as a primal-dual, saddle-point problem,
%
\begin{equation}\label{eq:saddle_qz}
\min_\uu\,\max_{\zz,\qq} \,\, \inner{D\uu,\zz} + \inner{\qq,A\uu-\ff} - \delta_{\mathcal{S}/N}(\zz) - \epsilon\|W^{-1/2}\qq\|_2,
\end{equation}
%
where the indicator function is either $\delta_{\mathcal{S}}$ for the channel-by-channel $\TV_S$ or $\delta_N$ for the proposed $\TNV$. We can directly apply the CPPD update equations of Algorithm \ref{alg:CPPD}  by making the following assignments.
%
\begin{equation}
\begin{split}
x &\to \uu,\quad
y \to \begin{pmatrix} \qq \\ \zz \end{pmatrix}, \quad
K \to \begin{pmatrix} A \\ D \end{pmatrix}\\
F(\qq,\zz) &= \delta_{\mathcal{S}/N}(\zz) - \epsilon\|W^{-1/2}\qq\|_2-\inner{\qq,\ff},\quad
G(\uu) = 0\\
\end{split}
\end{equation}
%
The resulting update equations are given by Algorithm \ref{alg:VTVupdates}.
\begin{algorithm}
\caption{Data-constrained VTV Update Equations}\label{alg:VTVupdates}
\begin{algorithmic}[1]
\State Initialize: $\theta=1$, $\,\, \sigma\tau\left\|\begin{pmatrix}A\\D\end{pmatrix}\right\|_2^2 = 1$, $\,\,(\uu,\uubar, \qq,\zz)\to \bm{0}$ 
\Repeat
\State $\zz^{(k+1)} = \Pi_{\mathcal{S}/N} \left(\zz^{(k)} + \sigma D \uubar^{(k)} \right) $
\State $\qq^{(k+1)} = \prox{\sigma\epsilon F_q}{\qq^{(k)} + \sigma \left(A\uubar^{(k)}-\ff\right)}$
\State $\uu^{(k+1)} = \uu^{(k)} + \tau\left(\Div\zz^{(k+1)} - A^T\qq^{(k+1)} \right)$
\State $\uubar^{(k+1)} = \uu^{(k+1)} + \theta\left(\uu^{(k+1)}-\uu^{(k)} \right)$
\Until  \Comment{stop when convergence criteria met}
\end{algorithmic}
\end{algorithm}
%
The operator $\Pi_{\mathcal{S}/N}$ is a Euclidean projection onto the set $\mathcal{S}$ or $N$, and the $\Div$ operator is the negative transpose of the discrete derivative operator $D$. Note that the $\qq$-update is written in terms of the proximal map of the function $F_q$, which we define as $F_q = \|W^{(-1/2)}\qq\|_2$. While there does not appear to be a closed-form expression for this proximal map, it can be reduced to a very simple scalar optimization problem. This is detailed in the appendix along with an efficient method for performing the Euclidean projection in the primal variable update. 

\subsection{Simulation studies}

To investigate the impact of the proposed vectorial TV regularization, we performed two numerical simulation studies with the computerized XCAT phantom. In the first experiment we simulate an ideal, photon-counting system with 5 energy windows and directly reconstruct images corresponding to the log-normalized bin data. We will refer to this image basis as the "color" basis because each image channel corresponds directly to the measurements of one energy window. 

In the second experiment, we simulate an ideal, dual-kVp scan and perform what we will refer to as a "hybrid" reconstruction. As we will detail later, we first decompose the 80/140 kVp data into a bone/soft-tissue "material" basis and then we co-reconstruct this synthetic data with the raw 80/140 kVp sinograms. The $\TNV$ penalty couples all four image channels so that the relatively noisy basis-material channels may benefit from the higher SNR color channels. 

\subsubsection{Data generation}

\paragraph{The XCAT shoulder phantom}

All simulations used the same 2D pixelized shoulder phantom, which was generated from an axial slice of the XCAT phantom and is pictured in Figure \ref{fig:XCAT}. The XCAT software package was used to generated a set of bone and soft-tissue density maps on a $2048\times2048$ pixel grid. These density maps were then used as input to a polyenergetic, distance-driven projector model to generate the simulated raw data that would be inputted to our reconstruction algorithm. 

\paragraph{Ideal photon counting model}

To generate the 5-bin photon-counting data, we used a realistic 120 kVp simulated x-ray tube spectrum with hard thresholds at 40,60,80,100, and 120 keV. We did not model any physical factors in the detector response, so our bin response functions are perfect rect functions. Our simulated spectrum had virtually no emissions below 20 keV, so we can think of these bins as being evenly spaced. Using a distance-driven projector model, polyenergetic projections were computed with 896 detector elements per view and 400 views, and Poisson noise was added. 

\paragraph{Dual kVp hybrid model}
The other configuration we looked at was an 80/140 dual kVp acquisition with the same number of detectors and views as the photon counting simulation. We simulated realistic 80 and 140 kVp x-ray tube spectra to generate two sets of consistent projection data. This type of consistent, dual-kVp data can be acquired on many current scanners by performing two back-to-back scans or can be estimated from a fast kV-switching geometry. Gaussian noise was added to approximate a compound Poisson noise model. From the noisy 80 and 140 kVp projection data, we synthesized bone and soft-tissue basis sinograms using a maximum-likelihood material decomposition. In the "hybrid" study, we will co-reconstruct these 4 channels of data, consisting of our log-normalized, dual kVp data (color basis) and the synthesized, material-basis data. 

\begin{figure}[htbp]
\centering
\newcommand{\pwidth}{0.31\linewidth}
\subcaptionbox{60 keV reference image}{
	\includegraphics[width=\pwidth,trim=0 350 0 500,clip=true]{Figures/xcat_60kev.png}
} \hfill
\subcaptionbox{bone density map}{
	\includegraphics[width=\pwidth,trim=0 350 0 500,clip=true]{Figures/xcat_bone.png}
} \hfill
\subcaptionbox{soft-tissue density map}{
	\includegraphics[width=\pwidth,trim=0 350 0 500,clip=true]{Figures/xcat_tissue.png}
} \hfill
\caption{The 2D XCAT shoulder phantom used for all simulations studies, depicted at 60 keV (a) for reference. The projection data were generated using bone (b) and soft-tissue (c) density maps with the appropriate mass-attenuation curves.}
\label{fig:XCAT}
\end{figure}

\subsubsection{Preprocessing and reconstruction}

Both the 5-channel photon-counting data and the 4-channel dual kVp hybrid data were reconstructed using the reconstruction model outlined in (\ref{eq:constrainedVTV}). For reconstruction, we used a $512\times512$ grid of $1$ mm pixels and a matched projector/backprojector pair based on Joseph's method. 

\paragraph{Computing data-weights}

For the raw-basis projection data, we will use the data-weighting approach described by Bouman \cite{Bouman1996}, which results in a diagonal $W$ matrix, where the diagonal elements equal the pre-logged projection data. This is a quadratic approximation to the log-likelihood. For the synthetic, bone/soft-tissue sinograms, we estimate the variances using the Cramer-Rao lower bound \cite{Roessl2009} and weight by their inverse. This is similar to the approach described by Schirra \cite{Schirra2013} and Sawatzky \cite{Sawatzky2014}, but we ignore the off-diagonal terms in this work for simplicity. We found that having a non-diagonal weighting matrix made it immensely more complicated to derive the algorithm update equations when using this data-constrained form. However, in practice, it is often easier to solve the regularized (unconstrained) problem, and incorporating the off-diagonal covariance terms poses no challenge in this case. We focus on the data-constrained form in this work, simply because it provides a straightforward way to compare the vectorial TV to the channel-by-channel TV.

\paragraph{Tuning the data-constraint parameter}
The only parameter that we will vary in our reconstruction model (\ref{eq:constrainedVTV}) is $\epsilon$, which controls the trade-off between data fidelity and smoothness. In general, smaller values of $\epsilon$ will result in noisier images, while higher values allow the $\VTV$ term to find a smoother solution. At the extreme case of $\epsilon=0$, only images with projections that exactly match the measured data are allowed. There is also some finite value $\epsilon=\epsilon_\mathrm{max}$ for which the algorithm will return an image of all zeros. In order to determine an interesting range for $\epsilon$, we first compute a ground truth image $\utrue$ by performing FBP on noiseless, non-sparse projection data. Then we define a new parameter $\epsilon^*$ which is defined by equation \ref{eq:epsilonstar},
%
\begin{equation}\label{eq:epsilonstar}
\epsilon^* \equiv \|A\utrue-\ff\|_W
\end{equation}
%
where $\ff$ are the noisy projection data we wish to reconstruct. In this study we will select values of $\epsilon$ that satisfy $\epsilon = \alpha\epsilon^*$, where $\alpha \in \left(0,1\right)$. We subjectively selected a range of $\alpha$ values that represent a range of solutions from under-smoothed to over-smoothed in order to demonstrate how the TNV compares to the channel-by-channel TV in various regimes. As described earlier, the projection data were re-scaled prior to reconstruction, so that the average noise levels were approximately the same in every spectral channel.

\begin{comment}
\paragraph{Converting to a noise-balanced coordinate space}
As described earlier, we scale the data so that the average noise-level is approximately equal among all image channels. 
%
\begin{align}
\min_{\bm{u}} &\quad \VTV(\bm{u})\\
\mathrm{s.t.} &\quad \|A\bm{u}-\bm{f}\|_W \leq \epsilon
\end{align}
%
When the weighting matrix $W$ is diagonal, we can write this as
%
\begin{align}
\min_{\bm{u}} &\quad \VTV(\bm{u})\\
\mathrm{s.t.} &\quad \sum_j\left\|W_j^{1/2}(A_j\bm{u}_j-\bm{f_j})\right\|^2_2 \leq \epsilon^2.
\end{align}
%
where $j$ indexes the spectral channel, and $W^{1/2}_j$ is a diagonal matrix where the diagonal elements correspond to the approximate the standard deviations of the noisy projection data $f_j$. By defining the average noise-level for each channel by the scalar, $\bar{W}^{1/2}_j$, as the mean value of $\mathrm{diag}\left(W^{1/2}_j\right)$, we can separate out the ray-by-ray noise variations from the global channel-by-channel noise variations. We rewrite the optimization problem as
%
\begin{align}
\min_{\bm{u}} &\quad \VTV(\bm{u})\\
\mathrm{s.t.} &\quad \sum_j\left\|\hat{W}_j^{1/2}\bar{W}^{1/2}_j(A_j\bm{u}_j-\bm{f_j})\right\|^2_2 \leq \epsilon^2,
\end{align}
%
where $\hat{W}_j^{1/2} \equiv W^{1/2}_j/\bar{W}^{1/2}_j$. We observe that by making the following change of variables
%
\begin{align}
\bm{u}'_j &\gets \bar{W}_j^{1/2}\bm{u}_j\\
\bm{f}'_j &\gets \bar{W}_j^{1/2}\bm{f}_j\\
\end{align}
%
and solving the new optimization problem over $\bm{u'}$, given by (\ref{eq:weightedproblem}),
%
\begin{equation}\label{eq:weightedproblem}
\begin{split}
\min_{\bm{u'}} &\quad \VTV(\bm{u'})\\
\mathrm{s.t.} &\quad \sum_j\left\|\hat{W}_j^{1/2}(A_j\bm{u'}_j-\bm{f'_j})\right\|^2_2 \leq \epsilon^2.
\end{split}
\end{equation}
%
we achieve better overall noise suppression. We simply convert back to the original image space after the algorithm is converged by undoing the scalar, per-channel weighting. This change of basis has no impact on the channel-by-channel TV reconstruction, but it does non-trivially change the behavior of the $\TNV$. We do not currently have a theoretical understanding of why this data-scaling is advantageous. 

\end{comment}

\section{Results}

\begin{figure}[htbp]
\centering
\includegraphics[width=2.5in]{Tikz/xcat_rois.pdf}
\caption{XCAT reference image, with ROI's indicated in yellow}
\label{fig:rois}
\end{figure}

\subsection{Photon counting study}

Here we present the results of performing a simultaneous reconstruction of the five photon-counting bins after log-normalization, using both the channel-by-channel $\TV_S$ and the proposed $\TNV$. This is an example of performing a joint reconstruction on data in the color basis because we did not perform a material decomposition. In this setup, bin 1 (0-40 keV) was the noisiest, due to significant attenuation below 40 keV. We find that the inter-channel coupling introduced by the total nuclear variation has the greatest impact on the noisiest channels, so we present reconstructed images from this energy bin, focusing on the ROI's indicated in Figure \ref{fig:rois}. The resulting images are depicted in figures \ref{fig:bin1bone} through \ref{fig:bin1tissue}. The image window was manually adjusted for each ROI to highlight the relevant structures but is fixed within a particular figure. In general, we find that the images reconstructed with $\TNV$ regularization suffer from less edge blurring as the smoothing parameter $\epsilon$ is increased. The profile plot in Figure \ref{fig:profile_plot} gives a closer look at how the TNV better preserves bony structures in the lowest energy bin image. We confirmed that these profiles were extracted from images of similar noise levels by measuring the sample variance in a nearby, uniform muscle ROI. This is expected because our data-constrained, reconstruction model allows us to select the noise level directly by tuning the $\epsilon$ parameter. 

\begin{figure}[htbp]
\centering
\newcommand{\pwidtha}{1in}
% 0013,0016,001,002,003
\begin{minipage}[t]{0.85\linewidth}
\subcaptionbox{${\TV_S}$\\${\epsilon=0.001\epsilon^*}$\\RMSE$ = 0.0444$}{\includegraphics[width=\pwidtha]{Figures/Results/PC/u_0-001_channel0_TVS_bone_ROI.png}}\hfill
\subcaptionbox{$\TV_S$,\\${\epsilon=0.0013\epsilon^*}$\\RMSE$ = 0.0436$}{\includegraphics[width=\pwidtha]{Figures/Results/PC/u_0-0013_channel0_TVS_bone_ROI.png}}\hfill
\subcaptionbox{$\TV_S$,\\${\epsilon=0.0016\epsilon^*}$\\RMSE$ = 0.0450$}{\includegraphics[width=\pwidtha]{Figures/Results/PC/u_0-0016_channel0_TVS_bone_ROI.png}}\hfill
\subcaptionbox{$\TV_S$,\\${\epsilon=0.002\epsilon^*}$\\RMSE$ = 0.0484$}{\includegraphics[width=\pwidtha]{Figures/Results/PC/u_0-002_channel0_TVS_bone_ROI.png}}\hfill
\subcaptionbox{$\TV_S$,\\${\epsilon=0.003\epsilon^*}$\\RMSE$ = 0.0609$}{\includegraphics[width=\pwidtha]{Figures/Results/PC/u_0-003_channel0_TVS_bone_ROI.png}}\hfill

\subcaptionbox{$\TNV$,\\${\epsilon=0.001\epsilon^*}$\\RMSE$ = 0.0390$}{\includegraphics[width=\pwidtha]{Figures/Results/PC/u_0-001_channel0_TVN_bone_ROI.png}}\hfill
\subcaptionbox{$\TNV$,\\${\epsilon=0.0013\epsilon^*}$\\RMSE$ = 0.0347$}{\includegraphics[width=\pwidtha]{Figures/Results/PC/u_0-0013_channel0_TVN_bone_ROI.png}}\hfill
\subcaptionbox{$\TNV$,\\${\epsilon=0.0016\epsilon^*}$\\RMSE$ = 0.0332$}{\includegraphics[width=\pwidtha]{Figures/Results/PC/u_0-0016_channel0_TVN_bone_ROI.png}}\hfill
\subcaptionbox{$\TNV$,\\${\epsilon=0.002\epsilon^*}$\\RMSE$ = 0.0327$}{\includegraphics[width=\pwidtha]{Figures/Results/PC/u_0-002_channel0_TVN_bone_ROI.png}}\hfill
\subcaptionbox{$\TNV$,\\${\epsilon=0.003\epsilon^*}$\\RMSE$ = 0.0456$}{\includegraphics[width=\pwidtha]{Figures/Results/PC/u_0-003_channel0_TVN_bone_ROI.png}}\hfill
\end{minipage}
\caption{Bin 1 (0-40keV) image, bone ROI/window comparison with channel-by-channel TV (top) and $\TNV$ (bottom). Grayscale window in cm$^{-1}$ $\left[0.30,0.85 \right]$.}
\label{fig:bin1bone}
\end{figure}
\begin{figure}[htbp]
\centering
\newcommand{\pwidtha}{1in}
\begin{minipage}[t]{0.85\linewidth}
\subcaptionbox{$\TV_S$,\\${\epsilon=0.001\epsilon^*}$\\RMSE$ = 0.0219$}{\includegraphics[width=\pwidtha]{Figures/Results/PC/u_0-001_channel0_TVS_soft_tissue_ROI.png}}\hfill
\subcaptionbox{$\TV_S$,\\${\epsilon=0.0013\epsilon^*}$\\RMSE$ = 0.0200$}{\includegraphics[width=\pwidtha]{Figures/Results/PC/u_0-0013_channel0_TVS_soft_tissue_ROI.png}}\hfill
\subcaptionbox{$\TV_S$,\\${\epsilon=0.0016\epsilon^*}$\\RMSE$ = 0.0210$}{\includegraphics[width=\pwidtha]{Figures/Results/PC/u_0-0016_channel0_TVS_soft_tissue_ROI.png}}\hfill
\subcaptionbox{$\TV_S$,\\${\epsilon=0.002\epsilon^*}$\\RMSE$ = 0.0223$}{\includegraphics[width=\pwidtha]{Figures/Results/PC/u_0-002_channel0_TVS_soft_tissue_ROI.png}}\hfill
\subcaptionbox{$\TV_S$,\\${\epsilon=0.003\epsilon^*}$\\RMSE$ = 0.0266$}{\includegraphics[width=\pwidtha]{Figures/Results/PC/u_0-003_channel0_TVS_soft_tissue_ROI.png}}\hfill

\subcaptionbox{$\TNV$,\\${\epsilon=0.001\epsilon^*}$\\RMSE$ = 0.0216$}{\includegraphics[width=\pwidtha]{Figures/Results/PC/u_0-001_channel0_TVN_soft_tissue_ROI.png}}\hfill
\subcaptionbox{$\TNV$,\\${\epsilon=0.0013\epsilon^*}$\\RMSE$ = 0.0176$}{\includegraphics[width=\pwidtha]{Figures/Results/PC/u_0-0013_channel0_TVN_soft_tissue_ROI.png}}\hfill
\subcaptionbox{$\TNV$,\\${\epsilon=0.0016\epsilon^*}$\\RMSE$ = 0.0181$}{\includegraphics[width=\pwidtha]{Figures/Results/PC/u_0-0016_channel0_TVN_soft_tissue_ROI.png}}\hfill
\subcaptionbox{$\TNV$,\\${\epsilon=0.002\epsilon^*}$\\RMSE$ = 0.0189$}{\includegraphics[width=\pwidtha]{Figures/Results/PC/u_0-002_channel0_TVN_soft_tissue_ROI.png}}\hfill
\subcaptionbox{$\TNV$,\\${\epsilon=0.003\epsilon^*}$\\RMSE$ = 0.0175$}{\includegraphics[width=\pwidtha]{Figures/Results/PC/u_0-003_channel0_TVN_soft_tissue_ROI.png}}\hfill
\end{minipage}
\caption{Bin 1 (0-40keV) image,  soft-tissue ROI comparison with channel-by-channel TV (top) and $\TNV$ (bottom). Grayscale window in cm$^{-1}$ $\left[0.30,0.35 \right]$.}
\label{fig:bin1tissue}
\end{figure}

\begin{comment}
\begin{figure}[htbp]
\centering
\newcommand{\pwidtha}{1in}
\begin{minipage}[t]{0.85\linewidth}
\subcaptionbox{$\TV_S$,\\${\epsilon=0.001\epsilon^*}$\\RMSE$ = 0.00713$}{\includegraphics[width=\pwidtha]{Figures/Results/PC/u_0-001_channel4_TVS_bone_ROI.png}}\hfill
\subcaptionbox{$\TV_S$,\\${\epsilon=0.0013\epsilon^*}$\\RMSE$ = 0.00499$}{\includegraphics[width=\pwidtha]{Figures/Results/PC/u_0-0013_channel4_TVS_bone_ROI.png}}\hfill
\subcaptionbox{$\TV_S$,\\${\epsilon=0.0016\epsilon^*}$\\RMSE$ = 0.00424$}{\includegraphics[width=\pwidtha]{Figures/Results/PC/u_0-0016_channel4_TVS_bone_ROI.png}}\hfill
\subcaptionbox{$\TV_S$,\\${\epsilon=0.002\epsilon^*}$\\RMSE$ = 0.00388$}{\includegraphics[width=\pwidtha]{Figures/Results/PC/u_0-002_channel4_TVS_bone_ROI.png}}\hfill
\subcaptionbox{$\TV_S$,\\${\epsilon=0.003\epsilon^*}$\\RMSE$ = 0.00640$}{\includegraphics[width=\pwidtha]{Figures/Results/PC/u_0-003_channel4_TVS_bone_ROI.png}}\hfill

\subcaptionbox{$\TNV$,\\${\epsilon=0.001\epsilon^*}$\\RMSE$ = 0.00760$}{\includegraphics[width=\pwidtha]{Figures/Results/PC/u_0-001_channel4_TVN_bone_ROI.png}}\hfill
\subcaptionbox{$\TNV$,\\${\epsilon=0.0013\epsilon^*}$\\RMSE$ = 0.00515$}{\includegraphics[width=\pwidtha]{Figures/Results/PC/u_0-0013_channel4_TVN_bone_ROI.png}}\hfill
\subcaptionbox{$\TNV$,\\${\epsilon=0.0016\epsilon^*}$\\RMSE$ = 0.00435$}{\includegraphics[width=\pwidtha]{Figures/Results/PC/u_0-0016_channel4_TVN_bone_ROI.png}}\hfill
\subcaptionbox{$\TNV$,\\${\epsilon=0.002\epsilon^*}$\\RMSE$ = 0.00405$}{\includegraphics[width=\pwidtha]{Figures/Results/PC/u_0-002_channel4_TVN_bone_ROI.png}}\hfill
\subcaptionbox{$\TNV$,\\${\epsilon=0.003\epsilon^*}$\\RMSE$ = 0.00562$}{\includegraphics[width=\pwidtha]{Figures/Results/PC/u_0-003_channel4_TVN_bone_ROI.png}}\hfill
\end{minipage}
\caption{Bin 5 (100-120keV) image, bone ROI/window comparison with channel-by-channel TV (top) and $\TNV$ (bottom). Grayscale window in cm$^{-1}$ $\left[0.18,0.24 \right]$.}
\label{fig:bin5bone}
\end{figure}
%
\begin{figure}[htbp]
\centering
\newcommand{\pwidtha}{1in}
\begin{minipage}[t]{0.85\linewidth}
\subcaptionbox{$\TV_S$,\\${\epsilon=0.001\epsilon^*}$\\RMSE$ = 0.00876$}{\includegraphics[width=\pwidtha]{Figures/Results/PC/u_0-001_channel4_TVS_soft_tissue_ROI.png}}\hfill
\subcaptionbox{$\TV_S$,\\${\epsilon=0.0013\epsilon^*}$\\RMSE$ = 0.00726$}{\includegraphics[width=\pwidtha]{Figures/Results/PC/u_0-0013_channel4_TVS_soft_tissue_ROI.png}}\hfill
\subcaptionbox{$\TV_S$,\\${\epsilon=0.0016\epsilon^*}$\\RMSE$ = 0.00700$}{\includegraphics[width=\pwidtha]{Figures/Results/PC/u_0-0016_channel4_TVS_soft_tissue_ROI.png}}\hfill
\subcaptionbox{$\TV_S$,\\${\epsilon=0.002\epsilon^*}$\\RMSE$ = 0.00637$}{\includegraphics[width=\pwidtha]{Figures/Results/PC/u_0-002_channel4_TVS_soft_tissue_ROI.png}}\hfill
\subcaptionbox{$\TV_S$,\\${\epsilon=0.003\epsilon^*}$\\RMSE$ = 0.00493$}{\includegraphics[width=\pwidtha]{Figures/Results/PC/u_0-003_channel4_TVS_soft_tissue_ROI.png}}\hfill

\subcaptionbox{$\TNV$,\\${\epsilon=0.001\epsilon^*}$\\RMSE$ = 0.00939$}{\includegraphics[width=\pwidtha]{Figures/Results/PC/u_0-001_channel4_TVN_soft_tissue_ROI.png}}\hfill
\subcaptionbox{$\TNV$,\\${\epsilon=0.0013\epsilon^*}$\\RMSE$ = 0.00705$}{\includegraphics[width=\pwidtha]{Figures/Results/PC/u_0-0013_channel4_TVN_soft_tissue_ROI.png}}\hfill
\subcaptionbox{$\TNV$,\\${\epsilon=0.0016\epsilon^*}$\\RMSE$ = 0.00667$}{\includegraphics[width=\pwidtha]{Figures/Results/PC/u_0-0016_channel4_TVN_soft_tissue_ROI.png}}\hfill
\subcaptionbox{$\TNV$,\\${\epsilon=0.002\epsilon^*}$\\RMSE$ = 0.00656$}{\includegraphics[width=\pwidtha]{Figures/Results/PC/u_0-002_channel4_TVN_soft_tissue_ROI.png}}\hfill
\subcaptionbox{$\TNV$,\\${\epsilon=0.003\epsilon^*}$\\RMSE$ = 0.00538$}{\includegraphics[width=\pwidtha]{Figures/Results/PC/u_0-003_channel4_TVN_soft_tissue_ROI.png}}\hfill
\end{minipage}
\caption{Bin 5 (100-120keV) image,  soft-tissue ROI/window comparison with channel-by-channel TV (top) and $\TNV$ (bottom). Grayscale window in cm$^{-1}$ $\left[0.17,0.19 \right]$.}
\label{fig:bin5tissue}
\end{figure}
\end{comment}

\begin{figure}[htbp]
\centering
\begin{tikzpicture}
\node at (0,0) {\includegraphics[width=3in]{Figures/bin1_bone_profile_0016.pdf}};
\node at (0,-1.3in) {Pixel Index};
\begin{scope}
\node[rotate=90] at (-1.7in,0in) {Attenuation (cm$^{-1}$)};
\end{scope}
%\node at (5,0) {\includegraphics[width=1in]{Figures/bone_profile_indicator.png}};
\end{tikzpicture}
\caption{This is a vertical profile through the 0-40 keV, bone-ROI image with $\epsilon = 0.0016 \epsilon^*$. The TNV regularized reconstruction shows slightly better preservation of bony structures. The noise levels were estimated from a nearby ROI in a uniform muscle region ($\sigma_{TV_S} = 0.00245$, $\sigma_{\TNV} = 0.00221)$.}
\label{fig:profile_plot}
\end{figure}

\FloatBarrier

\subsection{Dual kVp Hybrid Study}

In this study, we co-reconstruct the synthetic soft-tissue/bone image channels with the raw 80/140 kVp data. In general, the material decomposition greatly amplifies noise, so that the soft-tissue and bone images tend to have a much lower CNR than the raw 80 and 140 kVp images. This noise amplification is due to the ill-conditioning of the inversion step in the basis change, which is caused by the relatively poor spectral separation between basis materials. This poor spectral separation also explains the strong negative correlation between the synthesized material channels. We hypothesize that by coupling the high CNR raw image channels and the low CNR synthetic image channels, we may be able to improve noise suppression in the synthetic data. We call this technique of reconstructing the synthetic and raw data simultaneously "hybrid" reconstruction.  We present ROI's from both the bone density image and the soft-tissue density image over a range of different $\epsilon$ values in Figures \ref{fig:boneimage} and \ref{fig:tissueimage}. 
%
\begin{figure}[htbp]
\centering
\newcommand{\pwidtha}{1in}
\begin{minipage}[t]{0.85\linewidth}
\subcaptionbox{$\TV_S$,\\${\epsilon=0.06\epsilon^*}$\\RMSE$ = 0.0417$}{\includegraphics[width=\pwidtha]{Figures/Results/Hybrid/u_0-06_channel2_TVS_bone_ROI.png}}\hfill
\subcaptionbox{$\TV_S$,\\${\epsilon=0.07\epsilon^*}$\\RMSE$ = 0.0421$}{\includegraphics[width=\pwidtha]{Figures/Results/Hybrid/u_0-07_channel2_TVS_bone_ROI.png}}\hfill
\subcaptionbox{$\TV_S$,\\${\epsilon=0.08\epsilon^*}$\\RMSE$ = 0.0447$}{\includegraphics[width=\pwidtha]{Figures/Results/Hybrid/u_0-08_channel2_TVS_bone_ROI.png}}\hfill
\subcaptionbox{$\TV_S$,\\${\epsilon=0.09\epsilon^*}$\\RMSE$ = 0.0505$}{\includegraphics[width=\pwidtha]{Figures/Results/Hybrid/u_0-09_channel2_TVS_bone_ROI.png}}\hfill
\subcaptionbox{$\TV_S$,\\${\epsilon=0.1\epsilon^*}$\\RMSE$ = 0.0564$}{\includegraphics[width=\pwidtha]{Figures/Results/Hybrid/u_0-1_channel2_TVS_bone_ROI.png}}\hfill

\subcaptionbox{$\TNV$,\\${\epsilon=0.06\epsilon^*}$\\RMSE$ = 0.0410$}{\includegraphics[width=\pwidtha]{Figures/Results/Hybrid/u_0-06_channel2_TVN_bone_ROI.png}}\hfill
\subcaptionbox{$\TNV$,\\${\epsilon=0.07\epsilon^*}$\\RMSE$ = 0.0398$}{\includegraphics[width=\pwidtha]{Figures/Results/Hybrid/u_0-07_channel2_TVN_bone_ROI.png}}\hfill
\subcaptionbox{$\TNV$,\\${\epsilon=0.08\epsilon^*}$\\RMSE$ = 0.0392$}{\includegraphics[width=\pwidtha]{Figures/Results/Hybrid/u_0-08_channel2_TVN_bone_ROI.png}}\hfill
\subcaptionbox{$\TNV$,\\${\epsilon=0.09\epsilon^*}$\\RMSE$ = 0.0392$}{\includegraphics[width=\pwidtha]{Figures/Results/Hybrid/u_0-09_channel2_TVN_bone_ROI.png}}\hfill
\subcaptionbox{$\TNV$,\\${\epsilon=0.1\epsilon^*}$\\RMSE$ = 0.0405$}{\includegraphics[width=\pwidtha]{Figures/Results/Hybrid/u_0-1_channel2_TVN_bone_ROI.png}}\hfill
\end{minipage}
\caption{Bone basis image, bone ROI comparison with channel-by-channel TV (top) and $\TNV$ (bottom). Grayscale window in g/mL$^{-1}$ $\left[0.00,0.52 \right]$.}
\label{fig:boneimage}
\end{figure}
\begin{figure}[htbp]
\centering
\newcommand{\pwidtha}{1in}
\begin{minipage}[t]{0.85\linewidth}
\subcaptionbox{$\TV_S$,\\${\epsilon=0.06\epsilon^*}$\\RMSE$ = 0.0445$}{\includegraphics[width=\pwidtha]{Figures/Results/Hybrid/u_0-06_channel3_TVS_soft_tissue_ROI.png}}\hfill
\subcaptionbox{$\TV_S$,\\${\epsilon=0.07\epsilon^*}$\\RMSE$ = 0.0410$}{\includegraphics[width=\pwidtha]{Figures/Results/Hybrid/u_0-07_channel3_TVS_soft_tissue_ROI.png}}\hfill
\subcaptionbox{$\TV_S$,\\${\epsilon=0.08\epsilon^*}$\\RMSE$ = 0.0383$}{\includegraphics[width=\pwidtha]{Figures/Results/Hybrid/u_0-08_channel3_TVS_soft_tissue_ROI.png}}\hfill
\subcaptionbox{$\TV_S$,\\${\epsilon=0.09\epsilon^*}$\\RMSE$ = 0.0376$}{\includegraphics[width=\pwidtha]{Figures/Results/Hybrid/u_0-09_channel3_TVS_soft_tissue_ROI.png}}\hfill
\subcaptionbox{$\TV_S$,\\${\epsilon=0.1\epsilon^*}$\\RMSE$ = 0.0400$}{\includegraphics[width=\pwidtha]{Figures/Results/Hybrid/u_0-1_channel3_TVS_soft_tissue_ROI.png}}\hfill

\subcaptionbox{$\TNV$,\\${\epsilon=0.06\epsilon^*}$\\RMSE$ = 0.0431$}{\includegraphics[width=\pwidtha]{Figures/Results/Hybrid/u_0-06_channel3_TVN_soft_tissue_ROI.png}}\hfill
\subcaptionbox{$\TNV$,\\${\epsilon=0.07\epsilon^*}$\\RMSE$ = 0.0396$}{\includegraphics[width=\pwidtha]{Figures/Results/Hybrid/u_0-07_channel3_TVN_soft_tissue_ROI.png}}\hfill
\subcaptionbox{$\TNV$,\\${\epsilon=0.08\epsilon^*}$\\RMSE$ = 0.0391$}{\includegraphics[width=\pwidtha]{Figures/Results/Hybrid/u_0-08_channel3_TVN_soft_tissue_ROI.png}}\hfill
\subcaptionbox{$\TNV$,\\${\epsilon=0.09\epsilon^*}$\\RMSE$ = 0.0386$}{\includegraphics[width=\pwidtha]{Figures/Results/Hybrid/u_0-09_channel3_TVN_soft_tissue_ROI.png}}\hfill
\subcaptionbox{$\TNV$,\\${\epsilon=0.1\epsilon^*}$\\RMSE$ = 0.0372$}{\includegraphics[width=\pwidtha]{Figures/Results/Hybrid/u_0-1_channel3_TVN_soft_tissue_ROI.png}}\hfill
\end{minipage}
\caption{soft-tissue basis image, soft-tissue ROI comparison with channel-by-channel TV (top) and $\TNV$ (bottom). Grayscale window in g/mL$^{-1}$ $\left[0.80,1.07 \right]$}
\label{fig:tissueimage}
\end{figure}
%
We find that using the TNV allows for a high degree of noise suppression (high values of $\epsilon$) without significantly deteriorating bone or soft-tissue structures. However, when the na\"{i}ve channel-by-channel TV is used, these same $\epsilon$ values eventually lead to significant blurring artifacts. We also point out that even though the soft-tissue, 80 kVp, and 140 kVp images have edges that are not present in the bone density image, the $\TNV$ coupling does not falsely propagate these edges into the bone channel. 

\FloatBarrier

\section{Summary}

We have described a framework for jointly reconstructing multi-channel spectral CT images using a generalized vectorial regularizer. Specifically, we presented a multi-channel generalization of the total variation which couples the different image channels by encouraging their gradient vectors to point in a parallel or anti-parallel direction. Preliminary results suggest that this coupling may allow for greater noise suppression with less blurring of image structures. This regularization strategy can be used to reconstruct raw energy bin data, synthetic basis-material images, or even both simultaneously, such as in our hybrid reconstruction study. This hybrid approach may allow for better noise suppression in the synthetic material images, which typically suffer from elevated noise levels. We believe that the proposed $\TNV$ provides an elegant way to regularize multi-channel images and may be a useful complement to previously published statistical, spectral CT reconstruction methods. 

\bibliography{VTV2}

% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 
% % % % % % % % % % % % %        APPENDIX               % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %

\clearpage
\appendix
\appendixpage

\section{Deriving the saddle-point problem}

To derive the saddle-point formulation of our optimization problem, we used two fundamental results of convex analysis to "dualize" the data fidelity constraint as well as the VTV term. The transformation of the $\VTV$ term follows straightforwardly from the definition of the so-called "dual-norm," which is defined by
%
\begin{equation}\label{eq:dualnorm}
\begin{split}
\|x\|' = \sup_z &\,\, \inner{x,z} \,\,\, \mathrm{s.t.} \,\, \|z\|\leq 1.
\end{split}
\end{equation}
%
Every norm $\|.\|$ has an associated dual-norm $\|.\|'$, that obeys this relationship. The $\ell_2$ norm utilized in the scalar TV is self-dual, while the nuclear norm and spectral norm form a dual pair. 

\paragraph{Dualizing the channel-by-channel TV}
Recall that the channel-by-channel TV can be written as 
\begin{equation}
\TV_S(\uu) = \sum_{m,\ell} \left\|\left[D\uu_m \right]_\ell \right\|_2.
\end{equation}
The definition of the dual-norm allows us to rewrite this as 
\begin{equation}
\begin{split}
\TV_S(\uu) &= \sum_{m,\ell} \sup_\zz \,\, \inner{\left[D\uu_m\right]_\ell,\left[\zz_m\right]_\ell} \,\, \mathrm{s.t.} \,\, \|\left[\zz_m\right]_\ell\|_2 \leq 1\\
&= \sup_\zz \,\, \inner{D\uu,\zz} \,\, \mathrm{s.t.} \,\, \left\|\left[\zz_m\right]_\ell \right\|_2 \leq 1 \,\,\, \forall m \in 1,\ldots,\mathcal{M} \,\,\,; \ell \in 1,\ldots,\mathcal{L}
\end{split}
\end{equation}
%
which is equivalent to (\ref{eq:vtvsdual}).

\paragraph{Dualizing the proposed TNV}

The proposed TNV is 
\begin{equation}
\TNV(\bm{u}) = \sum_{\ell=1}^{\mathcal{L}} \|\left[D\bm{u}\right]_\ell\|_\star.
\end{equation}
%
Substituting the definition of the dual-norm, we get
%
\begin{equation}
\begin{split}
\TNV(\uu) &= \sum_{\ell} \sup_\zz \,\, \inner{\left[D\uu\right]_\ell,\left[\zz\right]_\ell} \,\, \mathrm{s.t.} \,\, \sigma_{\mathrm{max}}\left(\left[\zz\right]_\ell\right) \leq 1\\
&= \sup_\zz \,\, \inner{D\uu,\zz} \,\, \mathrm{s.t.} \,\, \sigma_{\mathrm{max}}\left(\left[\zz\right]_\ell\right) \leq 1 \,\,\, \forall \ell \in 1,\ldots,\mathcal{L},
\end{split}
\end{equation}

which is equivalent to (\ref{eq:vtvndual}). We have used the fact that the spectral norm (maximum singular value) is dual to the nuclear norm. 
%
\paragraph{Dualizing the data-fidelity term using the Fenchel conjugate}

To dualize the data-fidelity term, we make use of the Fenchel conjugate. For a convex, lower semi-continuous function $f$, the conjugate $f^*$ is defined by 
%
\begin{equation}\label{eq:fenchel}
f^*(x) = \sup_{x'} \,\, \inner{x,x'} - f(x'),
\end{equation}
%
where it is also true that $\left(f^*\right)^* = f$. Using this definition, it is easy to show that the following functions form a conjugate pair,
\begin{equation}
\delta_\epsilon(x) \Leftrightarrow \epsilon \|W^{-1/2}x\|_2
\end{equation}
where $\delta_\epsilon(x)$ is the indicator function, defined by (\ref{eq:dataconstraint}). This along with the definition of the Fenchel conjugate leads directly to Eqn. \ref{eq:dataconjugate}.

\section{The proximal map of $\epsilon\sigma\|W^{{-1/2}}\qq\|_2$}
The update equations for Algorithm \ref{alg:VTVupdates} involve evaluating the proximal map of $\epsilon\sigma F_q \equiv \epsilon\sigma\|W^{-1/2}\qq\|_2$, which does not admit a closed form. However, we will now show how it can still be efficiently evaluated using a simple root-finding procedure. First we explicitly write out the proximal mapping as 
%
\begin{equation}
\prox{\epsilon\sigma F_q}{\qq'} = \argmin_\qq \,\, \|W^{-1/2}\qq\|_2 + \frac{1}{2\epsilon\sigma}\|\qq-\qq'\|_2^2.
\end{equation}
%
Next, we compare this optimization problem to a slightly easier problem, 
%
\begin{equation}
\argmin_\qq \,\, \frac{1}{2}\|W^{-1/2}\qq\|_2^2 + \frac{1}{2\lambda\sigma}\|\qq-\qq'\|_2^2,
\end{equation}
%
and note that for some choice of $\epsilon$ and $\lambda$, these objectives both have the same set of level curves. We can find the solution to this problem by setting the gradient equal to zero, and for a symmetric weighting matrix $W$, it is given by 
\begin{equation}
\qq^* = \left(I + \sigma\lambda W^{-1}\right)^{-1}\qq'.
\end{equation}
If the gradient of this easier problem is equal to that of the original problem, then $\qq^*$ must be an optimizer of the original problem as well. Momentarily, we ignore the non-differentiable point $\qq=\bm{0}$, and set the gradients equal, which yields
%
\begin{equation}\label{eq:qroot}
\begin{split}
\epsilon &= \lambda \|W^{-1/2}\qq^*\|_2.\\
&= \lambda \|W^{-1/2}\left(I + \sigma \lambda W^{-1}\right)^{-1}\qq'\|_2.
\end{split}
\end{equation}
%
In this work, we will only consider the case where $W$ is a diagonal matrix, which allows us to simplify this expression.
%
\begin{equation}\label{eq:lambdaroot}
\epsilon^2 - \lambda ^2\sum_{j=1}^J \frac{W_{jj}}{\left(W_{jj} + \sigma\lambda \right)^2}[\qq']^2_j  = 0.
\end{equation}
%
The index $j$ is a linear index over every element of $\qq$. Now, we just need to perform a 1D search for a positive value of $\lambda$ that satisfies this equation, which can be done very efficiently using Newton's method or a variety of other algorithms. Once this root, $\lambda^*$, is found, we simply plug it back into (\ref{eq:qroot}):
%
\begin{equation}
\prox{\epsilon\sigma F_q}{\qq'} = \left\{\begin{array}{lr} \left(I + \sigma\lambda^* W^{-1}\right)^{-1}\qq' & \|W^{1/2}\qq'\|_2 > \sigma\epsilon \\
0 & \|W^{1/2}\qq'\|_2 \leq \sigma\epsilon \end{array}\right..
\end{equation}
%
Note that because of the non-differentiability of $\|W^{-1/2}\qq\|_2$ at $\qq=0$, there will not always be a solution to equation \ref{eq:lambdaroot}. Specifically, if $\sigma\epsilon > \|W^{1/2}\qq'\|_2$, then there is no viable root. In this case, we need not perform the root-finding procedure, because the optimum must occur at $\bm{q} = \bm{0}$. 

\section{Implementation of the projection operators $\Pi_{\mathcal{S}/N}$}

\paragraph{$\bm{\Pi_\mathcal{S}}$}

Now, we will describe the Euclidean projection onto the set $\mathcal{S}$ defined in (\ref{eq:setS}). Consider projecting $\zz\in\mathbb{R}^{(\mathcal{M}\times 2)\cdot\mathcal{L}}$ onto $\mathcal{S}$. We can define this operation element-wise on each vector $[\zz_m]_\ell \in \mathbb{R}^2$ corresponding to the $m^\text{th}$ spectral channel and the $\ell^\text{th}$ pixel-location. The projection of this element is given by 
%
\begin{equation}
\Pi_{\mathcal{S}}\left([\zz_m]_\ell \right) = \left\{ \begin{array}{lr} [\zz_m]_\ell  & \|[\zz_m]_\ell\|_2 \leq 0 \\
\frac{[\zz_m]_\ell}{\|[\zz_m]_\ell\|_2} & \text{otherwise}   \end{array}\right..
\end{equation}
For every image channel $m$ and pixel location $\ell$, we simply project the vector $[\zz_m]_\ell$ onto the unit ball. 

\paragraph{$\bm{\Pi_N}$}
The projection onto set $N$ defined by (\ref{eq:setN}) is slightly more complicated. We define this operation element-wise on each $\mathcal{M}\times2$ matrix  $[\zz]_\ell$. This time we need to threshold the eigenvalues of $[\zz]_\ell$. This projection can be succinctly described by
%
\begin{equation}\label{eq:svd}
\Pi_N\left([\zz]_\ell \right) = U\Sigma_pV^T,
\end{equation}
%
where $U\Sigma V^T$ is the SVD of $[\zz]_\ell$, and 
%
\begin{equation}
\Sigma_p = \mathrm{sgn}\left(\Sigma \right)\mathrm{min}\left(|1|,|\Sigma| \right)
\end{equation}
is the thresholded version of $\Sigma$. To form $\Sigma_p$ the singular values on the diagonal of $\Sigma$ are simply thresholded so that their magnitude does not exceed 1. An equivalent projection formula that leads to a much more computationally efficient solution is given by 
%
\begin{equation}\label{eq:eigenproject}
\Pi_N\left([\zz]_\ell \right) = [\zz]_\ell V \Sigma^\dagger\Sigma_pV^T.
\end{equation}
%
The quantity $\Sigma^\dagger$ is the pseudo-inverse of $\Sigma$. Since we are only working with 2 spatial dimensions the matrix $V$ will by a $2\times2$ square matrix. Therefore, to compute the projection according to (\ref{eq:eigenproject}) we only need to compute the eigenvalues and eigenvectors of a $2\times 2$ matrix, for which a very simple closed form is available. This can also be done very efficiently for 3D images, where $V$ will be $3\times 3$. Because of this, the update equations that result from the proposed vectorial TV are only trivially more expensive than the channel-by-channel TV, and the projection/backprojection operations are likely to swamp this difference. 

\begin{comment}
\section{A graphical explanation of the noise-balancing procedure}


\paragraph{VTV denoising of a 1D multi-channel signal without noise-balancing.}

Consider a noisy two-channel time signal $\bm{g(t)}$, that we wish to denoise. For illustrative purposes, we consider a scenario where the second signal channel has 10$\times$ the standard deviation.
%
 
\begin{figure}[htbp][!ht]
\centering
\begin{tikzpicture}[scale=3.0]
\begin{scope}
\draw[->] (-1,0) -- (1,0); \node[below right] at (1,0) {$t$};
\node[red] at (0,0.6) {$g_1(t)$};
\draw[red,line width=0.03cm] plot [domain=-1:1,samples=500] (\x,{0.4*(\x>=0.0)});
\draw[red,line width=0.03cm,opacity=0.3] plot [domain=-1:1,samples=500] (\x,{0.4*(\x>=0.0)+0.04*rand});
\end{scope}
\begin{scope}[shift = {(2.5,0)}]
\draw[->] (-1,0) -- (1,0); \node[below right] at (1,0) {$t$};
\node[blue] at (0,0.9) {$g_2(t)$};
\draw[blue,line width=0.03cm] plot [domain=-1:1,samples=500] (\x,{0.4*(\x>=0.0)});
\draw[blue,line width=0.03cm,opacity=0.3] plot [domain=-1:1,samples=500] (\x,{0.4*(\x>=0.0)+0.4*rand});
\end{scope}
\end{tikzpicture}
\begin{tikzpicture}[scale=3.0]
\begin{scope}
\draw[->] (-1,0) -- (1,0); \node[below right] at (1,0) {$t$};
\node[red] at (0,0.6) {$g'_1(t)$};
\draw[red,line width=0.03cm,opacity=0.3] plot [domain=-1:1,samples=500] (\x,{0.04*rand});
\draw[red,line width=0.03cm] (0,0)--(0,0.4); \fill[red] (0,0.4) circle (0.03cm);
\end{scope}
\begin{scope}[shift = {(2.5,0)}]
\draw[->] (-1,0) -- (1,0); \node[below right] at (1,0) {$t$};
\node[blue] at (0,0.6) {$g'_2(t)$};
\draw[blue,line width=0.03cm,opacity=0.3] plot [domain=-1:1,samples=500] (\x,{0.4*rand});
\draw[blue,line width=0.03cm] (0,0)--(0,0.4); \fill[blue] (0,0.4) circle (0.03cm);
\end{scope}
\end{tikzpicture}
\caption{Plot of a 2-channel, noisy time signal (top) and its derivative (bottom). The dark lines indicate the true signal, while the translucent lines indicate the noise-corrupted measurements.}
\end{figure}
%
\noindent We can use the VTV to denoise this signal by minimizing the following functional.
\begin{equation}\label{eq:tvdenoise}
\min_{\bm{f(t)}} \quad \frac{1}{2}\int \|\bm{f}(t)-\bm{g}(t)\|_2^2 \, \mathrm{dt} + \lambda \int \|D \bm{f}\|_\star \, \mathrm{dt}
\end{equation}
%
where $\bm{f(t)} = \begin{bmatrix} f_1(t) & f_2(t)\end{bmatrix}^T$, $\bm{g(t)} = \begin{bmatrix} g_1(t) & g_2(t)\end{bmatrix}^T$, $D\bm{f} = \begin{bmatrix} f_1'(t) & f'_2(t)\end{bmatrix}^T$. Note that for 1D signals, all vectorial TV's based on Schatten norms of the Jacobian are equivalent, so we will simply refer to this regularizer as \emph{the VTV}. A technique that is sometimes employed is to perform the denoising in the derivative domain, where the signal is inherently sparse, and then to recover the signal via integration. Fairly trivial assumptions are required in order to transform between the signal and its time-derivative. The optimization over the derivative $\bm{f'(t)}$ is given by (\ref{eq:tvdenoisederivative}).
%
\begin{equation}\label{eq:tvdenoisederivative}
\min_{\bm{f'(t)}} \quad \frac{1}{2}\int \|\bm{f'(t)} - \bm{g'(t)}\|_2^2 \,\mathrm{dt}+ \lambda \int \|\bm{f'(t)}\|_2\,\mathrm{dt}
\end{equation}
%
The reason for working with the derivatives is that the above problem admits a closed form solution. This variational problem is identical to the well-known \emph{group lasso}, and its solution is given by the multi-dimensional shrinkage/thresholding operator (\ref{eq:vectorshrink}).

\begin{equation}\label{eq:vectorshrink}
\begin{split}
\bm{f'(t)} &= \frac{\bm{g'(t)}}{\|\bm{g'(t)}\|_2}\mathrm{shrink}_\lambda(\|\bm{g'(t)}\|_2)\\
\mathrm{shrink}_\lambda(x) &= \left\{ \begin{array}{lr} x-\lambda & x < -\lambda\\
												-x+\lambda & x>\lambda \\
												0 & -\lambda \leq x \leq \lambda \end{array} \right.
\end{split}
\end{equation}


\noindent In words, vectors $\bm{g'(t)}$ that have magnitude less than $\lambda$ are thresholded to zero, while all other vectors are decreased in length by an amount $\lambda$.

In order to better visualize the shrinkage/thresholding operator, we make a parametric plot of the 2-channel signal. The mean signal is depicted by the dark purple line, while the envelope of the noise-magnitude is characterized by the translucent, purple ellipse. For an arbitrary choice of $\lambda$, the dotted gray circle indicates the vectors that will be thresholded to zero. In this case, the VTV denoising does help preserve some of the signal from the noisier (blue) channel, but we are still faced with the problem of either over-smoothing the cleaner channel or under-smoothing the noisier channel. 

\begin{figure}[htbp][!ht]
\centering
\begin{tikzpicture}[scale=4.0]
\begin{scope}
\draw[->,opacity=0.5] (-0.7,0) -- (0.7,0); \node[below right,red] at (0.7,0) {$g'_1(t)$};
\draw[->,opacity=0.5] (0,-0.7) -- (0,0.7); \node[below right,blue] at (0,0.7) {$g'_2(t)$};
\draw[red!50!blue,line width=0.04cm] (0,0) -- (0.4,0.4);\fill[red!50!blue] (0.4,0.4) circle (0.03cm);
\fill[red!50!blue,opacity=0.25] (0,0) circle [x radius=0.04, y radius=0.4];
\draw[->] (-115:0.5) -- (-95:0.35);\node[left] at (-115:0.5) {noise envelope};
\begin{scope}[opacity=0.3]
%\draw[-,blue,line width=0.04cm] (0,0) -- (0,0.4);\fill[blue] (0,0.4) circle (0.03cm);
%\draw[-,red,line width=0.04cm] (0,0) -- (0.4,0); \fill[red] (0.4,0) circle (0.03cm);
\end{scope}
\draw[dashed,opacity=0.5] (0,0) circle (0.3); \node[opacity=0.5] at (-45:0.35) {$\lambda$};
\end{scope}
\end{tikzpicture}
\caption{The true signal is denoted by the arrow in purple, while the ellipse characterizes the magnitude of the noise in both channels. The y-radius of the ellipse is larger because the second channel is noisier. The dotted gray circle indicates the amount of shrinkage that will occur with the denoising.}
\end{figure}


\FloatBarrier

\paragraph{Improving VTV denoising by balancing the per-channel noise.}

Note that this VTV denoising works better, if we first scale the signal channels in such a way that the magnitude of the noise is the same in both channels. The signals and their derivatives after noise-scaling are depicted below.
%
\begin{figure}[htbp][!ht]
\centering
\begin{tikzpicture}[scale=3.0]
\begin{scope}
\draw[->] (-1,0) -- (1,0); \node[below right] at (1,0) {$t$};
\node[red] at (0,0.6) {$g_1(t)$};
\draw[red,line width=0.03cm] plot [domain=-1:1,samples=500] (\x,{0.4*(\x>=0.0)});
\draw[red,line width=0.03cm,opacity=0.3] plot [domain=-1:1,samples=500] (\x,{0.4*(\x>=0.0)+0.04*rand});
\end{scope}
\begin{scope}[shift = {(2.5,0)}]
\draw[->] (-1,0) -- (1,0); \node[below right] at (1,0) {$t$};
\node[blue] at (0,0.3) {$g_2(t)$};
\draw[blue,line width=0.03cm] plot [domain=-1:1,samples=500] (\x,{0.04*(\x>=0.0)});
\draw[blue,line width=0.03cm,opacity=0.3] plot [domain=-1:1,samples=500] (\x,{0.04*(\x>=0.0)+0.04*rand});
\end{scope}
\end{tikzpicture}
\begin{tikzpicture}[scale=3.0]
\begin{scope}
\draw[->] (-1,0) -- (1,0); \node[below right] at (1,0) {$t$};
\node[red] at (0,0.6) {$g'_1(t)$};
\draw[red,line width=0.03cm,opacity=0.3] plot [domain=-1:1,samples=500] (\x,{0.04*rand});
\draw[red,line width=0.03cm] (0,0)--(0,0.4); \fill[red] (0,0.4) circle (0.03cm);
\end{scope}
\begin{scope}[shift = {(2.5,0)}]
\draw[->] (-1,0) -- (1,0); \node[below right] at (1,0) {$t$};
\node[blue] at (0,0.6) {$g'_2(t)$};
\draw[blue,line width=0.03cm,opacity=0.3] plot [domain=-1:1,samples=500] (\x,{0.04*rand});
\draw[blue,line width=0.03cm] (0,0)--(0,0.04); \fill[blue] (0,0.04) circle (0.03cm);
\end{scope}
\end{tikzpicture}
\caption{Time signals and their derivatives after noise-balancing}
\end{figure}
%
To see why the noise-balancing is advantageous, we look again at the parametric representation of the scaled signal.
%
\begin{figure}[htbp][!ht]
\centering
\begin{tikzpicture}[scale=5.0]
\begin{scope}
\draw[->,opacity=0.5] (-0.3,0) -- (1,0); \node[below right,red] at (1,0) {$g'_1(t)$};
\draw[->,opacity=0.5] (0,-0.3) -- (0,0.3); \node[below right,blue] at (0,0.4) {$g'_2(t)$};
\draw[red!50!blue,line width=0.04cm] (0,0) -- (0.8,0.08); \fill[red!50!blue] (0.8,0.08) circle (0.02cm);
\fill[red!50!blue,opacity=0.25] (0,0) circle [x radius=0.08, y radius=0.08];
\draw[->] (-135:0.25) -- (-135:0.08);\node[below left] at (-135:0.25) {noise envelope};
\begin{scope}[opacity=0.3]
%\draw[blue,line width=0.04cm] (0,0) -- (0,0.08);\fill[blue] (0,0.08) circle (0.02cm);
%\draw[red,line width=0.04cm] (0,0) -- (0.8,0); \fill[red] (0.8,0.0) circle (0.02cm);
\end{scope}
\draw[dashed,opacity=0.5] (0,0) circle (0.2); \node[opacity=0.5] at (-45:0.25) {$\lambda$};
\end{scope}
\end{tikzpicture}
\caption{Parametric signal plot after noise-balancing. Now the noise distribution is isotropic, so we can easily pick a $\lambda$ to suppress the noise without severely biasing the actual signal. }
\end{figure}
%
Now that the noise is evenly distributed between the 2 channels we can easily select a $\lambda$ value that gets rid of most of the noise while also preserving the majority of the signal. Even the extremely noisy second channel is well-recovered because of the tight coupling to the less noisy (red) channel. This suggests that noise suppression may be improved if, prior to reconstruction, the data channels are scaled so that the noise magnitude is approximately, uniformly distributed among all channels.
\end{comment}

\end{document}
paper_draft4_revised.bbl
\begin{thebibliography}{10}

\bibitem{Rigie2014}
D.~S. Rigie and P.~J. {La Rivi\`{e}re}, ``{A generalized vectorial
  total-variation for spectral CT reconstruction},'' {\em Proceedings of the
  Third International Conference on Image Formation in X-ray Computed
  Tomography}, vol.~1, no.~1, pp.~9--12, 2014.

\bibitem{Holt2014}
K.~Holt, ``{Total Nuclear Variation and Jacobian Extensions of Total Variation
  for Vector Fields.},'' {\em IEEE transactions on image processing : a
  publication of the IEEE Signal Processing Society}, vol.~23, pp.~3975--3989,
  June 2014.

\bibitem{Rudin1992}
L.~I. Rudin, S.~Osher, and E.~Fatemi, ``{Nonlinear total variation based noise
  removal algorithms},'' {\em Physica D: Nonlinear Phenomena}, vol.~60,
  pp.~259--268, 1992.

\bibitem{Chan1998}
T.~F. Chan and C.~K. Wong, ``{Total variation blind deconvolution.},'' {\em
  IEEE transactions on image processing : a publication of the IEEE Signal
  Processing Society}, vol.~7, pp.~370--375, 1998.

\bibitem{Komatsu2005}
T.~Komatsu and T.~Saito, ``{Super-resolution sharpening-demosaicking method for
  removing image blurs caused by An optical low-pass filter},'' in {\em
  Proceedings - International Conference on Image Processing, ICIP}, vol.~1,
  pp.~845--848, 2005.

\bibitem{Candes2006}
E.~J. Cand\`{e}s, J.~Romberg, and T.~Tao, ``{Robust Uncertainty Principles :
  Exact Signal Frequency Information},'' {\em IEEE Transactions on Information
  Theory}, vol.~52, no.~2, pp.~489--509, 2006.

\bibitem{Sidky2006}
E.~Sidky, C.~Kao, and X.~Pan, ``{Accurate image reconstruction from few-views
  and limited-angle data in divergent-beam CT},'' {\em Journal of X-ray Science
  and Technology}, vol.~14, pp.~119--139, 2006.

\bibitem{Bian2013}
J.~Bian, J.~Wang, X.~Han, E.~Y. Sidky, L.~Shao, and X.~Pan,
  ``{Optimization-based image reconstruction from sparse-view data in
  offset-detector CBCT.},'' {\em Physics in medicine and biology}, vol.~58,
  pp.~205--30, 2013.

\bibitem{Song2007}
J.~Song, Q.~H. Liu, G.~A. Johnson, and C.~T. Badea, ``{Sparseness prior based
  iterative image reconstruction for retrospectively gated cardiac
  micro-CT.},'' {\em Medical physics}, vol.~34, pp.~4476--4483, 2007.

\bibitem{Tian2011}
Z.~Tian, X.~Jia, and K.~Yuan, ``{Low-dose CT reconstruction via edge-preserving
  total variation regularization},'' {\em Physics in medicine and \ldots},
  vol.~56, no.~18, pp.~5949--5967, 2011.

\bibitem{Sapiro1997}
G.~Sapiro, ``{Color Snakes},'' {\em Computer Vision and Image Understanding},
  vol.~68, pp.~247--253, 1997.

\bibitem{Blomgren1998}
P.~Blomgren and T.~F. Chan, ``{Color TV: Total variation methods for
  restoration of vector-valued images},'' {\em IEEE Transactions on Image
  Processing}, vol.~7, pp.~304--309, 1998.

\bibitem{Goldluecke2012}
B.~Goldluecke, E.~Strekalovskiy, and D.~Cremers, ``{The Natural Vectorial Total
  Variation Which Arises from Geometric Measure Theory},'' {\em SIAM Journal on
  Imaging Sciences}, vol.~5, pp.~537--563, 2012.

\bibitem{Lefkimmiatis2013}
S.~Lefkimmiatis, A.~Roussos, M.~Unser, and P.~Maragos, ``{Convex
  Generalizations of Total Variation Based on the Structure Tensor with
  Applications to Inverse Problems},'' in {\em Scale Space and Variational
  Methods in Computer Vision 4th International Conference, SSVM 2013, Schloss
  Seggau, Leibnitz, Austria, June 2-6, 2013. Proceedings} (A.~Kujiper,
  K.~Bredies, T.~Pock, and H.~Bischof, eds.), pp.~48--60, Springer Berlin
  Heidelberg, 4th~ed., 2013.

\bibitem{keren1998}
D.~Keren and A.~Gotlib, ``{Denoising Color Images Using Regularization and
  "Correlation Terms"},'' {\em Journal of Visual Communication and Image
  Representation}, vol.~9, no.~4, pp.~352--365, 1998.

\bibitem{Holt2011}
K.~M. Holt, ``{Angular regularization of vector-valued signals},'' in {\em
  ICASSP, IEEE International Conference on Acoustics, Speech and Signal
  Processing - Proceedings}, pp.~1105--1108, 2011.

\bibitem{Fessler2002}
J.~A. Fessler, I.~A. Elbakri, P.~Sukovic, and N.~H. Clinthorne,
  ``{Maximum-likelihood dual-energy tomographic image reconstruction},'' {\em
  SPIE Medical Imaging}, vol.~4684, pp.~38--49, 2002.

\bibitem{Long2014}
Y.~Long and J.~A. Fessler, ``{Multi-Material Decomposition Using Statistical
  Image Reconstruction for Spectral CT.},'' {\em IEEE transactions on medical
  imaging}, vol.~33, pp.~1614--26, Aug. 2014.

\bibitem{Alvarez1976}
R.~E. Alvarez and A.~Macovski, ``{Energy-selective Reconstructions in X-ray
  Computerized Tomography},'' {\em Phys Med Biol}, vol.~21, pp.~733--744, 1976.

\bibitem{Chambolle2011}
A.~Chambolle and T.~Pock, ``{A first-order primal-dual algorithm for convex
  problems with applications to imaging},'' {\em Journal of Mathematical
  Imaging and Vision}, vol.~40, pp.~120--145, 2011.

\bibitem{Segars2010}
W.~P. Segars, G.~Sturgeon, S.~Mendonca, J.~Grimes, and B.~M.~W. Tsui, ``{4D
  XCAT phantom for multimodality imaging research.},'' {\em Medical physics},
  vol.~37, pp.~4902--4915, 2010.

\bibitem{Han2012}
X.~Han, J.~Bian, E.~L. Ritman, E.~Y. Sidky, and X.~Pan, ``{Optimization-based
  reconstruction of sparse images from few-view projections},'' {\em Physics in
  Medicine and Biology}, vol.~57, pp.~5245--5273, 2012.

\bibitem{Xu2014}
Q.~Xu, A.~Sawatzky, M.~A. Anastasio, and C.~O. Schirra, ``{Sparsity-regularized
  image reconstruction of decomposed K-edge data in spectral CT.},'' {\em
  Physics in medicine and biology}, vol.~59, pp.~N65--79, May 2014.

\bibitem{Cai2010}
J.~Cai, E.~Cand\`{e}s, and Z.~Shen, ``{A singular value thresholding algorithm
  for matrix completion},'' {\em SIAM Journal on Optimization}, pp.~1--28,
  2010.

\bibitem{Candes2011}
E.~Cand\`{e}s, X.~Li, Y.~Ma, and J.~Wright, ``{Robust principal component
  analysis?},'' {\em Journal of the ACM (JACM)}, 2011.

\bibitem{Sidky2012}
E.~Y. Sidky, J.~H. J{\o}rgensen, and X.~Pan, ``{Convex optimization problem
  prototyping for image reconstruction in computed tomography with the
  Chambolle–Pock algorithm},'' {\em Physics in Medicine and Biology},
  vol.~57, pp.~3065--3091, 2012.

\bibitem{Parikh2013}
N.~Parikh and S.~Boyd, ``{Proximal algorithms},'' {\em Foundations and Trends
  in Optimization}, vol.~1, no.~3, pp.~123--231, 2013.

\bibitem{Bouman1996}
C.~A. Bouman and K.~Sauer, ``{A unified approach to statistical tomography
  using coordinate descent optimization},'' {\em IEEE Transactions on Image
  Processing}, vol.~5, pp.~480--492, 1996.

\bibitem{Roessl2009}
E.~Roessl and C.~Herrmann, ``{Cram\'{e}r-Rao lower bound of basis image noise
  in multiple-energy x-ray imaging.},'' {\em Physics in medicine and biology},
  vol.~54, pp.~1307--18, Mar. 2009.

\bibitem{Schirra2013}
C.~O. Schirra, E.~Roessl, T.~Koehler, B.~Brendel, A.~Thran, D.~Pan, and
  R.~Anastasio, Mark A.and~Proksa, ``{Statistical reconstruction of material
  decomposed data in spectral CT.},'' {\em IEEE transactions on medical
  imaging}, vol.~32, pp.~1249--57, July 2013.

\bibitem{Sawatzky2014}
A.~Sawatzky, Q.~Xu, C.~O. Schirra, and M.~A. Anastasio, ``{Proximal ADMM for
  Multi-Channel Image Reconstruction in Spectral X-ray CT.},'' {\em IEEE
  transactions on medical imaging}, vol.~33, pp.~1657--1668, Aug. 2014.

\end{thebibliography}
