\section{Introduction}

In recent papers it has been argued that grand unification based on groups larger than $SU(5)$ has several features that make it a natural context for theories of dark matter \cite{unified-DM, BC-sphaleron}.
(1) The fermion multiplets of such unified groups typically contain Standard Model singlets, which could play the role of dark matter. (2) Unified theories can have accidental global symmetries (a well-known example is $B-L$ in minimal $SU(5)$), and one of these could be the quantum number $X$ that stabilizes dark matter. And (3), unification groups of rank 6 or larger can contain additional non-abelian interactions, which can play a role in generating dark matter. One way this can happen is through the sphaleron processes \cite{sphaleron} of the new non-abelian interaction, which can convert baryons and leptons into dark matter particles, thus giving a scenario in which dark matter and baryonic matter are ``co-generated."

Many papers have explored the idea that sphalerons are involved in the genesis of dark matter. Usually, these scenarios involve only the sphalerons of the electroweak interactions \cite{BCF, ADM-sphaleron}. However, as emphasized in \cite{ADM-nonSMsph} (and later in \cite{BC-sphaleron}), the idea that the sphaleron processes of a new non-Standard Model interaction do the co-generation has several attractive features. First, if electroweak sphalerons produced the dark matter asymmetry, the dark matter particles have to be chiral under electroweak $SU(2)$, which makes them more visible both directly and through their loop effects, severely constraining such models. By contrast, to be produced by the sphalerons of a new interaction, the dark matter particles need only be chiral under the new interaction and can be vectorlike under the electroweak $SU(2)$, which allows them to be heavy and to have negligible loop effects on Standard Model predictions. Such models are therefore less severely constrained and much easier to build. Second, having two kinds of sphalerons --- those of the electroweak $SU(2)$ and those of a new gauge interaction --- means that there are {\it two} sphaleron equilibrium conditions, which is enough to determine the ratios of $B$, $L$, and $X$. That means (under the assumption that the dark matter is almost purely asymmetric, i.e. $n_{\overline{DM}} \ll n_{DM}$) that the number density of dark matter particles today can be calculated, giving a definite (though of course model-dependent) prediction for the mass of the dark matter particle. As sphaleron processes typically make $B$, $L$, and $X$ comparable, and because $\Omega_B$ and $\Omega_{DM}$ are comparable, one typically finds that the mass of the dark mass particle is GeV-scale .

To have the dark matter be asymmetric \cite{nussinov, BCF, ADM-sphaleron, ADM-nonSMsph, ADM-other, Graesser, LYZ}, there must have been some processes in the early universe that annihilated the ``symmetric component" of the dark matter. That is, most of the anti-dark matter particles must have annihilated leaving just the asymmetry (as was the case with baryons). This is a non-trivial issue if the dark matter particles have masses of order a GeV, as then
annihilation via electroweak interactions is not efficient enough and would have left a relic abundance of anti-dark matter particles that was far too large. One way to solve this problem \cite{BC-sphaleron} is to posit the existence of a light field (typically it must have mass less than the weak scale) that mediates an interaction {\it by which} dark matter particles and their anti-particles can annihilate, as well as some new massless (or extremely light) particles {\it into which} they can annihilate. Thus the existence of an entire sector of light particles, besides the dark matter particles themselves, would be postulated {\it ad hoc}.

In this paper, we point out another way: The same new non-abelian gauge interactions whose sphalerons co-generate dark matter in the first place, are also able to annihilate the symmetric component of the dark matter. This leads to a very economical, tightly integrated scenario, in which everything comes from unification in larger groups: the dark matter fermions themselves, and the new interaction that generates them and annihilates the symmetric component.

At first, it might seem that the new non-abelian gauge interactions whose sphaleron processes generate the dark matter asymmetry cannot do the annihilation efficiently enough, any more than the weak interactions could, for the scale characterizing their breaking must be even larger than the weak scale. However, as we show in this paper, they can do it through a mechanism that we call``pre-annihilation." This mechanism was proposed in the context of other scenarios by Kitano and Low \cite{KitanoLow} and by Kang and Li \cite{KangLi}. (It can also be seen as a limiting case of the ``exodus" mechanism of Unwin \cite{Unwin}.) Because the scenario we are describing is tightly integrated and to a large extent determined by the grand unified structure, one can use cosmological observables to constrain it.

The idea of pre-annihilation is that an asymmetry in $X$ is first produced among some heavy metastable dark particles, which are able to annihilate efficiently enough to make negligible the number of their anti-particles. These heavy dark particles then decay late, out of equilibrium, to stable dark particles that ``inherit" their asymmetry and form the dark matter that is observed today. (We shall call the heavy dark particles $\chi$ and the light stable dark particles $\zeta$.) We shall show that the lifetime of these late decays must fall within a certain window. A nice feature of the scenario we are discussing is that these decays can result from effective dim-5 operators that would give the right lifetime if they were suppressed by $O(1/M_{P\ell})$ or $O(1/M_{GUT})$.

The basic cosmological scenario we shall discuss has the following stages. When the temperature is above the mass scale
of the heavy dark particles, $m_{\chi}$ (which is much larger than the weak scale but less than $3 \times 10^6$ GeV), the
sphaleron processes of the new non-Standard Model gauge interaction, whose gauge group we shall call $G_*$, produce an
$X$ asymmetry comparable to the baryon asymmetry. This asymmetry is in the form of an excess of $\chi$ over their
anti-particles $\chi^c$. As the temperature drops below $m_{\chi}$, the $\chi$ and $\chi^c$ annihilate by means of the
$G_*$ gauge interaction, leaving only the excess $\chi$ due to the $X$ asymmetry. At a later epoch (when $T$ is between
$10^{-4} m_{\chi}$ and $2 \times 10^{-7} m_{\chi}$), each heavy dark particle $\chi$ decays via interactions suppressed by $1/M_{GUT}$ or $1/M_{P \ell}$ into one light dark matter particle $\zeta$ plus some Standard Model quarks and leptons that equilibrate. The $\zeta$ are singlets under all gauge interactions --- except possibly ones broken at $M_{GUT}$ --- and are never in thermal equilibrium. (At some very early time, their density was driven to exponentially small values by inflation, and the only such particles present today are those caused by the decay of the heavy dark particles.) Thus the stable light dark matter particles inherit the asymmetry produced at very large scales. Because they have never been in equilibrium the dark matter particles are produced with a mean momentum significantly higher than the photon temperature. Depending on the parameters of the model, this scenario can give either cold or warm dark matter.

As noted above, some of the ingredients of the scenario we are proposing and studying here
have been considered elsewhere, such as the use of non-Standard Model sphalerons to co-generate dark matter \cite{ADM-nonSMsph, BC-sphaleron}, and the pre-annihilation of the symmetric component of dark matter \cite{KitanoLow,KangLi,Unwin}. What is new here is that we combine these elements into a relatively simple scenario and show that it emerges naturally in the context of grand unification based on groups of rank 6 or greater. And within models of the type we are proposing, we derive significant new limits coming from several considerations. These include (a) the condition that the decays of the heavy unstable dark matter not be so slow as to make the resulting stable dark matter hot, and (b) that the decays of the heavy unstable dark matter not be so rapid as to regenerate a large symmetric component of dark matter. These limits give a relatively narrow window of decay lifetimes of the heavy dark matter, which we show turns out to be nicely consistent with decay via GUT- or Planck-suppressed operators that one might expect to exist in realistic models.

All of this we illustrate in a specific grand unified model based on $SU(7)$. In that model, we show how the various critical ingredients tend naturally to exist in grand unified models with large gauge groups. The model illustrates, in particular, how the dark matter and ordinary matter particles are unified in GUT multiplets; how the $G_*$ interactions are unified with the Standard Model interactions; how in a typical case there are quarks, leptons and dark matter particles that transform under $G_*$, as needed for the co-generation mechanism; and how the quantum number $X$ can arise as an accidental symmetry in a grand unified model.

We shall first discuss the cosmology of the scenario in more detail in section 2, and then the details of the illustrative $SU(7)$ model in section 3.

\section{The details of pre-annihilation}

The new gauge interaction based on the group $G_*$ is broken at a scale $M_*$ by the VEV of a Higgs field that we shall call $\Omega$. It is the sphalerons of the $G_*$ interaction that will violate $X$ and relate the $X$ asymmetry to the $B$ and $L$ asymmetries. It is also the $G_*$ gauge interaction that will mediate the annihilation of the heavy dark particles ($\chi$) with their anti-particles ($\chi^c$). The heavy dark particles ($\chi$) which have $X\neq 0$, are neutral under the Standard Model gauge interactions, but transform in a chiral way under $G_*$. (It is because of this chirality, of course, that the $G_*$ sphalerons violate $X$.) The particles $\chi$ and $\chi^c$ obtain masses of order $M_*$ from the VEV of $\langle \Omega \rangle$. We can therefore denote the scale of the breaking of $G_*$ by $m_{\chi}$ (instead of $M_*$) and shall do so from now on. Some of the Standard Model quarks and leptons will also transform under $G_*$, so that $G_*$ sphalerons also violate $B$ and $L$. (A simple model having all these features will be presented in section 3.)

The light dark particles ($\zeta$) and anti-particles ($\zeta^c$) are singlets under $SU(3)_c \times SU(2)_L \times U(1)_Y \times G_*$. They couple to the other particles of the model through dimension-5 operators suppressed by $1/M_{GUT}$ or $1/M_{P \ell}$.

We assume that at very early times inflation wiped out $B$, $L$, and $X$. After
(or during) reheating an asymmetry of one or more of these quantum numbers was generated.
It will not matter for our discussion which asymmetry was generated or what process generated it. We assume that no significant number of $\zeta$ and $\zeta^c$ particles appeared during reheating, because they have no coupling to the inflaton and only couple to other light particles through GUT-scale-suppressed operators. (This will put an upper limit on the reheat temperature that depends on the GUT scale, as will be discussed later.)

When $T$ was large compared to $m_{\chi}$, the $G_*$ sphaleron processes were happening and produced an equilibrium among $B$, $X$, and $L$, making these
asymmetries comparable.  (Of course, as is common, we use the term``sphaleron process" loosely to mean any anomalous $G_*$ interaction, not just the ones involving tunneling via the sphaleron configuration.) Let us define $\epsilon_B \equiv (n_B - n_{\overline{B}})/s$, $\epsilon_L \equiv (n_L - n_{\overline{L}})/s$, and $\epsilon_X \equiv (n_X - n_{\overline{X}})/s$. In order to have $\Omega_{DM} \approx 5 \Omega_B$, one requires that $\epsilon_X \approx (5 m_p/m_{\zeta}) \epsilon_B$, where $m_p$ is the proton mass.

When $T$ dropped below $m_{\chi}$, the number of
$\chi$ and $\chi^c$ became Boltzmann suppressed as they annihilated into lighter states.
These annihilations would have happened primarily through $G_*$ gauge interactions (and partly through $\Omega$ exchange). (Therefore,
they would not have annihilated into the light dark particles $\zeta + \zeta^c$.)   The relic
abundance of an asymmetric particle can be expressed in terms of the asymmetry parameter $r$,
which is defined to be
\begin{equation}
r = \frac{n_\chi^c}{n_\chi}.
\end{equation}
For an annihilation cross section scaling as $\sigma_0 (m_{\chi}/T)^{-n}$,
the value of $r$ at late times, $r_\infty$, is given by \cite{Graesser,LYZ} (see also earlier calculations in Refs. \cite{Earlierfreezeout}):
\begin{equation}
\label{rinf1}
r_\infty \approx \exp\left(-\frac{\epsilon_\chi \lambda \sqrt{g_{*f}}}{x_f^{n+1}(n+1)} \right),
\end{equation}
where $\lambda = 0.264~ M_{P\ell} m_{\chi} \sigma_0$ (with $M_{P\ell}$ denoting the Planck mass), $g_{*f}$ is the effective
number of relativistic degrees of freedom when the annihilations freeze out at temperature $T_f$, and $x_f = m_{\chi}/T_f$.
(We will use $g_*$ throughout to denote relativistic degrees of freedom, with the subscript labeled
to denote the temperature at which this quantity is evaluated in each case).  Taking the annihilation
in our model to be $s$-wave ($n=0$) with $\sigma_0 = \alpha_*/m_{\chi}^2$, Eq. (\ref{rinf1}) simplifies to
\begin{equation}
\label{rinf2}
r_\infty \approx \exp\left(- \frac{0.264~\epsilon_\chi \alpha_* M_{P\ell} \sqrt{g_{*f}}}{m_{\chi} x_f} \right)
\end{equation}

Our goal is to achieve $\epsilon_\chi \sim 10^{-10}$ (the same order of magnitude as for baryons).
Somewhat arbitrarily, we will take the dividing line between symmetric ($r_\infty = 1$) and
asymmetric ($r_\infty \ll 1$) freeze-out to be the point at which the argument in the exponential of Eq. (\ref{rinf2}) is $-1$;
asymmetric freeze-out will occur when the absolute value of this argument is $> 1$.
For the model parameters considered here,
$\sqrt{g_{*f}} \sim 20$ and $x_f \sim 20$. (As noted in Ref. \cite{Graesser}, $x_f$ is little changed from its symmetric
value unless $r_\infty \ll 1$).  Assuming that $\alpha_* \sim 10^{-2}$, we find that the $\chi$ freeze-out will be asymmetric
as long as $m_{\chi} < 3 \times 10^6$ GeV = 3000 TeV.  When this limit is satisfied, to a good approximation the only heavy dark particles
$\chi$ that remain after freeze-out are those that existed due to the $X$ asymmetry.

As mentioned, we assume each $\chi$ decays
via a dimension-5 GUT suppressed operator into a light dark particle $\zeta$ plus
Standard Model particles. Let us call the coefficient of this dimension-5 operator
$c_5/M_{GUT}$, where $M_{GUT}$ is the scale where the Standard Model gauge couplings unify, and $c_5$ is a dimensionless number. The $\chi$ decays will happen when the temperature of the universe is $T_{dec}$, given by

\begin{equation}
\left( \frac{c_5}{M_{GUT}} \right)^2 m_{\chi}^3 \sim \sqrt{g_{*dec}} \frac{T_{dec}^2}{M_{P\ell}}
\;\; \Longrightarrow \;\; \frac{T_{dec}}{m_{\chi}} \sim g_{*dec}^{-1/4} \frac{\sqrt{m_{\chi} M_{P \ell}}}{(M_{GUT}/c_5)}.
\end{equation}

So far, we have been giving a fairly standard analysis of asymmetric freeze out. We now show that there are two considerations that yield
both an upper bound and a lower bound on $T_{dec}$ that must be satisfied to achieve a satisfactory asymmetric dark matter scenario. The upper bound comes from considering the $\chi^c$ that decay into $\zeta^c$ when $T \sim m_{\chi}$, when the
$\chi^c$ are not Boltzmann suppressed. Of course, only a small fraction do, since
these decays have a lifetime much larger than the age of the universe when $T \sim m_{\chi}$.
But this fraction must be smaller than $\epsilon_X \cong \frac{ 5 m_p}{m_{\zeta}} \epsilon_B \sim 10^{-10}$ if one is not to produce
a number of $\zeta^c$ that is comparable to the $\zeta$ asymmetry. This gives a lower bound on the lifetime of the decay, and thus an upper bound on $T_{dec}$.
One has

\begin{equation}
\left( \frac{n_{\zeta^c}}{s} \right)_{T_{dec}} \sim
\frac{t_{m_{\chi}}}{t_{dec}} \left( \frac{n_{\chi^c}}{s} \right)_{T \sim m_{\chi}}
\sim \frac{t_{m_{\chi}}}{t_{dec}} \frac{1}{g_{*{m_{\chi}}}},
\end{equation}

\noindent where $t_{m_{\chi}}$ is the age of the universe when $T \sim m_{\chi}$,
and $t_{dec}$ is the age of the universe when $T = T_{dec}$. Demanding that this be less than
$\epsilon_X \cong \frac{ 5 m_p}{m_{\zeta}} \epsilon_B$, gives

\begin{equation}
\frac{t_{m_{\chi}}}{t_{dec}} \frac{1}{g_{*{m_{\chi}}}} =
\frac{(\sqrt{g_{*{m_{\chi}}}} m_{\chi}^2/M_{P \ell})^{-1}}{(\sqrt{g_{*dec}} T_{dec}^2/M_{P \ell})^{-1}} \frac{1}{g_{*{m_{\chi}}}} < \frac{ 5 m_p}{m_{\zeta}} \epsilon_B
\;\; \Longrightarrow  \frac{T_{dec}}{m_{\chi}} < \left( \frac{5 m_p}{m_{\zeta}} \right)^{1/2}
\frac{g_{*{m_{\chi}}}^{3/4}}{g_{*dec}^{1/4}} \epsilon_B^{1/2} \sim 10^{-4} .
\end{equation}
We have neglected two additional effects that will tend to suppress the decays of the $\chi^c$.
At $T \sim m_{\chi}$, the decay rate will be Lorentz suppressed, and inverse decays
will be energetically allowed.
Both of these will decrease the number of decay-produced $\zeta^c$, so our limit is a conservative one.

The lower bound on $T_{dec}$ comes from requiring that the dark matter streaming length be small enough to allow for the formation of the observed
large-scale structure. When the dark matter particles $\zeta$ are produced by decays of the
$\chi$, they have a momentum that is of order $m_{\chi}/3$ (assuming three-body decays).  This momentum redshifts as the inverse of the scale factor, so
a lower $T_{dec}$ allows for less redshift prior to the epoch of structure formation, resulting in a larger velocity and
corresponding larger
free-streaming
length.

The comoving free-streaming length $\lambda_{FS}$ up to the time $t_{eq}$ when the matter and radiation
densities become equal and density perturbations can begin to grow, is given by \cite{KolbTurner}
\begin{equation}
\label{lambdaFS}
\lambda_{FS} = \int_0^{t_{eq}} v(t^\prime) \frac{R_0}{R(t^\prime)} dt^\prime,
\end{equation}
where $R$ is the scale factor, and the zero subscript refers to present-day quantities.
Using the approximation that the velocity is $v=1$ when $\zeta$ is relativistic and scales as
$1/R$ when it becomes nonrelativistic at $t_{NR}$, we obtain \cite{KolbTurner}
\begin{equation}
\lambda_{FS} \approx t_{NR}\left(\frac{R_0}{R_{NR}}\right)\left[2+ \ln\left(\frac{t_{eq}}{t_{NR}}\right)\right].
\end{equation}

The momentum of $\zeta$ at scale factor $R$ is given by
$p = (m_{\chi}/3)(R_{dec}/R)$, and $\zeta$ becomes nonrelativistic when $p \approx m_\zeta$, so we have
\begin{equation}
T_{NR} = T_{dec} \left(\frac{3 m_\zeta}{m_{\chi}}\right)\left(\frac{g_{*dec}}{g_{*NR}}\right)^{1/3},
\end{equation}
where we have used the fact that
$g_* (RT)^3$ is constant \cite{KolbTurner}.  In the radiation-dominated era, $t = 0.30 g_*^{-1/2} M_{P\ell} T^{-2}$\cite{KolbTurner},
so we obtain
\begin{equation}
t_{NR} = 0.03 M_{P\ell}\left(\frac{m_{\chi}}{m_\zeta}\right)^2 T_{dec}^{-2}~ g_{*NR}^{1/6} ~g_{*dec}^{-2/3}.
\end{equation}
Further, $R_0/R_{NR}$ in Eq. (\ref{lambdaFS}) can be written in terms of $T_{dec}$:
\begin{equation}
\frac{R_0}{R_{NR}}
= \frac{3 m_\zeta}{m_{\chi}}\left(\frac{T_{dec}}{T_0}\right)\left(\frac{g_{*dec}}{43/4}\right)^{1/3}\left(
\frac{11}{4}\right)^{1/3},
\end{equation}
where the two factors on the right-hand side include the heating from the annihilation of all particles heavier than electrons, and the
heating due to electron-positron annihilation, respectively.

Assembling all of these factors and using the measured values for $t_{eq}$ and $T_0$, the free-streaming length can be written as
\begin{equation}
\lambda_{FS} = 0.2 ~{\rm pc}~ \left(\frac{5 m_p}{m_\zeta}\right)\left(\frac{m_{\chi}}{T_{dec}}\right)~ g_{*NR}^{1/6}~ g_{*dec}^{-1/3}
\left(1 - 0.04~\ln\left[\left(\frac{5 m_p}{m_\zeta}\right)\left(\frac{m_{\chi}}{T_{dec}}\right)~g_{*NR}^{1/12}
~g_{*dec}^{-1/3} \right]\right).
\end{equation}

A reasonable order-of-magnitude upper bound on $\lambda_{FS}$ to allow structure formation is $\lambda_{FS} \la$ 100 kpc.
Using this limit and appropriate values for the degrees of freedom, we obtain:
\begin{equation}
\label{FSlimit}
\frac{T_{dec}}{m_{\chi}} > 2 \times 10^{-7}\left(\frac{5 m_p}{m_\zeta}\right).
\end{equation}
Roughly speaking, values of ${T_{dec}}/{m_{\chi}}$ no more than an order of magnitude above this bound will
produce $\zeta$ that behaves as
warm dark matter, while larger values of ${T_{dec}}/{m_{\chi}}$ will yield
cold dark matter.

Clearly, there is a broad range between the limits given in Eqs. (6) and (\ref{FSlimit}). Substituting Eq. (4) into those inequalities gives a range for $M_{GUT}/c_5$

\begin{equation}
2 \times 10^{16} \; {\rm GeV} \left( \frac{m_{\zeta}}{5 m_p} \right)^{1/2} \left( \frac{m_{\chi}}{3 \times
10^6 ~{\rm GeV}} \right)^{1/2} < \frac{M_{GUT}}{c_5}
< 8 \times 10^{18} \; {\rm GeV} \left( \frac{m_{\zeta}}{5 m_p} \right)
\left( \frac{m_{\chi}}{3 \times 10^6~ {\rm GeV}} \right)^{1/2}.
\end{equation}

We have assumed that a negligible number of $\zeta^c$ are generated by scattering of thermal particles after reheating. In particular, we require that
$n_{\zeta^c}/s < \epsilon_X = (5 m_p/m_{\zeta}) \epsilon_B$.  This gives the following upper limit on the reheat temperature:

\begin{equation}
\left( \frac{n_{\zeta^c}}{s} \right)_{T_{rh}}
\sim \frac{c_5^{\;2} T_{rh}^3}{M_{GUT}^2} \; \frac{M_{P\ell}}{\sqrt{g_{*rh}} T_{rh}^2}
\; \frac{1}{g_{*rh}} < (5 m_p/m_{\zeta}) \epsilon_B
\; \Longrightarrow \; T_{rh} < \frac{(M_{GUT}/c_5)^2}{3 \times 10^{25} {\rm GeV}}
\left( \frac{5 m_p}{m_{\zeta}} \right).
\end{equation}

All of the bounds we have derived in this section are based on the assumption that there is a strong asymmetry in the dark matter, i.e. $n_{\overline{DM}} \ll n_{DM}$. We regard this as an interesting limit. Moreover, in models where the dark matter asymmetry is generated by sphalerons of a new non-abelian interaction, assuming that limit typically allows a prediction  for the mass of the dark matter particle, as noted in \cite{ADM-nonSMsph}.
But it is interesting also to consider the more general case, where there is a weaker asymmetry in dark matter. This would allow smaller masses for the light stable dark matter particle $\zeta$. It would also allow $M_*$ (and $m_{\chi}$) to be larger, as less efficient pre-annihilation is allowed. It would also relax the upper bound on
$T_{dec}/m_{\chi}$. We leave a detailed analysis of the more general case to future work
\cite{BCS}.

\section{An $SU(7)$ example}

We now present a simple $SU(7)$ grand unified model that illustrates the ideas discussed above. As noted in the Introduction, grand unification based on large groups has several features that make it a natural context for models of dark matter. (1) The large unified fermion multiplets typically contain SM singlets that could be dark matter candidates, so that such models would {\it unify} dark matter with quarks and leptons. (2) The symmetry that stabilizes the dark matter can arise as an accidental global symmetry in a manner analogous to the way $U(1)_{B-L}$ arises as an accidental symmetry in minimal $SU(5)$. And (3), the large unified groups (if they have rank  $\geq 6$) contain extra non-abelian subgroups that can play the role of $G_*$ in co-generating and pre-annihilating dark matter. That is, the $G_*$ interactions can be unified with the Standard Model interactions. All of these features are realized in the $SU(7)$ model we present in this section.

The group $SU(7)$ contains the subgroups $SU(7) \supset SU(5) \times SU(2)_* \times U(1)_{T_7} \supset G_{SM} \times SU(2)_* \times U(1)_{T_7}$, where $G_{SM}$ is the Standard Model gauge group, and where the generator $T_7$ of the extra abelian factor
can be written (in the fundamental representation of $SU(7)$) as $T_7 =
diag (\frac{2}{7}, \frac{2}{7}, \frac{2}{7}, \frac{2}{7}, \frac{2}{7}, -\frac{5}{7}, - \frac{5}{7})$. We assume that $SU(7)$ breaks to $G_{SM} \times SU(2)_* \times U(1)_{T_7}$ at superlarge scales.

The fermions of one family consist of ${\bf 21} + 3 \times \overline{{\bf 7}}$, which we denote
$\psi^{[AB]}$, $\psi_A$, $\psi^{(k)}_A, k = 1,2$. The indices $A, B, ... $ etc. are $SU(7)$ indices, while $k$ is just a label distinguishing two of the anti-fundamental fermion multiplets. This is the smallest possible anomaly-free set that gives one SM family. In addition, each family is assumed to have three $SU(7)$-singlet fermions that will be denoted simply $\psi^{(m)}$, where $m = 1,2,3$ is a label distinguishing them.
Besides the Higgs fields that break $SU(7)$ at superlarge scales, there are Higgs multiplets $h^{[ABC]}$, $h'_A$, and $\Omega^{(k)}_A$, $k = 1,2$. These will break the electroweak $SU(2)_L$ and the $SU(2)_*$ and give mass to the fermions, as will be seen. (Here too the index $k$ is just a label distinguishing the two $\Omega_A$ multiplets.)

We shall use lower-case Greek letters for indices of the $SU(5)$ subgroup of $SU(7)$, so $\alpha, \beta, \gamma, ... = 1, ..., 5$. We shall use capital Latin letters from the middle of the alphabet for indices of the $SU(2)_*$ subgroup of $SU(7)$, so $I, J, ... = 6,7$. We suppress all family indices. The fermion multiplets decompose under $SU(5) \times SU(2)_*$ as shown in Table I. In that table, the $SU(7)$ fermion multiplets of one family are listed across the top row. In the column under each $SU(7)$ multiplet are given the $SU(5) \times SU(2)_*$ multiplets contained in it. The superscripts $0, \pm 1, \pm 2$ on the $SU(5) \times SU(2)_*$ multiplets are the values of the global charge $X$, which will be discussed shortly. The three rows denote three types of fermion: ``SM" stands for Standard Model fermions, which have $X= 0$; ``VL" stands for vectorlike fermions, which have $X = \pm 1$; and ``DM" stands for dark matter fermions, which have $X = \pm 2$.

\noindent {\large\bf Table I:} The decomposition under $SU(5) \times SU(2)_*$ of the fermion multiplets of a family. The superscripts $0$, $\pm 1$, $\pm 2$ are the values of $X$.

\vspace{0.2cm}
\[
\begin{tabular}{|l|l|l|l|l|}
\hline & $\psi^{[AB]} = {\bf 21}$ & $\psi_A = \overline{{\bf 7}}$ & $\psi^{(k)}_A = 2 \times \overline{{\bf 7}}$ & $\psi^{(m)} = 3 \times {\bf 1}$ \\
\hline
SM & $\psi^{\alpha \beta}$ & $\psi_{\alpha}$ & &  \\
& $= (10, 1)^0$ & $=(\overline{5}, 1)^0$ & & \\
\hline
VL & $\psi^{\alpha I}$ & $\psi_I$ & $\psi^{(k)}_{\alpha}$ &  \\
& $= (5,2)^{-1}$ & $=(1,2)^{+1}$ & $= 2\times (\overline{5}, 1)^{+1}$ & \\
\hline
DM & $\psi^{[IJ]}$ & & $\psi^{(k)}_I$ & $\psi^{(m)}$ \\
& $= (1,1)^{-2}$ & & $=2 \times (1,2)^{+2}$ & $=3 \times (1,1)^{-2}$ \\
\hline
\end{tabular}
\]

\noindent Note that, as is required for co-generation of dark matter and ordinary matter by $SU(2)_*$ sphalerons, there are dark matter particles that transform under $SU(2)_*$
(namely those in $(1,2)^{+2}$)
as well as quarks and leptons that do so (namely those in $(5,2)^{-1}$). This is not special to $SU(7)$, but would be a generic feature of models in which $SU(2)_*$ is
unified with $G_{SM}$ in a simple group.

The Higgs multiplets are assumed to have the vacuum expectation values (VEVs)
$\langle h^{267} \rangle = v_u \sim M_W$, $\langle h'_2 \rangle = v_d \sim M_W$,
and $\langle \Omega^{(k)}_6 \rangle \sim \langle \Omega^{(k)}_7 \rangle \sim M_*$. (The indices of the $SU(2)_L$ electroweak subgroup of $SU(5)$ are $\alpha = 1,2$, while the color indices are $\alpha = 3,4,5$.)
The VEV $v_u$ will give mass to the up quarks; the VEV $v_d$ will give mass to the down quarks and charged leptons; and the VEVs of $\Omega^{(k)}_I $ will give mass to the vectorlike fermions and to the dark matter, as well as breaking $SU(2)_*$. To align these VEVs there must be a coupling connecting these Higgs multiplets. A coupling sufficient to do this is
$h^{[ABC]} h'_A \Omega^{(1)}_B \Omega^{(2)}_C$.

The Yukawa couplings are as follows

\begin{equation}
{\cal L}_{Yukawa} \supset (\psi^{[AB} \psi^{CD} ) h^{EFG]}, \;\; (\psi^{[AB]} \psi_A)
h'_B, \;\; (\psi^{[AB]} \psi^{(k)}_A ) \Omega^{(k')}_B, \;\;
(\psi^{(k)}_A \psi^{(m)}) (\Omega^{(k')})^{*A}.
\end{equation}

\noindent Assumed {\it not} to be present are couplings similar to those in
Eq. (16), except with $h'_A \leftrightarrow \Omega^{(k)}_A$. It is easy to check that there is an accidental $U(1) \times U(1)$ global symmetry consistent with all the Yukawa couplings in Eq. (16) and the Higgs self-coupling $h^{[ABC]} h'_A \Omega^{(1)}_B \Omega^{(2)}_C$. We shall be interested in a subgroup of this, which we will call $U(1)_R \times K$, where $K$ is a $Z_2$ symmetry. The $R$ charges and $K$ parities of the various fermion and Higgs multiplets defined so far are given by

\begin{equation}
R^K(\psi^{[AB]}, \psi_A, \psi^{(k)}_A, \psi^{(m)}, h^{[ABC]}, h'_A, \Omega^{(k)}_A)
= (-4^-, 2^+, 9^+, -14^-, 8^+, 2^-, -5^-).
\end{equation}

\noindent When the $\Omega^{(k)}_I = (1,2)$ get non-zero VEVs of order $M_*$, the group $SU(2)_*$ is spontaneously broken, and the group $U(1)_{T_7} \times U(1)_R$ is broken down to a global
symmetry that we shall call $U(1)_X$, whose generator is

\begin{equation}
X= T_7 + \frac{1}{7}R.
\end{equation}

\noindent One easily sees from Eqs. (17) and (18) and the definition $T_7 \cong diag (\frac{2}{7}, \frac{2}{7}, \frac{2}{7}, \frac{2}{7}, \frac{2}{7}, -\frac{5}{7}, - \frac{5}{7})$ that the $SU(5) \times SU(2)_*$ multiplets in Table I have the $X$ values shown there.
One also sees that the Higgs VEVs $\langle h^{267} \rangle, \langle h'_2 \rangle$
and $\langle \Omega^{(k)}_6 \rangle \sim \langle \Omega^{(k)}_7 \rangle$ do not break $U(1)_X$.

\noindent The masses of the Standard Model up quarks come from the coupling
$(\psi^{[\alpha \beta} \psi^{\gamma \delta}) \langle h^{267]} \rangle$. Those of the Standard Model down quarks and charged leptons come from $(\psi^{[\alpha 2]} \psi_\alpha) \langle h'_2 \rangle$, with $\alpha$ being a color index (3,4,5) for the down quarks and a weak index (1) for the charged leptons. (Of course, in a realistic model of quark and lepton masses there must always be other terms to avoid the ``bad" unification relation that the charged leptons and down quarks have equal mass at the GUT scale. But there are several ways of adding such terms, and this does not affect anything else we say here.)
The heavy vectorlike fermions  $\psi^{[AI]} = (5,2)^{-1}$ and $\psi^{(k)}_A = 2 \times
(5,1)^{+1}$ (see the row VL in Table I)) obtain $O(M_*)$ masses from $(\psi^{[AI]} \psi^{(k)}_A) \Omega^{(k')}_I \rangle$. The Yukawa couplings so far mentioned leave the vectorlike fermions $\psi_I = (1,2)^{+1}$ massless. (The coupling $(\psi^{[2I} \psi_I)
\langle h'_2 \rangle$ mixes them with the heavy neutrinos, but leaves the resulting mixture massless.) The masslessness of these fermions would not be a phenomenological problem, but they can be given Dirac masses by introducing some gauge singlet fermions
into the model.

The heavy dark matter particles get mass of $O(M_*)$ from
$(\psi^{[IJ]} \psi^{(k)}_I) \langle \Omega^{(k')}_J \rangle$ and
$(\psi^{(k)}_I \psi^{(m)}) \langle (\Omega^{(k')})^{*I} \rangle$. The $\psi^{(k)}_I$
(of which there are four per family, as $k = 1,2$ and $I = 6,7$) are what we called $\chi$ in earlier sections. The $\psi^{[IJ]}$ (of which there is one per family) and the
$\psi^{(m)}$ (of which there are three per family) are what we called $\chi^c$.

In order for the pre-annihilation mechanism discussed in previous sections to work, the $\chi$ and $\chi^c$ must be able to decay to light stable dark matter particles through dim-5 operators that are suppressed by $O(1/M_{GUT})$ or $O(1/M_{P\ell})$. One such operator is $(\psi^{IJ} \zeta) h^*_{[2IJ]} h'^{*2}/M_{GUT}$. Since $\psi^{IJ}$ is one of the $\chi^c$, as noted in the previous paragraph, and $h^{[2IJ]}$ and $h'_2$ are just the weak Higgs doublets $h_u$ and $h_d$ whose VEVs do the electroweak breaking, we may write this operator as $(\chi^c \zeta) h^* h'^*/M_{GUT}$. The $SU(7)$-singlet fermion field $\zeta$ can be assigned $R^K= 14^+$ (implying by Eq. (18) that it has $X = 2$ as the other dark matter particles $\psi^{(k)}_I$ do), in which case this dim-5 operator conserves $R$, $K$ and $X$.
In order to give a Dirac mass to $\zeta$, let us introduce an SU(7)-singlet fermion $\zeta^c$ that has $R^K = -14^+$ and therefore $X= -2$. (Notice that $K$ parity prevents
the $\zeta$ from having a mass term with $\psi^{(m)}$, and $\zeta^c$ from mixing with  $\psi^{(m)}$.)

The Dirac fermion made up of $\zeta$ and $\zeta^c$ must have mass of order 1 GeV. The smallness of this mass can be explained in a way that is technically natural \cite{naturalness} by assuming it arises from the dynamical breaking of a symmetry. Let us suppose, therefore that there a $Z_2$ symmetry, which we will call $P$, under which $\zeta^c$ is odd, while all the other fields defined up to this point are even. (That is, all the fermions in Table I, the fermion $\zeta$, and the Higgs fields $h^{[ABC]}, h'_A, \Omega^{(k)}_A$ are even under $P$.)
The $P$ prevents an explicit (or bare) mass term for $\zeta, \zeta^c$. But such a mass can arise dynamically as follows. Suppose that there are terms
$\zeta^c \zeta S$, $\overline{\Psi} \Psi S^*$, and $M_{GUT}^2 S^* S$, where $S$ is a singlet scalar and $\overline{\Psi}$ and $\Psi$ are fermions that transform under a force that confines at a scale $\Lambda$ and acquire a condensate $\langle \overline{\Psi} \Psi \rangle \sim \Lambda^3$. Assume that $S$ and $\Psi^c$ are odd under $P$ while $\Psi$ is even, so as to allow these couplings.
Then $m_{\zeta} \sim \langle S \rangle \sim \Lambda^3/M_{GUT}^2$, which can be of order 1 GeV if the confinement scale $\Lambda \sim 10^{11}$ GeV.

The mechanism just described for making the mass of $\zeta$ small may seem somewhat contrived, but it is ``natural" in the technical sense \cite{naturalness}, {\it i.e.} in the sense that no coefficient of a renormalizable term in the Lagrangian has to be ``fine-tuned" to be very small. Perhaps a somewhat more elegant way could be found to
make $m_{\zeta}$ small. But this example does serve to illustrate the general fact that the mass of a fermion can be much smaller than the breaking scale of the symmetry that ``protects" it (here the symmetry $P$), without violating the condition of ``technical naturalness". Another example of this (which in fact is closely analogous to the mechanism just described) is the type II see-saw mechanism for neutrino masses, where neutrino masses are protected by the electroweak $SU(2)_L$, but the neutrino masses that result from the breaking of $SU(2)_L$ at the scale
Fermi scale $v$ are only of order $v^2/M_{GUT} \ll v$.

With all the particles of the model having the $R$ charges and $K$ and $P$ parities given above, it is easy to see that the lowest dimension operator that connects the light dark matter particles $\zeta$ and $\zeta^c$ to the other particles of the low energy theory (i.e. those with mass less than or of order $M_*$) is the dim-5 operator $(\psi^{IJ} \zeta) h^*_{2IJ} h'^{*2}$, which is suppressed by a large ($O(M_{P\ell})$ or $O(M_{GUT})$) mass. That is important, as seen in the previous sections, as otherwise there would be a danger that interactions of $\zeta$ and $\zeta^c$ with other matter would produce a thermal abundance of them and and lead to a large symmetric component of dark matter.

One can see that the quantum number $X$ serves to stabilize not only the dark matter, but also the heavy vectorlike $5 + \overline{5}$ fermions. This is easiest to see in the case (which we shall henceforth assume) that all the vectorlike (VL) fermions, which have odd values of $X$, are heavier than all the DM fermions, which have even values of $X$. Then the lightest particle with odd $X$ will be absolutely stable, and the lightest particle with $X = \pm 2$ will also be stable. That $X$ conservation stabilizes both the dark matter and the vectorlike fermions seems to be a generic feature of the kind of model we are discussing if it is unified in a simple group. If $X$ remained exactly unbroken, there would be residual heavy quarks and leptons from these vectorlike multiplets that would come to dominate the energy density of the universe. It is thus necessary that some small spontaneous breaking of $X$ in the model. It must be large enough for the decay of the heavy vectorlike fermions to be sufficiently fast not to cause cosmological problems, but small enough that (a) it does not cause too rapid decay of the heavy dark matter $\chi, \chi^c$ (i.e. dominating over the decays caused by the dim-5 operator), and (b) $X$-violating processes do not wipe out any asymmetries in dark matter.

Fortunately, there are simple ways that such $X$ violation can be introduced into the model that satisfy these conditions. One way is to suppose that there are $X$-violating VEVs of the Higgs multiplets $h'_A$ and $\Omega^{(k)}_A$ that are of order $\epsilon \ll 1$ times the $X$-conserving VEVs:
$\langle h'_I \rangle \sim \epsilon v_d \sim \epsilon M_W$, and
$\langle \Omega^{(k)}_2 \rangle \sim \epsilon M_*$. Then $\Delta X = \pm 1$ effects are suppressed by $O(\epsilon)$. This would have the following effects. (1) There would be an $O(\epsilon)$ mixing of the SM down quarks $d = \psi ^{[a2]}$ with the heavy vectorlike down quarks $D = \psi^{[aI]}$, so that $D$ could decay into
$u + W_-$ with an amplitude suppressed by a factor of $\epsilon$. (2) There would be similar mixing of the SM charged leptons $\ell^-$ with the heavy vectorlike charged leptons $L^-$, allowing the latter also to decay weakly with an amplitude suppressed by a factor
of $\epsilon$.
(3) There would be mixing among all the neutral fermions, namely the SM neutrinos ($\nu$), the heavy vectorlike neutrinos ($N$), and the heavy dark matter particles ($\chi$). This would allow weak decay of the heavy dark matter particles into $ \ell^- + W^+$. However, this amplitude would be suppressed by order $\epsilon^2$. The reason is that the diagonalization of the mass matrices of the charged and neutral fermions turns out to give electroweak lepton doublets of the following form:

\begin{equation}
\left( \begin{array}{l} \nu + O(\epsilon) N + O(\epsilon^2) \chi \\
\ell^- + O(\epsilon) L^- \end{array} \right), \;\;\;
\left( \begin{array}{l} N + O(\epsilon) \nu + O(\epsilon) \chi \\
L^- + O(\epsilon) \ell^- \end{array} \right).
\end{equation}

\noindent The powers of $\epsilon$ are directly traceable to the number of units by which $X$ is violated. Because a decay of a heavy dark matter particle into Standard Model particles violates $X$ by two units, its amplitude is suppressed by $\epsilon^2$.

The condition that the $X$-violating decay of the heavy dark matter be slow compared to that caused by dim-5 operator is that $\epsilon^4 M_* \ll M_*^3/M_{GUT}^2 \Rightarrow
\epsilon^2 \ll M_*/M_{P \ell}$, which is also the condition that $X$-violating scattering be slow enough not to wipe out the dark matter asymmetries. But that condition still allows the weak decays of the heavy vectorlike quarks and leptons to happen very rapidly; indeed with a lifetime short compared to the Hubble time when $T \sim M_*$. (These issues will be
analyzed in more detail in a future paper \cite{BCS}.)

It should be noted that despite the spontaneous violation of $X$, the $\zeta, \zeta^c$ are extremely stable on the time scale of the present age of the universe, since the lowest-dimension operators that couple $\zeta, \zeta^c$ to the light particles of the model have dimension 5 and are therefore suppressed by $O(1/M_{GUT})$ or $O(1/M_{P \ell})$, and also because their decay amplitude to Standard Model particles would be further suppressed by two powers of $\epsilon$. For example, if their decay amplitudes are of order $\epsilon^2/M_{GUT}$, they have a lifetime of $\tau_{\zeta} \sim \left( m_{\zeta}^3 \frac{\epsilon^4}{M_{P\ell}^2}
\right)^{-1}$, which for $m_{\zeta} \sim 5$ GeV and $\epsilon^2 \ll M_*/M_{P \ell}$  gives $\tau_{\zeta} \gg 10^{29}$ yrs.

In the $SU(7)$ model we have presented, both the dark matter particles and heavy vectorlike quarks and leptons carry non-zero $X$. (This is a generic feature of such models when unified in a simple group.) As one can see from Table I, one can write $X = 2N_{DM} + N_{VL}$, where $N_{DM}$ is dark matter number (the number of dark matter particles minus the number of anti-dark matter particles) and $N_{VL}$ is the vectorlike particle number. The contribution of the dark matter fermions and vectorlike fermions to $\Delta X$ in $SU(2)_*$ sphaleron processes cancel, as can be seen from Table I. In other words, $SU(2)_*$ sphaleron processes violate $N_{DM}$ and $N_{VL}$, but do not violate $X$.
(This is also a generic feature of models of the type we have been discussing when unified in a simple group.) Thus, whatever primordial asymmetry is generated after reheating must have $X = 2N_{DM} + N_{VL} \neq 0$, as otherwise $SU(2)_*$ sphalerons will wipe out the asymmetries in both $N_{DM}$ and $N_{VL}$. (This is analogous to the fact that electroweak sphalerons do not violate $B-L$, so that unless a primordial asymmetry in $B-L$ exists electroweak sphalerons will wipe out both $B$ and $L$.)

\section{Conclusions}

We have presented a scenario in which a new non-abelian gauge interaction, whose group we called $G_*$, is responsible for both co-generating a dark matter asymmetry and pre-annihilating the symmetric component of the dark matter. The co-generation of dark matter with ordinary matter is done by the sphaleron processes of $G_*$, while scattering mediated by the exchange
of $G_*$ gauge bosons accomplishes the pre-annihilation of dark matter with anti-dark matter, so as to leave almost purely asymmetric dark matter.

We showed that sufficient pre-annihilation in this scenario requires the scale of $G_*$ breaking $M_*$ to be below about $3 \times 10^6$ GeV.
The pre-annihilation mechanism requires that there be heavy unstable dark matter particles, an asymmetry in which is initially generated and which later decay to a light stable form of dark matter, which inherits this asymmetry. (These heavy unstable dark matter particles, which we call $\chi$, have mass $m_{\chi} \sim M_*$. The light stable dark matter particles, which we call $\zeta$, should have mass of order 5 GeV to obtain $\Omega_{DM} \sim 5 \Omega_B$.)
We have shown that these $\chi$ decays must have lifetimes that fall within a certain relatively narrow range.
Shorter lifetimes would cause these decays to regenerate the symmetric component
of dark matter, while longer lifetimes would result in the stable dark matter being too hot. The range for $T_{dec}$ (the temperature when the $\chi$ particles decay) is roughly between $2 \times 10^{-7}$ and $10^{-3}$ times $m_{\chi}$.
We have shown that the $\chi$ particles can have decay lifetimes in the required range,
if they come from dim-5 operators suppressed by $1/M_{GUT}$ or $1/M_{P \ell}$.

Finally, we have shown that this scenario can be implemented within the context of grand unified models based on groups of rank $\geq 6$. We presented a fairly simple $SU(7)$ model as an example. This model illustrates how several features of the scenario can arise naturally within grand unification: (1) Standard-Model-singlet fermions that can play the role of the heavy dark matter are contained in the fermion multiplets of the large unified group, thus unifying ordinary quarks and leptons with dark matter; (2) global symmetries that stabilize the dark matter can arise as accidental global symmetries, analogously to how $B-L$ arises as an accidental symmetry in minimal $SU(5)$; and (3) new (non-Standard Model) non-abelian gauge symmetries are left after the breaking of the grand unified group, which can play the role of $G_*$.

There are several model-building issues that deserve further study. One has to do with the breaking of the global symmetry $X$ that stabilizes dark matter. We have found that in fully gauge unified models that realize the scenario we discuss here, the global symmetry $X$ not only stabilizes dark matter but also stabilizes certain heavy vectorlike fermions that generically exist in such models. $X$ must therefore be slightly broken in such a way that the heavy vectorlike particles decay rapidly, while keeping the dark matter sufficiently stable. We showed one way to do this, but it would interesting to see if there are other better ways. Similarly, the way that the light dark matter particles $\zeta$ were given small mass in a ``technically natural" way could possibly be improved upon.

\section*{Acknowledgements} S.M.B. was
supported in part by the Department of Energy (DE-FG02-12ER41808).  R.J.S. was supported in part by the Department of Energy
(DE-SC0011981).


\section{Introduction}
Image denoising as a fundamental problem in image restoration has been extensively studied in image processing and computer vision (see [1]-[12] for few various representative works).
In this paper, we address the problem of image restoration from \emph{mixed Gaussian-impulse noise} data. Suppose the clean image $u\in \Bbb{R}^n$ is corrupted by both additive white Gaussian noise (AWGN) of $\eta\in \Bbb{R}^n$ (with variance $\sigma^2$) and an impulse noise. Also, let ${\Bbb{N}}_{imp}$ represent the process of image degradation with impulse noise. Accordingly, the overall degradation by mixed Gaussian-impulse noise can be modeled as the following:
\begin{equation}
f= {\Bbb{N}}_{imp}(u+\eta),
\label{eq.1}
\end{equation}
where $f\in \Bbb{R}^n$ is the observed image (noisy measurements). In the literature, there are two common types of impulse noise: salt-and-pepper (SP) noise and random-valued (RV) noise. Let denote $u_{a,b}$ as the pixel value of image $u$ at location $(a,b)\in {\cal{I}}$, and $[d_{\min}, d_{\max}]$ be the dynamic range of $u$. The operator ${\Bbb{N}}_{imp} (\cdot)$  is defined as follows:
\begin{equation}
{\tilde{f}}_{a,b}={\Bbb{N}}_{imp}(u_{a,b})=\left \{\begin{array}{ccc} d_{a,b},& {\rm{with~probability}}~{{r}/{2}}, \\ d_{a,b}^{\prime}, &  {\rm{with~probability}}~{{r}/{2}}, \\ u_{a,b},& {\rm{with~probability}}~{1-r}, \end{array}\right.
\label{eq.2}
\end{equation}
where $d_{a,b}$ and $d_{a,b}^{\prime}$ are in the range of $[d_{\min}, d_{\max}]$ and independent from $u_{a,b}$. When $d_{a,b}=d_{\min}$ and $d_{a,b}^{\prime}=d_{\max}$, the model is SP noise; however, if $d_{a,b}=d_{a,b}^{\prime}$ is identically and uniformly distributed random numbers in $[d_{\min}, d_{\max}]$ with probability $r$, the model is called RV noise. Median-type filters, e.g. adaptive median filter (AMF) [13] and adaptive center-weighted median filter (ACWMF) [14], are widely used in the literature to remove impulse noise.

Let $\cal{N}$ denote the set of locations of noisy pixel candidates corrupted by impulse noise, and accordingly, the set of pixels that are likely to be uncorrupted by impulse noise is defined as $\cal{A}=\cal{I}/\cal{N}$. Therefore, the degradation matrix $\Phi$ can be defined as $\Phi_{a,b}:=\{1:(a,b)\in{\cal{A}}~or~0:(a,b)\in{\cal{N}}\}$ \footnote{In this paper, similar to [5], [7]-[12] and for the sake of a fair comparison, $\Phi$ is determined by AMF and ACWMF for SP and RV noise detection, respectively.}.
More precisely, the pixels probably being corrupted by impulse noise are assigned with $0$ while the rest of pixel values are more likely to be corrupted by AWGN.
Delving further into this issue reveals this fact that, $\Phi\circ f$ follows an approximate Gaussian distribution, where $\circ$ denotes the element-wise multiplication between two matrices.
Also, due to the \emph{ill-posed} nature of (1), regularization-based methods are often used to convert (1) into a \emph{well-posed} problem by imposing \emph{a priori} information of the underlying signal [7], [9]-[12].
Hence, from the Bayesian inference, the following regularization-based framework for mixed Gaussian-impulse noise removal is written as:
 \begin{align}
\mathop{\min }_{u} {1 \over 2}\|\Phi\circ(f-u)\|_{\ell_2}^2 +\lambda{\cal{J}}(u),
\label{eq.33}
\end{align}
where the first term is a penalty that represents the closeness of the solution to the observed scene; the second term is a regularization term which represents a priori information of the original scene; and $\lambda$ is a regularization parameter that balances the contribution of both terms. The term ${\cal{J}}(u)$ could be various choices, e.g., Tikhonov regularization [15], Geman and Reynolds' half quadratic variational model [16], Rudin, Osher and Fatemi's total variation model [17], Mumford-Shah model [18], and framelet based model [7].
These approaches do not need to find the damaged pixels and perform well in impulse noise removal. However, for the case of images corrupted by mixed Gaussian impulse noise, the Gaussian noise is not treated properly [12]. In recent works, the local smoothness and the nonlocal self-similarity properties (or image sparsity prior) exhibited in natural images are combined into the final cost functional of image restoration solution to achieve better performance, e.g., see [10]-[12].

The main contributions of this paper are listed as follows. First, we propose an \emph{adaptive curvelet thresholding} (ACT) criterion, in a statistical manner, for adaptively characterizing the perturbations appeared during noise removal process---exploiting local consistency. Second, we establish a new statistical regularization term, called \emph{joint adaptive statistical prior} (JASP) which enforces both the local and nonlocal statistical consistencies, simultaneously, in a unified manner. Third, we propose a novel technique for mixed Gaussian plus impulse noise removal using JASP in a variational scheme, we refer to it as De-JASP. To solve the above variational scheme, an efficient alternating minimization algorithm is developed based on split Bregman iterative framework [19].

The rest of this paper is organized as follows. Section II provides a brief background on a recently proposed modeling for nonlocal self-similarity. The proposed ACT is introduced in Section III. The proposed regularization term, along with its incorporation into the variational framework of mixed noise removal, and implementation details are introduced in Section IV. Numerical results are given in Section V and finally, Section VI concludes the paper.
\section{Background on Nonlocal Statistical Modeling (NLSM)}
Motivated by the success of BM3D [4] in image denoising, Zhang \emph{et al}. [11] proposed a model for nonlocal self-similarity prior information, which is employed efficiently in image restoration.
The NLSM explores the nonlocal self-similarity by means of the distribution of the transform coefficients, which are obtained by transforming the 3D array generated by stacking similar image patches.
To elucidate on, First, the image $u$ (of size $\sqrt{n}\times \sqrt{n}$) is divided into $P$ overlapped patches of equal sizes, i.e., $u_{p_{{l}^{\prime}}}$ (of size $\sqrt{B_p}\times \sqrt{B_p}$) at location ${l}^{\prime}$, where $l^{\prime}=1,2,\cdots,P$. Then, for each patch the best $c$ matched patches are found within a searching window (of size $\omega\times\omega$). Next, these similar patches are stacked into a 3D array, $G_{u_{p_{{l}^{\prime}}}}$, which is called a group. Now, a 3D transform, ${\cal{T}}^{3D}(\cdot)$, is conducted on the group to obtain its transform coefficients, followed by arranging them in a lexicographic order.
The mathematical formulation of the NLSM for self-similarity in transform domain can be written as:
\begin{equation}
\Psi_{{}_{\rm{NLSM}}}(u)=\|\Theta_u\|_{\ell_1}={\sum_{l^{\prime}=1}^{P}}\|{\cal{T}}^{3D}(G_{u_{p_{l^{\prime}}}})\|_{\ell_1},
\label{eq.patch}
\end{equation}
where $\Psi_{{}_{\rm{NLSM}}}(\cdot)$ corresponds to the NLSM operator.
To put it simply and briefly, after obtaining $\Theta_u$, the new estimate of $u$ is achieved by $\hat{u}=\Omega_{{}_{\rm{NLSM}}}(\Theta_u)$, where $\Omega_{{}_{\rm{NLSM}}}(\cdot)$ corresponds the inverse operator of $\Psi_{{}_{\rm{NLSM}}}$.
\section{The Proposed Adaptive Curvelet Thresholding (ACT)}
\begin{figure*}[t!]
{\hspace{-5.6cm}
         {\centering{
            \begin{subfigure}[b]{1\textwidth}
                 \includegraphics[scale=0.068]{histogram-eps-converted-to.pdf}
                  \centering

         \end{subfigure}\\}}}%
         {\footnotesize{Fig. 1: (a)-(b) The histograms of two successive orientations at one scale of the curvelet coefficients of the clean image ``Fence" (which is shown in \emph{Fig. 2(a)}). (c)-(d) The histograms of the same orientations at the same scale of the curvelet coefficients of the initiated image $u_{init}$ (which is shown in \emph{Fig. 2(c)}), where, it is reconstructed by employing AMF on noisy image of \emph{Fig. 2(b)}. (e)-(f) PDFs of noise distribution (empirical perturbations) in natural log domain at the same orientations and scale (corresponding to the histograms of \emph{Figs. 1(a)-1(d)}) along with their best fitted PDFs: the noise distribution (blue), the Gaussian fit (red), and the Laplacian fit (black).}}
\end{figure*}%
The curvelet transform [20] is a multi-scale directional transform which provides a new multi-resolution representation with several features that are superior to existing representation, i.e., wavelets and steerable pyramids, thereby gaining good
performance in image denoising [1], [21]. For these reasons, here, an adaptive denoising/shrinkage procedure in curvelet domain via estimation of clean curvelet coefficients from available rough estimated data is proposed

We assume that the initially estimated image $u_{init}$ achieved by the median filter (AMF or ACWMF) is perturbed and model that perturbation by additive noise, with no prejudging on being AWGN or any specific species of noise. To elucidate on, this modeling of perturbations by additive noise helps us to model the noise probability density function (PDF) effectively and being able to solve the ensuing optimization problem efficiently.
Suppose, the initiated image coefficients in curvelet domain (${\hat{\vartheta}}=\Psi u_{init}$)\footnote{In matrix notation, we can write $\vartheta=\Psi u$, and $u=\Psi^{\dag}\vartheta$, where $\vartheta$ is the curvelet coefficient of function $u$, $\Psi$ is a forward curvelet transform and $\Psi^{\dag}=(\Psi^{T}\Psi)^{-1}\Psi^{T}$ is its Moore-Penrose pseudo-inverse transform (corresponding to the minimal dual synthesis frame). Due to the Parseval tight frame property of curvelet transform (i.e., $\Psi^{T}\Psi=I$), we have $\Psi^{\dag}=\Psi^{T}$.} can be written as:
\begin{equation}
{\hat{\vartheta}}_{j,l,{\rm{\bf{k}}}}={\vartheta}_{j,l,{\rm{\bf{k}}}}+n_{j,l,{\rm{\bf{k}}}},
\label{eq.curvecoeff}
\end{equation}
where ${\hat{\vartheta}}$ and ${{\vartheta}}$ are the initiated image curvelet coefficient (available rough coefficient which is assumed to be noisy) and the noise-free (or clean) curvelet coefficients, respectively, and $n$ is unknown additive noise (or \emph{perturbation}) vector. Also, the three indices of $j,l,{\rm{\bf{k}}}$ are curvelet triple index of \emph{scale}, \emph{rotation} and \emph{position}, respectively. The problem of image denoising in curvelet domain can be expressed as estimation of clean coefficients from noisy data with Bayesian estimation techniques. Either minimum mean-squared error (MMSE) or maximum \emph{a posteriori} (MAP) estimator is used for solving this problem, the solution requires a prior knowledge about the distribution of curvelet coefficients. More specifically, these methods are optimized w.r.t. the marginal statistics of the coefficients within each sub-band, by imposing a prior distribution on the clean transform coefficients [6].

Figs. 1(a)-1(b) show the histograms of two successive orientations at one scale of the curvelet coefficients of the clean image ``Fence" ($256\times 256$)---see Fig. 2(a).
Obviously, it can be observed that the distributions are characterized by a very sharp peak at zero amplitude and extended tails to both sides of the peak (\emph{leptokurtic distribution}). This leptokurtic property is observed in all histograms of each scale and orientation of all test images.
In this paper, Laplacian distribution is chosen to model the marginal distributions of curvelet coefficients of the clean image, i.e., ${\Pr}_{{}_{{{\vartheta}}}}({{\vartheta}}_{j,l,{\rm{\bf{k}}}})=({{1}/{{\sqrt{2}}\sigma_{j,l,{\rm{\bf{k}}}}}})\cdot\exp{\big({{-{\sqrt{2}}~|{{\vartheta}}_{j,l,{\rm{\bf{k}}}}|}/{{\sigma_{j,l,{\rm{\bf{k}}}}}}}\big)}$,
where $\sigma_{j,l,{\rm{\bf{k}}}}$ is the standard deviation of clean coefficients.
This choice is despite the fact that marginal distributions of curvelet coefficients have significantly heavier tails than a Laplacian distribution; tend to go zero slower than a Gaussian distribution; and can be well modelled by a hyper-Laplacian distribution.
But, actually, this choice makes a trade-off between modeling the coefficients statistics accurately and being able to solve the ensuing optimization problem efficiently.

The histograms of the same two successive orientations at the same scale of the coefficients of the initially estimated image $u_{init}$---see Fig. 2 (c)---are shown in Figs. 1(c)-1(d). More specifically, these histograms are obtained using the curvelet coefficients of the initiated image $u_{init}$, which is achieved by the initial median filtering step. As can be seen, here, the leptokurtic behavior and heavy tails do not exist in the shown distributions. More accurately, these two properties do not exist in all distributions of each scale and orientation of all initiated images achieved by this step.

Considering (5), the noise distribution can be obtained using the curvelet coefficients of the clean image and the initiated image. Figs. 1(e)-1(f) show the PDFs of noise distribution at the same two successive orientations of the same scale corresponding to those of Figs. 1(a)-1(d). It can be observed that the empirical distributions of noise coefficients $n_{j,l,{\rm{\bf{k}}}}$ can be well characterized by Gaussian distributions, while the Laplacian distributions have much larger fitting errors.
This observation motivates us to model the PDF of $n_{j,l,{\rm{\bf{k}}}}$ by a standard Gaussian distribution: ${\Pr}_{{}_{n}}({n}_{j,l,{\rm{\bf{k}}}})=({{1}/{{\sqrt{2\pi}}\sigma_{{n}_{j,l}}}})\cdot\exp{\big({{-{n}_{j,l,{\rm{\bf{k}}}}^{2}}/{{2\sigma_{{n}_{j,l}}^{2}}}}\big)}$, where $\sigma_{{n}_{j,l}}^{2}$ is the noise variance in curvelet domain, which can be estimated using robust median estimator [22] and Monte-Carlo simulations [20].

\begin{figure}[t!]
{\hspace{1.2cm}
         {\centering
         \begin{subfigure}[b]{0.21\textwidth}
                 \includegraphics[scale=0.35]{fence-eps-converted-to.pdf}
                  \centering
                 \caption{}
                 \label{fig:1}
         \end{subfigure}%

         \begin{subfigure}[b]{0.201\textwidth}
                 \includegraphics[scale=0.35]{F_noisy-eps-converted-to.pdf}
                 \centering
                 \caption{}
                 \label{fig:1}
         \end{subfigure}

         \begin{subfigure}[b]{0.201\textwidth}
                 \includegraphics[scale=0.35]{F_init-eps-converted-to.pdf}
                  \centering
                 \caption{}
                 \label{fig:1}
         \end{subfigure}%

         \begin{subfigure}[b]{0.21\textwidth}
                 \includegraphics[scale=0.35]{F_ACT-eps-converted-to.pdf}
                 \centering
                 \caption{}
                 \label{fig:1}
         \end{subfigure}}}\\
{\footnotesize{Fig. 2: (a) Original \emph{Fence} image (of szie $256\times 256$). (b) The noisy image which is corrupted by Gaussian plus salt-and-pepper impulse noise with $\sigma=30$ and $r_{sp}=30\%$ (PSNR=9.92 dB, SSIM=0.143). (c) The initiated image using AMF (PSNR=18.81 dB, SSIM=0.432). (d) The recovered image via ACT by 11 iterations (PSNR=22.11 dB, SSIM=0.589).}}
\end{figure}

Consider the marginal distribution of the clean curvelet coefficients. Since in image denoising problems, we do not have access to the clean image, the variance of the clean image curvelet coefficients ${{\sigma_{{j,l,{\rm{\bf{k}}}}}^{2}}}$ can be estimated by \emph{maximum likelihood} (ML) estimator, through averaging over the neighboring initiated image coefficients (available ${\hat{\vartheta}}_{j,l,{\rm{\bf{k}}}}$)  within a squared window [23].

In order to estimate ${{\vartheta}_{j,l,{\rm{\bf{k}}}}}$ from the available observation ${\hat{\vartheta}_{j,l,{\rm{\bf{k}}}}}$, the MAP estimator is employed:
\begin{equation}
{\vartheta}_{j,l,{\rm{\bf{k}}}}=\arg\max_{{\vartheta}_{j,l,{\rm{\bf{k}}}}}\bigg\{{{\Pr}_{{}_{{{\vartheta}}|\hat\vartheta}}}\big({\vartheta_{j,l,{\rm{\bf{k}}}}}|{{\hat{\vartheta}}_{j,l,{\rm{\bf{k}}}}}\big)\bigg\}.
\label{eq.map}
\end{equation}
By considering (5) and using Bayes' rule, we get:
\begin{align}
{\vartheta}_{j,l,{\rm{\bf{k}}}}=\arg\max_{{\vartheta}_{j,l,{\rm{\bf{k}}}}}\bigg\{{{\Pr}_{{}_{{\hat{\vartheta}}|\vartheta}}}\big({\hat\vartheta_{j,l,{\rm{\bf{k}}}}}|{\vartheta_{j,l,{\rm{\bf{k}}}}}\big)\cdot{{\Pr}_{{}_{{\vartheta}}}\big({\vartheta_{j,l,{\rm{\bf{k}}}}}\big)}\bigg\} \nonumber \\
=\arg\max_{{\vartheta}_{j,l,{\rm{\bf{k}}}}}\bigg\{{{\Pr}_{{}_{n}}}\big({\hat\vartheta_{j,l,{\rm{\bf{k}}}}}-{\vartheta_{j,l,{\rm{\bf{k}}}}}\big)\cdot{{\Pr}_{{}_{{\vartheta}}}\big({\vartheta_{j,l,{\rm{\bf{k}}}}}\big)}\bigg\}.
\label{eq.bays}
\end{align}
Therefore, (7) allows us to write this estimation in terms of the PDF of the noise coefficient ($\Pr_{{}_{n}}$), and the PDF of the clean image coefficient ($\Pr_{{}_{\vartheta}}$).
By substituting $\Pr_{{}_{n}}$ and $\Pr_{{}_{\vartheta}}$ into (7) and doing some manipulations, the estimation of $\vartheta_{j,l,{\rm{\bf{k}}}}$ is achieved as:
\begin{equation}
{\vartheta}_{j,l,{\rm{\bf{k}}}}={\rm{sign}}(\hat\vartheta_{j,l,{\rm{\bf{k}}}})\max\big({0,|\hat\vartheta_{j,l,{\rm{\bf{k}}}}|-{\frac{\sqrt{2} \sigma_{{n}_{j,l}}^{2}}{\sigma_{j,l,{\rm{\bf{k}}}}}}}\big),
\label{eq.varshrink}
\end{equation}
where, in fact, (8) is the classical soft shrinkage function [24] defined as ${\cal{S}}_{\rho}(y):={\rm{sgn}}(y)\cdot \max(0, |y|-\rho)$.
Not only can the proposed ACT applied for the initiated image achieved by the initial median filtering step, but also it can be employed in further iterations---by using the last approximated image achieved in previous iteration as the noisy image for the current iteration---which leads to a simple and effective mixed Gaussian-impulse noise removal scheme, e.g., see Fig. 2(d).

\section{The Proposed Image Denoising via JASP (De-JASP)}
The proposed \emph{joint adaptive statistical prior} (JASP) is defined by integrating both the local information prior (which depicts the local smoothness and geometric regularity of image structures achieved by discrete curvelet transform (DCuT) [20]) and the nonlocal self-similarity prior (corresponding to the NLSM in 3D transform domain achieved by the method introduced in Section II). In a mathematical expression, the proposed JASP is expressed as:
\begin{align}
{\cal{J}}_{{}_{\rm{JASP}}}(u)&=\Psi_{{}_{\rm{DCuT}}}(u)+{{\tau}}\Psi_{{}_{\rm{NLSM}}}(u) \nonumber \\
&=\|\Psi u\|_{\ell_1}+{{\tau}}\|\Theta_u\|_{\ell_1},
\label{eq.jasr}
\end{align}
where the first and the second terms represent the image local smoothness prior and nonlocal self-similarity prior, respectively. Also, ${{\tau}}$ is a regularization parameter which controls the trade-off between two competing (statistical) terms.

By substituting the proposed JASP (9) for the regularization term ${\cal{J}}(u)$ in the regularization-based framework of (3), and by introducing two auxiliary variables $\vartheta$ and $w$, the proposed optimization problem for image recovery is expressed as follow:
 \begin{align}
\mathop{\min }_{u,\vartheta, w} \{{1 \over 2}\|\Phi\circ(f-u)\|_{\ell_2}^2 +\lambda\|\vartheta\|_{\ell_1}+{\tau}^{\prime} \|\Theta_w\|_{\ell_1}\} \nonumber \\
s.t.~~~~\vartheta=\Psi u,~ u=w,~ w=\Omega_{{}_{\rm{NLSM}}}(\Theta_w).
\label{eq.constrained}
\end{align}
In order to solve the above minimization problem, an alternating split Bregman iterative algorithm [19] is invoked. We finally obtain the following schemes:
   \begin{subequations}\label{eq:subeqns}
     \begin{align}
(&u^{k+1}, \vartheta^{k+1}, w^{k+1})\leftarrow\arg\min_{u, \vartheta, w}\Big\{{{1\over 2}{\|\Phi\circ(f-u)\|_{\ell_2}^2}+{\lambda}\|\vartheta\|_{\ell_1}}\nonumber\\&{+{\tau}^{\prime} \|\Theta_w\|_{\ell_1}+{\frac{\mu_1}{2}}{\|\vartheta-\Psi u-b^{k}\|_{\ell_{2}}^2}}
+{\frac{\mu_2}{2}}{\|u-w-c^{k}\|_{\ell_{2}}^2}\Big\} \label{eq:subeq1}\\
&b^{k+1}\leftarrow b^{k}-\vartheta^{k+1}+\Psi u^{k+1}\label{eq:subeq4}\\
&c^{k+1}\leftarrow c^{k}-u^{k+1}+w^{k+1}.\label{eq:subeq5}
     \end{align}
    \end{subequations}
Here, $\mu_1$ and $\mu_2$ are fixed value parameters for improving the numerical stability of the algorithm.

Given $\vartheta^{k}$ and $w^k$, the $u$ \emph{sub-problem} of (11a) consists of minimizing a strictly convex quadratic function that can be solved easily, i.e., $u^{k+1}={\big(\Phi^{T}\Phi+\mu I\big)}^{-1}\circ {\cal{D}}$, where $\mu=\mu_1+\mu_2$ and ${\cal{D}}=\big(\mu_1\Psi^{T}(\vartheta^{k}-b^{k})+\mu_2(w^{k}+c^{k})+\Phi^{T}\Phi\circ f\big)$. In order to avoid computing the matrix inversion in the first term, the matrix inversion lemma of the Sherman-Morrison formula is applied to that term. Thus, the $u$ \emph{sub-problem} can be obtained as:
\begin{equation}
u^{k+1}={\frac{1}{\mu}}\big({I-\frac{\Phi^{T}\Phi}{\mu+\Phi\Phi^{T}}}\big)\circ {\cal{D}}.
\label{eq.usubp}
\end{equation}

By given $u$ in hand (according to (12)) and considering $\hat\vartheta=\Psi u+b$ (for simplicity, the superscript $k$ is dropped without confusion), the $\vartheta$ \emph{sub-problem} of (11a) becomes:
\begin{equation}
\min_{\vartheta}{{1\over 2}{\|\vartheta-\hat\vartheta\|_{\ell_{2}}^{2}+{\lambda \over \mu_1}{\|\vartheta\|_{\ell_{1}}}}}.
\label{eq.coeffsubp}
\end{equation}
By these transformations, we regard $\vartheta^{\prime}$ as some type of the noisy unknown coefficient $\vartheta$. More precisely, the sub-problem of (13) can be interpreted as the denoising in curvelet coefficient domain.
Considering the fact that the unknown variable $\vartheta$ is component-wise separable in curvelet domain, each of its component ${\vartheta}_{j,l,{\rm{\bf{k}}}}$ can be obtained by a component-wise shrinkage procedure:
\begin{equation}
{\vartheta}_{j,l,{\rm{\bf{k}}}}=\arg\min_{{\vartheta}_{j,l,{\rm{\bf{k}}}}}\bigg\{{\frac{1}{2}}{\big|}{\vartheta}_{j,l,{\rm{\bf{k}}}}-{\hat\vartheta}_{j,l,{\rm{\bf{k}}}}{\big|}^{2}+{\frac{\lambda_{j,l,{\rm{\bf{k}}}}}{\mu_1}}{\big|}{\vartheta}_{j,l,{\rm{\bf{k}}}}{\big|}\bigg\},
\label{eq.cws}
\end{equation}
where its solution is the same as the one introduced in (8) with regarding ${\frac{\lambda_{j,l,{\rm{\bf{k}}}}}{\mu_1}}={\frac{{\sqrt{2}}\sigma_{n_{j,l}}^{2}}{\sigma_{j,l,{\rm{\bf{k}}}}}}$. To elucidate on, the $\vartheta$ \emph{sub-problem} of (13) can be interpreted as the proposed ACT introduced in Section III. Thus, the parameter $\lambda$ is self-adaptive.

Given $u$, analogously, the $w$ \emph{sub-problem} of (11a) can be written as:
\begin{equation}
\min_{w}{{1\over 2}{\|w-z\|_{\ell_{2}}^{2}+{{\tau}^{\prime} \over \mu_2}{\|\Theta_w\|_{\ell_{1}}}}},
\label{eq.wsubp}
\end{equation}
where $z=u-c$ (for simplicity the superscript $k$ is omitted without confusion). Owing to the complicated definition of $\Theta_w$, it seems difficult to solve (15) directly. In order to solve (15) amenable, according to [11], a reasonable assumption is used which leads to obtain a closed-form solution of (15). By assuming that all elements of $w-z$ (such that $w, z\in \Bbb{R}^{n}$) are independent and identically distributed (i.i.d.) with zero-mean and variance ${\sigma^{\prime}}^{2}$ \footnote{It is worth emphasizing that the above assumption does not need to be Gaussian, Laplacian or generalized Gaussian process.}, and by invoking the \emph{law of large numbers} in probability theory, for any $\epsilon^{\prime}>0$, it leads to:
\begin{equation}
\lim_{n\rightarrow\infty}\Pr\{{\big|}{{1}\over{n}}\|w-z\|_{\ell_{2}}^{2}-{\sigma^{\prime}}^{2}{\big|}<{{\epsilon^{\prime}}\over{2}}\}=1.
\label{eq.lln1}
\end{equation}
\begin{table}[!t]
\renewcommand{\arraystretch}{1.2}
{\footnotesize{{~~~~~~~~~~~~~~~~~TABLE I: T\sc he Proposed De-JASP Framework}}}\\\\
\label{table_example}
\centering
\begin{tabular}{l}
\hline\hline
 {\bfseries{Input}}: $f$, $B_p, \omega, {\tau}^{\prime}, \mu_1, \mu_2, Tol$, $k_{max}^{\prime}$: maximum iteration number\\
 {\bfseries{Output}}: $u^*$: Recovered image;\\
 {\bfseries{Initialization}}: Set $k=0$; $(b^{0},c^{0})=(0,0)$\\
$u_{init}=u^{0}$= median-filtering$(f)$; compute $\Phi$ \\
 {\bfseries{While}} a stop criterion is not satisfied  {\bfseries{do}}\\
1.~~~~~$u^{k+1}=\tilde{u}=u^{k}$;~${\hat{\vartheta}}^{k+1}=\Psi u^{k+1}+b^{k}$;~$z^{k+1}=u^{k+1}-c^{k}$\\
2.~~~~~solve $\vartheta$ \emph{sub-problem} using (14) and concatenating all coefficients\\
3.~~~~~solve $w$ \emph{sub-problem} using (19)\\
4.~~~~~solve $u$ \emph{sub-problem} using (12)\\
5.~~~~~compute $s^{k+1}$= SSIM($u^{k+1}, \tilde{u}$)\\
6.~~~~~compute \emph{diff} = {\rm {abs}}$(s^{k+1}-s^{k})$\\
7.~~~~~update $b^{k}, c^{k}$ using (11b) and (11c), respectively\\
8.~~~~~$k\leftarrow k+1$\\
{\bfseries{end}} While\\
\hline
{\bfseries{stopping criterion}}: $k=k_{max}^{\prime}$ or {\emph{diff}} $\leq Tol$$$.\\ \hline
\end{tabular}
\end{table}
Due to the orthogonal property of 3D transform ${\cal{T}}^{3D}$, and denoting the error vector in 3D transform domain by $\Theta_e=\Theta_w-\Theta_z$ (such that $\Theta_w, \Theta_z\in \Bbb{R}^{K_{\Theta}}$, where ${K_{\Theta}}=B_p*c*P$), it is concluded that all elements of $\Theta_e$ are i.i.d. with zero-mean and variance ${\sigma^{\prime}}^{2}$. Thanks to the \emph{law of large numbers}, by doing the same manipulations with (16) to $\Theta_e$, for any $\epsilon^{\prime}>0$, it yields:
\begin{equation}
\lim_{K_{\Theta}\rightarrow\infty}\Pr\{{\big|}{{1}\over{K_{\Theta}}}\|\Theta_w-\Theta_z\|_{\ell_{2}}^{2}-{\sigma^{\prime}}^{2}{\big|}<{{\epsilon^{\prime}}\over{2}}\}=1.
\label{eq.lln2}
\end{equation}
Therefore, according to (16) and (17) the following relationship is described (almost surely): ${\frac{1}{n}}\|w-z\|_{\ell_{2}}^{2}={\frac{1}{K_\Theta}}\|\Theta_w-\Theta_z\|_{\ell_{2}}^{2}$, where incorporating it into (15) leads to:
\begin{equation}
\arg\min_{\Theta_w}{{1\over 2}{\|\Theta_w-\Theta_z\|_{\ell_{2}}^{2}+{\frac{K_{\Theta} {\tau}^{\prime}}{\mu_2~n}}{\|\Theta_w\|_{\ell_{1}}}}}.
\label{eq.newwsub}
\end{equation}
Since the unknown variable $\Theta_w$ is component-wise separable in (18), each of its component $\Theta_w(j^{\prime}), j^{\prime}=1,\cdots,K_{\Theta}$ can be independently obtained by a component-wise (soft) shrinkage procedure. Thus, the closed-form solution of (15) is written as:
\begin{equation}
w=\Omega_{{}_{\rm{NLSM}}}\big({{\Theta_w}}\big)=\Omega_{{}_{\rm{NLSM}}}\big({{{\cal{S}}_{{\frac{K_{\Theta} {\tau}^{\prime}}{\mu_2n}}}\big(\Theta_z\big)}}\big).
\label{eq.omega}
\end{equation}

The proposed De-JASP algorithm is delineated in Table I.

\section{Experimental Results}
In this section, we evaluate the performance of the proposed method---De-JASP.
The performances of our experiments are evaluated on 4 gray-scale test images of \emph{Barbara}, \emph{Boat}, \emph{Fence}, and \emph{Parrot} (of size $256\times 256$). To evaluate our simulation results, we use two applicable quality assessors, the pick signal-to-noise ratio (PSNR) in dB and the structural similarity (SSIM) [25].
All the experimental results reported in this paper are averaged over 10 independent trials. Also, the second generation of DCuT via wrapping [20] is employed in our implementations. To demonstrate the effectiveness of the proposed method in image denoising, we have compared it with two competitive mixed noise removal techniques: (i) iterative framelet-based approximation/sparsity deblurring algorithm IFASDA [7], and (ii) image denoising via joint statistical modeling (JSM) [11].
It is worth emphasizing that these two recently proposed denoising methods are among the state-of-the-art techniques  in mixed-noise removal which achieve the best performance so far.
Note that, to make a fair comparison among the competing methods, we have carefully tuned their parameters to achieve the best performance. Also, for the sake of fair comparison, the same test conditions are used in all experiments; i.e., the same noisy measurements are applied for each method. Due to the space limitations, only parts of the experimental results are shown in this paper. Interested readers may contact the corresponding author for all the images and the source code.

In the following, the effectiveness of De-JASP is investigated. In our implementation, in process of finding NLSM for self-similarity in 3D transform domain, the size of each patch is set to $8\times 8$. The size of training window for searching matched blocks is set to $20\times 20$, and the number of best matched blocks is set to $10$. The other main parameters of the proposed algorithm are empirically set: $k_{max}^{\prime}=20$, $Tol=10^{-3}$, $\mu_1=0.16\mu$, $\mu_2=0.84\mu$ and $\mu=2.7\times 10^{-3}$. Also, the orthogonal 3D transform denoted by ${\cal{T}}^{3D}$ is composed of 2D DCT and 1D Haar transform.
\begin{table*}[!t]
\renewcommand{\arraystretch}{1}
{\footnotesize{\footnotesize{~~~~~~~~~~~~~~~~~TABLE II:~\sc{PSNR/SSIM Comparisons for Gaussian plus Salt-and-Pepper Noise Removal.}}}}\\~\\
\label{table_example}
\centering
\begin{tabular}{|c|c|c|c|c|c|c|c|}\cline{6-8}
\multicolumn{5}{c|}{}&\multicolumn{3}{c|}{Algorithm}\\ \cline{6-8}\hline
Test Image&$r_{sp}$&$\sigma$&Noisy Image& Init. Image&IFASDA [7]&JSM [11]& De-JASP\\ \hline
{}&{0.2}& 20&12.10/0.147&22.48/0.529 &24.92/0.756&27.73/0.795&\bf{28.06/0.835}\\
{Barbara}&{0.2} &30&11.69/0.133 &19.57/0.385& 20.46/0.514&22.76/0.537&\bf{24.83/0.676}\\
{}&{0.3} &30& 10.21/0.091& 19.10/0.364& 20.34/0.502& 22.63/0.533 &\bf{24.43/0.655}\\ \hline\hline
{}&{0.2}& 20&12.14/0.134 &23.08/0.510 &27.81/0.771 & 28.40/0.784&\bf{29.22/0.833}\\
{Boat}&{0.2} & 30 &11.72/0.122 &19.93/0.373 & 22.74/0.506&23.22/0.520 &\bf{25.88/0.679}\\
{}&{0.3}& 30& 10.24/0.085& 19.47/0.355&22.67/0.499 & 23.15/0.518&\bf{25.47/0.659}\\ \hline\hline
{}&{0.2}& 20&11.77/0.212 & 21.58/0.560&24.99/0.782 & 26.60/0.791&\bf{26.79/0.812}\\
{Fence}&{0.2} & 30 &11.39/0.197 &19.13/0.443 & 21.03/0.577&22.35/0.583 &\bf{24.15/0.693}\\
{}&{0.3}& 30& 9.92/0.143& 18.81/0.432& 20.82/0.569&22.16/0.577 &\bf{23.80/0.679}\\ \hline\hline
{}&{0.2}& 20&11.89/0.095 & 23.15/0.407& 27.78/0.759& 28.83/0.765&\bf{29.65/0.837}\\
{Parrot}&{0.2} & 30 &11.50/0.087 &19.97/0.277 & 22.65/0.424& 23.51/0.438&\bf{26.21/0.640}\\
{}&{0.3}& 30& 10.02/0.063&19.51/0.264 &22.46/0.411 &23.31/0.435 & \bf{25.61/0.616}\\ \hline
\end{tabular}
\end{table*}

Table II lists the PSNR and SSIM values of the recovered noisy test images obtained by the proposed De-JASP compared to those obtained by JSM and IFASDA in different noise level hypotheses (best results are emphasized in bold face). From Table II, it can be inferred that the performance of the proposed De-JASP is superior to those of compared algorithms in terms of both employed objective quality assessors. The average PSNR (SSIM) gain of De-JASP over JSM and IFASDA, in mixed \emph{AWGN+RV} noise removal, can be up to 2.31 dB (0.154) and 3.55 dB (0.167), respectively.
For visual comparison, some recovered images by the competing methods, in different noise levels, are shown in Figs. 4-7. It can be observed that De-JASP outperforms the other competing methods and reproducing clearer images. Better visual comparison can be made by zooming the images on the screen.

\begin{figure}
{\hspace{1.2cm}
         {\centering
         \begin{subfigure}[b]{0.21\textwidth}
                 \includegraphics[scale=0.35]{B_noisy-eps-converted-to.pdf}
                  \centering
                 \caption{}
                 \label{fig:1}
         \end{subfigure}%

         \begin{subfigure}[b]{0.201\textwidth}
                 \includegraphics[scale=0.35]{B_IF-eps-converted-to.pdf}
                 \centering
                 \caption{}
                 \label{fig:1}
         \end{subfigure}

         \begin{subfigure}[b]{0.201\textwidth}
                 \includegraphics[scale=0.35]{B_JSM-eps-converted-to.pdf}
                  \centering
                 \caption{}
                 \label{fig:1}
         \end{subfigure}%

         \begin{subfigure}[b]{0.21\textwidth}
                 \includegraphics[scale=0.35]{B_JASP-eps-converted-to.pdf}
                 \centering
                 \caption{}
                 \label{fig:1}
         \end{subfigure}}}\\
{\footnotesize{Fig. 4: Visual quality comparison of mixed Gaussian-impulse noise removal on image ``Barbara". (a) noisy image corrupted by AWGN+SP noise with $\sigma=20$ and $r_{sp}=20\%$. Denoised results by (b) IFASDA [7]; (c) JSM [11]; and (d) the proposed De-JASP. For numerical comparison, see Table II.}}
\end{figure}
\begin{figure}
{\hspace{1.2cm}
         {\centering
         \begin{subfigure}[b]{0.21\textwidth}
                 \includegraphics[scale=0.35]{Bo_noisy-eps-converted-to.pdf}
                  \centering
                 \caption{}
                 \label{fig:1}
         \end{subfigure}%

         \begin{subfigure}[b]{0.201\textwidth}
                 \includegraphics[scale=0.35]{Bo_IF-eps-converted-to.pdf}
                 \centering
                 \caption{}
                 \label{fig:1}
         \end{subfigure}

         \begin{subfigure}[b]{0.201\textwidth}
                 \includegraphics[scale=0.35]{Bo_JSM-eps-converted-to.pdf}
                  \centering
                 \caption{}
                 \label{fig:1}
         \end{subfigure}%

         \begin{subfigure}[b]{0.21\textwidth}
                 \includegraphics[scale=0.35]{Bo_JASP-eps-converted-to.pdf}
                 \centering
                 \caption{}
                 \label{fig:1}
         \end{subfigure}}}\\
{\footnotesize{Fig. 5: Visual quality comparison of mixed Gaussian-impulse noise removal on image ``Boat". (a) noisy image corrupted by AWGN+SP noise with $\sigma=10$ and $r_{sp}=60\%$. Denoised results by (b) IFASDA [7] (PSNR=27.58 dB, SSIM=0.769); (c) JSM [11] (PSNR=28.06 dB, 0.793); and (d) the proposed De-JASP (PSNR=28.37 dB, SSIM=0.814).}}
\end{figure}

\begin{figure}
{\hspace{1.2cm}
         {\centering
         \begin{subfigure}[b]{0.21\textwidth}
                 \includegraphics[scale=0.35]{P_noisy-eps-converted-to.pdf}
                  \centering
                 \caption{}
                 \label{fig:1}
         \end{subfigure}%

         \begin{subfigure}[b]{0.201\textwidth}
                 \includegraphics[scale=0.35]{P_IF-eps-converted-to.pdf}
                 \centering
                 \caption{}
                 \label{fig:1}
         \end{subfigure}

         \begin{subfigure}[b]{0.201\textwidth}
                 \includegraphics[scale=0.35]{P_JSM-eps-converted-to.pdf}
                  \centering
                 \caption{}
                 \label{fig:1}
         \end{subfigure}%

         \begin{subfigure}[b]{0.21\textwidth}
                 \includegraphics[scale=0.35]{P_JASP-eps-converted-to.pdf}
                 \centering
                 \caption{}
                 \label{fig:1}
         \end{subfigure}}}\\
{\footnotesize{Fig. 6: Visual quality comparison of mixed Gaussian-impulse noise removal on image ``Parrot". (a) noisy image corrupted by AWGN+RV noise with $\sigma=20$ and $r_{rv}=20\%$. Denoised results by (b) IFASDA [7] (PSNR=26.34 dB, SSIM=0.726); (c) JSM [11] (PSNR=26.96 dB, 0.789); and (d) the proposed De-JASP (PSNR=27.42 dB, SSIM=0.816)}}
\end{figure}
\begin{figure}
{\hspace{1.2cm}
         {\centering
         \begin{subfigure}[b]{0.21\textwidth}
                 \includegraphics[scale=0.35]{F_noisy-eps-converted-to.pdf}
                  \centering
                 \caption{}
                 \label{fig:1}
         \end{subfigure}%

         \begin{subfigure}[b]{0.201\textwidth}
                 \includegraphics[scale=0.35]{F_IF-eps-converted-to.pdf}
                 \centering
                 \caption{}
                 \label{fig:1}
         \end{subfigure}

         \begin{subfigure}[b]{0.201\textwidth}
                 \includegraphics[scale=0.35]{F_JSM-eps-converted-to.pdf}
                  \centering
                 \caption{}
                 \label{fig:1}
         \end{subfigure}%

         \begin{subfigure}[b]{0.21\textwidth}
                 \includegraphics[scale=0.35]{F_JASP-eps-converted-to.pdf}
                 \centering
                 \caption{}
                 \label{fig:1}
         \end{subfigure}}}\\
{\footnotesize{Fig. 7: Visual quality comparison of mixed Gaussian-impulse noise removal on image ``Fence". (a) noisy image corrupted by AWGN+SP+RV noise with $\sigma=10$, $r_{sp}=10\%$ and $r_{rv}=10\%$. Denoised results by (b) IFASDA [7] (PSNR=20.11 dB, SSIM=0.516); (c) JSM [11] (PSNR=21.79 dB, 0.622); and (d) the proposed De-JASP (PSNR=23.19 dB, SSIM=0.698)}}
\end{figure}
\section{Conclusion}
In this paper, first, an adaptive curvelet thresholding (ACT) criterion is introduced to adaptively characterize and abate the perturbations which are appeared during denoising process. Then, a novel joint adaptive statistical prior (JASP) and a new strategy for mixed Gaussian-impulse noise removal via JASP, called De-JASP, is proposed. Extensive experimental results clearly confirm that De-JASP significantly outperforms the current state-of-the-art mixed Gaussian-impulse noise removal techniques.

For future work, we would like to propose an effective nonlocal self-similarity modeling, and pursue this direction to extend our proposed statistical prior to other image restoration applications such as image deblurring and inpainting.

\section*{Acknowledgment}
The authors would like to express their gratitude to Prof. J. Fadili (CNRSENSICAEN-Universit\`{e} de Caen) and Dr. J. Zhang (Peking
University) for many fruitful discussions. They also would like to thank the authors of [7] and [11] for sharing the source code of their papers used in Section V.

\ifCLASSOPTIONcaptionsoff
  \newpage
\fi