text:
Home > African & African American Studies > Senior Theses > 10
African & African American Studies Senior Theses
Narrative of Flight: Camden, New Jersey in the Post-War Era
Chuck Hovanic, Fordham UniversityFollow
Commuters boarding the Speed Line in the small town of Lindenwold, New Jersey witness a rapid change in landscape as their train whisks them towards Philadelphia. Each stop on the train represents a distinct stage of suburban development ranging from the sprawling McMansions of Voorhees to the tightly packed two-family houses of Collingswood. The passengers that board a the various stations accurately represents the demographics of each town. This broad range of people is treated to a memorable, yet brief glimpse of poverty when the train enter the City of Camden. Gliding above decaying row homes and empty streets, passengers stare at the mot dangerous city in America from the relative safety of the train. This fleeting encounter with Camden that is witnessed daily by thousands of South Jersey residents, represents the extant to which mot people interact with this embattled city. Only memories, told by grandparents and parents of these commuters, serve to illustrate the Camden that once was.
Hovanic, Chuck, "Narrative of Flight: Camden, New Jersey in the Post-War Era" (2013). African & African American Studies Senior Theses. 10.
https://research.library.fordham.edu/aaas_senior/10
African American Studies Commons
Undergraduate Research Commons

redpajama_set_name: RedPajamaCommonCrawl


text:
\section{Introduction}
The field of sparse recovery has a wide range of applications in many different areas such as medical imaging \cite{LusDonPau,CTRef}, astronomy \cite{Astronomy1,Astronomy2}, electron microscopy \cite{EM1} etc. One of the great successes in this area are new developments of multiscale sparsifying transforms since the invention of wavelets and the wavelet transform. Indeed wavelets are known to compress natural images very effectively since most natural images are sparse in a wavelet domain. However, they lack in directional sensitivity and are therefore not optimal for images that are governed by curvilinear structures. For precisely that reason, almost a decade ago directional systems such as curvelets \cite{CanDon} and shearlets \cite{KitKutLim, GuoKutLab2006} have been created to overcome this deficit. We will briefly recall the concepts of wavelets and shearlets in Section \ref{sec:Section2} as both systems will play a significant role in the upcoming content of this paper.

The field of sparse recovery is predominantly influenced by the development of  \emph{compressed sensing} \cite{CanRomTao, Don}, a theory that guarantees the recovery of \emph{sparse signals} from incomplete measurements under the assumptions of \emph{sparsity} and \emph{incoherence}. These sparse signals are typically obtained by solving a convex optimization problem of the form
\begin{align}
    \min_u \| \Psi  u \|_1 \quad \text{ subject to } \quad Au = y, \label{eq:AnaMin}
\end{align}
where $u$ is the object of interest, $\Psi$ is a sparsifying transform, $A$ is a matrix representing the measurement process, and $y$ are the resulting measurements. If $u$ is already \emph{sparse} it suffices to let $\Psi$ to be the identity. Otherwise \eqref{eq:AnaMin} is called \emph{Basis Pursuit in the analysis formulation}.
 Very recently, Ahmad and Schniter considered a generalized variant of \eqref{eq:AnaMin} in \cite{AhmSch}, where they have used not only one sparsifying transform, but a composition of several sparsifying transforms. Furthermore, in order to improve the reconstruction quality the authors have combined their ideas with reweighted $\ell^1$-minimization into their framework. Indeed, the concept of \emph{reweighted $\ell^1$-minimization} introduced by Cand\`{e}s et al. in \cite{CanWakBoy} further promotes the sparsity of the recovered signal by iteratively updating a weighting matrix in the minimization problem. More precisely, one solves
 \begin{align}
\min_u \| W_k\Psi u \|_1 \quad \text{ subject to } \quad Au = y, \label{eq:AnaMinRe}
\end{align}
iteratively for $k=1,2,\dots$ where after each iteration $k$ the diagonal weighting matrix $W_k$ is updated according to the sparsity structure of the current solution $u_{k+1}$ of \eqref{eq:AnaMinRe}. The role of $W_k$ is to mimic the actual sparsity structure of the true signal that one wishes to recover,. We will recap the ideas of reweighted $\ell^1$ later in Section \ref{sec:Section2} in more detail as this will also be one of the main ingredients of the algorithm that we propose in this work.

In order to find solutions of \eqref{eq:AnaMin} or approximations of such solutions,
many different possibilities available in the literature. In this work, we focus on the \emph{split Bregman} algorithm \cite{SB,SB2}, which besides minor differences, is a reinvention of the \emph{alternating direction method of multipliers} (ADMM) \cite{GabMer, EckBer, BoyADMM}. It can also be seen as \emph{Douglas Rachford splitting} of the dual problem. For details on the relationship of these algorithms we refer the interested reader to \cite{Setzer,Esser}.  The split Bregman algorithm transforms the constrained problem \eqref{eq:AnaMin} into an unconstrained formulation and by introducing splitting variables they break the original problem down into, hopefully, easier ones. A key ingredient is then the so-called \emph{soft-thresholding}, \emph{shrinking} or \emph{proximal mapping} which gives a closed-form solution to some of the subproblems. One of the great advantages of ADMM and in particular split Bregman is that they are easily derived and very flexible in terms of  multiple regularizers. Furthermore, as we shall explain in this paper, it can be greatly combined with the idea of reweighted $\ell^1$.

 It is often beneficial to consider a second sparsity promoting regularizer such as total variation (TV) \cite{SB,RecPF} to reduce artifacts coming from the sparsifying transform $\Psi$ in the
solutions of sparse imaging problems  that are obtained via  \eqref{eq:AnaMin}. TV was initially proposed for denoising problems \cite{RudOshFat} and is very well established in image processing by now. However, it has been noticed that severe staircasing artifacts may  appear in images recovered by TV regularized reconstructions for large noise levels. In order to overcome this issue the authors of \cite{BreKunPoc} have introduced \emph{total generalized variation (TGV)} which is a generalization of TV to higher order derivatives. Since then TGV has been used in many applications \cite{BloTobUecFra,KnoBrePocRud} as a regularizer for inverse problems.

 The consideration of a joint regularization scheme using shearlets and TGV has already been done by Guo et al. in \cite{GuoQinYin} in order to solve problems of the form \eqref{eq:AnaMin}. Our main contribution is to further exploit the general structure of the (multilevel) sparsity by combining iterative reweighting and adaptive multilevel weights associated to any multiscale transform, not only shearlets. This is also related to the approach in \cite{AhmSch}, but the algorithm and the conclusions derived in this work are different. In order to solve the resulting multilevel reweighted $\ell^1$ problem we make use of the flexibility of split Bregman by directly incorporating adaptive multilevel thresholds into the subproblem that is solved by a simple thresholding step. This results in an algorithm that comes with very little additional computational cost and almost automatically chosen regularization parameters. In particular, as we will show in this paper, the results are greatly improved compared to the non-reweighted analogue. We wish to mention, that the reweighting approach that we are considering in this paper is not to be confused with weighted $\ell^1$-minimization. Weighted $\ell^1$-minimization requires addition a priori knowledge to carefully design effetive weights. whereas in reweighting the weights are adaptive and do not require a priori knowledge. Furthermore, reweighted $\ell^1$-minimization is an iterative scheme and the minimization problem has to be solved several times (with updated weights) which is not the case for weighted $\ell^1$.

 In Figure \ref{fig:abstract} we give a motivation and a first glance for the possible benefit of combining reweighting methods with, in this particular case, the wavelet transform. The reconstructions shown in Figure \ref{fig:abstract}are obtained from partial Fourier measurements of an synthetic test image which has parts that are certainly sparse in a wavelet dictionary and parts that are very well suited for TGV. In Section \ref{sec:Section4} we will present all details of the numerical implementation. Note that in Figure \ref{fig:abstract} the proposed method reduces the artifacts while still being able to reconstruct fine details using only $10.28\%$ of Fourier measurements.

 \begin{figure*}
 \includegraphics[width=.24\textwidth]{./truthA.png}
 \includegraphics[width=.24\textwidth]{./zeroFillA.png}
 \includegraphics[width=.24\textwidth]{./29WL1TGV.png}
 \includegraphics[width=.24\textwidth]{./29WIRL1TGV.png}
\\[.5ex]
 \includegraphics[width=.24\textwidth]{./ztruthA.png}
 \includegraphics[width=.24\textwidth]{./zzeroFillA.png}
 \includegraphics[width=.24\textwidth]{./z29WL1TGV.png}
 \includegraphics[width=.24\textwidth]{./z29WIRL1TGV.png}
 \caption{Reconstructions from 25 radial lines ($10.28\%$) through the k-space origin. \textbf{First column}: Original image with zoom. \textbf{Second column}: Reconstruction with inverse Fourier transform. Relative error: 0.146. Structured similarity index: 0.735. \textbf{Third column}: Reconstruction with redundant Daubechies 4 wavelets and total generalized variation regularizer without reweighting. Relative error: 0.060. Structured similarity index: 0.896. \textbf{Fourth column}: Reconstruction with the proposed method. Relative error: 0.031. Structured similarity index: 0.951.}\label{fig:abstract}
 \end{figure*}

\subsection{Outline}

In Section \ref{sec:Section2} we will give a compact overview of all methods that are needed in order to derive our proposed algorithm. In Section \ref{sec:Section3} we then present our ideas and the final method. The last section, Section \ref{sec:Section4} contains carefully conducted experiments with method, including important comparisons to other classical and novel algorithms that are known to work well for the recovery from incomplete Fourier measurements. Furthermore we will demonstrate the successful performance of our algorithm in the context of computed tomography.

\section{Sparse recovery, convex optimization, and sparsifying transforms}\label{sec:Section2}

In this section we present a short overview of current concepts and methods that are standard in the area of sparse recovery and are necessary to follow the rest of this work. For more details we refer the interested reader to the indicated literature.

\subsection{Compressed sensing and reweighted $\ell^1$-minimization}

The problem considered in compressed sensing can be explained by solving a system of underdetermined linear equations using prior information. Indeed, one is interested in solutions $u$ of the equation
\begin{align}
 Au = y, \label{eq:Axy}
\end{align}
where $y$ is a vector representing the acquired data, $A$ is a sensing matrix, and $u$ is the object of interest. The assumption that makes this problem in particular interesting is that $y$ should be of very
small dimension compared to $u$ that lives in a much higher dimensional space, i.e., $A \in \mathbb{C}^{m \times n}$ with $m \ll n$. Furthermore, the aforementioned prior information that enables us to solve such an underdetermined
system is \emph{sparsity}, i.e., although $u$ might be drawn from a much higher dimensional space only very few of its entries are nonzero. More precisely, a signal $u \in \mathbb{C}^n$ is called \emph{$s$-sparse} if
\begin{align*}
\| u \|_0 := \#\{ i \in \mathbb{N} \, : \, u(i) \neq 0, 1\leq i \leq n\}\leq s.
\end{align*}

A common approach to obtain sparse solutions of \eqref{eq:Axy} is to solve the following constrained convex optimization problem
\begin{align}
\label{min:standardl1}
  \min_u \| u\|_1 \quad \text{ subject to } \quad Au = y,
\end{align}
see for example \cite{CanRomTao,FouRau}.

One of the possibilities to improve the recovery model \label{min:l1minimization} is to strengthen the effect of sparsity in the minimization problem. This can be done, for instance, by using the idea of reweighted $\ell^1$ introduced by Cand\`{e}s et al. in \cite{CanWakBoy} which can be described as follows.

Suppose we are given measurements $y = Au \in \mathbb{C}^m$ of an $s$-sparse signal $$u = (u(1),\dots,u(n))^T \in \mathbb{C}^n$$ for a measurement matrix $A \in \mathbb{C}^{m \times n}$ for $m \ll n$.
When solving the minimization problem \eqref{min:standardl1} iteratively, one would ask for the following effect: large coefficients should be quickly identified and hence become ``cheaper'' in the minimization of the  objective function in \eqref{min:standardl1}, whereas very small coefficients
should be neglected in the minimization since they are most likely going to be zero in the true signal. More precisely, let $u_0$ be the true signal and define a  diagonal weighting matrix  $W$ by
\begin{align}
 W(i,i) = \begin{cases}
                                    \frac{1}{|u_{0}(i)|},&  u_{0}(i) \neq 0 \\
                                    \infty, & u_{0}(i) = 0.
                                   \end{cases} \label{weights}
\end{align}
Now, if the signal $u_0$ was $s$-sparse, then under some assumptions \cite{CanWakBoy} the \emph{weighted} $\ell^1$-minimization problem
\begin{align*}
 \min_u \| W u\|_1 \quad\text{ subject to } \quad Au = y,
\end{align*}
will find the exact solution. However, since $u_0$ is usually unknown such weights are practically infeasible. Therefore Cand\`{e}s et al. proposed adaptive weights in \cite{CanWakBoy} that change at each iteration depending on the previously computed solution $u^k$ which is an approximation to $u_0$. More precisely, the following sequence of minimization problems are then considered
\begin{align*}
u_{k+1} =  \argmin_u \| W_k u\|_1 \quad \text{ subject to } \quad Au = y,
\end{align*}
with a weighting matrix
\begin{align*}
 W_{k}(i,i) =   \frac{1}{|u_{k}(i)|+\varepsilon},
\end{align*}
where $\varepsilon > 0$ is a stability parameter and the initial weighting matrix $W^0$ is set to be the identity. In a series of numerical experiments it was shown in \cite{CanWakBoy} that such reweighting methods find sparse solution much faster with significantly reduced errors.

\subsection{Sparsifying multilevel transforms}

In the previous section, we discussed how the concept of sparsity is used to recover signals from possibly highly undersampled data. However, in many applications the signals are not directly sparse, but only  after the application of certain transforms. Such so-called \emph{sparsifying transforms} are often build upon systems that are equipped with a multiscale structure. Indeed, very recently a new direction of compressed sensing has been developed that is very much motivated by the sparsity structure of multiscale systems \cite{AdcHanPooRom}. Typical examples of such multiscale transforms are the \emph{wavelet} transform \cite{Dau}, the \emph{shearlet} transform \cite{GuoKutLab2006, KitKutLim, Lim}, and the \emph{curvelet} transform \cite{CanDon}. For the numerical results of this paper we will only consider the first two transforms depending on the particular signals that are to be recovered.  In particular, both, wavelet reconstructions as well as shearlet reconstructions from an incomplete amount of  Fourier measurements have been analyzed in the literature, for instance, in \cite{AdcHanKutMa} and \cite{Ma2015a, KutLimFourier}, respectively. However, the methodology of algorithm applies to any other multiscale transform other than wavelets and shearlets.

The multiscale structure of a wavelet basis comes from the use of dyadic scaling matrices of the form
\begin{align*}
    D_j = \begin{pmatrix} 2^j & 0 \\ 0 & 2^{j} \end{pmatrix}, \quad j =0,1,\ldots.
\end{align*}
Using these matrices together with  simple translations one can eventually obtain \emph{orthonormal bases} for $L^2(\mathbb{R}^2)$, the space of square integrable functions, of the form
\begin{align*}
\{ \phi(\cdot - m ) \, : \, m \in \mathbb{Z}^2 \} &\cup\{ \psi^1(D_j \cdot - m ) \, : \, j \geq 0, m \in \mathbb{Z}^2 \}\\
&\cup\{ \psi^2(D_j \cdot - m ) \, : \, j \geq 0, m \in \mathbb{Z}^2 \}\\
&\cup\{ \psi^3(D_j \cdot - m ) \, : \, j \geq 0, m \in \mathbb{Z}^2 \},
\end{align*}
where $\phi, \psi^1, \psi^2, \psi^3 \in L^2(\mathbb{R}^2)$ with certain regularity properties, see \cite{Dau} for more details.

The multiscale structure of shearlets are obtained by the use of parabolic scaling matrices of the form
\begin{align*}
A_{2^j} = \begin{pmatrix} 2^j & 0 \\ 0 & 2^{j/2} \end{pmatrix}, \quad \widetilde{A}_{2^j} = \begin{pmatrix} 2^{j/2} & 0 \\ 0 & 2^{j} \end{pmatrix}, \quad j \in \mathbb{N}.
\end{align*}
In addition to the parabolic scaling matrix, a shearlet system is equipped with a directional component that can be obtained by using shear matrices
\begin{align*}
S_k = \begin{pmatrix} 1 & k \\ 0 & 1 \end{pmatrix}, \quad k \in \mathbb{Z}.
\end{align*}
The so-called \emph{cone-adapted shearlet system} is then defined as follows.
\begin{deff}
Let $\varphi, \psi, \widetilde{\psi} \in L^2(\mathbb{R}^2)$. Then we call $\Phi(\phi,c) \cup \Psi(\phi, c) \cup \widetilde{\Psi}(\widetilde{\psi},c)$ a \emph{cone-adapted shearlet system}, where
\begin{align*}
\Phi(\phi,c) &= \{ \phi_m \, : \, m \in \mathbb{Z}^2\}, \\
\Psi(\psi,c) &= \{ \psi_{j,k,m} \, : \, j \geq 0, |k| \leq 2^{j/2}, m \in \mathbb{Z}^2\}, \\
\widetilde{\Psi}(\widetilde{\psi},c) &= \{ \widetilde{\psi}_{j,k,m} \, : \, j \geq 0, |k| \leq 2^{j/2}, m \in \mathbb{Z}^2\},
\end{align*}
and
\begin{align*}
\phi_m &= \phi(\cdot - c_1m), \\
\psi_{j,k,m} &= 2^{3/4j} \psi(S_k A_{2^j} \cdot - cm), \\
 \widetilde{\psi}_{j,k,m} &= 2^{3/4j} \widetilde{\psi}(S_k^T \widetilde{A}_{2^j} \cdot - \widetilde{c}m),
\end{align*}
with $c = (c_1,c_2)^T \in \mathbb{R}^2_+, \widetilde{c} = (c_2,c_1)$ and the multiplication of $c$ and $\widetilde{c}$ with $m$ to be understood componentwise.
\end{deff}

The shearlet coefficients $\{ (\langle u , \phi_m\rangle )_m\} \cup \{ (\langle u, \psi_{j,k,m} \rangle)_{j,k,m} \} \cup \{ (\langle u, \widetilde{\psi}_{j,k,m} \rangle)_{j,k,m} \}$  computed from a shearlet transform are known to have a fast decay which ensures within the model of so-called \emph{cartoon-like images} \cite{DonSparse} an optimal sparse approximation rate \cite{KutLim}. Wavelets on the other hand, do not fulfill this optimal approximation rate. In Figure \ref{Sparsity2} we depicted some of the shearlet coefficients of each scale for an MRI test image.

\begin{figure}[H]
\includegraphics[width=0.24\textwidth]{./Im0.png}
\includegraphics[width=0.24\textwidth]{./Im1c.png}
\includegraphics[width=0.24\textwidth]{./Im2c.png}
\includegraphics[width=0.24\textwidth]{./Im3c.png}
\caption{Shearlet coefficients at different scales of the GLPU Brain phantom \cite{GLPU}. For better visual difference the contrast at scale 1,2, and 3 has been changed.
\textbf{First}: Shearlet coefficients at scale $j = 0$. \textbf{Second}: Shearlet coefficients at scale $j =1$. \textbf{Third}: Shearlet coefficients at scale $j =2$. \textbf{Fourth}: Shearlet coefficients at scale $j =3$. } \label{Sparsity2}
\end{figure}

\subsection{Example: Multilevel iterative reweighting for inpainting}\label{sec:inpainting}

As already outlined above, reweighted $\ell^1$-minimization can significantly improve the reconstruction of certain sparse signals. The reader might wonder at this point, why an adaptive strategy for choosing the weights is beneficial if multilevel transforms are used. We demonstrate the basic idea of multilevel iterative reweighting for a simple \emph{inpainting} or \emph{image restoration} example using shearlets. Suppose $A$ denotes a masking operator, $y = Au$ is the masked version of the original image $u$, cf. Figure \ref{fig:inpainting}, and $\Psi, \Psi^{-1}$ denote the forward and backward shearlet transform, respectively.  We compute two reconstructions using the iterative hard thresholding  Algorithm \ref{algo:IHT} for two different thresholding strategies. Note that thresholding algorithms are classical and well known in the literature for such image restoration tasks. It has also been used in combination with shearlets in
\cite{KutLimRei}.

In the following, let $\delta = (\delta(1), \ldots, \delta(N))^T \in \mathbb{R}^N$, where $N$ is the length of the transform coefficient vector $\Psi u$. Further, let $T_\delta: \mathbb{R}^N \longrightarrow \mathbb{R}^N$ be the hard thresholding operator applied entrywise, that is
\begin{align*}
T_\delta(c) = h, \quad \text{where} \quad h(k) = \begin{cases} c(k) & |c(k)| > \delta(k), \\ 0 & |c(k)| \leq \delta(k) \end{cases}
\end{align*}
for $k = 1, \ldots, N$.
\begin{algorithm}
\caption{Iterative hard thresholding}\label{algo:IHT}
\begin{algorithmic}
\STATE{\underline{Input}:
Measurements $y$, sampling operator $A$, sparsifying transform $\Psi$, initial values for $\delta$ and $\lambda$, thresholding strategy \eqref{eq:Strat1}  or \eqref{eq:Strat2} denoted by $f$, factor $\sigma <1$.\\
\underline{Initialization}:\\
$u_{\text{rec}},u_{\text{res}} \gets 0$;
}
\FOR{$i =1:N$}
\STATE{$u_{\text{res}} = A(y - u_{\text{rec}})$;\\
$u_{\text{rec}} = {\Psi}^{-1}(T_\delta(\Psi(u_{\text{res}} + u_{\text{rec}})))$; \\
$\delta = \sigma \cdot f(u_{\text{res}} + u_{\text{rec}}, \lambda)$;
}
\ENDFOR
\end{algorithmic}
\end{algorithm}

The shearlet coefficients $\{ (\langle u, \phi_m\rangle )_m\} \cup \{ (\langle u, \psi_{j,k,m} \rangle)_{j,k,m} \} \cup \{ (\langle u, \widetilde{\psi}_{j,k,m} \rangle)_{j,k,m} \}$  can be divided into corresponding levels $j$. For the first cone this looks as follows
\begin{align*}
 (\langle u, \psi_{j,k,m} \rangle)_{j,k,m} = ((\langle u, \psi_{j,k,m} \rangle)_{(j,k,m) \in I_j})_j
\end{align*}
where $I_j$ contains all indices at scale $j$. Two different thresholding strategies involving reweighted $\ell^1$ are now for example
\begin{align}
 f_1(u,\lambda) &= \left (\frac{\lambda}{|\langle u, \psi_{j,k,m} \rangle| + \varepsilon}\right )_{j,k,m}  \label{eq:Strat1}
\end{align}
and
\begin{align}
 f_2(u,\mu) &= \left (\frac{\mu\max\{(|\langle u, \psi_{j,k,m}\rangle|)_{(j,k,m) \in I_j}\}}{|\langle u, \psi_{j,k,m} \rangle| + \varepsilon}\right )_{j,k,m},   \label{eq:Strat2}
\end{align}
where $\lambda, \mu$ and $\varepsilon>0$ are fixed parameters. Notice that the first strategy \eqref{eq:Strat1} is a direct application of reweighted $\ell^1$ to the thresholding based recovery algorithm Algorithm \ref{algo:IHT}, whereas strategy \eqref{eq:Strat2} involves an additional compensation factor to adapt the idea of reweighting to the multilevel structure. Such adaptations are useful as the transform coefficients of a multilevel system have a natural decrease due to scaling. Independently of the usage of reweighting, such types of compensation factors should be used in thresholding based algorithms. This is, for instance, used in \cite{HauMa} and we will discuss other choices in Section \ref{sec:Section3}.

In Figure \ref{fig:inpainting}  one can observe that by introducing level adapted weights, the information across all scales will be treated equally important. Without such weights, the fine detail coefficients are underrated in terms of their importance which might lead to crucial quality loss of detail information.

\begin{figure*}
 \includegraphics[width=0.247\textwidth]{./masked.png}
\includegraphics[width=0.37\textwidth]{./IRzoom2.png}
\includegraphics[width=0.37\textwidth]{./MIRzoom2.png}
\caption{\textbf{Left:} Sampled data of an image. \textbf{Middle:} Reconstructed image using iterative reweighting (Strategy \eqref{eq:Strat1}). \textbf{Right:} Reconstructed image using multilevel iterative reweighting (Strategy \eqref{eq:Strat2}). Note that the fine details are much better restored in the second case.} \label{fig:inpainting}
\end{figure*}

Algorithm \ref{algo:IHT} is simple and effective, however, it is not very sophisticated and also not sufficient for our purposes as in this work we also want to consider an additional regularizer. Hence, we continue with an introduction to \emph{split Bregman}.

\subsection{Split Bregman}
\label{SB}
Split Bregman is an effective algorithm to solve constrained optimization problems by introducing a split variable and solving the resulting decoupled problem with \emph{Bregman Iterations}. It was proposed in 2009 by Goldstein and Osher and became a popular method since then \cite{SB,SB2,Wu,Plonka,Steidl}. However, as we have already mentioned above, the algorithm was already invented under the name ADMM in \cite{GabMer,EckBer}; see also \cite{Setzer,Esser}.

Even though split Bregman can handle general convex regularizers it is sufficient for our purposes to focus on $\ell^1$-regularized problems in the analysis formulation. We briefly follow the steps in \cite{Burger,SB,Yin} to derive the basic form of the algorithm.

Consider the basis pursuit problem in analysis formulation
\begin{align}
    \label{BPinAna}
    \min_u  \norm{W \Psi u}_1 \quad \mbox{ subject to } \quad  Au = y,
\end{align}
for a possibly redundant dictionary $\Psi \in \mathbb{R}^{N}\times n, N, n \in \mathbb{N}$, a measurement matrix $A \in \mathbb{C}^{m\times n}$, and a diagonal weighting matrix with entries
$
W(l,l)$ for $l=1,\dots,N$.
Then instead of using a continuation method for enforcing the constraint, i.e., taking $\beta \to \infty$ in
\begin{align*}
    u= \argmin_u \norm{W\Psi u}_1 + \frac{\beta}{2} \norm{y -Au}_2^2,
\end{align*}
problem \eqref{BPinAna} is transformed into a sequence of unconstrained problems using \emph{Bregman iterations}
\begin{align}
\label{BPinAnaBreg}
\begin{cases}
    u_{k+1} = \argmin \limits_u \norm{W \Psi u}_1 + \frac{\beta}{2} \norm{Au - y_k}_2^2,\\
    y_{k+1} = y_k + y - Au_{k+1},
\end{cases}
\end{align}
for a suitable $\beta >0$. For keeping the presentation in this theoretical study as concise as possible we focus  on  noiseless measurements and remark that if noise is present, model \eqref{BPinAna} and the  iterations \eqref{BPinAnaBreg} have to be adapted carefully, cf. \cite{Yin, Burger}. A performance of the presented framework in this situation is demonstrated in an upcoming work \cite{MaMaeFunSchKutSchKol}.
We continue by introducing a split variable $d = \Psi u$ for the $\ell^1$-part of the minimization problem in \eqref{BPinAnaBreg} and executing an additional Bregman iteration step to obtain
\begin{alignat*}{1}
\begin{cases}
(u_{k+1},d_{k+1}) & = \argmin \limits_{u,d}  \norm{Wd}_1   + \frac{\beta}{2} \norm{Au - y_k}_2^2  \\
& \qquad  \qquad + \frac{\mu}{2} \norm{d - \Psi u - b_k}_2^2, \\
b_{k+1} & =   b_k + \Psi u_{k+1} - d_{k+1},\\
y_{k+1} & =  y_k + y - Au_{k+1}.
\end{cases}
\end{alignat*}
To solve the $(u,d)$-minimization problem one or multiple \emph{nonlinear block Gauss-Seidel} iterations are used, which alternate between minimizing with respect to $u$ and $d$. This yields the \emph{split Bregman Algorithm}
\begin{algorithm}
\caption{Split Bregman algorithm}\label{SBAlgo}
\begin{algorithmic}

\FOR{$i =1:N$}
\STATE{$u_{k+1} =\argmin_{u}  \frac{\beta}{2} \norm{Au - y_k}_2^2 + \frac{\mu}{2} \norm{d^{\text{cur}} - \Psi u - b_k}_2^2$\\
$d_{k+1}  = \argmin_d \norm{Wd}_1 + \frac{\mu}{2} \norm{d - \Psi u^{\text{cur}} - b_k}_2^2,$
}
\ENDFOR
\STATE{$b_{k+1}  =   b_k + \Psi u_{k+1} - d_{k+1}$,\\
$y_{k+1}  =  y_k + y - Au_{k+1},$}
\end{algorithmic}
\end{algorithm}

where $x^{\text{cur}}$ denotes the latest available stage of the variable $x \in \{ d, u \}$. Note that the solution of the $d$-subproblem  is explicitly given by \emph{soft-thresholding}
\[
d_{k+1}(l) = \text{shrink}\left(\left(\Psi u^{\text{cur}}\right)(l) + b_{k}(l), \frac{W(l,l)}{\mu} \right),
\]
for $l=1,\dots,N$ where
\[
\text{shrink}\left(z,\lambda\right) = \begin{cases} \max\left(\norm{z} - \lambda,0\right) \frac{z}{\norm{z}}, & z \neq 0, \\
0, & z = 0. \end{cases}
\]
In \cite{SB} it was furthermore observed that the minimization with respect to $u$ in \eqref{SBAlgo} does not have to be solved to full precision and in many applications only few steps of an iterative method are sufficient.

\subsection{TV and TGV}\label{TGV}

\emph{Total Variation (TV)} based methods were initially proposed by Rudin, Osher, and Fatemi in 1992  for image denoising \cite{RudOshFat} and are now widely used for image reconstruction and compressed sensing, see for example \cite{RecPF,NeedellWard}. TV is based on the assumption that the reconstructed image is piecewise constant and therefore gradient sparse. This results in preserving sharp edges. But for realistic images, which are usually not piecewise constant, this can lead to severe \emph{oil painting artifacts} or  \emph{staircasing effects} leading to unnatural looking reconstructed images.

\emph{Total Generalized Variation} (TGV) is a generalization of TV and has been proposed to improve on these issues by involving higher order derivatives \cite{BreKunPoc}. We will now briefly give a definition of the second order TGV regularizer in $\mathbb{R}^2$ together with some basic facts. Its general derivation and more details on this subject can be found in \cite{TGVinverse,Bredies2015,BreKunPoc,Holler,Holler2}. Please note that the discretized TGV-model, which is solved numerically, was already introduced in \cite{SS08}.

The so-called \emph{pre-dual} formulation of second order TGV is given by
\begin{align}
\TGV_{\alpha}^2(u) =  \sup  \bigg\{ & \int_{\Omega}  u \divergenz^2 v \, dx : v \in C_c^2(\Omega,\mathcal{S}^{2\times 2}), \nonumber \\
& \norm{v}_{\infty} \leq \alpha_0, \norm{\divergenz v}_{\infty} \leq \alpha_1 \bigg\},
\end{align}
for $\alpha = (\alpha_0,\alpha_1) \in \mathbb{R}_+^2$, $\Omega \subseteq \mathbb{R}^2$ a bounded Lipschitz domain, $\mathcal{S}^{2 \times 2}$ the space of symmetric $2\times 2$ matrices, and $u \in L^1(\Omega,\mathbb{C})$. Thereby the divergences are defined as
\[
(\divergenz v)_i =\sum_{j=1}^2  \frac{\partial v_{ij}}{\partial x_j}, \qquad i=1,2,
\]
and
\[
\divergenz^2 v = \sum_{i,j=1}^2 \frac{\partial^2 w_{ij}}{\partial x_i \partial x_j},
\]
together with the norms
\[
\norm{v}_{\infty} = \sup_{l\in\Omega} \left(\sum_{i,j=1}^2 |v_{ij}(l)|^2 \right)^{1/2},
\]
and
\[
\norm{\divergenz v}_{\infty} = \sup_{l\in\Omega} \left(\sum_{i=1}^2 |(\divergenz v)_i(l)|^2 \right)^{1/2}.
\]
Under certain conditions, an equivalent and more convenient form of $\TGV_{\alpha}^2$ is given by the \emph{minimum representation} as
\begin{align}
\TGV_{\alpha}^2(u) = \inf_{v \in \text{BD}(\Omega,\mathbb{C}^2)} \alpha_1 \norm{\nabla u -v}_1  + \alpha_0 \norm{\mathcal{E}(v)}_1,
\end{align}
where $\text{BD}(\Omega, \mathbb{C}^2)$ is the space of symmetric tensor fields of bounded deformation and $\mathcal{E}$ the symmetrized derivative defined as
\[
\mathcal{E}(v) = \left(\begin{matrix}\partial_x v_1 &  \frac{1}{2}(\partial_y v_1 +  \partial_x v_2) \\ \frac{1}{2}(\partial_y v_1 + \partial_x v_2) & \partial_y v_2  \end{matrix}\right).
\]
In this form $\TGV_{\alpha}^2$ can be interpreted as balancing the first and second derivatives of $u$ controlled by the ratio of $\alpha_0$ and $\alpha_1$. In \cite{BreKunPoc} and \cite{KnoBrePocRud} it was observed that the use of TGV as a regularizer indeed leads to reconstructed images with sharp edges but without the staircaising effects of TV. Being furthermore convex and lower semi-continuous with respect to $L^1$-convergence, this makes TGV a numerical feasible and therefore a suitable alternative for TV \cite{TGVinverse}.

\section{Proposed multilevel based reweighting algorithm with TGV}\label{sec:Section3}

\subsection{Model and discretization}
In this section we aimto develop an algorithm for solving the multilevel reweighted $\ell^1$-problem. In order to do so, the split Bregman approach introduced in Section \ref{SB} shall be equipped with an appropriate iteratively reweighted soft-thresholding procedure.  For a further reduction of artifacts and improving the reconstruction of piecewise constant as well as smooth regions we not only use the regularizer belonging to the reweighted multilevel decomposition but also TGV as a second regularization term.

Let $A$ be a measurement operator, let $y$ the measurements of our signal of interest $u$, and let $\sigma>0$ be a fidelity parameter. The recovery problem can then be stated as
\begin{align*}
    \min_u \sum_{j=1}^{\infty} \lambda_j \norm{W_j \Psi_j u}_1 + \TGV_{\alpha}^2(u)  \quad \mbox{ subject to } \quad Au = y,
\end{align*}
where $\Psi_j$ corresponds to the $j$-th subband of the multilevel transform $\Psi $, $\lambda_j$ are regularization parameters accounting for the multilevel structure of $\Psi $ and $W_j$ are diagonal matrices containing the weights. For the sake of clearness we have assumed that there is only one subband per level otherwise an additional index has to be attached to $\Psi_j$ to specify the current subband.
Note that after we have established a basic split Bregman framework for solving the minimization problem we will aim to update $\lambda_j$ and $W_j$ iteratively.
Using the characterization of $\TGV_{\alpha}^2$ presented in Section \ref{TGV} the objective can be rewritten as
\begin{align}
\label{Objective}
\min_{u,v} \sum_{j=1}^{\infty} \lambda_j \norm{W_j \Psi_j u}_1 + \alpha_1 \norm{\nabla u - v}_1 + \alpha_0 \norm{\mathcal{E}(v)}_1.
\end{align}
For the \emph{discretization}, let $u\in\mathbb{C}^{n^2}$ be the vectorized finite-dimensional image of interest which is for simplicity assumed to be of square size. Let $A\in \mathbb{C}^{m\times n^2}$ be the finite dimensional measurement matrix and $y\in \mathbb{C}^m$ the observed data. Furthermore, let $\nabla^f$ and $\nabla^b$ denote a discrete gradient operator with periodic boundary conditions using forward and respectively backward differences. Following \cite{BreKunPoc,Holler2} we  approximate the derivatives in \eqref{Objective} by
\[
\nabla u \approx \nabla^f u = \left(\begin{matrix} \nabla_x^f u\\ \nabla_y^f u \end{matrix}\right)
\]
and
\[
\mathcal{E}(v) \approx \mathcal{E}^b v =  \left(\begin{matrix} \nabla_x^b v_x & \frac{1}{2}(\nabla_y^b v_x + \nabla_x^b v_y )\\ \frac{1}{2}(\nabla_y^b v_x + \nabla_x^b v_y ) & \nabla_y^b v_y \end{matrix} \right).
\]
A finite dimensional approximation of \eqref{Objective} is then given by
\begin{align}
\label{discrObj}
\min_{u,v} \sum_{j=1}^{J} \lambda_j \norm{W_j \Psi_j u}_1 + \alpha_1 \norm{\nabla^f u - v}_1 + \alpha_0 \norm{\mathcal{E}^bv}_1,
\end{align}
where $J$ is some fixed a priori chosen maximum scale and $\Psi$ is the discrete transform acting on the vectors.
Let us furthermore introduce the notation
\begin{align}
\begin{split}
	\Psi u  = \left(\Psi_j  u \right)_{j=0,\dots,J}  = (\langle \psi_{j, l}, u \rangle)_{j=0,\ldots, J, l =1, \ldots, N_j} \label{eq:CoeffNj}
	\end{split}
\end{align}
for dividing the analysis coefficients into $J$ subbands, each consisting of $N_j \in \mathbb{N}$ elements.
For a documentation of the discrete transforms, we refer the interested reader to Chapter 8 in \cite{Walnut} for wavelets and \cite{KutLimRei} for shearlets.
Note that the $\ell^1$-norm in the second summand of \eqref{discrObj} is thereby defined as
\[
\norm{v}_1 = \sum_{l=1}^{n^2} \left(|v_x(l)|^2 +|v_y(l)|^2 \right)^{1/2},
\]
and for the third summand as
\[
\norm{e}_1 = \sum_{l=1}^{n^2} \norm{e(l)}_F = \sum_{l=1}^{n^2} \norm{\left(\begin{matrix}
 e(l)_1 & e(l)_2 \\
 e(l)_2 & e(l)_3
\end{matrix} \right)}_F,
\]
where $\norm{\cdot}_F$ is the Frobenius norm of a $2 \times 2$ matrix.
\subsection{Split Bregman framework}
The proposed constrained optimization problem can be casted into the form given in  \eqref{BPinAna} by introducing the variable ${{\bf u} = (u,v)^T}$, together with the matrix
\[
{\bf \Psi } = \left(\begin{matrix} \Psi  & 0 \\  \nabla^f & - I \\ 0 &  \mathcal{E}^b  \end{matrix} \right).
\]
In order to  come up with the explicit form of the resulting split Bregman algorithm as given in Section \ref{SB}, let us  split  as follows:
\[
\left( \begin{matrix}
w\\d \\t
\end{matrix}
\right)
= \left(
\begin{matrix}
 \Psi u \\ \nabla^f  u -v  \\ \mathcal{E}^bv
\end{matrix}  \right).
\]
The $(u,v)$-subproblem of Algorithm \ref{SBAlgo} is then given by
\begin{align}
\begin{split}
\label{uProb}
 &(u_{k+1},v_{k+1}) \\
 &=  \argmin_{u,v}  \frac{\beta}{2} \norm{Au - y_k}_2^2   + \frac{\mu_1}{2} \norm{w^{\cur}- \Psi u - b^w_k}_2^2 \\
   &+ \frac{\mu_2}{2}\norm{d^{\cur} - (\nabla^f u - v) - b^d_k }_2^2   + \frac{\mu_3}{2}\norm{t^{\cur} - \mathcal{E}^b v - b^t_k}_2^2.
\end{split}
\end{align}
We furthermore obtain the subproblems
\begin{align}
\label{subw}
w^j_{k+1} = \argmin_{ w^j} \lambda_j  \norm{W_j  w^j}_1
+ \frac{\mu_1}{2}  \norm{w^j - \Psi_j u^{\cur} - b^{w,j}_k}_2^2,
\end{align}
for each subband $j=1,\dots,J$, as well as
\begin{align}
\label{subd}
d_{k+1} = \argmin_d \alpha_1 \norm{d}_1 + \frac{\mu_2}{2} \norm{d - (\nabla^f u^{\cur} - v^{\cur}) - b^{d}_k}_2^2,
\end{align}
and
\begin{align}
\label{subt}
t_{k+1} = \argmin_t \alpha_0 \norm{t}_1 + \frac{\mu_3}{2} \norm{t - \mathcal{E}^b v^{\cur} - b^t_k}_2^2.
\end{align}
Note that the regularization parameters $\lambda_j, \alpha_0,$ and $\alpha_1$ have thereby been subsumed into a weighting matrix $W$. Also we are allowing some more flexibility by incorporating different values for $\mu_i$ for $i=1,2,3$.
Furthermore we obtain the following \emph{Bregman updates:}
\begin{align*}
\begin{cases}
b^w_{k+1} & = b^w_k + \Psi u_{k+1} - w_{k+1}. \\
b^d_{k+1} & = b^d_k + (\nabla^f u_{k+1} -v_{k+1}) -d_{k+1},\\
b^t_{k+1} & = b^t_k + \mathcal{E}^b v_{k+1} - t_{k+1},
\end{cases}
\end{align*}
as well as
\begin{align*}
    y_{k+1}  = y_k + y - Au_{k+1}.
\end{align*}

\subsection{Solutions of the subproblems}
The solution of the subproblem \eqref{uProb} can be obtained by setting the first derivatives with respect to $u$, $v_x$, and $v_y$ to zero. We then obtain the linear system
\begin{align}
\label{sys}
   \left(\begin{matrix}
    b_1 & b_4^* & b_5^* \\
    b_4 & b_2 & b_6^* \\
    b_5 & b_6 & b_3
   \end{matrix} \right)
   \left(
   \begin{matrix}
    u \\ v_x \\ v_y
   \end{matrix}
   \right)
   = \left(
   \begin{matrix}
    R_1 \\ R_2 \\ R_3
   \end{matrix}
   \right),
\end{align}
where $b_i$ are $n^2\times n^2$ block matrices defined as
\begin{align*}
     &  b_1 = \beta A^*A + \mu_1 \Psi^* \Psi + \mu_2 (\nabla^f)^* \nabla^f, \\
    & b_2 = \mu_3 (\nabla_x^b)^* \nabla_x^b + \frac{\mu_3}{2}   (\nabla_y^b)^* \nabla_y^b + \mu_2 I,\\
    & b_3 = \mu_3 (\nabla_y^b)^* \nabla_y^b + \frac{\mu_3}{2} (\nabla_x^b)^* \nabla_x^b
     + \mu_2 I, \\
     & b_4 = - \mu_2  \nabla_x^f, \\
     & b_5 = -\mu_2 \nabla_y^f,\\
     & b_6 = \frac{\mu_3}{2} (\nabla_x^b)^* \nabla_y^b,
\end{align*}
and the components of the right hand side are given by
\begin{alignat*}{2}
     R_1 &= \beta A^*(y + y_k)  + \mu_1 \Psi^*(w^{\cur} - b_k^w) \\
     &\qquad + \mu_2 (\nabla^f)^* (d^{\cur} -b_k^d), \\
     R_2 &=   \mu_2 ( b_{k,x}^d - d_x^{\cur}) + \mu_3 \big[  (\nabla_x^b)^* (t_1^{\cur}-b^t_{k,1})\\
          &\qquad +(\nabla_y^b)^* (t_2^{\cur}-b^t_{k,2})\big], \\
  R_3 &=   \mu_2 ( b_{k,y}^d - d_y^{\cur})   + \mu_3 \big[  (\nabla_x^b)^* (t_2^{\cur}-b^t_{k,2})\\
 & \qquad+(\nabla_y^b)^* (t_3^{\cur}-b^t_{k,3})\big].
\end{alignat*}
Similar to \cite{SB}, it was observed in \cite{GuoQinYin}, that in many cases the linear system in \eqref{sys} can be efficiently solved by using the $2$D-Fourier transform $\mathcal{F} \in \mathbb{C}^{n^2 \times n^2}$. Note that $\nabla^f$ and $\nabla^b$ are circulant, since they correspond to periodic boundary conditions. Therefore
\[
\mathcal{F} \nabla^* \nabla \mathcal{F}^*
\]
is a diagonal matrix. For a tight frame $\Psi$ we have
\begin{align}
\Psi^* \Psi = a I,\label{eq:TightFrame}
\end{align}
where $a\in\mathbb{R}$ is the frame bound of $\Psi$. This is for example the case in \cite{GuoQinYin}, where the Fast Finite Shearlet Transform (\cite{Hae13,Hae14}) was used, which forms a Parseval frame for $\mathbb{R}^{n\times n}$, i.e., Equation \eqref{eq:TightFrame} with $a = 1$. Note that also for the non-tight shearlet system of ShearLab \cite{KutLimRei}  the matrix $\mathcal{F}\Psi^*\Psi\mathcal{F}^*$ is diagonal and can be explicitly computed.

In the case of  subsampled Fourier measurements, the measurement matrix can be written as
\begin{align*}
A = P \mathcal{F},
\end{align*}
where $P \in \left\{0,1\right\}^{m\times n^2}$ is selecting or discarding measurements. In this case
\[
A^* A = \mathcal{F}^* P \mathcal{F}
\]
is naturally diagonalized by the 2D Fourier transform.

If all blocks $b_i$ for $i=1,\dots,6$ can be diagonalized in this way, the authors of \cite{GuoQinYin} proposed to multiply with a preconditioner matrix from the left to obtain the system
\begin{align}
\label{diagsys}
\left( \begin{matrix}
\widehat{b_1} & \widehat{b_4}^* &\widehat{b_5}^*\\
\widehat{b_4} & \widehat{b_2} & \widehat{b_6}^* \\
\widehat{b_5} & \widehat{b_6} & \widehat{b_3} \end{matrix} \right)
\left(
\begin{matrix}
\mathcal{F}u \\ \mathcal{F} v_x \\ \mathcal{F} v_y \end{matrix} \right)
 = \left( \begin{matrix}\mathcal{F}R_1 \\ \mathcal{F}R_2 \\ \mathcal{F} R_3 \end{matrix} \right),
\end{align}
where each $\widehat{b_j} = \mathcal{F} b_j\mathcal{F}^*$ is a $n^2\times n^2$ diagonal matrix. A closed form solution can then be obtained by applying Cramer's rule.

In other applications, such as inpainting or reconstruction from the subsampled cosine transform, the measurement matrix can be written as
\begin{align}
A = P \mathcal{T}, \label{eq:AUnit}
\end{align}
for a unitary matrix $\mathcal{T} \in \mathbb{C}^{n^2 \times n^2}$. In the ADMM model of \cite{GuoQinYin} it is proposed to include an additional split to deal with the fact that $A^* A$ cannot be diagonalized under $\mathcal{F}$.

\noindent However, note that a representation of the form \eqref{eq:AUnit} is not possible for instance in
\begin{enumerate}
\item \emph{Partial Parallel Imaging (PPI)} \cite{PPI,PPI2} in MRI,  where the image is to be reconstructed by using subsampled Fourier data from parallel scans of multiple coils.
\item \emph{Computed tomography}, where the image has to be reconstructed from \emph{fan beam} projections, see Section \ref{sec:Radon} for more details.
\end{enumerate}
In such cases we therefore propose to solve the linear system of equation \eqref{sys} only approximately. As it was observed in \cite{SB}, usually only a few steps of an iterative solver are sufficient for the convergence of the resulting algorithm. A proof of this fact can be found in \cite{Nien}. In order to save memory and computation time, we still want to use  that all $d_j$ are diagonal under $\mathcal{F}$ except for $\widehat{d_1}$.  Therefore we multiply with the same preconditioner matrix as above and obtain a similar system as in equation \eqref{diagsys}, with the difference that $\widehat{d_1}$ is not diagonal anymore. With a few explicitly given steps the block matrix can be brought to \emph{lower triangular form}. In this way
we only need to apply an iterative solver to one $n^2 \times n^2$ system involving $A^* A$ instead of solving the entire block system.
In both of the above listed cases, we are using a projected conjugate gradient method together with a warm start (obtained through the previous iterations of the  split Bregman algorithm), so that only very few steps are necessary for sufficient precision; see also Section \ref{sec:Radon} and \cite{MaMaeFunSchKutSchKol}.

Let us now briefly discuss the solution of the other subproblems:
As described in Section \ref{SB} a closed-form solution of \eqref{subw} is given by
\begin{align}
\label{shrinkw}
w_{k+1}^j(l) = \text{shrink}\left(\left(\Psi_j u^{\cur}\right)(l) + b_k^{w,j}(l), \frac{\lambda_j  {W_j}(l)}{\mu_1} \right),
\end{align}
for $l=1,\dots,N_j-N_{j-1}+1$.
For equation \eqref{subd} we obtain
\begin{align*}
d_{k+1}(l) = \text{shrink}_2\left( \nabla^f u^{\cur}(l) - v^{\cur}(l) + b_k^d(l),\frac{\alpha_1}{\mu_2} \right),
\end{align*}
for $l=1,\dots,n^2$ and the shrinkage rule
\begin{align*}
\text{shrink}_2\left(x,\lambda\right) = \begin{cases}
\max(\norm{(x}_2 - \lambda),0) \frac{x}{\norm{x}_2}, & x \neq 0,\\
0, & x= 0.
\end{cases}
\end{align*}
Similarly, the solution of \eqref{subt} is given by
\begin{align*}
t_{k+1}(l) = \text{shrink}_F\left(\left(\mathcal{E}^b v^{\cur}\right)(l) + b^t_k(l), \frac{\alpha_0}{\mu_3} \right),
\end{align*}

for $l=1,\dots,n^2$ and
\begin{align*}
\text{shrink}_F\left(x,\lambda\right) = \begin{cases}
\max(\norm{(x}_F - \lambda),0) \frac{x}{\norm{x}_F}, & x\neq 0,\\
0, & x = 0.
\end{cases}
\end{align*}

\subsection{Combining reweighted $\ell^1$ with multiscale transforms}\label{Sec:IIIb)}

As we already outlined in Section \ref{sec:Section2}, the idea of \emph{reweighted $\ell^1$} as proposed by Cand\`{e}s et al. in \cite{CanWakBoy} is to improve the $\ell^1$-norm as a sparsity regularizer by an iterative reweighting of the nonzero coefficients. In contrast to $\ell^1$-minimization, $\ell^0$ only counts the number of nonzero coefficients and thus spares larger coefficients more than $\ell^1$ does.
Hence the idea of reweighted $\ell^1$ is that since larger coefficients of an iterative solution are more likely to be nonzero in the true signal, they should be multiplied with a smaller weight during the optimization process.
Put differently, the guiding principle of reweighted $\ell^1$ is that small coefficients of an iterative solution are likely going to be zero in the true signal. However, this principle is not valid for multiscale sparse signals, i.e., signals that can be sparsely represented using a multiscale transform. The magnitudes of multiscale coefficients naturally decrease with increasing scales, but the high scale nonzero coefficients of an iterative solution are not necessarily less important or more likely zero in the actual signal, if compared to low scale coefficients which are intrinsically larger. In the following section we are aiming to compensate for this misfit by including additional weighting parameters for each level in the transformation.

Suppose $u \in \mathbb{R}^{n^2}$ is the true signal. Let us recall the notation
\begin{align*}
\begin{split}
	\Psi u  = \left(\Psi_j  u \right)_{j=0,\dots,J}  = (\langle \psi_{j, l}, u \rangle)_{j=0,\ldots, J, l =1, \ldots, N_j},
	\end{split}
\end{align*}
for dividing the analysis coefficients into $J$ subbands, each consisting of $N_j$ elements. In \cite{AhmSch} a multi dictionary reweighting algorithm (Co-L1) was proposed  which iteratively updates $\lambda_j^{k}$ in
\begin{align}
\label{prob}
u_{k+1} = \argmin_{u} \sum_{j=0}^{J} \lambda_j^{k} \norm{W_j \Psi_j u}_1 \quad \text{ subject to } \quad Au = y,
\end{align}
by setting
\begin{align}
\label{Schniter}
\lambda_j^{k} = \frac{N_j}{\varepsilon + \norm{\Psi_j u^{k}}_1},
\end{align}
and $W_j = I$ for all iterations of solving \eqref{prob}. It was shown therein that the resulting algorithm can be interpreted as applying a Majorization-Minimization algorithm to the unconstrained formulation of \eqref{prob} with the regularizer
\begin{align*}
    \sum_{j=0}^J N_j \log\left(\varepsilon + \norm{\Psi_j u}_1 \right).
\end{align*}
This update rule was proposed in \cite{AhmSch} for a composition of multiple different dictionaries instead of just one multiscale dictionary divided into its subbands. In the latter case it is less likely to expect that $\log \left(\varepsilon + \norm{\Psi_j u}_1 \right)$ promotes the sparsity structure of $u$ within each of the subbands $\Psi_j$ sufficiently. Indeed, as it was argued in \cite{CanWakBoy}, the log-sum penalty is more sparsity enforcing than the $\ell^1$-norm by putting a larger penalty on small nonzero coefficients. In the case of the $\ell^1$-norm of an entire subband this approach seems to be less effective in promoting the sparsity within each level. It was furthermore proposed in \cite{AhmSch} (Co-IRW-L1, Algorithm 4) to combine the update rule \eqref{Schniter} with the classical elementwise reweighting update
\begin{align}
\label{IRL1}
W_{j} = \diag \left( \frac{1}{\varepsilon + \left| \langle \psi_{j,l},u^{k} \rangle \right| } \right),
\end{align}
for $j=1,\dots,J$.
However, note that this combination is very different to what we are aiming for, since there is even more emphasize put on penalizing the smaller coefficients in higher levels which can happen to delete too many highscale coefficients.

\begin{figure*}
\includegraphics[width=0.49\textwidth]{./Schniter.png}
\includegraphics[width=0.49\textwidth]{./Quantile.png}
\\[1ex]
\includegraphics[width=0.49\textwidth]{./New.png}
\includegraphics[width=0.49\textwidth]{./NewZoom.png}
\caption{\textbf{Upper left}: Coefficients (in blue) with a multi-level reweighting curve as in \eqref{Schniter} (in orange). \textbf{Upper right}: Coefficients (in blue) with multi-level curve (in orange) using \eqref{our} only. \textbf{Lower left}: Coefficients with the proposed reweighting strategy (in orange). \textbf{Lower right}: Zoom of proposed multilevel reweighting strategy.} \label{Reweighting}
\end{figure*}

This fact is visualized from a different point of view in Figure \ref{Reweighting}, where we have depicted the shearlet coefficients of a MRI phantom introduced in \cite{GLPU} together with the reweighting rule we have just discussed in the top-left of the figure.  The shearlet coefficients are depicted in blue and the values of $\frac{\lambda_j}{\mu_1}$ for a realistic value $\mu_1$ are shown in orange. Note that according to the update rule  \eqref{shrinkw} of the split Bregman algorithm everything below the orange curve would be thresholded.

\subsection{Multilevel adapted reweighting}{Sec:IIIc)}
Considering the previous discussion one of the disadvantages is that the weights corresponding to higher levels might become too large. This can be prevented, for instance, by choosing the regularization parameters as
\begin{align}
\label{our}
\lambda_j =  \max\left\{\left| \langle \psi_{j,l}, u \rangle \right| : l = 1,\dots,N_j \right\},
\end{align}
for $j =1, \ldots, J$ and zero otherwise, i.e., if $j =0$. Notice that the necessity of such weights was already elaborated in Section \ref{sec:inpainting}.
Note that we set $\lambda_0=0$, since for real life signals the low frequency part is usually not sparse, see also Figure \ref{Sparsity2}. This was also observed in  \cite{SelFig}, where it was shown that this idea can be accomplished more effectively using an analysis instead of a synthesis prior. The attentive reader might wonder that the magnitude of the analysis coefficients could be very irregular per level, in particular, one could have strong outliers. For such cases it is better to take a quantile instead of the maximum over all coefficients. Note that this choice of $\lambda_j$ is a heuristic substitute accounting for the unknown constant in the theoretical decay of the multilevel coefficients.  A schematic representation from the thresholding perspective of split Bregman can be found in the second image of Figure \ref{Reweighting}.

\begin{figure*}
\includegraphics[width=0.495\textwidth]{./RelErrTrue.png}
\includegraphics[width=0.495\textwidth]{./SSIMTrue.png}
\caption{Convergence plot: Behavior of the error and structured similarity index for reweighting with true shearlet coefficients with respect to increasing number of iterations. Used signal: phantom from \cite{GLPU}. Reconstructed from 6\% of Fourier data with proposed algorithm without TGV (see section \ref{sec:Section4} for details). In IRL1 we are choosing $W_j$ as in \eqref{IRL1} and $\lambda_j$ constant and for the proposed emthod we additionally define $\lambda_j$ as in \eqref{our}. \textbf{Left figure}: Relative error of in each iteration.  \textbf{Right figure}: SSIM of each iteration.} \label{true}
\end{figure*}

Our proposed method combines the classical reweighting of \eqref{IRL1} with the above choice for $\lambda_j$. The idea behind this is that we are still using the power of pointwise iterative reweigthing, but since our multiscale coefficients naturally come in a decreasing order in magnitude, we apply it to each level separately weighted with $\lambda_j$. That means that within each level we follow the democratic philosophy of reweighting which is that small analysis coefficients of the current iterate $\Psi_j u^{\cur}$ are likely to be zero in $\Psi_j u$. By multiplying with the latter choice of $\lambda_j$ we  also gain more control and account for the multilevel structure of $\Psi u$. An artificial experiment using the (in reality unknown) true analysis coefficients demonstrates that this update rule does seem to perform better than standard reweighting without such a compensation of multilevel weights, see Figure \ref{true}. For the signal $u$ we are choosing the phantom of \cite{GLPU}. Using the
usually unknown shearlet coefficients for the construction of $\lambda_j$ and $W_j$ for $j=1,\dots,J$ as explained above we are reconstructing $u$ from only $6 \%$ of its Fourier measurements obtained by radial lines through the k-space origin by using the proposed algorithm of the last section, but without the additional TGV regularizer. Note that both reconstructions approximately start with the same error, which indicates that the set of regularization parameters is chosen equally good. However, it can be said that tuning the iterative reweighted shrinking method within the split Bregman framework without the automatic choice of the subband-weights $\lambda_j$ is rather difficult and highly signal dependent. This example further shows the potential of reweighting as the data is significantly undersampled.

 \subsection{Proposed algorithm}

Having explained the idea of our method, we now state the final resulting algorithm that is a composition of the split Bregman framework for solving the constrained optimization problem \eqref{discrObj} and the previously explained idea of multilevel weighting and iteratively reweighting, respectively.
In contrast to the traditional reweighted $\ell^1 $ approaches we do not iterate between solving the $\ell^1$-problem up to convergence and updating the weights. We propose to incorporate the multilevel adapted reweighting rule directly into the split Bregman algorithm. This is done in such a way that only the shrinking of the $w$-subproblem is changed to a \emph{multilevel adapted, iteratively reweighted shrinking rule}. Note that by choosing the level-weights $\lambda_j$ depending on the magnitude of the signal coefficients the resulting method appears to stable towards the alternation of signals.

\begin{algorithm*}
\caption{Proposed algorithm}
\label{alg:buildtree}
\begin{algorithmic}
\STATE{\underline{Input}: \\
Measurement operator $A$, multilevel transform $\Psi$,\\
regularization parameters: $\alpha_0, \alpha_1,\mu_1,\mu_2,\mu_3,\beta,$\\
iteration numbers $N$ and maxIter.\\
\underline{Data}:\\
Measured data $y$.\\
\underline{Initialization}:\\
$ k \gets 0$;\\
$ u_0 \gets A^*y$;\\
$ y_0,v_0, d_0, w_0,t_0,b_0^t,b_0^d,b_0^w \gets 0$;
}
\WHILE{$k \leq \maxIter$}
\FOR{$i=1,\ldots,N$}
\STATE{$(u_{k+1},v_{k+1}) \gets$ \text{solve linear system \eqref{diagsys}};}
\FOR{$j=1,\ldots,J$}
\STATE{$\lambda_j =  \max\left\{\left| \langle \psi_{j,l}, u \rangle \right| : l = 1,\dots,N_j  \right\}$;\\
$W_{j} = \diag \left( \frac{1}{\varepsilon + \left| \langle \psi_{j,l},u^{\cur} \rangle \right| } \right)$;\\
$w_{k+1}^j(l) \gets \text{shrink}\left(\left(\Psi_j u^{\cur}\right)(l) +  b_k^{w,j}(l), \frac{\lambda_j W_j(l)}{\mu_1} \right), \quad l=1,\dots,N_j$};\\
\ENDFOR
\STATE{$d_{k+1}(l) \gets \text{shrink}_2\left( \nabla^f u^{\cur}(l) - v^{\cur}(l) + b_k^d(l),\frac{\alpha_1}{\mu_2} \right),\quad l=1,\dots,N_j$;\\
${t_{k+1}(l) \gets \text{shrink}_F\left(\left(\mathcal{E}^b v^{\cur}\right)(l) + b^t_k(l), \frac{\alpha_0}{\mu_3} \right),\quad l=1,\dots,N_j}$;
}
\ENDFOR
\STATE{
$b^w_{k+1}  \gets b^w_k + \Psi u_{k+1} - w_{k+1}$;\\
$b^d_{k+1}  \gets b^d_k + (\nabla^f u_{k+1} -v_{k+1}) -d_{k+1}$;\\
$b^t_{k+1}  \gets b^t_k + \mathcal{E}^bv_{k+1} - t_{k+1}$;\\
$y_{k+1}  \gets y_k + y - Au_{k+1}$;\\
$k \gets k+1$;}
\ENDWHILE
\RETURN Reconstruction $u_{\text{maxIter}}$.
\end{algorithmic}
\end{algorithm*}

We like to comment on two things regarding the algorithm above. First, the weighting matrix $W_j$ depends on the initialization of an $\varepsilon >0$. The choice of $\varepsilon$ is rather empirical and in many cases does not effect the solution. This was already noticed in the beginning of reweighted $\ell^1$ in \cite{CanWakBoy}. This is also the case for our algorithm. As we explained in Section \ref{sec:Section2} the role $\varepsilon$ is essentially just the maximum threshold that our algorithm will do. It provides stability by preventing a division by zero, but the magnitude of the coefficient is mostly determined by the respective analysis coefficient. Let us also mention that for obvious reasons, adapting the $\epsilon$ with respect to the levels cannot achieve the same effect as the proposed choice of the parameters $\lambda_j$.

Furthermore, we have not incorporated an additional stopping criterion besides a maximum number of iterations. In all our test cases the algorithm shows convergence, see Figure \ref{true}, \ref{GLPU}, \ref{Pepper} and \ref{CT}.

\section{Numerics} \label{sec:Section4}

In this section we will recover numerous signals from their Fourier measurements and  Radon measurements. These are classical problems in applied mathematics where the image of interest is sparse under a multiscale transform such was the wavelet transform or the shearlet transform. One of the most known applications for the recovery of Fourier measurements is magnetic resonance imaging (MRI) where data is collected in the so-called \emph{k-space} which is the Fourier domain, i.e. every point in the k-space can be interpreted as a Fourier coefficient of the object of interest. This is also one of the very first areas where compressed sensing has had a great impact on, see, for instance \cite{LusDonPau}. Our implementation is, however, not restricted to these two types of measurements or sparsifying transforms and other cases can be directly tested.

The numerical performance of our algorithm will be judged by the following three criteria that will be the main points of our analysis:
\begin{itemize}
\item[(N1)] Quality of the reconstruction,
\item[(N2)] Convergence of the algorithm,
\item[(N3)] Stability towards change of signals.
\end{itemize}
 For (N1) we compare our algorithm with different existing and established methods that are known to perform well in the recovery problem from Fourier measurements these methods are shown in Table \ref{Abbrev}. The results are then shown in Section \ref{Sec:Compare}. Moreover, in some cases we will only use a Haarwavelet-based wavelet transform and in some other cases the shearlet transform. This is because we are not interested in promoting a particular transform, but rather show the effectiveness of our algorithm. Therefore we have always chosen the transform that performs best on the specific data set. Only in Figure \ref{Compare} we will show both methods in a direct comparison.

 For (N2) we have chosen two quality measurements. First, the relative error which is computed by the formula
 \begin{align*}
    \text{RE} = \frac{\| u^{\text{ref}} - u^{\text{rec}} \|_2}{\| u^{\text{ref}}\|_2},
 \end{align*}
 where $u^{\text{ref}}$ is the reference image and $u^{\text{rec}}$ the reconstructed image, both in a vectorized form. Second, we use the structured similarity index that as introduced in \cite{SSIM}.

 The stability (N3) is verified by the fact that we have chosen the same parameters for each multiscale transform across all experiments for each problem. Although an extensive tuning of all parameters for different images might yield superior results we have chosen not to do so. The reason behind is two-fold: First, our algorithm already performance very well with a fixed choice of parameters for all different images used in this section. Second, iterative reweighting combined with the proposed multilevel weighting strategy already suggests the level of thresholds for all coefficients and should therefore be less sensitive to the choice of additional parameters.

\subsection{Recovery from Fourier measurements}

In this section we give a description of which multiscale transforms we used precisely and the parameters.

For \emph{wavelets} we have used the undecimated 2D wavelet transform of the Spot package available at
\begin{center}
 \texttt{http://www.cs.ubc.ca/labs/scl/spot/}
\end{center}
Otherwise differently stated we have used a 4 scale Daubechies-2 wavelet system.
The chosen parameters for the proposed algorithm (WIRL1 + TGV) are
\begin{itemize}
\item[-] $\mu$: [6e2 1e1 2e1],
\item[-] $\alpha$: [1 2],
\item[-] $\beta$: 1e4,
\item[-] $\varepsilon$: 1e-4.
\end{itemize}

For \emph{shearlets} we have used the shearlet transform available from
\begin{center}
 \texttt{http://www.shearlab.org/}
\end{center}
The discrete shearlet system is generated by using 4 scales and $[1 \, 1 \,  2 \, 2]$ for the directional parameters. The chosen parameters for the proposed algorithm (SIRL1 + TGV) are
\begin{itemize}
\item[-] $\mu$: [5e3 1e1 2e1],
\item[-] $\alpha$: [1 1],
\item[-] $\beta$: 1e5,
\item[-] $\varepsilon$: 1e-5.
\end{itemize}
For all experiments we chose 2 number of block Gauss-Seidel iterations  and performed 4 inner iterations before updating  $y_k$.  All experiments were conducted in MATLAB R2015b with an Intel i3 CPU with 8GB memory.

\subsection{Comparison with other methods}\label{Sec:Compare}

\begin{table*}
\centering
\begin{tabularx}{1\textwidth}{|X|X|}
\hline
Fourier inverse & Fourier inversion of data\\
\hline
RecPF & Total variation and wavelet regularization from \cite{RecPF}\\
\hline
FFST+TGV & Shearlets with TGV from \cite{GuoQinYin}\\
\hline
Co-IRL1 & Composite iterative reweighting from \cite{AhmSch}\\
\hline
PANO & Patch-based nonlocal operator, \cite{QuHouLamGuoZhongChen}\\
\hline
TV & Total variation with our algorithm \\
\hline
TGV &  Total generalized variation with our algorithm\\
\hline
WL1 & Wavelets without reweighting\\
\hline
WIRL1 & Wavelets with proposed reweighting\\
\hline
WIRL1+TGV & Wavelets with proposed reweighting and TGV\\
\hline
SL1 & Shearlets without reweighting\\
\hline
SR + TGV & Shearlets with standard reweighting and TGV\\
\hline
SIRL1 & Shearlets with proposed reweighting\\
\hline
SIRL1+TGV & Shearlets with proposed reweighting and TGV\\
\hline
\end{tabularx}
\caption{Table for abbreviations}\label{Abbrev}
\end{table*}

Our first experiment shows the recovery from a $256 \times 256$ rose image available from the open source framework
\begin{center}
 \texttt{http://aforgenet.com/framework/}
\end{center}
We took 30 radial lines through the k-space origin ($\approx 12,2\%$ of Fourier data) representing the measurement data.  We then compared our results obtained by our algorithm using wavelets as well as shearlets.  For both results we have used the proposed iterative multi-level reweighting and a generalized total variation regularizer. The results are compared to RecPF by Yang et al. \cite{RecPF}, Co-IRL1 by Ahmad and Schniter \cite{AhmSch}, PANO by Qu et al. \cite{QuHouLamGuoZhongChen}, and FFST+TGV by Guo et al. \cite{GuoQinYin}.
In order to make the experiments comparable we have used the same scales and number of directions in \cite{GuoQinYin}. Furthermore, for Co-IRL1 we have used two redundant Daubechies wavelet dictionaries with the same number of scales. Moreover, one dictionary consists of Haar wavelets and the second one of Daubechies 2 wavelets.

The final result comparing all methods can be found in Figure \ref{Compare}. It can be observed that the recovery obtained
by the proposed method shows the least amount of artifacts and overall the highest quality of reconstruction.

\begin{figure*}
\hspace*{-.9cm}
\begin{tabular}{cccc}
\hspace*{0.4cm}\includegraphics[width=0.235\textwidth]{./truthR} & \hspace*{-0.4cm}
\includegraphics[width=0.235\textwidth]{./ztruthR.png}
&
\includegraphics[width=0.235\textwidth]{./zeroFR.png}
&
\hspace*{-0.4cm}
\includegraphics[width=0.235\textwidth]{./zzeroFR.png}
\\
\multicolumn{2}{c}{\footnotesize{Original}}& \multicolumn{2}{c}
{\footnotesize{Fourier inverse: \quad RE = 0.167, SSIM = 0.739 }}\\[1ex]
 \hspace*{0.4cm}\includegraphics[width=0.235\textwidth]{./RecPFR} & \hspace*{-0.4cm}
\includegraphics[width=0.235\textwidth]{./zRecPFR.png}
&
\includegraphics[width=0.235\textwidth]{./CoL1R.png}
&
\hspace*{-0.4cm}
\includegraphics[width=0.235\textwidth]{./zCoL1R.png}
\\
\multicolumn{2}{c}{\footnotesize{RecPF: \quad RE = 0.105, SSIM = 0.841}} & \multicolumn{2}{c}{\footnotesize{Co-IRL1: \quad RE = 0.076, SSIM = 0.898}}\\[1ex]
\hspace*{0.4cm}\includegraphics[width=0.235\textwidth]{./PanoR.png}
&
\hspace*{-0.4cm}
\includegraphics[width=0.235\textwidth]{./zPanoR.png}
&
\includegraphics[width=0.235\textwidth]{./FFSTR.png}
&
\hspace*{-0.4cm}
\includegraphics[width=0.235\textwidth]{./zFFSTR.png}
\\
\multicolumn{2}{c}{\footnotesize{PANO: \quad RE = 0.086, SSIM = 0.881 }} & \multicolumn{2}{c}{\footnotesize{FFST+TGV: \quad RE = 0.095, SSIM = 0.855}} \\[1ex]
\hspace*{0.4cm}\includegraphics[width=0.235\textwidth]{./TGVR.png}
&
\hspace*{-0.4cm}
\includegraphics[width=0.235\textwidth]{./zTGVR.png}
&
\includegraphics[width=0.235\textwidth]{./WIRL1TGVR.png}
&
\hspace*{-0.4cm}
\includegraphics[width=0.235\textwidth]{./zWIRL1TGVR.png}
\\
 \multicolumn{2}{c}{\footnotesize{TGV: \quad RE = 0.113, SSIM = 0.874}} &  \multicolumn{2}{c}{\footnotesize{WIRL1+TGV: \quad RE = 0.058, SSIM = 0.936}}
\end{tabular}
\caption{Different reconstructions from 30 radial lines ($\approx 12,2\%$) through the k-space origin with relative error and structured similarity index. 50 iterations are used for the reconstruction. See Table \ref{Abbrev} for used abbreviations.} \label{Compare}
\end{figure*}

We will next extensively study our algorithm in terms of the effectiveness and stability.

\subsection{Convergence, signal independence, and the effect of reweighting}\label{Sec:Convergence}

In this section we analyze (N2) for our algorithm. We do this by considering two images, one that is well suited for wavelets and the other one where shearlets perform better. We start with a $256 \times 256$ phantom that was designed by Guerquin-Kern et al. in \cite{GLPU} for MRI studies. As this image is piecewise constant we have chosen a 4 scale wavelet transform generated by Haar wavelets. We reconstructed the image using our algorithm for TV only, WL1, WIRL1 and WIRL1+TGV. Moreover, as this image is very compressible in a Haar wavelet basis, the recovery allows a much lower sampling rate. In fact, we have only used 21 radial lines which corresponds to 8.73\%. It is interesting to mention that the exact solution is returned after almost 80 iterations when WIRL1+TGV. Note that, for 24 lines ($\approx 9.83 \%$) wavelets with the proposed iterative reweighting step (WIRL1) will eventually also return the exact solution, see Figure \ref{GLPU}. Also note that the TV reconstruction performs worst, although
this image should be well suited for TV. Our explanation is that at these extremely low sampling rates a highly redundant transformation is needed to still guarantee recovery.

\begin{figure*}
\footnotesize
\begin{tabular}{ccccc}
\includegraphics[width=0.17\textwidth]{./truthG.png}
&\includegraphics[width=0.17\textwidth]{./TV.png}
&\includegraphics[width=0.17\textwidth]{./WL1G.png}
&\includegraphics[width=0.17\textwidth]{./WIRL1G.png}
&\includegraphics[width=0.17\textwidth]{./WIRL1TGVG.png}
\\
\footnotesize{Original} & \footnotesize{TV} & \footnotesize{WL1} & \footnotesize{WIRL1} & \footnotesize{WIRL1+TGV} \\
\includegraphics[width=0.17\textwidth]{./ztruthG.png}
&\includegraphics[width=0.17\textwidth]{./zTV.png}
&\includegraphics[width=0.17\textwidth]{./zWL1G.png}
&\includegraphics[width=0.17\textwidth]{./zWIRL1G.png}
&\includegraphics[width=0.17\textwidth]{./zWIRL1TGVG.png}
\\[2ex]
\multicolumn{5}{c}{\includegraphics[width=0.4\textwidth]{./RelErrGLPU.png}
\includegraphics[width=0.4\textwidth]{./RelErrGLPUp1.png}
}
\end{tabular}
\caption{Different reconstructions from 24 radial lines through the k-space origin with relative error and structured similarity index. The lower left graphics corresponds to 21 radial lines ($\approx 8.73 \%$) and the lower right to 24 radial lines ($\approx 9.83 \%$). 100 iterations are used for the reconstruction. See Table \ref{Abbrev} for used abbreviations.} \label{GLPU}
\end{figure*}

Our third numerical example for Fourier measurements concerns the $256 \times 256$ pepper image, see Figure \ref{Pepper}. It has many more structures than the previously considered GLPU phantom. More importantly, it does not consist of piecewise constant areas. This image is particularly well suited for shearlets and thus we have chosen the shearlet transform  with four scales. The result for a fixed threshold (SL1), i.e. without the proposed reweighting is significantly worse than the one obtained by the multilevel iterative reweighting method (SIRL1). In this case adding total generalized variation as an additional regularizer does not improve the image quality much. Indeed, the improvement of adding TGV as a regularizer depends strongly on the image. It improves the quality if more piecewise constant areas are present in the image (especially if the background is constant as this is typically the case for MRI images), see for instance Figure \ref{GLPU}.

\begin{figure*}
\begin{tabular}{cccc}
\includegraphics[width=0.22\textwidth]{./truthP.png}
&\includegraphics[width=0.22\textwidth]{./SL1TGVP.png}
&\includegraphics[width=0.22\textwidth]{./SIRL1TGVP.png}
&\includegraphics[width=0.22\textwidth]{./SNewIRL1TGVP.png}
\\
\footnotesize{Original} & \footnotesize{SL1} & \footnotesize{SR+TGV} & \footnotesize{SIRL1+TGV}
\\
\includegraphics[width=0.22\textwidth]{./ztruthP.png}
&\includegraphics[width=0.22\textwidth]{./zSL1TGVP.png}
&\includegraphics[width=0.22\textwidth]{./zSIRL1TGVP.png}
&\includegraphics[width=0.22\textwidth]{./zSNewIRL1TGVP.png}
\\[2ex]
\multicolumn{4}{c}{\includegraphics[width=0.4\textwidth]{./RelErrP.png}
\includegraphics[width=0.4\textwidth]{./SSIMP.png}
}
\end{tabular}
\caption{Different reconstruction from 30 radial lines ($\approx 12,2\%$) through the k-space origin with relative error and structured similarity index. See Table \ref{Abbrev} for used abbreviations. 50 iterations are used for the reconstruction.} \label{Pepper}
\end{figure*}

\subsection{Timings}

Although our code is not optimized at all, reconstructions are obtained in reasonable time. Below we display the timings that our algorithm needs in order to compute a reconstruction. The timings are recorded from the experiment shown in Figure \ref{Compare} (image size: $256\times 256$). As our algorithm is also capable of computing reconstructions using only a TV- or TGV-regularizer we additionally recorded these timings.

\begin{table}[H]
\centering
\begin{tabular}{|c|c|}
\hline
 Method & Time in seconds\\
\hline
TV & 16.97\\
\hline
TGV &  47.01\\
\hline
WL1 & 254.77 \\
\hline
WIRL1 & 310.78\\
\hline
WIRL1+TGV & 350.61 \\
\hline
SL1 & 286.57\\
\hline
SIRL1 & 419.20\\
\hline
SIRL1+TGV & 479.74\\
\hline
\end{tabular}
\caption{Table of timings.}\label{Timings}
\end{table}

Note that TGV is more time consuming than TV, mostly due to the larger system \eqref{diagsys}. Moreover, as we already mentioned, the wavelet transform as well as the shearlet transform used in this paper are redundant transforms, i.e. the number of coefficients is in this case significantly larger than the number of pixels available. This explains the much slower performance in comparison to the variational methods.

We also like to mention that the additional cost of reweighting (mostly due to a rather unsophisticated implementation) seems manageable, especially given the benefit of quality that one gets in the recovered images.

\subsection{Recovery from Radon measurements} \label{sec:Radon}
In order to demonstrate the versatility of our proposed method, we furthermore consider an example of a recovery from few Radon measurements. This is another classical linear inverse problem, which forms the underlying model for \emph{computed tomography}. As a test signal we again use the brain phantom of \cite{GLPU}. For the measurement system we use  \emph{fan-beam} projections with 45 equally spaced rotation angles. For the reconstruction of the piecewise constant test signal we make use of the same 4 scale redundant Daubechies-2 wavelet system as above, together with the following parameter setup:
\begin{itemize}
\item[-] $\mu$: [1e3 1e2 2e3],
\item[-] $\alpha$: [1e-3 2e-3],
\item[-] $\beta$: 1e1,
\item[-] $\varepsilon$: 1e-6.
\end{itemize}
Since the matrix $A^*A$ is not diagonalized by $\mathcal{F}$ anymore, we follow the the steps described in Section \ref{sec:Section3} for solving \eqref{diagsys}. The resulting linear system is solved iteratively using 75 \texttt{pcg} iterations.

We compare the result with WL1, WIRL1, TV, TGV, Tikhonov regularization and filtered backprojection using a \texttt{`Hann'} filter as provided in Matlab.
The results are displayed in Figure \ref{CT} and the timings (for 150 iterations) are displayed in Table \ref{TimingsCT}. It can be observed that the convergence is slower as in the Fourier example. The fact that the timings are much higher than before is mostly due to the fact that \eqref{diagsys} is not diagonalizable. Although the Radon-Wavelet measurement matrix is known to be coherent \cite{Sampling,GuaCzaAroLea} (and thus classical CS theorems do not apply), we find it quite remarkable that the measurements contain enough information for an (almost) perfect recovery.

\begin{figure*}

\begin{tabular}{cccc}
\hspace*{0.4cm}\includegraphics[width=0.215\textwidth]{truth} & \hspace*{-0.4cm}
\includegraphics[width=0.215\textwidth]{truthZ.png}
&
\includegraphics[width=0.215\textwidth]{Filtered.png}
&
\hspace*{-0.4cm}
\includegraphics[width=0.215\textwidth]{FilteredZ.png}
\\
\multicolumn{2}{c}{\footnotesize{Original}}& \multicolumn{2}{c}
{\footnotesize{Filtered backprojection: \quad RE = 0.507, SSIM = 0.277 }}\\[1ex]
 \hspace*{0.4cm}\includegraphics[width=0.215\textwidth]{Tikhonov.png} & \hspace*{-0.4cm}
\includegraphics[width=0.215\textwidth]{TikhonovZ.png}
&
\includegraphics[width=0.215\textwidth]{CTTV.png}
&
\hspace*{-0.4cm}
\includegraphics[width=0.215\textwidth]{TVZ.png}
\\
\multicolumn{2}{c}{\footnotesize{Tikhonov: \quad RE = 0.264, SSIM = 0.479}} & \multicolumn{2}{c}{\footnotesize{TV: \quad RE = 0.094, SSIM = 0.966}}\\[1ex]
\hspace*{0.4cm}\includegraphics[width=0.215\textwidth]{TGV.png}
&
\hspace*{-0.4cm}
\includegraphics[width=0.215\textwidth]{TGVZ.png}
&
\includegraphics[width=0.215\textwidth]{Haar.png}
&
\hspace*{-0.4cm}
\includegraphics[width=0.215\textwidth]{HaarZ.png}
\\
\multicolumn{2}{c}{\footnotesize{TGV: \quad RE = 0.094, SSIM = 0.967}} & \multicolumn{2}{c}{\footnotesize{WL1: \quad RE = 0.101, SSIM = 0.956}} \\[1ex]
\hspace*{0.4cm}\includegraphics[width=0.215\textwidth]{HaarIRL1.png}
&
\hspace*{-0.4cm}
\includegraphics[width=0.215\textwidth]{HaarIRL1Z.png}
&
\includegraphics[width=0.215\textwidth]{HaarTGVIRL1.png}
&
\hspace*{-0.4cm}
\includegraphics[width=0.215\textwidth]{HaarTGVIRL1Z.png}
\\
 \multicolumn{2}{c}{\footnotesize{WIRL1: \quad RE = 0.058 SSIM = 0.985}} &  \multicolumn{2}{c}{\footnotesize{WIRL1+TGV: \quad RE = 0.009, SSIM = 0.999}}
 \\[2ex]
\multicolumn{4}{c}{\includegraphics[width=0.42\textwidth]{CTRelErr.png}
}
\end{tabular}
\caption{Different reconstructions from 45 angles fan-beam projections with relative error and structured similarity index. 150 iterations are used for the last five methods. See Table \ref{Abbrev} for used abbreviations.} \label{CT}
\end{figure*}

\begin{table}
\centering
\begin{tabular}{|c|c|}
\hline
 Method & Time in seconds\\
\hline
Filtered backprojection & 0.5\\
\hline
Tikhonov &  0.67\\
\hline
TV & 718 \\
\hline
TGV & 1208\\
\hline
WL1 & 533 \\
\hline
WIRL1 & 703\\
\hline
WIRL1+TGV & 1117\\
\hline
\end{tabular}
\caption{Table of timings for Radon measurements.}\label{TimingsCT}
\end{table}

\section{Conclusion}\label{sec:Section5}

In this paper we presented a novel split Bregman based algorithm that incorporates an iteratively reweighted shrinkage step that is specifically adapted to the structure of multiscale coefficients in order to enhance the quality of image reconstruction. We presented an extensive study of the algorithm thereby focusing on the quality of the image that is to be reconstructed as well as the stability of the parameters towards the change of different signals. Although our focus in the numerical experiments was on the reconstruction problem from Fourier measurements, our algorithm can also be applied to other common problems such as denoising and inpainting.

Even though all experiments where based on a 2D multiscale transform our findings can be carried over to higher dimensions. This is demonstrated in case of real-world 3D MR data with noise in a practical paper \cite{MaMaeFunSchKutSchKol} using 3D shearlets, \cite{3DShearlets}.

\section*{Acknowledgements}
This work was done when J. Ma was affiliated with the Technische Universit\"{a}t Berlin, he acknowledges support from the DFG Collaborative Research Center TRR 109 "Discretization in Geometry and Dynamics". JM and MM acknowledge support from the Berlin Mathematical School. The authors would further like to thank Gitta Kutyniok for helpful discussions. Moreover, the authors would like to thank Professor Rizwan Ahmad and Professor Philip Schniter for providing a code for the Co-IRL1 algorithm.

\bibliographystyle{abbrv}

redpajama_set_name: RedPajamaArXiv


text:
When one talks about the budget or mid-range Android smartphones invariably one tends to gravitate towards a Chinese phenomenon called Xiaomi. Over the last few years, it has made a mockery of what one should expect in a smartphone that costs a third of an iPhone. Ironically, this comes from a company which has been dubbed the Apple of China. The Mi A1 launch also confirms that Google is still committed to its Android One program. Xiaomi's new Android One phone could be the first of many in the pipeline. Apart from its stock Android experience, the Mi A1 is also Xiaomi's first dual-camera device in India. It packs a decent set of innards and is being marketed as the new flagship. We received a Xiaomi Mi A1 unit ahead of the launch, and here are our first impressions.
At first look, the Mi A1's finish and overall look instantly tell us that it belongs to the more premium Mi range and not the budget Redmi series from Xiaomi. The metal unibody and 2.5D curved-edge glass at the front make it look good. The rounded corners are comfortable in the palm.
In a nutshell, the Mi A1's design can be called a lovechild between the OnePlus 5 and the Xiaomi Redmi Note 4. The facade of the phone will remind you of the Redmi Note 4 which shouldn't surprise anyone as they are the same family of device, however, its corners are more rounded which is kind of amazing as it looks nicer than the Redmi. On the back, it looks a lot like the OnePlus 5. It has a similar style of a camera visor/hump, and the antenna lines are also similarly designed.
The front is bare of any branding, which we think looks good. You'll see the Mi logo and the Android One tag towards the bottom of the rear. The Xiaomi Mi A1 will eventually be manufactured locally and will sport a "Made in India" tag, but initial units will be shipped in from different countries. The antenna lines are discreet and the design is different from what we have seen on the Xiaomi Redmi Note 4 (Review). In our opinion, the Xiaomi Mi A1 feels the best in hand out of the entire current Xiaomi product line.
The Mi A1 features a 5.5-inch full-HD display with 2.5D glass. It is powered by an octa-core Qualcomm Snapdragon 625 processor coupled with 4GB of RAM. It also has 64GB of storage and supports expansion using a microSD card. The Mi A1 has a hybrid dual-SIM design, and support two Nano-SIMs if you aren't using a microSD card. There's a fingerprint scanner on the back. The smartphone also has a USB Type-C port and packs a 3080mAh non-removable battery.
One of the most marketed features of the Mi A1 is its cameras. It sports two 12-megapixel sensors at the back, one with a wide-angle lens and f/2.2 aperture, and the other with a telephoto lens and f/2.6 aperture. This is similar to the setup of the flagship Xiaomi Mi 6.
The Mi A1 supports 4G with VoLTE (Voice over LTE) as well as Wi-Fi, Bluetooth, and a 3.5mm audio jack.
In terms of software, Xiaomi's Mi A1 runs pure Android 7.1.2 Nougat which is the most recent Nougat version available for non-Google-branded Android smartphones.Google's Android One program was created to do two things — one to reduce software fragmentation of Android wherein there are many devices running outdated versions of the operating system which could also be potential security hazards and secondly for reaching a baseline benchmark for what an Android smartphone should be like.
Apart from the Google suite of apps, the Mi A1 also has some software from Xiaomi such as the Mi Remote and Mi Store apps. At least Xiaomi's fans won't completely miss out on features seen on other Xiaomi products. The Settings app also has a Mi Services section that can be enabled to send the company diagnostic data about how you use your phone.
Stock Android is not the only USP for this device. It's also the only Android One phone to sport a dual-camera. There are two 12MP sensors on the back of this device.Essentially, Xiaomi has taken the dual 12-megapixel camera stack from the Mi 6 and outfitted it on the Mi A1. This, however, is inferior as it has a shallower aperture of f/2.2. That said, it does a 2x optical zoom via the 12-megapixel telephoto lens like the iPhone 7 Plus and the OnePlus 5. The main wide-angle camera is also great and when used in tandem there's a creamy and effective portrait mode which is better than the one on the OnePlus 5.
Camera Samples
The portrait mode implementation is also a very good one. It has great edge detection and you get this natural creamy bokeh effect on the background. Sure, it isn't on the same level as the iPhone 7 Plus, but it is better than the OnePlus 5 which praiseworthy as this is a phone that costs half as much.
Mi A1 Picture Perfect Dual Camera
Xiaomi's Mi A1 not only has a lot riding on it for the Chinese company perspective but also for Google. It will be interesting to see the response to this phone in India and other markets where it's launching.
download nulled
Tags Mobile
Nvidia is going to release a new generation of graphics card named Nvidia GeForce GTX …


Here at Cleveland Cascades Ltd we produce high quality goods with in-house processes for the fabrication, manufacturing and assembly of our chutes and their component parts.
Our fabrication workshop produces bespoke design components using specialist equipment such as punching, bending and welding rigs. Some of our fabrication processes can be seen in the Gallery below.
As well as fabrication, we also produce our own GRP cones for use with our Cascade and FFC chutes.
Using a mould we Spray the GRP cones to size before trimming, cleaning, polishing and finally lining. the pictures below show the cone spraying booth and lining workshop with cones in the process of being manufactured and lined with UHMWPE.
We also apply lining to many our other manufactured components such as our head chutes and trimmer extensions as seen below in the gallery with a completed ceramic lined trimmer.
Once each of the sections of shroud are manufactured they are then moved into our assembly workshop where they are finished into a full assembly and finally quality inspected before packing.
Once all of our items have been manufactured they are then quality inspected and where necessary completed as an assembly in our workshop by our trained cascade chute engineers. We also complete torque tightening on all required areas whilst assembling as seen in the gallery below with our engineers working on one of our 1700 carrier assemblies and on the electrical components of one of our 900 size head chutes.
Once we have all items assembled we then complete trial fitting and function testing to ensure a smooth installation process once the chute is on site for the customer/end user.


text:
Email this to a friend:
Obama's War
FRONTLINE is calling Afghanistan "Obama's War." The PBS documentary series kicks off its new season by playing out the wartime scenarios and the diplomatic frustrations in a war that has three fronts – the battlefield, the afghan leadership, and the Pakistani government. It's a firsthand look at how the president's new strategy is taking shape, delivering vivid, on-the-ground reporting from this war's many fronts and the internal debates over President Obama's grand attempt to combat terrorism at its roots.
Send This To:
Terms & Conditions | Privacy | Website Feedback | Copyright 2017 © Southern California Public Radio


the chilling effect on job creation from the NLRB's actions.
organizing elections. The bottom-line is that the NLRB's actions undermine employer-employee relations.
impossible for anyone to challenge the bargaining unit chosen by the union. The decision will divide employees.


How to survive the land down under's record-breaking heatwave. Source: Shutterstock.
AUSTRALIA'S summer is so abnormally hot right now that a slab of porterhouse steak could cook in a car parked in a shade.
In fact, the land down under's heatwave is so extreme, it is reaching record-breaking temperatures of 50c in some parts. And according to Australia's Bureau of Meteorology (BoM), it is expected to blaze on for a few more weeks.
The intense heat is brought on by the delayed onset of monsoon rains in the north, and the lack of typical cool fronts spinning up from the Antarctic in the south have exacerbated the warmth, according to the BoM.
The high temperatures are taking an increasing toll on the country's flora and fauna, CNN reported. More than a million dead fish have washed up on the banks of the Murray-Darling River Basin.
Colonies of bats and an estimated 23,000 spectacled flying foxes have also succumbed to the heat and dropping from trees.
In South Australia, the extreme heat is causing stone fruit, including peaches and nectarines, to cook from the inside out, according to ABC.
While a cool-down is expected, it will be shortlived.
The rest of South Australia is also headed up this week; into the mid-40s probably from around Tuesday, Narramore added.
The BoM provides timely updates of the heatwave, urging people to stay cool, check on loved ones, and follow advice from health authorities.
The New South Wales government has released a health alert warning people to stay indoors during the hottest part of the day, minimize physical activity and keep hydrated.
"Heat puts a lot of strain on the body and can cause dehydration, heat stress and heat stroke. It can also make underlying health conditions worse. People over 75, people with chronic medical conditions and people who live alone are particularly vulnerable," the statement said.
If the sweltering heat can cook fruits from inside out, it can also do damage to the human body. So remember to drink lots of water, slap on heaps of sun protection, and stay indoors as much as possible.


'Tis the season to overspend. Each year, we promise ourselves that we're not going to go overboard with holiday spending. Each year, it seems it's difficult to pay the bills in January and beyond. But you can be festive and save money. Here's how.
Go branch hunting outside. Spray-paint the branches in a silver or gold, pop them in a vase and voilà, elegant holiday decor. Take it a step further and dangle some ornaments on them or tie a few bows.
Glue leftover cards or those from previous years on a board as holiday artwork. The boards could be trimmed down to the size of a single card or left as a whole board with numerous cards attached to it lined up against a wall or on top of a coffee table.
Take your throw pillows and tie them with holiday coloured ribbons either once around the middle or like a Christmas gift. Add a rhinestone pin for an extra kick!
Fill large tubular vases with red cranberries to add a beautiful splash of Christmas red to any room. If you don't have vases, take three or more clear glass wine bottles or other household bottles and fill them with cranberries. You can even take it a step further and add flowers, twigs, or greenery with water. Another easy tip is to fill candle trays with cranberries or even pour them into the centre of your dining table with candles.
Take clear vases or bottles, fill them with water, add coloured food colouring — a few drops will do — and float candles in them or just line several up in a row to add a simple, cool holiday accent.


text:
GREAT AMERICAN BITES
Memphis' historic music room returns with soul food
Larry Olmsted
Special for USA TODAY
The scene: Memphis is renowned as one of the top cities in the world for enjoying live music, and Lafayette's Music Room combines this with a historical trip back in time. Before the famous Beale Street Entertainment District was created as a pedestrianized cluster of bars and venues, Overton Square was the city's beating heart spot for live music and theater. After years of decline, the neighborhood is suddenly hot again, arguably the trendiest in Memphis, with new restaurants, shops and brewpubs. Back in the '70s, the last heyday of Overton Square, Lafayette's was the place that up-and-coming acts wanted to play and the likes of Billy Joel, KISS, Barry Manilow and Styx all took the stage here. The bar closed for 38 years before reopening in 2014 amidst the neighborhood's resurgence. It has since become one of the city's hot spots for live music fans again.
The new Lafayette's is in the heart of Overton Square and spans two levels with outdoor seating and dining on both. The interior of the ground floor is one large room with the main stage on one side, a long bar on the other and floor tables across the room. Upstairs there's a horseshoe-shaped balcony with seating wrapped around and views of the stage below. The stage itself is a serious setup with state of the art sound and light systems. The entire place has exposed wooden walls, dark wooden tables and chairs, and a worn but upscale club feel. Lafayette's' features two live acts nightly, with a focus on promoting local talent.
Reason to visit: Live music, Pimento cheese waffle fries, po' boys, southern pizzas
The food: Not surprisingly for Memphis, Lafayette's Music Room serves a modern and slightly upscale version of traditional southern and soul food, with both upgraded classics and creative inventions. For an example of the latter, one of the more popular appetizers is a sort of southern spin on poutine, the classic Montreal comfort food of French fries with gravy and cheese curds, in turn a potato-based riff on nachos. Here they top crispy waffle fries with house-made pimento cheese, a southern staple, and your choice of crawfish or bacon. The crawfish was too unusual for me to pass up and the use of waffle fries, and their crispiness, helps a lot, as one of the major problems with any poutine-style dish is sogginess. Not here — this is a winner; it's really good, with the pimento cheese just spicy enough to have a kick but not so hot as to make it a challenge. The crawfish is a nice touch, as is serving it in a cast iron skillet, though I think it would be just as good with bacon. Another top starter with lots of regional flair is the fried green tomatoes, topped with corn maque choux, sort of a Cajun take on creamed corn with corn kernels cooked in lots of butter and cream along with chopped bell peppers, onions, garlic, herbs and spices, all drizzled with a Cajun remoulade — a very New Orleans spin on the classic fried green tomatoes. The tomatoes themselves were good, not the best I've had, but the topping was delicious and made the dish. Other notable starters include sweet corn and jalapeno hushpuppies, char broiled oysters and New Orleans-style barbecue shrimp.
Memphis' fast food barbecue chain is 'Tops'
The shrimp and grits is a good example of the entrees here: a haute version of the classic, a mound of dense grits in the middle of a large bowl surrounded by a ring of brothy, almost soup-like sauce brimming with shrimp, slices of andouille sausage and mushrooms. It's not the most traditional take on the dish but it's quite tasty and very moist, with a splash of sauce on every bite of grits. Even the simplest dishes seek authentic touches, like the bags of New Orleans' beloved Zapp's chips (Voodoo flavor) that come with the po' boy sandwiches, served on butcher paper. One real standout in this category is the fried shrimp po' boy — tasty shrimp perfectly breaded and cooked till crispy on very good bread. Even the extensive list of wood fired pizza has several southern takes on the genre, including a shrimp and andouille version, and my favorite, one with big cubes of barbecue pork belly, pickled red onions, cheddar cheese and a smoky chipotle barbecue sauce. The crust won't win any pizza awards, but the creative and delicious toppings make up for it.
The menu is very large and includes something for every taste from entrée salads to basics like grilled tuna, but to me the appeal is the southern spin, such as fried chicken on a waffle, pasta Jambalaya, or a unique main course comprised of an assortment of side dishes (you choose three) such as mac and cheese, smoked cheddar grits, roasted corn risotto, Spanish rice or corn maque choux, all served with a heathy side of Abita beer-flavored bread.
Lunch hour! Quirky Memphis eatery evokes elementary school
The main appeal of Lafayette's Music Room is the live music every night, and it is a big appeal, but the food, which could have been much more of an afterthought, is fun and tasty, and the setting is very pleasant. It doesn't hurt that you can start with a choice from a very creative craft cocktail list or enjoy numerous local beers on draught, and end with tempting southern dessert specialties like a layered mason jar take on Bananas Foster or classic New Orleans beignets (fried dough fritters) with confectioners sugar and decadent caramel sauce.
Pilgrimage-worthy?: Yes — for the music with good elevated southern comfort fare under one roof, it's a great place to get a feel of Memphis.
Rating: Yum! (Scale: Blah, OK, Mmmm, Yum!, OMG!)
Price: $$ ($ cheap, $$ moderate, $$$ expensive)
Details: 2119 Madison Ave., Memphis, TN; 901-207-5097; lafayettes.com/memphis/
Larry Olmsted has been writing about food and travel for more than 15 years. An avid eater and cook, he has attended cooking classes in Italy, judged a barbecue contest and once dined with Julia Child. Follow him on Twitter, @TravelFoodGuy, and if there's a unique American eatery you think he should visit, send him an e-mail at travel@usatoday.com. Some of the venues reviewed by this column provided complimentary services.


The good news: this is not a post about bad golfers like myself who routinely score quadruple bogeys on any particular hole (i.e. take 4 more strokes than par).
However, it sort of is, as three of the "Bogeys" are related to golf. Let me explain: five pictures will be worth a total of the order of a thousand words.
Depending on where you are, a drive of around fifty miles can provide much different time frames.
Much of that can vary, depending on speed limits and traffic. This post, in some sense, reflects on why I am sometimes hesitant to drive to Milwaukee or Chicago for JUST an evening.