# Influence of primary hair and plasma on intensity distribution of black hole
shadows

Vitalii Vertogradov  Physics Department, Herzen State Pedagogical University
of Russia, 48 Moika Emb., Saint Petersburg 191186, Russia SPb branch of SAO
RAS, 65 Pulkovskoe Rd, Saint Petersburg 196140, Russia
<EMAIL_ADDRESS>Parth Bambhaniya  International Center for Space and
Cosmology, Ahmedabad University, Ahmedabad 380009, Gujarat, India
<EMAIL_ADDRESS>Maxim Misyura  Department of High Energy and Elementary
Particles Physics, Saint Petersburg State University, University Emb. 7/9,
Saint Petersburg, 199034, Russia<EMAIL_ADDRESS>
###### Abstract

In this paper, we investigate the influence of primary hair ($l$) on the
shadows of hairy Schwarzschild and Reissner-Nordström black holes obtained
through gravitational decoupling. In the context of hairy Schwarzschild black
holes, $l$ either has no effect or consistently enlarges the photon sphere
radius. Notably, even when it violates the strong energy condition, it can
decrease the radius. For Reissner-Nordström black holes, an additional matter
field consistently expands the photon sphere radius, potentially reaching
$3M$, akin to the pure Schwarzschild case. Remarkably, we demonstrate that
black holes can exist even when overcharged ($Q^{2}>M^{2}$), casting shadows.
Specific intensity calculations reveal $l$ consistently reduces it in both
scenarios. Furthermore, we investigate the impact of pressureless plasma,
finding $l$ exerts a stronger influence on visible size than plasma. These
results can help in our understanding of theoretical models of black hole
shadows and can be tested by comparison with the images obtained by EHT
collaboration.

keywords: Black holes, Shadows, EHT, Plasma, No-hair theorem.

##  I Introduction

Black holes possess distinctive features, notably an event horizon, a one-way
membrane in spacetime where objects, including light, can only enter but not
escape. However, the existence of an event horizon lacks conclusive evidence.
The curvature of spacetime near a black hole is so strong that the light can
move along a circle. These light rings around black hole form a unstable
photon sphere which forms a black hole shadow Synge:1966okc , the images of
which have been obtained by the Event Horizon Telescope Collaboration for the
central object in M87 galaxy EventHorizonTelescope:2019dse ;
EventHorizonTelescope:2019uob  and our own galaxy
EventHorizonTelescope:2022xnr . These images show that a black hole can’t be
perceived as an exotic theoretical object anymore but it is a real
astrophysical one.

Spacetime around a supermassive compact object is generally believed to abide
by vacuum solutions of the Einstein field equations, which include
Schwarzschild or Kerr solutions. On the other hand, observations indicate
toward a dense concentration of matter surrounding the Galactic Center. As
consequence, the existence of a vacuum solution in this situation seems less
probable. Therefore, it becomes crucial to speculate about a matter
distribution around the center while looking at the physical signatures that
go along with it. This approach allows it conceivable to test the no-hair
theorem and the influence of plasma on black holes in the context of the
Galactic Center. A shadow is a very useful tool to test different black hole
models and modified theories of gravitation Vagnozzi:2022moj . It can also be
used as a standard cosmological ruler Tsupko:2019pzg ; Vagnozzi:2020quf .

When light from a distant source approaches a black hole, it can be captured
by the black hole. The border of escaping to the distant observer is called a
photon sphere. Thus, the observer should see the black spot against the light
source distribution. A comprehensive review of a black hole shadow can be
found in the paper Perlick:2021aok . Various properties of a black hole shadow
have been extensively studied by numerous authors Claudel:2000yi ; Hod:2012ax
; Khoo:2016xqv ; Decanini:2010fz ; Shoom:2017ril ; Cederbaum:2015fra ;
Johannsen:2013vgc ; Abdujabbarov:2015xqa ; Younsi:2016azx ; Perlick:2018iye ;
Tsupko:2021yca ; Bambhaniya:2021ybs ; Bambhaniya:2021ugr ; bib:vaidya1 ;
bib:vaidya2 ; Heydarzade:2023gmd ; bib:tsupko_first ; bib:tsupko_plasma ;
bib:tsupko_lensing ; bib:understanding ; bib:joshi_shadow . Recently, it has
been understood that naked singularities can also cast a shadow
bib:joshi_shadow ; Saurabh:2023otl ; Saurabh:2022jjv ; Solanki:2021mkt .
Therefore, it is not possible to distinguish between a naked singularity and a
black hole based on their shadows.

In general, individuals believe that a photon sphere is what causes the
shadows to form. But recently, the authors of Joshi:2020tlq  introduced a new
spherically symmetric null naked singularity solution to the Einstein field
equation that casts a shadow albeit not having a photon sphere. They then
derive the general conditions under which a shadow could form in the absence
of a photon sphere for both timelike and null naked singularities, where both
sorts of singularities fulfil all the energy conditions Joshi:2020tlq .

A famous fact regarding black holes is the ‘no-hair’ theorem, which states
that a black hole can be described using only three parameters, i.e. mass $M$,
electric charge $Q$ and spin $a$. However, it has shown that a black hole may
possess a soft hair bib:hok_hair . Recently, it has been understood that the
‘no-hair’ theorem can be evaded by using minimal geometrical deformations
bib:mgd1 ; bib:mgd2 ; bib:mgd3  and extended gravitational decoupling methods
bib:gd1 ; bib:gd2 ; bib:gd3 . Ovalle et all bib:bh1  have introduced a
generalization of a Schwarzschild black hole surrounded by an anisotropic
fluid, which can have primary hair due to extended gravitational decoupling.
This new solution has provoked a significant number of papers that have
further generalized this solution to hairy Kerr bib:hairy_kerr , Vaidya and
generalized Vaidya bib:vermax , regular hairy black holes bib:vermax2 ;
bib:ovalle_regular  and many others. Gravitational decoupling represents a
novel and powerful tool for obtaining new solutions to the Einstein equations.

Gravitational decoupling enables the derivation of new solutions of the
Einstein field equations that describe the exterior geometry of a black hole.
These compact objects have fascinating physical and geometrical properties,
and it would be observationally significant to distinguish them. Therefore, it
is crucial to study the observational aspects of these solutions and test them
with observational data to confirm or reject a model. This paper explores the
properties of the black hole shadow of hairy Schwarzschild and Reissner-
Nordstrom black holes obtained by gravitational decoupling. The impact of a
primary hair on the radius of the photon sphere and black hole shadow is
investigated. Additionally, the specific intensity is calculated, and it is
shown that a primary hair always reduces it compared to the usual black hole
solution. The influence of pressureless plasma on the size of a shadow is also
evaluated.

This paper is organized as follows. In sec. II, we consider the influence of a
primary hair on the black hole shadow in hairy Schwarzschild and Reissner-
Nordstrom spacetimes. In sec. III, we prove that a primary hair can lead to an
overcharged black hole with $Q^{2}>M^{2}$ ,and this black hole can also cast a
shadow. In sec. IV, we evaluate the specific intensity of these solutions and
show that a primary hair always decreases it. We investigate the influence of
pressureless plasma on a shadow in sec. V, and sec. VI is the discussions and
conclusions. The system of units which we use throughout the paper, is
geometrized one in which $G=c=1$. The metric signature is considered as
$-+++$.

##  II Shadow of hairy black holes

In this section, we will consider several hairy black hole models suggested in
the paper bib:bh1 . We express the line element in terms of the lapse function
$f=f(r)$

$ds^{2}=-fdt^{2}+f^{-1}dr^{2}+r^{2}d\Omega^{2},$ (1)

here $d\Omega^{2}=d\theta^{2}+\sin^{2}\theta d\varphi^{2}$ is the metric on
two-sphere.

The spacetime (1) admits two constants of motion, i.e. energy per unit mass
$E$

$E=f\frac{dt}{d\lambda}\,,$ (2)

and angular momentum per unit mass $L$

$L=r^{2}\sin^{2}\theta\frac{d\varphi}{d\lambda},$ (3)

where $\lambda$ is an affine parameter along the null geodesic. The considered
models are spherically symmetric. Therefore, without loss of generality, we
can consider the motion only in the equatorial plane $\theta=\frac{\pi}{2}$.
When discussing the black hole shadow, the primary objective is to locate the
unstable photon orbit, also known as the photon sphere. By using the eqn. (2),
(3) and the null geodesic condition, $g_{ik}u^{i}u^{k}=0$, one can obtains the
radial part of geodesic equation

$\left(\frac{dr}{d\lambda}\right)^{2}=-V_{eff},$ (4)

where

$V_{eff}=f\frac{L^{2}}{r^{2}}-E^{2}\,,$ (5)

is an effective potential. The radius $r=r_{ph}$ is the radius of photon
sphere, if the following conditions are held:

$V_{eff}(r_{ph})=0,~{}~{}V^{\prime}_{eff}(r_{ph})=0\,,$ (6)

where a prime denotes the derivative with respect to the radial coordinate
$r$. Using the second condition from (6), it can determine the radius of the
possible photon sphere:

$V^{\prime}_{eff}(r_{ph})=0\rightarrow
f^{\prime}(r_{ph})r_{ph}-2f(r_{ph})=0\,.$ (7)

By using the first condition (6) and introducing the impact parameter
$b=\frac{L}{E}$, one can obtain the visible size of the shadow,

$V_{eff}(r_{ph})=0\rightarrow b=\frac{r_{ph}}{\sqrt{f(r_{ph})}}\,.$ (8)

The observer will see the black spot of the radius $b$, where coordinates are,

$\displaystyle x$ $\displaystyle=$ $\displaystyle b\cos\psi,$ $\displaystyle
y$ $\displaystyle=$ $\displaystyle b\sin\psi,$ $\displaystyle 0$
$\displaystyle\leq$ $\displaystyle\psi<2\pi.$ (9)

###  II.1 model 1

The first model which we examine is the hairy Schwarzschild black hole, the
lapse function $f$ of which is given by

$f(r)=1-\frac{2m}{r}+\alpha e^{-\frac{r}{M}}\,,$ (10)

where the Schwarzschild mass $M$ is related to the new mass parameter $m$ by
the relation

$m=M+\frac{\alpha l}{2}\,.$ (11)

Here, $\alpha$ is a coupling constant, $l$ is related to a primary hair, and
it has the length dimension. The influence of these parameters on the hawking
temperature, geodesic motions, gravitational lensing and energy extraction has
been investigated in bib:thermo ; bib:kudr_thermo ; bib:geod ; bib:lens ;
bib:energy . To find out what impact these parameters have on the shadow of a
black hole, we use eqn. (10) in the second condition of eqn. (6),

$V^{\prime}_{eff}(r_{ph})=0\rightarrow\frac{6m}{r}-\alpha
e^{-\frac{r}{M}}\left(\frac{r}{M}+2\right)-2=0\,.$ (12)

If we assume that the radius of the photon sphere is not influenced by the new
matter source then from eqn. (12), we obtain the following value for a primary
hair $l$ regardless of the value of $\alpha$.

$r_{ph}=3M\rightarrow l=\frac{5M}{e^{3}}\,.$ (13)

However, the visible size of the black hole shadow will change as,

$b=\frac{3\sqrt{3}M}{\sqrt{1-2\alpha e^{-3}}}\,.$ (14)

Now, as shown in eqn. (14), if there is no additional matter source
$\alpha\rightarrow 0$, the visible size of the black hole shadow remains
unchanged $b=3\sqrt{3}M$ \- as it would be. However, the presence of an extra
matter source increases the visible size of the black hole shadow if

$0<\alpha<\frac{e^{3}}{2}\,.$ (15)

One should note that the coupling constant is assumed to be less than unity,
so the case $\alpha=\frac{e^{3}}{2}\rightarrow b=\infty$ corresponds to
unphysical situation. If we consider the eqn. (12) in general case then the
impact parameter is

$\begin{split}b^{2}=&\frac{r_{ph}^{3}(r_{ph}+2M)}{r_{ph}^{2}-\beta
r_{ph}+M\beta},\\\ \\\\[-10.0pt] \beta=&\alpha l+2M.\end{split}$ (16)

One should note that both $l$ and $\alpha$ are supposed to be positive. So
$\beta>0$ and the black hole shadow exists for $r_{ph}>M$, as can be seen from
the denominator of eqn. (16). Also, if we assume that $\alpha=0$ then
$\beta=2M,r_{ph}=3M$ and $b=3\sqrt{3}M$ \- Schwarzschild case. However, the
influence of a primary hair $l$ on the radius of a photon sphere is not clear
from eqn. (12). We know that a parameter $\alpha$ is supposed to be small
$\alpha\ll 1$, so we can consider the radius of a photon sphere in the form,

$r_{ph}=r_{ph}^{(0)}+\alpha r_{ph}^{(1)}\,,$ (17)

where $r_{ph}^{(0)}=3M$ is the radius of a photon sphere for the Schwarzschild
black hole. If we apply the first condition of eqn. (6) to the effective
potential and write it in the form:

$V_{eff}=0\rightarrow b^{2}F(r)\left(1+\alpha G(r)\right)=1\,,$ (18)

where,

$\displaystyle F(r)$ $\displaystyle\equiv$ $\displaystyle\frac{r-2m}{r^{3}},$
$\displaystyle G(r)$ $\displaystyle\equiv$
$\displaystyle\frac{e^{-\frac{r}{M}}-\frac{l}{r}}{1-\frac{2M}{r}}.$ (19)

Now, applying the second condition of eqn.(6) into (18), one can obtains

$\frac{dV_{eff}}{dr}=0\rightarrow F^{\prime}(1+\alpha G)+\alpha
FG^{\prime}=0,$ (20)

where a prime denotes the derivative with respect to the $r$ coordinate. From
eqn. (20), one can easily find that $F^{\prime}(r_{ph}^{(0)})=0$. Substituting
eqn. (17) into (20) and expand on small parameter $\alpha$, we can finds

$\displaystyle F^{\prime\prime}(r_{ph}^{(0)})\alpha r_{ph}^{(1)}+\alpha
G^{\prime}(r_{ph}^{(0)})F(r_{ph}^{(0)})=0\rightarrow$ $\displaystyle
r_{ph}^{(1)}=-\frac{G^{\prime}(r_{ph}^{(0)})F(r_{ph}^{(0)})}{F^{\prime\prime}(r_{ph}^{(0)})}.$
(21)

Now, by using eqn. (II.1) into (II.1), we can evaluate the sign of
$r_{ph}^{(1)}$ in order to find out whether the primary hair will decreases or
increases the radius of a photon sphere. For this purpose, let’s find
$F(r_{ph}^{(0)})\,,F^{\prime\prime}(r_{ph}^{(0)})\,,G^{\prime}(r_{ph}^{(0)})$:

$\displaystyle F\left(r_{ph}^{(0)}\right)$ $\displaystyle=$
$\displaystyle\frac{1}{27M^{2}},$ $\displaystyle
F^{\prime\prime}\left(r_{ph}^{(0)}\right)$ $\displaystyle=$
$\displaystyle-\frac{2}{81M^{4}},$ $\displaystyle
G^{\prime}\left(r_{ph}^{(0)}\right)$ $\displaystyle=$
$\displaystyle\frac{l}{M^{2}}-\frac{5}{M}e^{-3}.$ (22)

Substituting these results into (II.1), one obtains for $r_{ph}^{(1)}$

$r_{ph}^{(1)}=\frac{3l}{2}-\frac{15M}{2}e^{-3}\,.$ (23)

From this result we can note the following points:

  * •

$l=5Me^{-3}\rightarrow r_{ph}^{(1)}=0$. We arrive at the condition (13), when
a primary hair doesn’t influence the radius of a photon sphere for the
Schwarzschild black hole;

  * •

$l>5Me^{-3}$. In this case, $r_{ph}>r_{ph}^{(0)}$, i.e. a primary hair $l$
increases the radius of a photon sphere in comparison with the Schwarzschild
black hole;

  * •

$l<5Me^{-3}$. This case leads to the smaller radius of photon sphere than in
Schwarzschild case $r_{ph}<r_{ph}^{(0)}$. However, in this case, the energy-
momentum tensor of extra matter fields violates the strong energy condition
bib:bh1 .

Thus, if one demands the fulfilment of the strong energy condition, then the
radius of the photon sphere always increases. Note that all these conditions
don’t depend upon the coupling constant $\alpha$. The only constraint we
impose on it is $\alpha\ll 1$.

###  II.2 model 2

The second model is given by the lapse function $f(r)$

$f(r)=1-\frac{2M}{r}-\frac{\alpha l}{r}+\frac{Q^{2}}{r^{2}}-\frac{\alpha
M}{r}e^{-\frac{r}{M}}\,.$ (24)

To fulfill the dominant energy condition, one should impose the following
restrictions on $Q$ and $l$ bib:geod

$Q^{2}\geq\frac{\alpha M^{2}}{e^{2}},~{}~{}l\geq\frac{M}{e^{2}}\,.$ (25)

If we substitute the lapse function (24) into the second condition of eqn.
(6), then one obtains

$\begin{split}V^{\prime}_{eff}(r_{ph})=0&\rightarrow\frac{6M+3\alpha
l}{r}-\frac{4Q^{2}}{r^{2}}-2\\\ &+\left(\frac{3\alpha
M}{r}+\alpha\right)e^{-\frac{r}{M}}\,.\end{split}$ (26)

In usual Reissner-Nordstrom solution ($\alpha=0$), one has

$\displaystyle 2r_{ph}$ $\displaystyle=$ $\displaystyle
3M+\sqrt{9M^{2}-8Q^{2}},$ $\displaystyle b$ $\displaystyle=$
$\displaystyle\frac{r^{2}_{ph}}{\sqrt{Mr_{ph}-Q^{2}}}.$ (27)

Note that in the hairy Schwarzschild case, one might have the same $r_{ph}$ as
in the Schwarzschild case. The comparison of eqn. (II.2) with eqn. (26) shows
that there is no photon sphere at $r=r_{ph}$ for the considered hairy model.
If we substitute (26) into the first condition of eqn. (6) then one obtains

$\displaystyle\beta=3+\frac{r_{ph}}{M}\,,$
$\displaystyle\kappa=\beta(r_{ph}^{2}-2mr_{ph}+Q^{2})+$
$\displaystyle+6Mr_{ph}+3\alpha lr_{ph}-4Q^{2}-2r_{ph}^{2}\,,$ $\displaystyle
b=\frac{\sqrt{\beta}\,r_{ph}^{2}}{\sqrt{\kappa}}\,.$ (28)

Again, like in the previous model, the expression (26) is too cumbersome to
find out what an impact a primary hair $l$ has on the radius of a photon
sphere. Thus, To evaluate this influence, we utilize the same method as in the
previous subsection. We assume that the radius of a photon sphere has the form
(17) and, using the first condition of eqn. (6), we write down the effective
potential in the form of eqn. (18). The functions $F(r)$ and $G(r)$ for the
spacetime with lapse function (24) are given by,

$\displaystyle r_{ph}^{(0)}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left(3M+\sqrt{9M^{2}-8Q^{2}}\right),$ $\displaystyle
F(r)$ $\displaystyle=$ $\displaystyle\frac{r^{2}-2Mr+Q^{2}}{r^{4}},$
$\displaystyle G(r)$ $\displaystyle=$
$\displaystyle-\frac{r\left(l+Me^{-\frac{r}{M}}\right)}{r^{2}-2Mr+Q^{2}}.$
(29)

Now applying the same method as described above, we evaluate functions $F$,
$F^{\prime\prime}$ and $G$ at the radius $r_{ph}^{(0)}$, and then we have

$\displaystyle F\left(r_{ph}^{(0)}\right)$ $\displaystyle=$
$\displaystyle\frac{Mr_{ph}^{(0)}-Q^{2}}{\left(r_{ph}^{(0)}\right)^{4}},$
$\displaystyle F^{\prime\prime}\left(r_{ph}^{(0)}\right)$ $\displaystyle=$
$\displaystyle\frac{8Q^{2}-6Mr_{ph}^{(0)}}{\left(r_{ph}^{(0)}\right)^{6}},$
$\displaystyle G^{\prime}\left(r_{ph}^{(0)}\right)$ $\displaystyle=$
$\displaystyle\frac{3l+3Me^{-\frac{r_{ph}^{(0)}}{M}}+r_{ph}^{(0)}e^{-\frac{r_{ph}^{(0)}}{M}}}{Mr_{ph}^{(0)}-Q^{2}},$
(30)

substituting these results into the definition $r_{ph}^{(1)}$ (II.1), we can
obtain

$r_{ph}^{(1)}=-\left(r_{ph}^{(0)}\right)^{2}\frac{3l+3Me^{-\frac{r_{ph}^{(0)}}{M}}+r_{ph}^{(0)}e^{-\frac{r_{ph}^{(0)}}{M}}}{8Q^{2}-6Mr_{ph}^{(0)}}\,.$
(31)

As one can see that the sign of $r_{ph}^{(1)}$ depends on the sign of the
expression

$8Q^{2}-6Mr_{ph}^{(0)}=8Q^{2}-9M^{2}-3M\sqrt{9M^{2}-8Q^{2}}\,.$ (32)

Introducing a dimensionless parameter $\xi$ by

$\xi^{2}=\frac{Q^{2}}{M^{2}},~{}~{}0\leq\xi\leq 1\,,$ (33)

one can write eqn. (32) as,

$-9M^{2}\left(1-\frac{8}{9}\xi^{2}+\sqrt{1-\frac{8}{9}\xi^{2}}\right)<0\,.$
(34)

Thus, we can come to the conclusion that in the hairy Reissner-Nordstrom case,
the primary hair $l$ always increases the radius of a photon sphere in
comparison with the RN black hole. The radius of a photon sphere can be equal
to $3M$ like in Schwarzschild spacetime. In this case, from eqn. (26), we
obtain

$r=3M\rightarrow l=\frac{4Q^{2}}{9\alpha M}-2e^{-3}\,.$ (35)

##  III Shadow of overcharged black hole

In pure Reissner-Nordstrom spacetime, if a black hole is overcharged, i.e.
$Q^{2}>M^{2}$, then this spacetime contains a naked singularity. It doesn’t
cast a shadow. However, the primary hair $l$ affects the spacetime and
overcharged black holes with $Q^{2}>M^{2}$ can exist. We will prove that
overcharged black holes can also cast a shadow. We will show it implicitly
since the lapse function (24) and the second condition of eqn. (26) are too
cumbersome to obtain an analytical solution. First, we will show that there
can be a black hole with $Q^{2}>M^{2}$, and then we will prove that one of the
radii of the photon spheres is outside the event horizon and the other is
inside. Here, we consider the charge $Q$ as

$Q^{2}=M^{2}+\delta,~{}~{}\delta\ll M^{2},$ (36)

then the lapse function (24) can be written as,

$f(r)=\left(1-\frac{M}{r}\right)^{2}+\frac{\delta}{r^{2}}-\frac{\alpha
l}{r}-\frac{\alpha M}{r}e^{-\frac{r}{M}}\,.$ (37)

From eqn. (37), it is obvious that we can choose parameters
$\delta\,,\alpha\,,l$ such that $f(M)<0$ and $f(2M)>0$, which indicates the
existence of the event horizon in the region $M<r<2M$. Now, by substituting
the charge (36) into the second condition of eqn. (26), one obtains:

$\begin{split}G(r)&\equiv\frac{6M+3\alpha
L}{r}-\frac{4Q^{2}+4\delta}{r^{2}}-\\\ &2+\left(\frac{3\alpha
M}{r}+\alpha\right)e^{-\frac{r}{M}}=0\,.\end{split}$ (38)

Neglecting $\frac{\delta}{M^{2}}$ in comparison with $(\frac{3\alpha
l}{2M}+\left(\frac{3\alpha l}{2M}+\alpha\right)e^{-2})$111To obtain this
expression, we put $r=2M$ into the expression of eqn. (38)., it is easy to
show that $G(M)$ and $G(2M)$ are both positive indicating that the first
photon sphere radius is the inside the event horizon. However, if one
satisfies the condition

$\frac{4}{9}\geq\frac{\alpha l}{M}+\left(\frac{\alpha
l}{M}+\alpha\right)e^{-3}\,,$ (39)

then the radius of the second photon sphere is located outside the event
horizon in the region $2M<r<3M$ and can form a black hole shadow with an
impact parameter of the form

$b=\frac{r_{ph}}{\sqrt{f\left(r_{ph}\right)}}.$ (40)

Thus, we implicitly showed that overcharged hairy Reissner-Nordstrom black
hole can cast a shadow.

##  IV Specific intensity of a shadow

In the previous section, we have considered an idealized scenario where the
light source is located far behind a black hole, and there are no light
sources between the observer and the compact object. A real astrophysical
black hole is typically surrounded by an accretion disc. We will compute the
influence of plasma on the radius of a photon sphere in the next section. In
this section, we consider an optically thin, radiating, accreting flow that
surrounds the hairy black hole and intensity distributions as a function of
the impact parameter. The radiation can be emitted in the vicinity of a black
hole, including the region inside the photon sphere. This radiation will
affect the specific intensity $I_{\nu_{0}}$ at the observed photon frequency
$\nu_{0}$ in the observer’s sky at the point $(x,y)$. It has the following
form bib:jaros ; bib:bambi

$I_{\nu_{0}}\left(x,y\right)=\int_{\gamma}G^{\prime
3}J\left(\nu_{e}\right)dl_{prop},$ (41)

here, $\nu_{e}$ is emitted frequency, $G^{\prime}=\frac{\nu_{0}}{\nu_{e}}$ is
the redshift factor, $J\left(\nu_{e}\right)$ is the emitter’s rest frame
emissivity per unit volume, $dl_{prop}=-K_{i}u^{i}_{(e)}d\lambda$ is the
infinitesimal proper length in the rest frame of the emitter,
$K^{i}\equiv\frac{dx^{i}}{d\lambda}$ is the photon four-momentum,
($K^{i}K_{i}=0$), $u^{i}_{(e)}$ is the four-velocity of an emitter and
$\lambda$ is an affine parameter. The integration is evaluated along the
observed photon path $\gamma$. The redshift $G^{\prime}$ in terms of the four-
velocities $K^{i},u^{i}_{(0)},u^{i}_{(e)}$ is given by bib:hartle

$G^{\prime}=\frac{K_{i}u^{i}_{(0)}}{K_{j}u^{j}_{(e)}}\,,$ (42)

where $u^{i}=\left\\{1\,,0\,,0\,,0\right\\}$ is the four-velocity of an
observer. We consider accreting gas which is in radial free fall with four-
velocity $u^{i}_{(e)}$ of the form bib:bambi

$\displaystyle u^{0}_{(e)}$ $\displaystyle=$ $\displaystyle\frac{1}{f},$
$\displaystyle u^{1}_{(e)}$ $\displaystyle=$ $\displaystyle-\sqrt{1-f},$
$\displaystyle u^{2}_{(e)}$ $\displaystyle=$ $\displaystyle u^{3}_{(e)}=0.$
(43)

The four-vector $K^{i}\equiv\frac{dx^{i}}{d\lambda}$ is given by:

$\displaystyle K^{0}$ $\displaystyle=$ $\displaystyle\frac{E}{f},$
$\displaystyle K^{1}$ $\displaystyle=$
$\displaystyle\pm\sqrt{E^{2}-f\frac{L^{2}}{r^{2}}},$ $\displaystyle K^{2}$
$\displaystyle=$ $\displaystyle 0,$ $\displaystyle K^{3}$ $\displaystyle=$
$\displaystyle\frac{L^{2}}{r^{2}}.$ (44)

Here the sign $K^{1}$ depends on the photon motion, i.e. $+$ if it goes away
from a black hole and $-$ if approaches it. Substituting (IV) and (IV) into
(42), we find the redshift $G^{\prime}$ in the form

$G^{\prime}=\frac{K_{0}f}{K_{0}-fK_{1}\sqrt{1-f}}\,.$ (45)

For the specific emissivity $J$, we assume the following simple model
bib:bambi  in which the emission is monochromatic with emitter’s rest-frame
frequency $\nu^{*}$ and the emission has a $\frac{1}{r^{2}}$ radial profile:

$J(\nu_{e})\propto\frac{\delta(\nu_{e}-\nu^{*})}{r^{2}}\,.$ (46)

The proper length $dl_{prop}$ is given by

$dl_{prop}=-\frac{-K_{0}}{G^{\prime}K^{1}}dr\,.$ (47)

Now substituting (42), (46) and (47) into (41) and integrating over all
observed frequencies, we find the observed intensity

$I_{o}(x,y)\sim-\int_{\gamma}\frac{G^{\prime 3}K_{0}}{r^{2}K^{1}}dr\,.$ (48)

As we have shown, a primary hair increases the radius of a photon sphere. As a
result, the specific intensity decreases, as shown in figures 1 (c) and (d).

Figure 1: These plots show observed intensity distributions as a function of
the impact parameter $b$ for $M=1$. Plot (a) shows the function $I_{0}$ versus
the impact parameter $b$ for the Schwarzschild black hole. Plot (b) shows the
function $I_{0}$ versus the impact parameter $b$ for the Reissner-Nordstrom
black hole with $Q=0.98$. Plot (c) shows the function $I_{0}$ versus the
impact parameter $b$ for model $1$ if $\alpha=0.5$ and $l=2$. Plot (d) shows
the function $I_{0}$ versus the impact parameter $b$ for model $2$ if
$Q=0.98$, $\alpha=0.5$ and $l=2$.

##  V The influence of a plasma on a radius of a photon sphere

In an idealized situation, we assume that there is no matter between the
observer and the black hole. However, as we showed above, the accretion disc
has an impact on the specific intensity. Here we consider the influence of
plasma for a radius of a photon sphere of a hairy black hole. Our analysis
goes close to the method described in bib:tsupko_plasma . We consider low-
density, non-magnetized pressureless plasma. In the absence of the pressure,
the particles move along timelike geodesics for which $g_{ik}u^{i}u^{k}=-1$.
In order to estimate the influence of plasma, one needs to calculate the
plasma frequency $\omega_{p}$ which, for our model, depends on the plasma
density $\rho$. Let’s consider the continuity equation

$\left(\sqrt{-g}\rho u^{i}_{(p)}\right)_{,i}=0\,.$ (49)

Here $g=\det g_{ik}$ and $u^{i}_{(p)}$ is the four-velocity of the particles.
We simplify our model by assuming that electron and proton four-velocities are
coincides. The plasma is a neutral hydrogen one. We neglect the electron mass
$m_{e}$ in comparison with proton mass $m_{p}$. So, the plasma density $\rho$
can be written as

$\rho\sim m_{p}N\,.$ (50)

Where $N$ is the number density of protons. Again, we assume the equatorial
geodesic motion i.e. $\theta=\frac{\pi}{2}$. As a result the metric tensor
components depend only upon the radial coordinate $r$ and the eqn. (49)
becomes

$\left(r^{2}\rho u^{1}\right)_{,r}=0\,.$ (51)

An integration gives,

$4\pi\rho r^{2}u^{1}=-\dot{m}=const.\,,$ (52)

where $\dot{m}$ is a stationary mass flux. Here, we assume that particles
starts its radial motion from rest at infinity. It means $E=1$, then $u^{1}$
reads

$u^{1}=-\sqrt{1-f}\,.$ (53)

The density $\rho$ can be expressed through the ratio of plasma $\omega_{p}$
and photon $\omega_{0}$ frequencies

$\frac{\omega_{p}^{2}}{\omega_{0}}=\frac{4\pi
e^{2}N}{m_{e}\omega_{0}^{2}}=\frac{4\pi
e^{2}\rho}{m_{p}m_{e}\omega_{0}^{2}}=\beta_{0}\frac{1}{r^{2}\sqrt{1-f}}\,,$
(54)

where,

$\beta_{0}=\frac{e^{2}\dot{m}}{m_{e}m_{p}\omega_{0}^{2}}.$ (55)

As we have mentioned above, we consider a low density plasma i.e.
$\omega_{p}\ll\omega_{0}$. In this case, the angular size $(\theta_{sh})$ of
the shadow observed at the radius $r_{o}$ is given by (See (56) in
bib:tsupko_plasma )

$\begin{split}\sin^{2}\theta_{sh}&\approx\frac{r_{ph}^{2}}{\left(r_{o}^{(0)}\right)^{2}}\frac{f\left(r_{o}\right)}{f\left(r_{ph}^{(0)}\right)}\biggl{(}1-\beta_{0}f\left(r_{ph}^{(0)}\right)\times\\\
&\sqrt{1-f\left(r_{ph}^{(0)}\right)}+\beta_{0}f(r_{o})\sqrt{1-f(r_{o})}\biggr{)}\,.\end{split}$
(56)

Where $r_{ph}^{(0)}$ is the radius of a photon sphere without plasma. For
estimating the plasma correction numerically, we re-express all quantities in
Gaussian CGS units, i.e., we restore factors of $c$ and $G$ . Then $\beta_{0}$
reads

$\beta_{0}=\frac{e^{2}\dot{m}c^{3}}{\sqrt{2GM}(GM)^{3/2}m_{e}m_{p}\omega_{0}^{2}}\,.$
(57)

The accretion rate $\dot{m}$ can be estimated through observed luminosity as,

$\lambda\sim\eta\dot{m}c^{2}\,,$ (58)

where $\eta$ is dimensionless quantity related to accretion efficiency. An
observed photon frequency at infinity is related to the wavelength
$\lambda_{0}$ by relation,

$\omega_{0}=\frac{2\pi c}{\lambda_{0}}\,.$ (59)

We consider two models $\eta=10^{-4}$ and $\eta=0.1$. and two massive black
holes Sgr $A^{*}$ and M87*.

$\begin{split}\text{Sgr $A^{*}$}\quad,\quad M=4.3\times
10^{6}\,M_{\odot}\,,\\\ r_{o}=8.3\,Kpc\,,\\\ \text{M87*}\quad,\quad M=3\times
10^{9}M_{\odot}\,,\\\ r_{o}=18\,Mpc\,,\\\
\lambda=10^{6}\lambda_{\odot}\,,\lambda_{0}=1\,\text{mm}\,,\lambda_{0}=10\,\text{cm}\,,\end{split}$
(60)

without plasma, the angular radius of the shadow of black hole Sgr $A^{*}$ has
the value $\theta_{sh}\simeq 26.572$ $\mu\text{as}$. In general, as one can
see from tables 1 and 2, the influence of a primary hair on the visible size
of a shadow is much stronger than the influence of plasma except for several
particular cases of large $\lambda_{0}$ and small $\eta$. For pure M87*, the
angular radius of the shadow has the following value $\theta_{sh}\simeq 8.548$
$\mu\text{as}$, and the effect of plasma gives a change in the size of the
shadow of the order $10^{-6}\mu\text{as}$ for $\eta=10^{-4}$ and
$\lambda_{0}=10$, whereas primary hair for $\alpha=0.5$ and $l=1.5$ gives the
value $\theta_{sh}\simeq 11.602\mu\text{as}$.

Table 1:  This table shows the value of the angular radius of the shadow for Sgr $A^{*}$ black hole for different quantities of parameters $\alpha$ and $l$. The Upper part corresponds to the Schwarzschild case, the middle part to the Reissner-Nordstrom black hole with $Q=0.1$ and the bottom of the table for $Q=0.9$. The units of measure are $\mu\text{as}$ Schwarzschild | $l=0.5$ | $l=1$ | $l=1.5$
---|---|---|---
$\alpha=0.01$ | $26.618$ | $26.685$ | $26.751$
$\alpha=0.1$ | $27.044$ | $27.718$ | $28.391$
$\alpha=0.5$ | $29.052$ | $32.583$ | $36.066$
RN $Q=0.1$ | $l=0.5$ | $l=1$ | $l=1.5$
$\alpha=0.01$ | $26.601$ | $26.667$ | $26.734$
$\alpha=0.1$ | $27.254$ | $27.915$ | $28.576$
$\alpha=0.5$ | $30.074$ | $33.332$ | $36.61$
RN $Q=0.9$ | $l=0.5$ | $l=1$ | $l=1.5$
$\alpha=0.01$ | $22.190$ | $22.274$ | $22.359$
$\alpha=0.1$ | $23.078$ | $23.887$ | $24.684$
$\alpha=0.5$ | $26.594$ | $30.248$ | $33.839$
Table 2:  This table shows the influence of plasma on the angular radius of the shadow for Sgr $A^{*}$ black hole for different parameters $\eta$ and $\lambda_{0}$. The Upper part corresponds to the pure Schwarzschild case, the middle part to the pure Reissner-Nordstrom black hole with $Q=0.1$, and the bottom of table for $Q=0.9$. The units of measure are $\mu\text{as}$ Schwarzschild | $\eta=10^{-4}$ | $\eta=0.1$
---|---|---
$\lambda_{0}=10$ | $22.4$ | $26.568$
$\lambda_{0}=0.1$ | $26.571$ | $26.572$
RN $Q=0.1$ | $\eta=10^{-4}$ | $\eta=0.1$
$\lambda_{0}=10$ | $22.367$ | $26.524$
$\lambda_{0}=0.1$ | $26.527$ | $26.527$
RN $Q=0.9$ | $\eta=10^{-4}$ | $\eta=0.1$
$\lambda_{0}=10$ | $19.077$ | $22.085$
$\lambda_{0}=0.1$ | $22.087$ | $22.087$

##  VI Discussions and Conclusions

In this paper, we have considered the influence of a primary hair $l$ on the
shadow of hairy Schwarzschild and Reissner-Nordstrom black holes obtained by
gravitational decoupling. We have found out that in the hairy Schwarzschild
case, a primary hair $l$ either doesn’t have any impact or always increases
the radius of the photon sphere. It can decrease the radius of a photon
sphere, but in this case, the energy-momentum tensor $\Theta_{ik}$ violates
the strong energy condition. For the Reissner-Nordstrom case, we have shown
that an additional matter field $\Theta_{ik}$ always increases the radius of a
photon sphere, and this radius can even be equal to $3M$ like in the pure
Schwarzschild case. We also proved that a black hole can exist even in an
overcharged case, i.e. when $Q^{2}>M^{2}$ and it also casts a shadow. We have
calculated the specific intensity and showed that a primary hair $l$ always
decreases it in both Schwarzschild and Reissner-Nordstrom cases. The influence
of a pressureless plasma on the black hole has been also considered. We have
found out that, in general, a primary hair $l$ has a stronger influence on the
visible size in comparison with pressureless plasma. However, the real
astrophysical black holes are supposed to be rotating ones. In bib:geod , it
has been shown that under some conditions on $\alpha$ and $l$, the hairy
Schwarzshild solution can be a mimicker of a Kerr black hole. In the case of a
shadow, it can be viewed as a Kerr black hole mimicker only if we observe a
shadow on a certain angle. Because, in general the shadow in Kerr spacetime is
flattened, but for the black holes considered in this paper, a shadow is
always circular in shape. The results obtained in this paper, can be further
used for the construction of a black hole shadow for rotating hairy Kerr and
Kerr-Newman spacetimes obtained by gravitational decoupling bib:hairy_kerr .
After this, all these models can be tested by comparison with the images
obtained by EHT collaboration.

Acknowledgments: V. Vertogradov thanks the Basis Foundation (grant number
23-1-3-33-1) for the financial support.

## References

  * (1) J. L. Synge, Mon. Not. Roy. Astron. Soc. 131, no.3, 463-466 (1966).
  * (2) K. Akiyama et al. [Event Horizon Telescope], Astrophys. J. Lett. 875, L1 (2019).
  * (3) K. Akiyama et al. [Event Horizon Telescope], Astrophys. J. Lett. 875, no.1, L2 (2019).
  * (4) K. Akiyama et al. [Event Horizon Telescope], Astrophys. J. Lett. 930, no.2, L12 (2022).
  * (5) S. Vagnozzi, R. Roy, Y. D. Tsai, L. Visinelli, M. Afrin, A. Allahyari, P. Bambhaniya, D. Dey, S. G. Ghosh and P. S. Joshi, et al. Class. Quant. Grav. 40, no.16, 165007 (2023).
  * (6) O. Y. Tsupko, Z. Fan and G. S. Bisnovatyi-Kogan, Class. Quant. Grav. 37, no.6, 065016 (2020).
  * (7) S. Vagnozzi, C. Bambi and L. Visinelli, Class. Quant. Grav. 37, no.8, 087001 (2020).
  * (8) V. Perlick and O. Y. Tsupko, Phys. Rept. 947, 1-39 (2022).
  * (9) C. M. Claudel, K. S. Virbhadra and G. F. R. Ellis, J. Math. Phys. 42, 818-838 (2001).
  * (10) S. Hod, Phys. Lett. B 718, 1552-1556 (2013).
  * (11) F. S. Khoo and Y. C. Ong, Class. Quant. Grav. 33, no.23, 235002 (2016).
  * (12) Y. Decanini, A. Folacci and B. Raffaelli, Phys. Rev. D 81, 104039 (2010).
  * (13) A. A. Shoom, Phys. Rev. D 96, no.8, 084056 (2017).
  * (14) C. Cederbaum and G. J. Galloway, Class. Quant. Grav. 33, 075006 (2016).
  * (15) T. Johannsen, Astrophys. J. 777, 170 (2013).
  * (16) A. A. Abdujabbarov, L. Rezzolla and B. J. Ahmedov, Mon. Not. Roy. Astron. Soc. 454, no.3, 2423-2435 (2015).
  * (17) Z. Younsi, A. Zhidenko, L. Rezzolla, R. Konoplya and Y. Mizuno, Phys. Rev. D 94, no.8, 084025 (2016).
  * (18) V. Perlick, O. Y. Tsupko and G. S. Bisnovatyi-Kogan, Phys. Rev. D 97, no.10, 104062 (2018).
  * (19) O. Y. Tsupko, Phys. Rev. D 103, no.10, 104019 (2021).
  * (20) P. Bambhaniya, D. Dey, A. B. Joshi, P. S. Joshi, D. N. Solanki and A. Mehta, Phys. Rev. D 103, no.8, 084005 (2021).
  * (21) P. Bambhaniya, S. K, K. Jusufi and P. S. Joshi, Phys. Rev. D 105, no.2, 023021 (2022).
  * (22) J. Solanki, V. Perlick  Phys. Rev. D 105 ,064056, (2022).
  * (23) Y. Koga, N. Asaka, M. Kimura and K. Okabayashi,  Phys. Rev. D 105, 104040 (2022).
  * (24) Y. Heydarzade and V. Vertogradov, arXiv:2311.08930 [gr-qc].
  * (25) O. Yu. Tsupko, G. S. Bisnovatyi-Kogan,  Int. J. Mod. Phys. D 29, 2050062 (2020).
  * (26) V. Perlick, O. Yu. Tsupko, G. S. Bisnovatyi-Kogan. Phys. Rev. D 92, 104031 (2015).
  * (27) G. S. Bisnovatyi-Kogan, O. Yu. Tsupko, Plasma Phys. Rep. 41, 562 (2015).
  * (28) A. K Mishra, S. Chakraborty and S. Sarkar.  Phys. Rev. D 99, 104080 (2019).
  * (29) R. Shaikh, P. Kocherlakota, R. Narayan and P. S. Joshi. Mon. Not. R. Astron. Soc. 482, 52 (2019).
  * (30) S. Saurabh, P. Bambhaniya and P. S. Joshi, Astron. Astrophys. 682, A113 (2024).
  * (31) S. Saurabh, P. Bambhaniya and P. S. Joshi, arXiv:2202.00588 [gr-qc].
  * (32) D. N. Solanki, P. Bambhaniya, D. Dey, P. S. Joshi and K. N. Pathak, Eur. Phys. J. C 82, no.1, 77 (2022).
  * (33) A. B. Joshi, D. Dey, P. S. Joshi and P. Bambhaniya, Phys. Rev. D 102, no.2, 024022 (2020).
  * (34) S. W. Hawking, M. J. Perry, A. Strominger Phys. Rev. Lett. 116, 231301 (2016).
  * (35) J. Ovalle, Int. J. Mod. Phys. Conf. Ser., 41, 1660132 (2016).
  * (36) R. Casadio, J. Ovalle, R. da Rocha  Class. Quantum Grav. 32, 215020 (2015).
  * (37) J. Ovalle, R. Casadio, R. da Roch, A. Sotomayor and Z. Stuchlik Eur. Phys. J. C 78, 960. (2018).
  * (38) J. Ovalle Phys. Rev. D 95, 104019 (2017).
  * (39) J. Ovalle Phys. Lett. B 788, 213. (2019).
  * (40) E. Contreras, J. Ovalle, R. Casadio Phys. Rev. D 103, 044020. (2021).
  * (41) J. Ovalle, R. Casadio, E. Contreras and A. Sotomayor Phys. Dark Universe 31, 100744 (2021).
  * (42) S. Mahapatra, I. Banerjee Phys. Dark Universe 39, 101172 (2023).
  * (43) V. Vertogradov, M. Misyura Universe 8 567, (2022).
  * (44) V. Vertogradov, M. Misyura Phys. Sci. Forum 7, 27 (2023).
  * (45) J. Ovalle, R. Casadio, A. Giusti Phys. Lett. B 844, 138085 (2023).
  * (46) R. T. Cavalcanti, K. d. S. Alves,J. M. .H. d. Silva  Universe 8, 363 (2022).
  * (47) V. D. Vertogradov, D. A. Kudryavcev  Phys. Complex Syst. 4,2 (2023).
  * (48) A. Ramos, C. Arias, R. Avalos, E. Contreras Annals Phys 431 431, 168557 (2021).
  * (49) S. K. Jha, A. Rahaman arXiv:2205.06052 [gr-qc].
  * (50) Z. Li, F. Yuan Phys. Rev. D 108, 024039 (2023).
  * (51) M. Jaroszynski, A. Kurpiewski  A.A 326, 419 (1997).
  * (52) C. Bambi  Phys. Rev. D 87, 107501 (2013).
  * (53) J. B. Hartle, _Gravity: An introduction to Einstein’s general relativity_ (Pearson, San Francisco, 2003).


# LoMAE: Low-level Vision Masked Autoencoders for Low-dose CT Denoising

Dayang Wang  Student Member, IEEE  Yongshun Xu  Shuo Han  Zhan Wu
Li Zhou  Bahareh Morovati  Hengyong Yu  Senoir Member, IEEE The paper was
submitted on Oct. 1, 2023. This work was supported in part by NIH/NIBIB under
grants R01EB032807 and R01EB3473, and NIH/NCI under grant R21CA264772.D. Wang,
Y. Xu, S. Han, L. Zhou, B. Morovati, and H. Yu are with the Department of
Electrical and Computer Engineering, University of Massachusetts Lowell,
Lowell, MA, USA, 01854.Z. Wu is with the Laboratory of Image Science and
Technology, Southeast University, Nanjing 210096, China, and also with the Key
Laboratory of Computer Network and Information Integration, Southeast
University, Ministry of Education, Nanjing 210096, China.Dr. H. Yu is the
corresponding author (e-mail: hengyong-yu@ieee.org).

###### Abstract

Low-dose computed tomography (LDCT) offers reduced X-ray radiation exposure
but at the cost of compromised image quality, characterized by increased noise
and artifacts. Recently, transformer models emerged as a promising avenue to
enhance LDCT image quality. However, the success of such models relies on a
large amount of paired noisy and clean images, which are often scarce in
clinical settings. In the fields of computer vision and natural language
processing, masked autoencoders (MAE) have been recognized as an effective
label-free self-pretraining method for transformers, due to their exceptional
feature representation ability. However, the original pretraining and fine-
tuning design fails to work in low-level vision tasks like denoising. In
response to this challenge, we redesign the classical encoder-decoder learning
model and facilitate a simple yet effective low-level vision MAE, referred to
as LoMAE, tailored to address the LDCT denoising problem. Moreover, we
introduce an MAE-GradCAM method to shed light on the latent learning
mechanisms of the MAE/LoMAE. Additionally, we explore the LoMAE’s robustness
and generability across a variety of noise levels. Experiments results show
that the proposed LoMAE can enhance the transformer’s denoising performance
and greatly relieve the dependence on the ground truth clean data. It also
demonstrates remarkable robustness and generalizability over a spectrum of
noise levels.

Low-dose CT, Masked Autoencoder, Self-pretraining, Transformer.

##  1 Introduction

In recent years, Low-dose computed tomography (LDCT) has become the mainstream
in clinical applications due to the potential risks from x-ray radiation in
normal-dose CT (NDCT). Nonetheless, LDCT images exhibit compromised image
quality and reduced diagnostic value, presenting a significant barrier for
their widespread application. To overcome this issue, numerous deep learning
models have been explored along this direction [1, 2, 3, 4].

Figure 1: Pipelines of the original MAE (a) and LoMAE (b) in LDCT denoising.

Recently, masked autoencoders (MAE) have emerged as an excellent self-
supervised learning strategy for various computer vision tasks [5, 6]. He et.
al. revealed that small fraction of an image can infer complex and holistic
visual concepts like semantics [5]. Feichtenhofer et. al. applied MAE on
videos and found that the MAE method can almost eliminate inductive bias on
spacetime with strong representations [7]. Germain et. al. modified MAE to
data distribution estimation and interpreted the MAE output as conditional
probabilities and the full joint probability, thus prototyping a strong
distribution estimator with faster speed and higher scalability [8]. Gao et.
al. showed that more discriminative representations can be further learned via
a MAE multi-scale hybrid convolution-transformer architecture [9].
Furthermore, they proposed a block-wise masking strategy to prevent
information leakage and guarantee computational efficiency. Besides these
explorations on deep architecture design and general natural images, the MAE
was also explored in the medical imaging field. Zhou et al. were among the
first to demonstrate that MAE can improve the performance of medical image
segmentation and classification [10]. However, these works are on high-level
vision tasks. While masked image modeling has been explored as a robust
denoiser [11, 12, 13], currently, there lacks sufficient exploration of MAE as
self-pretraining on low-level vision tasks like LDCT denoising.

Hence, we are inspired to delve into the potential of the MAE in LDCT
denoising. We believe this study is significant for two compelling reasons: i)
clinical applications often face the challenge of limited ground truth data
availability. The self-pretraining paradigm of the MAE can reap the benefits
of the unlabeled data, making it an ideal choice for LDCT denoising task. ii)
LDCT denoising places significant importance on structural preservation and
enhancement. But, the anatomical structures within a CT image are supposed to
be intricately connected both mechanically and functionally. Here, the MAE’s
innate ability to aggregate contextual information and infer masked structures
has great potential to fortify the interdependence between anatomical regions,
thus complementing structural loss. Consequently, not only does this
reinforcement bolster the model’s robustness, but it also potentially augments
its capacity to generalize across a diverse spectrum of noise types and
levels, rendering it invaluable in LDCT denoising endeavors.

The contributions of this paper are fourfolds. 1) Tailored low-level vision
MAE: a low-level vision masked autoencoder (LoMAE) framework is designed and
specifically tailored to the demands of LDCT denoising tasks. This novel
approach ensures that our method is well-aligned with the unique challenges
presented by LDCT image enhancement. 2) Reduced ground truth dependency: our
LoMAE represents a breakthrough in reducing the reliance on ground truth data.
Remarkably, even with access to just two patients’ clean data, LoMAE
outperforms fully supervised baseline models. This achievement underscores the
practicality and effectiveness of LoMAE in scenarios where annotated data is
scarce. 3) Robustness and generalizability: the exceptional robustness and
generalizability of the LoMAE model are demonstrated across a spectrum of
noise levels. This finding underscores the model’s resilience and its ability
to deliver consistent performance in different real-world settings. 4) MAE-
GradCAM insight: for the first time, the MAE-GradCAM method is introduced,
shedding light on the intricate pretraining process of the Masked Autoencoder.
This novel approach unveils hidden patterns and provides valuable insights
into the inner workings of MAE-based models.

This manuscript is a significant extension of the IEEE International Symposium
on Biomedical Imaging (ISBI) 2023 paper entitled ”Masked Autoencoders for Low-
dose CT Denoising” [14, 15], which is the first work to adapt MAE to the low-
level vision tasks like denoising. In this paper, more extensive experiments
are conducted to demonstrate the properties of the proposed LoMAE strategy on
the LDCT denoising problem.

##  2 Related Works

There are two lines of research which are highly related to this manuscript on
LDCT denoising.

Convolutional Models. In the past few years, deep learning has been prevailing
in various science and engineering fields [16, 17, 18]. Its applications in
LDCT have also achieved many state-of-the-art performances. Chen et al. are
the pioneers in designing a residual encoder-decoder convolutional model (RED-
CNN) for the LDCT denoising [19]. Shan proposed a modularized deep neural
network (MAPNN) for LDCT reconstruction and achieved better qualities than
commercial iterative reconstruction methods with higher structural fidelity
and faster speed [3]. Fan et al. introduced a quadratic encoder-decoder LDCT
denoising network with more efficiency and robustness. Despite these
achievements, the CNN models are limited in capturing global contextual
information, thus less powerful to model long-stride structure dependence
among different anatomical regions in the CT images.

Transformer Models. In recent years, transformer models have gained
significant attention for their proficiency in capturing global contextual
information. Several researches also extended their applications to the realm
of LDCT denoising. Wang et al. advanced the field by introducing a
convolution-free encoder-decoder transformer, , providing insights into the
static and dynamic latent learning behaviors of the model through attention
map analysis [4]. Meanwhile, Yang et al. pioneered the development of a
sinogram inner-structure transformer, leveraging sinogram inner-structures to
suppress LDCT noise effectively [20]. Recently, the Swin transformer has
gained prominence as a versatile backbone architecture for a variety of
downstream tasks [21]. SwinIR is its notable variant, specifically tailored
for tasks such as image denoising, super-resolution, and artifact reduction
[22]. Notably, some SwinIR-based transformer models or modules have also been
introduced to CT image enhancement/denoising and demonstrating exceptional
performances [23]. Fan et al. also proposed another U-shaped Denoising
transformer, namely SUNet [24], which was based on Swin transformer modules
and achieved the SOTA performance in image denoising [24].

##  3 Methods

###  3.1 Noise Model

In LDCT images, two primary types of noise are prevalent. The first type is
thermal noise, which emanates from the detection system. This noise adheres to
a Gaussian distribution and is structurally independent. The second type is
photon statistical noise due to the inherent quantum fluctuations during the
X-ray emission. This noise typically conforms to a Poisson distribution that
is associated with the object structures. The LDCT noise is a complex
interplay of both Gaussian and Poisson components, as expressed by the
equation:

$\boldsymbol{\psi}\thicksim\mathcal{N}(0,\sigma^{2})+\mathcal{P}(\xi(\boldsymbol{\rho})),$
(1)

where $\boldsymbol{\psi}$ is the noise variable, $\mathcal{N}(0,\sigma^{2})$
denotes the Gaussian distribution whose variance is $\sigma^{2}$, and
$\mathcal{P}(\xi(\boldsymbol{\rho}))$ is Poisson distribution in which
$\xi(\boldsymbol{\rho})$ represents the latent mapping from projections
$\boldsymbol{\rho}$ of the clean image structure to the distribution
expectancy.

The proposed LDCT denoising model directly learns from paired noisy LDCT
images to their corresponding clean NDCT images $\boldsymbol{y}$ directly.
Experimentally, $\mathbf{L}_{1}$ and $\mathbf{SSIM}$ losses are employed to
suppress these noises during training. Mathematically, the $\mathbf{L}_{1}$
and $\mathbf{SSIM}$ losses are calculated as follows:

$\mathbf{L}_{1}(\tilde{\boldsymbol{y}},\boldsymbol{y})=||\tilde{\boldsymbol{y}}-\boldsymbol{y}||_{1},$
(2)
$\mathrm{SSIM}({\tilde{\boldsymbol{y}}},{\boldsymbol{y}})=\frac{(2u_{\tilde{\boldsymbol{y}}}u_{\boldsymbol{y}}+c_{1})(2\sigma_{{\tilde{\boldsymbol{y}}}{\boldsymbol{y}}}+c_{2})}{(u_{\tilde{\boldsymbol{y}}}^{2}+u_{\boldsymbol{y}}^{2}+c_{1})(\sigma_{\tilde{\boldsymbol{y}}}^{2}+\sigma_{\boldsymbol{y}}^{2}+c_{2})},$
(3)

where $\tilde{\boldsymbol{y}}$ denote the prediction. $u_{a}$ and
$\sigma_{a}^{2}$ represent the mean and variance of variable $a$,
respectively. $c_{1}$ and $c_{2}$ are two small values to avoid zero
denominator. The final loss is the combination of the two losses:

$\mathcal{L}_{f}({\tilde{\boldsymbol{y}}},{\boldsymbol{y}})=\alpha\cdot\mathbf{L}_{1}(\tilde{\boldsymbol{y}},\boldsymbol{y})+\beta\cdot\mathbf{SSIM}(\tilde{\boldsymbol{y}},\boldsymbol{y}),$
(4)

where $\alpha$ and $\beta$ are hyperparameters to balance the losses. The
$\mathbf{L}_{1}$ loss serves to suppress structural-independent Gaussian noise
while preserving the image’s brightness and color fidelity. Concurrently, the
$\mathbf{SSIM}$ loss plays a pivotal role in mitigating Poisson noise by
maintaining contrast in the high-frequency regions, thereby enhancing crucial
structural features for accurate denoising.

###  3.2 Low-level Masked Autoencoder Design

The model flowchart in Figure 1 demonstrates the straightforward
implementation of MAE learning paradigm for LDCT denoising. In the self-
pretraining, the LDCT images are subjected to patch-wise masking at a high
rate, typically $75\%$ as proposed in previous work [5]. These masked images
are then processed through an encoder and an arbitrary decoder to reconstruct
the original image. In fine-tuning, the lightweight decoder is replaced by a
specialized decoder tailored to match the target task, such as LDCT denoising
in our context. While this approach demonstrates its applicability in various
high-level vision tasks such as classification, segmentation, and object
detection, its effectiveness diminishes in the context of low-level vision
tasks like denoising.

To address this challenge, we introduce a novel masked image-based self-
pretraining approach, termed LoMAE, tailored for low-level vision tasks like
LDCT denoising. Unlike conventional methods, LoMAE (illustrated in Figure
1(b)) adopts a more streamlined architecture, consolidating the encoder and
decoder into a single MAE denoiser. This transformation is grounded in the
shared image-to-image training paradigm between LDCT denoising and MAE
pretraining. However, it is essential to emphasize the fundamental distinction
between these two tasks. Particularly, MAE pretraining focuses on
supplementing and restoring missing structures while the denoiser primarily
aims at Gaussian/Poisson noise removal. To bridge this gap and enhance
structural preservation, we redesign the pretraining model by adding a front-
to-end residual shortcut within the pretraining model. This redesign empowers
the model to better tackle denoising tasks while maintaining information
consistency and avoiding structural loss.

Particularly, in the pertaining stage, the model learns from masked LDCT to
whole LDCT images $\boldsymbol{x}$,

$\underset{\mathbf{W_{o}}}{\min}\ \
\mathcal{J}_{p}(\mathbf{W_{o}};\boldsymbol{x})=\mathbf{L}_{1}(\mathcal{D}^{-}(\mathbf{W_{o}};\boldsymbol{x}\cdot\boldsymbol{\mathcal{M}}),\boldsymbol{x}),$
(5)

where $\mathcal{D}^{-}(\mathbf{W_{o}};\boldsymbol{x})$ is the denoising model
without front-to-end shortcut, and $\mathbf{W_{o}}$ is a collection of initial
parameters before pretraining. After well-trained, let’s assume the parameters
become $\mathbf{W_{p}}$. $\boldsymbol{\mathcal{M}}$ represents the mask on the
patch-wise tokens from the images.

During finetuning, the model is trained from LDCT to NDCT image
$\boldsymbol{y}$ to learn the noise removal process. Here, we remove the image
masks and transfer the model weights from the pretraining stage, maintaining
an identical architecture. This finetuning stage can be mathematically
expressed as:

$\underset{\mathbf{W_{p}}}{\min}\ \
\mathcal{J}_{f}(\mathbf{W_{p}};\boldsymbol{x})=\mathcal{L}_{f}(\mathcal{D}^{-}(\mathbf{W_{p}};\boldsymbol{x})+\boldsymbol{x},\boldsymbol{y}).$
(6)

Figure 2: The architecture of the SwinIR and SUNet as MAE denoiser. Front-to-
end shortcuts are removed in the MAE pretraining stage and then are
reconnected in the finetuning stage.

###  3.3 Swin Transformer

The Swin transformer introduced a hierarchical attention mechanism featuring
shifted windows, enhancing contextual information fusion while controlling
computational costs. Initially, the image is divided into non-overlapping
patches, followed by linear mapping into tokens as transformer block input. In
a Swin transformer module (Figure 2(c)), input token sequence
$\mathbf{z}^{l-1}\in R^{b\times w\times n\times d_{o}}$ first undergoes
initial layer normalization (LN). Here, $b$ is batch size, $w$ is window
count, $n$ is token count, and $d_{o}$ represents token embedding dimension.
Then, $\mathbf{z}^{l-1}$ goes through linear transformations to derive query
(Q), key (K), and value (V) components to compute the window-based multi-head
attention (WMHA):

$\mathrm{WMHA}(\mathbf{Q},\mathbf{K},\mathbf{V})=\mathrm{softmax}(\frac{\mathbf{Q}\mathbf{K^{\top}}}{\sqrt{d_{k}}}+\mathbf{B})\mathbf{V},$
(7)

where $\frac{1}{\sqrt{{{d}_{k}}}}$ is the scaling factor. $\mathbf{B}$ is a
bias matrix. Subsequently, an additional LN layer, a multiple-layer perceptron
(MLP), and two residual shortcuts are employed to enhance feature inference.
As shown in Fig. 2 (c), a successive Swin transformer module is characterized
as

$\begin{cases}&\hat{\mathbf{z}}^{l}=\mathrm{WMHA}(\mathrm{LN}(\mathbf{z}^{l-1}))+\mathbf{z}^{l-1}\\\
&\mathbf{z}^{l}=\mathrm{MLP}(\mathrm{LN}(\hat{\mathbf{z}}^{l}))+\hat{\mathbf{z}}^{l}\\\
&\hat{\mathbf{z}}^{l+1}=\mathrm{S\text{-}WMHA}(\mathrm{LN}(\mathbf{z}^{l}))+\mathbf{z}^{l}\\\
&\mathbf{z}^{l+1}=\mathrm{MLP}(\mathrm{LN}(\hat{\mathbf{z}}^{l+1}))+\hat{\mathbf{z}}^{l+1}.\\\
\end{cases}$ (8)

Compared with the $\mathrm{WMHA}(\cdot)$ function,
$\mathrm{S\text{-}WMHA}(\cdot)$ incorporates an extra cyclic shift operation
to further enhance rich spatial information fusion. In addition, to reduce
computational complexity and support hierarchical representation, the patch
merging layers are interposed between successive transformer blocks to reduce
the feature map size by half.

In recent years, an abundance of deep learning models has harnessed Swin
transformer modules as their foundation. Similarly, the transformer models we
utilize to study the performance of the LoAME are also built upon Swin
modules, namely SwinIR [22] and SUNet [24], which represent two distinct
denoising architectures, as shown in Fig. 2(a) and (b). Both variants have
demonstrated formidable performance in addressing low-level denoising
challenges.

###  3.4 MAE Interpretation

While the MAE has been effective across diverse areas, interpreting its inner
workings remains a formidable challenging task as the majority of other deep
models. Naturally, our curiosity extends to understand how the model
complements the missing patches in the pretraining stage, thus facilitating
robust structural preservation. Among the handful methods for peering into the
workings of deep learning models, GradCAM [25] stands out as a favored choice
since it provides easy-to-understand saliency maps to elucidate what element
models are relying on to make a decision. GradCAM has found widespread utility
in a multitude of general computer vision tasks, spanning classification,
regression, segmentation, and more. Motivated by this, we have endeavored to
adapt the GradCAM methodology, here termed as MAE-GradCAM, to decode the
latent learning behavior of MAE. Please note that this method can be applied
to both the proposed LoMAE and original MAE since they share the same
pretraining procedure. Specifically, the traditional GradCAM generates the
saliency map $\boldsymbol{\mathcal{S}}$ through the following formula:

$\begin{cases}&\boldsymbol{\mathcal{S}}=\mathrm{ReLU}(\sum_{k}\boldsymbol{\mathcal{G}}^{k}_{c}\mathbf{\mathcal{A}}^{k})\\\
&\boldsymbol{\mathcal{G}}^{k}_{c}=\frac{1}{\mathrm{N}}\sum_{u,v}\frac{\partial
y^{c}}{\partial\mathbf{\mathcal{A}}^{k}_{uv}},\end{cases}$ (9)

where $y^{c}$ is the logits for the selected class. The selected layer is
typically the last convolution layer prior to the classification head.
However, these do not work in MAE because there is no concept of ’logits’ in
MAE and the last layer aims to restore the original pixels thus representing
low-level features with detailed structures.

To unravel the concealed mechanisms governing how MAE can fill the masked
patches, the MAE-GradCAM is distinct from the original GradCAM in two crucial
aspects. First, $\partial y^{c}$ is replaced by
$\partial\mathbf{L_{1}}(\hat{y}_{p},y_{p})$, where $\hat{y}_{(i,j)\in
p}\subset\hat{y}$ and $y_{(i,j)\in p}\subset y$ are prediction and target
pixels within certain interested patches $p$. In this sense, the gradients are
calculated within a defined region, allowing the MAE-GradCAM to unveil which
areas positively contribute to the generation of the missing patches. Second,
we choose the middle layer of the feature maps to visualize the gradients. In
contrast to classification models where the last convolution usually encodes
rich high-level semantic information, the last layer in the MAE training
typically encapsulates low-level features devoid of semantic content.
Consequently, we select the middle layer $\mathbf{\mathcal{A}}_{m}$ of base
network as our area of interest since it includes both high-level features and
spatial information. Therefore, the saliency map $\boldsymbol{\mathcal{S}}$
generated by MAE-GradCAM is defined as follows:

$\begin{cases}&\boldsymbol{\mathcal{S}}=\mathrm{ReLU}(\sum_{k}\boldsymbol{\mathcal{G}}^{k}\mathbf{\mathcal{A}}^{k}_{m})\\\
&\boldsymbol{\mathcal{G}}^{k}=\frac{1}{\mathrm{N}}\sum_{u,v}\frac{\partial\mathbf{L_{1}}(\hat{y}_{p},y_{p})}{\partial\mathbf{\mathcal{A}}^{k}_{uv}}.\end{cases}$
(10)

##  4 Experiments and Results

###  4.1 Dataset

Low-dose dataset: We employed the publicly available dataset from the 2016
NIH-AAPM-Mayo Clinic LDCT Grand
Challenge111https://www.aapm.org/GrandChallenge/LowDoseCT/ [26] for both model
training and evaluation. This dataset encompasses 2,378 CT images with a slice
thickness of 3.0 mm, collected from ten distinct patients. To ensure
robustness, we implemented a ten-fold cross-validation strategy, where each
fold reserved one patient for testing and the remaining for training. Due to
GPU limitations, the original $512\times 512$ images were uniformly resized to
$256\times 256$ using inter_area interpolation. Furthermore, techniques such
as image rotation and flipping are utilized for data augmentation.

Simulated dataset: To test the robustness and generalizability of the proposed
model, we simulate multiple LDCT scans with different dose level based on
[26]. To mimic the realistic clinical environment, we first obtain the
projection datasets of Radon transform from the Manifold and Graph Integrative
Convolutional Network (MAGIC) [27]. Then, the CT images of different dose
level are synthesized by adding Poisson and electronic noise to the clean
projections. Here, the geometry of the scan is listed as follows: The distance
from the x-ray source to the isocenter of image is 595 mm. The distance from
the detector to the source is 1085.6 mm. The size of the image pixel is 0.6641
mm. The number of detector units is 1024, with a unit pitch of 1.2854 mm. With
the obtained noise-free projections, the noisy projections are obtained by the
following formula:

$\boldsymbol{y}=\mathrm{ln}\frac{\mathrm{I}_{0}}{\mathrm{Poisson}(\mathrm{I}_{0}\mathrm{exp}(-\hat{\boldsymbol{y}}))+\mathrm{Normal}(0,\sigma_{e}^{2})},$
(11)

where $\mathrm{I}_{0}$ represents the x-ray intensity. $\hat{y}$ is the noise-
free sinogram and $\sigma_{e}^{2}$ is the variance in the normal distribution.
In our experiment, the electronic noise variance is set to 10 according to
[28, 29]. Then, the noisy CT images are reconstructed using filtered
backprojection algorithm in MAGIC.

###  4.2 Experiment Settings

Our experiments are conducted on Ubuntu 18.04.5 LTS system equipped with an
Intel(R) Core(TM) i9-9920X CPU @ 3.50GHz. All models are implemented with
PyTorch 1.5.0 and CUDA 10.2.0 on four NVIDIA 2080TI 11G GPUs for accelerated
computations. Here are the experiment settings:

  * •

The masked patch size is 8 with a masked rate of $75\%$ for both SwinIR and
SUNet during pretraining.

  * •

The SwinIR model was configured with block lengths of [4,4,4,4] with a
constant embedding dimension of 60 and 6 heads for the multiple head
attention.

  * •

The model depths for the SUNet are [2,2,2,2] for both encoder and decode with
a uniform embedding dimension of 80 and 8 heads for attention.

  * •

In both the LoMAE pretraining and finetuning phases, we optimize the model
over 50 epochs using the Adam optimizer. We initiate the learning rate at
$1.5\times 10^{-4}$ and implement a scheduled learning rate decay of 0.5 every
3,000 iterations.

  * •

For both LoMAE pretraining and finetuning processes, the batch size is set to
1.

###  4.3 Result Comparison

To assess the effectiveness of LoMAE pretraining for enhancing transformer
denoising model, the following state-of-the-art models are investigated:
REDCNN [1], MAPNN [3], SUNet [24], and SwinIR [22]. Specifically, all these
methods are rigorously implemented and optimized according to the officially
disclosed codebases. Furthermore, we explore the integration of LoMAE as a
pretraining strategy into SwinIR (SwinIR+LoMAE) and SUNet (SUNet+LoMAE).

Figure 3: The denoised results of Lesion 575 from patient L506. (a) LDCT, (b)
REDCNN, (c) MAPNN, (d) SUNet, (e) SUNet+LoMAE, (f) SwinIR, (g) SwinIR+LoMAE,
(h) NDCT. The display window is [-160,240] HU.

Results in Fig. 3 show that all these methods are effective in removing the
noises from LDCT images. The SwinIR and SUNet can generate clearer noise-free
images than the REDCNN and MAPNN. Upon closer examination of the regions of
interest (ROIs) in Fig. 3, it’s evident that our proposed approaches,
SwinIR+LoMAE and SUNet+LoMAE, exhibit superior performance in generating more
accurate and smoother lesion structures while retaining essential details when
compared to the original models. Moving on to Fig. 4, we provide
visualizations of horizontal and vertical intensity profiles for a
representative slice. Results indicate that the overall trend of the intensity
profiles of the studied methods are close, whereas the proposed methods
(SwinIR+LoMAE and SUNet+LoMAE) generate the most accurate structures, closely
aligning with the NDCT image. Furthermore, the zoomed-in view and bar plot in
Fig. 4 highlight that our models exhibit the lowest mean absolute error in
terms of intensities compared to others. In addition, as illustrated in Fig.
5, the loss curves of the SwinIR+LoMAE and SUNet+LoMAE display a smoother
trajectory compared to SwinIR and SUNet, respectively. This indicates that the
LoMAE methods benefit from a much more stable and consistent training process
than the original approaches.

Figure 4: The intensity profiles of different methods and the corresponding
mean absolute errors on a representative slice. Notably, the left two figures
feature a zoomed-in view that indicates the pixel-wise mean absolute error
when compared to NDCT within the region of interest.  Figure 5: The loss
curves of the SwinIR and SUNet with their corresponding LoMAE methods.

In addition to qualitative assessments, our model undergoes quantitative
evaluation using two metrics including structural similarity (SSIM) and root
mean square error (RMSE). Results in Tab. 1 indicate that the transformer
models (SwinIR and SUNet) can achieve better performance than the
convolutional models (REDCNN and MAPNN). Notably, SwinIR achieves the highest
scores in both metrics. Furthermore, the incorporation of LoMAE into SwinIR
and SUNet leads to substantial quantitative improvements. Specifically, we
observe an increase of 0.0063 and 0.0046 in SSIM, as well as reductions of
0.4877 and 0.4795 in RMSE for SwinIR+LoMAE and SUNet+LoMAE, respectively. In
summary, our proposed LoMAE significantly enhances the denoising performance
of transformer-based models, both qualitatively and quantitatively.

Table 1: Quantitative results of different methods on L$506$. The bold-faced
numbers are the best results.

Method | SSIM$\uparrow$ | RMSE$\downarrow$
---|---|---
LDCT | 0.9359 | 10.4833
REDCNN | 0.9501$\pm{0.0012}$ | 7.8016$\pm{0.0834}$
MAPNN | 0.9504$\pm{0.0019}$ | 7.5357$\pm{0.1342}$
SUNet | 0.9555$\pm{0.0010}$ | 7.3147$\pm{0.1084}$
SUNet+LoMAE* | 0.9601$\pm{0.0009}$ | 6.8352$\pm{0.1022}$
SwinIR | 0.9546$\pm{0.0011}$ | 7.2232$\pm{0.0981}$
SwinIR+LoMAE* | 0.9609$\pm{0.0007}$ | 6.7355$\pm{0.0800}$

Figure 6: The quantitative results of SSIM and RMSE using different numbers of
patients for training. Figure 7: The quantitative results of different methods
on all testing slices from patient L506. Please note that the solid line ’-’
indicates the SwinIR+LoMAE method, while the dashed line ’- -’ denotes SwinIR.

###  4.4 Ground Truth Data Dependency

In clinical applications, obtaining well-organized paired data for supervised
deep models can be exceedingly challenging, if not impossible. This difficulty
arises from the impracticality and ethical concerns associated with subjecting
patients to multiple radiation-based tests. Even if such practices were
feasible, the dynamic nature of the human body renders this approach far from
ideal. Consequently, there is a growing need for models that require less
ground truth data, aligning more closely with clinical realities. Our model is
tailored to meet this demand.

In the pretraining stage, the proposed method relies solely on noisy data.
Only in the fine-tuning stage, both noisy and ground truth clean data become
necessary to support a supervised training paradigm. In this context, the
proposed model can leverage the noisy data to encode the structural dependence
across various anatomical regions and thus potentially demands less ground
truth data to learn meaningful representations.

To empirically validate this distinctive feature, SwinIR is studied since it
exhibits superior denoising performance as demonstrated in Subsection C. We
randomly select different numbers of patients for model training and then
evaluate the models performances under the varying data configurations. As
shown in Fig. 6, as the number of patients in the training data increases,
both SSIM and RMSE metrics show noticeable improvement for both SwinIR and
SwinIR+LoMAE models. Generally speaking, the curve of SwinIR+LoMAE method lies
above SwinIR’s in SSIM and under in RMSE. To simplify notation, the number
following each method indicates the quantity of patients with ground-truth
data used for model training. For instance, ’LoMAE2’ denotes that the model
was trained with data from only two patients with clean images. It’s
particularly noteworthy that even with only two patients’ worth of clean data,
SwinIR+LoMAE consistently outperforms SwinIR, clearly demonstrating its
advantage. Furthermore, t-tests are performed between LoMAE2 and SwinIR
methods. The LoMAE2 method shows significant improvement ($p<0.05$) over
SwinIR8 ($p=6.07\times 10^{-3}$), SwinIR6 ($p=5.99\times 10^{-3}$), SwinIR4
($p=2.56\times 10^{-3}$), and SwinIR2 ($p=1.22\times 10^{-3}$).

Figure 8: The MAE reconstruction results on the testing cases. The display
window for the gray images is [-160,240] HU. Figure 9: The MAE-GradCAM results
from four representative slices. Each row represents one slice, while the
three columns denote three different locations. The red arrows indicate the
interested $8\times 8$ patch areas for the MAE-GradCAM. It is better to view
this image in a zoomed version.

Moreover, results of all slices in Fig. 7 provide a visual representation of
these findings. It illustrates that the majority of curves of the SwinIR+LoMAE
are better than that of SwinIR, which indicates that the SwinIR+LoMAE method
can greatly reduce its reliance on ground truth data. This attribute makes it
a highly practical choice for clinical scenarios where noisy data is abundant
while labeled data remains scarce.

Figure 10: Representative testing results of the same trained models under
different noise levels (x-ray intensity). The display window is [-160,240] HU.

###  4.5 MAE Interpretation

In the pretraining phase, the MAE adeptly learns the mapping from corrupted
images to their original counterparts. The results, as depicted in Figure 8,
demonstrate the remarkable capability of the MAE, indicating that even when
only 25% of the original image is available, it can effectively reconstruct
the entire legitimate image with remarkable accuracy. While minor details are
missing, the overall structural integrity and anatomical regions remain
largely preserved post-pretraining, underscoring the efficacy of this
pretraining process.

Furthermore, with the proposed MAE-GradCAM, some latent behaviors of the MAE
pretraining can be decoded. As shown from the representative slices in Fig. 9,
although the exact causal relationship of the patch generation process is
inevident, several common behaviors can be unveiled:

  * •

The MAE can clearly discriminate the foreground body tissue and the background
areas. In all exemplary slices, the attention is distinctly focused either on
the foreground or background.

  * •

The saliency map with similar weights indicates similar tissues or organs,
suggesting that MAE potentially possesses remarkable semantic segmentation
capability. For example in Fig. 9 (b1, b2, c3, d1, d3), the attention on the
liver and kidney is consistently distributed across the entire organs.

  * •

The completeness of the missing background is typically highly affected by
other existing background patches. For instance in Fig. 9 (c1), the missing
background patches rely on the presence of other background areas.

  * •

The boundaries between low-intensity and high-intensity are predominantly
determined by the high-intensity organs and partially determined by the low-
intensity areas. For example in Fig. 9 (b1, b2, b3, c3, d1, d2, d3), the
intersected area is primarily determined by the high-intensity area like
aorta, rib, and ilium, with a less influence from low-intensity areas like
liver, kidney, and spleen.

###  4.6 Model Robustness and Generalizability

Deep learning models excel at learning complicated latent functions. However,
their performance often hinges on the assumption that testing data follows the
same distribution as the training data. Controversially, the data in the real
world is much more diverse and complex with many unpredictable perturbations,
thus, posing a severe challenge to a model’s robustness and generalizability.
In our case, a clinically favorable denoising model should deliver constant
good denoising performance across various noise distributions. Naturally, we
are motivated to study the robustness and generalizability of the proposed
method. Specifically, we evaluate how the quarter dose-trained model performs
over different noise distributions.

Figure 11: The quantitative curves of different methods under different noise
levels. Figure 12: Comparison of the NPS maps SwinIR and SUNet with respect to
the corresponding LoMAE methods under different noise levels. The image center
of each map represents low frequency structural noise while the edges and
corners indicate the high frequency noise such as Poisson and Gaussian.
Figure 13: The cumulative outcomes of the NPS curves along the vertical axis.
The left figure displays the results of SwinIR and SwinIR+LoMAE while the
right figure indicates SUNet and SUNet+LoMAE.

Experiments are conducted on the Simulated Dataset where a variety of noise
settings is performed. The noise level is mainly determined by the x-ray
intensities (dose level), _e.g._ , $1\times 10^{5}$, $8\times 10^{4}$,
$6\times 10^{4}$, $4\times 10^{4}$, and $2\times 10^{4}$, which correspond to
$10\%$, $8\%$, $6\%$, $4\%$, and $2\%$ doses, respectively. As shown in Fig.
10, the studied methods exhibit denoising capabilities across different noise
levels. Unsurprisingly, all methods’ performance drops when the noise level
increases, which is expected since all models are only trained on 25% dose
images. Notably, our method exhibits a slower decline in denoising capability
and constantly delivers the best results across all tests as evidenced by both
qualitative results in Fig. 10 and quantitative curves in Fig. 11. For example
in Fig. 10(a), the LoMAE methods capture the hepatic vein structure with high
fidelity. In Fig. 10(b) where there exists a solid cystic mass in the left
kidney, the LoMAEs generate smoother interior structures. Next in Fig. 10(c),
LoMAE can deliver the most accurate colon image with minimal distortion.
Finally, Fig. 10(d) and Fig. 10(e) illustrate that our models render clearer
muscle and kidney tissues compared to others. Particularly at the lowest
intensity of $\mathrm{I}_{0}=2\times 10^{4}$, all other methods fail to work
but our method, especially the SwinIR+LoMAE, can still suppress the major
noises and retain the key structures.

Moreover, the noise power spectrum (NPS) maps of SwinIR and SUNet are
visualized along with the corresponding LoMAE-enhanced variants in Fig. 12.
These maps indicate that the LoMAE methods are subject to less noise across
different noise frequencies, especially on high-frequency noise. The
horizontally summed results in Fig. 13 further underscore the superior noise
reduction and structure preservation capabilities of LoMAE methods compared to
the raw approaches. Particularly, the low-frequency noise (middle of Fig. 13
on x axis) is most pronounced in SwinIR ($2\times 10^{4}$) and SUNet ($2\times
10^{4}$), indicating significant structural losses in these cases.
Controversially, the corresponding low-frequency noise is significantly lower
in the LoMAE-enhanced methods.

Figure 14: The CKA results on SwinIR and SwinIR+LoMAE. CKAs are calculated
between the last-layer feature maps of input images with different noise
levels.

To further demonstrate the generalizability of the proposed LoMAE, centered
kernel alignment (CKA) [30, 31] is performed on the feature maps of SwinIR and
LoMAE. Specifically, let’s consider $\mathbf{X}_{1}\in R^{m\times d}$ and
$\mathbf{X}_{2}\in R^{m\times d}$ as two token layers to compute the CKA. The
Gram matrices is denoted as $\mathbf{K}=\mathbf{X}_{1}\mathbf{X}_{1}^{\top}$
and $\mathbf{L}=\mathbf{X}_{2}\mathbf{X}_{2}^{\top}$, the CKA is expressed as:

$\mathrm{CKA}(\mathbf{K},\mathbf{L})=\frac{\mathrm{HSIC}(\mathbf{K},\mathbf{L})}{\sqrt{\mathrm{HSIC(\mathbf{K},\mathbf{K})}\mathrm{HSIC}(\mathbf{L},\mathbf{L})}}.$
(12)

Here, $\mathrm{HSIC}(\cdot)$ represents the Hilbert-Schhmidt Independence
Criterion, which is calculated as below:

$\mathrm{HSIC}(\mathbf{K},\mathbf{L})=\frac{vec(\mathbf{K}^{\prime})\cdot
vec(\mathbf{L}^{\prime})}{(m-1)^{2}},$ (13)

where $\mathbf{K}^{\prime}=\mathbf{H}\mathbf{K}\mathbf{H}$ and
$\mathbf{L}^{\prime}=\mathbf{H}\mathbf{L}\mathbf{H}$ are the centered Gram
metrics of $\mathbf{K}$ and $\mathbf{L}$, respectively.
$\mathbf{H}=\mathbf{I}_{n}-\frac{1}{n}\mathbf{1}_{n}\mathbf{1}_{n}^{\top}$ is
the centering metric. And $vec(\cdot)$ denotes the vectorization of a matrix.

In this analysis, CKA is calculated among the feature maps over various noise
settings. Higher CKA indicates higher generalizability across multiple noise
levels. As shown in Fig. 14, the LoMAE method demonstrates significantly
higher CKAs than raw SwinIR. Particularly, the CKA between $2\times 10^{4}$
and $1\times 10^{5}$ is 0.24 with SwinIR, but 0.76 with LoMAE. Therefore, it
can be concluded that the proposed LoMAE possesses strong robustness and
generalizability across different dose levels.

###  4.7 Ablation Studies

On the impact of base model: In this study, we primarily utilize SwinIR and
SUNet for LOMAE applications. Nevertheless, this strategy can also be applied
to non-swin based transformers. As an example, the LoMAE is tested on
Restormer [32], another SOTA transformer model for denoising. The original
Restormer achieves SSIM and RMSE values of $0.9552\pm{0.0010}$ and
$7.4543\pm{0.1084}$, respectively. With Restormer+LoMAE, the SSIM and RMSE
improve to $0.9582\pm{0.0009}$ and $6.9881\pm{0.1022}$, respectively. Those
results further confirm the efficacy of the LoMAE on raw transformer models.

On the impact of the front-to-end shortcut: Front-to-end shortcut is absent in
LoMAE pretaining but utilized in finetuning to ensure information consistency.
To evaluate its usefulness, we remove the shortcut in finetuning as a
comparison. Quantitative results on SUNet show that without the shortcut, the
SSIM and RMSE are $0.9571\pm{0.0012}$ and $7.0627\pm{0.1000}$, which are
inferior to the results of $0.9601\pm{0.0008}$ and $6.8352\pm{0.0900}$
obtained with the shortcut.

##  5 Conclusion

This paper explored the novel application of Masked Autoencoder in LDCT
denoising, marking the pioneering use of MAE as a self-pretraining strategy in
low-level denoising tasks. Specifically, we proposed a LoMAE design and
developed an MAE-GradCAM method to interpret the latent learning behavior.
Experimental results on the AAPM and simulated dataset underscored LoMAE’s
exceptional ability to enhance noise removal and structural preservation,
producing qualitatively and quantitatively superior CT images. Remarkably,
LoMAE proved effective in semi-supervised contexts, significantly reducing
reliance on ground truth data. Furthermore, our findings demonstrated LoMAE’s
high robustness and generalizability across a spectrum of noise levels, paving
the way for potential applications in various medical image modalities. Our
future research endeavors will focus on further exploring LoMAE’s
applicability in these domains.

## References

  * [1] H. Chen, Y. Zhang, M. K. Kalra, F. Lin, Y. Chen, P. Liao, J. Zhou, and G. Wang, “Low-dose ct with a residual encoder-decoder convolutional neural network,” IEEE transactions on medical imaging, vol. 36, no. 12, pp. 2524–2535, 2017.
  * [2] F. Fan, H. Shan, M. K. Kalra, R. Singh, G. Qian, M. Getzin, Y. Teng, J. Hahn, and G. Wang, “Quadratic autoencoder (q-ae) for low-dose ct denoising,” IEEE transactions on medical imaging, vol. 39, no. 6, pp. 2035–2050, 2019.
  * [3] H. Shan, A. Padole, F. Homayounieh, U. Kruger, R. D. Khera, C. Nitiwarangkul, M. K. Kalra, and G. Wang, “Competitive performance of a modularized deep neural network compared to commercial algorithms for low-dose ct image reconstruction,” Nature Machine Intelligence, vol. 1, no. 6, pp. 269–276, 2019.
  * [4] D. Wang, F. Fan, Z. Wu, R. Liu, F. Wang, and H. Yu, “Ctformer: Convolution-free token2token dilated vision transformer for low-dose ct denoising,” arXiv preprint arXiv:2202.13517, 2022.
  * [5] K. He, X. Chen, S. Xie, Y. Li, P. Dollár, and R. Girshick, “Masked autoencoders are scalable vision learners,” in Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp. 16000–16009, 2022.
  * [6] Z. Xie, Z. Zhang, Y. Cao, Y. Lin, J. Bao, Z. Yao, Q. Dai, and H. Hu, “Simmim: A simple framework for masked image modeling,” in Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp. 9653–9663, 2022.
  * [7] C. Feichtenhofer, Y. Li, K. He, et al., “Masked autoencoders as spatiotemporal learners,” Advances in neural information processing systems, vol. 35, pp. 35946–35958, 2022.
  * [8] M. Germain, K. Gregor, I. Murray, and H. Larochelle, “Made: Masked autoencoder for distribution estimation,” in International conference on machine learning, pp. 881–889, PMLR, 2015.
  * [9] P. Gao, T. Ma, H. Li, Z. Lin, J. Dai, and Y. Qiao, “Convmae: Masked convolution meets masked autoencoders,” arXiv preprint arXiv:2205.03892, 2022.
  * [10] L. Zhou, H. Liu, J. Bae, J. He, D. Samaras, and P. Prasanna, “Self pre-training with masked autoencoders for medical image analysis,” arXiv preprint arXiv:2203.05573, 2022.
  * [11] S. Tihon, M. U. Javaid, D. Fourure, N. Posocco, and T. Peel, “Daema: Denoising autoencoder with mask attention,” in International Conference on Artificial Neural Networks, pp. 229–240, Springer, 2021.
  * [12] Q. Wu, H. Ye, Y. Gu, H. Zhang, L. Wang, and D. He, “Denoising masked autoencoders are certifiable robust vision learners,” arXiv preprint arXiv:2210.06983, 2022.
  * [13] F. Zhao, D. Li, R. Luo, M. Liu, X. Jiang, and J. Hu, “Self-supervised deep learning for joint 3d low-dose pet/ct image denoising,” Computers in Biology and Medicine, p. 107391, 2023.
  * [14] D. Wang, Y. Xu, S. Han, and H. Yu, “Masked autoencoders for low-dose ct denoising,” in 2023 IEEE 20th International Symposium on Biomedical Imaging (ISBI), pp. 1–4, IEEE, 2023.
  * [15] D. Wang, Y. Xu, S. Han, and H. Yu, “Masked autoencoders for low dose ct denoising,” arXiv preprint arXiv:2210.04944, 2022.
  * [16] F.-L. Fan, D. Wang, H. Guo, Q. Zhu, P. Yan, G. Wang, and H. Yu, “On a sparse shortcut topology of artificial neural networks,” IEEE Transactions on Artificial Intelligence, vol. 3, no. 4, pp. 595–608, 2021.
  * [17] Z. Jia, Y. Luo, D. Wang, Q. N. Dinh, S. Lin, A. Sharma, E. M. Block, M. Yang, T. Gu, A. J. Pearlstein, et al., “Nondestructive multiplex detection of foodborne pathogens with background microflora and symbiosis using a paper chromogenic array and advanced neural network,” Biosensors and Bioelectronics, vol. 183, p. 113209, 2021.
  * [18] B. Morovati, R. Lashgari, M. Hajihasani, and H. Shabani, “Reduced deep convolutional activation features (r-decaf) in histopathology images to improve the classification performance for breast cancer diagnosis,” arXiv preprint arXiv:2301.01931, 2023.
  * [19] H. Chen, Y. Zhang, M. K. Kalra, F. Lin, Y. Chen, P. Liao, J. Zhou, and G. Wang, “Low-dose ct with a residual encoder-decoder convolutional neural network,” IEEE transactions on medical imaging, vol. 36, no. 12, pp. 2524–2535, 2017.
  * [20] L. Yang, D. Zhang, et al., “Low-dose ct denoising via sinogram inner-structure transformer,” arXiv preprint arXiv:2204.03163, 2022.
  * [21] Z. Liu, Y. Lin, Y. Cao, H. Hu, Y. Wei, Z. Zhang, S. Lin, and B. Guo, “Swin transformer: Hierarchical vision transformer using shifted windows,” arXiv preprint arXiv:2103.14030, 2021.
  * [22] J. Liang, J. Cao, G. Sun, K. Zhang, L. Van Gool, and R. Timofte, “Swinir: Image restoration using swin transformer,” in Proceedings of the IEEE/CVF International Conference on Computer Vision, pp. 1833–1844, 2021.
  * [23] M. Puttaguntaa, R. Subbanb, and N. K. B. Cc, “Swinir transformer applied for medical image super-resolution,” Procedia Computer Science, vol. 204, pp. 907–913, 2022.
  * [24] C.-M. Fan, T.-J. Liu, and K.-H. Liu, “Sunet: swin transformer unet for image denoising,” in 2022 IEEE International Symposium on Circuits and Systems (ISCAS), pp. 2333–2337, IEEE, 2022.
  * [25] R. R. Selvaraju, M. Cogswell, A. Das, R. Vedantam, D. Parikh, and D. Batra, “Grad-cam: Visual explanations from deep networks via gradient-based localization,” in Proceedings of the IEEE international conference on computer vision, pp. 618–626, 2017.
  * [26] C. H. McCollough, A. C. Bartley, R. E. Carter, B. Chen, T. A. Drees, P. Edwards, D. R. Holmes III, A. E. Huang, F. Khan, and S. Leng, “Low‐dose ct for the detection and classification of metastatic liver lesions: results of the 2016 low dose ct grand challenge,” Medical physics, vol. 44, no. 10, pp. e339–e352, 2017.
  * [27] W. Xia, Z. Lu, Y. Huang, Z. Shi, Y. Liu, H. Chen, Y. Chen, J. Zhou, and Y. Zhang, “Magic: Manifold and graph integrative convolutional network for low-dose ct reconstruction,” IEEE Transactions on Medical Imaging, vol. 40, no. 12, pp. 3459–3472, 2021.
  * [28] W. Xia, Z. Lu, Y. Huang, Y. Liu, H. Chen, J. Zhou, and Y. Zhang, “Ct reconstruction with pdf: Parameter-dependent framework for data from multiple geometries and dose levels,” IEEE Transactions on Medical Imaging, vol. 40, no. 11, pp. 3065–3076, 2021.
  * [29] S. Niu, Y. Gao, Z. Bian, J. Huang, W. Chen, G. Yu, Z. Liang, and J. Ma, “Sparse-view x-ray ct reconstruction via total generalized variation regularization,” Physics in Medicine & Biology, vol. 59, no. 12, p. 2997, 2014.
  * [30] C. Cortes, M. Mohri, and A. Rostamizadeh, “Algorithms for learning kernels based on centered alignment,” The Journal of Machine Learning Research, vol. 13, no. 1, pp. 795–828, 2012.
  * [31] M. Raghu, T. Unterthiner, S. Kornblith, C. Zhang, and A. Dosovitskiy, “Do vision transformers see like convolutional neural networks?,” Advances in Neural Information Processing Systems, vol. 34, pp. 12116–12128, 2021.
  * [32] S. W. Zamir, A. Arora, S. Khan, M. Hayat, F. S. Khan, and M.-H. Yang, “Restormer: Efficient transformer for high-resolution image restoration,” in Proceedings of the IEEE/CVF conference on computer vision and pattern recognition, pp. 5728–5739, 2022.


This research was supported in part by the Israel Science Foundation (grant
No. 1122/22). Xin Huang was supported in part at the Technion Israel Institute
of Technology by an Aly Kaufman Fellowship.

## References

  * [1] Amanatidis, G., Markakis, E., Nikzad, A., and Saberi, A. Approximation algorithms for computing maximin share allocations.  In International Colloquium on Automata, Languages, and Programming (2015), pp. 39–51.
  * [2] Aziz, H., Caragiannis, I., Igarashi, A., and Walsh, T. Fair allocation of indivisible goods and chores.  In Proceedings of the Twenty-Eighth International Joint Conference on Artificial Intelligence, IJCAI 2019, Macao, China, August 10-16, 2019 (2019), S. Kraus, Ed., ijcai.org, pp. 53–59.
  * [3] Aziz, H., Li, B., Moulin, H., and Wu, X. Algorithmic fair allocation of indivisible items: A survey and new questions.  arXiv preprint arXiv:2202.08713 (2022).
  * [4] Aziz, H., Li, B., and Wu, X. Approximate and strategyproof maximin share allocation of chores with ordinal preferences.  Mathematical Programming (2022), 1–27.
  * [5] Aziz, H., Rauchecker, G., Schryen, G., and Walsh, T. Algorithms for max-min share fair allocation of indivisible chores.  In Thirty-First AAAI Conference on Artificial Intelligence (2017).
  * [6] Babaioff, M., Ezra, T., and Feige, U. Best-of-both-worlds fair-share allocations.  CoRR abs/2102.04909 (2021).
  * [7] Babaioff, M., Ezra, T., and Feige, U. Fair-share allocations for agents with arbitrary entitlements.  In Proceedings of the 22nd ACM Conference on Economics and Computation (2021), pp. 127–127.
  * [8] Barman, S., and Krishnamurthy, S. K. Approximation algorithms for maximin fair division.  ACM Trans. Economics and Comput. 8, 1 (2020), 5:1–5:28.
  * [9] Barman, S., Krishnamurthy, S. K., and Vaish, R. Finding fair and efficient allocations.  In Proceedings of the 19th ACM Conference on Economics and Computation (EC) (2018), pp. 557–574.
  * [10] Berger, B., Cohen, A., Feldman, M., and Fiat, A. Almost full EFX exists for four agents.  In Thirty-Sixth AAAI Conference on Artificial Intelligence, AAAI 2022 (2022), AAAI Press, pp. 4826–4833.
  * [11] Bogomolnaia, A., and Moulin, H. A new solution to the random assignment problem.  J. Econ. Theory 100, 2 (2001), 295–328.
  * [12] Bouveret, S., and Lemaître, M. Characterizing conflicts in fair division of indivisible goods using a scale of criteria.  Autonomous Agents and Multi-Agent Systems 30, 2 (2016), 259–290.
  * [13] Budish, E. The combinatorial assignment problem: Approximate competitive equilibrium from equal incomes.  Journal of Political Economy 119, 6 (2011), 1061–1103.
  * [14] Caragiannis, I., Kurokawa, D., Moulin, H., Procaccia, A. D., Shah, N., and Wang, J. The unreasonable fairness of maximum nash welfare.  In Proceedings of the 17th ACM Conference on Economics and Computation (EC) (2016), pp. 305–322.
  * [15] Chaudhury, B. R., Cheung, Y. K., Garg, J., Garg, N., Hoefer, M., and Mehlhorn, K. Fair division of indivisible goods for a class of concave valuations.  Journal of Artificial Intelligence Research 74 (2022), 111–142.
  * [16] Chaudhury, B. R., Garg, J., and Mehlhorn, K. EFX exists for three agents.  In EC ’20: The 21st ACM Conference on Economics and Computation, Virtual Event, Hungary, July 13-17, 2020 (2020), P. Biró, J. Hartline, M. Ostrovsky, and A. D. Procaccia, Eds., ACM, pp. 1–19.
  * [17] Chaudhury, B. R., Garg, J., Mehlhorn, K., Mehta, R., and Misra, P. Improving efx guarantees through rainbow cycle number.  In Proceedings of the 22nd ACM Conference on Economics and Computation (2021), pp. 310–311.
  * [18] Feige, U., Sapir, A., and Tauber, L. A tight negative example for mms fair allocations.  In International Conference on Web and Internet Economics (2021), Springer, pp. 355–372.
  * [19] Garg, J., McGlaughlin, P., and Taki, S. Approximating maximin share allocations.  In 2nd Symposium on Simplicity in Algorithms, SOSA@SODA 2019, January 8-9, 2019 - San Diego, CA, USA (2019), J. T. Fineman and M. Mitzenmacher, Eds., vol. 69 of OASICS, Schloss Dagstuhl - Leibniz-Zentrum für Informatik, pp. 20:1–20:11.
  * [20] Garg, J., and Taki, S. An improved approximation algorithm for maximin shares.  In EC ’20: The 21st ACM Conference on Economics and Computation, Virtual Event, Hungary, July 13-17, 2020 (2020), P. Biró, J. Hartline, M. Ostrovsky, and A. D. Procaccia, Eds., ACM, pp. 379–380.
  * [21] Ghodsi, M., Hajiaghayi, M., Seddighin, M., Seddighin, S., and Yami, H. Fair allocation of indivisible goods: Improvements and generalizations.  In Proceedings of the 19th ACM Conference on Economics and Computation (EC) (2018), pp. 539–556.
  * [22] Ghodsi, M., Hajiaghayi, M., Seddighin, M., Seddighin, S., and Yami, H. Fair allocation of indivisible goods: Beyond additive valuations.  Artificial Intelligence 303 (2022), 103633.
  * [23] Huang, X., and Lu, P. An algorithmic framework for approximating maximin share allocation of chores.  In Proceedings of the 22nd ACM Conference on Economics and Computation (2021), pp. 630–631.
  * [24] Kurokawa, D., Procaccia, A. D., and Wang, J. Fair enough: Guaranteeing approximate maximin shares.  Journal of the ACM 65(2) (2018), 8:1–27.
  * [25] Li, B., Wang, F., and Zhou, Y. Fair allocation of indivisible chores: Beyond additive valuations.  arXiv preprint arXiv:2205.10520 (2022).
  * [26] Livanos, V., Mehta, R., and Murhekar, A. (Almost) Envy-free, proportional and efficient allocations of an indivisible mixed manna.  In Proceedings of the 21st International Conference on Autonomous Agents and Multiagent Systems (2022), pp. 1678–1680.