would be a substantial overcorrection for masks made of many tiny holes that
are partially re-filled by the optimal filtering. Similar results have been
found for all the other foreground masks.

Figure 18:  Relative differences between the reconstructed cross
$\hat{\kappa}\kappa$-spectrum and the true Websky $C_{L}^{\kappa\kappa}$ for
three different uncorrelated $W_{\rm halo}$ masks normalizing with $f_{\rm
sky}^{-1}$ to account for the mask area but no multiplicative Monte Carlo
correction. Dashed lines show the same relative differences without the
$f_{\rm sky}^{-1}$ normalization. Dash-dotted and dotted lines show the
theoretical expectations for the model of Eq. (25) with and without the
$f_{\rm sky}^{-1}$ normalization respectively. In the theoretical calculation,
we assumed that all the disks have the same radius equal the mean of the halo
sizes (in terms of $2\theta_{500,c}$) in each mass-limited sample. For small
holes, the $f_{\rm sky}^{-1}$ normalization over-corrects for the loss of
power as discussed in Appendix A.

##  Appendix B Effects of masking non-Gaussian fields

|  |
---|---|---

|  |
---|---|---

Figure 19:  Same as Fig. 4 for the sign-flipped masks $W_{-\rm CIB}$,
$W_{-y}$, and $W_{-\kappa}$ that remove minima of the foreground fields.
Simulation measurements are shown as data points, the semi-analytic Gaussian
theory predictions in solid and the fully-analytic Gaussian model in dashed
(evaluating the expectation values in Eq. (6) empirically from our simulations
and analytically respectively). The top and bottom rows show the masks
obtained by smoothing the foreground maps with a Gaussian beam of
$\theta_{1/2}=5.1^{\prime},1.7^{\prime}$ respectively.

In this Appendix, we explore the disagreement between the theory predictions
based on simple Gaussian models and the simulation results for the threshold
masks, $W_{\kappa}$, $W_{\rm CIB}$ and $W_{y}$, as shown in Fig. 4.

To test higher-order non-Gaussian effects, as in Paper I , we constructed
another set of foreground masks from the $y$, ${\rm CIB}$ and $\kappa$ fields
by inverting their sign prior to the thresholding operation, $W_{-\rm CIB}$,
$W_{-y}$ and $W_{-\kappa}$. This is equivalent to masking the minima instead
of the peaks of the foreground fields. This helps isolating higher-order
effects that involve a higher even power of the perturbations compared to the
expansion considered in our analytical model (where $X$ and $-X$ masks give
identical results). The results are shown in Fig. 19. We found a good
agreement only for $W_{-\kappa}$ while we recover at least a good qualitative
trend for $W_{-\rm CIB}$. However for $W_{y}$, the predictions instead lie
away from the analytic curves, consistent with strongly non-Gaussian higher-
order effects (clearly visible in the mask shape difference shown in Fig. 20)
not being captured in our Gaussian model. Here the semi-empirical results
suggest only a very small bias when masking $W_{-y}$.

Figure 20:  Websky $\kappa$ field masked with $W_{y}$ and $W_{-y}$ masks
obtained smoothing the foreground field with a $\theta_{1/2}=5.1^{\prime}$
Gaussian beam prior to thresholding. The masked fraction of the sky is the
same for both masks ($f_{\rm sky}=93.3\%$). The tSZ signal is highly non-
Gaussian and highly skewed, so that the mask covers large clusters in the
standard case, and a much larger number of small underdensities in the $-y$
case (for Gaussian fields the two masks would be statistically identical).

The theory predictions are based on full-sky spectra of non-Gaussian
foreground maps, and hence are sensitive to the peaks of the field inside the
sky areas that get masked. Such features are enhanced for non-Gaussian fields
compared to a purely Gaussian field. To test the impact of these non-Gaussian
features, and to understand if they are responsible for the disagreement with
the simulation results, we tried creating a Gaussianized version of the non-
Gaussian foreground fields. To do this, we remapped the values in each pixel
to the value that they would have in same the corresponding percentile in a
Gaussian distribution having the same variance. Ideally, the masked foreground
spectra after this procedure should be similar to the one obtained from non-
Gaussian maps, but the full-sky ones used in Eq. (6) may be rather different,
as now all the very bright extreme pixels have had their values scaled down.
In Fig. 21 we show the analytic prediction based on these Gaussianized maps.
We still do not find a perfect agreement between simulations and theory
predictions, but at least we are now able to reproduce the trends seen in the
simulations: in particular, the $W_{y}$ mask predictions now recover, at
intermediate scales, distinctly different curves for masks retaining different
$f_{\rm sky}$s. Moreover, as expected, the semi-analytic theory predictions
and the pure Gaussian model obtained using the full-sky power spectra from
Gaussianized maps are closer together (see dashed and solid lines).

|
---|---

|
---|---

Figure 21: Effect of LSS-correlated masking on the CMB lensing convergence
power spectrum $C_{L}^{\kappa\kappa}$ as fractional difference between the
$C_{L}^{\kappa\kappa}$ computed on the masked sky and its value computed on
the full sky. We compare a subset of the simulation results shown in Fig. 4
with alternative analytical predictions. The semi-analytical and analytical
predictions of the Gaussian model obtained using the full-sky power spectra of
Gaussianized foreground maps are shown in solid and in dashed, respectively.
The baseline theory predictions (both semi-analytical and fully analytical)
based on non-Gaussian foreground maps and shown in Fig. 4 appear as faded
lines. Top and bottom panels show results for the $W_{\rm CIB}$ and $W_{y}$
masks constructed on foregrounds field smoothed with different Gaussian beams.

Finally, to test the accuracy of the analytical model and to convince
ourselves that discrepancies we observe between the analytical models and the
simulations are simply due to non-Gaussianity of the foreground fields, we
built correlated Gaussian map of the foreground fields, $y^{\rm G}$,
$\mathrm{CIB}^{\rm G}$ and $\kappa^{\rm G}$, that preserve the cross-
correlation between $y$, CIB and $\kappa$ as measured from the full-sky Websky
maps. We used these Gaussian realizations of the foreground fields to
construct the Gaussian foreground $W^{\rm G}_{y}$, $W^{\rm G}_{\rm CIB}$ and
$W^{\rm G}_{\kappa}$ masks that we later applied on the corresponding Gaussian
$\kappa^{\rm G}$ map. The results are shown in Fig. 22 and are in excellent
agreement with the analytic model described in Sec. II.2, confirming that the
differences we see in the more realistic simulations are due to the non-
Gaussian statistics.

|  |
---|---|---

|  |
---|---|---

Figure 22: Same as Fig. 19 for the $W^{\rm G}_{y}$, $W^{\rm G}_{\rm CIB}$ and
$W^{\rm G}_{\kappa}$ masks derived thresholding a Gaussian realization of the
foreground fields and of CMB lensing keeping the same correlation between the
components observed in the Websky simualations.

## References

  * Ade  _et al._ (2016a) P.A.R. Ade _et al._ (Planck), “Planck 2015 results. XXVII. The Second Planck Catalogue of Sunyaev-Zeldovich Sources,” Astron. Astrophys. 594, A27 (2016a), arXiv:1502.01598 [astro-ph.CO] .
  * Bleem  _et al._ (2020) L.E. Bleem _et al._ (SPT, DES), “The SPTpol Extended Cluster Survey,” Astrophys. J. Suppl. 247, 25 (2020), arXiv:1910.04121 [astro-ph.CO] .
  * Hilton  _et al._ (2021) M. Hilton _et al._ (ACT, DES), “The Atacama Cosmology Telescope: A Catalog of $>$4000 Sunyaev–Zel’dovich Galaxy Clusters,” Astrophys. J. Suppl. 253, 3 (2021), arXiv:2009.11043 [astro-ph.CO] .
  * Holder (2002) Gilbert P. Holder, “Radio point sources and the thermal sz power spectrum,” Astrophys. J. 580, 36–41 (2002), arXiv:astro-ph/0205467 .
  * Shirasaki (2019) Masato Shirasaki, “Impact of radio sources and cosmic infrared background on thermal Sunyaev–Zel’dovich – gravitational lensing cross-correlation,” Mon. Not. Roy. Astron. Soc. 483, 342–351 (2019), arXiv:1807.09412 [astro-ph.CO] .
  * Allison  _et al._ (2015) Rupert Allison _et al._ (ACT), “The Atacama Cosmology Telescope: measuring radio galaxy bias through cross-correlation with lensing,” Mon. Not. Roy. Astron. Soc. 451, 849–858 (2015), arXiv:1502.06456 [astro-ph.CO] .
  * Dwek and Barker (2002) Eli Dwek and Michael K. Barker, “The cosmic radio and infrared backgrounds connection,” Astrophys. J. 575, 7–11 (2002).
  * Ade  _et al._ (2014a) P.A.R. Ade _et al._ (Planck), “Planck 2013 results. XVIII. The gravitational lensing-infrared background correlation,” Astron. Astrophys. 571, A18 (2014a), arXiv:1303.5078 [astro-ph.CO] .
  * Ade  _et al._ (2016b) P.A.R. Ade _et al._ (Planck), “Planck 2015 results. XXIII. The thermal Sunyaev-Zeldovich effect–cosmic infrared background correlation,” Astron. Astrophys. 594, A23 (2016b), arXiv:1509.06555 [astro-ph.CO] .
  * Aghanim  _et al._ (2016) N. Aghanim _et al._ (Planck), “Planck 2015 results. XXII. A map of the thermal Sunyaev-Zeldovich effect,” Astron. Astrophys. 594, A22 (2016), arXiv:1502.01596 [astro-ph.CO] .
  * Song  _et al._ (2003) Yong-Seon Song, Asantha Cooray, Lloyd Knox, and Matias Zaldarriaga, “The far-infrared background correlation with cmb lensing,” Astrophys. J. 590, 664–672 (2003), astro-ph/0209001 .
  * Osborne  _et al._ (2014) Stephen J. Osborne, Duncan Hanson, and Olivier Doré, “Extragalactic Foreground Contamination in Temperature-based CMB Lens Reconstruction,” JCAP 1403, 024 (2014), arXiv:1310.7547 [astro-ph.CO] .
  * van Engelen  _et al._ (2014) A. van Engelen, S. Bhattacharya, N. Sehgal, G. P. Holder, O. Zahn, and D. Nagai, “CMB Lensing Power Spectrum Biases from Galaxies and Clusters using High-angular Resolution Temperature Maps,” Astrophys. J. 786, 13 (2014), arXiv:1310.7023 [astro-ph.CO] .
  * Namikawa  _et al._ (2013) Toshiya Namikawa, Duncan Hanson, and Ryuichi Takahashi, “Bias-hardened CMB lensing,” MNRAS 431, 609–620 (2013), arXiv:1209.0091 [astro-ph.CO] .
  * Schaan and Ferraro (2019) Emmanuel Schaan and Simone Ferraro, “Foreground-immune CMB lensing with shear-only reconstruction,” Phys. Rev. Lett. 122, 181301 (2019), arXiv:1804.06403 [astro-ph.CO] .
  * Darwish  _et al._ (2020) Omar Darwish _et al._ , “The Atacama Cosmology Telescope: A CMB lensing mass map over 2100 square degrees of sky and its cross-correlation with BOSS-CMASS galaxies,” MNRAS 500, 2250–2263 (2020), arXiv:2004.01139 [astro-ph.CO] .
  * Madhavacheril and Hill (2018) Mathew S. Madhavacheril and J. Colin Hill, “Mitigating Foreground Biases in CMB Lensing Reconstruction Using Cleaned Gradients,” Phys. Rev. D98, 023534 (2018), arXiv:1802.08230 [astro-ph.CO] .
  * Fabbian  _et al._ (2019) Giulio Fabbian, Antony Lewis, and Dominic Beck, “CMB lensing reconstruction biases in cross-correlation with large-scale structure probes,” JCAP 1910, 057 (2019), arXiv:1906.08760 [astro-ph.CO] .
  * Naess  _et al._ (2020) Sigurd Naess _et al._ , “The Atacama Cosmology Telescope: arcminute-resolution maps of 18 000 square degrees of the microwave sky from ACT 2008–2018 data combined with Planck,” JCAP 12, 046 (2020), arXiv:2007.07290 [astro-ph.IM] .
  * Fabbian  _et al._ (2021) G. Fabbian, J. Carron, A. Lewis, and M. Lembo, “Lensed CMB power spectrum biases from masking extragalactic sources,” Phys. Rev. D 103, 043535 (2021), arXiv:2011.08841 [astro-ph.CO] .
  * Aghanim  _et al._ (2020a) N. Aghanim _et al._ (Planck), “Planck 2018 results. VIII. Gravitational lensing,” A&A 641, A8 (2020a), arXiv:1807.06210 [astro-ph.CO] .
  * Hu and Okamoto (2002) Wayne Hu and Takemi Okamoto, “Mass reconstruction with CMB polarization,” Astrophys. J. 574, 566–574 (2002), arXiv:astro-ph/0111606 [astro-ph] .
  * Kesden  _et al._ (2003) Michael Kesden, Asantha Cooray, and Marc Kamionkowski, “Lensing reconstruction with CMB temperature and polarization,” Phys. Rev. D 67, 123507 (2003), astro-ph/0302536 .
  * Böhm  _et al._ (2016) Vanessa Böhm, Marcel Schmittfull, and Blake D. Sherwin, “Bias to CMB lensing measurements from the bispectrum of large-scale structure,” Phys. Rev. D 94, 043519 (2016), arXiv:1605.01392 [astro-ph.CO] .
  * Pratten and Lewis (2016) Geraint Pratten and Antony Lewis, “Impact of post-Born lensing on the CMB,” JCAP 1608, 047 (2016), arXiv:1605.05662 [astro-ph.CO] .
  * Böhm  _et al._ (2018) Vanessa Böhm, Blake D. Sherwin, Jia Liu, J. Colin Hill, Marcel Schmittfull, and Toshiya Namikawa, “Effect of non-Gaussian lensing deflections on CMB lensing measurements,” Phys. Rev. D98, 123510 (2018), arXiv:1806.01157 [astro-ph.CO] .
  * Beck  _et al._ (2018) Dominic Beck, Giulio Fabbian, and Josquin Errard, “Lensing Reconstruction in Post-Born Cosmic Microwave Background Weak Lensing,” Phys. Rev. D98, 043512 (2018), arXiv:1806.01216 [astro-ph.CO] .
  * Liu  _et al._ (2016) Jia Liu, J. Colin Hill, Blake D. Sherwin, Andrea Petri, Vanessa Böhm, and Zoltán Haiman, “CMB lensing beyond the power spectrum: Cosmological constraints from the one-point probability distribution function and peak counts,” Phys. Rev. D 94, 103501 (2016), arXiv:1608.03169 [astro-ph.CO] .
  * Vallinotto (2012) Alberto Vallinotto, “Using CMB lensing to constrain the multiplicative bias of cosmic shear,” Astrophys. J. 759, 32 (2012), arXiv:1110.5339 [astro-ph.CO] .
  * Schaan  _et al._ (2017) Emmanuel Schaan, Elisabeth Krause, Tim Eifler, Olivier Doré, Hironao Miyatake, Jason Rhodes, and David N. Spergel, “Looking through the same lens: Shear calibration for LSST, Euclid, and WFIRST with stage 4 CMB lensing,” Phys. Rev. D 95, 123512 (2017), arXiv:1607.01761 [astro-ph.CO] .
  * Cawthon (2020) Ross Cawthon, “Effects of Redshift Uncertainty on Cross-Correlations of CMB Lensing and Galaxy Surveys,” Phys. Rev. D 101, 063509 (2020), arXiv:1809.09251 [astro-ph.CO] .
  * Schaan  _et al._ (2020) Emmanuel Schaan, Simone Ferraro, and Uroš Seljak, “Photo-z outlier self-calibration in weak lensing surveys,” JCAP 12, 001 (2020), arXiv:2007.12795 [astro-ph.CO] .
  * Ilić  _et al._ (2021) S. Ilić _et al._ (Euclid), “$Euclid$ preparation: XV. Forecasting cosmological constraints for the $Euclid$ and CMB joint analysis,” (2021), arXiv:2106.08346 [astro-ph.CO] .
  * Sailer  _et al._ (2021) Noah Sailer, Emanuele Castorina, Simone Ferraro, and Martin White, “Cosmology at high redshift – a probe of fundamental physics,” (2021), arXiv:2106.09713 [astro-ph.CO] .
  * Hill and Spergel (2014) J.Colin Hill and David N. Spergel, “Detection of thermal SZ-CMB lensing cross-correlation in Planck nominal mission data,” JCAP 02, 030 (2014), arXiv:1312.4525 [astro-ph.CO] .
  * Battaglia  _et al._ (2015) N. Battaglia, J. C. Hill, and N. Murray, “Deconstructing Thermal Sunyaev–Zel’dovich—gravitational lensing cross-correlations: Implications for the Intracluster Medium,” Astrophys. J. 812, 154 (2015), arXiv:1412.5593 [astro-ph.CO] .
  * Hojjati  _et al._ (2015) Alireza Hojjati, Ian G. McCarthy, Joachim Harnois-Deraps, Yin-Zhe Ma, Ludovic Van Waerbeke, Gary Hinshaw, and Amandine M. C. Le Brun, “Dissecting the thermal Sunyaev-Zeldovich-gravitational lensing cross-correlation with hydrodynamical simulations,” JCAP 10, 047 (2015), arXiv:1412.6051 [astro-ph.CO] .
  * Van Waerbeke  _et al._ (2014) Ludovic Van Waerbeke, Gary Hinshaw, and Norman Murray, “Detection of warm and diffuse baryons in large scale structure from the cross-correlation of gravitational lensing and the thermal Sunyaev-Zeldovich effect,” Phys. Rev. D 89, 023508 (2014), arXiv:1310.5721 [astro-ph.CO] .
  * Osato  _et al._ (2018) Ken Osato, Samuel Flender, Daisuke Nagai, Masato Shirasaki, and Naoki Yoshida, “Investigating cluster astrophysics and cosmology with cross-correlation of the thermal Sunyaev–Zel’dovich effect and weak lensing,” Mon. Not. Roy. Astron. Soc. 475, 532–542 (2018), arXiv:1706.08972 [astro-ph.CO] .
  * Bethermin  _et al._ (2012) Matthieu Bethermin, Olivier Dore, and Guilaine Lagache, “Where stars form and live at high redshift: clues from the infrared,” Astron. Astrophys. 537, L5 (2012), arXiv:1201.0546 [astro-ph.CO] .
  * Maniyar  _et al._ (2021) A. Maniyar, M. Béthermin, and G. Lagache, “Simple halo model formalism for the cosmic infrared background and its correlation with the thermal Sunyaev Zel’dovich effect,” A&A 645, A40 (2021), arXiv:2006.16329 [astro-ph.CO] .
  * Ade  _et al._ (2014b) P. A. R. Ade _et al._ (Planck), “Planck 2013 results. XVIII. The gravitational lensing-infrared background correlation,” Astron. Astrophys. 571, A18 (2014b), arXiv:1303.5078 [astro-ph.CO] .
  * Aguirre  _et al._ (2019) James Aguirre _et al._ (Simons Observatory), “The Simons Observatory: Science goals and forecasts,” JCAP 1902, 056 (2019), arXiv:1808.07445 [astro-ph.CO] .
  * Aravena  _et al._ (2019) Manuel Aravena _et al._ , “The CCAT-Prime Submillimeter Observatory,” (2019), arXiv:1909.02587 [astro-ph.IM] .
  * McCarthy and Madhavacheril (2021) Fiona McCarthy and Mathew S. Madhavacheril, “Improving models of the cosmic infrared background using CMB lensing mass maps,” Phys. Rev. D 103, 103515 (2021), arXiv:2010.16405 [astro-ph.CO] .
  * Sherwin and Schmittfull (2015) Blake D. Sherwin and Marcel Schmittfull, “Delensing the CMB with the Cosmic Infrared Background,” Phys. Rev. D 92, 043005 (2015), arXiv:1502.05356 [astro-ph.CO] .
  * Yu  _et al._ (2017) Byeonghee Yu, J. Colin Hill, and Blake D. Sherwin, “Multitracer CMB delensing maps from Planck and WISE data,” Phys. Rev. D 96, 123511 (2017), arXiv:1705.02332 [astro-ph.CO] .
  * Beck  _et al._ (2020) Dominic Beck, Josquin Errard, and Radek Stompor, “Impact of Polarized Galactic Foreground Emission on CMB Lensing Reconstruction and Delensing of B-Modes,” JCAP 06, 030 (2020), arXiv:2001.02641 [astro-ph.CO] .
  * Abazajian  _et al._ (2019) Kevork Abazajian _et al._ , “CMB-S4 Science Case, Reference Design, and Project Plan,” (2019), arXiv:1907.04473 [astro-ph.IM] .
  * Cooray and Sheth (2002) Asantha Cooray and Ravi K. Sheth, “Halo Models of Large Scale Structure,” Phys. Rept. 372, 1–129 (2002), arXiv:astro-ph/0206508 .
  * Navarro  _et al._ (1997) J. F. Navarro, C. S. Frenk, and S. D. M. White, “A Universal Density Profile from Hierarchical Clustering,” Astrophys. J. 490, 493 (1997), astro-ph/9611107 .
  * Duffy  _et al._ (2008) Alan R. Duffy, Joop Schaye, Scott T. Kay, and Claudio Dalla Vecchia, “Dark matter halo concentrations in the Wilkinson Microwave Anisotropy Probe year 5 cosmology,” MNRAS 390, L64–L68 (2008), arXiv:0804.2486 [astro-ph] .
  * Sheth and Tormen (1999) Ravi K. Sheth and Giuseppe Tormen, “Large-scale bias and the peak background split,” MNRAS 308, 119–126 (1999), arXiv:astro-ph/9901122 [astro-ph] .
  * Rotti  _et al._ (2021) Aditya Rotti, Boris Bolliet, Jens Chluba, and Mathieu Remazeilles, “Removing the giants and learning from the crowd: a new SZ power spectrum method and revised Compton $y$-map analysis,” Mon. Not. Roy. Astron. Soc. 503, 5310–5328 (2021), arXiv:2010.07797 [astro-ph.CO] .
  * Sheth  _et al._ (2001) Ravi K. Sheth, H. J. Mo, and Giuseppe Tormen, “Ellipsoidal collapse and an improved model for the number and spatial distribution of dark matter haloes,” MNRAS 323, 1–12 (2001), arXiv:astro-ph/9907024 [astro-ph] .
  * Sheth and Tormen (2002) Ravi K. Sheth and Giuseppe Tormen, “An excursion set model of hierarchical clustering: ellipsoidal collapse and the moving barrier,” MNRAS 329, 61–75 (2002), arXiv:astro-ph/0105113 [astro-ph] .
  * Tinker  _et al._ (2010) Jeremy L. Tinker, Brant E. Robertson, Andrey V. Kravtsov, Anatoly Klypin, Michael S. Warren, Gustavo Yepes, and Stefan Gottlöber, “The Large-scale Bias of Dark Matter Halos: Numerical Calibration and Model Tests,” Astrophys. J. 724, 878–886 (2010), arXiv:1001.3162 [astro-ph.CO] .
  * Smith  _et al._ (2018) Kendrick M. Smith, Mathew S. Madhavacheril, Moritz Münchmeyer, Simone Ferraro, Utkarsh Giri, and Matthew C. Johnson, “KSZ tomography and the bispectrum,” (2018), arXiv:1810.13423 [astro-ph.CO] .
  * Stein  _et al._ (2020) George Stein, Marcelo A. Alvarez, J. Richard Bond, Alexander van Engelen, and Nicholas Battaglia, “The Websky Extragalactic CMB Simulations,” JCAP 10, 012 (2020), arXiv:2001.08787 [astro-ph.CO] .
  * Stein  _et al._ (2019) George Stein, Marcelo A. Alvarez, and J. Richard Bond, “The mass-Peak Patch algorithm for fast generation of deep all-sky dark matter halo catalogues and its N-Body validation,” Mon. Not. Roy. Astron. Soc. 483, 2236–2250 (2019), arXiv:1810.07727 [astro-ph.CO] .
  * Blot  _et al._ (2019) Linda Blot _et al._ , “Comparing approximate methods for mock catalogues and covariance matrices II: Power spectrum multipoles,” Mon. Not. Roy. Astron. Soc. 485, 2806–2824 (2019), arXiv:1806.09497 [astro-ph.CO] .
  * Colavincenzo  _et al._ (2019) Manuel Colavincenzo _et al._ , “Comparing approximate methods for mock catalogues and covariance matrices – III: bispectrum,” Mon. Not. Roy. Astron. Soc. 482, 4883–4905 (2019), arXiv:1806.09499 [astro-ph.CO] .
  * Lippich  _et al._ (2019) Martha Lippich _et al._ , “Comparing approximate methods for mock catalogues and covariance matrices – I. Correlation function,” Mon. Not. Roy. Astron. Soc. 482, 1786–1806 (2019), arXiv:1806.09477 [astro-ph.CO] .
  * Sehgal  _et al._ (2010) Neelima Sehgal, Paul Bode, Sudeep Das, Carlos Hernandez-Monteagudo, Kevin Huffenberger, Yen-Ting Lin, Jeremiah P. Ostriker, and Hy Trac, “Simulations of the Microwave Sky,” Astrophys. J. 709, 920–936 (2010), arXiv:0908.0540 [astro-ph.CO] .
  * Li  _et al._ (2020) Zack Li, Giuseppe Puglisi, Mathew Madavacheril, and Marcelo Alvarez, “Websky: simulated catalogs and maps of radio galaxies,” in prep. (2020).
  * Böhm  _et al._ (2020) Vanessa Böhm, Chirag Modi, and Emanuele Castorina, “Lensing corrections on galaxy-lensing cross correlations and galaxy-galaxy auto correlations,” JCAP 03, 045 (2020), arXiv:1910.06722 [astro-ph.CO] .
  * Fabbian  _et al._ (2018) Giulio Fabbian, Matteo Calabrese, and Carmelita Carbone, “CMB weak-lensing beyond the Born approximation: a numerical approach,” JCAP 2018, 050 (2018), arXiv:1702.03317 [astro-ph.CO] .
  * Hivon  _et al._ (2002) E. Hivon _et al._ , “Master of the cmb anisotropy power spectrum: A fast method for statistical analysis of large and complex cmb data sets,” Astrophys. J. 567, 2–17 (2002), astro-ph/0105302 .
  * Alonso  _et al._ (2019) David Alonso, Javier Sanchez, Anže Slosar, and LSST Dark Energy Science Collaboration, “A unified pseudo-Cℓ framework,” MNRAS 484, 4127–4151 (2019), arXiv:1809.09603 [astro-ph.CO] .
  * Makiya  _et al._ (2018) Ryu Makiya, Shin’ichiro Ando, and Eiichiro Komatsu, “Joint analysis of the thermal Sunyaev–Zeldovich effect and 2MASS galaxies: probing gas physics in the local Universe and beyond,” Mon. Not. Roy. Astron. Soc. 480, 3928–3941 (2018), arXiv:1804.05008 [astro-ph.CO] .
  * Gupta  _et al._ (2019) N. Gupta _et al._ (SPT), “Fractional polarization of extragalactic sources in the 500 deg2 SPTpol survey,” Mon. Not. Roy. Astron. Soc. 490, 5712–5721 (2019), arXiv:1907.02156 [astro-ph.CO] .
  * Aghanim  _et al._ (2020b) N. Aghanim _et al._ (Planck), “Planck 2018 results. VIII. Gravitational lensing,” Astron. Astrophys. 641, A8 (2020b), arXiv:1807.06210 [astro-ph.CO] .
  * Okamoto and Hu (2003) Takemi Okamoto and Wayne Hu, “CMB lensing reconstruction on the full sky,” Phys. Rev. D 67, 083002 (2003), arXiv:astro-ph/0301031 [astro-ph] .
  * Raghunathan  _et al._ (2021) Srinivasan Raghunathan, Nathan Whitehorn, Marcelo A. Alvarez, Han Aung, Nicholas Battaglia, Gilbert P. Holder, Daisuke Nagai, Elena Pierpaoli, Christian L. Reichardt, and Joaquin D. Vieira, “Constraining Cluster Virialization Mechanism and Cosmology using Thermal-SZ-selected clusters from Future CMB Surveys,” (2021), arXiv:2107.10250 [astro-ph.CO] .
  * Smith  _et al._ (2007) Kendrick M. Smith, Oliver Zahn, and Olivier Dore, “Detection of Gravitational Lensing in the Cosmic Microwave Background,” Phys. Rev. D 76, 043510 (2007), arXiv:0705.3980 [astro-ph] .
  * Mirmelstein  _et al._ (2019) Mark Mirmelstein, Julien Carron, and Antony Lewis, “Optimal filtering for CMB lensing reconstruction,” Phys. Rev. D100, 123509 (2019), arXiv:1909.02653 [astro-ph.CO] .
  * Hanson  _et al._ (2011) Duncan Hanson, Anthony Challinor, George Efstathiou, and Pawel Bielewicz, “CMB temperature lensing power reconstruction,” Phys. Rev. D 83, 043005 (2011), arXiv:1008.4403 [astro-ph.CO] .
  * Lewis  _et al._ (2011) Antony Lewis, Anthony Challinor, and Duncan Hanson, “The shape of the CMB lensing bispectrum,” JCAP 1103, 018 (2011), arXiv:1101.2234 [astro-ph.CO] .
  * Namikawa  _et al._ (2013) Toshiya Namikawa, Duncan Hanson, and Ryuichi Takahashi, “Bias-Hardened CMB Lensing,” Mon. Not. Roy. Astron. Soc. 431, 609–620 (2013), arXiv:1209.0091 [astro-ph.CO] .
  * Adachi  _et al._ (2020) S. Adachi _et al._ (POLARBEAR), “Internal delensing of Cosmic Microwave Background polarization $B$-modes with the POLARBEAR experiment,” Phys. Rev. Lett. 124, 131301 (2020), arXiv:1909.13832 [astro-ph.CO] .
  * Han  _et al._ (2021) Dongwon Han _et al._ (ACT), “The Atacama Cosmology Telescope: delensed power spectra and parameters,” JCAP 01, 031 (2021), arXiv:2007.14405 [astro-ph.CO] .
  * Fabbian and Stompor (2013) Giulio Fabbian and Radek Stompor, “High precision simulations of weak lensing effect on Cosmic Microwave Background polarization,” Astron. Astrophys. 556, A109 (2013), arXiv:1303.6550 [astro-ph.CO] .
  * Manzotti  _et al._ (2017) A. Manzotti _et al._ (Herschel, SPT), “CMB Polarization B-mode Delensing with SPTpol and Herschel,” Astrophys. J. 846, 45 (2017), arXiv:1701.04396 [astro-ph.CO] .
  * Larsen  _et al._ (2016) Patricia Larsen, Anthony Challinor, Blake D. Sherwin, and Daisy Mak, “Demonstration of cosmic microwave background delensing using the cosmic infrared background,” Phys. Rev. Lett. 117, 151102 (2016), arXiv:1607.05733 [astro-ph.CO] .
  * Ade  _et al._ (2016c) P. A. R. Ade _et al._ (Planck), “Planck 2015 results. XXIII. The thermal Sunyaev-Zeldovich effect–cosmic infrared background correlation,” Astron. Astrophys. 594, A23 (2016c), arXiv:1509.06555 [astro-ph.CO] .
  * Madhavacheril  _et al._ (2017) Mathew S. Madhavacheril, Nicholas Battaglia, and Hironao Miyatake, “Fundamental physics from future weak-lensing calibrated Sunyaev-Zeldovich galaxy cluster counts,” Phys. Rev. D 96, 103525 (2017), arXiv:1708.07502 [astro-ph.CO] .
  * Raghunathan  _et al._ (2017) Srinivasan Raghunathan, Sanjaykumar Patil, Eric J. Baxter, Federico Bianchini, Lindsey E. Bleem, Thomas M. Crawford, Gilbert P. Holder, Alessandro Manzotti, and Christian L. Reichardt, “Measuring galaxy cluster masses with CMB lensing using a Maximum Likelihood estimator: Statistical and systematic error budgets for future experiments,” JCAP 08, 030 (2017), arXiv:1705.00411 [astro-ph.CO] .
  * Miyatake  _et al._ (2019) Hironao Miyatake _et al._ , “Weak-lensing Mass Calibration of ACTPol Sunyaev–Zel’dovich Clusters with the Hyper Suprime-Cam Survey,” Astrophys. J. 875, 63 (2019), arXiv:1804.05873 [astro-ph.CO] .
  * Raghunathan  _et al._ (2018) Srinivasan Raghunathan, Federico Bianchini, and Christian L. Reichardt, “Imprints of gravitational lensing in the Planck cosmic microwave background data at the location of WISExSCOS galaxies,” Phys. Rev. D 98, 043506 (2018), arXiv:1710.09770 [astro-ph.CO] .
  * Raghunathan  _et al._ (2019a) S. Raghunathan _et al._ (DES), “Mass Calibration of Optically Selected DES clusters using a Measurement of CMB-Cluster Lensing with SPTpol Data,” Astrophys. J. 872, 170 (2019a), arXiv:1810.10998 [astro-ph.CO] .
  * Seljak and Zaldarriaga (2000) Uros Seljak and Matias Zaldarriaga, “Lensing induced cluster signatures in cosmic microwave background,” Astrophys. J. 538, 57–64 (2000), astro-ph/9907254 .
  * Dodelson (2004) Scott Dodelson, “CMB-cluster lensing,” Phys. Rev. D 70, 023009 (2004), astro-ph/0402314 .
  * Holder and Kosowsky (2004) Gilbert P. Holder and Arthur Kosowsky, “Gravitational lensing of the microwave background by galaxy clusters,” Astrophys. J. 616, 8–15 (2004), astro-ph/0401519 .
  * Lewis and King (2006) Antony Lewis and Lindsay King, “Cluster masses from cmb and galaxy weak lensing,” Phys. Rev. D 73, 063006 (2006), astro-ph/0512104 .
  * Hu  _et al._ (2007) Wayne Hu, Simon DeDeo, and Chris Vale, “Cluster mass estimators from cmb temperature and polarization lensing,” New J. Phys. 9, 441 (2007), astro-ph/0701276 .
  * Raghunathan  _et al._ (2019b) Srinivasan Raghunathan, Gilbert P. Holder, James G. Bartlett, SanjayKumar Patil, Christian L. Reichardt, and Nathan Whitehorn, “An Inpainting Approach to Tackle the Kinematic and Thermal SZ Induced Biases in CMB-Cluster Lensing Estimators,” JCAP 11, 037 (2019b), arXiv:1904.13392 [astro-ph.CO] .
  * Maturi  _et al._ (2005) Matteo Maturi, Matthias Bartelmann, Massimo Meneghetti, and Lauro Moscardini, “Gravitational lensing of the cmb by galaxy clusters,” A&A 436, 37–46 (2005), astro-ph/0408064 .
  * Yoo and Zaldarriaga (2008) Jaiyul Yoo and Matias Zaldarriaga, “Improved estimation of cluster mass profiles from the cosmic microwave background,” Phys. Rev. D 78, 083002 (2008), arXiv:0805.2155 [astro-ph] .
  * Ade  _et al._ (2016d) P. A. R. Ade _et al._ (Planck), “Planck 2015 results. XXVII. The Second Planck Catalogue of Sunyaev-Zeldovich Sources,” Astron. Astrophys. 594, A27 (2016d), arXiv:1502.01598 [astro-ph.CO] .
  * Patil  _et al._ (2020) Sanjaykumar Patil, Srinivasan Raghunathan, and Christian L. Reichardt, “Suppressing the Thermal SZ-induced Variance in CMB-cluster Lensing Estimators,” Astrophys. J. 888, 9 (2020), arXiv:1905.07943 [astro-ph.CO] .
  * Reuter  _et al._ (2020) C. Reuter, J. D. Vieira, J. S. Spilker, A. Weiss, M. Aravena, M. Archipley, M. Béthermin, S. C. Chapman, C. De Breuck, C. Dong, W. B. Everett, J. Fu, T. R. Greve, C. C. Hayward, R. Hill, Y. Hezaveh, S. Jarugula, K. Litke, M. Malkan, D. P. Marrone, D. Narayanan, K. A. Phadke, A. A. Stark, and M. L. Strandet, “The Complete Redshift Distribution of Dusty Star-forming Galaxies from the SPT-SZ Survey,” Astrophys. J. 902, 78 (2020), arXiv:2006.14060 [astro-ph.GA] .
  * Cai  _et al._ (2013) Zhen-Yi Cai, Andrea Lapi, Jun-Qing Xia, Gianfranco De Zotti, Mattia Negrello, Carlotta Gruppioni, Emma Rigby, Guillaume Castex, Jacques Delabrouille, and Luigi Danese, “A hybrid model for the evolution of galaxies and Active Galactic Nuclei in the infrared,” Astrophys. J. 768, 21 (2013), arXiv:1303.2335 [astro-ph.CO] .
  * Lapi  _et al._ (2012) A. Lapi, M. Negrello, J. Gonzalez-Nuevo, Z. Y. Cai, G. De Zotti, and L. Danese, “Effective Models for Statistical Studies of Galaxy-Scale Gravitational Lensing,” Astrophys. J. 755, 46 (2012), arXiv:1206.1142 [astro-ph.CO] .
  * Alvarez  _et al._ (2021) Marcelo Alvarez, Federico Bianchini, Giulio Fabbian, Mattia Negrello, and Giuseppe Puglisi, “Private communication,” (2021).
  * Choi  _et al._ (2020) Steve K. Choi _et al._ , “Sensitivity of the Prime-Cam Instrument on the CCAT-prime Telescope,” J. Low Temp. Phys. 199, 1089–1097 (2020), arXiv:1908.10451 [astro-ph.IM] .
  * Hirata and Seljak (2003a) Christopher M. Hirata and Uros Seljak, “Analyzing weak lensing of the cosmic microwave background using the likelihood function,” Phys. Rev. D 67, 043001 (2003a), arXiv:astro-ph/0209489 [astro-ph] .
  * Hirata and Seljak (2003b) Christopher M. Hirata and Uros Seljak, “Reconstruction of lensing from the cosmic microwave background polarization,” Phys. Rev. D 68, 083002 (2003b), arXiv:astro-ph/0306354 [astro-ph] .
  * Carron and Lewis (2017) Julien Carron and Antony Lewis, “Maximum a posteriori CMB lensing reconstruction,” Phys. Rev. D 96, 063510 (2017), arXiv:1704.08230 [astro-ph.CO] .
  * Ponce  _et al._ (2019) Marcelo Ponce, Ramses van Zon, Scott Northrup, Daniel Gruner, Joseph Chen, Fatih Ertinaz, Alexey Fedoseev, Leslie Groer, Fei Mao, Bruno C. Mundim, Mike Nolta, Jaime Pinto, Marco Saldarriaga, Vladimir Slavnic, Erik Spence, Ching-Hsing Yu, and W. Richard Peltier, “Deploying a Top-100 Supercomputer for Large Parallel Workloads: the Niagara Supercomputer,” PEARC ’19: Proceedings of the Practice and Experience in Advanced Research Computing on Rise of the Machines (learning) , 24 (2019), arXiv:1907.13600 [cs.DC] .
  * Wu  _et al._ (2019) W. L. K. Wu _et al._ , “A Measurement of the Cosmic Microwave Background Lensing Potential and Power Spectrum from 500 deg2 of SPTpol Temperature and Polarization Data,” Astrophys. J. 884, 70 (2019), arXiv:1905.05777 [astro-ph.CO] .
  * Carron (2019) Julien Carron, “Position-space curved-sky anisotropy quadratic estimation,” (2019), arXiv:1908.02016 [astro-ph.CO] .
  * Wandelt  _et al._ (2001) Benjamin D. Wandelt, Eric Hivon, and Krzysztof M. Gorski, “The Pseudo-$C_{l}$ method: Cosmic microwave background anisotropy power spectrum statistics for high precision cosmology,” Phys. Rev. D 64, 083003 (2001), astro-ph/0008111 .
  * Challinor and Chon (2002) Anthony Challinor and Gayoung Chon, “Geometry of weak lensing of CMB polarization,” Phys. Rev. D 66, 127301 (2002), arXiv:astro-ph/0301064 [astro-ph] .


capbtabboxtable[][] Send correspondence to Ipek<EMAIL_ADDRESS>
# Unsupervised Denoising of Retinal OCT with Diffusion Probabilistic Model

Dewei Hu  Vanderbilt University, Dept. of Electrical and Computer Engineering
Yuankai K. Tao  Vanderbilt University, Dept. of Biomedical Engineering  Ipek
Oguz  Vanderbilt University, Dept. of Computer Science, Nashville, TN, USA

###### Abstract

Optical coherence tomography (OCT) is a prevalent non-invasive imaging method
which provides high resolution volumetric visualization of retina. However,
its inherent defect, the speckle noise, can seriously deteriorate the tissue
visibility in OCT. Deep learning based approaches have been widely used for
image restoration, but most of these require a noise-free reference image for
supervision. In this study, we present a diffusion probabilistic model that is
fully unsupervised to learn from noise instead of signal. A diffusion process
is defined by adding a sequence of Gaussian noise to self-fused OCT b-scans.
Then the reverse process of diffusion, modeled by a Markov chain, provides an
adjustable level of denoising. Our experiment results demonstrate that our
method can significantly improve the image quality with a simple working
pipeline and a small amount of training data. The implementation is available
at https://github.com/DeweiHu/OCT_DDPM.

###### keywords:

OCT, denoising, unsupervised, diffusion probabilistic model

##  1 INTRODUCTION

Due to limited spatial-frequency bandwidth, optical coherence tomography (OCT)
[1]has an inherent characteristic of speckle [2]. The speckle can severely
degrade the image quality by occluding essential anatomical structures such as
retinal vessels and thin layers (e.g., the external limiting membrane (ELM)).
Hence, despeckling is an urgent pre-processing step for both clinical
diagnoses [3] and further image analysis [4]. The traditional way to reduce
speckle noise is to average multiple noisy b-scan acquisitions at the same
location [5]. This approach requires a large number of repetitions, and the
prolonged acquisition time can be problematic for patient comfort;
additionally, registration artifacts caused by eye movement can be an issue.
Therefore, a denoising algorithm that does not require repeated acquisitions
is desirable.

Many deep learning methods have been investigated for the OCT denoising
problem [6, 7, 8, 9, 10]. Most of these train the model using low-noise
reference images, which is not always available in practice. To tackle this
lack of ground truth, some unsupervised approaches including Noise2Noise [11]
(N2N) and deep image prior [12] (DIP) have been applied to OCT image
restoration. Mao et al. [9] implement N2N on OCT denoising by using two
repeated b-scans as input and target for model training. Fan et al. [10]
demonstrate that DIP can have the drawback of serious over-smoothing in
despeckling results.

In recent research, Ho et al. proposed a diffusion probabilistic model [13]
that is used for image synthesis with a performance superior to generative
adversarial networks (GANs) [14]. The general idea of the diffusion model is
straightforward: Given a natural image, a Markov chain can be formed by adding
a small amount of Gaussian noise at each step. Given a sufficient number of
steps, a complex data distribution will finally be transformed to a Gaussian.
Conversely, given a Gaussian noise image, a meaningful image can be
synthesized in a reverse process. In this work, we propose to leverage the
diffusion probabilistic model to denoise retinal OCT b-scans. Since the model
is learning the speckle pattern instead of the retina appearance, the
reference image used for training is not required to be the true noise-free
image. In our experiments we apply the self-fusion [15, 16] method to obtain
the clean reference image and we train the parameterized Markov chain by
variational inference. As the number of reverse steps is adjustable, the
algorithm is able to produce different levels of denoising results. This is an
advantage since different tasks may require different levels of fine detail
retention in the images.

##  2 METHOD

Figure 1: General workflow in this study. In (a) the target b-scan is
highlighted in red. The surrounding b-scans within a radius $r$ are used as
atlases for reconstructing the target b-scan with higher SNR. In (b) the black
arrows indicate the diffusion direction while the red represent sampling.

We leverage our self-fusion method [15] as a pre-processing step in the
training stage (Fig. 1a). Self-fusion regards b-scans in a small vicinity of a
given target b-scan as ‘atlases’ for that b-scan, because of their structural
similarity. After registering the neighbors to the target b-scan, a pixel-wise
weighted average of these ‘atlases’ will result in an image with high signal-
to-noise ratio (SNR). This approach is easy to implement and robust for
retinal layer enhancement, but finer features like vessels and texture can be
over-smoothed. Nevertheless, since the diffusion probabilistic model aims to
learn the speckle pattern instead of the signal, the self-fusion output
$x_{0}$ can still be used as the clean image for our training purposes. Fig.
1b shows a Markov chain in forward (diffuse) and reverse (denoise) directions:

$q(x_{1:T}|x_{0})=\prod_{t=1}^{T}q(x_{t}|x_{t-1}),\qquad
p_{\theta}(x_{0:T})=p(x_{T})\prod_{t=1}^{T}p_{\theta}(x_{t-1}|x_{t})$ (1)

where $q(x_{0})$ represents the data distribution while
$p(x_{T})=\mathcal{N}(x_{t};\bm{0},\bm{I})$. $\theta$ is the parameters of the
model. The diffusion and sampling describe a transition between these two
distributions with $T$ discretized steps. Our goal is to train a deep model
$p_{\theta}$ to restore a noisy image $x$ with an adjustable parameter $t$.
Intuitively, an image with stronger speckle require a larger $t$ value that
indicates more denoising steps.

The image sequence $x_{0},x_{1},\ldots,x_{T}$ is created by gradually adding
small Gaussian noise with a variance schedule
$\\{\beta_{1},\ldots,\beta_{T}\\}$, where $\beta_{t}\in(0,1)$, $\forall
t\in(1,T)$.

$q(x_{t}|x_{t-1}):=\mathcal{N}(x_{t};\sqrt{\alpha_{t}}x_{t-1},\beta_{t}\bm{I})\quad\textrm{where}\quad\alpha_{t}=1-\beta_{t}$
(2)

Eq. 2 approximates the posterior distribution in the forward process assuming
that there is a small mean shift after one step of diffusion. Denote
$\bar{\alpha}_{t}=\prod_{s=1}^{t}\alpha_{s}$, then $x_{t}$ is acquired by
adding $t$ different Gaussian random variables to $x_{0}$.

$q(x_{t}|x_{0})=\mathcal{N}(x_{t};\sqrt{\bar{\alpha}_{t}}x_{0},(1-\bar{\alpha}_{t})\bm{I})$
(3)

As the sum of $t$ Gaussians is a Gaussian with variance
$\sum_{t=1}^{T}\beta_{t}\approx 1-\bar{\alpha}_{t}$. The high order terms with
regard to $\beta_{t}$ in $1-\bar{\alpha}_{t}$ are negligible because
$\beta_{t}$ is a small value in range $(0,1)$. In practice, Eq. 3 enables
sampling of $x_{t}$ with reparameterization:

$x_{t}(x_{0},\epsilon)=\sqrt{\bar{\alpha}_{t}}x_{0}+\sqrt{1-\bar{\alpha}_{t}}\epsilon,\quad\epsilon\in\mathcal{N}(\bm{0},\bm{I})$
(4)

Similar to Eq. 2, the reverse step is also modeled as a Gaussian since the
noise added in each step is small:

$p_{\theta}(x_{t-1}|x_{t})=\mathcal{N}(x_{t-1};\bm{\mu}_{\theta}(x_{t},t),\bm{\Sigma}_{\theta}(x_{t},t))$
(5)

In this work, we set the variance to be a fixed parameter
$\bm{\Sigma}_{\theta}(x_{t},t)=\beta_{t}$ and only learn to predict the mean.

To perfectly recover the image in the reverse process, the ideal solution is
to minimize the distance between $q(x_{t-1}|x_{t})$ and $p(x_{t-1}|x_{t})$.
However, according to Ho et al. [13], the direct KL-divergence between these
two distributions, $D_{KL}(q(x_{t-1}|x_{t})||p_{\theta}(x_{t-1}|x_{t}))$, is
not tractable. Alternatively, they introduce a tractable constraint
$D_{KL}(q(x_{t-1}|x_{t},x_{0})||p_{\theta}(x_{t-1}|x_{t}))$. By leveraging the
property of Markov chain, the following equation should hold:

$q(x_{t-1}|x_{t},x_{0})=\frac{q(x_{t}|x_{t-1},x_{0})q(x_{t-1}|x_{0})}{q(x_{t}|x_{0})}=q(x_{t}|x_{t-1})q(x_{t-1}|x_{0})$
(6)

From Eq. 2 and Eq. 3, we know that both terms in the product are Gaussian;
then $q(x_{t-1}|x_{t},x_{0})$ should also have the form
$\mathcal{N}(x_{t-1};\bm{\tilde{\mu}}_{t}(x_{t},x_{0}),\tilde{\beta_{t}}\bm{I})$:

$q(x_{t-1}|x_{t},x_{0})=\mathcal{N}\left(x_{t-1};\underbrace{\frac{\sqrt{\bar{\alpha}_{t-1}}\beta_{t}}{1-\bar{\alpha}_{t}}x_{0}+\frac{\sqrt{\alpha_{t}}(1-\bar{\alpha}_{t-1})}{1-\bar{\alpha}_{t}}x_{t}(x_{0},\epsilon)}_{\bm{\tilde{\mu}}_{t}(x_{t},x_{0})},\underbrace{\frac{1-\bar{\alpha}_{t-1}}{1-\bar{\alpha}_{t}}\beta_{t}}_{\tilde{\beta_{t}}}\bm{I}\right),\quad\epsilon\in\mathcal{N}(\bm{0},\bm{I})$
(7)

We use the negative log likelihood $-\log p_{\theta}(x_{0})$ as the objective
function. This can be optimized by minimizing its variational upper bound
$\mathcal{L}$ given by the Jensen’s inequality; the detailed derivation of
this is included in Appendix A:

$\mathcal{L}=\underbrace{D_{KL}(q(x_{T}|x_{0})\ \|\
p(x_{T}))}_{\mathcal{L}_{T}}+\sum_{t=2}^{T}\underbrace{D_{KL}(q(x_{t-1}|x_{t},x_{0})\
\|\
p_{\theta}(x_{t-1}|x_{t}))}_{\mathcal{L}_{t-1}}+\underbrace{\mathcal{H}(p_{\theta}(x_{0}|x_{1}))}_{\mathcal{L}_{0}}$
(8)

Because the variance schedule is fixed in this implementation,
$\mathcal{L}_{T}$ turns out to be a constant, so we only need to consider
$\mathcal{L}_{t-1}$ and $\mathcal{L}_{0}$ as loss function. Given Eq. 7 and
Eq. 5, $\mathcal{L}_{t-1}$ is the KL divergence of two Gaussian distributions
and can be reduced to Eq. 9. The derivation details can be found in Appendix
B.

$\mathcal{L}_{t-1}=\frac{1}{2\beta_{t}}\left\|\bm{\tilde{\mu}}_{t}(x_{t},x_{0})-\bm{\mu}_{\theta}(x_{t}(x_{0},\epsilon),t)\right\|^{2}$
(9)

Obviously, to minimize $\mathcal{L}_{t-1}$ we can set the mean prediction
equal to $\bm{\tilde{\mu}}_{t}(x_{t},x_{0})$ which is derived from Eq. 7 and
Eq. 4:

$\bm{\mu}_{\theta}(x_{t},t)=\frac{1}{\sqrt{\alpha_{t}}}\left(x_{t}-\frac{\beta_{t}}{\sqrt{1-\bar{\alpha_{t}}}}\epsilon_{\theta}(x_{t},t)\right)$
(10)

Combining Eq. 9 and Eq. 10 it is easy to see that the model is learning to
predict the sample $\epsilon$ drawn from a normal distribution with mean
square error (MSE). The expression of $\bm{\mu}_{\theta}(x_{t},t)$
demonstrates that the image restoration is actually done by adding a Gaussian
noise.

##  3 EXPERIMENTS

###  3.1 Dataset

We train our model on 6 optic nerve head (ONH) volumes of the human retina,
and test it on 6 fovea volumes. Each volume contains 500 b-scans of $512\times
500$ pixels. To test the performance of the model at different speckle levels,
the data is acquired with three different levels of signal-to-noise ratio
(92dB, 96dB, 101dB). As mentioned in Sec. 1, the ground truth used for
evaluation is often obtained by averaging repeated b-scan frames at the same
spatial location. In our dataset, we have 5 repeated acquisitions for each
b-scan which are averaged together to form the ground truth.

###  3.2 Implementation details

For self-fusion, we take adjacent b-scans within radius of 3 as candidates to
compute the weighted mean. We use the
greedy111https://github.com/pyushkevich/greedy software for diffeomorphic
image registration [17]. Images are padded to $512\times 512$ pixels, and the
intensity is normalized to $[-1,1]$.

The variance schedule is set to linearly increase from $10^{-4}$ to $6\times
10^{-3}$ in $T=100$ steps. The architecture of our model is a simplified
version of that in Ho et al. [13]. It is trained on an NVIDIA RTX 2080TI 11GB
GPU for 500 epochs, with a batch size of 2 using Adam optimizer. The starting
learning rate is $10^{-4}$ and decay by half every 5 epochs.

| SNR=92dB | SNR=96dB | SNR=101dB
---|---|---|---
t=0 |  |  |
t=20 |  |  |
t=41 |  |  |
t=46 |  |  |
t=51 |  |  |
t=70 |  |  |
Figure 2: Fovea denoising results for different input SNR levels and for
different t values. The best result (determined visually) for each SNR level
is highlighted with a red box. The red arrows point to the ELM, which becomes
visible with our denoising model. (Excess background trimmed.)

##  4 RESULTS

###  4.1 Qualitative results

In Fig. 2 we show the denoising results of our model for a range of $t$
values. $t=0$ is the original input image for each SNR level. Increasing $t$
values indicate more denoising steps. The best result (determined visually)
for each noise level, as highlighted by the red box, coincide with the
intuition that the noisier images benefit from a larger $t$. For the third
column, where the input noise level is relatively low, we can see that as $t$
increases from 41 to 51, retinal layers gradually become over-smoothed and
fine texture features fade away. As discussed in Sec. 2, Eq. 10 explains that
the denoising process is done by adding Gaussian noise to compensate for the
speckle pattern. When the $t$ is too large (e.g. $t=70$ in Fig. 2), the added
noise becomes excessive and produces poor results.

We further note that our proposed method performs well for vessel and layer
preservation. For example, in the second column, our result reveals the very
thin external limiting membrane (ELM) (marked by red arrows) which is hardly
visible in the noisy input.

###  4.2 Comparison to baseline denoising model

Hu et al. [8] present a pseudo-modality fusion network (PMFN) that improves
the feature preservation compared to a former method developed by Devalla et
al. [7]. We use the PMFN as the baseline in this study. In Fig. 3, we observe
that the retinal layers are more homogeneous in our proposed approach than in
PMFN for all input SNR levels. Downstream analysis tasks such as layer
segmentation would likely benefit from this improvement. We also note that
small features like vessels are not sacrificed, even though other regions of
the layers become denoised. To quantitatively confirm these observations, we
use the average of 5 repeated frames (5-mean) as the reference ground truth
image and we report several metrics in Table 1. Comparing with PMFN, the
proposed method improve the denoising performance in terms of SNR, CNR and
ENL. The results are significantly different in a paired, two-tailed t-test
with a significance threshold of 0.05.

| SNR | PSNR | CNR | ENL
---|---|---|---|---
PMFN | $29.18\pm 2.03$ | $81.51\pm 0.69$ | $1.89\pm 0.57$ | $10.91\pm 2.80$
proposed | $\bm{40.94\pm 1.78}$ | $74.67\pm 0.58$ | $\bm{2.12\pm 0.71}$ | $\bm{54.66\pm 15.84}$

Table 1: Quantitative evaluation. SNR: signal-to-noise ratio, PSNR: peak signal-to-noise ratio, CNR: contrast-to-noise ratio, ENL: equivalent number of looks. Improvements over the baseline are highlighted in boldface. | SNR=92dB | SNR=96dB | SNR=101dB
---|---|---|---
5-mean |  |  |
PMFN |  |  |
Proposed |  |  |
Figure 3: Result comparison with baseline model. 5-mean refers to the average
image of 5 repeated b-scans.

##  5 CONCLUSION

We propose an OCT restoration method leveraging a diffusion probabilistic
model. The model is unsupervised and requires a small amount of data to train.
Moreover, the parameter $t$ provides control over the output smoothness level.
Our results show that our model can efficiently suppress the speckle;
furthermore, the features including layers and vessels are not only preserved
but present higher visibility. A limitation of this approach is its Gaussian
assumption on the speckle pattern. The generalization to other noise
distribution types will be a potential direction of future work for better
speckle modeling.

##  6 ACKNOWLEDGEMENTS

This work is supported by NIH R01EY031769, NIH R01EY030490 and the Vanderbilt
University Discovery Grant Program.

## References

  * [1] Huang, D., Swanson, E. A., Lin, C. P., Schuman, J. S., Stinson, W. G., Chang, W., Hee, M. R., Flotte, T., Gregory, K., Puliafito, C. A., et al., “Optical coherence tomography,” science 254(5035), 1178–1181 (1991).
  * [2] Schmitt, J. M., Xiang, S., and Yung, K. M., “Speckle in optical coherence tomography: an overview,” in [Saratov Fall Meeting’98: Light Scattering Technologies for Mechanics, Biomedicine, and Material Science ], 3726, 450–461, International Society for Optics and Photonics (1999).
  * [3] Chiu, S. J., Allingham, M. J., Mettu, P. S., Cousins, S. W., Izatt, J. A., and Farsiu, S., “Kernel regression based segmentation of optical coherence tomography images with diabetic macular edema,” Biomedical optics express 6(4), 1172–1194 (2015).
  * [4] Xiang, D., Tian, H., Yang, X., Shi, F., Zhu, W., Chen, H., and Chen, X., “Automatic segmentation of retinal layer in oct images with choroidal neovascularization,” IEEE Transactions on Image Processing 27(12), 5880–5891 (2018).
  * [5] Sander, B., Larsen, M., Thrane, L., Hougaard, J. L., and Jørgensen, T. M., “Enhanced optical coherence tomography imaging by multiple scan averaging,” British Journal of Ophthalmology 89(2), 207–212 (2005).
  * [6] Ma, Y., Chen, X., Zhu, W., Cheng, X., Xiang, D., and Shi, F., “Speckle noise reduction in optical coherence tomography images based on edge-sensitive cgan,” Biomedical optics express 9(11), 5129–5146 (2018).
  * [7] Devalla, S. K., Subramanian, G., Pham, T. H., Wang, X., Perera, S., Tun, T. A., Aung, T., Schmetterer, L., Thiéry, A. H., and Girard, M. J., “A deep learning approach to denoise optical coherence tomography images of the optic nerve head,” Scientific reports 9(1), 1–13 (2019).
  * [8] Hu, D., Malone, J. D., Atay, Y., Tao, Y. K., and Oguz, I., “Retinal oct denoising with pseudo-multimodal fusion network,” in [International Workshop on Ophthalmic Medical Image Analysis ], 125–135, Springer (2020).
  * [9] Mao, Z., Miki, A., Mei, S., Dong, Y., Maruyama, K., Kawasaki, R., Usui, S., Matsushita, K., Nishida, K., and Chan, K., “Deep learning based noise reduction method for automatic 3d segmentation of the anterior of lamina cribrosa in optical coherence tomography volumetric scans,” Biomedical optics express 10(11), 5832–5851 (2019).
  * [10] Fan, W., Yu, H., Chen, T., and Ji, S., “Oct image restoration using non-local deep image prior,” Electronics 9(5), 784 (2020).
  * [11] Lehtinen, J., Munkberg, J., Hasselgren, J., Laine, S., Karras, T., Aittala, M., and Aila, T., “Noise2noise: Learning image restoration without clean data,” arXiv preprint arXiv:1803.04189 (2018).
  * [12] Ulyanov, D., Vedaldi, A., and Lempitsky, V., “Deep image prior,” in [Proceedings of the IEEE conference on computer vision and pattern recognition ], 9446–9454 (2018).
  * [13] Ho, J., Jain, A., and Abbeel, P., “Denoising diffusion probabilistic models,” arXiv preprint arXiv:2006.11239 (2020).
  * [14] Dhariwal, P. and Nichol, A., “Diffusion models beat gans on image synthesis,” arXiv preprint arXiv:2105.05233 (2021).
  * [15] Oguz, I., Malone, J. D., Atay, Y., and Tao, Y. K., “Self-fusion for OCT noise reduction,” in [SPIE Medical Imaging 2020: Image Processing ], 11313, 113130C (2020).
  * [16] Hu, D., Cui, C., Li, H., Larson, K. E., Tao, Y. K., and Oguz, I., “Life: A generalizable autodidactic pipeline for 3d oct-a vessel segmentation,” arXiv preprint arXiv:2107.04282 (2021).
  * [17] Yushkevich, P. A., Pluta, J., Wang, H., Wisse, L. E., Das, S., and Wolk, D., “Fast automatic segmentation of hippocampal subfields and medial temporal lobe subregions in 3 tesla and 7 tesla t2-weighted mri,” Alzheimer’s & Dementia 7(12), P126–P127 (2016).

##  Appendix A variational upper bound

According to the Jensen’s inequality, for a concave function $f$, we should
have:

$f(\mathbb{E}[x])\geq\mathbb{E}[f(x)]$ (11)

Hence the log likelihood $\log p_{\theta}(x_{0})$ has a evidence lower bound
(ELBO), as logarithm is a concave function:

$\displaystyle\log p_{\theta}(x_{0})$
$\displaystyle=\log\int_{x_{1:T}}p_{\theta}(x_{0}|x_{1:T})p(x_{1:T})$
$\displaystyle=\log\int_{x_{1:T}}p_{\theta}(x_{0:T})\frac{q(x_{1:T}|x_{0})}{q(x_{1:T}|x_{0})}$
$\displaystyle=\log\int_{x_{1:T}}\left[\frac{p_{\theta}(x_{0:T})}{q(x_{1:T}|x_{0})}\right]q(x_{1:T}|x_{0})$
$\displaystyle=\log\mathbb{E}_{x_{1:T}\sim
q(x_{1:T}|x_{0})}\left[\frac{p_{\theta}(x_{0:T})}{q(x_{1:T}|x_{0})}\right]$
$\displaystyle\geq\mathbb{E}_{x_{1:T}\sim
q(x_{1:T}|x_{0})}\left[\log\frac{p_{\theta}(x_{0:T})}{q(x_{1:T}|x_{0})}\right]$

The negative log likelihood will then have an upper bound:

$\displaystyle-\log p_{\theta}(x_{0})$
$\displaystyle\leq\mathbb{E}_{x_{1:T}\sim
q(x_{1:T}|x_{0})}\left[-\log\frac{p_{\theta}(x_{0:T})}{q(x_{1:T}|x_{0})}\right]$
$\displaystyle=\mathbb{E}_{x_{1:T}\sim
q(x_{1:T}|x_{0})}\left[-\log\frac{p(x_{T})\prod_{t=1}^{T}p_{\theta}(x_{t-1}|x_{t})}{\prod_{t=1}^{T}q(x_{t}|x_{t-1})}\right]$
$\displaystyle=\mathbb{E}_{x_{1:T}\sim q(x_{1:T}|x_{0})}\left[-\log
p(x_{T})\prod_{t=1}^{T}\frac{p_{\theta}(x_{t-1}|x_{t})}{q(x_{t}|x_{t-1})}\right]$
$\displaystyle=\mathbb{E}_{x_{1:T}\sim q(x_{1:T}|x_{0})}\left[-\log
p(x_{T})-\sum_{t=1}^{T}\log\frac{p_{\theta}(x_{t-1}|x_{t})}{q(x_{t}|x_{t-1})}\right]$

Then minimizing the negative log likelihood is equivalent to minimize the
upper bound $\mathcal{L}$.

$\displaystyle\mathcal{L}$ $\displaystyle=\mathbb{E}_{q}\left[-\log
p(x_{T})-\sum_{t=1}^{T}\log\frac{p_{\theta}(x_{t-1}|x_{t})}{q(x_{t}|x_{t-1})}\right]$
$\displaystyle=\mathbb{E}_{q}\left[-\log
p(x_{T})-\sum_{t=2}^{T}\log\frac{p_{\theta}(x_{t-1}|x_{t})}{q(x_{t}|x_{t-1})}-\log\frac{p_{\theta}(x_{0}|x_{1})}{q(x_{1}|x_{0})}\right]$

Intuitively $q(x_{t}|x_{t-1})$ can be approximated by
$q(x_{t}|x_{t-1},x_{0})$.

$q(x_{t}|x_{t-1})\approx
q(x_{t}|x_{t-1},x_{0})=\frac{q(x_{t-1}|x_{t},x_{0})q(x_{t}|x_{0})}{q(x_{t-1}|x_{0})}$

hence the bound $\mathcal{L}$ can be expressed as:

$\displaystyle\mathcal{L}$ $\displaystyle=\mathbb{E}_{q}\left[-\log
p(x_{T})-\sum_{t=2}^{T}\log\frac{p_{\theta}(x_{t-1}|x_{t})}{q(x_{t-1}|x_{t},x_{0})}\frac{q(x_{t-1}|x_{0})}{q(x_{t}|x_{0})}-\log\frac{p_{\theta}(x_{0}|x_{1})}{q(x_{1}|x_{0})}\right]$
$\displaystyle=\mathbb{E}_{q}\left[-\log
p(x_{T})-\sum_{t=2}^{T}\log\frac{p_{\theta}(x_{t-1}|x_{t})}{q(x_{t-1}|x_{t},x_{0})}-\log\frac{q(x_{1}|x_{0})}{q(x_{T}|x_{0})}-\log\frac{p_{\theta}(x_{0}|x_{1})}{q(x_{1}|x_{0})}\right]$
$\displaystyle=\mathbb{E}_{q}\left[\log\frac{q(x_{T}|x_{0})}{p(x_{T})}+\sum_{t=2}^{T}\log\frac{q(x_{t-1}|x_{t},x_{0})}{p_{\theta}(x_{t-1}|x_{t})}-\log
p_{\theta}(x_{0}|x_{1})\right]$
$\displaystyle=\underbrace{D_{KL}(q(x_{T}|x_{0})\ \|\
p(x_{T}))}_{\mathcal{L}_{T}}+\sum_{t=2}^{T}\underbrace{D_{KL}(q(x_{t-1}|x_{t},x_{0})\
\|\
p_{\theta}(x_{t-1}|x_{t}))}_{\mathcal{L}_{t-1}}+\underbrace{\mathcal{H}(p_{\theta}(x_{0}|x_{1}))}_{\mathcal{L}_{0}}$

##  Appendix B KL divergence of Gaussian variables

Let $p(x)$ and $q(x)$ be two Gaussian distributions:

$\displaystyle p(x)\sim\mathcal{N}(\mu_{1},\sigma_{1}^{2}),\
q(x)\sim\mathcal{N}(\mu_{2},\sigma_{2}^{2})$

The KL divergence between $p(x)$ and $q(x)$ can be derived as follows:

$\displaystyle D_{KL}(p(x)\ \|\ q(x))$
$\displaystyle=\mathbb{E}_{p}\left[\log\frac{p(x)}{q(x)}\right]=\int p(x)\log
p(x)dx-\int p(x)\log q(x)dx$ $\displaystyle=-\frac{1}{2}(1+\log
2\pi\sigma_{1}^{2})+\frac{1}{2}\log
2\pi\sigma_{2}^{2}+\frac{1}{2\sigma_{2}^{2}}\left[\sigma_{1}^{2}+(\mu_{1}-\mu_{2})^{2}\right]$
$\displaystyle=\frac{1}{2\sigma_{2}^{2}}\left[\sigma_{1}^{2}+(\mu_{1}-\mu_{2})^{2}\right]+\log\frac{\sigma_{2}}{\sigma_{1}}-\frac{1}{2}$

Apply this result on $\mathcal{L}_{t-1}$, if the ratio of variances is either
cancelled or regarded as constant

$\mathcal{L}_{t-1}=\frac{1}{2\sigma_{t}^{2}}\|\bm{\tilde{\mu_{t}}}(x_{t},x_{0})-\bm{\mu}_{\theta}(x_{t},t)\|^{2}+C$
(12)

From the reparameterization (Eq. 4) it is easy to get

$x_{0}=\frac{1}{\sqrt{\bar{\alpha_{t}}}}\left(x_{t}(x_{0},\epsilon)-\sqrt{1-\bar{\alpha_{t}}}\epsilon\right)$

Then plug this in $\bm{\tilde{\mu_{t}}}$ to get:

$\displaystyle\bm{\tilde{\mu_{t}}}(x_{t},x_{0})$
$\displaystyle=\bm{\tilde{\mu_{t}}}\left(x_{t}(x_{0},\epsilon),\frac{1}{\sqrt{\bar{\alpha_{t}}}}(x_{t}(x_{0},\epsilon)-\sqrt{1-\bar{\alpha_{t}}}\epsilon)\right)$
$\displaystyle=\frac{\sqrt{\bar{\alpha}_{t-1}}\beta_{t}}{1-\bar{\alpha}_{t}}\cdot\frac{1}{\sqrt{\bar{\alpha}}_{t}}(x_{t}(x_{0},\epsilon)-\sqrt{1-\bar{\alpha_{t}}}\epsilon)+\frac{\sqrt{\alpha_{t}}(1-\bar{\alpha}_{t-1})}{1-\bar{\alpha}_{t}}x_{t}(x_{0},\epsilon)$
$\displaystyle=\left[\frac{\beta_{t}}{(1-\bar{\alpha}_{t})\sqrt{\alpha_{t}}}+\frac{\sqrt{\alpha_{t}}(1-\bar{\alpha}_{t-1})}{1-\bar{\alpha}_{t}}\right]x_{t}(x_{0},\epsilon)-\frac{\beta_{t}}{(1-\bar{\alpha}_{t})\sqrt{\alpha_{t}}}\sqrt{1-\bar{\alpha_{t}}}\epsilon$
$\displaystyle=\frac{1}{\sqrt{\alpha_{t}}}\left(\frac{\beta_{t}+\alpha_{t}(1-\bar{\alpha}_{t-1})}{1-\bar{\alpha}_{t}}x_{t}(x_{0},\epsilon)-\frac{\beta_{t}}{\sqrt{1-\bar{\alpha_{t}}}}\epsilon\right)$
$\displaystyle=\frac{1}{\sqrt{\alpha_{t}}}\left(x_{t}(x_{0},\epsilon)-\frac{\beta_{t}}{\sqrt{1-\bar{\alpha_{t}}}}\epsilon\right)$

Then apply this result in the expression of $\mathcal{L}_{t-1}$, we get:

$\mathcal{L}_{t-1}-C=\frac{1}{2\sigma_{t}^{2}}\left\|\frac{1}{\sqrt{\alpha_{t}}}\left(x_{t}(x_{0},\epsilon)-\frac{\beta_{t}}{\sqrt{1-\bar{\alpha_{t}}}}\epsilon\right)-\bm{\mu}_{\theta}(x_{t}(x_{0},\epsilon),t)\right\|^{2}$
(13)

To minimize $\mathcal{L}_{t-1}$, given fixed $x_{t}$, the following should
hold:

$\displaystyle\bm{\mu}_{\theta}(x_{t},t)$
$\displaystyle=\frac{1}{\sqrt{\alpha_{t}}}\left(x_{t}-\frac{\beta_{t}}{\sqrt{1-\bar{\alpha_{t}}}}\epsilon_{\theta}(x_{t},t)\right)$
(14)