# Investigation of Reflection-Based Measurements of Microwave Kinetic
Inductance Detectors in the Optical Bands

Jie Hu ID GEPI, Observatoire de Paris, PSL Université, CNRS, 75014 Paris,
France  Université de Paris, CNRS, Astroparticule et Cosmologie, F-75013
Paris, France  Faouzi Boussaha ID GEPI, Observatoire de Paris, PSL Université,
CNRS, 75014 Paris, France  Paul Nicaise ID GEPI, Observatoire de Paris, PSL
Université, CNRS, 75014 Paris, France  Christine Chaumont ID GEPI,
Observatoire de Paris, PSL Université, CNRS, 75014 Paris, France  Maria
Appavou  GEPI, Observatoire de Paris, PSL Université, CNRS, 75014 Paris,
France  Viet Dung Pham ID Université de Paris, CNRS, Astroparticule et
Cosmologie, F-75013 Paris, France  Michel Piat ID Université de Paris, CNRS,
Astroparticule et Cosmologie, F-75013 Paris, France<EMAIL_ADDRESS>
###### Abstract

In this paper, we investigate the single photon response from the reflection
of the Microwave Kinetic Inductance Detector (MKID) array. Reflection
measurements are carried out using two configurations: one is measured
simultaneously with the transmission, and the other is obtained with a single-
ended MKID array terminated with an open load. Compared with the transmission,
reflection measurements significantly reduce the readout noise of the single-
ended MKID array. This is also reflected in the improvement of the median
energy resolving power by around 20-30% under pulsed photon illumination at
$\lambda=405$ nm, mainly due to an increase in the size of the resonance
circle on the IQ plane. This method has the potential to be used to read out
large MKID arrays.

††preprint: APS/123-QED

The Microwave Kinetic Inductance DetectorDay  _et al._ (2003) (MKID) has
emerged as a groundbreaking technology with wide-ranging applications in
astrophysics, spanning from millimeter-waveAdam, R.  _et al._ (2018) and
opticalMeeker  _et al._ (2018); Jie  _et al._ (2023) to X-ray
observationsSwimmer (2023) for its wideband operation, intrinsic frequency
domain multiplexity and relatively simple fabrication. In the optical and
near-infrared bands, MKIDs demonstrate a distinctive blend of excellent energy
resolution and the ease of scalability into large arraysMazin  _et al._
(2012); Gao  _et al._ (2012); de Visser  _et al._ (2021); Zobrist  _et al._
(2022), providing a significant advantage in applications like extreme faint
object detection and low-resolution spectroscopy like SPIAKIDJie  _et al._
(2023).

MKIDs typically consist of high-quality lumped superconducting resonators in
the optical and near-infrared bands. These resonators, as discussed in Beldi
_et al._ (2019); de Visser  _et al._ (2021); Hu  _et al._ (2021); Zobrist  _et
al._ (2022); Nicaise  _et al._ (2022), are constructed with an interdigitated
capacitor (IDC) and a meander inductor. When photons with energies
$h\nu>2\Delta$ (2$\Delta$ the Cooper pair binding energy) are absorbed in the
meander, they generate phonons, disrupt Cooper pairs, and produce
quasiparticles (QPs). These additional QPs increase the kinetic inductance,
lowering the resonator’s resonance frequency and quality factor. This response
can be detected through changes in the phase and amplitude of the resonator,
respectively, by applying a probing tone at the resonance frequency.

The performance of optical MKIDs is characterized by their energy-resolving
power $R$, which currently faces domination by three distinct noise sources.
The first is the noise caused by the escape of hot phonons from the meander to
the substrate, as discussed in the literaturede Visser  _et al._ (2021);
Zobrist  _et al._ (2022). The second is the two-level system noise (TLS)
generated within the substrate and at the superconductor/substrate interface
Gao  _et al._ (2008). The last is the readout noise from the system, primarily
stemming from the low noise amplifier (LNA), which plays a significant role.
Efforts have been dedicated to reducing this readout noise, notably through
implementing superconducting parametric amplifiersZobrist  _et al._ (2019,
2022), which are still not commercially available for reading out MKIDs.

Figure 1: MKIDs T-junction circuit

Presently, MKIDs are typically read out via the transmission of the resonator.
However, due to the inherent nature of MKIDs, strong reflections can occur at
resonance, especially in the case of high-quality superconducting resonators
where the reflected signal ($S_{11}$) can surpass the transmitted signal
($S_{21}$). Consequently, by utilizing $S_{11}$ for resonator readout, there
is potential to diminish the sensitivity of the resonator to readout noise
within the system. This signal, however, is not traditionally measured in the
conventional homodyne configuration.

Figure 2: Measurement Setup with reflection ($S_{11}$) and transmission
($S_{21}$). The $S_{11}$-open is obtained by terminating the output port of
the MKIDs with an open load. The reflected signal is obtained by putting an
isolator at the input of the resonator. The isolator is currently placed on
the 4K stage due to the limited space on the 0.1K stage.

In this paper, we investigate the single-photon response of TiN/Ti/TiN
MKIDsJie  _et al._ (2023) measured by reflection. To capture the reflected
signal, a cryogenic circulator is positioned at the input of the MKID array,
with the reflected signal ($S_{11}$) read from the circulator’s third port—a
methodology commonly employed for qubit readout. Wang  _et al._ (2014, 2021).
Two distinct configurations are compared. In the first configuration, both
$S_{11}$ and $S_{21}$, taken as a reference, are measured simultaneously. In
contrast, only $S_{11}$ is measured in the second configuration when the MKID
array is terminated with an open load. The latter is referred to as single-
ended MKIDs. Across all configurations, we observed a single-photon response,
with the single-ended MKID array demonstrating the most favorable median $R$,
approximately 20-30%

The equivalent circuit of the MKID shunted-coupled to a transmission line with
a characteristic impedance $Z_{0}$ can be generally represented in Fig. 1.
$Z_{1}$ and $Z_{2}$ are general descriptions of the port impedance, which can
either be the impedance mismatch or the impedance of the transmission lines.
$Z_{3}$ is the MKID impedance. $S_{11}$ can be obtained as

$\displaystyle S_{11}=\frac{AZ_{3}+\Gamma_{0}}{1+2Z_{3}/\hat{Z}_{0}}$ (1)

with

$\displaystyle\hat{Z}_{0}$
$\displaystyle=\frac{2(Z_{1}+Z_{0})(Z_{2}+Z_{0})}{2Z_{0}+Z_{1}+Z_{2}}$ (2)
$\displaystyle A$
$\displaystyle=\frac{Z_{1}+Z_{2}}{(Z_{0}+Z_{1})(Z_{0}+Z_{2})}$ (3)
$\displaystyle\Gamma_{0}$ $\displaystyle=\frac{Z_{1}-Z_{0}}{Z_{1}+Z_{0}}$ (4)

$\hat{Z}_{0}$ and $\Gamma_{0}$ are the equivalent characteristic impedance of
the transmission line and the reflection coefficient at the MKID array input.
$S_{21}$ of the MKIDs can be expressed asZobrist  _et al._ (2022)

$\displaystyle
S_{21}=\frac{2Z_{0}}{2Z_{0}+Z_{1}+Z_{2}}(1-\frac{1}{1+2Z_{3}/\hat{Z}_{0}})$
(5)

withKhalil  _et al._ (2012); Zobrist  _et al._ (2022)

$\displaystyle\frac{1}{1+2Z_{3}/\hat{Z_{0}}}=\frac{Q/Q_{c}-2jQx_{a}}{1+2jQx_{g}}$
(6)

where $Q$ and $Q_{c}$ are, respectively, the quality factor and the coupling
quality factor of the resonator, $x_{g}=(f-f_{r})/f_{r}$ is the detuning of
the generator frequency from resonance, and $x_{a}=\delta f_{a}/f_{r}$ is the
fractional detuning of the original resonance frequency due to the impedance
mismatches. By substituting Eq. (6) into Eq. (A),

$\displaystyle S_{11}$
$\displaystyle=\gamma-(\gamma-\Gamma_{0})\frac{Q/Q_{c}-2jQx_{a}}{1+2jQx_{g}}$
(7)

Then, $S_{11}$ can be fitted with

$\displaystyle
S_{11}=(\gamma-\Gamma_{0})(\gamma_{0}-\frac{Q/Q_{c}-2jQx_{a}}{1+2jQx_{g}})$
(8)

with $\gamma_{0}=\gamma/(\gamma-\Gamma_{0})$ and
$\gamma=(Z_{1}+Z_{2})/(2Z_{0}+Z_{1}+Z_{2})$ as a complex fitting parameter.
The amplitude of $S_{11}$ and the location of its center of the resonance
circle in the IQ plane will depend on the value of $\gamma_{0}$. For a matched
resonator, i.e. $Z_{1}=Z_{2}\rightarrow 0$,  $\gamma_{0}\rightarrow 0$.
$S_{11}$ will show a maximum at the resonance frequency. When
$Z_{2}\rightarrow\infty$ and $Z_{1}\rightarrow 0$, $\gamma_{0}\rightarrow
1/2$, the magnitude of $S_{11}$ will be similar to $S_{21}$ and shows a
minimum at the resonance. A more detailed mathematical procedure for Eq. (S26)
can be found in Appendix A.

Figure 3: Comparison of the single-photon response of the MKID at 4.018 GHz
from the transmission and reflection at $T_{bath}=200~{}\text{mK}$ (A)-(C):
Measured pulse statistics from $S_{21}$, $S_{11}$ and $S_{11}$-Open.
$\lambda_{p}$ is the fitted mean value of the Poisson distribution.
$\sigma_{0}$ and $\sigma_{1}$ are the standard deviations of the first and
second peaks. (D)-(E): Comparison of the amplitude of the MKID and its
trajectory on the IQ plane. The black dots indicate the readout frequency. The
dashed lines are the fitted curves with Eq. (S21) and (S26). (F): The averaged
single-photon response with hundreds of pulse events with $n=1$ as shown in
(A)-(C). (G): Comparison between the noise spectrum and the pulse shown in
(F), obtained by Fourier transformation. The dashed line ($--$) is the noise
in phase, and the dash-dot line ($-\cdot$) is the noise in amplitude.

The detailed measurement setup for the MKIDs is illustrated in Fig. 2. The
MKIDs undergo characterization within a pre-cooled two-stage pulse tube
Adiabatic Demagnetization Refrigerator (ADR). A 1.5mm-thick niobium cylinder,
sheets of metglas 2714a surrounding the MKIDs, and a 2mm-thick cylinder of
Cryophy are employed on the 4K stage to shield from magnetic interference.
Each MKID is individually read out using a standard homodyne mixing scheme.

The input signal, generated by a signal generator, is attenuated by a
programmable attenuator at room temperature and a 20 dB attenuator at 4 K. A
circulator is positioned at the input of the MKIDs to capture the $S_{11}$
signal. Due to space limitations on the 100 mK stage, attenuators on the 1 K
and 100 mK stages are temporarily removed to acquire the $S_{11}$ signal. The
thermal noise above 100 GHz is well-attenuated thanks to the thermal filters
on the 100 mK stage. Thus, the background noise at the MKIDs should remain
lowBaselmans  _et al._ (2012), although the thermal noise at the resonance
frequency band remains. It should also be noted that placing the circulator on
the 100 mK stage will reduce the system noise, thus further improving $R$.

The transmitted and reflected signals from MKIDs are initially amplified by an
LNA on the 4 K stage and further amplified by two room-temperature amplifiers.
Subsequently, the transmitted and reflected signals are down-converted to DC
by two IQ mixers and then sampled by an oscilloscope. The readout power at the
detector’s input is estimated to be approximately -100 dBm.

Illumination of the MKIDs array is achieved through an optical fiber
positioned 35 mm above the pixels. A 405 nm laser, modulated by a 250 Hz pulse
from a pulse generator with a width of around 50 ns, is employed for
illumination. The laser output power is adjusted outside the cryostat using a
digital step attenuator. The pulse response of the MKID is sampled by an
oscilloscope at a frequency of 5-10 MHz.

Two different configurations are employed for the measurement of the MKIDs.
First, $S_{11}$ and $S_{21}$ are measured simultaneously. Second, the
responses of the MKIDs terminated with an open load at the output port of the
array are measured, as depicted in Fig. 2, and the corresponding result is
labeled as $S_{11}$-open.

Our array comprises 30 $\times$ 30 MKIDs. Each MKID consists of a 10 nm/10
nm/10 nm-thick TiN/Ti/TiN trilayer meander connected to a niobium capacitor.
The meander size is $36\times 36~{}\mu\text{m}^{2}$. The measured critical
temperature $T_{c}$ of the trilayer film is around 1.75 K. A more detailed
description of the design of the pixels can be found in our previous
publicationJie  _et al._ (2023).

As depicted in Fig.3, we first conducted a comparison of the measured single-
photon response of the MKID obtained from $S_{11}$, $S_{21}$, and
$S_{11}$-Open. In all cases, the MKID resonating at around 4.018 GHz was
measured with a readout of approximately 2 dB below saturation at the
frequency showing the maximum response on the IQ plane. Each measurement
recorded around 10,000 pulse events and was processed by the standard wiener
filteringGao  _et al._ (2012). The resulting pulse height statistics are
illustrated in Fig.3-(A)-(C), and these data were fitted using the convolution
between the Gaussian and Poisson distributions, asMezzena  _et al._ (2019).

$\displaystyle
p(\phi)=\sum_{n=0}^{\infty}\frac{\lambda_{p}^{n}}{n!}e^{-(\phi-n\phi_{1}-\phi_{0})^{2}/(2\sigma_{n}^{2})}$
(9)

Here, $\lambda_{p}$ is the mean value of the Poisson distribution, $\phi_{0}$
is the mean value of the 0-th peak that is the background response from the
chip, $\sigma_{n}$ is the variance of the n-th peak, and $\phi_{1}$ is the
difference between the peaks, which is the response of a single photon.

Figure 4: (A)-(B): Comparison of amplitudes of $S_{21}$, $S_{11}$ and
$S_{11}$-Open respectively. There are two transmission zeros around 4.25 GHz
and 5.35 GHz.

Here, we consider the width of the 0-photon peak is dominated by the system
noise $\sigma_{s}$, i.e., the two-level system (TLS) noise $\sigma_{TLS}$ from
the MKID itself and the readout noise $\sigma_{r}$.

$\displaystyle\sigma_{0}^{2}$ $\displaystyle=\sigma_{TLS}^{2}+\sigma_{r}^{2},$
(10)

The width of the 1-photon peak is dominated by the readout noise and the noise
related from the photon absorption $\sigma_{p}$, i.e., the hot phonon escapede
Visser  _et al._ (2021) and the current non-uniformity in the meanderZobrist
_et al._ (2019) as

$\displaystyle\sigma_{1}^{2}$ $\displaystyle=\sigma_{0}^{2}+\sigma_{p}^{2}$
(11)

Here, $\sigma_{p}$ is dominated by the phonon escape noise. Since the TLS is
mainly in the phase response of MKID, $\sigma_{r}$ can be obtained by
convolving the amplitude response in the 0-photon peak with the wiener filter.
We list the contribution of each noise of the resonator at 4.0175 GHz in Table
1. It can be seen that the $\sigma_{TLS}$ and $\sigma_{p}$ are consistent with
each other. The difference in the $R$s of $S_{21}$, $S_{11}$ and $S_{11}$-open
is due to the different levels of the readout noise. The $R$ of $S_{11}$-Open
is dominated by the TLS and the phonon escape noise.

Table 1: Fitted parameters and contribution of each noise to $R$ for the resonator at 4.0175 GHz. |  $Q_{c}$ |  $Q_{i}$ |  $\gamma_{0}$ |  $\sigma_{TLS}$ |  $\sigma_{r}$ |  $\sigma_{p}$
---|---|---|---|---|---|---
|  (k) |  (k) |  |  (O. F.) |  (O. F.) |  (O. F.)
$S_{21}$ |  32k |  21k |  1.0 |  0.61 |  0.56 |  0.61
$S_{11}$ |  25k |  17k |  0.06 + 1.0j |  0.60 |  0.72 |  0.69
$S_{11}$-Open |  16k |  19k |  0.27+0.35j |  0.58 |  0.38 |  0.59

The increase in the readout noise in $S_{11}$ can be attributed to two
possible reasons. First, the presence of the reflection chain especially the
circulator can lead to an increase in readout noise. Second, the MKID is read
out at the frequency exhibiting maximum response in the IQ plane, resulting in
similar output powers at the output port of the MKID across all cases, as
depicted in Fig. 3-(D)-(E), where the black dots represent the position of the
readout signals. Consequently, the signal obtained from the reflection may not
be significantly higher than the $S_{21}$’s. The dashed lines are the fitted
results with Eq. (S21) and Eq. (S26) and the fitted parameters are listed in
Table 1, which shows the internal quality factor $Q_{i}$ are in good agreement
with each other. The initial values strongly influence the fitted parameters,
yet their impact on estimating the pulse response is minimal.

The reduction of the readout noise in $S_{11}$-open is due to the increase of
the radius of the resonance circle $r$, which is shown in Fig.3-(E). as there
is no significant difference in the pulse response shown in Fig.3-(F). The
readout noise is related to $r$ asBARENDS (2009)

$\displaystyle S_{DSB}=\frac{k_{B}T_{r}}{r^{2}P_{r}},$ (12)

where $k_{B}$ is the Boltzmann constant $T_{r}$ is the equivalent input noise
temperature of the readout system. $P_{r}$ is the level of the readout power.
Thus, the readout noise would be reduced when there is an increase in the
radius of the resonance circle. This is also supported by the noise spectrum
shown in Fig. 3-(G). Eq. (12) is qualitatively in agreement with the
measurement by taking into account that the readout power of $S_{11}$-open is
2 dB lower than that of $S_{21}$. For matched resonators (i.e, $Z_{1}=0$ and
$Z_{2}=0$ for the case of $S_{11}$ and $S_{21}$, and $Z_{1}=0$ and
$Z_{2}=\infty$ for $S_{11}$-Open), the radius of the resonance circle of
$S_{11}$, $S_{21}$ and $S_{11}$-open can be expressed as

$\displaystyle
r=\begin{cases}\frac{Q_{i}}{2(Q_{i}+Q_{c})}\sqrt{4(Q_{c}x_{a})^{2}+1}&S_{11}~{}\&~{}S_{21}\\\
\frac{Q_{i}}{Q_{i}+Q_{c}^{\prime}}\sqrt{4(Q_{c}^{\prime}x_{a})^{2}+1}&S_{11}-\text{Open}\end{cases}$
(13)

$Q_{i}$ is the internal quality factor of the resonator and
$1/Q=1/Q_{c}+1/Q_{i}$. The coupling quality factor $Q_{c}^{\prime}$ of the
$S_{11}$-Open will change due as

$\displaystyle
Q_{c}=\frac{C}{\omega_{r}C_{c}^{2}}\Re{\frac{2Z_{0}+Z_{1}+Z_{2}}{(Z_{0}+Z_{1})(Z_{0}+Z_{2})}}$
(14)

where $C$ and $C_{c}$ are, respectively, the capacitors of the resonance and
the coupling capacitor to the feedline. $\omega_{r}$ is the angular resonance
frequency. $Q_{c}^{\prime}$ will reduce when $Z_{0}+Z_{2}\rightarrow\infty$,
which is agreement with the $Q_{c}$ fitted from the $S_{21}$ and
$S_{11}$-open. In this case, for a matched resonator ($x_{a}=0$), the radius
of the resonance circle will be increased by more than a factor of 2. A
comprehensive mathematical derivation for determining the radius of the
resonance circle can be found in Appendix B.

Figure 5: Comparison of the energy resolution between the $S_{11}$, $S_{21}$,
and $S_{11}$-open over a bandwidth of around 100 MHz at 4.19-4.29 GHz,
4.75-4.95 GHz and 5.30-5.40 GHz. (A): The measured amplitude of the $S_{11}$,
$S_{21}$ and $S_{11}$-open. (B): The measured energy resolution versus
frequency. The input power is optimized for each pixel to maximize $R$. (C):
Statistics of $R$ in the measured bands, around 10 pixels in the measured
bands, do not show a single photon response.

The readout scheme utilizing $S_{11}$-Open presents a notable advantage for
systems featuring significant ripples in the baseline particularly due to the
absence of bridges that balance the CPW groundMeeker (2017). In Fig. 4, we
compare the amplitudes of $S_{11}$, $S_{21}$, and $S_{11}$-Open across the
frequency range from 4.00 GHz to 6.00 GHz. Notably, where a minimum exists in
$S_{21}$, there is a corresponding maximum in $S_{11}$, and vice versa.
Additionally, as the frequency surpasses 5.50 GHz, the magnitude of $S_{11}$
begins to exceed that of $S_{21}$, elucidating the unexpected additional
reduction in the magnitude of $S_{21}$ within our system. It can be seen that
terminating the array with an open load remarkably diminishes the grounding
effect to approximately $\pm 3$ dB, a level comparable to state-of-the-art
MKID arrays Baselmans  _et al._ (2017); Shibo  _et al._ (2018). Note the
amplitude of the $S_{11}$ and $S_{11}$-Open depends strongly on the
$\gamma_{0}$ in Eq. S26, which is determined by the matching condition of the
resonator to the circuit. When $\gamma_{0}$ is real positive, $S_{11}$ and
$S_{11}$-Open will have similar amplitude with $S_{21}$. The amplitude of
$S_{11}$-Open shown in Fig.3 is a special case with the center of the
resonance circle close to the origin of the IQ plane.

Further investigations will be performed for larger MKIDs with multiple
readout lines to check if the bridges that balance the ground can be
eliminated.

We conducted a comparison of the energy resolution ($R$) across three
frequency bands: 4.19-4.29 GHz, 4.75-4.95 GHz, and 5.30-5.40 GHz, where two of
them exhibit a minimum in the amplitude, as depicted in Fig. 5-(A). The
measured $R$s within the corresponding bands are illustrated in Fig. 5-(B).
Despite the amplitude of $S_{21}$ being approximately 10 dB lower than that of
$S_{11}$, the measured $R$ from $S_{21}$ remains comparable to $S_{11}$, while
$R$ from $S_{11}$-Open demonstrates a significant improvement. This trend is
consistent for the resonances observed in the 4.75-4.95 GHz and 5.30-5.40 GHz
bands.

The statistical analysis of the $R$s depicted in Fig. 5-(B) is presented in
Fig. 5-(C). Notably, the median $R$ is enhanced from 2.0 with $S_{21}$ to 2.6
with $S_{11}$-Open, and the maximum $R$ reaches approximately 3.5.

The impedance transformation of the readout line, denoted as ($Z_{0}$ \+
$Z_{2}$) in Eq. (14), could potentially approach zero, leading to an increase
in readout noise as $r$ decreases, as described in Eq. (12) and Eq. (13). To
mitigate this issue, one effective approach is to position the resonators at
multiples of half wavelength to the open load. Typically, an MKID array
exhibits resonances between 4 GHz and 8 GHz, with the wavelength of an 8 GHz
signal measuring around 16 mm on sapphire. By arranging the pixels within $\pm
1~{}$mm around the half wavelength, approximately 10 pixels can be
accommodated. However, achieving MKIDs in the optical band with resonance
frequencies at desired positions presents challenges due to the small size of
the pixels. Further investigations are warranted to optimize the performance
of MKIDs with a reflection readout scheme. Moreover, it’s intriguing to note
that despite our MKIDs not being designed to be terminated by an open load,
we’re observing an increase in the median $R$-value, suggesting the
effectiveness of our method. Additionally, we believe that the proposed
readout scheme remains applicable for MKID arrays in the millimeter bands, as
indicated by the validity of Eq. (12) and Eq. (13). Even if there’s a slight
deviation in the resonance frequency from the design, it should still be
manageableShu  _et al._ (2021).

In conclusion, we have implemented an MKID readout scheme based on reflection.
Our observations indicate no significant improvement in the median $R$ when
$S_{11}$ is measured together with $S_{21}$. Adopting a scheme to read out the
single-end MKID array yields an enhancement in the median $R$ of the single-
end MKID array by 20-30%. This improvement is primarily attributed to the
increase in the radius of the resonance circle. It is worth mentioning that
this method can also be implemented in systems utilizing superconducting
parametric amplifiers Zobrist  _et al._ (2019, 2022).

## Availability of data

The data supporting the findings of this study are available from the
corresponding authors upon reasonable request.

## AUTHOR DECLARATIONS

## Conflict of Interest

The authors have no conflicts to disclose.

###### Acknowledgements.

The authors would like to thank Shibo Shu from The Institute of High Energy
Physics of the Chinese Academy of Sciences for the suggestion for the MKIDs
measurement with reflection, as well as Florent Reix, Josiane Firminy, and
Thibaut Vacelet from Paris Observatory for assembly and mounting the devices.
This work is supported by the European Research Council (ERC) through Grant
835087 (SPIAKID) and the UnivEarthS Labex program.

##  Appendix A Reflection of the Resonator with Mis-match

The $Z_{in}$ of the system can be calculated as

$\displaystyle Z_{in}$
$\displaystyle=Z_{1}+\frac{Z_{3}(Z_{2}+Z_{0})}{Z_{3}+(Z_{2}+Z_{0})}$ (S15)

The $S_{11}$ can be obtained as

$\displaystyle S_{11}$
$\displaystyle=\frac{Z_{in}-Z_{0}}{Z_{in}+Z_{0}}=\frac{Z_{3}(Z_{1}+Z_{2})+(Z_{1}-Z_{0})(Z_{0}+Z_{2})}{Z_{3}(2Z_{0}+Z_{1}+Z_{2})+(Z_{0}+Z_{1})(Z_{0}+Z_{2})}$
(S16)

with

$\displaystyle\hat{Z}_{0}=\frac{2(Z_{1}+Z_{0})(Z_{2}+Z_{0})}{2Z_{0}+Z_{1}+Z_{2}}$
(S17)

Then, $S_{11}$ is

$\displaystyle S_{11}=\frac{AZ_{3}+\Gamma_{0}}{1+2Z_{3}/\hat{Z_{0}}}$ (S18)

with

$\displaystyle A$
$\displaystyle=\frac{Z_{1}+Z_{2}}{(Z_{0}+Z_{1})(Z_{0}+Z_{2})}$ (S19)
$\displaystyle\Gamma_{0}$ $\displaystyle=\frac{Z_{1}-Z_{0}}{Z_{1}+Z_{0}}$
(S20)

The $S_{21}$ can be expressed as

$\displaystyle
S_{21}=\frac{2Z_{0}}{2Z_{0}+Z_{1}+Z_{2}}(1-\frac{1}{1+2Z_{3}/\hat{Z_{0}}})$
(S21)

with

$\displaystyle\frac{1}{1+2Z_{3}/\hat{Z}_{0}}=\frac{Q/Q_{c}-2jQx_{a}}{1+2jQx_{g}}$
(S22)

$S_{11}$ can be segmented into two parts as

$\displaystyle
S_{11}=\gamma(1-\frac{1}{1+2Z_{3}/\hat{Z}_{0}})+\frac{\Gamma_{0}}{1+2Z_{3}/\hat{Z}_{0}}$
(S23)

with $\gamma=(Z_{1}+Z_{2})/(2Z_{0}+Z_{1}+Z_{2})$.

By taking into account of Eq. (S21), Eq. (S23) becomes

$\displaystyle
S_{11}=\frac{Z_{1}+Z_{2}}{2Z_{0}}S_{21}-\Gamma\frac{Q/Q_{c}-2jQx_{a}}{1+2jQx_{g}}$
(S24)

or

$\displaystyle S_{11}$
$\displaystyle=\gamma-(\gamma-\Gamma_{0})\frac{Q/Q_{c}-2jQx_{a}}{1+2jQx_{g}}$
(S25)

Then, the $S_{11}$ can be fitted with

$\displaystyle
S_{11}=(\gamma-\Gamma_{0})(\gamma_{0}-\frac{Q/Q_{c}-2jQx_{a}}{1+2jQx_{g}})$
(S26)

with $\gamma_{0}=\gamma/(\gamma-\Gamma_{0})$ as a complex fitting parameter.
For the single-end resonator, $Z_{2}\rightarrow\infty$ and $\gamma\rightarrow
1$. If the resonator is perfectly matched, $Z_{1}=0$ and $\Gamma_{0}=-1$.
Thus, $S_{11}$ can be expressed as

$\displaystyle S_{11}=1-2\frac{Q/Q_{c}-2jQx_{a}}{1+2jQx_{g}}$ (S27)

which is similar with the expression of $S_{11}$ found in previous
publications.

##  Appendix B The Coupling Quality Factor and Resonance circle

The coupling quality factor can be calculated as follows. The energy stored in
the capacitor and inductor are the same at the resonance frequency, which is

$\displaystyle E=2\times\frac{1}{2}CV^{2}$ (S28)

where $V$ is the voltage of the capacitor. The energy dissipated outside the
resonator is

$\displaystyle
P=|I|^{2}\Re{R}\approx\omega^{2}C_{c}^{2}V^{2}\Re{\frac{(Z_{0}+Z_{1})(Z_{0}+Z_{2})}{2Z_{0}+Z_{1}+Z_{2}}}$
(S29)

$Q_{c}$ can then be calculated by definition as

$\displaystyle Q_{c}=\omega\frac{E}{P}=\frac{C}{\omega
C_{c}^{2}}\Re{\frac{2Z_{0}+Z_{1}+Z_{2}}{(Z_{0}+Z_{1})(Z_{0}+Z_{2})}}.$ (S30)

We compare the radius of the resonance circle between the $S_{11}$, $S_{21}$
and $S_{11}$-Open. The term in the bracket of Eq. (S26)

$\displaystyle S_{11}^{b}=\gamma_{0}-\frac{Q/Q_{c}-2jQx_{a}}{1+2jQx_{g}}$
(S31)

is a circle on the IQ plane. The three points on the resonance circle can be
determined as

$\displaystyle S_{11}^{b}=\begin{cases}\gamma_{0}&x\rightarrow\infty\\\
\gamma_{0}-Q/Q_{c}-2jQx_{a}&x=0\\\
\gamma_{0}-\frac{Q/Qc-2jQx_{a}}{1+jQ}&x=1/2\\\ \end{cases}$ (S32)

Here we assume $\gamma_{0}=\gamma_{r}+j\gamma_{i}$, then the center of and the
radius of the resonance circle is

$\displaystyle[x_{0}^{b},~{}y_{0}^{b}]$
$\displaystyle=[\gamma_{r}-\frac{Q}{2Q_{c}},\gamma_{i}+2Qx_{a}]$ (S33)
$\displaystyle R_{b}$
$\displaystyle=\frac{Q}{2Q_{c}}\sqrt{4(Q_{c}x_{a})^{2}+1}$ (S34)

Taking into account of the coefficients, the radius of the resonator circle of
the $S_{11}$ is

$\displaystyle
R_{r}=|\gamma_{21}-\Gamma_{0}|\frac{Q}{2Q_{c}}\sqrt{4(Q_{c}x_{a})^{2}+1}$
(S35)

If the resonator is perfectly matched, $\gamma_{21}=0$. $\Gamma_{0}=-1$. In
the case that the $Z_{2}\rightarrow\infty$, which is the case that the MKIDs
is terminated with an open load. $\gamma_{21}=1$ In the case

$\displaystyle
R_{r}=\begin{cases}\frac{Q_{i}}{2(Q_{i}+Q_{c})}\sqrt{4(Q_{c}x_{a})^{2}+1}&Z_{1}=0~{}~{}\&~{}~{}Z_{2}=0\\\
\frac{Q_{i}}{Q_{i}+Q_{c}^{\prime}}\sqrt{4(Q_{c}^{\prime}x_{a})^{2}+1}&Z_{1}=0~{}~{}\&~{}~{}Z_{2}=\infty\end{cases}$
(S36)

with $1/Q=1/Q_{i}+1/Q_{c}$. Thus, there is an increase of a factor of 2 in the
resonator terminated with an open load. Similarly, the radius of $S_{21}$

$\displaystyle
R_{t}=|\frac{2Z_{0}}{2Z_{0}+Z_{1}+Z_{2}}|\frac{Q}{2Q_{c}}\sqrt{4(Q_{c}x_{a})^{2}+1}$
(S37)

In the situation that $Z_{1}=Z_{2}=0$,

$\displaystyle R_{t}=\frac{Q_{i}}{2(Q_{i}+Q_{c})}\sqrt{4(Q_{c}x_{a})^{2}+1}$
(S38)

## References

  * Day  _et al._ (2003) P. K. Day, H. G. LeDuc, B. A. Mazin, A. Vayonakis, and J. Zmuidzinas, Nature 425, 817 (2003).
  * Adam, R.  _et al._ (2018) Adam, R., Adane, A., Ade, P. A. R., André, P., Andrianasolo, A., Aussel, H., Beelen, A., Benoît, A., Bideaud, A., Billot, N., Bourrion, O., Bracco, A., Calvo, M., Catalano, A., Coiffard, G., Comis, B., De Petris, M., Désert, F.-X., Doyle, S., Driessen, E. F. C., Evans, R., Goupy, J., Kramer, C., Lagache, G., Leclercq, S., Leggeri, J.-P., Lestrade, J.-F., Macías-Pérez, J. F., Mauskopf, P., Mayet, F., Maury, A., Monfardini, A., Navarro, S., Pascale, E., Perotto, L., Pisano, G., Ponthieu, N., Revéret, V., Rigby, A., Ritacco, A., Romero, C., Roussel, H., Ruppin, F., Schuster, K., Sievers, A., Triqueneaux, S., Tucker, C., and Zylka, R., A&A 609, A115 (2018).
  * Meeker  _et al._ (2018) S. R. Meeker, B. A. Mazin, A. B. Walter, P. Strader, N. Fruitwala, C. Bockstiegel, P. Szypryt, G. Ulbricht, G. Coiffard, B. Bumble, G. Cancelo, T. Zmuda, K. Treptow, N. Wilcer, G. Collura, R. Dodkins, I. Lipartito, N. Zobrist, M. Bottom, J. C. Shelton, D. Mawet, J. C. van Eyken, G. Vasisht, and E. Serabyn, Publications of the Astronomical Society of the Pacific 130, 065001 (2018).
  * Jie  _et al._ (2023) H. Jie, N. Paul, B. Faouzi, M. Jean-Marc, C. Christine, M. Alexine, R. Florent, F. Josiane, V. Thibaut, P. V. Dung, P. Michel, C. Elisabetta, and B. Piercarlo, Journal of Low Temperature Physics (2023), 10.1007/s10909-023-03018-5.
  * Swimmer (2023) N. J. Swimmer, Ph. d thesis, University of California, Santa Barbara (2023).
  * Mazin  _et al._ (2012) B. A. Mazin, B. Bumble, S. R. Meeker, K. O’Brien, S. Mchugh, and E. Langman, Optics Express 20, 1503 (2012).
  * Gao  _et al._ (2012) J. Gao, M. R. Vissers, M. O. Sandberg, F. C. S. da Silva, S. W. Nam, D. P. Pappas, D. S. Wisbey, E. C. Langman, S. R. Meeker, B. A. Mazin, H. G. Leduc, J. Zmuidzinas, and K. D. Irwin, Applied Physics Letters 101, 142602 (2012).
  * de Visser  _et al._ (2021) P. J. de Visser, S. A. de Rooij, V. Murugesan, D. J. Thoen, and J. J. Baselmans, Phys. Rev. Applied 16, 034051 (2021).
  * Zobrist  _et al._ (2022) N. Zobrist, W. H. Clay, G. Coiffard, M. Daal, N. Swimmer, P. Day, and B. A. Mazin, Physical Review Letters 129, 017701 (2022).
  * Beldi  _et al._ (2019) S. Beldi, F. Boussaha, J. Hu, A. Monfardini, A. Traini, F. Levy-Bertrand, C. Chaumont, M. Gonzales, J. Firminy, F. Reix, M. Rosticher, S. Mignot, M. Piat, and P. Bonifacio, Opt Express 27, 13319 (2019).
  * Hu  _et al._ (2021) J. Hu, F. Boussaha, J.-M. Martin, P. Nicaise, C. Chaumont, S. Beldi, M. Piat, and P. Bonifacio, Applied Physics Letters 119 (2021), 10.1063/5.0074103.
  * Nicaise  _et al._ (2022) P. Nicaise, J. Hu, J.-M. Martin, S. Beldi, C. Chaumont, P. Bonifacio, M. Piat, H. Geoffray, and F. Boussaha, Journal of Low Temperature Physics 209, 1242 (2022).
  * Gao  _et al._ (2008) J. Gao, J. Zmuidzinas, A. Vayonakis, P. Day, B. Mazin, and H. Leduc, Journal of Low Temperature Physics 151, 557 (2008).
  * Zobrist  _et al._ (2019) N. Zobrist, G. Coiffard, B. Bumble, N. Swimmer, S. Steiger, M. Daal, G. Collura, A. B. Walter, C. Bockstiegel, N. Fruitwala, I. Lipartito, and B. A. Mazin, Applied Physics Letters 115, 213503 (2019).
  * Wang  _et al._ (2014) C. Wang, Y. Y. Gao, I. M. Pop, U. Vool, C. Axline, T. Brecht, R. W. Heeres, L. Frunzio, M. H. Devoret, G. Catelani, L. I. Glazman, and R. J. Schoelkopf, Nat Commun 5, 5836 (2014).
  * Wang  _et al._ (2021) H. Wang, S. Singh, C. R. H. McRae, J. C. Bardin, S. X. Lin, N. Messaoudi, A. R. Castelli, Y. J. Rosen, E. T. Holland, and D. P. Pappas, Quantum Science and Technology 6, 035015 (2021).
  * Khalil  _et al._ (2012) M. S. Khalil, M. J. A. Stoutimore, F. C. Wellstood, and K. D. Osborn, Journal of Applied Physics 111, 054510 (2012).
  * Baselmans  _et al._ (2012) J. Baselmans, S. Yates, P. Diener, and P. de Visser, Journal of Low Temperature Physics 167, 360 (2012).
  * Mezzena  _et al._ (2019) R. Mezzena, M. Faverzani, E. Ferri, A. Giachero, B. Margesin, A. Nucciotti, A. Puiu, and A. Vinante, Journal of Low Temperature Physics 199, 73 (2019).
  * BARENDS (2009) R. BARENDS, Ph. d thesis, Tu Delft (2009).
  * Meeker (2017) S. R. Meeker, Doctor thesis, University of California Santa Barbara (2017).
  * Baselmans  _et al._ (2017) J. J. A. Baselmans, J. Bueno, S. J. C. Yates, O. Yurduseven, N. Llombart, K. Karatsu, A. M. Baryshev, L. Ferrari, A. Endo, D. J. Thoen, P. J. de Visser, R. M. J. Janssen, V. Murugesan, E. F. C. Driessen, G. Coiffard, J. Martin-Pintado, P. Hargrave, and M. Griffin, Astronomy & Astrophysics 601, A89 (2017).
  * Shibo  _et al._ (2018) S. Shibo, M. Calvo, J. Goupy, S. Leclercq, A. Catalano, A. Bideaud, A. Monfardini, and E. F. C. Driessen, Applied Physics Letters 113, 082603 (2018).
  * Shu  _et al._ (2021) S. Shu, M. Calvo, J. Goupy, S. Leclercq, A. Catalano, A. Bideaud, A. Monfardini, and E. F. C. Driessen, Applied Physics Letters 119 (2021), 10.1063/5.0058335.


# Pixel-weighted Multi-pose Fusion for Metal Artifact Reduction in X-ray
Computed Tomography

Diyu Yang1, Craig A. J. Kemp2, Soumendu
Majee${}^{3}\textsuperscript{\textsection}$, Gregery T. Buzzard1, and Charles
A. Bouman1
1Purdue University-Main Campus, West Lafayette, IN 47907.
2Eli Lilly and Company, Indianapolis, IN 46225.
3Samsung Research America, Mountain View, CA 94043.
Diyu Yang was supported by Eli Lilly and Company. Gregery T. Buzzard was
supported by NSF grant CCF-1763896. Charles A. Bouman was supported by the
Showalter Trust.

###### Abstract

X-ray computed tomography (CT) reconstructs the internal morphology of a three
dimensional object from a collection of projection images, most commonly using
a single rotation axis. However, for objects containing dense materials like
metal, the use of a single rotation axis may leave some regions of the object
obscured by the metal, even though projections from other rotation axes (or
poses) might contain complementary information that would better resolve these
obscured regions.

In this paper, we propose pixel-weighted Multi-pose Fusion to reduce metal
artifacts by fusing the information from complementary measurement poses into
a single reconstruction. Our method uses Multi-Agent Consensus Equilibrium
(MACE), an extension of Plug-and-Play, as a framework for integrating
projection data from different poses. A primary novelty of the proposed method
is that the output of different MACE agents are fused in a pixel-weighted
manner to minimize the effects of metal throughout the reconstruction. Using
real CT data on an object with and without metal inserts, we demonstrate that
the proposed pixel-weighted Multi-pose Fusion method significantly reduces
metal artifacts relative to single-pose reconstructions.

###### Index Terms:

Inverse problems, Computed tomography, Model based reconstruction, Plug-and-
play, Multi-agent consensus equilibrium, Metal artifact reduction.

§§footnotetext: This work was done while Soumendu Majee was employed at Purdue
University.

##  I Introduction

X-ray computed tomography (CT) imaging is widely used in industrial [1] and
medical [2] applications for non-destructive visualization of internal sample
morphology. X-ray CT uses a series of projection images from various angles to
reconstruct a 3D array of attenuation coefficients that describe the sample
[3]. Traditional CT reconstruction methods, such as Filtered Back Projection
(FBP), use projection images acquired around a single rotation axis. However,
for some objects, projection images from different rotation axes may contain
complementary information. In such cases, more accurate reconstructions could
in principle be achieved by collecting projection images from multiple
rotation axes (or poses), and then performing a joint or multi-pose
reconstruction.

Multi-pose reconstruction is motivated by samples with dense materials such as
metal, which produce substantial, spectrally-dependent attenuation in the
projection data. Standard reconstructions from such data exhibit bright or
dark streaks, commonly known as metal artifacts, radiating out from the metal
object [4]. These metal artifacts are strongly influenced by the measurement
pose of the object [4, 5, 6, 7], and specific regions of the object may be
corrupted in one pose but artifact-free in another. Therefore, multi-pose
reconstruction holds the potential to mitigate metal artifacts by effectively
taking the best quality information from each pose.

One approach to multi-pose reconstruction is to transform separate
reconstructions to a common pose and fuse them using a post-processing
technique such as averaging [8], pixel-wise minimum/maximum [9], or weighted
sum [10, 11]. However, this approach operates in the image domain and hence
does not fully exploit the complementary information in sinogram measurements.
In contrast, Kostenko et al. [12] and Herl et al. [13] adapt the scanner
geometry of multiple poses to a common pose to produce a joint reconstruction
from sinogram data. However, this method is complex to implement, since it
requires reprogramming the system matrix and integrating pose information into
the reconstruction software.

Metal artifacts can also be reduced by the use of Model-based Iterative
Reconstruction (MBIR) [14, 15, 16, 17, 18, 19, 20], which leverages the
regularity within reconstructions to compensate for noisy or limited
measurements. More recently, MBIR has evolved to include Plug-and-Play (PnP)
[21, 22], which uses a denoiser to model the prior distribution of an image,
and Multi-Agent Consensus Equilibrium (MACE), which allows multiple agents to
represent various objectives in an inverse problem [23, 24, 25, 26, 27]. One
notable advantage of MACE is its formulation as an equilibrium problem, making
it applicable even in cases where there is no well-defined cost function to
minimize. In this context, Yang et al. introduced a Multi-pose Fusion approach
to perform a joint tomographic reconstruction using multiple poses of a single
object [27].

Figure 1: Multi-pose Fusion Overview: Multiple sets of CT measurements are
acquired from different poses of the same object. The measurement data from
different poses are then fused in a MACE framework to form a joint
reconstruction.

In this paper, we build on the prior work of [27], and propose a pixel-
weighted Multi-pose Fusion algorithm, which reduces metal artifacts by
integrating CT measurements from multiple poses of a single object (see Figure
1). As illustrated in Fig. 1, the proposed algorithm uses MACE to integrate
information from the various poses, with each pose being represented by a
single agent. The final reconstruction is then computed as the fixed point of
a set of equilibrium equations. This allows for a straightforward, modular
implementation using standard CT reconstruction software. A primary novelty of
the proposed method is that the output of different MACE agents are fused in a
pixel-weighted manner to minimize the effects of metal throughout the
reconstruction. We present experimental results for the problem of non-
destructive evaluation (NDE) using measured CT data, and demonstrate that
pixel-weighted Multi-pose Fusion method is effective in reducing metal
artifacts and improving image quality.

##  II Problem Formulation

In multi-pose CT imaging, multiple sets of CT scans are taken using different
poses of the object, as illustrated in Fig. 1. Notice that the imaging
geometry is the same for each pose. However, in practice some poses provide
more useful information, particularly when the object contains dense or even
opaque components that may obscure portions of the object. The objective of
multi-pose CT reconstruction is then to perform a joint MBIR reconstruction
from scans acquired from multiple poses of the object.

Let $y_{k}\in\mathbb{R}^{M_{k}}$ be the sinogram measurements for the $k^{th}$
pose, where $k\in\left\\{0,...,K-1\right\\}$. Then our goal is to recover
$x\in\mathbb{R}^{N}$, the image vector containing attenuation coefficients in
the reconstruction coordinate system.

For each pose $k$, we also define a transformation function $x_{k}=T_{k}x$
where $x$ is the object represented in the common reconstruction coordinate
system and $x_{k}$ is the object represented in the $k^{th}$ pose. So
intuitively, $T_{k}$ transforms the raster sampled object from the common
reconstruction coordinate system to the posed coordinate system. In practice,
$T_{k}$ typically implements a rigid body transformation [28], so it requires
that the discretized function be resampled on the transformed sampling grid.
This process requires some form of typically spline-based interpolation
algorithm [29]. We will also require an approximate inverse transformation
$T_{k}^{-1}$. Since both transforms require interpolation and resampling, we
note that they will not in general be exact inverses of each other.

Using this notation, the forward model for each pose has the form

$y_{k}=A_{k}T_{k}x+w_{k}\ ,$ (1)

where $A_{k}\in\mathbb{R}^{N\times M_{k}}$ is the scanner system matrix and
$w_{k}\sim N(0,\alpha\Lambda_{k}^{-1})$ is independent additive noise, each
for the $k^{th}$ pose.

The joint MBIR reconstruction for the multi-pose problem is then given by

$x^{*}=\arg\min_{x}\left\\{\sum_{k=0}^{K-1}f_{k}(x)+h(x)\right\\}$ (2)

with data fidelity terms given by $f_{k}(x)=-\log p(y_{k}|x)+\mbox{const}$
where

$f_{k}(x)=\frac{1}{2}\lVert y_{k}-A_{k}T_{k}x\rVert_{\Lambda_{k}}^{2}\ ,$ (3)

and a prior term given by $h(x)=-\log p(x)$ that imposes regularity.

Notice that direct implementation of (2) is difficult since it requires that
software be written to minimize a sum of complex tomographic reconstruction
terms each with a different transformation $T_{k}$.

Alternatively, one can compute the MBIR reconstruction by using consensus ADMM
[30] to minimize the sum of $K+1$ terms consisting of $h$ and the $K$ terms in
(3). However, this approach has two serious disadvantages. First, the proximal
map terms required for each pose will be very computationally expensive to
compute. Second, we can improve reconstruction quality by replacing the prior
term $h(x)$ with a PnP denoiser.

##  III MACE formulation of Pixel-weighted Multi-pose Fusion

In this section, we introduce the MACE framework for solving the multi-pose
reconstruction problem [23, 26].

###  III-A Agent formulation

We start by describing the agents in our approach. Let $x^{\prime}=F_{k}(x)$
be an agent for the $k^{th}$ pose. Intuitively, this agent should take a
reconstruction $x$ and return a reconstruction $x^{\prime}$ that better fits
the measurements $y_{k}$ associated with the data from pose $k$. One approach
would be to use (3) directly in a proximal map. However, this would be
computationally expensive and difficult to compute since it requires that the
transformation $T_{k}$ be integrated into the reconstruction software.

Alternatively, we propose to use a _Conjugate Proximal Map_ given by

$F_{k}(v)=T_{k}^{-1}F(T_{k}v;y_{k})\ ,$ (4)

where $F(v;y)$ is the standard proximal map in reconstruction coordinates
given by

$F(v;y)=\arg\min_{x}\left\\{\frac{1}{2}\lVert
y-A_{k}x\rVert_{\Lambda_{k}}^{2}+\frac{1}{2\sigma^{2}}\lVert
x-v\rVert^{2}\right\\}.$ (5)

Notice that the conjugate proximal map of (4) can be computed easily since it
requires only the computation of the standard proximal map of (5) in the
standard coordinates and pre- and post-composition with the spline-based maps
$T_{k}$ and $T_{k}^{-1}$. In fact, software for computing the proximal map in
(5) is openly available [31, 32].

The conjugate proximal map is exactly equivalent to the conventional proximal
map based on (3) when the rigid body transformations $T_{k}$ and $T_{k}^{-1}$
are exact inverses of each other (see Appendix A). In practice, we note that
$T_{k}$ and $T_{k}^{-1}$ will generally not be exact inverses due to the
nonlinear interpolation, hence the solution to the conjugate proximal map is
an approximation to the conventional proximal map.

For the prior model, we will use a BM4D denoiser [33]. We denote this agent by
$F_{K}$.

The MACE agents are concatenated together to form a single operator ${\bf
F}:\mathbb{R}^{(K+1)N}\rightarrow\mathbb{R}^{(K+1)N}$, defined as:

${\bf F}({\bf x})=\left[F_{0}(x_{0}),\cdots,F_{K}(x_{K})\right],$ (6)

where ${\bf x}=\left[x_{0},...,x_{K}\right]$ denotes the full MACE state.

###  III-B Pixel-weighted MACE

The MACE state vector ${\bf x}$ contains multiple, potentially inconsistent
reconstructions $x_{0},...,x_{K}$. In order to produce a single coherent MACE
reconstruction, we define a pixel-weighted averaging operator

${\bf G_{M}}({\bf x})=\left[\bar{x}_{\bf M}({\bf x}),\cdots,\bar{x}_{\bf
M}({\bf x})\right]\ ,$ (7)

where $\bar{x}_{\bf M}({\bf x})$ is a pixel-weighted average of the input
vector components given by

$\bar{x}_{\bf M}({\bf
x})=\frac{1}{1+\beta}\sum_{k=0}^{K-1}M_{k}x_{k}+\frac{\beta}{1+\beta}x_{K}.$
(8)

Here, $\beta>0$ controls the amount of regularization relative to the data-
fitting agents, and $M_{k}\in\mathbb{R}^{N\times N}$ is a diagonal weight
matrix specific to each data-fitting agent satisfying the property

$\begin{cases}\left[M_{k}\right]_{ii}\geq 0\\\ \sum_{k=0}^{K-1}M_{k}=I\\\
\end{cases}\ .$

Intuitively, ${\bf G}_{\bf M}({\bf x})$ computes a weighted average of the
components in ${\bf x}$ according to the weight matrices $M_{0},...,M_{K-1}$
as well as the regularization parameter $\beta$, and then returns a state
vector formed by replicating this average $K+1$ times.

Notice that $M_{k}$ provides a mechanism to weight each pixel in each pose
separately, which could be valuable in the metal artifact scenario, where
certain components from specific poses are corrupted. The design of this
weight matrix is discussed in Section IV.

Using this notation, the MACE equilibrium equation is

${\bf F}({\bf x^{*}})={\bf G_{M}}({\bf x^{*}})\,$ (9)

where ${\bf x^{*}}$ solves the equation, and the final reconstruction is then
given by $x^{*}=\bar{x}_{\bf M}({\bf x^{*}})$. This equation enforces that all
agents have the same output value (consensus) and that the vectors
$\delta^{*}_{k}=x^{*}_{k}-F_{k}(x^{*}_{k})$ satisfy $\bar{x}_{\bf
M}({\bf\delta})=0$ (equilibrium) [23].

###  III-C Computing the MACE Solution

It can be shown that the solution to (9) is also the fixed point of the
operator ${\bf T}=(2{\bf G_{M}}-I)(2{\bf F}-I)$ (see Appendix B). One popular
method of finding such a fixed point is Mann iteration

${\bf w}\leftarrow(1-\rho){\bf w}+\rho{\bf T}{\bf w},$ (10)

where $\rho\in(0,1)$ controls the convergence speed.

Algorithm 1 Pixel-weighted MACE algorithm

Input: Initial Reconstruction: $x^{(0)}\in\mathbb{R}^{N}$
Input: Weight matrices $M_{0},...,M_{K}\in\mathbb{R}^{N\times N}$
Output: Final Reconstruction: $x^{*}\in\mathbb{R}^{N}$

1:${\bf w}\leftarrow[x^{(0)},...,x^{(0)}]$

2:while not converged do

3: ${\bf x}\leftarrow{\bf F}({\bf w})$

4: ${\bf z}\leftarrow{\bf G_{M}}(2{\bf x}-{\bf w})$

5: ${\bf w}\leftarrow{\bf w}+2\rho({\bf z}-{\bf x})$

6:return $x^{*}\leftarrow\bar{x}_{\bf M}({\bf x})$

Algorithm 1 shows the general method of solving MACE with Mann iterations. The
algorithm starts from an initial reconstruction $x^{(0)}$, and uses Mann
iterations to find the equilibrium point between the prior and forward model
terms. From [23], when the agents $F_{0},...,F_{K}$ are all proximal maps of
associated cost functions $f_{0},...,f_{K-1},h$, and all agents are equally
weighted ($M_{k}=\frac{1}{K+1}I$ for $k\in\left\\{0,...,K-1\right\\}$), this
equilibrium point is exactly the solution to the consensus optimization
problem of (2).

##  IV Pixel-weighted Averaging for Metal Artifact Reduction

We leverage the weight matrix $M_{k}$ in (8), which provides a mechanism of
applying a separate weight to each pixel in each pose, to propose a pixel-
weighted averaging algorithm for metal artifact reduction.

Given an initial reconstruction $x^{(0)}$ in the reconstruction coordinates,
we locate the metal and object components through the binary masks
$b^{metal},b^{object}\in\mathbb{R}^{N}$, where 1 indicates a metal/object
pixel, and 0 elsewhere:

$\left[b^{metal}\right]_{i}=\begin{cases}1,&\left[x^{(0)}\right]_{i}>\tau_{metal}\\\
0,&\left[x^{(0)}\right]_{i}\leq\tau_{metal}\\\ \end{cases}$ (11)

$\left[b^{object}\right]_{i}=\begin{cases}1,&\left[x^{(0)}\right]_{i}>\tau_{object}\\\
0,&\left[x^{(0)}\right]_{i}\leq\tau_{object}\\\ \end{cases}$ (12)

where $\tau_{metal},\tau_{object}\in\mathbb{R}$ are threshold values to
identify the metal and object pixels. We then transform the image masks from
the reconstruction pose to each of the measurement poses
$k\in\left\\{0,...,K-1\right\\}$:

$b_{k}^{metal}=T_{k}b^{metal}$ (13) $b_{k}^{object}=T_{k}b^{object}$ (14)

For each pose $k$, we compute a distortion image $D_{k}\in\mathbb{R}^{N}$,
which estimates the level of distortion at each pixel in the associated
measurement pose:

$\left[D_{k}\right]_{i}=\frac{\left[A_{k}^{t}A_{k}b_{k}^{metal}\right]_{i}}{\left[A_{k}^{t}A_{k}b_{k}^{object}\right]_{i}+\epsilon},$
(15)

where small $\epsilon>0$ prevents division by zero. Here, $D_{k}$ uses both
the location of the metal components ($b_{k}^{metal}$) and the scanner
geometry ($A_{k}$) to predict the pixel-wise level of distortion. Intuitively,
a larger entry in $D_{k}$ indicates that the associated projections traverse a
longer distance through metal, which could lead to more corruptions in the
reconstructed pixel.

The weight matrices $M_{0},...,M_{K-1}$ are then calculated using a softmax
function across all distortion images:

$\left[M_{k}\right]_{ii}=\frac{\exp{-\alpha\left[T_{k}^{-1}D_{k}\right]_{i}}}{\sum_{m=0}^{K-1}\exp{-\alpha\left[T_{m}^{-1}D_{m}\right]_{i}}},$
(16)

where $T_{k}^{-1}$ transforms the distortion image from the measurement pose
to the common reconstruction pose, and $\alpha\geq 0$ controls the range of
the weight matrices.

Subsequently, pixel-weighted Multi-pose Fusion (Algorithm 1) uses these weight
matrices to produce a joint reconstruction that selectively fits the
informative measurements from different poses.

Notice that the proposed algorithm can also work as a post-processing method
to directly integrate the images from different poses. We call this method the
pixel-weighted post-processing. The pseudo-code for this method is depicted in
Algorithm 2. The algorithm takes standard CT reconstructions
$x_{0},...,x_{K-1}$ as inputs, each from a distinct pose. For each pose, the
distortion image $D_{k}$ and subsequently the weight matrix $M_{k}$ are
computed. The post-processed image is formed by taking the pixel-weighted
average among the input reconstructions $x_{0},...,x_{K-1}$.

Algorithm 2 Pixel-weighted Post-processing for Metal Artifact Reduction

Input: Standard reconstructions: $x_{0},...,x_{K-1}\in\mathbb{R}^{N}$
Input: Transformation mappings: $T_{0},...,T_{K-1}$
Output: Pixel-weighted average: $x^{*}\in\mathbb{R}^{N}$

1:for $k\leftarrow 0$ to $K-1$ do

2: Compute $b_{k}^{metal}$, $b_{k}^{object}$ $\triangleright$ image masks Eq.
(13), (14)

3: Compute $D_{k}$ $\triangleright$ distortion image, Eq. (15)

4: Compute $M_{k}$ $\triangleright$ weight matrix, Eq. (16)

5:return $x^{*}=\sum_{k=0}^{K-1}M_{k}T_{k}^{-1}x_{k}$

##  V Experimental Results

We evaluate the effectiveness of pixel-weighted Multi-pose Fusion on a real
cone beam CT dataset featuring two distinct measurement poses.

Figure 2: Demonstration of the object of interest. The object is composed of a plastic component with four removable metal disks. (a) A photo of the object of interest. (b) and (c) vertical and horizontal slices without the metal components. In the scans below, metal disks are inserted in the regions indicated by the red arrows.  Figure 3: 3D Visualization of the measurement poses. The CT scans are acquired from two measurement poses: (a) vertical pose. (b) horizontal pose.  TABLE I: Experimental Specifications for Plastic Part with Metal Disk Dataset labels |  Vertical Pose Scan |  Horizontal Pose Scan
---|---|---
Source-Detector Distance |  569.363 mm |  580.247 mm
Magnification |  3.87 |  3.47
Number of Views |  900 |  1150
Detector Array Size |  960$\times$768 |  768$\times$960
Detector Pixel Pitch |  0.254 mm |  0.254 mm

The object of interest (Fig. 2) is a plastic component with four removable
metal disks. The object is scanned from two different poses (Fig. 3) using a
North Star Imaging X50 X-ray CT system. The experimental specifications are
detailed in Table I.

Figure 4: Comparison of different reconstruction methods and different poses.
(a) PnP, vertical pose. (b) PnP, horizontal pose. (c) Averaging of a and b.
(d) Pixel-weighted averaging of a and b. (e) Baseline Multi-pose Fusion. (f)
Pixel-weighted Multi-pose Fusion (proposed).
Pixel-weighted Multi-pose Fusion produces the best image quality with
significantly reduced metal artifacts.

We perform a pixel-weighted Multi-pose Fusion with Algorithm 1, and compare
the results with various reconstruction and post-processing methods listed
below:

  1. (a)

PnP, vertical pose: PnP recon from vertical pose.

  2. (b)

PnP, horizontal pose: PnP recon from horizontal pose.

  3. (c)

Averaging: Averaging of the PnP results.

  4. (d)

Pixel-weighted averaging: Pixel-weighted averaging of the PnP results.

  5. (e)

Baseline MPF: Baseline Multi-pose Fusion. The agents are equally weighted.

  6. (f)

Pixel-weighted MPF (proposed): Pixel-weighted Multi-pose Fusion. A separate
weight is assigned to each pixel in each agent.

For a fair comparison, we use the same BM4D denoiser in PnP and Multi-pose
Fusion algorithms.

Fig. 4 shows the results with different methods. The comparison of PnP results
reveals varying metal artifact characteristics across different poses. The
artifacts mainly manifest as horizontal streaks in the vertical pose (in red
boxes), while appearing as shadowy artifacts between adjacent metal disks in
the horizontal pose (in yellow boxes). Compared to simple averaging, pixel-
weighted averaging further reduces metal artifacts by selectively integrating
the informative components from different poses. Pixel-weighted Multi-pose
Fusion further reduces metal artifacts by incorporating this pixel-weighted
averaging mechanism into MACE framework, which produces the best image quality
with significantly reduced metal artifacts.

##  VI Conclusion

In this paper, we have introduced a novel CT reconstruction method called
pixel-weighted Multi-pose Fusion. Our method uses Multi-Agent Consensus
Equilibrium (MACE), an extension of Plug-and-Play, as a framework for
integrating projection data from different poses. A primary novelty of the
proposed method is that the output of different MACE agents are fused in a
pixel-weighted manner to minimize the effects of metal throughout the
reconstruction. Our experiment demonstrated that pixel-weighted Multi-pose
Fusion delivers a significant reduction in metal artifacts compared to single-
pose reconstruction and post-processing methods.

## Appendix

###  VI-A Relationship between Conventional Proximal Map and Conjugate
Proximal Map

The conjugate proximal map for the $k^{th}$ pose is given by

$F_{k}(v)=T_{k}^{-1}\arg\min_{x}\left\\{\frac{1}{2}\lVert
y_{k}-Ax\rVert_{\Lambda_{k}}^{2}+\frac{1}{2\sigma^{2}}\lVert
x-T_{k}v\rVert^{2}\right\\}.$ (17)

We show that the conjugate proximal map (17) is exactly equivalent to the
conventional proximal map based on (3) when the transformation mappings
$T_{k},T_{k}^{-1}$ satisfy the following conditions for all
$x\in\mathbb{R}^{N}$:

  1. (a)

$x=T_{k}^{-1}T_{k}x$ (inverse condition).

  2. (b)

$\lVert x\rVert=\lVert T_{k}x\rVert$ (rigid body transformation)

With a change of variable $x=T_{k}x_{k}$, we may rewrite the conjugate
proximal map (17) as

$\displaystyle F_{k}(v)$
$\displaystyle=T_{k}^{-1}\arg\min_{T_{k}x_{k}}\left\\{\frac{1}{2}\lVert
y_{k}-AT_{k}x_{k}\rVert_{\Lambda_{k}}^{2}+\frac{1}{2\sigma^{2}}\lVert
T_{k}x_{k}-T_{k}v\rVert^{2}\right\\}$ (18)
$\displaystyle=\arg\min_{x_{k}}\left\\{\frac{1}{2}\lVert
y_{k}-AT_{k}x_{k}\rVert_{\Lambda_{k}}^{2}+\frac{1}{2\sigma^{2}}\lVert
T_{k}(x_{k}-v)\rVert^{2}\right\\}$ (19)
$\displaystyle=\arg\min_{x_{k}}\left\\{\frac{1}{2}\lVert
y_{k}-AT_{k}x_{k}\rVert_{\Lambda_{k}}^{2}+\frac{1}{2\sigma^{2}}\lVert
x_{k}-v\rVert^{2}\right\\},$ (20)

where (19) requires the inverse condition (a), and (20) requires the rigid
body transformation condition (b).

Notice that after substituting the dummy variable $x_{k}$ with $x$ in (20),
the conjugate proximal map is exactly equivalent to the conventional proximal
map:

$\tilde{F}_{k}(v)=\arg\min_{x}\left\\{\frac{1}{2}\lVert
y-A_{k}T_{k}x\rVert_{\Lambda_{k}}^{2}+\frac{1}{2\sigma^{2}}\lVert
x-v\rVert^{2}\right\\}$ (21)

###  VI-B Reformulating MACE as a Fixed-point Problem

We follow [23] and show that MACE equation (9) can be reformulated as a fixed
point problem.

For notation simplicity, we rewrite the weight matrices as

$M_{k}^{\prime}=\begin{cases}\frac{1}{1+\beta}M_{k},&0\leq k\leq K-1\\\
\frac{\beta}{1+\beta}I,&k=K\end{cases}\ .$

With this notation, the weighted averaging operator (8) can be rewritten as:

$\bar{x}_{\bf M}({\bf x})=\sum_{k=0}^{K}M_{k}^{\prime}x_{k}.$ (22)

By definitions of $\bar{x}_{\bf M}$ and ${\bf G}_{\bf M}$, we have

$\displaystyle\bar{x}_{\bf M}({\bf G}_{\bf M}({\bf x}))$
$\displaystyle=\sum_{k=0}^{K}M_{k}^{\prime}\left[{\bf G}_{\bf M}({\bf
x})\right]_{k}$ (23) $\displaystyle=\sum_{k=0}^{K}M_{k}^{\prime}\bar{x}_{\bf
M}({\bf x})$ (24) $\displaystyle=\bar{x}_{\bf M}({\bf x}),$ (25)

where (25) holds because $\sum_{k=0}^{K}M_{k}^{\prime}=I$. This gives the
following property:

${\bf G}_{\bf M}({\bf G}_{\bf M}({\bf x}))={\bf G}_{\bf M}({\bf x})$ (26)

With this linear property, we may further show that ${\bf G}_{\bf M}$
satisfies the property

$(2{\bf G}_{\bf M}-{\bf I})^{-1}=(2{\bf G}_{\bf M}-{\bf I}).$ (27)

From this, we may reformulate the MACE equation (9) as a fixed point problem:

$(2{\bf G}_{\bf M}-{\bf I})(2{\bf F}-{\bf I})({\bf x^{*}})={\bf x^{*}},$ (28)

or ${\bf Tx^{*}}={\bf x^{*}}$, where ${\bf T}=(2{\bf G}_{\bf M}-{\bf I})(2{\bf
F}-{\bf I})$.

## References

  * [1] W. Sun, S. B. Brown, and R. K. Leach, “An overview of industrial x-ray computed tomography.” NPL Report, January 2012.
  * [2] H. P. Hiriyannaiah, “X-ray computed tomography for medical imaging,” _IEEE signal Processing magazine_ , vol. 14, no. 2, pp. 42–59, 1997.
  * [3] P. J. Withers, C. Bouman, S. Carmignato, V. Cnudde, D. Grimaldi, C. K. Hagen, E. Maire, M. Manley, A. Du Plessis, and S. R. Stock, “X-ray computed tomography,” _Nature Reviews Methods Primers_ , vol. 1, no. 1, p. 18, 2021.
  * [4] L. Gjesteby, B. De Man, Y. Jin, H. Paganetti, J. Verburg, D. Giantsoudi, and G. Wang, “Metal artifact reduction in CT: where are we after four decades?” _Ieee Access_ , vol. 4, pp. 5826–5849, 2016.
  * [5] J. H. Brown, E. Lustrin, M. H. Lev, C. S. Ogilvy, and J. M. Taveras, “Reduction of aneurysm clip artifacts on CT angiograms: a technical note.” _AJNR. American journal of neuroradiology_ , vol. 20 4, pp. 694–6, 1999.
  * [6] M. Lewis, A. Toms, K. Reid, and W. Bugg, “CT metal artefact reduction of total knee prostheses using angled gantry multiplanar reformation.” _The Knee_ , vol. 17 4, pp. 279–82, 2010.
  * [7] M. Luckow, H. Deyhle, F. Beckmann, D. Dagassan-Berndt, and B. Müller, “Tilting the jaw to improve the image quality or to reduce the dose in cone-beam computed tomography,” _European journal of radiology_ , vol. 80, no. 3, pp. e389–e393, 2011.
  * [8] H. Ballhausen, M. Reiner, U. Ganswindt, C. Belka, and M. Söhn, “Post-processing sets of tilted ct volumes as a method for metal artifact reduction,” _Radiation Oncology_ , vol. 9, no. 1, pp. 1–10, 2014.
  * [9] J. Guhathakurta, S. Kiess, J. Hillebrand, J. Denecke, I. Effenberger, and S. Simon, “Reducing computed tomography reconstruction and beam hardening artefacts by data fusion,” in _International Symposium on Digital Industrial Radiology and Computed Tomography (DIR)_ , 2015, pp. 330–337.
  * [10] K. Kobayashi, Y. Ohtake, Y. Nagai, and H. Suzuki, “Multi-pass ct scanning with orientation selection and volume integration for metal artifact reduction,” in _15th Asia Pacific Conference for Non-Destructive Testing (APCNDT2017)_ , 2017.
  * [11] G. Herl, S. Rettenberger, J. Hiller, and T. Sauer, “Metal artifact reduction by fusion of ct scans from different positions using the unfiltered backprojection,” in _Proceedings of 8th Conference on Industrial Computed Tomography_ , 2018, pp. 1–8.
  * [12] A. Kostenko, V. Andriiashen, and K. J. Batenburg, “Registration-based multi-orientation tomography,” _Optics express_ , vol. 26, no. 22, pp. 28 982–28 995, 2018.
  * [13] G. Herl, J. Hiller, and T. Sauer, “Artifact reduction in x-ray computed tomography by multipositional data fusion using local image quality measures,” in _Proceedings of the 9th Conference on Industrial Computed Tomography (iCT 2019), Padova, Italy_ , 2019, pp. 13–15.
  * [14] C. Bouman and K. Sauer, “A generalized gaussian image model for edge-preserving MAP estimation,” _IEEE Transactions on Image Processing_ , vol. 2, no. 3, pp. 296–310, 1993.
  * [15] J.-B. Thibault, K. D. Sauer, C. A. Bouman, and J. Hsieh, “A three-dimensional statistical approach to improved image quality for multislice helical CT,” _Medical physics_ , vol. 34, no. 11, pp. 4526–4544, 2007.
  * [16] R. Zhang, J.-B. Thibault, C. A. Bouman, K. D. Sauer, and J. Hsieh, “Model-based iterative reconstruction for dual-energy X-ray CT using a joint quadratic likelihood model,” _IEEE transactions on medical imaging_ , vol. 33, no. 1, pp. 117–134, 2013.
  * [17] P. Jin, C. A. Bouman, and K. D. Sauer, “A model-based image reconstruction algorithm with simultaneous beam hardening correction for X-ray CT,” _IEEE Transactions on Computational Imaging_ , vol. 1, no. 3, pp. 200–216, 2015.
  * [18] S. J. Kisner, E. Haneda, C. A. Bouman, S. Skatter, M. Kourinny, and S. Bedford, “Model-based CT reconstruction from sparse views,” in _Second International Conference on Image Formation in X-Ray Computed Tomography_ , 2012, pp. 444–447.
  * [19] S. J. Kisner, P. Jin, C. A. Bouman, K. Sauer, W. Garms, T. Gable, S. Oh, M. Merzbacher, and S. Skatter, “Innovative data weighting for iterative reconstruction in a helical ct security baggage scanner,” in _2013 47th International Carnahan Conference on Security Technology (ICCST)_. IEEE, 2013, pp. 1–5.
  * [20] T. Balke, S. Majee, G. T. Buzzard, S. Poveromo, P. Howard, M. A. Groeber, J. McClure, and C. A. Bouman, “Separable models for cone-beam mbir reconstruction,” _electronic imaging_ , vol. 30, pp. 1–1817, 2018.
  * [21] S. V. Venkatakrishnan, C. A. Bouman, and B. Wohlberg, “Plug-and-play priors for model based reconstruction,” in _2013 IEEE Global Conference on Signal and Information Processing_. IEEE, 2013, pp. 945–948.
  * [22] S. Sreehari, S. V. Venkatakrishnan, B. Wohlberg, G. T. Buzzard, L. F. Drummy, J. P. Simmons, and C. A. Bouman, “Plug-and-play priors for bright field electron tomography and sparse interpolation,” _IEEE Transactions on Computational Imaging_ , vol. 2, no. 4, pp. 408–423, 2016.
  * [23] G. T. Buzzard, S. H. Chan, S. Sreehari, and C. A. Bouman, “Plug-and-play unplugged: Optimization-free reconstruction using consensus equilibrium,” _SIAM Journal on Imaging Sciences_ , vol. 11, no. 3, pp. 2001–2020, 2018.
  * [24] V. Sridhar, X. Wang, G. T. Buzzard, and C. A. Bouman, “Distributed iterative ct reconstruction using multi-agent consensus equilibrium,” _IEEE Transactions on Computational Imaging_ , vol. 6, pp. 1153–1166, 2020.
  * [25] S. Majee, T. Balke, C. A. Kemp, G. T. Buzzard, and C. A. Bouman, “Multi-slice fusion for sparse-view and limited-angle 4D CT reconstruction,” _IEEE Transactions on Computational Imaging_ , vol. 7, pp. 448–461, 2021.
  * [26] C. A. Bouman, _Foundations of Computational Imaging: A Model-Based Approach_. Philadelphia: Society for Industrial and Applied Mathematics, 2022.
  * [27] D. Yang, C. A. Kemp, G. T. Buzzard, and C. A. Bouman, “Multi-pose fusion for sparse-view ct reconstruction using consensus equilibrium,” in _2022 58th Annual Allerton Conference on Communication, Control, and Computing (Allerton)_. IEEE, 2022, pp. 1–5.
  * [28] D. W. Eggert, A. Lorusso, and R. B. Fisher, “Estimating 3-d rigid body transformations: a comparison of four major algorithms,” _Machine vision and applications_ , vol. 9, no. 5, pp. 272–290, 1997.
  * [29] C. De Boor and C. De Boor, _A practical guide to splines_. springer-verlag New York, 1978, vol. 27.
  * [30] S. Boyd, N. Parikh, E. Chu, B. Peleato, J. Eckstein _et al._ , “Distributed optimization and statistical learning via the alternating direction method of multipliers,” _Foundations and Trends® in Machine learning_ , vol. 3, no. 1, pp. 1–122, 2011.
  * [31] SVMBIR Development Team, “Super-Voxel Model Based Iterative Reconstruction (SVMBIR),” Software library available from https://github.com/cabouman/svmbir, 2020.
  * [32] MBIRCONE Development Team, “Model-Based Iterative Reconstruction for cone-beam (MBIRCONE),” Software library available from https://github.com/cabouman/mbircone, 2020.
  * [33] M. Maggioni, V. Katkovnik, K. Egiazarian, and A. Foi, “Nonlocal transform-domain filter for volumetric data denoising and reconstruction,” _IEEE transactions on image processing_ , vol. 22, no. 1, pp. 119–133, 2012.