# Transforming photon statistics through zero-photon subtraction

C. M. Nunn  S. U. Shringarpure  T. B. Pittman  Physics Department, University
of Maryland Baltimore County, Baltimore, MD 21250

###### Abstract

Zero-photon subtraction (ZPS) is a conditional measurement process that can
reduce the mean photon number of quantum optical states without physically
removing any photons. Here we show that ZPS can also be used to transform
certain super-Poissonian states into sub-Poissonian states, and vice versa.
Combined with a well-known “no-go” theorem on conditional measurements, this
effect leads to a new set of non-classicality criteria that can be
experimentally tested through ZPS measurements.

##  I Introduction

The investigation of super-Poissonian and sub-Poissonian light sources plays a
significant role in the history of quantum optics [1, 2]. The photon-number
distributions of these sources display variances that are, respectively, wider
or narrower than the benchmark Poissonian statistics of a coherent state with
the same average photon number [3]. Experimentally, these properties can be
conveniently characterized by Mandel’s $Q$ parameter, with $Q>0$ for super-
Poissonian sources and ${-1\leq Q<0}$ for sub-Poissonian sources [4].

Given either type of source, it is also interesting to consider the physical
processes that can be used to actively transform the emitted super-Poissonian
light into sub-Poissonian light, and vice versa. Quintessential examples
include the use of optical nonlinearities such as two-photon absorption [5] or
photon blockades [6] to transform super- into sub-Poissonian light, and
amplification [7, 8] or phase randomization [9] to transform sub- into super-
Possonian light. More recently, it has been shown that conditional measurement
processes such as photon addition and photon subtraction can also be used to
implement these statistical transformations [10, 11, 12, 13]. Here we show
that the relatively new conditional measurement process of zero-photon
subtraction (ZPS) can also transform certain sub-Poissonian states into super-
Poissonian states (and vice versa), despite no photons being added to or
subtracted from the system. Other related work includes performing statistical
transformations by manipulating the initial parameters in atomic fluorescence
[14], sideband squeezing [15], photon catalysis [16], and displacement
operation [17] experiments.

Figure 1: Overview of a typical zero-photon subtraction (ZPS) setup. An input
state $\hat{\rho}_{in}$ enters one input port of a beamsplitter with variable
reflectivity $R$, and the detection of zero reflected photons at $D_{1}$
heralds a noiselessly attenuated output state $\hat{\rho}_{out}$. Despite no
photons being removed from the system, ZPS can dramatically transform the
photon statistics of certain input states. These effects can be observed and
quantified by measuring trends in a relative attenuation parameter $K(R)$
using an auxiliary detector $D_{2}$.

Figure 1 shows an overview of a typical ZPS setup. ZPS was originally proposed
in the context of noiseless attenuation [18], and its counterintuitive ability
to reduce the mean photon number of input states has now been demonstrated in
several experiments [19, 20, 21, 22, 23, 24]. As illustrated in the figure, an
input state $\hat{\rho}_{in}$ with average photon number
$\langle\hat{n}\rangle_{in}$ passes through a beamsplitter with variable
reflectance $R$, and the output state is only accepted when no photons are
reflected. This is accomplished by actively “heralding on zero” photons (HoZ)
using detector $D_{1}$ [25], which occurs with a probability of success
$P_{s}$ that decreases as $R$ is increased. Upon successful operation, the
properties of the heralded output state $\hat{\rho}_{out}$ can then be
measured with an auxiliary detection system represented by $D_{2}$.

It can be shown that this conditional measurement process implements the basic
transformation $\ket{n}\rightarrow t^{n}\ket{n}$ in the Fock state basis
(here, transmittivity $|t|^{2}\equiv T=1-R$), which results in an overall
attenuation (i.e., $\langle\hat{n}\rangle_{out}<\langle\hat{n}\rangle_{in}$)
for all but Fock state inputs [26]. Roughly speaking, this attenuation effect
can be understood by considering that larger-$n$ terms in the Fock state
expansion of a general input state are less likely to “survive” the HoZ
process, resulting in a conditional re-weighting of the Fock state
coefficients towards smaller-$n$ terms in the heralded output state. Here we
extend this idea to show how ZPS can also transform the statistical properties
of input states, and then use these transformations to define a set of
nonclassicality criteria based on ZPS measurements.

As a simple example of the ability of ZPS to perform statistical
transformations, consider an input superposition state
$\ket{\psi}_{in}=\frac{1}{\sqrt{2}}(\ket{1}+\ket{5})$, which has a mean photon
number $\langle\hat{n}\rangle_{in}=3$ and a positive Mandel $Q$ parameter
$Q_{in}\approx 0.33$, denoting super-Poissonian photon statistics. After ZPS
with a standard 50:50 beamsplitter $(R=\frac{1}{2})$, the output state becomes
$\ket{\psi}_{out}=\frac{4}{\sqrt{17}}(\ket{1}+\frac{1}{4}\ket{5})$, which has
an attenuated mean photon number $\langle\hat{n}\rangle_{out}\approx 1.2$ and
a negative Mandel $Q$ parameter $Q_{out}\approx-0.28$. This negative value
denotes sub-Poissonian photon statistics and thus non-classical properties of
the heralded transmitted light [3].

Combined with a well-known “no-go” theorem which states that conditional
photon-counting measurements at a beamsplitter cannot produce nonclassical
output states unless the input is also nonclassical [27, 28], this simple
example shows how even a super-Poissonian input state may possess other
“hidden” non-classical properties that could be revealed by ZPS measurements.
In this paper, we formalize this argument by characterizing the ZPS
measurement process through a relative attenuation parameter $K(R)$ and
deriving nonclassicality criteria based on its behavior. While $K(R)$
typically changes monotonically as beamsplitter reflectance $R$ is increased,
we find that non-monotonic behavior arises for “transformable” states that
change their sub- or super-Poissonian character, and that this
“transformability” is always nonclassical.

The remainder of the paper is structured as follows: In Section II we define
the ZPS relative attenuation parameter $K(R)$ and explore its relationship
with $Q$. We illustrate the key features of $K(R)$ by comparing its behavior
for three basic (non-transformable) states with two simple examples of
transformable states. In Section III, we extend the concept of
“transformability” to define ZPS-based nonclassicality criteria, and in
Section IV we describe a method to predict which input states will
statistically transform under ZPS. In Section V we consider two more rich
examples of ZPS input states– the displaced squeezed state [29] and the
catalyzed coherent state [30]– which illustrate the concept of
transformability over a limited range of input state parameter space. Finally,
in Section VI we consider the effects of realistic detectors on experimentally
observing these ZPS-based statistical transformations.

##  II Relative attenuation parameter

When exactly zero photons are reflected by the beamsplitter in the setup of
Fig. 1, ZPS is successful and noiselessly attenuates the input state
$\hat{\rho}_{in}$ [18]. This process is described by the nonunitary operator
$t^{\hat{n}}$, producing a modified state $\hat{\rho}_{out}$ with resulting
mean photon number:

$\langle\hat{n}\rangle_{out}=\frac{\sum_{n}p_{n}nT^{n}}{\sum_{n}p_{n}T^{n}},$
(1)

where $\\{p_{n}\\}$ are the diagonal elements of $\hat{\rho}_{in}$ in the Fock
basis. With the exception of pure Fock state inputs, this always results in
$\langle\hat{n}\rangle_{out}<\langle\hat{n}\rangle_{in}$ [26]. This can be
seen by differentiating Eq. 1 with respect to $T$:

$\frac{d\langle\hat{n}\rangle_{out}}{dT}=\frac{1}{T}\left[\frac{\sum_{n}(n-\langle\hat{n}\rangle_{out})^{2}p_{n}T^{n}}{\sum_{n}p_{n}T^{n}}\right]=\frac{\langle(\Delta
n)^{2}\rangle_{out}}{T}\geq 0,$ (2)

and noting that as $T$ decreases ($R$ increases),
$\langle\hat{n}\rangle_{out}$ is nonincreasing and only remains constant in
the case of zero photon number variance $\langle(\Delta
n)^{2}\rangle_{out}=0$.

We quantify this photon number reduction through the relative attenuation
parameter $K(R)$, which is defined as the ratio of the mean output photon
number with HoZ (ZPS attenuation) to the mean output photon number without HoZ
(“ordinary” beamsplitter attenuation) [24]:

$K(R)\equiv\frac{\langle\hat{n}\rangle_{out}}{(1-R)\langle\hat{n}\rangle_{in}}$
(3)

This particular definition facilitates the use of coherent states (Poissonian
statistics) as a benchmark in experimental measurements. Because a coherent
state impinging on the beamsplitter will produce two completely uncorrelated
coherent states in the output ports [28], HoZ in the reflected port makes no
difference on $\langle\hat{n}\rangle_{out}$ in the transmitted port.
Consequently, $K(R)=1$ for all $R$ for any coherent state input
$\ket{\alpha}$.

In contrast, for any pure Fock state input $\ket{n}$ (sub-Poissonian
statistics), ZPS will not reduce $\langle\hat{n}\rangle_{out}$ while
“ordinary” beamsplitter attenuation grows with $R$, causing $K(R)$ to increase
monotonically. Conversely, for any thermal state input $\hat{\rho}_{th}$
(super-Poissonian statistics), ZPS will reduce $\langle\hat{n}\rangle_{out}$
more than “ordinary” beamsplitter attenuation [20, 24], causing $K(R)$ to
decrease monotonically with $R$.

These trends in the behavior of $K(R)$ for different types of input state
statistics can be formally seen by noting from Eq. 2 that the effects of ZPS
are highly dependent on photon number variance, and by extension Mandel’s $Q$
parameter, defined $Q=\frac{\langle(\Delta
n)^{2}\rangle}{\langle\hat{n}\rangle}-1$ [4]. Combining the results of Eqs.
1-3, one can obtain expressions for $Q$ in terms of the measurable quantity
$K$:

$Q_{out}(R^{\prime})=-\frac{(1-R^{\prime})}{K(R^{\prime})}\frac{dK}{dR}\Big{|}_{R=R^{\prime}}$
(4)

for the output state, and for the input state $(R=0)$:

$Q_{in}=-\frac{dK}{dR}\Big{|}_{R=0}$ (5)

In both cases, the slope of $K(R)$ is determined by the sign of $Q$, with
negative slopes denoting super-Poissonian statistics ($Q>0$) and positive
slopes denoting sub-Poissonian statistics ($Q<0$). Consequently, the initial
slope of $K(R)$ near $R=0$ determines whether the input state was super- or
sub-Poissonian [24], while any changes in the sign of the slope as $R$ is
increased denote the statistical transformations of interest here.

Figure 2: Plots of the ZPS relative attenuation parameter $K(R)$ for five
different input states. Non-monotonic behavior is a signature of a
nonclassical input state, with local extrema denoting transformations between
super- and sub-Poissonian photon statistics in the output state. Plots (a) -
(c) represent standard baseline examples of “non-transformable” input states,
while plots (d) and (e) represent simple examples of states that are
“transformable” by the ZPS process. For comparative purposes,
$\langle\hat{n}\rangle_{in}=3$ in examples (a) - (c).

Figure 2 illustrates these ideas by showing plots of relative attenuation
$K(R)$ for five different input states. By definition, each state begins with
$K=1$ at $R=0$ ($\langle\hat{n}\rangle_{out}=\langle\hat{n}\rangle_{in}$), and
the plots diverge from this common point as reflectance increases. Plots (a) -
(c) show the baseline examples of the three non-transformable states mentioned
above ($|\alpha\rangle$, $|n\rangle$, and $\hat{\rho}_{th}$) which show
monotonic behavior with initial slopes determined by Eq. 5. Plot (d)
corresponds to the transformable state described in Section I,
$|\psi\rangle_{in}=\frac{1}{\sqrt{2}}(|1\rangle+|5\rangle)$, which shows non-
monotonic behavior with a local minimum near $R\sim 0.4$. Past this point, we
see the state becomes sub-Poissonian according to Eq. 4.

The super- to sub-Poissonian statistical transformation revealed by this non-
monotonic behavior in $K(R)$ can be loosely understood in the following way:
when $R=0$, the equal weighting of $|1\rangle$ and $|5\rangle$ in
$|\psi\rangle_{in}$ yields a photon number variance larger than that of the
coherent state benchmark (i.e., super-Poissonian), resulting in a negative
initial slope for $K(R)$. As $R$ is increased, ZPS increasingly drives the
equally weighted superposition towards the Fock state $|1\rangle$, causing a
reduction in the photon number variance below that of the coherent state
benchmark (i.e., sub-Poissonian) and thus a transition to a positive slope for
$K(R)$.

Finally, plot (e) in Figure 2 displays a reverse transformation from sub- to
super-Poissonian statistics. Here we use the input state
$|\psi^{\prime}\rangle_{in}=\frac{1}{\sqrt{10}}|1\rangle+\frac{3}{\sqrt{10}}|5\rangle$,
which is heavily weighted towards $|5\rangle$ and thus starts with a small
(sub-Poissonian) photon-number variance and a positive initial slope for
$K(R)$. As $R$ is increased, ZPS drives the state towards a more equally
weighted superposition of $|1\rangle$ and $|5\rangle$, which has a larger
(super-Poissonian) photon-number variance, and causes the transition to a
negative slope for $K(R)$. As $R$ is further increased (past $R\sim 0.7$), we
replicate the behavior seen in plot (d): ZPS drives the superposition closer
to the Fock state $|1\rangle$ and $K(R)$ once again has a positive slope.
Overall, this state undergoes two statistical transformations– from sub- to
super-Poissonian statistics and back again– over the full range of $R$.

##  III Nonclassicality criteria

It is well known that conditional measurements at a beamsplitter (such as ZPS)
cannot produce a nonclassical output unless the input is also nonclassical
[27, 28]. Here, nonclassical states are defined in the usual way as those
which cannot be expressed as a mixture of coherent states [31]. The Glauber-
Sudarshan $P$-representation for such states cannot be interpreted as a valid
probability distribution, becoming highly singular or taking on negative
values [32, 33]. It is generally accepted that sub-Poissonian statistics are
nonclassical by this definition [3, 31]. It follows, therefore, that any input
state which becomes sub-Poissonian after ZPS must have been nonclassical to
begin with.

This useful restriction can be formulated in terms of the relative attenuation
parameter $K(R)$. Any classical (i.e., not nonclassical) input state
$\hat{\rho}_{in}$ is bound by the following inequality:

$\frac{dK}{dR}\leq 0\quad\text{for all {\hbox{R}}}$ (6)

Any state which is sub-Poissonian before or after ZPS will violate the above
inequality and certify $\hat{\rho}_{in}$ as nonclassical. Given that $K=1$ at
$R=0$ by definition, it immediately follows that for classical states:

$K(R)\leq 1\quad\text{for all {\hbox{R}}}$ (7)

As seen in Figure 2(a), the coherent state input $|\alpha\rangle$ saturates
these bounds with $K(R)=1$ and $dK/dR=0$ for all $R$.

To summarize, violations of Eqs. 6 and/or 7 represent nonclassicality criteria
that can be easily measured with the ZPS setup of Fig. 1, and are based on the
ability of ZPS to generate sub-Poissonian statistics in the output. This
transformation occurs if and only if the slope of $K(R)$ is positive for some
value of $R$. Alternatively, a single measurement of $K(R)>1$ for some $R$ is
sufficient (though not necessary) to identify such a nonclassical state. Of
course, sub-Poissonian input states are guaranteed to trigger both of these
nonclassicality criteria.

##  IV Predicting transformability

From the simple examples in Section II, it is clear that not all super-
Poissonian states can transform into sub-Poissonian states under ZPS, or vice
versa. We have shown that “transformable” states must be nonclassical in Sec.
III, but the nonclassicality criteria in Eqs. 6 and 7 provide little a priori
insight into which states will actually transform. In this section, we derive
sufficient conditions for predicting transformability based on the first few
terms of the photon number distribution.

The key insight is that if a state only transforms once over the full range of
$R$, we need only to consider the behavior of $K(R)$ at maximum attenuation to
verify if the statistics have changed. In the limit $R\rightarrow 1$, $K(R)$
and its derivative are given by:

$\displaystyle\lim_{R\to
1}K(R)=\frac{1}{\langle\hat{n}\rangle_{in}}\left(\frac{p_{1}}{p_{0}}\right)$
(8) $\displaystyle\lim_{R\to
1}\frac{dK}{dR}=\frac{1}{\langle\hat{n}\rangle_{in}}\left(\frac{p_{1}^{2}-2p_{0}p_{2}}{p_{0}^{2}}\right)$
(9)

Remarkably, these values are entirely determined by the mean photon number
$\langle\hat{n}\rangle_{in}$ and first three photon number probabilities
$(p_{0},p_{1},p_{2})$ of the input state.

For a super-Poissonian input state, transformability amounts to violating the
classical bounds in Eqs. 6 and 7, meaning $dK/dR>0$ or $K>1$ for some value of
$R$. Combining this with Eqs. 8 and 9, we obtain the following
transformability criteria:

$\displaystyle p_{0}$ $\displaystyle=0$ (10) $\displaystyle
p_{0}\langle\hat{n}\rangle_{in}$ $\displaystyle<p_{1}$ (11) $\displaystyle
2p_{0}p_{2}$ $\displaystyle<p_{1}^{2}$ (12)

Satisfying any of the above is sufficient to show that a super-Poissonian
input state will transform for sufficiently large $R$.

We note that Eq. 10 is equivalent to Lee’s theorem, which states that any
input state with zero vacuum probability $p_{0}=0$ is nonclassical [34], and
is often invoked to understand the nonclassical nature of photon-added states
[10, 13]. In the case of ZPS, it is clear that if $p_{0}=0$, $K(R)$ and its
derivative will diverge to infinity as $R\rightarrow 1$, triggering both
nonclassicality criteria in Eqs. 6 and 7. In addition, Eq. 12 is equivalent to
Klyshko’s well-known nonclassicality criterion [35]. Thus, all super-
Poissonian states which trigger Lee’s or Klyshko’s nonclassicality criteria
are transformable.

Conversely, a sub-Poissonian input state will transform if $dK/dR<0$ or $K<1$
for some value of $R$. It follows that satisfying either of Eqs. 11 or 12,
with the inequalities reversed, is sufficient to show that a sub-Poissonian
state is transformable. Lee’s criteria (Eq. 10) cannot be used to predict
transformability for sub-Poissonian states.

The above methods for predicting transformability rely only on comparing
statistics at $R=0$ and $R\rightarrow 1$. Thus, states which transform an even
number of times over the full range of $R$ may not trigger these criteria, as
in the example of Figure 2(e) which is sub-Poissonian at both extremes. This
also holds for super-Poissonian input states, so satisfying Klyshko’s
inequality (Eq. 12) is sufficient but not necessary to be transformable. For
example, any mixture/superposition of 0, 2 and 6 photons does not satisfy
Klyshko’s inequality, yet it can be shown that the specific case of
$p_{0}=0.04$, $p_{2}=0.48$, $p_{6}=0.48$ corresponds to a transformable super-
Poissonian state.

##  V Input state parameters

For many types of quantum optical states, the initial statistical character is
determined by a set of input parameters. A trivial example is the two-term
Fock state superposition considered in Figs. 2(d) and (e),
$|\psi\rangle_{in}=\gamma|1\rangle+\sqrt{(1-\gamma^{2})}|5\rangle$, which was
initially super- or sub-Poissonian based on the value of $\gamma$. More
complex examples include squeezed and/or displaced states, in which the
magnitude and angle of squeezing and displacement can result in dramatically
different photon statistics [17, 15], and conditional state preparation
techniques like photon catalysis [16].

It is interesting to extend this idea to address the following question: for a
given type of input state, what regions of its input parameter space will
enable it to transform under ZPS? Here we consider two rich examples within
this context– displaced squeezed states, as considered by Dodonov et al. [29],
and catalyzed coherent states [30].

###  V.1 Displaced squeezed state

Noiseless attenuation of general gaussian states was previously considered in
Ref. [26], and it was shown that ZPS preserves gaussianity. However, we show
here that ZPS may not preserve the sub- or super-Poissonian character of these
states for a significant portion of their parameter space. We consider a
displaced squeezed state with real squeezing parameter $r$ and displacement
magnitude $z$. To restrict our space to two variables, the angles of squeezing
and displacement are chosen so the two quadrature variances are equal
$\sigma_{x}=\sigma_{p}$, and their means opposite
$\langle\hat{x}\rangle=-\langle\hat{p}\rangle$, following Dodonov et al. [29].

Starting with the photon statistics derived in [29], the relative attenuation
$K(R)$ can be shown to be:

$\begin{split}K(R)=\frac{1}{\sinh^{2}{r}+z^{2}}\left[\frac{1}{2\left(\coth{r}-(1-R)\right)}\right.\\\
-\left.\frac{1}{\coth{r}+\left(1-R\right)}\left(\frac{1}{2}-2z^{2}\frac{\coth{2r}+1}{1+(1-R)\tanh{r}}\right)\right]\end{split}$
(13)

Figure 3 shows plots of $K(R)$ for multiple combinations of displacement and
squeezing parameters $z$ and $r$. For $z=1$ and $r=1$ in Fig. 3(c), a local
minimum appears. This indicates that after sufficient noiseless attenuation,
the super-Poissonian input state undergoes a transition to sub-Poissonian
statistics in the output.

Figure 3: Relative attenuation $K(R)$ for a displaced squeezed state with
parameters (a) $z=1$, $r=0.05$, (b) $z=1$, $r=0.35$, (c) $z=1$, $r=1$, and (d)
$z=0.1$, $r=0.1$. Cases (a), (b) and (d) do not exhibit statistical
transformations, while case (c) illustrates a transformation from super- to
sub-Poissonian statistics under ZPS.

There is a wide portion of parameter space for which this transition is
possible, as shown by the blue shaded region in Fig. 4(a). The bounding curves
are calculated analytically in Mathematica. The bottom curve corresponds to
the bound of Klyshko’s inequality (Eq. 12), while the top curve divides sub-
from super-Poissonian input states. This leaves three total regions: (i) the
blank bottom-most region which contains super-Poissonian states with low
squeezing, including squeezed vacuum states; (ii) the transformable super-
Poissonian states, with moderate displacement and squeezing; and (iii) the
upper-most region of sub-Poissonian states, including those with the highest
squeezing. There are no cases in which ZPS transforms a sub-Poissonian state
of this kind into a super-Poissonian one.

Figure 4: Regions shaded blue and red correspond to parameter spaces for which
the input of (a) a displaced squeezed state and (b) a catalyzed coherent
state, has a local minimum and maximum respectively in $K(R)$. For (b), these
regions overlap. An example from within the overlapping region is shown in
Fig. 5(c).

As discussed in Sec. III, the transition in Fig. 3(c) is clearly nonclassical
according to Eq. 6, and likewise satisfies Klyshko’s inequality (Eq. 12).
However, it is interesting to note that the input has no “hidden” higher-order
sub-Poissonian statistics as defined in Ref. [36]. In other words, it can be
shown that the normalized correlations functions for this state satisfy
$g^{(n)}(0)>1$ for all $n$. In this sense, ZPS does not “reveal” sub-
Poissonian statistics, but rather generates them from other nonclassical
aspects of the initial number distribution.

###  V.2 Catalyzed coherent state

Catalyzed coherent states (CCS) can be generated by mixing a coherent state
$\ket{\alpha}$ with a single-photon Fock state at a beamsplitter with
reflectance $\Lambda$, and conditioning the output on the detection of exactly
one photon in the auxiliary output mode [30]. As first proposed by Xu and
Yuan, a similar catalysis procedure can be implemented with the signal and
idler modes of an optical parametric amplifier (OPA) with an equivalent
catalysis parameter $\Lambda=1-1/g^{2}$, where $g$ is the gain of the OPA [37,
38, 39].

Here, we consider an input CCS with coherent state amplitude $\alpha$ and
catalysis parameter $\Lambda$. It can be shown that the photon number
distribution for this state is:

Figure 5: Relative attenuation $K(R)$ for a catalyzed coherent state with
parameters (a) $\Lambda=0.3$, $\alpha=1$, (b) $\Lambda=0.75$, $\alpha=1$, (c)
$\Lambda=0.2$, $\alpha=4$, and (d) $\Lambda=0.9$, $\alpha=1$. Cases (a) and
(b) exhibit a super-to-sub- and sub-to-super-Poissonian transformation,
respectively, while case (c) shows both and (d) shows neither.

$p_{n}=\frac{e^{-|\alpha|^{2}(1-\Lambda)}|\alpha|^{2n}(1-\Lambda)^{n-1}(\Lambda+\Lambda
n-1)^{2}}{n!\left[1-\Lambda-|\alpha|^{2}\Lambda(\Lambda^{2}-4\Lambda+2)\right]}$
(14)

Using Eqs. 1 and 3, $K(R)$ is then given by:

$K(R)\propto\frac{\displaystyle\begin{multlined}1-4\Lambda(1-\Lambda)-\alpha^{2}\Lambda(2-5\Lambda)(1-R)\\\
+\alpha^{4}\Lambda^{2}(1-\Lambda)^{2}(1-R)^{2}\end{multlined}1-4\Lambda(1-\Lambda)-\alpha^{2}\Lambda(2-5\Lambda)(1-R)\\\
+\alpha^{4}\Lambda^{2}(1-\Lambda)^{2}(1-R)^{2}}{\displaystyle\begin{multlined}1-\Lambda-\alpha^{2}\Lambda(2-3\Lambda)(1-R)\\\
+\alpha^{4}\Lambda^{2}(1-\Lambda)(1-R)^{2}\end{multlined}1-\Lambda-\alpha^{2}\Lambda(2-3\Lambda)(1-R)\\\
+\alpha^{4}\Lambda^{2}(1-\Lambda)(1-R)^{2}}$ (15)

with the necessary normalization such that $K=1$ at $R=0$.

Like the displaced squeezed state, relative attenuation of the CCS can exhibit
a local minimum for particular values of the input parameters. In Fig. 5(a), a
minimum is found for $\Lambda=0.3$ and $\alpha=1$. Just like the example in
Fig. 3(c), this state is super-Poissonian to arbitrarily high order
($g^{(n)}(0)>1$) and satisfies Klyshko’s nonclassicality criteria. Unlike the
previous example, however, $K(R)$ can also have a local maximum for different
parameters, as in Fig. 5(b). In this example, all higher-order statistics of
the input state are sub-Poissonian ($g^{(n)}(0)<1$), and yet attenuation by
ZPS produces super-Poissonian output statistics for $R\gtrsim 0.52$.

The parameter space for local minima and maxima is shown in Fig. 4(b), again
calculated analytically with Mathematica. Interestingly, the two regions
overlap, allowing for input cases with both a minimum and maximum, as in Fig.
5(c). Additionally, Klyshko’s inequality is satisfied when $\Lambda\lesssim
0.365$ and when $\Lambda\gtrsim 0.775$. The former region contains all super-
to-sub transformations, while the latter contains no transforming states.

##  VI Realistic detectors

In a real experiment, the ability to observe the statistical transformations
of interest will critically depend on the performance of detectors $D_{1}$ and
$D_{2}$ in Fig. 1. Here we consider the effects of dark counts and
inefficiency in both detectors, as well as the lack of photon number
resolution (PNR) for $D_{2}$ in the transmitted output port of the
beamsplitter.

If heralding detector $D_{1}$ has a reduced effective efficiency $\eta_{1}$
(including losses in the reflected mode or before the beamsplitter) then it
can be shown that [24]:

$\displaystyle K_{\text{exp}}(R)=K_{\text{ideal}}(R\eta_{1})$ (16)

If $\eta_{1}$ is known, then $K_{\text{ideal}}(R)$ can be recovered over the
smaller domain $R\in[0,\eta_{1})$. However, if any statistical transformations
or other interesting behavior only exist beyond $R=\eta_{1}$, they can no
longer be generated in the transmitted output port. In the limit of
$\eta_{1}\rightarrow 0$, ZPS is equivalent to ordinary attenuation, which can
never change the sign of Mandel’s $Q$ parameter [40].

Meanwhile, $D_{2}$ efficiency $\eta_{2}$ or loss in the output mode has no
effect on $K(R)$ as defined in Eq. 3. The numerator and denominator of this
ratio would be reduced by the same efficiency factor $\eta_{2}$, which then
cancels. However, the efficiency $\eta_{2}$ must be carefully considered if
using a non-PNR detector $D_{2}$ (see Fig. 6). In this case, mean intensity or
expected photon number at $D_{2}$ is replaced with the probability of a click
event:

$\displaystyle K_{\text{click}}(R)$
$\displaystyle\equiv\frac{P(C_{2}|NC_{1})}{P(C_{2})}$
$\displaystyle=\frac{\text{tr}\left\\{\hat{B}\hat{\rho}_{in}\hat{B}^{\dagger}\hat{\Pi}^{(NC)}_{1}\hat{\Pi}^{(C)}_{2}\right\\}}{\text{tr}\left\\{\hat{B}\hat{\rho}_{in}\hat{B}^{\dagger}\hat{\Pi}^{(C)}_{2}\right\\}\text{tr}\left\\{\hat{B}\hat{\rho}_{in}\hat{B}^{\dagger}\hat{\Pi}^{(NC)}_{1}\right\\}}$
$\displaystyle=\frac{\sum_{n}p_{n}\left[(1-R\eta_{1})^{n}-(1-R\eta_{1}-T\eta_{2})^{n}\right]}{\left[\sum_{n}p_{n}(1-R\eta_{1})^{n}\right]\left[1-\sum_{n}p_{n}(1-T\eta_{2})^{n}\right]}$
(17)

where the subscripts $1$ and $2$ indicate detection channels $D_{1}$ and
$D_{2}$, $\hat{B}$ is the unitary beam-splitter operator [3], and we have used
the standard POVMs for non-PNR detectors for “click” $(C)$ and “no-click”
$(NC)$ events [41]:

$\displaystyle\hat{\Pi}^{(C)}_{i}$ $\displaystyle=\mathbbm{1}-\Pi^{(NC)}_{i}$
$\displaystyle\hat{\Pi}^{(NC)}_{i}$
$\displaystyle=\sum_{n}(1-\eta_{i})^{n}\ket{n}\bra{n}$ (18)

for $i=1,2$. We find that $\lim_{\eta_{2}\to
0}K_{\text{click}}(\eta_{1},\eta_{2},R)=K(\eta_{1}R)$ after applying
L’Hôpital’s rule. Alternatively, the binomial approximation can be applied so
long as $|\eta_{2}n|\ll 1$ and higher-number terms are insignificant. This is
similar to the requirements for measuring normalized correlation functions
$g^{(n)}$. In these multiphoton coincidence measurements, non-PNR detectors
are effective for low efficiency and low counting rates, when the detector
response to photon number is approximately linear [42]. These measurements are
considered loss-tolerant exactly like the transmitted mode in ZPS.

Figure 6: ZPS measurement with two non-PNR detectors. Detector $D_{2}$
measures the rate of clicks with and without conditioning on a no-click at
$D_{1}$. Detector inefficiencies $\eta_{1}$ and $\eta_{2}$ are modeled with
beamsplitters that reflect photons into the environment. To approximate $K(R)$
without PNR, the effective efficiency of $D_{2}$ must satisfy $|\eta_{2}n|\ll
1$.

Dark counts at the heralding detector $D_{1}$ have no effect on $K(R)$, only
reducing the probability of success for HoZ [25]. These false click events can
be safely ignored so long as $D_{1}$ dark counts are completely uncorrelated
with the counts at the other detector $D_{2}$.

In contrast, dark counts at $D_{2}$ do have an effect on the measured relative
attenuation. For some dark count probability $d$ at $D_{2}$, Eq. VI is
modified as follows:

$K_{\text{dark}}(R)=\frac{P(C_{2}|NC_{1})+d\cdot
P(NC_{2}|NC_{1})}{P(C_{2})+d\cdot P(NC_{2})}$ (19)

If we assume the probability of a simultaneous dark and “true” count is
negligible– i.e., counting rates are sufficiently low with and without HoZ,
then this becomes:

$K_{\text{dark}}(R)\approx\frac{P(C_{2}|NC_{1})+d}{P(C_{2})+d}$ (20)

If the dark count rate $d$ is simply subtracted from both the numerator and
denominator, we regain the original ratio $K_{\text{click}}$. In any case,
lower $D_{2}$ dark counts are desirable to improve the signal-to-noise ratio.

In summary, the function $K(R)$ can be reliably measured with currently
available single-photon detectors, even without PNR capability. The heralding
detector $D_{1}$ requires high efficiency to generate the full range of
possible output statistics, but is robust to dark counts and other
uncorrelated background noise. The output detector $D_{2}$ has essentially
opposite requirements, requiring low effective efficiency (if non-PNR) and
relatively low dark counts.

##  VII Summary & Conclusions

We have shown how the conditional measurement process of zero-photon
subtraction (ZPS) can transform certain super-Poissonian states into sub-
Poissonian states, and vice versa. These effects can be experimentally
observed by measuring the relative attenuation parameter $K(R)$, which is
simply defined as the ratio of ZPS-based attenuation to ordinary beamsplitter
attenuation. Because the slope of $K(R)$ is proportional to Mandel’s
$Q$-parameter, an observed local minimum or maximum as the beamsplitter
reflectivity is tuned from $R=0\rightarrow 1$ is the signature of a
statistical transformation of interest.

We described how input states which transform in this way are necessarily
nonclassical, which allowed us to establish nonclassicality criteria based on
ZPS measurements. The connection between these criteria and Klyshko and Lee’s
theorems [35, 34] showed that certain restrictions on photon number
probabilities provide sufficient, but not necessary, conditions to predict a
priori which input states will transform through the ZPS process.

We considered several simple examples of input states that illustrated the
basic physics of these ZPS-based statistical transformations, as well as two
more complex example states that showed interesting parameter-dependent
behavior. In all cases, the effects of realistic detector parameters (low
efficiency and dark counts) were found to only moderately reduce the ability
to observe the statistical transformations of interest. Consequently, an
experimental demonstration could be feasible with currently available detector
technology [43], although high-fidelity preparation of the various quantum
optical input states considered here would remain a challenge.

## Acknowledgements

We would like to thank J. D. Franson for many valuable discussions pertaining
to this research. This work was supported by the National Science Foundation
under Grant No. 2013464.

## References

  * Hanbury Brown and Twiss [1956] R. Hanbury Brown and R. Q. Twiss, Nature 178, 1046 (1956).
  * Short and Mandel [1983] R. Short and L. Mandel, Phys. Rev. Lett. 51, 384 (1983).
  * Agarwal [2012] G. S. Agarwal,  _Quantum Optics_ (Cambridge University Press, 2012).
  * Mandel [1979] L. Mandel, Opt. Lett. 4, 205 (1979).
  * Gilles and Knight [1993] L. Gilles and P. L. Knight, Phys. Rev. A 48, 1582 (1993).
  * Birnbaum  _et al._ [2005] K. M. Birnbaum, A. Boca, R. Miller, A. D. Boozer, T. E. Northup, and H. J. Kimble, Nature 436, 87 (2005).
  * Hong  _et al._ [1985] C. K. Hong, S. Friberg, and L. Mandel, J. Opt. Soc. Am. B 2, 494 (1985).
  * Boyd  _et al._ [2008] R. W. Boyd, G. S. Agarwal, K. W. Clifford Chan, A. K. Jha, and M. N. O’Sullivan, Opt. Commun. 281, 3732 (2008).
  * Li  _et al._ [2020] S.-W. Li, F. Li, T. Peng, and G. S. Agarwal, Phys. Rev. A 101, 063806 (2020).
  * Agarwal and Tara [1991] G. S. Agarwal and K. Tara, Phys. Rev. A 43, 492 (1991).
  * Ourjoumtsev  _et al._ [2006] A. Ourjoumtsev, R. Tualle-Brouri, J. Laurat, and P. Grangier, Science 312, 83 (2006).
  * Zavatta  _et al._ [2008] A. Zavatta, V. Parigi, M. S. Kim, and M. Bellini, New J. Phys. 10, 123006 (2008).
  * Barnett  _et al._ [2018] S. M. Barnett, G. Ferenczi, C. R. Gilson, and F. C. Speirits, Phys. Rev. A 98, 013809 (2018).
  * Arnoldus and Nienhuis [1983] H. Arnoldus and G. Nienhuis, Opt. Acta 30, 1573 (1983).
  * Grosse  _et al._ [2007] N. B. Grosse, T. Symul, M. Stobińska, T. C. Ralph, and P. K. Lam, Phys. Rev. Lett. 98, 153603 (2007).
  * Bartley  _et al._ [2012] T. J. Bartley, G. Donati, J. B. Spring, X. M. Jin, M. Barbieri, A. Datta, B. J. Smith, and I. A. Walmsley, Phys. Rev. A 86, 043820 (2012).
  * de Oliveira  _et al._ [1990] F. A. M. de Oliveira, M. S. Kim, P. L. Knight, and V. Buek, Phys. Rev. A 41, 2645 (1990).
  * Mičuda  _et al._ [2012] M. Mičuda, I. Straka, M. Miková, M. Dušek, N. J. Cerf, J. Fiurášek, and M. Ježek, Phys. Rev. Lett. 109, 180503 (2012).
  * Allevi  _et al._ [2010] A. Allevi, A. Andreoni, M. Bondani, M. G. Genoni, and S. Olivares, Phys. Rev. A 82, 013816 (2010).
  * Zhai  _et al._ [2013] Y. Zhai, F. Becerra-Chavez, B. Glebov, J. Fan, S. W. Nam, and A. Migdall, Opt. Lett. 38, 2171 (2013).
  * Bogdanov  _et al._ [2017] Y. I. Bogdanov, K. G. Katamadze, G. V. Avosopiants, L. V. Belinsky, N. A. Bogdanova, A. A. Kalinkin, and S. P. Kulik, Phys. Rev. A 96, 063803 (2017).
  * Magaña-Loaiza  _et al._ [2019] O. S. Magaña-Loaiza, R. de J. León-Montiel, A. Perez-Leija, A. B. U’Ren, C. You, K. Busch, A. E. Lita, S. W. Nam, R. P. Mirin, and T. Gerrits, npj Quantum Inf. 5, 1 (2019).
  * Katamadze  _et al._ [2020] K. G. Katamadze, G. V. Avosopiants, N. A. Bogdanova, Y. I. Bogdanov, and S. P. Kulik, Phys. Rev. A 101, 013811 (2020).
  * Nunn  _et al._ [2022] C. M. Nunn, J. D. Franson, and T. B. Pittman, Phys. Rev. A 105, 033702 (2022).
  * Nunn  _et al._ [2021] C. M. Nunn, J. D. Franson, and T. B. Pittman, Phys. Rev. A 104, 033717 (2021).
  * Gagatsos  _et al._ [2014] C. N. Gagatsos, J. Fiurášek, A. Zavatta, M. Bellini, and N. J. Cerf, Phys. Rev. A 89, 062311 (2014).
  * Ban [1996] M. Ban, J. Mod. Opt. 43, 1281 (1996).
  * Kim  _et al._ [2002] M. S. Kim, W. Son, V. Buzžek, and P. L. Knight, Phys. Rev. A 65, 032323 (2002).
  * Dodonov  _et al._ [1994] V. V. Dodonov, I. M. Dremin, P. G. Polynkin, and V. I. Man’ko, Phys. Lett. A 193, 209 (1994).
  * Lvovsky and Mlynek [2002] A. I. Lvovsky and J. Mlynek, Phys. Rev. Lett. 88, 250401 (2002).
  * Dodonov [2002] V. V. Dodonov, Journal of Optics B 4, R1 (2002).
  * Glauber [1963] R. J. Glauber, Physical Review 130, 2529 (1963).
  * Sudarshan [1963] E. C. Sudarshan, Phys. Rev. Lett. 10, 277 (1963).
  * Lee [1995] C. T. Lee, Phys. Rev. A 52, 3374 (1995).
  * Klyshko [1996] D. Klyshko, Phys. Lett. A 213, 7 (1996).
  * Erenso  _et al._ [2002] D. Erenso, R. Vyas, and S. Singh, J. Opt. Soc. Am. B 19, 1471 (2002).
  * Xu and Yuan [2016] X. X. Xu and H. C. Yuan, Phys. Lett. A 380, 2342 (2016).
  * Shringarpure and Franson [2019] S. U. Shringarpure and J. D. Franson, Phys. Rev. A 100, 043802 (2019).
  * Barnett  _et al._ [1999] S. M. Barnett, D. T. Pegg, and J. Jeffers, Opt. Commun. 172, 55 (1999).
  * Alléaume  _et al._ [2004] R. Alléaume, F. Treussart, J. M. Courty, and J. F. Roch, New J. Phys. 6, 85 (2004).
  * Stevens [2013] M. J. Stevens,  _Photon Statistics, Measurements, and Measurements Tools_, edited by A. Migdall, S. Polyakov, J. Fan, and J. Bienfang, Vol. 45 (Academic Press, 2013) pp. 25–68.
  * Avenhaus  _et al._ [2010] M. Avenhaus, K. Laiho, M. V. Chekhova, and C. Silberhorn, Phys. Rev. Lett. 104, 063602 (2010).
  * Eisaman  _et al._ [2011] M. D. Eisaman, J. Fan, A. Migdall, and S. V. Polyakov, Rev. Sci. Instrum. 82, 071101 (2011) .


* [Ion02] Eleny-Nicoleta Ionel  “Topological recursive relations in $H^{2g}(M_{g,n})$”  In _Invent. Math._ 148.3, 2002, pp. 627–658  DOI: 10.1007/s002220100205
  * [IP04] Eleny-Nicoleta Ionel and Thomas H. Parker  “The symplectic sum formula for Gromov-Witten invariants”  In _Ann. of Math. (2)_ 159.3, 2004, pp. 935–1025  DOI: 10.4007/annals.2004.159.935
  * [Kra18] Thomas Kragh  “The Viterbo transfer as a map of spectra”  In _J. Symplectic Geom._ 16.1, 2018, pp. 85–226  DOI: 10.4310/JSG.2018.v16.n1.a3
  * [Lip06] Robert Lipshitz  “A cylindrical reformulation of Heegaard Floer homology”  In _Geom. Topol._ 10, 2006, pp. 955–1096  DOI: 10.2140/gt.2006.10.955
  * [McD06] Dusa McDuff  “Groupoids, branched manifolds and multisections”  In _J. Symplectic Geom._ 4.3, 2006, pp. 259–315  URL: http://projecteuclid.org/euclid.jsg/1180135649
  * [McD19] Dusa McDuff  “Constructing the virtual fundamental class of a Kuranishi atlas”  In _Algebr. Geom. Topol._ 19.1, 2019, pp. 151–238  DOI: 10.2140/agt.2019.19.151
  * [Moc95] Shinichi Mochizuki  “The geometry of the compactification of the Hurwitz scheme”  In _Publ. Res. Inst. Math. Sci._ 31.3, 1995, pp. 355–441  DOI: 10.2977/prims/1195164048
  * [MS21] Cheuk Yu Mak and Ivan Smith  “Non-displaceable Lagrangian links in four-manifolds”  In _Geom. Funct. Anal._ 31.2, 2021, pp. 438–481  DOI: 10.1007/s00039-021-00562-8
  * [MW17] Dusa McDuff and Katrin Wehrheim  “Smooth Kuranishi atlases with isotropy”  In _Geom. Topol._ 21.5, 2017, pp. 2725–2809  DOI: 10.2140/gt.2017.21.2725
  * [Par16] John Pardon  “An algebraic approach to virtual fundamental cycles on moduli spaces of pseudo-holomorphic curves”  In _Geom. Topol._ 20.2, 2016, pp. 779–1034  DOI: 10.2140/gt.2016.20.779
  * [Prü17] Sven Prüfer  “Symplectic Geometry of Moduli Spaces of Hurwitz Covers”, 2017  arXiv:1711.07048
  * [Rit10] Alexander F. Ritter  “Deformations of symplectic cohomology and exact Lagrangians in ALE spaces”  In _Geom. Funct. Anal._ 20.3, 2010, pp. 779–816  DOI: 10.1007/s00039-010-0074-7
  * [Rit13] Alexander F. Ritter  “Topological quantum field theory structure on symplectic cohomology”  In _J. Topol._ 6.2, 2013, pp. 391–489  DOI: 10.1112/jtopol/jts038
  * [Sei08] Paul Seidel  “A biased view of symplectic cohomology”  In _Current developments in mathematics, 2006_ Int. Press, Somerville, MA, 2008, pp. 211–253
  * [Sei08a] Paul Seidel  “Fukaya categories and Picard-Lefschetz theory”, Zurich Lectures in Advanced Mathematics  European Mathematical Society (EMS), Zürich, 2008, pp. viii+326  DOI: 10.4171/063
  * [Sie97] Bernard Siebert  “Gromov-Witten invariants of general symplectic manifolds” preprint, 1997
  * [Vit95] Claude Viterbo  “Generating functions, symplectic geometry, and applications”  In _Proceedings of the International Congress of Mathematicians, Vol. 1, 2 (Zürich, 1994)_ Birkhauser, Basel, 1995, pp. 537–547
  * [Vit99] C. Viterbo  “Functors and computations in Floer homology with applications. I”  In _Geom. Funct. Anal._ 9.5, 1999, pp. 985–1033  DOI: 10.1007/s000390050106


# Preserving Patient Privacy in MRI Scans: A Comprehensive Approach with 3D
Masked Autoencoders

Lennart Alexander Van der Goten
KTH Royal Institute of Technology
SciLifeLab
Stockholm, SWEDEN
<EMAIL_ADDRESS>Kevin Smith
KTH Royal Institute of Technology
SciLifeLab
Stockholm, SWEDEN
<EMAIL_ADDRESS>
###### Abstract

MRI scans provide valuable medical information, however they also contain
sensitive and personally identifiable information (PII) that needs to be
protected. Whereas MRI metadata is easily sanitized, MRI image data is a
privacy risk because it contains information to render highly-realistic 3D
visualizations of a patient’s head, enabling malicious actors to possibly
identify the subject by cross-referencing a database. Data anonymization and
de-identification is concerned with ensuring the privacy and confidentiality
of individuals’ personal information. Traditional MRI de-identification
methods remove privacy-sensitive parts (_e.g_. eyes, nose etc.) from a given
scan. This comes at the expense of introducing a domain shift that can throw
off downstream analyses. Recently, a GAN-based approach was proposed to de-
identify a patient’s scan by remodeling it (_e.g_. changing the face) rather
than by removing parts. In this work, we propose CP-MAE, a model that de-
identifies the face using masked autoencoders and that outperforms all
previous approaches in terms of downstream task performance as well as de-
identification. With our method we are able to synthesize scans of resolution
up to $256^{3}$ (previously $128^{3}$) which constitutes an eight-fold
increase in the number of voxels. Using our construction we were able to
design a system that exhibits a highly robust training stage, making it easy
to fit the network on novel data.

---
Figure 1: MRI scans pose a privacy risk since highly-realistic face renderings
can be crafted and misused for malicious purposes. Our model aims to advance
the so-called _remodeling-based_ subclass of MRI de-identification which
retains the brain and remodels all other features. (_top_ : 3D view, _bottom_
slice-view)

##  1 Introduction

Magnetic resonance imaging (MRI) is a non-invasive, high-resolution imaging
technology that enables the analysis of anatomical structures such as brain,
spine, joints, abdomen, pelvis, the cardiovascular system and associated
diseases. While MRI scans are usually visualized as 2D slices over one of the
three axes, it is also possible to render a high-quality 3D model using
techniques such as volumetric raytracing, yielding a realistic depiction of a
patient’s face. This is problematic in terms of data privacy: Given a face
database with associated identites, a malicious actor can find the closest
match to a given face rendering, _allowing them to infer the patient’s
identity_.

Various de-identification methods have appeared over the years that address
the removal of privacy-sensitive parts111In the most extreme case just
retaining the brain and mostly differ in their level of aggressiveness: While
FACE MASK [24] merely blurs out the face, skull-stripping methods such as MRI
WATERSHED [35] only retain the brain and remove everything else. However,
these traditional approaches are potentially problematic. Simply put, many
tools in the medical workflow require the presence of certain landmarks as
preconditions. If those landmarks are absent, the analysis might be impaired
or, in the worst case, it might be impossible to perform the analysis at all.
On the other hand, if a method is not aggressive enough it might still be
feasible to infer a patient’s identity, a fact that proves especially harmful
when such a tool is to be used automatically and without oversight on a large
corpus of scans.

Consequently, _removal-based_ methods seek to find a sweet spot which ensures
acceptable levels of both privacy and downstream task compatibility. A
recently proposed GAN-based approach named CP-GAN [11] introduced a new
paradigm: to remodel the face and leave the medically-relevant information
intact. This _remodeling-based_ class of models aims to remodel the skull and
facial features (staying true to the overall data distribution of MRI scans)
while leaving the brain _untouched_. Put differently, a patient’s MRI scan is
de-identified by copying their brain and remodeling the remaining parts such
that they “look like” as if they belong to an actual MRI scan but do not give
away information about the patient’s real face (or identity). This approach
eliminates the need for the aforementioned trade-off, as it ensures accurate
placement of landmarks and thus minimizes domain shift, leading to more
effective de-identification performance.

Recently, the emergence of _diffusion-probabilistic_ models (DPMs) and _masked
autoencoders_ (MAEs) has propelled generative approaches forward, deviating
from the previously dominant _generative adversarial network_ (GAN)
methodology. Instead of relying on the definition of an adversarial game to
learn the underlying data distribution, these new model classes adopt
significantly different principles. Masked autoencoders draw inspiration from
natural language processing (NLP) models such as BERT [12], and focus on
predicting stochastically-masked segments of the data while conditioning on
the unmasked portions. Data points can be sampled by starting with randomly
initialized data and allowing the masked autoencoder to fill in the missing
parts based on its highest confidence. MAEs are generally considered to be
more efficient [30, 41] than DPMs during inference as significantly fewer
steps are needed to synthesize a high-fidelity sample.

MAEs also offer two fundamental advantages over GANs: (i) Higher training
stability and (ii) typically a lower memory footprint222Due to the absence of
a discriminator at training time. Considering that MRI de-identification with
generative models involves high memory requirements and often suffers from
instability due to small batch sizes, MAEs seem to be a promising alternative.

In this work, we adapt masked autoencoders for MRI de-identfication. Our
contributions are as follows:

  1. 1.

We suggest a novel _masked autoencoder_ -based model named CP-MAE

  2. 2.

We show that CP-MAE can deal with higher-resolution MRI scans than previous
remodeling-based methods, improving the state-of-the-art from a resolution of
$128^{3}$ to $256^{3}$, effectively octupling the number of synthesized voxels

  3. 3.

To the best of our knowledge, our model is the first instantiation of
combining a volumetric VQ-VAE with an MAE for 3D MR image synthesis

  4. 4.

We demonstrate that CP-MAE features superior de-identification performance
compared to other methods and can be robustly trained on modestly-sized
datasets (_e.g_. ADNI, OASIS-3)

  5. 5.

We show that de-identification with CP-MAE introduces minimal effects on brain
tissue and subcortical segmentation tasks compared to other approaches.

Figure 2: Preprocessing/Vector Quantization Stage. We first apply standard
preprocessing algorithms to a given MRI scan to make the data easier to model.
When then execute two independent stages: (i) We use ROBEX [20] to extract the
(continuous) brain of the scan, and inversely (ii) we remove the brain from
the full skull. Both representations are then _compressed_ independently by
two VQ-VAEs into 3D integer volumes of much lower resolution. To help
generalization we apply various augmentation techniques. After the training
phase, a “paired dataset” is created that contains the fixed discrete codes
derived from both representations.

##  2 Related Work

MRI De-Identification. MRI de-identification is a critical pre-processing step
in neuroimaging, and many methods have been proposed to effectively perform
this task. Smith _et al_. introduced the Brain Extraction Tool (BET) [36],
which is widely used for its simplicity. The algorithm creates a binary brain
mask by thresholding the input image and uses a deformable model to refine the
initial mask. Despite its effectiveness, it sometimes fails to exclude non-
brain tissues, such as dura and eyes.

To address some of the limitations of BET, Iglesias _et al_. proposed Robust
Brain Extraction (ROBEX) [20]. ROBEX employs a trained random forest
classifier to distinguish brain and non-brain voxels. MRI WATERSHED [35] is
perhaps the most widely used algorithm due to its presence in the _FreeSurfer_
[13] library, as the name suggests it employs a watershed transformation
followed by a deformable model. Milchenko _et al_. [24] propose FACE MASK,
focusing on the protection of sensitive facial information in MRI data. This
method combines the benefits of defacing and skull-stripping by generating a
mask that can be used to blur out facial features while preserving brain
anatomy. Schimke _et al_. introduced a tool named QUICKSHEAR [34], which is
intended to compute a hyperplane in 3D space that separates the facial
features from the brain. DEFACE [5] is a deformable model that estimates which
voxels belong to the brain, said voxels can then be cut out for de-
identification purposes.

Finally, Van der Goten _et al_. suggest a remodeling-based approach using
GANs, CP-GAN [11]. This framework anonymously remodels non-brain tissues as
opposed to removing them. In CP-GAN the de-identification process is performed
by a conditional GAN that takes as input the patient’s brain and a convex hull
hinting at the metric extent of the skull to be generated. Our approach
follows this general framework of remodeling the face and skull, as opposed to
opting for something that would remove these parts altogether.

Vector Quantization. Vector-quantized autoencoders (VQ-VAE) [40, 31] have
become a staple in the computer vision community over the last years [33, 29]
and have also found widespread application in other fields [15, 9, 44]. VQ-
VAEs are neural networks that marry the concept of autoencoding with vector
quantization. The latter concept can best be described as _snapping_ a
continuous vector $c\in\mathbb{R}^{d}$ to a so-called codebook vector
$v_{i}\in\mathbb{R}^{d}$ ($i=1,\ldots,K$) according to $i=\text{argmax}_{j}\
d(c,v_{j})$ where $d(\cdot,\cdot)$ denotes some suitable distance function.
This concept enables efficient data compression and reconstruction. The
impossibility to derive gradients through the discrete argmax operation is
addressed by the _straight-through_ estimator trick[4].

A recurrent theme within recent works in computer vision is a two stage
approach: Firstly, a VQ-VAE (or a variant thereof) is trained in an
_unsupervised_ fashion on the whole dataset $\mathcal{D}$. The encoder of the
VQ-VAE is used to _translate_ each dataset element $X$ (image, spectrogram
etc.) into a discrete333It is also possible to translate to _continuous_
vectors which is done by taking the _inputs_ of the discretization operation
rather than its outputs representation $e(X)$ (_i.e_. a matrix of integers)
constituting a (fixed) dataset $\mathcal{D}_{e}$. Secondly, the actual task is
formulated. For image synthesis this could mean to learn how to stochastically
sample from $\mathcal{D}_{e}$. Said samples can then be fed through the
decoder to obtain actual samples mimicking the dataset $\mathcal{D}$. The
distribution can be learned via a multitude of model classes. The most common
ones are autoregressive [31], DPM-based [1, 7] or MAE-based [6, 30].

Image Synthesis. Image synthesis stands as a potent subset of machine learning
models, boasting a wide range of uses such as the generation of synthetic
data, identifying anomalies, and facilitating semi-supervised learning. A
diverse array of generative models exists, encompassing Variational
Autoencoders (VAEs) [22], Generative Adversarial Networks (GANs) [16], and,
more recently, diffusion-probabilistic models (DPMs) [37] and masked
autoencoders (MAEs) [12, 27, 28].

While diffusion-probabilistic models often produce state-of-the-art results,
they are challenged by long execution times during inference because for a
fixed number of steps $N$ (typically $1000$) the image needs to be
sequentially denoised $N$ times. DDIM [38] has shown that a non-Markovian
reinterpretation of the diffusion process can reduce the number of steps
significantly (_e.g_. to merely 50 steps). Yet, we argue that this is still a
lengthy process if one is concerned with high-resolution 3D scans, especially
if consumer-grade hardware is to be used.

Masked autoencoder (MAE) approaches are related to diffusion models in that
they aim to predict an image that was previously degraded. While DPMs use
noise to achieve this, MAEs stochastically mask out certain portions of the
data. Recently, MaskGIT [6] demonstrated that state-of-the-art performance on
ImageNet [10] can be achieved with a remasking scheme, Paella [30] uses a
similar training objective but simplifies the inference stage and furthermore
allows the model to re-predict tokens. Given said properties we opt to adapt
the Paella model to 3D in order to model the latent integer codes emitted by
the 3D VQ-VAE networks.

Figure 3: Latent Modeling Stage (MAE). Given the “paired dataset” featured in
Figure 2 we first extract the two components derived from the (continuous)
brain and the full skull without the brain. The former serves as a
_conditioning_ variable that informs the network about the metric constraints
(_e.g_. the size of the to-be-generated skull) of the synthesis process.
During training, the latter (a volume of integers encoding the full skull) is
stochastically-perturbed according to a _noise schedule_. The MAE’s task is
then to “unmask” the tokens. Since the targets are categorical we opt for the
_cross entropy_ loss.

##  3 Method

We follow the problem definition of CP-GAN [11]: Given a set of 3D scans
$(X^{(i)})_{i=1,\ldots,N}$ following a data distribution $\mathcal{P}_{X}$ and
having a resolution of $S^{3}$, we aim to find a mapping
$Y=G_{\Phi}(\gamma(X))$ parameterized by $\Phi$ that de-identifies a _raw_
scan $X$ and yields $Y$. The purpose of the _privacy transform_
$\gamma(\cdot)$ is to provide $G_{\Phi}(\cdot)$ with a minimal blueprint to
guide the skull synthesis without leaking information about the patient’s
facial features444As the skull is synthesized around the brain,
$\gamma(\cdot)$ should at least contain a binary brain mask. Fundamentally,
$\gamma(X)$ should at least include a representation of the brain to inform
the model about the metric constraints of the to-be-synthesized skull.

Figure 4: Test-time De-Identification. An _unseen_ MRI scan is de-identified
by applying the same preprocessing as during train time. Afterwards, the brain
is extracted yielding both a continuous and a binary brain representation. The
former is _vector-quantized_ to a compact volume of integers and used as the
_condition_ in the inference stage of the masked autoencoder. Starting from a
randomly-initialized $\hat{X}$ the network refines its estimate of how a skull
around the given brain could look like. A certain number of _renoising_ steps
is necessary as the network was not tasked to do one-step reconstructions
during the training phase. The final _de-identified_ scan is obtained by
_blending_ the original scan with the last estimate $\hat{X}$ where the
_binary_ brain representation acts as a _mask_. This step ensures that the
brain is preserved.

Overview. Synthesizing (3D) MRI scans is a challenging endeavor in terms of
memory as a volume contains $S^{3}$ (_e.g_. $S\in\\{128,256\\})$ voxels which
are ideally predicted by a single pass. VQ-VAEs paired with DPMs have already
been successfully shown to accurately model high-resolution (2D) imagery, most
prominently in the case of Stable Diffusion [33] that manages555Without using
additional super resolution techniques to synthesize (RGB) images of
resolution up to $256^{2}$. Applying a VQ-VAE on MRI scans is therefore
particularly promising as the number of voxels could in the most extreme case
be reduced from $256^{3}$ to a resolution of $64^{3}$, reducing the overall
number of voxels by a factor of 64. This is equivalent to a grayscale image of
resolution $512^{2}$, which is still pushing the limits of modern hardware and
algorithms.

One key objective of our work is to be able to run the training on consumer-
grade hardware with modest training/inference times. Having this in mind, a
bird’s eye view of our method can be described as follows: We first generalize
the VQ-VAE model to three dimensions. One instance of the VQ-VAE is tasked to
model the brains only, a second complementary instance models the full skulls
_without_ the brain. Applying both trained instances on an MRI dataset
produces a paired dataset where each item is a tuple of the two latent integer
volumes coming from each VQ-VAE encoder (see Fig. 2).

In a second step, we employ a Paella-style MAE that conditions on the
previously derived latents from the brain and models the distribution of the
latents pertaining to the full skulls (see Fig. 3). Unseen MRI scans can then
be de-identified by computing the latents associated to the brain and letting
the MAE synthesize a realistic skull using the brain latent as conditioning
variable (see Fig. 4).

Vector Quantization Stage. In order to reduce the memory requirements of high-
resolution MR imagery, we leverage the VQ-VAE framework to compress (parts of)
MRI scans into latent codes.

As a preparatory step we use ROBEX [20] to compute the (binary) brain mask
$B(X)$ of an MRI scan $X\in\mathcal{D}$. We need the brain representation for
two reasons: Firstly, it serves as a conditioning variable that informs the
synthesis about the proportions of the skull, and secondly, we require it to
later copy the original brain into the de-identified scan.

We then train a (3D) VQ-VAE for each of the two representations $B(X)\odot X$
(_“brain”_) and $\overline{B}(X)\odot X$ (_“skull”_) _in isolation_ where
$\odot$ denotes the Hadamard product and $\overline{B}(X)\ \hat{=}\ 1-B(X)$.
Both representations are complementary666By adding both one recovers $X$ to
each other in that $B(X)\odot X$ contains solely the brain (_i.e_. has non-
zero brain intensities) without featuring any of the remaining parts (_i.e_.
zero non-brain intensities) whereas the inverted properties hold for
$\overline{B}(X)\odot X$.

Training both models yields two encoder/decoder pairs $(e_{1},d_{1})$ and
$(e_{2},d_{2})$. After training, we deploy the two encoders to translate the
dataset of MRI scans $\mathcal{D}$ into a highly-compressed _paired_ dataset
of integer codes $\mathcal{D}_{e}=\\{(e_{1}[B(X)\odot
X)],e_{2}[\overline{B}(X)\odot X]\mid X\in\mathcal{D}\\}$. For the sake of
simplity, we assume that both $e_{1}[\cdot]$ and $e_{2}[\cdot]$ share the same
output shape $s\times s\times s$ where $s$ is a power of two that evenly
divides $S$. A conceptual depiction of this and the applied preprocessing can
be found in Fig. 2.

Architecture. We adapt the original VQ-VAE [40] model from (2D) images to (3D)
MRI scans. Given a tensor777$Z$ either corresponds to $B(X)\odot X$ or
$\overline{B}(X)\odot X$ $Z\in\mathbb{R}^{S\times S\times S}$, its resolution
$S\in\\{128,256\\}$ and a desired output resolution $s\ll S$ (a power of two),
we first feed $Z$ through an encoder. The encoder consists of a stack of
$\log_{2}{\nicefrac{{S}}{{s}}}$ convolutional blocks (the _“feature block”_)
followed by a vector quantizer.

We apply layer normalization [2] along the channel dimension before each
convolution. Let us express the final output of the feature block as
$f(Z)\in\mathbb{R}^{s\times s\times s\times L}$ where $L$ denotes an arbitrary
embedding dimension.

The vector quantizer, part of the encoder, operates independently on the
$s^{3}$ vectors of its input $f(Z)$, producing a volumetric output
$e(Z)\in\\{0,\ldots,N_{\text{CV}}-1\\}^{s\times s\times s}$ of integers where
each integer encodes the index of the closest codebook vector in a (trainable)
codebook of size $N_{\text{CV}}$. The decoder is a mirror image of the
encoder, with its own set of parameters. We employ a mean square loss between
$Z$ and the output of the decoder. To account for space constraints we moved
specifics as well as the used augmentation regimen to the _Appendix_.

Latent Modeling Stage. Having the compressed latent integer representations of
$E_{1}\hat{=}e_{1}[B(X)\odot X]$ (_“brain”_) and
$E_{2}\hat{=}e_{2}[\overline{B}(X)\odot X]$ (_“skull”_) now at our disposal,
we describe a generative MAE $H_{\Phi}(E_{2}|E_{1})$ that predicts the latter
(integer) skull representation while conditioning on the former (integer)
brain representation.

Architecture. Recall that both $E_{1}$ and $E_{2}$ are 3D volumes/grids of
integers from the set $\\{0,\ldots,N_{\text{CV}}-1\\}^{s\times s\times s}$. In
order to model the distribution of $E_{2}$ the MAE methodology requires that
one defines a (time-dependent) scheme to _perturb_ the integers in the $E_{2}$
representation. We decide to use _Paella_ -style [30] perturbations, _i.e_. a
value $v\in\\{0,\ldots,N_{\text{CV}}-1\\}$ is kept with probability
$\alpha_{t}$ and resampled from the same set with a probability
$1-\alpha_{t}$, where the time $t\sim\mathcal{U}(0,1)$ is sampled
independently for each element in the batch. This can be efficiently
implemented by bundling the point-wise masking decisions in a tensor $M$ and
the three-dimensional uniform variates in $U$:

$\tilde{E_{2}}=(1-M)\odot E_{2}+M\odot U$ (1)

As the spatial shape of the input agree with the spatial shape of the output,
we decide to use a U-Net architecture (with layer normalization) to denoise
the perturbed representations. During training the U-Net is given three
inputs: Namely, the perturbed input $\tilde{E_{2}}$, the condition
$E_{1}\hat{=}e_{1}[B(X)\odot X]$ (also called $\gamma(X)$) and the time $t$ to
provide a general sense of how corrupted $\tilde{E_{2}}$ is. We choose to
merge $E_{1}$ and $\tilde{E_{2}}$ by using a simple (channel-wise)
concatenation before feeding it into the network. The scalar $t$ is fed to an
MLP (w. GeLU activations [18]) at each layer deriving two parameters
$\delta_{1},\delta_{2}$ which are used to shift and rescale an input $Z$ as
$Z^{\prime}=Z\cdot(\delta_{1}+1)+\delta_{2}$. The U-Net finally produces a
logit tensor $Y_{\text{logits}}\in\mathbb{R}^{s\times s\times s\times
N_{\text{CV}}}$ that can be used to compute the cross entropy loss w.r.t. the
(unperturbed) $E_{2}$. For further details refer to Fig. 3 and the _Appendix_.

| _User-based_ | _Model-based_
---|---|---
| ADNI | OASIS-3 | ADNI | OASIS-3
| $128^{3}$ | $256^{3}$ | $128^{3}$ | $256^{3}$ | $128^{3}$ | $256^{3}$ | $128^{3}$ | $256^{3}$
BLACK | $20.17\pm 25.45$ | $20.30\pm 26.33$ | $18.37\pm 10.90$ | $18.54\pm 11.08$ | $19.66\pm 2.25$ | $18.75\pm 2.08$ | $19.86\pm 0.89$ | $19.98\pm 1.19$
BLURRED | $45.05\pm 28.91$ | $46.29\pm 29.42$ | $41.63\pm 15.57$ | $42.93\pm 18.06$ | $87.46\pm 7.65$ | $86.54\pm 8.61$ | $97.05\pm 2.33$ | $97.31\pm 1.67$
ORIGINAL | $55.33\pm 30.70$ | $58.93\pm 29.24$ | $61.86\pm 15.00$ | $59.27\pm 16.94$ | $100.00\pm 0.00$ | $100.00\pm 0.00$ | $100.00\pm 0.00$ | $100.00\pm 0.00$
MRI WATERSHED | $19.03\pm 25.31$ | $21.48\pm 27.25$ | $22.56\pm 13.11$ | $22.20\pm 14.41$ | $44.75\pm 4.09$ | $47.03\pm 5.06$ | $67.76\pm 4.39$ | $67.50\pm 7.39$
DEFACE | $43.53\pm 30.16$ | $45.58\pm 29.53$ | $38.14\pm 9.58$ | $43.17\pm 16.04$ | $99.27\pm 0.86$ | $99.22\pm 1.51$ | $99.78\pm 0.36$ | $99.85\pm 0.20$
QUICKSHEAR | $38.51\pm 30.18$ | $39.81\pm 30.43$ | $40.70\pm 16.68$ | $35.85\pm 13.78$ | $98.79\pm 0.87$ | $95.66\pm 1.98$ | $99.81\pm 0.21$ | $99.83\pm 0.20$
FACE MASK v1 | $48.55\pm 30.49$ | $50.24\pm 30.54$ | $50.23\pm 17.39$ | $52.68\pm 14.32$ | $96.31\pm 2.71$ | $98.75\pm 1.59$ | $99.65\pm 0.50$ | $99.71\pm 0.23$
FACE MASK v2 | $38.05\pm 28.79$ | $43.60\pm 29.90$ | $33.02\pm 17.93$ | $35.85\pm 15.00$ | $94.42\pm 5.11$ | $98.06\pm 1.20$ | $99.78\pm 0.31$ | $99.63\pm 0.29$
CP-GAN | $28.46\pm 27.95$ | ✗ | $30.47\pm 14.30$ | ✗ | $56.11\pm 5.05$ | ✗ | $56.40\pm 2.89$ | ✗
CP-MAE | $\bm{23.14\pm 27.88}$ | $\bm{25.91\pm 27.65}$ | $\bm{22.56\pm 11.77}$ | $\bm{23.41\pm 15.27}$ | $\bm{39.91\pm 9.49}$ | $\bm{41.74\pm 6.91}$ | $\bm{48.82\pm 4.32}$ | $\bm{58.19\pm 4.04}$

Table 1: De-Identification Quality. We compare CP-MAE against traditional
methods, CP-GAN and four control methods (above ”- -” line) in terms of their
de-identification capabilities. For this, users on _Amazon Mechanical Turk_
were asked to map an original rendering onto its de-identified counterpart
(c.f.  _”User-based”_). Similarly, we harnessed a neural network trained
within a _contrastive learning_ framework to perform the very same task (c.f.
_”Model-based”_). Both scenarios involve finding the correct option among five
alternatives (of which four belong to a different subject) and the showcased
values correspond to the percentage of correct guesses ($\pm$ s.d.). The ”✗”
for CP-GAN denotes that the method does not support a resolution of $256^{3}$.

(a) ORIGINAL

(b) BLACK

(c) BLURRED

(d) MRI WATERSHED

(e) DEFACE

(f) QUICKSHEAR

(g) FACE MASK v1

(h) FACE MASK v2

(i) CP-GAN

(j) CP-MAE

Figure 5: De-Identification Quality: Renderings. A single subject (from
OASIS-3) de-identified by various methods (and control methods) featuring
differing levels of aggressiveness.

Sampling Stage. After having trained the MAE, we execute the following
sampling steps: (i) Compute the brain mask $B(X^{\prime})$, (ii) compute the
integer representation $E_{1}=e_{1}[B(X^{\prime})\odot X^{\prime}]$, (iii) use
$H_{\Phi}$ to sample $E_{2}$ using $E_{1}$ as a conditioning variable, and
finally (iv), pass $E_{2}$ through the decoder $d_{2}(\cdot)$ to recover a
skull that harmonizes well with the brain in $X^{\prime}$. To ensure that the
de-identified scan accurately represents the real brain and not a hallucinated
version thereof, we implement a simple blending scheme that is inspired by CP-
GAN in order to copy the brain from the original $X$:

$Y=(1-B(X))\odot d_{2}(\hat{E_{2}})+B(X)\odot X$ (2)

where $\hat{E_{2}}$ denotes the latent as estimated by letting the MAE
condition on $E_{1}\hat{=}e_{1}[B(X)\odot X]$ (_i.e_. the brain).

Privacy Considerations. During the training phase, the model can theoretically
”remember” or associate the original scan given only a representation of the
brain. However, when we split the data into training and testing sets, we find
that the de-identification property holds on the testing set. This is because
the network has never seen brains from the test set before and a test-set
brain hence does not provide enough information to accurately recreate the
skull or any distinctive facial features of an unseen individual. We validate
this finding within our user study in the _Experiments_ section.

##  4 Experiments

In this section, we conduct a comprehensive set of experiments to evaluate the
performance and effectiveness of our proposed methodology for MRI scan de-
identification. We aim to assess the capabilities of our approach in terms of
preserving privacy while maintaining the diagnostic value of the scans. To
achieve this, we carefully design experiments that encompass various aspects,
including downstream task performance and the degree of de-identification
achieved. To gauge how much a de-identification method affects a real-world
downstream task, we employ various segmentation tasks (_i.e_. SIENAX [21],
FIRST [26] & FASTSURFER [19]) and quantify how the segmentation changes in
comparison to running the methods on raw data. We also conduct a user study on
Amazon Mechanical Turk in which people are tasked with mapping de-identified
scans to their raw/original counterpart given 3D renderings of the respective
scans.

| ADNI | OASIS-3
---|---|---
| $128^{3}$ | $256^{3}$ | $128^{3}$ | $256^{3}$
MRI WATERSHED | ✗ | $0.229\pm 0.208$ | ✗ | ✗
DEFACE | $0.977\pm 0.016$ | $0.977\pm 0.022$ | $0.981\pm 0.007$ | $0.979\pm 0.013$
QUICKSHEAR | $0.982\pm 0.007$ | $0.980\pm 0.015$ | $0.982\pm 0.044$ | $0.981\pm 0.037$
FACE MASK v1 | $0.959\pm 0.007$ | $0.956\pm 0.011$ | $0.950\pm 0.022$ | $0.958\pm 0.012$
FACE MASK v2 | $0.959\pm 0.007$ | $0.952\pm 0.062$ | $0.950\pm 0.043$ | $0.958\pm 0.012$
CP-GAN | $0.980\pm 0.013$ | ✗ | $0.984\pm 0.003$ | ✗
CP-MAE | $\bm{0.985\pm 0.010}$ | $\bm{0.986\pm 0.011}$ | $\bm{0.994\pm 0.005}$ | $\bm{0.981\pm 0.012}$

Table 2: Downstream Tasks: Brain segmentation. We utilize the widely-
established SIENAX tool from the FSL suite to gauge how much de-identification
methods affect the outcome of a common downstream task. A perfect method gives
rise to the same segmentation. We quantify the difference between the
segmentation maps by computing the class-averaged Dice score over four
classes. Higher is better. (”✗” denotes that the model is not available for
this resolution, ”✗” means that the method is failing in more than 90% of all
cases)

Datasets and implementation details. In this study, we use the well-
established and publicly available datasets ADNI [42, 43] (2,172 T1-weighted
MRI scans) and OASIS-3 [23] (2,556 T1-weighted MRI scans). We use a per-
patient train ($80\%$)-test ($20\%$) split to partition the data.

For training our networks, we leverage two NVIDIA RTX 4090 GPUs, each equipped
with 24 GiB of GPU memory. Our implementation relies on the PyTorch framework
[25]. To optimize memory consumption, we utilize the _inductor_ backend for
model compilation, which has proven to be highly advantageous. Moreover, we
adopt the recently introduced Lion optimizer [8], which enables us to scale up
the model size compared to other available optimizers. These implementation
details contribute to the efficiency and effectiveness of our training process
and facilitate robust experimentation.

Preprocessing & Benchmarks. The MRI scans are standardized to RAS+
orientation, resampled to a common resolution ($128^{3}$ resp. $256^{3}$),
bias field corrected by the N4 method [39], intensity normalized using the
_Fuzzy C-means_ algorithm [32] and finally co-registered to a fixed scan with
an affine deformation.

We juxtapose our findings by comparing against six de-identification methods.
The first three belong to the removal-based class and are comprised of
QUICKSHEAR [34], FACE MASK [24] and DEFACE [5]. On the extreme end we also
include MRI WATERSHED as a control method which removes all non-brain tissue.
Apart from the published version of FACE MASK we also cover an
updated888https://github.com/mmilch01/mask_face2 version of the toolkit
created by the same author. For the sake of clarity we refer to the two
variants as FACE MASK v1 and FACE MASK v2. Furthermore, we benchmark against
CP-GAN [11] which is so far the only member of the remodeling-based class. A
visual depiction of the methods can be found in Figure 5. All procedures were
applied with default settings on images of resolution $128^{3}$ and $256^{3}$
(with the exception of CP-GAN which synthesizes only up to a resolution of
$128^{3}$). A detailed explanation of each method can be found in the
_Appendix_.

De-Identification Quality: User-based Study. In an effort to assess the
resilience of various de-identification approaches, we replicate the user
study of CP-GAN utilizing Amazon Mechanical Turk. Participants were exposed to
a series of quiz questions designed to bypass various de-identification
techniques. Each question was created by sampling an untouched (3D) rendering
of a scan, accompanied by five more (3D) renderings created by a single fixed
method (_e.g_. CP-MAE). Within these renderings, one was the de-identified
counterpart of the original patient scan. The goal for the participants was to
correctly match the de-identified rendering with the rendering of the
unaltered scan. The renderings were produced using volumetric ray tracing and
feature a $45^{\circ}$ side profile of the patients. This is in contrast to
the user study conducted by CP-GAN which uses a frontal view. We argue that a
side view presents more discernible information to the user as various
features such as skull size, ears and head posture can be factored into the
decision process.

For maintaining comparability, we incorporated control tasks with varying
degrees of de-identification: ORIGINAL, which signifies the untouched scan
with no de-identification applied; BLURRED, introducing a mild degree of
obfuscation by blurring the 2D renderings; BLACK, consistently presenting an
entirely black image; and MRI WATERSHED, a brain extraction tool which only
retains the brain tissue.

For each method (or control method) we posted $220$ distinct999One question
for each element in the test set questions and a single question was given to
five workers, meaning that $1100$ questions were answered in total for each
method. Taking the two resolutions into account, we collected $2\times
10\times 1100=22,000$ answers. Since the workers are presented with five
alternatives the optimal guessing rate is $20\%$ which corresponds to perfect
de-identification. More specifics, including a sample question, are available
in the _Appendix_. Illustrative renderings are presented in Figure 5 whereas
the results are showcased in Table 1.

We observe that CP-MAE attains almost optimal de-identification results,
outperforming the removal-based methods by a margin of $15$ to $25$ percent.
CP-MAE also outperforms CP-GAN albeit by a smaller margin. Although both
methods share a similar methodology, we explain the difference by the fact
that CP-GAN is additionally conditioning on a tightly-fit (filled) convex hull
around the skull. This effectively constrains the synthesis to produce skull
shapes that are more correlated to the skull shape of the raw/original scan as
compared to CP-MAE which solely conditions on the brain.

ADNIOASIS-3FIRSTFASTSURFERDEFACEFACE MASK v1CP-GANQUICKSHEARFACE MASK v2CP-
MAEDEFACEFACE MASK v1CP-GANQUICKSHEARFACE MASK v2CP-
MAE0.800.850.900.951.000.800.850.900.951.00$128^{3}$$256^{3}$

Figure 6: Downstream Tasks: Subcortical segmentation. In a similar vein to
Table 2 we analyze how de-identification affects downstream tasks. FIRST is an
integral part of the FSL suite and FASTSURFER is a deep learning-based
alternative to FREESURFER that runs orders of magnitude faster. The depicted
values are the _class-averaged_ Dice scores over $15$ classes for FIRST and
resp. $78$ classes for FASTSURFER. Due to the much lower performance of MRI
WATERSHED ($\approx 0.21$, averaged over all configurations) we decided to
exclude the method to improve visual discernability among the other methods.
Higher values are preferable.

De-Identification Quality: Model-based Study. Building on the methodology of
the previous experiment, we evaluated the efficacy of different de-
identification models by attempting to circumvent them. In this iteration, we
employed a neural network to determine the similarity between original
(unaltered) and de-identified renderings, an approach that is similar in
fashion to the one conducted for CP-GAN albeit with a much more powerful
similarity-quantification network.

To achieve this, we adopted a similar metric learning strategy by training a
Siamese network $h(\cdot,\cdot)$ . This network is designed to discern whether
two given inputs are associated with the same patient. Specifically, for an
original rendering $x$ and a de-identified rendering $y$, the network $h$
functions by independently embedding $x$ and $y$ through a sub-network
$\tilde{h}$—comprising a ResNet-18 block [17] (pre-trained on ImageNet-1K), a
flattening layer, and a fully-connected layer—and subsequently calculating the
Euclidean distance between the resulting embeddings. Thus, the network’s
output is given by
$h(x,y)\triangleq\lVert\tilde{h}(x)-\tilde{h}(y)\rVert_{2}$. We employed the
Triplet Margin loss function as outlined in [3], with the margin set to $1$.
The previously established hold-out dataset was partitioned, allocating $80\%$
of patients to the training set and the remaining $20\%$ to the test set. To
calculate standard deviations we repeat the aforementioned process to obtain
$k=10$ independent folds. Detailed training procedures are outlined in the
_Appendix_.

Table 1 illustrates the network’s proficiency in overcoming de-identification
methods. The values are obtained by selecting a patient from one of the test
sets (pertaining to a specific fold) and a corresponding original scan. A de-
identification method $m$ was then chosen, and the $m$-rendering of the scan
was designated as the correct option. The other four options were compiled by
randomly choosing $m$-renderings from scans of different patients. The
network’s decision is given by the option $y$ that exhibits the shortest
Euclidean distance to $\tilde{h}(x)$.

First and foremost, we observe that the Siamese approach is much more powerful
in terms of defeating the de-identification methods. Similar to the user-based
study, we observe that CP-MAE exhibits superior de-identification performance
in comparison to traditional de-identification methods and compares favorably
against CP-GAN.

Effect of De-Identification on Medical Analyses. It is of utmost importance
that a scan de-identified with a de-identification method does not negatively
affect the analysis performed by downstream tasks. To explore the level of
impairment we employ three established methods. Firstly, SIENAX [21] measures
the brain volume and provides a segmentation into four classes (total brain
volume, white matter, gray matter & ventricular CSF) and is the de-facto
standard to assess how individuals compare in terms of tissue composition.
Secondly, FIRST is a Bayesian (non deep-learning based) point distribution-
based approach to perform efficient and accurate subcortical segmentations.
Subcortical segmentations are an important means in the analysis of
neurodegenerative diseases. Thirdly, to gauge how deep learning-based
downstream tasks are effected, we further employ FASTSURFER [19] which is a
recently suggested tool that also addresses the subcortical segmentation task.
It is meant to offer a more efficient and robust alternative to the widely-
established FREESURFER [14] algorithm.

To quantify whether one of these segmentation methods is thrown off by a de-
identification method, we replicate the relativistic evaluation approach
already employed by CP-GAN: First a segmentation method is run on an original
MRI scan yielding a segmentation map $\mathcal{S}$, then the segmentation
method is run once again, this time on the de-identified counterpart of the
scan, producing a new segmentation map $\mathcal{S}^{\prime}$. Possible
differences are then quantified using the Dice coefficent w.r.t. each
segmentation class. We report the _class-averaged Dice score_ for SIENAX in
Table 2 and for FIRST resp. FASTSURFER in Figure 6 .

We observe that CP-MAE outperforms the other methods w.r.t. both SIENAX and
FIRST and exhibits slightly better Dice scores for FASTSURFER. Furthermore, we
observe that de-identification with MRI WATERSHED (full-brain extraction)
causes downstream tasks to either get impaired or to fail completely.

##  5 Conclusion

We present a novel MRI de-identification method called CP-MAE which robustly
anonymizes an MRI scan by employing a generative model that conditions only on
the brain. Within our work, we have shown how the challenge of handling high-
resolution MR imagery can be addressed by leveraging two instances of a VQ-VAE
together with an MAE. The MAE operates in the low-dimensional space defined by
the VQ-VAE pairs and is able to stochastically produce high-fidelity MR scans
that retain the brain of the raw scan but exhibit vastly different appearance
outside of the brain. Using this construction we are able to improve the
state-of-the-art 3D de-identification resolution of CP-GAN from $128^{3}$ to
$256^{3}$, representing an eight-fold increase in the number of voxels. We
validate the de-identification quality of our model by conducting a user study
in which CP-MAE outperforms all other methods by moderate to substantial
margins. Analyzing how de-identification affects common brain segmentation
tasks, we conclude that our model affects the segmentation maps only to a
negligible degree. Future work will include the CT domain.

## References

  * Austin et al. [2021] Jacob Austin, Daniel D. Johnson, Jonathan Ho, Daniel Tarlow, and Rianne van den Berg.  Structured denoising diffusion models in discrete state-spaces.  In _Advances in Neural Information Processing Systems_ , 2021.
  * Ba et al. [2016] Jimmy Lei Ba, Jamie Ryan Kiros, and Geoffrey E. Hinton.  Layer normalization, 2016.
  * Balntas et al. [2016] Vassileios Balntas, Edgar Riba, Daniel Ponsa, and Krystian Mikolajczyk.  Learning local feature descriptors with triplets and shallow convolutional neural networks.  In _British Machine Vision Conference_ , 2016.
  * Bengio et al. [2013] Yoshua Bengio, Nicholas Léonard, and Aaron Courville.  Estimating or propagating gradients through stochastic neurons for conditional computation, 2013.
  * Bischoff-Grethe et al. [2021] A. Bischoff-Grethe et al.  Defacing method for magnetic resonance imaging data.  _Journal of Neuroscience Methods_ , 350:109034, 2021.
  * Chang et al. [2022] Huiwen Chang, Han Zhang, Lu Jiang, Ce Liu, and William T. Freeman.  Maskgit: Masked generative image transformer, 2022.
  * Chen et al. [2022] Ting Chen, Ruixiang Zhang, and Geoffrey E. Hinton.  Analog bits: Generating discrete data using diffusion models with self-conditioning.  _ArXiv_ , abs/2208.04202, 2022.
  * Chen et al. [2023] Xiangning Chen, Chen Liang, Da Huang, Esteban Real, Kaiyuan Wang, Yao Liu, Hieu Pham, Xuanyi Dong, Thang Luong, Cho-Jui Hsieh, Yifeng Lu, and Quoc V. Le.  Symbolic discovery of optimization algorithms, 2023.
  * Cifka et al. [2021] Ondrej Cifka, Alexey Ozerov, Umut Simsekli, and Gael Richard.  Self-supervised VQ-VAE for one-shot music style transfer.  In _ICASSP 2021 - 2021 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2021.
  * Deng et al. [2009] Jia Deng, Wei Dong, Richard Socher, Li-Jia Li, Kai Li, and Li Fei-Fei.  Imagenet: A large-scale hierarchical image database.  In _2009 IEEE Conference on Computer Vision and Pattern Recognition_ , pages 248–255, 2009.
  * der Goten et al. [2021] Lennart Alexander Van der Goten, Tobias Hepp, Zeynep Akata, and Kevin Smith.  Conditional de-identification of 3d magnetic resonance images, 2021.
  * Devlin et al. [2018] Jacob Devlin, Ming-Wei Chang, Kenton Lee, and Kristina Toutanova.  Bert: Pre-training of deep bidirectional transformers for language understanding.  _arXiv preprint arXiv:1810.04805_ , 2018.
  * Fischl [2012a] Bruce Fischl.  Freesurfer.  _NeuroImage_ , 62(2):774–781, 2012a.  20 YEARS OF fMRI.
  * Fischl [2012b] Bruce Fischl.  Freesurfer.  Neuroimage, 2012b.
  * Garbacea et al. [2019] Cristina Garbacea, Aaron van den Oord, Yazhe Li, Felicia S C Lim, Alejandro Luebs, Oriol Vinyals, and Thomas C Walters.  Low bit-rate speech coding with VQ-VAE and a WaveNet decoder.  In _ICASSP 2019 - 2019 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2019.
  * Goodfellow et al. [2014] Ian Goodfellow, Jean Pouget-Abadie, Mehdi Mirza, Bing Xu, David Warde-Farley, Sherjil Ozair, Aaron Courville, and Yoshua Bengio.  Generative adversarial nets.  In _Advances in Neural Information Processing Systems_ , pages 2672–2680, 2014.
  * He et al. [2015] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Deep residual learning for image recognition.  _CoRR_ , abs/1512.03385, 2015.
  * Hendrycks and Gimpel [2020] Dan Hendrycks and Kevin Gimpel.  Gaussian error linear units (gelus), 2020.
  * Henschel et al. [2019] Leonie Henschel, Sailesh Conjeti, Santiago Estrada, Kersten Diers, Bruce R. Fischl, and Martin Reuter.  Fastsurfer - a fast and accurate deep learning based neuroimaging pipeline.  _NeuroImage_ , 219:117012 – 117012, 2019.
  * Iglesias et al. [2011] J. E. Iglesias et al.  Robust brain extraction across datasets and comparison with publicly available methods.  _IEEE Transactions on Medical Imaging_ , 30(9):1617–1634, 2011.
  * Jenkinson et al. [2002] Mark Jenkinson, Paul Bannister, Michael Brady, and Stephen Smith.  Improved optimization for the robust and accurate linear registration and motion correction of brain images.  _Neuroimage_ , 17(2):825–841, 2002.
  * Kingma and Welling [2013] Diederik P Kingma and Max Welling.  Auto-encoding variational bayes.  _arXiv preprint arXiv:1312.6114_ , 2013.
  * LaMontagne et al. [2019] Pamela J. LaMontagne, Tammie LS. Benzinger, John C. Morris, Sarah Keefe, Russ Hornbeck, Chengjie Xiong, Elizabeth Grant, Jason Hassenstab, Krista Moulder, Andrei G. Vlassenko, Marcus E. Raichle, Carlos Cruchaga, and Daniel Marcus.  Oasis-3: Longitudinal neuroimaging, clinical, and cognitive dataset for normal aging and alzheimer disease.  _medRxiv_ , 2019.
  * Milchenko and Marcus [2013] Mikhail Milchenko and Daniel Marcus.  Obscuring surface anatomy in volumetric imaging data.  _Neuroinformatics_ , 11(1):65–75, 2013.
  * Paszke et al. [2019] Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Köpf, Edward Z. Yang, Zach DeVito, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala.  Pytorch: An imperative style, high-performance deep learning library.  _CoRR_ , abs/1912.01703, 2019.
  * Patenaude et al. [2011] Brian Patenaude, Stephen Smith, David Kennedy, and Mark Jenkinson.  A bayesian model of shape and appearance for subcortical brain segmentation.  _NeuroImage_ , 56:907–22, 2011.
  * Radford et al. [2018] Alec Radford, Karthik Narasimhan, Tim Salimans, and Ilya Sutskever.  Improving language understanding by generative pre-training.  _OpenAI Blog_ , 2018.
  * Radford et al. [2019] Alec Radford, Jeff Wu, Rewon Child, David Luan, Dario Amodei, and Ilya Sutskever.  Language models are unsupervised multitask learners.  _OpenAI Blog_ , 1(8), 2019.
  * Ramesh et al. [2021] Aditya Ramesh, Mikhail Pavlov, Gabriel Goh, Scott Gray, Chelsea Voss, Alec Radford, Mark Chen, and Ilya Sutskever.  Zero-shot text-to-image generation, 2021.
  * Rampas et al. [2022] Dominic Rampas, Pablo Pernias, Elea Zhong, and Marc Aubreville.  Fast text-conditional discrete denoising on vector-quantized latent spaces, 2022.
  * Razavi et al. [2019] Ali Razavi, Aäron van den Oord, and Oriol Vinyals.  Generating diverse high-fidelity images with VQ-VAE-2.  _CoRR_ , abs/1906.00446, 2019.
  * Reinhold et al. [2019] Jacob C Reinhold, Blake E Dewey, Aaron Carass, and Jerry L Prince.  Evaluating the impact of intensity normalization on MR image synthesis.  In _Medical Imaging 2019: Image Processing_ , page 109493H. International Society for Optics and Photonics, 2019.
  * Rombach et al. [2022] Robin Rombach, Andreas Blattmann, Dominik Lorenz, Patrick Esser, and Björn Ommer.  High-resolution image synthesis with latent diffusion models, 2022.
  * Schimke et al. [2011] Nakeisha Schimke, Mary Kuehler, and John Hale.  Preserving privacy in structural neuroimages.  In _Lecture Notes in Computer Science (including subseries Lecture Notes in Artificial Intelligence and Lecture Notes in Bioinformatics)_ , pages 301–308. Springer, Berlin, Heidelberg, 2011.
  * Segonne et al. [2004] Florent Segonne, Anders Dale, Evelina Busa, Maureen Glessner, David Salat, Horst Hahn, and Bruce Fischl.  A hybrid approach to the skull stripping problem in mri.  _NeuroImage_ , 22:1060–75, 2004.
  * Smith et al. [2002] S. M. Smith et al.  Fast robust automated brain extraction.  _Human Brain Mapping_ , 17(3):143–155, 2002.
  * Sohl-Dickstein et al. [2020] Jascha Sohl-Dickstein, Eric A Weiss, Niru Maheswaranathan, and Anitha Martin.  Denoising diffusion probabilistic models.  _arXiv preprint arXiv:2006.11239_ , 2020.
  * Song et al. [2020] Jiaming Song, Chenlin Meng, and Stefano Ermon.  Denoising diffusion implicit models.  _CoRR_ , abs/2010.02502, 2020.
  * Tustison et al. [2010] Nicholas Tustison, Brian Avants, Philip Cook, Yuanjie Zheng, Alexander Egan, Paul Yushkevich, and James Gee.  N4itk: improved n3 bias correction.  _Medical Imaging, IEEE Transactions on_ , 29:1310 – 1320, 2010.
  * van den Oord et al. [2017] Aäron van den Oord, Oriol Vinyals, and Koray Kavukcuoglu.  Neural discrete representation learning.  _CoRR_ , abs/1711.00937, 2017.
  * Villegas et al. [2022] Ruben Villegas, Mohammad Babaeizadeh, Pieter-Jan Kindermans, Hernan Moraldo, Han Zhang, Mohammad Taghi Saffar, Santiago Castro, Julius Kunze, and Dumitru Erhan.  Phenaki: Variable length video generation from open domain textual description, 2022.
  * Weiner et al. [2017] Michael W. Weiner, Dallas P. Veitch, Paul S. Aisen, Laurel A. Beckett, Nigel J. Cairns, Robert C. Green, Danielle Harvey, Clifford R. Jack, William Jagust, John C. Morris, Ronald C. Petersen, Jennifer Salazar, Andrew J. Saykin, Leslie M. Shaw, Arthur W. Toga, and John Q. Trojanowski.  The Alzheimer’s Disease Neuroimaging Initiative 3: Continued innovation for clinical trial improvement, 2017.
  * Wyman et al. [2013] Bradley T. Wyman, Danielle J. Harvey, Karen Crawford, Matt A. Bernstein, Owen Carmichael, Patricia E. Cole, Paul K. Crane, Charles Decarli, Nick C. Fox, Jeffrey L. Gunter, Derek Hill, Ronald J. Killiany, Chahin Pachai, Adam J. Schwarz, Norbert Schuff, Matthew L. Senjem, Joyce Suhy, Paul M. Thompson, Michael Weiner, and Clifford R. Jack.  Standardization of analysis sets for reporting results from ADNI MRI data, 2013.
  * Xue et al. [2019] Yifan Xue, Michael Ding, and Xinghua Lu.  Supervised vector quantized variational autoencoder for learning interpretable global representations, 2019.