# STAIC regularization for spatio-temporal image reconstruction

Deepak G Skariah
Department of Electrical Engineering
IISc
Bengaluru, 560012
<EMAIL_ADDRESS>
&Muthuvel Arigovindan
Department of Electrical Engineering
IISc
Bengaluru,560012
<EMAIL_ADDRESS>

###### Abstract

We propose a regularization-based image restoration scheme for 2D images
recorded over time (2D+t). We design an infimal convolution-based
regularization function which we call spatio-temporal Adaptive Infimal
Convolution (STAIC) regularization. We formulate the infimal convolution in
the form of an additive decomposition of the 2D+t image such that the extent
of spatial and temporal smoothing is controlled in a spatially and temporally
varying manner. This makes the regularization adaptable to the local
characteristics of the motion leading to an improved ability to handle noise.
We also develop a minimization method for image reconstruction by using the
proposed form of regularization. We demonstrate the effectiveness of the
proposed regularization using TIRF images recorded over time and compare with
some selected existing regularizations.

_Keywords_ Spatio-Temporal $\cdot$ Regularization $\cdot$ Restoration

##  1 Introduction

Image restoration is an inverse problem Katsaggelos (1989) where a higher
quality image estimate is generated from a corrupted observation by exploiting
knowledge of image statistics. Regularization methods Engl et al. (1996)
constitute an important category among methods for biomedical image
restoration Rangayyan (2004). Imaging modalities where regularization methods
have been successful in restoration of corrupted images include MRI Ramani et
al. (2012); Viswanath et al. (2020), Widefield Microscopy Arigovindan (2013);
Li et al. (2018) and Total Internal Reflection (TIRF) Microscopy Fan et al.
(2019) among others. Regularization involves formulating the required restored
image ${g}_{opt}$ as a minimizer of a cost function involving the observed
image ${m}$. The cost is formulated as a sum of a data fitting term $G({g,m})$
and a regularization term $R({g})$ resulting in the following minimization
problem:

${{g}}_{opt}=\operatornamewithlimits{argmin}\limits_{{g}}{\;\;G({g,m})+\lambda
R({g})}$ (1)

where $\lambda$ is the regularization parameter which controls the relative
weighting of regularization term against the data fitting term in the overall
cost. The choice of data fitting term $G(g,m)$ is dependent on the image
formation forward-model. The regularization term enforces any prior we have
about the class of images we are trying to restore often in the form of some
image regularity condition.

Regularization functional $R(g)$ is primarily designed based on study of image
statistics, while some are data driven designs as observed in learning based
methods. Image statistics often depend on imaging modality and the type of
objects being imaged. One of the earliest successful regularization function
for restoration of natural images is Total Variation. Total Variation, more
specifically, first-order total variation (TV-1) is defined as the sum of
absolute values of image derivatives along image directions. It works on the
principle that discrete natural images have a limited amount of variation
which is captured via the first-order derivatives along the two image
dimensions. For bio-medical imaging modalities often a better choice is
second-order Total Variation (TV-2) where the sum of second-order derivatives
along the two image dimensions are employed instead of the first order ones.
Replacing first-order derivatives with the second-order ones leads to more
natural intensity variation in the reconstructed images Lefkimmiatis et al.
(2011). A generalization of TV-2 known as Hessian-Schatten regularization (HS)
was successfully employed for the restoration of biological images
Lefkimmiatis et al. (2013). HS regularization employed Schatten norm on the
image Hessian. HS regularization and its variants have been successfully
applied in a wide variety of imaging forward models including natural images
and biological images.

Further developments led to more complex forms of regularization. Among these,
two forms are prominent in the literature: (i) sum of norms regularization
Lindsten et al. (2011) and (ii) infimal decomposition-based regularization
Holler and Kunisch (2014). Both approaches multiple types of derivatives, but
they are fundamentally different in the way the multiple order are combined,
which leads to a significant difference in the performance. In the first
approach, $R(g)$ takes the form:

$R(g)=\lambda_{1}R_{1}(g)+\lambda_{2}R_{2}(g)$ (2)

where $R_{1}$ and $R_{2}$ are the functions capturing different priors based
on different types of derivatives. This will produce a reconstruction that
balances regularity assumptions enforced by $R_{1}$ and $R_{2}$ [COTV].
Unfortunately, this balancing is global in the sense that we can control only
the overall agreement of the solution to the priors captured by the terms
$R_{1}$ and $R_{2}$. Real images have non-stationary statistics, and hence it
will be more advantageous to combine different terms in a spatially adaptive
manner. This can be accomplished by, for example, making the weights
$\lambda_{1}$ and $\lambda_{2}$ spatially varying, i.e., by making it vary
from pixel to pixel. In this case, the weights $\lambda_{1}$ and $\lambda_{2}$
themselves become images and determining these images along with the required
image becomes a challenge. This problem can also be handled by an auxiliary
regularization Viswanath et al. (2020). The main problem in this approach is
that the overall cost becomes non-convex.

The second approach is based on infimal convolution Holler and Kunisch (2014).
Here the regularization functional $R(g)$ itself is defined by a minimization
problem. Infimal convolution based regularization functional $R(g)$ takes the
following form:

$R(g)=\operatornamewithlimits{min}\limits_{v_{1}+v_{2}=g}{\;\;R_{1}(v_{1})+R_{2}(v_{2})}$
(3)

Here $v_{1}$ and $v_{2}$ are auxiliary variables satisfying the condition that
$v_{1}+v_{2}=g$. Note that here too, $R(g)$ is defined as a sum of two
functions $R_{1}$ and $R_{2}$. But this time, $R_{1}$ acts on a variable
$v_{1}$ and $R_{2}$ acts on a variable $v_{2}$ where $v_{1}+v_{2}=g$. In the
presence of data-fitting cost that is defined on the sum $g=v_{1}+v_{2}$,
minimization with respect to $v_{1}$ and $v_{2}$ leads to spatial adaptivity;
the relative influence of $R_{1}$ and $R_{2}$ is determined adaptively in a
spatially varying manner depending on the spatial structure of the input image
$m$. Hence we get an effect that is similar to the method of Viswanath et al.
(2020) with an advantage that the cost here is convex. An additional advantage
is that this formulation has a higher extent of adaptivity: it is possible to
get a solution in which either $R_{1}$ or $R_{2}$ has no influence at all in
some spatial locations depending on the local structure in $m$ such that
overall distortion is minimized. Total Generalized Variation (TGV) Guo et al.
(2014); Bredies et al. (2010) is the most popular regularization designed
using this concept and is used successfully in restoration of MRI images.

More recently learning based methods have been applied successfully to
restoration of images from a wide variety of imaging modalities. These belong
to the class of deep learning based methods where the image prior is
represented by using Convolution Neural networks. These methods use end to end
learning based methods or use of CNN networks as a prior for regularization-
based image restoration Ulyanov et al. (2018). In microscopy it has been
applied in restoration of 2D images. It has also been applied to other image
analysis problems in microscopy as well. Xing et al. (2017); von Chamier et
al. (2021); Liu et al. (2021). The disadvantages of using learning in
microscopy image analysis has been discussed in Hoffman et al. (2021).

A related category of inverse problems is restoration of spatio-temporal 2D
images observed over time. Even though a spatio-temporal image is simply a
sequence of 2D images, the following are the disadvantages of applying 2D
image restoration for each time-point independently: (i) temporal correlation
is ignored leading to poorer quality of restored images; (ii) temporal
continuity is lost in the restored images Holler and Kunisch (2014). Due to
the inherent need for imposing regularity in both time and space, spatio-
temporal image restoration demands for the employment of sum of norms or
infimal decomposition based construction of the regularization functional. In
the sum-of-norms based construction, a weighted sum of two terms as given in
the equation (2) is typically used [ref]: (i) the first term applies spatial
regularization on each time point of the image sequence; (ii) the second term
enforces temporal smoothing for each pixel independently. We will call such
construction the combined spatio-temporal regularization (CST). In infimal
decomposition approach, a form as given in the equation Equation 3 is used.
Here both the terms use spatial and temporal derivatives, but, they differ by
the relative weight applied to the temporal derivatives. This method was
called the Infimal Convolution Total Variation (ICTV), and we will refer this
methods as ICTV-2DT. ICTV-2DT regularization has been successfully applied in
biomedical imaging modalities such as dynamic MRISchloegl et al. (2017) and
SPECT imaging Zhang et al. (2018). Among other regularization based
approaches, low rank sparse decomposition based model for dynamic MRI
reconstruction Trémoulhéac et al. (2014) is the most well known. In MRI more
recently deep learning based methods have been proposed. For restoration of
real world video signals deep neural network based algorithms have been
proposed. More recently RNN networks had been proposed for dynamic MRI
reconstruction Chen et al. (2022).

In this work, we develop a novel spatio-temporal regularization based on
infimal convolution employing two terms as done in ICTV-2DT. However, our
regularization significantly differs from ICTV-2DT in how spatial and temporal
derivatives as distributed among the two terms. We construct the form of each
term to make the overall regularization suitable for restoration of 2D
fluorescence images recorded as a function of time. We call the novel
regularization the Spatio-Temporal Adaptive Infimal Convolution (STAIC)
regularization. We also develop a computational algorithm for dynamic image
restoration by using the STAIC regularization, and demonstrates the usefulness
of our method using simulation experiments. The paper is organized as follows.
In section Section 2, we introduce the notations and mathematical
preliminaries. In section Section 3, we describe the important existing
regularization approaches to deal with spatio-temporal signal restoration
problems. In section Section 4, we introduce the proposed STAIC
regularization. In section Section 5, we present the STAIC regularized
reconstruction as a convex optimization problem, and describe how we solve the
optimization problem using the ADMM approach. In section Section 6, we solve
the ADMM sub-problems that arise out of the reconstruction problem in more
detail. Finally in Section 7, we validate the utility of STAIC regularization.

##  2 Notations and mathematical preliminaries

  1. 1.

Images are represented by lower case English alphabets. For example $g$.

  2. 2.

In the discussion we will use the idea of vector valued images often refereed
to as vector images. Vector images are discrete 2D arrays where each pixel
location has a vector quantity. It is denoted by lower-case bold-faced letter
with a bold-faced lower-case letter as an argument. For example,
${\mathbf{v}}({\mathbf{r}})$ is a vector image with
${\mathbf{r}}=[x\;y]^{\top}$ representing a 2D pixel location. Depending on
the context, the symbol denoting the pixel location may be omitted.

  3. 3.

For a vector image, ${\mathbf{v}}({\mathbf{r}})$, $\|{\mathbf{v}}\|_{1,2}$
denotes
$\|{\mathbf{v}}\|_{1,2}=\sum_{\mathbf{r}}\|{\mathbf{v}}({\mathbf{r}})\|_{2}$.
It is the sum of pixel-wise $l_{2}$ norms, where $\sum_{\mathbf{r}}$ denotes
the sum across pixel indices. The bound of sum is the first to last pixel
location in this notation. The norm $\|{\mathbf{v}}\|_{1,2}$ is a composition
of norms and is often refereed to as a mixed norm.

  4. 4.

In a scalar image having multiple frames, we use the subscript notation to
refer to a particular frame number. Example $g_{i}$ refers to frame number $i$
of the spatio-temporal image $g$.

  5. 5.

Index $\mathbf{r}$ is used to refer to a spatio-temporal image (2DT). Index
$\mathbf{k}$ is used to refer to a 2D image. Let $\delta(\mathbf{k})$ and
$\delta(\mathbf{r})$ represent 2D and 3D Kronecker delta respectively.

  6. 6.

$\|\mathbf{g}\|_{1,\kappa}=\sum_{\mathbf{r}}\|\mathbf{g}(\mathbf{r})\|_{\kappa}$
with the definition of
$\|\mathbf{y}\|_{\kappa}=\sqrt{\kappa^{2}(y_{1}^{2}+y_{2}^{2})+y_{3}^{2}}$.
Further we have that
$\|\mathbf{g}\|_{1,1/\kappa}=\sum_{\mathbf{r}}\|\mathbf{g}(\mathbf{r})\|_{1/\kappa}$
and
$\|\mathbf{x}\|_{1/\kappa}=\sqrt{y_{1}^{2}+y_{2}^{2}+\kappa^{2}y_{3}^{2}}$.
Here $\kappa$ is a parameter.

  7. 7.

$\bar{d}_{xx}(\mathbf{k}),\bar{d}_{xy}(\mathbf{k}),\bar{d}_{yy}(\mathbf{k})$
are discrete filters implementing 2D second order derivatives. In addition
$d_{xx}(\mathbf{r}),d_{xy}(\mathbf{r}),d_{xt}(\mathbf{r})$,$d_{yy}(\mathbf{r}),d_{yt}(\mathbf{r}),d_{tt}(\mathbf{r})$
are discrete filters implementing the 3D second order derivatives.

  8. 8.

$M(\mathbf{r})$ is the 3D Hessian operator encapuslating all the second order
derivative operators where the directions are $(x,y,t)$

##  3 Related methods in detail

Before discussing related work, it is important to understand the expected
behavior from an ideal spatio-temporal regularization function. In regions
with motion, regularization function must be able to discourage temporal
smoothing. In regions without motion, temporal smoothing must be promoted. All
existing approaches to regularization for this inverse problem is trying to
achieve this ideal behavior in some way or the other. This is possible only if
the regularization function $R(g)$ is powerful enough to differentiate between
regions with motion and regions without motion. This points towards the need
to incorporate both spatial and temporal derivatives in the design of
regularization function. The presence of temporal derivatives also ensures
that such functions operate differently to classical image regularization
methods when employed in the resultant optimization problem.

Early methods that consider temporal derivatives together with spatial
derivatives constructed the overall regularization as a simple weighted sum of
functionals constructed using spatial and temporal derivative independently.
This sum of norm based spatio-temporal regularization which we called Combined
Spatio-Temporal (CST) is defined as :

$R_{CST}(g)=\|\mathcal{M}*{g}\|_{1,2}+\|d_{tt}*{g}\|_{1,2}$ (4)

where $\mathcal{M}$ is the second order Hessian operator in 3D and $D_{tt}$ is
second order derivative along time dimension. In restoration of spatio-
temporal signals, $R_{CST}$ is known not to satisfy the desirable properties
of an ideal spatio-temporal regularization stated above. One disadvantage of
the above formulation is that it is blind to the fact that characteristics of
motion varies spatially and temporally, and there is no spatially and temporal
varying trade-off between spatial and temporal derivatives. In regions with
motion both components will have nontrivial contribution to overall cost
leading temporal smoothing which is undesirable. In regions without motion
spatial smoothing occurs in addition to temporal smoothing which is again
undesirable.This shortcoming prompted design of spatio-temporal regularization
employing the infimal decomposition approach as done in the ICTV-2DT method
introduced before. Let $g(\mathbf{r})$ be the discrete candidate spatio-
temporal (2D+time) image, where $\mathbf{r}=[x\;y\;t]^{\top}$ is the discrete
3D pixel index. Assume that image $g(\mathbf{r})$ is of dimensions
$n_{1}\times n_{2}\times n_{F}$ where $n_{F}$ is the number of frames (2D
images). Hence each frame is of dimension $n_{1}\times n_{2}$. The
regularization is constructed in the form of minimization as given below:

$R_{ICTV}(g)=\operatornamewithlimits{min}\limits_{v}\|{\cal
S}_{\kappa_{1}}\nabla(g-v)\|_{1,\kappa_{1}}+\|{\cal
S}_{\kappa_{2}}\nabla(v)\|_{1,\kappa_{2}},$ (5)

where ${\cal S}_{\kappa}=diag(\kappa,\kappa,1)$ and $\nabla$ is the 3D
gradient operator. Here $\kappa_{1}$ and $\kappa_{2}$ are user parameters that
determines relative weight between spatial and temporal derivatives. It may be
noted that the method requires tuning for $\kappa_{1}$ and $\kappa_{2}$ in
addition to the tuning of the overall weight $\lambda$. This form has the
advantage over CST in the same way as the form introduced in the equation (2)
has an advantage over the form in the equation (3). In other words, combining
two different types of sub-functional in the form of infimal convolution is
always advantageous over a simple weighted. However, ICTV-2DT still has a
disadvantage. This can be understood by analyzing its effect on regions with
motion and without motion. In regions with motion, it will be desirable to
eliminate temporal smoothing completely, but, ICTV-2DT does not do so.
Similarly, in regions without motion, it will be desirable to eliminate
spatial smoothing completely, but, ICTV-2DT does not do so. The reason is that
both sub-functional have both temporal and spatial derivatives and, they
differ only by the relative weight between spatial and temporal derivatives.
Further, it uses first order derivatives which is again undesirable.

##  4 Spatio-Temporal Adaptive Infimal Convolution (STAIC) Regularization for
2DT image restoration

We propose STAIC regularization for spatio-temporal images by removing the
disadvantages of ICTV-2DT while strengthening the design to cater to unique
features of 2DT fluorescence images. Let $g_{i}(\mathbf{k})$ represent the
$i^{th}$ frame of the spatio-temporal image $g(\mathbf{r})$ where
$\mathbf{k}=[x\;y]^{\top}$ is the discrete $2D$ pixel index. Consider two
linear operators $A(\mathbf{k})$ and $M(\mathbf{r})$ defined as follows :

$A(\mathbf{k})=\begin{bmatrix}\bar{d}_{xx}(\mathbf{k})\\\
\bar{d}_{xy}(\mathbf{k})\\\ \bar{d}_{yx}(\mathbf{k})\\\
\bar{d}_{yy}(\mathbf{k})\\\
\delta(\mathbf{k})\end{bmatrix},M(\mathbf{r})=\begin{bmatrix}d_{xx}(\mathbf{r})&d_{xy}(\mathbf{r})&d_{xt}(\mathbf{r})\\\
d_{yx}(\mathbf{r})&d_{yy}(\mathbf{r})&d_{yt}(\mathbf{r})\\\
d_{tx}(\mathbf{r})&d_{ty}(\mathbf{r})&d_{tt}(\mathbf{r})\end{bmatrix}$

Next, we define the following components, which act on 2DT images:

$\displaystyle\bar{R}_{1}(u_{i})$ $\displaystyle=\|A*u_{i}\|_{1,2}$ (6)
$\displaystyle R_{2}(v)$ $\displaystyle=\|(M*v)\|_{1,F}$ (7)

Here $u$ and $v$ are 2DT images, and $u_{i}$ is the $i^{th}$ 2D time frame of
of the 2DT image $u$. Note that $\|A*u_{i}\|_{1,2}=\sum_{\bf
k}\|(A*u_{i})({\bf k})\|_{2}=\sum_{\bf
k}\sqrt{((\bar{d}_{xx}*u_{i})(\mathbf{k}))^{2}+2((\bar{d}_{xy}*u_{i})(\mathbf{k}))^{2}+((\bar{d}_{yy}*u_{i})(\mathbf{k}))^{2}+(u_{i}({\bf
k}))^{2}}$. This is a modified 2D second order TV applied on $u_{i}$, where
the modification is simply the addition of the intensity term. Next, we note
that $R_{2}(v)$ computes the 3D second-order TV applied on the 2DT sequence
$v$.

The proposed STAIC regularization is defined as follows :

$S(g,\alpha_{s},\alpha_{t})=\operatornamewithlimits{min}\limits_{v}{\;\;\alpha_{s}R_{1}(g,v)+\alpha_{t}R_{2}(v)}$

where $R_{1}(g,v)=\sum_{i=1}^{n_{F}}\bar{R}_{1}(g_{i}-v_{i})$. The user
defined parameters $\alpha_{s}$, $\alpha_{t}$ helps control the relative
strength of $R_{1}(\cdot,\cdot)$ and $R_{2}(\cdot)$ in the regularization
function $S(g,\alpha_{s},\alpha_{t})$.

Despite the role of $v$ as an auxiliary variable in the definition of
regularization, the component $v$ is regularized to have finite amount of
variations in its 3D structure. This is achieved by constraining the Frobenius
norm of 3D second order derivatives. This reflects our desire to ensure that
$v$ also resemble a realistic spatio-temporal signal. In the first term
$R_{1}(g,v)$, we put a constraint on spatial regularity of the difference
signal by incorporating an intensity term in addition to the spatial
derivative term. This combination of derivatives and pixel intensity promotes
the presence of blob-like structure in the restored image. Such image features
are often encountered in fluorescence microscopy images which is desirable in
our setting due to the TIRF microscopy origin of our signal of interest.

To understand the effect of the two terms on the overall regularization
function, it is important to consider regions with motion and without motion
separately.

  * •

In regions with motion, $v$ will tend to be zero, and so only $R_{1}$ comes
into play; this means that there will be negligible temporal smoothing, which
will reduce motion blurring

  * •

In regions without motion, $v$ will tend to be close to $g$ so that only
$R_{2}$ come into play; since $R_{2}$ imposes 3D smoothing, this lead to
robust noise removal

It may be noted that the design of $R_{1}$ in STAIC removes an important
drawback of ICTV-2DT in regions with motion. STAIC due to the presence of only
spatial derivatives in its $R_{1}$ term discourages temporal smoothing
completely in regions with motion. This is an important feature as it avoids
motion blurring in the final image estimate. In addition, a crucial design
change is the presence of a pure intensity term $g-v$ in the $R_{1}$ term.
This is done due the following reasons :

  * •

We target our method of restoration fluorescence images; fluorescence
microscopy images often have images exhibiting a blob like pattern which is
promoted by the presence of intensity term

  * •

In addition the intensity term help promote sparsity of regions with motion in
turn helping remove micro motions which are often the artifacts caused by
Poisson nature of noise. This helps in recovering all large and small motions
while removing any spurious motion in the signal.

In addition, STAIC formulations also do not employ any parameter to control
the relative strength of spatial and temporal derivatives by virtue of its
design leaving us one less parameter to tune. In summary,(i) the redesigned
design of purely spatial term $R1$ (ii) intensity augmentation of $R1$ helps
us achieve the goal of designing a regularization scheme that has close to
ideal spatio-temporal behavior in regions with motion and without it with
limited number of parameters. This ensures that the reconstruction is closer
to the ground truth in terms of both spatial and temporal features.

##  5 Signal Reconstruction Using STAIC regularization

###  5.1 The Cost Function

Let $h$ denote the 2D impulse response of the microscope. The noise model in
fluorescence images is modeled to be mixed Poisson-Gaussian where a signal
corrupted by Poisson noise is subject to an additive Gaussian noise. The
generation of measured image $m$ is therefore modeled as follows:

$m=\mathcal{P}(h*g)+\eta$

where $\mathcal{P}$ is the operator representing Poisson process, $\eta$
represents the additive Gaussian noise, $m$ represents the observed noisy
blurred image and $g$ represents the ground truth image. The data fitting term
is often chosen to be the maximum likelihood estimate of the noise model. But
we choose the least squares data fitting term due to ease of optimization
Arigovindan et al. (2013). The data fitting term is

$G(h,g,{m})=\frac{1}{2}\sum_{i=1}^{n_{F}}\|(h*g_{i})-m_{i}\|_{F}^{2}$ (8)

It may be noted that the blurring happens frame wise as we are considering a
2D signal observed over time. The STAIC regularized image reconstruction
optimization problem can be expressed as

$g_{opt}=\operatornamewithlimits{argmin}\limits_{g}{\;\;G(h,g,m)+S(g,\alpha_{s},\alpha_{t})}$
(9)

By observing that the regularization functional $S(g,\alpha_{s},\alpha_{t})$
is defined as an optimization problem over variable $v$, we can express the
reconstruction problem as:

$\displaystyle(g_{opt},v_{opt})=\operatornamewithlimits{argmin}\limits_{g,v}{\;\;G(h,g,m)+\bar{S}(g,v,\alpha_{s},\alpha_{t})}$
(10) $\displaystyle\text{ where
}\bar{S}(g,v,\alpha_{s},\alpha_{t})=\alpha_{s}R_{1}(g,v)+\alpha_{t}R_{2}(v)$

resulting in an optimization problem in 2 variables $g$ and $v$. This is
possible because, the variable $v$ do not appear in the data fitting term
enabling it to be moved to the outer optimization step. We introduce a
variable (vector image)
$\mathbf{f(\mathbf{r})}=[g(\mathbf{r})\;\;v(\mathbf{r})]^{\top}$ concatenating
the two variables in the above optimization problem to enable a simpler
restatement of the cost of optimization. To this end we define the following
new linear operators :

$\displaystyle{\mathbf{h}}$ $\displaystyle=[h(\mathbf{k})\;\;0],$
$\displaystyle{\mathbf{T}}_{s}$
$\displaystyle=\left[\begin{array}[]{cc}d_{xx}(\mathbf{k})&0\\\
d_{xy}(\mathbf{k})&0\\\ d_{yx}(\mathbf{k})&0\\\ d_{yy}(\mathbf{k})&0\\\
\delta(\mathbf{k})&0\\\ 0&d_{xx}(\mathbf{k})\\\ 0&d_{xy}(\mathbf{k})\\\
0&d_{yx}(\mathbf{k})\\\ 0&d_{yy}(\mathbf{k})\\\
0&\delta(\mathbf{k})\end{array}\right]\text{ and
}\;\;{\mathbf{T}}_{t}=\left[\begin{array}[]{cc}0&d_{xx}(\mathbf{r})\\\
0&d_{yy}(\mathbf{r})\\\ 0&d_{xy}(\mathbf{r})\\\ 0&d_{yx}(\mathbf{r})\\\
0&d_{xt}(\mathbf{r})\\\ 0&d_{tx}(\mathbf{r})\\\ 0&d_{yt}(\mathbf{r})\\\
0&d_{ty}(\mathbf{r})\\\ 0&d_{tt}(\mathbf{r})\end{array}\right]$

$\displaystyle{A}_{s}=\left[\begin{array}[]{cccccccccc}1/\sqrt{2}&0&0&0&0&-1/\sqrt{2}&0&0&0&0\\\
0&1/\sqrt{2}&0&0&0&0&-1/\sqrt{2}&0&0&0\\\
0&0&1/\sqrt{2}&0&0&0&0&-1/\sqrt{2}&0&0\\\
0&0&0&1/\sqrt{2}&0&0&0&0&-1/\sqrt{2}&0\\\
0&0&0&0&1/\sqrt{2}&0&0&0&0&-1/\sqrt{2}\end{array}\right]$ (16)

The optimization problem in (10) can now be reformulated in terms of
$\mathbf{f}$ only as :

$\displaystyle J(\mathbf{f},\alpha_{s},\alpha_{t})=$
$\displaystyle\frac{1}{2}\sum_{i=1}^{n_{F}}\|\mathbf{h}*\mathbf{f}_{i}-m_{i}\|_{F}^{2}+\sqrt{2}\alpha_{s}\sum_{i=1}^{n_{F}}\sum_{\mathbf{k}}\|A_{s}\big{(}(\mathbf{T_{s}}*\mathbf{f}_{i})(\mathbf{k})\big{)}\|_{2}+\alpha_{t}\sum_{\mathbf{r}}\|(\mathbf{T_{t}}*\mathbf{f})(\mathbf{r})\|_{2}+\mathcal{B_{C}}(\mathbf{f})$
(17)

Here, the additional term $\mathcal{B_{C}}(\cdot)$ is added to impose bound
constraint on $\mathbf{f}$ which restrict the range of pixel values
$\mathbf{f}$ can take. This helps impose commonly used non-negativity
condition on pixel values. Here $\mathcal{B_{C}}(\cdot)$ is defined as follows

$\mathcal{B_{C}}(\mathbf{f})=\begin{cases}\text{0}&\quad\text{if
}\mathbf{f}\in\mathcal{C}\\\ \infty&\quad\text{otherwise}\end{cases}$ (18)

Here $\mathcal{C}=\\{\mathbf{f}\;|\;lb\leq\mathbf{f}(\mathbf{r})\leq
ub\;\;\forall\mathbf{r}\\}$ where $lb$ and $ub$ represents the smallest and
highest permitted pixel values. The new image restoration optimization problem
over variable $\mathbf{f}$ now takes the form :

$\mathbf{f}_{opt}=\operatornamewithlimits{argmin}\limits_{\mathbf{f}}{\;\;J(\mathbf{f},\alpha_{s},\alpha_{t})}$
(19)

###  5.2 Proposed ADMM Method

The cost $J(\mathbf{f},\alpha_{s},\alpha_{t})$ to be minimized is a convex
function in variable $\mathbf{f}$. We propose to design an algorithm to
minimize the cost by employing the Alternating Direction Method of Multipliers
(ADMM) Boyd et al. (2011) algorithmic framework. ADMM framework is most suited
for minimizing convex cost under linear equality constraints. Our original
optimization problem in (19) can be reformulated to get an equivalent linearly
constrained convex problem. This is done by by introduction of new variables
$w_{m},\mathbf{w}_{s},\mathbf{w}_{t}$ and $\mathbf{w}_{b}$ and as a result,
the updated problem becomes:

$\displaystyle\operatornamewithlimits{min}\limits_{\mathbf{f},w_{m},\mathbf{w}_{s},\mathbf{w}_{t},\mathbf{w}_{b}}{\;\;\frac{1}{2}\sum_{i=1}^{n_{F}}\|w_{m_{i}}-m_{i}\|_{F}^{2}}+\sqrt{2}\alpha_{s}\sum_{i=1}^{n_{F}}\sum_{\mathbf{k}}\|A_{s}\big{(}\mathbf{w}_{s_{i}}(\mathbf{k})\big{)}\|_{2}+\alpha_{t}\sum_{\mathbf{r}}\|\mathbf{w}_{t}(\mathbf{r})\|_{2}+\mathcal{B_{C}}(\mathbf{w}_{b})$
$\displaystyle\text{ subject to }\quad\quad
w_{m_{i}}=\mathbf{h}*\mathbf{f}_{i},\;\;\mathbf{w}_{s_{i}}=\mathbf{T}_{s}*\mathbf{f}_{i},\;\;\mathbf{w}_{t}=\mathbf{T}_{t}*\mathbf{f},\;\;\mathbf{w}_{b}=\mathbf{f}$

It may be noted that we have converted the unconstrained problem in (19) to a
constrained optimization form where all constraints are linear equality
constraints. To allow a simpler algorithm statement, we introduce a combined
operator $\mathbf{T}$ and a combined vector $\mathbf{w}$ defined as follows :

$\displaystyle\mathbf{T}=\begin{bmatrix}\mathbf{h}(\mathbf{k})\\\
\mathbf{T}_{s}(\mathbf{k})\\\ \mathbf{T}_{t}(\mathbf{r})\\\
\mathbf{e}(\mathbf{r})\end{bmatrix},\quad\mathbf{w}=\begin{bmatrix}w_{m}\\\
\mathbf{w}_{s}\\\ \mathbf{w}_{t}\\\ \mathbf{w}_{b}\end{bmatrix}$ (20)

where $\mathbf{e}(\mathbf{r})=[\delta(\mathbf{r})\;\;\delta(\mathbf{r})]$.
Under this definition, $\mathbf{e}*\mathbf{f}=\mathbf{f}$. The linearly
constrained problem can now be stated in a compact form as :

$\displaystyle(\mathbf{f}^{*},\mathbf{w}^{*})=\operatornamewithlimits{argmin}\limits_{\mathbf{f},\mathbf{w}}{\;\;R(\mathbf{w},\alpha_{s},\alpha_{t})}$
(21) $\displaystyle\text{ subject to
}\quad\quad\mathbf{T}*\mathbf{f}=\mathbf{w}$

where
$R(\mathbf{w},\alpha_{s},\alpha_{t})=G(w_{m},m)+\sqrt{2}\alpha_{s}\sum_{i=1}^{n_{F}}\sum_{\mathbf{k}}\|{A}_{s}(\mathbf{w}_{s_{i}}(\mathbf{k}))\|_{2}+\alpha_{t}\sum_{\mathbf{r}}\|\mathbf{w}_{t}(\mathbf{r})\|_{2}+\mathcal{B_{C}}(\mathbf{w}_{b})$
and $G(w_{m},m)=\frac{1}{2}\sum_{i=1}^{n_{F}}\|w_{m_{i}}-m_{i}\|_{F}^{2}$. The
next step in ADMM framework is to construct the Augmented Lagrangian Boyd et
al. (2011)
$\mathcal{L}(\mathbf{f},\mathbf{w},{\boldsymbol{\beta}},\alpha_{s},\alpha_{t})$
of the above linearly constrained cost as:

$\mathcal{L}(\mathbf{f},\mathbf{w},\boldsymbol{\beta},\alpha_{s},\alpha_{t})=R(\mathbf{w},\alpha_{s},\alpha_{t})+\langle{\boldsymbol{\beta}},\mathbf{T}*\mathbf{f}-\mathbf{w}\rangle+\frac{\rho}{2}\|\mathbf{T}*\mathbf{f}-\mathbf{w}\|_{2}^{2}$
(22)

where ${\boldsymbol{\beta}}$ is the Lagrangian multiplier and $\rho$ is a user
supplied ADMM parameter. Here the dimensions of $\boldsymbol{\beta}$ is same
as that of $\mathbf{w}$. Finally, ADMM algorithm involves collection of
alternative minimization of sub-problem with respect to variables $\mathbf{f}$
and $\mathbf{w}$ followed by an update step involving variable
$\boldsymbol{\beta}$. Assume that
$\mathbf{f}^{(k)},\mathbf{w}^{(k)},{\boldsymbol{\beta}}^{(k)}$ are the current
estimates , the ADMM algorithm involves the following steps

$\displaystyle\mathbf{w}^{(k+1)}$
$\displaystyle=\operatornamewithlimits{argmin}\limits_{\mathbf{f}}{{L}(\mathbf{f}^{(k)},\mathbf{w},{\boldsymbol{\beta}}^{(k)},\alpha_{s},\alpha_{t})}$
(23) $\displaystyle\mathbf{f}^{(k+1)}$
$\displaystyle=\operatornamewithlimits{argmin}\limits_{\mathbf{f}}{L(\mathbf{f},\mathbf{w}^{(k+1)},{\boldsymbol{\beta}}^{(k)},\alpha_{s},\alpha_{t})}$
(24) $\displaystyle\text{ and }\;\;{\boldsymbol{\beta}}^{(k+1)}$
$\displaystyle={\boldsymbol{\beta}}^{(k)}+\rho\big{(}\mathbf{T}*\mathbf{f}^{(k+1)}-\mathbf{w}^{(k+1)}\big{)}$
(25)

The first two equation involves solving two optimization problems over
variables $\mathbf{f}$ and $\mathbf{w}$ respectively which is discussed next.

##  6 Solving the Sub-problems of ADMM Algorithm

We will now discuss how the sub-problems in (23), and (24) are solved to
obtain the intermediate variables $\mathbf{w}^{(k+1)}$ and
$\mathbf{f}^{(k+1)}$ that appear in the ADMM iterative scheme.

###  6.1 The w problem

The sub-problem in (23) can be equivalently simplified to the following form

$\displaystyle
L_{w,k}(\mathbf{w},\alpha_{s},\alpha_{t})=R(\mathbf{w},\alpha_{s},\alpha_{t})+\frac{\rho}{2}\|\mathbf{w}-\bar{\mathbf{x}}^{(k)}\|_{2}^{2}$
(26) $\displaystyle\quad\text{ where
}\quad\bar{\mathbf{x}}^{(k)}=\mathbf{T}*{\mathbf{f}^{(k)}}+\frac{1}{\rho}{\boldsymbol{\beta}^{(k)}}$
(27)

For cleaner presentation of sub-problems, we introduce the notation
$\mathbf{x}=\bar{\mathbf{x}}^{(k)}$ and $\hat{\mathbf{w}}=\mathbf{w}^{(k+1)}$.
Since $\mathbf{w}$ is made up of sub vectors
$w_{m},\mathbf{w}_{b},\mathbf{w}_{t},\mathbf{w}_{s}$, we separate the above
problem into sub-problems involving constituent variables.

$w_{m}$-prob.:
$\displaystyle\hat{w}_{m}=\underset{w_{m}}{\operatorname{argmin}}\underbrace{\frac{\rho}{2}\|x_{m}-w_{m}\|_{2,2}^{2}+G({w}_{m},{m})}_{\bar{L}_{m}(w_{m},x_{m})}$
(28) $\mathbf{w}_{b}$-prob.:
$\displaystyle\hat{\mathbf{w}}_{b}=\underset{\mathbf{w}_{b}}{\operatorname{argmin}}\;\;\underbrace{\frac{\rho}{2}\|\mathbf{x}_{b}-\mathbf{w}_{b}\|_{2,2}^{2}+{\mathcal{B_{C}}}(\mathbf{w}_{b})}_{\bar{L}_{b}(\mathbf{w}_{b},\mathbf{x}_{b})}$
(29) ${\mathbf{w}}_{t}$-prob.:
$\displaystyle\hat{\mathbf{w}}_{t}=\underset{{\mathbf{w}}_{t}}{\operatorname{argmin}}\;\;\underbrace{\frac{\rho}{2}\|{\mathbf{x}}_{t}-{\mathbf{w}}_{t}\|_{2,2}^{2}+\alpha_{t}\|\mathbf{w}_{t}\|_{1,2}}_{\bar{L}_{t}({\mathbf{w}}_{t},{\mathbf{x}}_{t},\alpha_{t})}$
(30) ${\mathbf{w}}_{s}$-prob.:
$\displaystyle\hat{\mathbf{w}}_{s}=\underset{{\mathbf{w}}_{s}}{\operatorname{argmin}}\;\;\sum_{i=1}^{n_{F}}\underbrace{\frac{\rho}{2}\|{\mathbf{x}}_{s_{i}}-{\mathbf{w}}_{s_{i}}\|_{2,2}^{2}+\sqrt{2}\alpha_{s}\sum_{\mathbf{k}}\|{A}_{s}(\mathbf{w}_{s_{i}}(\mathbf{k}))\|_{2}}_{\bar{L}_{s_{i}}({\mathbf{w}}_{s_{i}},{\mathbf{x}}_{s_{i}},\alpha_{s})}$
(31)

####  6.1.1 Decomposing problems pixel-wise

The cost $\bar{L}_{m}(w_{m},x_{m})$ is separable across pixels as shown below:

$\displaystyle\bar{L}_{m}(w_{m},x_{m})=$
$\displaystyle\frac{\rho}{2}\|x_{m}-w_{m}\|_{2,2}^{2}+G({w}_{m},m)$ (32)
$\displaystyle=$
$\displaystyle\frac{\rho}{2}\|x_{m}-w_{m}\|_{2,2}^{2}+\frac{1}{2}\|w_{m}-m\|_{2}^{2}$
(33) $\displaystyle=$
$\displaystyle\sum_{\mathbf{r}}\underbrace{\frac{\rho}{2}(x_{m}(\mathbf{r})-w_{m}(\mathbf{r}))^{2}+\frac{1}{2}(w_{m}(\mathbf{r})-m(\mathbf{r}))^{2}}_{{L}_{m}(w_{m}(\mathbf{r}),x_{m}(\mathbf{r}))}$
(34)

Hence the pixel wise cost ${L}_{m}(w_{m}(\mathbf{r}),x_{m}(\mathbf{r}))$ is
given by:

${L}_{m}(w_{m}(\mathbf{r}),x_{m}(\mathbf{r}))=\frac{\rho}{2}(x_{m}(\mathbf{r})-w_{m}(\mathbf{r}))^{2}+\frac{1}{2}(w_{m}(\mathbf{r})-m(\mathbf{r}))^{2}$

The cost function $\bar{L}_{b}(w_{b},x_{b})$ is separable across 3D pixel
index $\mathbf{r}$ because
$\mathcal{B_{C}}(\mathbf{w}_{b})=\sum_{\mathbf{r}}\bar{\mathcal{B_{C}}}(\mathbf{w}_{b}(\mathbf{r}))$
where

$\bar{\mathcal{B_{C}}}(\mathbf{w}_{b}(\mathbf{r}))=\begin{cases}\text{0}&\quad\text{if
}\mathbf{w}_{b}(\mathbf{r})\geq 0\\\ \infty&\quad\text{otherwise}\end{cases}$
(35)

The cost function reformulated as a sum over pixels can be stated as:

$\bar{L}_{b}(w_{b},x_{b})=\sum_{\mathbf{r}}\underbrace{\frac{\rho}{2}(\mathbf{x}_{b}(\mathbf{r})-\mathbf{w}_{b}(\mathbf{r}))^{2}+\bar{\mathcal{B_{C}}}(\mathbf{w}_{b}(\mathbf{r}))}_{{L}_{b}(\mathbf{w}_{b}(\mathbf{r}),\mathbf{x}_{b}(\mathbf{r}))}$
(36)

Now $\bar{L}_{s_{i}}({\mathbf{w}}_{s_{i}},{\mathbf{x}}_{s_{i}},\alpha_{s})$
can be expanded across pixels as follows:

$\bar{L}_{s_{i}}({\mathbf{w}}_{s_{i}},{\mathbf{x}}_{s_{i}},\alpha_{s})=\sum_{\mathbf{k}}\underbrace{\frac{\rho}{2}\|{\mathbf{x}}_{s_{i}}(\mathbf{k})-{\mathbf{w}}_{s_{i}}(\mathbf{k})\|_{2}^{2}+\sqrt{2}\alpha_{s}\|A_{s}{\mathbf{w}}_{s_{i}}(\mathbf{k})\|_{{2}}}_{{L}_{s}({\mathbf{w}}_{s_{i}}(\mathbf{k}),{\mathbf{x}}_{s_{i}}(\mathbf{k}),\alpha_{s})},$
(37)

Finally, $\bar{L}_{t}({\mathbf{w}}_{t},{\mathbf{x}}_{t},\alpha_{t})$ can also
be expanded across pixels courtesy the use of mixed vector matrix norm

$\bar{L}_{t}({\mathbf{w}}_{t},{\mathbf{x}}_{t},\alpha_{t})=\sum_{\mathbf{r}}\underbrace{\frac{\rho}{2}\|{\mathbf{x}}_{t}(\mathbf{r})-{\mathbf{w}}_{t}(\mathbf{r})\|_{2}^{2}+\alpha_{t}\|{\mathbf{w}}_{t}(\mathbf{r})\|_{2}}_{{L}_{t}({\mathbf{w}}_{t}(\mathbf{r}),{\mathbf{x}}_{t}(\mathbf{r}),\alpha_{t})}$
(38)

We have shown so far that all the sub-problems are separable across pixels.
Hence the solution to the minimization problems of equations (28), (29), (31),
and (30), can be expressed as following:

$\displaystyle\hat{w}_{m}(\mathbf{r})=$
$\displaystyle\underset{z\in\mathbb{R}}{\operatorname{argmin}}\;\;{L}_{m}(z,x_{m}(\mathbf{r})),$
(39) $\displaystyle\hat{\mathbf{w}}_{b}(\mathbf{r})=$
$\displaystyle\underset{z\in\mathbb{R}}{\operatorname{argmin}}\;\;{L}_{b}(z,\mathbf{x}_{b}(\mathbf{r})),$
(40) $\displaystyle\hat{\mathbf{w}}_{s_{i}}(\mathbf{k})=$
$\displaystyle\underset{{\mathbf{z}}\in\mathbb{R}^{10}}{\operatorname{argmin}}\;\;{L}_{s}({\mathbf{z}},{\mathbf{x}}_{s_{i}}(\mathbf{k}),\alpha_{s})\text{
, and}$ (41) $\displaystyle\hat{\mathbf{w}}_{t}(\mathbf{r})=$
$\displaystyle\underset{{\mathbf{z}}\in\mathbb{R}^{9}}{\operatorname{argmin}}\;\;{L}_{t}({\mathbf{z}},{\mathbf{x}}_{t}(\mathbf{r}),\alpha_{t}).$
(42)

####  6.1.2 Solution to the pixel-wise sub-problems

The solution to the ${w}_{m}$ sub-problem is obtained by exploiting the fact
that the cost ${L}_{m}(z,x_{m}(\mathbf{r}))$ is a differentiable function. The
minima is obtained by finding the stationary point of the cost function and
the resultant optimal point $\hat{w}_{m}(\mathbf{r})$ is:

$\hat{w}_{m}(\mathbf{r})=\frac{\rho x_{m}(\mathbf{r})+m(\mathbf{r})}{\rho+1}$
(43)

The solution to the $w_{b}$-problem is also simple, and it is the clipping of
the pixels by bound that defines the set $\mathcal{C}$ Parikh et al. (2014).
The optimal point $\hat{w}_{b}(\mathbf{r})$ as given below:

$\hat{\mathbf{w}}_{b}(\mathbf{r})=\mathbb{P}_{\mathcal{C}}(\mathbf{x}_{b}(\mathbf{r})),$
(44)

where ${\mathbb{P}}_{\mathcal{C}}(\cdot)$ denotes the operation of clipping
the pixel values within the bounds in definition of $\mathcal{C}$.

The $\mathbf{w}_{t}$ sub-problem could be understood as evaluating the well
known proximal operator Parikh et al. (2014) of $\ell_{2}$ norm
$\frac{\alpha_{t}}{\rho}\|\cdot\|_{2}$ at the point
${\mathbf{x}}_{t}(\mathbf{r})$. Hence the solution to the $\mathbf{w}_{t}$
sub-problem is given by

$\displaystyle\hat{\mathbf{w}}_{t}(\mathbf{r})$
$\displaystyle=\underset{{\mathbf{z}}\in\mathbb{R}^{9}}{\operatorname{argmin}}\;\;\frac{\rho}{2}\|{\mathbf{x}}_{t}(\mathbf{r})-\mathbf{z}\|_{2}^{2}+\alpha_{t}\|{\mathbf{z}}\|_{2}$
(45)
$\displaystyle=\begin{cases}\big{(}1-\frac{\alpha_{t}}{\rho\|{\mathbf{x}}_{t}(\mathbf{r})\|_{2}}\big{)}{\mathbf{x}}_{t}(\mathbf{r})&\|{\mathbf{x}}_{t}(\mathbf{r})\|_{2}\geq\frac{\alpha_{t}}{\rho}\\\
\mathbf{0}&otherwise\end{cases}$ (46)

The $\mathbf{w}_{s}$ sub-problem is more complicated as its involves
composition of a linear operator with a norm function. The solution of this
sub-problem is given by the following lemma:

###### Lemma 1

Let $\mathbf{y}\in\mathbb{R}^{10}$. Let $\mathbf{y}_{1}\in\mathbb{R}^{5}$ be
its sub-vector with first five entries and $\mathbf{y}_{2}\in\mathbb{R}^{5}$
be its sub-vector with last five entries. The solution of optimization problem
$\;\;\operatornamewithlimits{argmin}\limits_{\mathbf{z}\in\mathbb{R}^{10}}{\;\;\frac{\rho}{2}\|{{\mathbf{y}}}-\mathbf{z}\|_{2}^{2}+\sqrt{2}\alpha_{s}\|A_{s}\mathbf{z}\|_{{2}}}\;\;$
where $A_{s}$ is defined in (16) is given by

$\mathbf{z}^{*}=P\left[\begin{array}[]{c}\gamma{{\mathbf{y}}}_{1}\\\
{{\mathbf{y}}}_{2}\end{array}\right]$ (47)

Here $P$ defined in (57) is the eigenvector matrix of $A_{s}$ and
$\gamma=\max(0,1-\frac{\sqrt{2}\alpha_{s}}{\rho\|{{\mathbf{y}}}_{1}\|_{2}})$.

The proof of Lemma 1 is given in Appendix A. By applying this lemma, we can
conclude that the solution of $\mathbf{w}_{s}$ sub-problem is as follows:

$\hat{\mathbf{w}}_{s_{i}}(\mathbf{k})=P\left[\begin{array}[]{c}\gamma\mathbf{x}_{s_{i1}}(\mathbf{k})\\\
\mathbf{x}_{s_{i2}}(\mathbf{k})\end{array}\right]$ (48)

where $\mathbf{x}_{s_{i1}}(\mathbf{k})\in\mathbb{R}^{5}$ is the sub-vector of
$\mathbf{x}_{s_{i}}(\mathbf{k})$ with first five entries and
$\mathbf{x}_{s_{i2}}(\mathbf{k})\in\mathbb{R}^{5}$ is the sub-vector of
$\mathbf{x}_{s_{i}}(\mathbf{k})$ with last five entries.

###  6.2 The f sub-problem

$\mathbf{f}^{(k+1)}=\operatornamewithlimits{argmin}\limits_{\mathbf{f}}{L(\mathbf{f},\mathbf{w}^{(k+1)},{\boldsymbol{\beta}}^{(k)},\alpha_{s},\alpha_{t})}$
(49)

The sub-problem in variable given in (24) has a simpler form once you ignore
all the terms not depending on $\mathbf{f}$ in the optimization problem. The
simpler from of $\mathbf{f}$ sub-problem may be stated as follows:

$\displaystyle\mathbf{f}^{(k+1)}=$
$\displaystyle\operatornamewithlimits{argmin}\limits_{\mathbf{f}}\frac{1}{2}\|\mathbf{T}*\mathbf{f}-\mathbf{y}^{(k)}\|_{2}^{2}$
(50) where
$\displaystyle\mathbf{y}^{(k)}=\mathbf{w}^{(k+1)}-\frac{1}{\rho}\boldsymbol{\beta}^{(k)}$

For notational convenience, we let $\mathbf{y}=\mathbf{y}^{(k)}$ and
$\hat{\mathbf{f}}=\mathbf{f}^{(k+1)}$. Recall that
$\mathbf{f(\mathbf{r})}=[g(\mathbf{r})\;\;v(\mathbf{r})]^{\top}$From the
definition of $\mathbf{T}(\mathbf{r})$ it can be observed that the cost is
separable across the components $g$ and $v$ of $\mathbf{f}(\mathbf{r})$.
Assume that
$\mathbf{y}(\mathbf{r})=[{y}_{m}(\mathbf{r})\;{y}_{s,1}(\mathbf{r})\;{y}_{s,2}(\mathbf{r})\;{y}_{s,3}(\mathbf{r})\;{y}_{s,4}(\mathbf{r})\;{y}_{s,5}(\mathbf{r})\;{y}_{s,6}(\mathbf{r})\;{y}_{s,7}(\mathbf{r})\;{y}_{s,8}(\mathbf{r})\;{y}_{s,9}(\mathbf{r})\;{y}_{s,10}(\mathbf{r})\;{y}_{t,1}(\mathbf{r})\;\\\
{y}_{t,2}(\mathbf{r})\;{y}_{t,3}(\mathbf{r})\;{y}_{t,4}(\mathbf{r})\;{y}_{t,5}(\mathbf{r})\;{y}_{t,6}(\mathbf{r})\;{y}_{t,7}(\mathbf{r})\;{y}_{t,8}(\mathbf{r})\;{y}_{t,9}(\mathbf{r})\;{y}_{b,1}(\mathbf{r})\;{y}_{b,2}(\mathbf{r})]^{\top}$.
This simplification is achieved by observing the structure $\mathbf{y}^{(k)}$
it inherits from $\mathbf{w}^{(k+1)}$ and $\boldsymbol{\beta}^{(k)}$. Now the
cost separated along $g$ and $v$ is given by:

$\displaystyle\bar{L}_{1}(g)=$
$\displaystyle\frac{1}{2}\big{(}\|h*g-y_{m}\|_{2}^{2}+\|d_{xx}*g-y_{s,1}\|_{2}^{2}+\|d_{yy}*g-y_{s,2}\|_{2}^{2}+$
$\displaystyle\|d_{xy}*g-y_{s,3}\|_{2}^{2}+\|d_{xy}*g-y_{s,4}\|_{2}^{2}+\|g-y_{s,5}\|_{2}^{2}+\|g-y_{b,1}\|_{2}^{2}\big{)}$
$\displaystyle\bar{L}_{2}(v)=$
$\displaystyle\frac{1}{2}\big{(}\|d_{xx}*v-y_{s,6}\|_{2}^{2}+\|d_{yy}*v-y_{s,7}\|_{2}^{2}+\|d_{xy}*v-y_{s,8}\|_{2}^{2}+\|d_{xy}*v-y_{s,9}\|_{2}^{2}+$
$\displaystyle\|v-y_{s,10}\|_{2}^{2}+\|d_{xx}*v-y_{t,1}\|_{2}^{2}+\|d_{yy}*v-y_{t,2}\|_{2}^{2}+\|d_{xy}*v-y_{t,3}\|_{2}^{2}+$
$\displaystyle\|d_{xy}*v-y_{t,4}\|_{2}^{2}\|d_{yt}*v-y_{t,5}\|_{2}^{2}+\|d_{yt}*v-y_{t,6}\|_{2}^{2}+\|d_{xt}*v-y_{t,7}\|_{2}^{2}+$
$\displaystyle\|d_{xt}*v-y_{t,8}\|_{2}^{2}+\|d_{tt}*v-y_{t,9}\|_{2}^{2}+\|v-y_{b,2}\|_{2}^{2}\big{)}$

For notational convenience, we let $\mathbf{y}=\mathbf{y}^{(k)}$ and
$\hat{\mathbf{f}}=\mathbf{f}^{(k+1)}$. Recall that
$\mathbf{f(\mathbf{r})}=[g(\mathbf{r})\;\;v(\mathbf{r})]^{\top}$. The function
$\bar{L}_{1}(g)$ and $\bar{L}_{2}(v)$ are quadratic in nature in the variables
$g$ and $v$ respectively. The minima of both these functions can be obtained
by solving the equations$\nabla_{g}\bar{L}_{1}(g)=\mathbf{0}$ and
$\nabla_{v}\bar{L}_{2}(v)=\mathbf{0}$ respectively. This requires evaluation
of the gradient expressions which are given below.

$\displaystyle\nabla_{g}\bar{L}_{1}(g)=$
$\displaystyle\tilde{h}*h*g+\tilde{d}_{xx}*d_{xx}*g+\tilde{d}_{yy}*d_{yy}*g+2\tilde{d}_{xy}*d_{xy}*g+g$
(51)
$\displaystyle-\tilde{h}*y_{m}-\tilde{d}_{xx}*y_{s,1}-\tilde{d}_{yy}*y_{s,2}-\tilde{d}_{xy}*y_{s,3}-\tilde{d}_{xy}*y_{s,4}-y_{s,5}-y_{b,1}$
$\displaystyle\nabla_{v}\bar{L}_{2}(v)=$
$\displaystyle\tilde{d}_{xx}*d_{xx}*v+\tilde{d}_{yy}*d_{yy}*v+2\tilde{d}_{xy}*d_{xy}*v+2v+\tilde{d}_{xt}*d_{xt}*v$
(52)
$\displaystyle+\tilde{d}_{yt}*d_{yt}*v+\tilde{d}_{tt}*d_{tt}*v-\tilde{d}_{xx}*y_{s,5}-\tilde{d}_{yy}*y_{s,6}-\tilde{d}_{xy}*y_{s,7}-\tilde{h}*y_{m}-y_{s,8}$
$\displaystyle-\tilde{d}_{xt}*y_{t,1}-\tilde{d}_{yt}*y_{t,2}-\tilde{d}_{tt}*y_{t,3}-\tilde{d}_{xt}*y_{t,4}-\tilde{d}_{yt}*y_{t,5}-\tilde{d}_{tt}*y_{t,6}$

This completes the solution of the $\mathbf{f}$ sub-problem which form part of
the ADMM iterates.

##  7 Experiments

To demonstrate the effectiveness of STAIC regularization, we consider
restoration of time varying TIRF images, and compare with the method ICTV,
3D-TV2 and CST. We selected five image sequences obtained from a high
Numerical Aperture (high bandwidth) TIRF microscope under nearly noise-free
conditions, and designate them as the ground truth models. The ground truth
models are given in Figure 1. We then simulate measured images by blurring
these models with PSF corresponding to low NA systems and by adding mixed
Poisson-Gaussian noise as shown below :

$m=\mathcal{P}(\gamma_{p}(h*g))+\eta$

Here $\gamma_{p}$ is a parameter to control the strength of Poisson noise. We
consider PSF corresponding to five NA values namely 0.8, 0.9, 1.0, 1.1, and
1.2. We kept the Gaussian noise level fixed, and considered two $\gamma_{p}$
values of 1 and 5 to generate the dataset. A higher value of $\gamma_{p}$
results in a lower level of Poisson noise corruption. This makes a set of 50
measured image sequences generated from our five ground truth images.

###  7.1 Experiment 1

In first experiment, we choose to demonstrate the advantages gained by
treating restoration of time varying images as a separate problem as opposed
to treating them under classical image restoration by treating time as the
third dimension. To this end we consider 3D-TV2 regularization based
restoration of our dataset of 50 images. The resultant restoration scheme
which we call 3D-TV2 is posed as the following optimization problem

$\hat{g}=\operatornamewithlimits{argmin}\limits_{g}{\;\;\frac{1}{2}\sum_{i=1}^{n_{F}}\|(h*g_{i})-m_{i}\|_{F}^{2}+\lambda\sum_{\mathbf{r}}\|(M*g)(\mathbf{r})\|_{F}}$
(53)

where $M(\mathbf{r})$ is the 3D Hessian operator defined earlier and $\lambda$
is the regularization parameter. The algorithm is tuned for $\lambda$ to
obtain the best SNR. The results are presented in Tables 1 and 2 for the two
Poisson noise levels $\gamma_{p}=1,5$ respectively along with results of the
STAIC scheme. The resultant SNR (in dB) was not found to be competitive in
comparison to our STAIC scheme or even other algorithms we used in subsequent
experiments. The SSIM scores also demonstrate the shortcomings of this
approach. It also demonstrates that we indeed need to treat the time dimension
differently by considering temporal variations as a different phenomenon in
relation to spatial variations in the TIRF image.

###  7.2 Experiment 2

In the second set of experiments we compare the performance of proposed STAIC
scheme against the Combined Spatio-Temporal Regularization (CST) and Infimal
Convolution TV (ICTV-2DT). In CST the measured image is subject to both
spatial and temporal regularization by employing a sum of norms
regularization. The corresponding image restoration optimization problems
using CST regularization takes the form:

$\operatornamewithlimits{argmin}\limits_{g}{\;\;\frac{1}{2}\sum_{i=1}^{n_{F}}\|(h*g_{i})-m_{i}\|_{F}^{2}+\lambda
R_{CST}(g)}$ (54)

In addition, we also generated the restoration results of applying ICTV-2DT
regularization scheme. As discussed earlier ICTV belongs to the family of
infimal convolution regularization where spatial and temporal components of
signal are modeled differently by the regularization. The ICTV-2DT
optimization problem takes the form:

$\operatornamewithlimits{argmin}\limits_{g}{\;\;\frac{1}{2}\sum_{i=1}^{n_{F}}\|(h*g_{i})-m_{i}\|_{F}^{2}+\lambda
R_{ICTV}(g)}$ (55)

The results of ICTV-2DT, CST and STAIC are aslo presented in Tables 1 and 2
for the two Poisson noise levels. Our proposed algorithm STAIC perform better
than both CST and ICTV-2DT regularized approaches to restoration for the TIRF
2DT dataset. STAIC was found to result in a higher restoration quality as
measured using SNR as well as SSIM which is another popular image quality
measure.. It may be noted that SSIM of restored spatio-temporal image
presented in both the tables is the average of SSIM values across the image
frames. Our STAIC achieves SNR improvement over the second best approach in
the range of $0.5$ dB to $2$ dB in the 50 image dataset considered. This
improvement can be attributed to the improved regularization design that
captures the prior of temporal component of the signal in comparison to the
other algorithms. In SSIM measure too we have a clear gain over the other
approaches considered in our experiments.

Figure 1:  The five original images  | STAIC | ICTV | CST | 3DTV2
---|---|---|---|---
Image | NA | ssim | snr | ssim | snr | ssim | snr | ssim | snr
5135 | 0.8 | 0.853 | 11.58 | 0.758 | 8.84 | 0.841 | 8.80 | 0.656 | 8.48
0.9 | 0.873 | 13.01 | 0.781 | 9.91 | 0.863 | 9.86 | 0.669 | 9.45
1 | 0.887 | 14.31 | 0.799 | 10.84 | 0.880 | 10.80 | 0.677 | 10.28
1.1 | 0.897 | 15.46 | 0.812 | 11.65 | 0.891 | 11.61 | 0.683 | 10.98
1.2 | 0.905 | 16.54 | 0.822 | 12.39 | 0.901 | 12.34 | 0.688 | 11.60
5142 | 0.8 | 0.912 | 11.88 | 0.814 | 10.17 | 0.886 | 10.18 | 0.620 | 9.96
0.9 | 0.926 | 13.26 | 0.831 | 11.30 | 0.905 | 11.31 | 0.629 | 11.01
1 | 0.934 | 14.52 | 0.848 | 12.31 | 0.918 | 12.32 | 0.636 | 11.94
1.1 | 0.942 | 15.68 | 0.865 | 13.21 | 0.928 | 13.22 | 0.640 | 12.75
1.2 | 0.946 | 16.78 | 0.876 | 14.04 | 0.936 | 14.05 | 0.643 | 13.49
5147 | 0.8 | 0.908 | 12.48 | 0.901 | 11.94 | 0.903 | 11.95 | 0.746 | 11.70
0.9 | 0.926 | 13.67 | 0.914 | 13.05 | 0.921 | 13.06 | 0.755 | 12.73
1 | 0.939 | 14.73 | 0.914 | 14.04 | 0.934 | 14.05 | 0.761 | 13.64
1.1 | 0.947 | 15.70 | 0.923 | 14.93 | 0.943 | 14.94 | 0.766 | 14.44
1.2 | 0.954 | 16.58 | 0.930 | 15.74 | 0.950 | 15.75 | 0.769 | 15.15
5157 | 0.8 | 0.874 | 8.04 | 0.794 | 7.28 | 0.814 | 7.29 | 0.615 | 6.94
0.9 | 0.901 | 9.46 | 0.835 | 8.59 | 0.850 | 8.60 | 0.628 | 8.11
1 | 0.919 | 10.75 | 0.854 | 9.76 | 0.877 | 9.76 | 0.638 | 9.11
1.1 | 0.932 | 11.90 | 0.844 | 10.81 | 0.897 | 10.79 | 0.644 | 9.96
1.2 | 0.941 | 12.96 | 0.853 | 11.75 | 0.912 | 11.73 | 0.649 | 10.70
5158 | 0.8 | 0.943 | 14.38 | 0.918 | 13.18 | 0.927 | 13.14 | 0.513 | 12.00
0.9 | 0.954 | 15.76 | 0.917 | 14.37 | 0.940 | 14.33 | 0.519 | 12.86
1 | 0.962 | 17.02 | 0.928 | 15.44 | 0.949 | 15.39 | 0.522 | 13.58
1.1 | 0.967 | 18.15 | 0.936 | 16.38 | 0.955 | 16.33 | 0.524 | 14.16
1.2 | 0.970 | 19.21 | 0.941 | 17.24 | 0.960 | 17.18 | 0.526 | 14.64
Table 1: The Restoration Quality for Noise Level 1 ($\gamma_{p}=1$) | STAIC | ICTV | CST | 3DTV2
---|---|---|---|---
Image | NA | ssim | snr | ssim | snr | ssim | snr | ssim | snr
5135 | 0.8 | 0.867 | 11.60 | 0.761 | 8.85 | 0.857 | 8.81 | 0.658 | 8.48
0.9 | 0.889 | 13.03 | 0.789 | 9.91 | 0.881 | 9.87 | 0.671 | 9.45
1 | 0.903 | 14.34 | 0.809 | 10.85 | 0.898 | 10.81 | 0.678 | 10.28
1.1 | 0.914 | 15.52 | 0.823 | 11.68 | 0.910 | 11.63 | 0.684 | 10.99
1.2 | 0.922 | 16.59 | 0.835 | 12.41 | 0.920 | 12.36 | 0.689 | 11.60
5142 | 0.8 | 0.929 | 11.89 | 0.892 | 10.19 | 0.892 | 10.18 | 0.621 | 9.96
0.9 | 0.943 | 13.28 | 0.908 | 11.32 | 0.911 | 11.31 | 0.629 | 11.01
1 | 0.953 | 14.55 | 0.920 | 12.33 | 0.925 | 12.32 | 0.636 | 11.94
1.1 | 0.960 | 15.72 | 0.929 | 13.24 | 0.936 | 13.22 | 0.640 | 12.75
1.2 | 0.966 | 16.84 | 0.935 | 14.08 | 0.944 | 14.06 | 0.643 | 13.49
5147 | 0.8 | 0.914 | 12.49 | 0.910 | 11.95 | 0.907 | 11.95 | 0.746 | 11.70
0.9 | 0.932 | 13.68 | 0.924 | 13.07 | 0.925 | 13.06 | 0.756 | 12.74
1 | 0.945 | 14.75 | 0.935 | 14.06 | 0.938 | 14.06 | 0.762 | 13.64
1.1 | 0.954 | 15.71 | 0.942 | 14.95 | 0.948 | 14.94 | 0.766 | 14.44
1.2 | 0.961 | 16.60 | 0.948 | 15.77 | 0.955 | 15.76 | 0.770 | 15.15
5157 | 0.8 | 0.881 | 8.04 | 0.795 | 7.28 | 0.815 | 7.29 | 0.615 | 6.94
0.9 | 0.908 | 9.47 | 0.836 | 8.59 | 0.852 | 8.60 | 0.628 | 8.11
1 | 0.927 | 10.76 | 0.855 | 9.77 | 0.879 | 9.77 | 0.638 | 9.12
1.1 | 0.940 | 11.91 | 0.869 | 10.81 | 0.899 | 10.80 | 0.644 | 9.96
1.2 | 0.949 | 12.97 | 0.854 | 11.76 | 0.914 | 11.73 | 0.649 | 10.70
5158 | 0.8 | 0.949 | 14.39 | 0.927 | 13.20 | 0.931 | 13.14 | 0.514 | 12.00
0.9 | 0.961 | 15.78 | 0.939 | 14.40 | 0.945 | 14.34 | 0.519 | 12.86
1 | 0.969 | 17.04 | 0.947 | 15.48 | 0.954 | 15.40 | 0.522 | 13.57
1.1 | 0.974 | 18.19 | 0.953 | 16.44 | 0.960 | 16.35 | 0.525 | 14.16
1.2 | 0.978 | 19.26 | 0.958 | 17.30 | 0.965 | 17.20 | 0.526 | 14.64
Table 2: The Restoration Quality for Noise Level 2 $(\gamma_{p}=5)$

##  8 Conclusion

We proposed a spatio-temporal regularization method for restoration of time
varying 2D images. Our design of the regularization terms led to an underlying
additive decomposition of the signal enabling differentiated regularization of
regions with motion and without motion. The effectiveness of this
decomposition is evident in the superior reconstruction quality evaluated by
structure similarity as well as signal to noise ratio measures of image
quality. Our method was found to perform better than other algorithms each
designed using different philosophy towards spatio-temporal signal
restoration. This better quality reconstruction was enabled through design of
a fast and efficient algorithm for minimizing the resultant regularized cost
using the ADMM framework. The proposed algorithms effectiveness can be
attributed to the advantages provided by the infimal convolution framework in
exploiting correlation across frames inherent to spatio-temporal signals.

##  Appendix A proof

###### Proof 1

Let us start by considering the optimization problem

$\underset{{\mathbf{z}}\in\mathbb{R}^{10}}{\operatorname{argmin}}\;\;{L}_{s}(\mathbf{z},\mathbf{y},\alpha_{s})$
(56)

where
${L}_{s}(\mathbf{z},\mathbf{y},\alpha_{s})=\frac{\rho}{2}\|\mathbf{y}-\mathbf{z}\|_{2}^{2}+\sqrt{2}\alpha_{s}\|A_{s}\mathbf{z}\|_{{2}}$.
Note that $\mathbf{y,z}\in\mathbb{R}^{10}$.

Consider symmetric matrix $A_{s}^{\top}A_{s}$ which is orthogonally
diagonalizable. In other words, we can find an orthogonal matrix of
eigenvectors $P$ and a diagonal matrix of eigenvalues $D$ such that
$A_{s}^{\top}A_{s}=PDP^{\top}$.

Since $A_{s}^{\top}A_{s}$ can be derived from given matrix $A_{s}$, we can
derive $P$ and $D$ by eigen decomposition and is given below:

$P=\left[\begin{array}[]{cccccccccc}\frac{1}{\sqrt{2}}&0&0&0&0&-\frac{1}{\sqrt{2}}&0&0&0&0\\\
0&\frac{1}{\sqrt{2}}&0&0&0&0&-\frac{1}{\sqrt{2}}&0&0&0\\\
0&0&\frac{1}{\sqrt{2}}&0&0&0&0&-\frac{1}{\sqrt{2}}&0&0\\\
0&0&0&\frac{1}{\sqrt{2}}&0&0&0&0&-\frac{1}{\sqrt{2}}&0\\\
0&0&0&0&\frac{1}{\sqrt{2}}&0&0&0&0&-\frac{1}{\sqrt{2}}\\\
-\frac{1}{\sqrt{2}}&0&0&0&0&\frac{1}{\sqrt{2}}&0&0&0&0\\\
0&-\frac{1}{\sqrt{2}}&0&0&0&0&\frac{1}{\sqrt{2}}&0&0&0\\\
0&0&-\frac{1}{\sqrt{2}}&0&0&0&0&\frac{1}{\sqrt{2}}&0&0\\\
0&0&0&-\frac{1}{\sqrt{2}}&0&0&0&0&\frac{1}{\sqrt{2}}&0\\\
0&0&0&0&-\frac{1}{\sqrt{2}}&0&0&0&0&\frac{1}{\sqrt{2}}\\\
\end{array}\right]\text{ and }D=diag(\left[\begin{array}[]{cccccccc}1\\\ 1\\\
1\\\ 1\\\ 1\\\ 0\\\ 0\\\ 0\\\ 0\\\ 0\end{array}\right])$ (57)

where $\operatorname{diag}$ is the operation that creates a diagonal matrix
with its input vector forming the diagonal entries. We define
$\hat{\mathbf{z}}=P^{\top}\mathbf{z}$ and
$\hat{\mathbf{y}}=P^{\top}\mathbf{y}$. Since $P$ is orthogonal, it is also
true that $\|P^{\top}\mathbf{x}\|_{2}=\|\mathbf{x}\|_{2}$. Since
$\|\mathbf{x}\|_{2}=\sqrt{\langle\mathbf{x},\mathbf{x}\rangle}$ where
$\langle\cdot,\cdot\rangle$ is the standard inner product, it can be seen that

$\|A_{s}\mathbf{z}\|_{2}=\sqrt{\langle
A_{s}\mathbf{z},A_{s}\mathbf{z}\rangle}=\sqrt{\langle\mathbf{z},A_{s}^{\top}A_{s}\mathbf{z}\rangle}=\sqrt{\mathbf{z}^{\top}A_{s}^{\top}A_{s}\mathbf{z}}=\sqrt{\mathbf{z}^{\top}PDP^{\top}\mathbf{z}}=\sqrt{\hat{\mathbf{z}}^{\top}D\hat{\mathbf{z}}}$
(58)

Also observe that:

$\frac{\rho}{2}\|\mathbf{y}-\mathbf{z}\|_{2}^{2}=\frac{\rho}{2}\|P^{\top}\mathbf{y}-P^{\top}\mathbf{z}\|_{2}^{2}=\frac{\rho}{2}\|\hat{\mathbf{y}}-\hat{\mathbf{z}}\|_{2}^{2}$
(59)

Let us define $\hat{\mathbf{z}}_{1},\hat{\mathbf{z}}_{2}\in\mathbb{R}^{5}$
such that $\hat{\mathbf{z}}=\left[\begin{array}[]{c}\hat{\mathbf{z}}_{1}\\\
\hat{\mathbf{z}}_{2}\end{array}\right]$. Also define
$\hat{\mathbf{y}_{1}},\hat{\mathbf{y}_{2}}\in\mathbb{R}^{4}$ such that
$\hat{\mathbf{y}}=\left[\begin{array}[]{c}\hat{\mathbf{y}_{1}}\\\
\hat{\mathbf{y}_{2}}\end{array}\right]$

The original optimization problem in terms of $\hat{\mathbf{z}}$ is given by

$\underset{{\mathbf{z}}\in\mathbb{R}^{10}}{\operatorname{argmin}}\;\;\frac{\rho}{2}\|\hat{\mathbf{y}}-\hat{\mathbf{z}}\|_{2}^{2}+\sqrt{2}\alpha_{s}\sqrt{\hat{\mathbf{z}}D\hat{\mathbf{z}}}$
(60)

By appealing to the sub vectors definition in $\hat{\mathbf{z}}$ and
$\hat{\mathbf{y}}$, the problem can be stated as

$\underset{\hat{\mathbf{z}}_{1},\hat{\mathbf{z}}_{2}\in\mathbb{R}^{5}}{\operatorname{argmin}}\;\;\frac{\rho}{2}\|\hat{\mathbf{y}_{1}}-\hat{\mathbf{z}}_{1}\|_{2}^{2}+\frac{\rho}{2}\|\hat{\mathbf{y}_{2}}-\hat{\mathbf{z}}_{2}\|_{2}^{2}+\sqrt{2}\alpha_{s}\|{\hat{\mathbf{z}_{1}}}\|_{2}$
(61)

By considering two separate optimization problems in $\hat{\mathbf{z}}_{1}$
and $\hat{\mathbf{z}}_{2}$, we get

$\displaystyle\underset{\hat{\mathbf{z}}_{1}\in\mathbb{R}^{5}}{\operatorname{argmin}}\;\;\frac{\rho}{2}\|\hat{\mathbf{y}_{1}}-\hat{\mathbf{z}}_{1}\|_{2}^{2}+\sqrt{2}\alpha_{s}\|{\hat{\mathbf{z}_{1}}}\|_{2}$
(62)
$\displaystyle\underset{\hat{\mathbf{z}}_{2}\in\mathbb{R}^{5}}{\operatorname{argmin}}\;\;\frac{\rho}{2}\|\hat{\mathbf{y}_{2}}-\hat{\mathbf{z}}_{2}\|_{2}^{2}$
(63)

The optimization problem in variable $\hat{\mathbf{z}}_{1}$ is simply the
proximal operator of $\ell_{2}$ norm which is well known Parikh et al. (2014)
and is given by

$\hat{z}_{1}^{*}=\gamma\hat{\mathbf{y}_{1}}$ (64)

where
$\gamma=\max(0,1-\frac{\sqrt{2}\alpha_{s}}{\rho\|\hat{\mathbf{y}_{1}}\|_{2}})$.
The solution of optimization problem in variable $\hat{\mathbf{z}}_{2}$ is
given by $\hat{\mathbf{z}}_{2}^{*}=\hat{\mathbf{y}_{2}}$. Now, we map back to
space of variable $\mathbf{z}$ by multiplying with $P$

${\mathbf{z}}^{*}=P\hat{\mathbf{z}}^{*}$ (65)

where $\hat{z}^{*}=\left[\begin{array}[]{c}\gamma\hat{\mathbf{y}_{1}}\\\
\hat{\mathbf{y}_{2}}\end{array}\right]$.

## References

  * Katsaggelos [1989] Aggelos K Katsaggelos.  Iterative image restoration algorithms.  _Optical engineering_ , 28(7):735–748, 1989.
  * Engl et al. [1996] Heinz Werner Engl, Martin Hanke, and Andreas Neubauer.  _Regularization of inverse problems_ , volume 375.  Springer Science & Business Media, 1996.
  * Rangayyan [2004] Rangaraj M Rangayyan.  _Biomedical image analysis_.  CRC press, 2004.
  * Ramani et al. [2012] Sathish Ramani, Zhihao Liu, Jeffrey Rosen, Jon-Fredrik Nielsen, and Jeffrey A Fessler.  Regularization parameter selection for nonlinear iterative image restoration and mri reconstruction using gcv and sure-based methods.  _IEEE Transactions on Image Processing_ , 21(8):3659–3672, 2012.
  * Viswanath et al. [2020] Sanjay Viswanath, Manu Ghulyani, Simon De Beco, Maxime Dahan, and Muthuvel Arigovindan.  Image restoration by combined order regularization with optimal spatial adaptation.  _IEEE Transactions on Image Processing_ , 29:6315–6329, 2020.
  * Arigovindan [2013] Muthuvel Arigovindan.  Image deconvolution research: its scope and importance in live cell microscopy.  _Current Science_ , pages 1501–1511, 2013.
  * Li et al. [2018] Jizhou Li, Feng Xue, and Thierry Blu.  Accurate 3d psf estimation from a wide-field microscopy image.  In _2018 IEEE 15th International Symposium on Biomedical Imaging (ISBI 2018)_ , pages 501–504. IEEE, 2018.
  * Fan et al. [2019] Junchao Fan, Xiaoshuai Huang, Liuju Li, Liangyi Chen, and Shan Tan.  One-step deconvolution for multi-angle tirf microscopy with enhanced resolution.  _Biomedical optics express_ , 10(3):1097–1110, 2019.
  * Lefkimmiatis et al. [2011] Stamatios Lefkimmiatis, Aurélien Bourquard, and Michael Unser.  Hessian-based norm regularization for image restoration with biomedical applications.  _IEEE Transactions on Image Processing_ , 21(3):983–995, 2011.
  * Lefkimmiatis et al. [2013] Stamatios Lefkimmiatis, John Paul Ward, and Michael Unser.  Hessian schatten-norm regularization for linear inverse problems.  _IEEE transactions on image processing_ , 22(5):1873–1888, 2013.
  * Lindsten et al. [2011] Fredrik Lindsten, Henrik Ohlsson, and Lennart Ljung.  Clustering using sum-of-norms regularization: With application to particle filter output computation.  In _2011 IEEE Statistical Signal Processing Workshop (SSP)_ , pages 201–204. IEEE, 2011.
  * Holler and Kunisch [2014] Martin Holler and Karl Kunisch.  On infimal convolution of tv-type functionals and applications to video and image reconstruction.  _SIAM Journal on Imaging Sciences_ , 7(4):2258–2300, 2014.
  * Guo et al. [2014] Weihong Guo, Jing Qin, and Wotao Yin.  A new detail-preserving regularization scheme.  _SIAM journal on imaging sciences_ , 7(2):1309–1334, 2014.
  * Bredies et al. [2010] Kristian Bredies, Karl Kunisch, and Thomas Pock.  Total generalized variation.  _SIAM Journal on Imaging Sciences_ , 3(3):492–526, 2010.
  * Ulyanov et al. [2018] Dmitry Ulyanov, Andrea Vedaldi, and Victor Lempitsky.  Deep image prior.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pages 9446–9454, 2018.
  * Xing et al. [2017] Fuyong Xing, Yuanpu Xie, Hai Su, Fujun Liu, and Lin Yang.  Deep learning in microscopy image analysis: A survey.  _IEEE transactions on neural networks and learning systems_ , 29(10):4550–4568, 2017.
  * von Chamier et al. [2021] Lucas von Chamier, Romain F Laine, Johanna Jukkala, Christoph Spahn, Daniel Krentzel, Elias Nehme, Martina Lerche, Sara Hernández-Pérez, Pieta K Mattila, Eleni Karinou, et al.  Democratising deep learning for microscopy with zerocostdl4mic.  _Nature communications_ , 12(1):1–18, 2021.
  * Liu et al. [2021] Zhichao Liu, Luhong Jin, Jincheng Chen, Qiuyu Fang, Sergey Ablameyko, Zhaozheng Yin, and Yingke Xu.  A survey on applications of deep learning in microscopy image analysis.  _Computers in Biology and Medicine_ , 134:104523, 2021.
  * Hoffman et al. [2021] David P Hoffman, Isaac Slavitt, and Casey A Fitzpatrick.  The promise and peril of deep learning in microscopy.  _Nature methods_ , 18(2):131–132, 2021.
  * Schloegl et al. [2017] Matthias Schloegl, Martin Holler, Andreas Schwarzl, Kristian Bredies, and Rudolf Stollberger.  Infimal convolution of total generalized variation functionals for dynamic mri.  _Magnetic resonance in medicine_ , 78(1):142–155, 2017.
  * Zhang et al. [2018] Jiahan Zhang, Si Li, Andrzej Krol, C Ross Schmidtlein, Edward Lipson, David Feiglin, and Yuesheng Xu.  Infimal convolution-based regularization for spect reconstruction.  _Medical physics_ , 45(12):5397–5410, 2018.
  * Trémoulhéac et al. [2014] Benjamin Trémoulhéac, Nikolaos Dikaios, David Atkinson, and Simon R Arridge.  Dynamic mr image reconstruction–separation from undersampled (k,t)-space via low-rank plus sparse prior.  _IEEE transactions on medical imaging_ , 33(8):1689–1701, 2014.
  * Chen et al. [2022] Eric Z Chen, Puyang Wang, Xiao Chen, Terrence Chen, and Shanhui Sun.  Pyramid convolutional rnn for mri image reconstruction.  _IEEE Transactions on Medical Imaging_ , 2022.
  * Arigovindan et al. [2013] Muthuvel Arigovindan, Jennifer C Fung, Daniel Elnatan, Vito Mennella, Yee-Hung Mark Chan, Michael Pollard, Eric Branlund, John W Sedat, and David A Agard.  High-resolution restoration of 3d structures from widefield images with extreme low signal-to-noise-ratio.  _Proceedings of the National Academy of Sciences_ , 110(43):17344–17349, 2013.
  * Boyd et al. [2011] Stephen Boyd, Neal Parikh, and Eric Chu.  _Distributed optimization and statistical learning via the alternating direction method of multipliers_.  Now Publishers Inc, 2011.
  * Parikh et al. [2014] Neal Parikh, Stephen Boyd, et al.  Proximal algorithms.  _Foundations and trends in Optimization_ , 1(3):127–239, 2014.