# Efficient detection of multidimensional single-photon time-bin
superpositions

Adam Widomski  Faculty of Physics, University of Warsaw, Pasteura 5, 02-093
Warszawa, Poland  The authors contributed equally to this work.
<EMAIL_ADDRESS>Maciej Ogrodnik  Faculty of Physics, University of
Warsaw, Pasteura 5, 02-093 Warszawa, Poland  The authors contributed equally
to this work.  Michał Karpiński  Faculty of Physics, University of Warsaw,
Pasteura 5, 02-093 Warszawa, Poland

###### Abstract

The ability to detect quantum superpositions lies at the heart of fundamental
and applied aspects of quantum mechanics. The time-frequency degree of freedom
of light enables encoding and transmitting quantum information in a multi-
dimensional fashion compatible with fiber and integrated platforms. However,
the ability to efficiently detect time-frequency superpositions is not yet
available. Here we show, that multidimensional time-bin superpositions can be
detected using a single time-resolved photon detector. Our approach uses off-
the shelf components and is based on the temporal Talbot effect – a time-
frequency counterpart of the well-know near field diffraction effect. We
provide experimental results and discuss the possible applications in quantum
communication, quantum information processing, and time-frequency quantum
state tomography.

††journal: opticajournal

##  1 Introduction

The ability to detect quantum superpositions is key for most of quantum
experiments, from quantum communication [1, 2, 3, 4], quantum state and
process tomography [5, 6, 7, 8, 9, 10, 11, 12], entanglement verification [13,
14, 15, 16, 17] to quantum computation [18]. Given the growing use of the
time-frequency degree of freedom of quantum light, which enables
multidimensional encoding of quantum information in a single spatial mode,
compatible with fiber optic and integrated optical setups [19, 20, 21, 22,
23], a necessity arises to efficiently detect single-photon superpositions in
time and frequency [24, 25, 20, 19].

The measurement of temporal mode superpositions was performed by means of the
optically nonlinear quantum pulse gate [26, 27], enabling quantum tomography
[28], and by means of electro-optic sideband generation in frequency bins,
enabling time-frequency entanglement detection [8, 7, 21]. However, those
approaches require active spectral modification, either electro-optically or
by nonlinear optical interactions. Franson interferometers [29, 30] are an
alternative for the detection of time-bin superpositions, but a single
interferometer detects the phase between only two time bins. Building a setup
enabling multidimensional state discrimination requires nesting the
interferometers [31, 32]. This implies higher complexity, cost, and in
particular, losses. In the Franson interferometer tree method the probability
of detecting a $d$-dimensional superposition scales as $1/d$, due to the
increasing number of possible paths in the nested interferometer [31].

Another possibility to approximately detect time-bin superpositions is by
directly resolving their spectrum. However, resolving the signals in time and
frequency simultaneously at the single-photon level is challenging due to the
limits given by the time-frequency uncertainty relation [20, 33]. The highest
resolution of efficient multiplexed single-photon spectral measurements is
offered by means of the dispersive Fourier Transform [34, 35, 36, 37, 38],
however, the large values of group dispersion delay needed to reach the
required spectral resolution can result in prohibitively high losses [20].

Here we show, that the temporal Talbot effect [39], a time-frequency
equivalent of the well-known near field diffraction effect [40], enables
discrimination between $d$ orthogonal states of a $d$-dimensional time-bin
superposition. Our method only requires a passive experimental setup
comprising a dispersive medium and a single time-correlated single photon
counter. We show, that whereas the discrimination is not perfect, the
information content increases with the dimension. The setup can be constructed
with off-the-shelf components, requires significantly less dispersion than a
conventional dispersive Fourier spectrometer, and offers constant detection
loss regardless of the dimension $d$.

##  2 Methods

###  2.1 Temporal Talbot effect

The Talbot effect is a phenomenon of self-imaging occurring for periodically
placed intensity distributions [40]. Those self-images are created due to
diffraction, and repeat themselves after a certain characteristic distance.
The optical space-time duality [41] lets us translate spatial effects within
the paraxial approximation to temporal effects in short pulse propagation in a
second-order dispersive medium, in the narrowband regime [42]. This means that
a superposition of equally spaced pulses is self-imaged by propagation in a
medium with a certain group dispersion delay (GDD) (which we will further
refer to as the Talbot GDD) [39], just as equally spaced pulses would be self-
imaged after some propagation distance due to diffraction (Fig. 1a). This
temporal Talbot effect enables multiple bright light applications such as
multiplication of the repetition rate of a periodic pulse train [43],
factorization of numbers [44] or real-time calculation of microwave
spectrograms [45].

In analogy to the spatial case, the temporal location of the self-images
depends on the relative phases of the sources. This enables us to employ the
temporal Talbot effect to detect different time bin superpositions. In
particular, we use the fact that the temporal Talbot effect can be seen as a
discrete Fourier transform of input pulse trains [46, 47, 48].

Figure 1: a) A Talbot carpet displaying probability of finding a photon at a
specified instance of time as a function of group delay dispersion (GDD).
Color intensity depicts the probability density. Horizontal dashed lines
represent GDD values necessary to measure the distribution in the temporal far
field and by means of the temporal Talbot effect. The carpet was generated
considering no phase difference between the pulses. b) An orthogonal Talbot
carpet representing the probability to detect a photon in time as a function
of symbol separation, calculated for $\textrm{GDD}=12900\textrm{ ps}^{2}$ and
$\tau=284$ ps corresponding to the first Talbot separation.

###  2.2 Detecting superpositions

Let us consider optical pulses in 4 time bins:

$\ket{t_{0}},\ket{t_{1}},\ket{t_{2}},\ket{t_{3}}.$ (1)

Those time bins are approximately Gaussian-shaped, $40\textrm{ ps}$ wide and
separated by $284\textrm{ ps}$ (the experimental details will be presented in
sec. 2c.).

Those states can be directly measured in time by means of time-correlated
single photon counting (TCSPC). We also created 4 superpositions of those
states that form a discrete Fourier transform [49] of the time bin basis:

$\ket{f_{n}}=\frac{1}{\sqrt{d}}\sum_{m=0}^{d-1}e^{\frac{-2\pi
inm}{d}}\ket{t_{m}},$ (2)

with the dimension $d=4$.

Let us now consider propagation through a medium with group delay dispersion
corresponding to the Talbot GDD ($\beta_{2}$). The Talbot GDD is linked to the
separation between the pulses ($\tau$) [43] as:

$\tau=\sqrt{\frac{2\pi\beta_{2}}{s}},s=1,2....$ (3)

Finally, the signals are detected by time-resolved single photon counting at
the output of the GDD medium. When the condition in Eq. 3 is met, a train of
pulses in time will be observed (cf. the dashed line in Fig. 1a due to the
temporal Talbot effect.

In Fig. 1b we show the dependence of the output time-resolved photon detection
probabilities on the symbol separation for a given GDD value of $12900$ ps2,
corresponding to the first Talbot separation in the experiment. Fringe pattern
can be observed for symbol separations matching the conditions.

In Fig. 2a we present the optical pulses representing the 4 time bins. Those
pulses were further used to form superpositions. Adequate phases are adjusted
with the voltage signals (blue rectangles) driving an electro-optic phase
modulator. The histograms corresponding to the photon times of arrival of the
states from the basis $\\{\ket{f_{i}}\\}$ are presented in the Fig. 2b.

Figure 2: a) Numerical calculation of signal’s temporal profiles. Yellow
rectangles represent temporal symbols generated by means of fast electro-optic
amplitude modulation. Blue rectangles represent voltage signals used for phase
modulation. b) Simulated single-photon time of arrival distributions for the
temporal Talbot effect for the experimental parameters provided in the text.

In our method only a single detector is used, and all detections are taken
into account, but for a fraction of detections the result may be ambiguous.
This is in contrast to the Franson interferometer tree method [32] where for a
fraction of detections we get unambiguous measurement. For single-shot state
discrimination in Talbot effect method we need a procedure to assign symbol to
a measured time of arrival. We want to do it with maximal _correctness_ , that
is a ratio of correct outcomes to errors. We could try to do that by some
probabilistic procedure based on conditional probabilities of each symbol for
each time of arrival. Next we show that the optimal strategy is to simply take
the most probable symbol at a given time.

With a fixed experimental setup this is a task of distinguishing classical
random variables which is a simple case of Bayesian decision theory [50]. Let
us denote the random variable associated with the symbol $i$ as $X_{i}$.
Symbol number $i$ is chosen randomly, described by the random variable
$\sigma$ with the uniform distribution on the $\\{0,1,\ldots,d-1\\}$. Next, a
random time of arrival value $X_{i}=t_{0}$ is drawn. Our task is to decide
from which $X_{i}$ a value $t_{0}$ was drawn. Here we consider case with
discrete time which corresponds to the realistic setup with a finite time
resolution. We want to guess the outcome of $\sigma$ given the outcome $t_{0}$
of $X_{\sigma}$. Based on probabilities $P(X_{i}=t_{0})$ we may decide the
outcome probabilistically. Given the result $t_{0}$ we give answer $i$ with a
probability $\alpha_{t_{0}}(i)$. We want to maximize the correctness, ratio of
correct results to wrong results, given by

$C=\sum_{t}\sum_{i}\alpha_{t}(i)P(\sigma=i|X_{\sigma}=t).$ (4)

By the Bayes rule we have

$P(\sigma=i|X_{\sigma}=t)=\frac{P(X_{\sigma}=t|\sigma=i)P(\sigma=i)}{P(X_{\sigma}=t)}=\frac{P(X_{i}=t)\frac{1}{d}}{P(X_{\sigma}=t)},$
(5)

so the function to optimize takes the form

$C=\sum_{t}\frac{1}{d\cdot P(X_{\sigma}=t)}\sum_{i}\alpha_{t}(i)P(X_{i}=t).$
(6)

We have

$\sum_{i}\alpha_{t}(i)P(X_{i}=t)\leq\max_{j}P(X_{j}=t)\sum_{i}\alpha_{t}(i)=\max_{j}P(X_{j}=t),$
(7)

and this bound is saturated for the deterministic strategy where
$\alpha_{t}(i)=1$ for $i$ such that $P(X_{i}=t)$ is maximal. For this game the
optimal strategy is to always choose the outcome $i$ such that

$P(X_{i}=t_{0})=\max_{j}P(X_{j}=t_{0}).$ (8)

Depending on the application it may be advantageous to have lower error rate
at the cost of lower efficiency. We can do that in the postprocessing by
taking only measurement times at which one symbol is much more probable than
others.

###  2.3 Experimental setup

Our experimental setup is shown in Fig. 3. A laser operating at the telecom
wavelength of 1560 nm in the continuous wave (CW) mode was used as source of
light. For technical reasons, we amplified the optical signal with an erbium-
doped fiber amplifier (EDFA, Pritel, HPP-PMFA-22-10) and employed a bandpass
filter to reduce the noise. Optical pulses were carved by means of fast
amplitude electro-optic modulation with a Mach-Zehnder modulator (MZM,
Thorlabs, LNA6213). Those pulses were used to form superpositions consisting
of $d=4$ single states in a predefined phase relation.The superpositions were
directed to the electro-optic phase modulator (EOPM, EOspace), which was used
to impose adequate phases on every component of the superposition. Both
modulators support 40 GHz of analog bandwidth, and were driven with an
amplified radio-frequency (RF) signal generated with the fast arbitrary
waveform generator (AWG, Keysight, M1896A) providing 33 GHz of the analog
bandiwdth and sampling rate up to 92.16 Gs/s. The full-width at half maximum
(FWHM) duration of a one-bit signal was therefore $\sim\\!12\textrm{ ps}$. The
phase factors were adjusted with the EOPM by programming four driving voltage
signals consisting of approximately $150$-ps-wide rectangular pulses (Fig. 2a)
such, that their amplitudes corresponded to fractions or multiples of the
half-wave voltage (${V_{\pi}}$). To generate optical signals we assigned 3
bits from the AWG memory per symbol. The separation between the symbols was
equal to 282 ps insted of 284 ps due to the finite sampling frequency. The AWG
was also a source of the 10 MHz clock (CLK) signal distributed over an
electrical cable to the detection section. We biased the MZM for minimal
transmission with the direct current (DC) voltage from the power supply
(Keysight, E36313A) using the feedback loop consisting of a 90:10 fiber optic
beam splitter and a power meter (Thorlabs, PM400). Finally, the pulses were
attenuated to the single photon level with the electronic variable optical
attenuator (EVOA, Thorlabs, EVOA1550F).

On the detection side, we reconnected the optical path in order to either
measure the photon counts directly in time, or to detect the superpositions.
Time-domain measurement was achieved by directing the signals to niobium
nitride superconducting nanowire single-photon detectors (SNSPDs, Single
Quantum) and histogramming the time tags using a time-to-digital converter
(Swabian Instruments, Time Tagger Ultra). The time tagger exhibited jitter
equal to $\sim\\!10\textrm{ ps}$ root mean square (RMS), and the SNSPDs of
$\sim\\!5\textrm{ ps}$ RMS. Thus, the total jitter of the receiver system was
$\sim\\!11\textrm{ ps}$ RMS. For the superposition measurement, photons were
additionally transmitted through a chirped-fiber-Bragg-grating-based
dispersion compensating module (DCM), resulting in dispersive stretching of
the pulses. We used a DCM providing group delay dispersion (GDD) equivalent to
$12900\textrm{ ps}^{2}$ (Proximion, DCMHDC-100H-P510). This amount of
dispersion was not enough to perform the dispersive Fourier transform, but
sufficed to observe the temporal Talbot effect (see Fig.1a), which we used to
distinguish superpositions of the optical signals with different phases.

Figure 3: Experimental setup. A continuous wave (CW) telecommunications laser
is modulated with an electro-optic Mach-Zehnder modulator (MZM) and electro-
optic phase modulator (EOPM) to generate optical pulses forming
superpositions. The optical signals are attenuated to the single photon level
with an electronic variable optical attenuator (EVOA). The signals are
detected either directly in the time domain (T) or in the superposition domain
(S) by means of the temporal Talbot effect performed with the dispersion
compensating module (DCM) providing group delay dispersion (GDD). Photon time-
of-arrival histograms were acquired by means of time-correlated single photon
counting (TCSPC). Yellow lines represent optical fiber connections.

##  3 Results

We demonstrated our technique with four four-dimensional superpositions. We
analysed histogram’s widths, shapes, and locations on the time axis to verify
our method. Temporal symbols were approximately-Gaussian, 46 ps wide, and
separated by approximately 284 ps, as expected. We examined these temporal
properties by utilising the $\ket{f_{0}}$ state, since all the components of
this superposition were equally modulated in phase (Fig 4a). Fig 4b presents
histograms of all the $\\{f_{i}\\}$ states after propagation through the DCM.
The obtained fringes were distinguishable in time yielding the possibility to
detect each state from the $\\{f_{i}\\}$ basis.

To display the essence of utilising the Talbot effect for efficient detection,
we measured the probability of successful measurement for symbol separations
lower and greater than the first Talbot separation. The outcomes of this
measurement are presented in Fig. 4c. We evaluated the correctness, according
to the aforementioned selection criterion, eq. 8. The result indicates that
any deviation from the Talbot separation results in the decrease of the
correctness. This indicates that the measurements are performed outside the
temporal far field regime.

Figure 4: a) Measured histogram of the $\ket{f_{0}}$ state in the time domain.
b) Measured histograms of the $\ket{f_{i}}$ states after frequency-to-time
mapping c) Measured correctness of the $\ket{f_{i}}$ state discrimination,
compared to the theoretical dependencies taking into account no jitter
(noiseless), and 11 ps root mean square (RMS) jitter (noisy) cases. The error
bars indicate sampling error due to the finite sampling frequency of the AWG.
Imperfections in symbol preparation (not included in the simulation) lead to a
small decrease in the measured correctness.

##  4 Conclusion and discussion

We propose a new method to detect time-bin states from discrete Fourier
transform basis. Those states could be also detected by a tree of Franson
interferometers as proposed in [32, 51]. For a system of dimension $d$ this
method requires $d-1$ aligned Franson interferometers and $d$ detectors and
uses postselection of results that discards $\frac{d-1}{d}$ of measurements.
Real efficiency is further limited by the insertion loss of interferometers in
the optical path. In our method only one dispersion module and one detector
are used for the measurement. It offers constant efficiency irrespective of
the dimension, at the cost of error rate increasing with the number of
different states to be distinguished. One can recognize a trade-off between
the Franson interferometer tree method with a constant error rate, decreasing
efficiency and increasing complexity, and our method with a constant
efficiency, constant complexity and error rate increasing with dimension. The
error rate in our method can also be decreased by postselection of the results
without the change in the experimental setup.

With larger alphabet we send more bits of information, but we also witness
higher error rates. To quantify that trade-off we may use mutual information
to evaluate the scaling with dimension of the amount of information that can
be sent using the Talbot effect method for detection. In particular, we
calculate mutual information between the random variable $\sigma$ of the
prepared state and the random variable of the time of arrival $X_{\sigma}$. In
Fig. 5 we can see how it increases with the increase of the dimension of the
prepared superposition.

Figure 5: Numerical simulation of the mutual information of input state and
(noisy) time of arrival random variables.

In summary, we show that the temporal Talbot effect in the single-photon-
counting regime enables efficient detection of time-bin superpositions in an
all-fiber setup. The proposed selection algorithm allowed us to distinguish
the superpositions for all detected photons with with error rate of $36.5\%$
for the Talbot separation between pulses. Compared to the Franson
interferometer tree method we have higher efficiency at the cost of higher
error rate. Other post selection strategies, including rejection of ambiguous
measurements, can alter that trade-off. Our technique can be easily expanded
to higher dimensions with the increase in bits of information per pulse. It
offers low complexity and high flexibility, which is key for quantum-photonic
applications. We expect it will enable development of robust and efficient
quantum communication methods, as well as techniques for time-frequency
entanglement detection and characterization and quantum state and process
tomography.

Funding A part of this work was carried out within the Project QuICHE,
supported by the National Science Centre of Poland project no.
2019/32/Z/ST2/00018, under QuantERA, which has received funding from the
European Union’s Horizon 2020 research and innovation programme under grant
agreement no 731473. The authors appreciate funding from the University of
Warsaw within the of the "Excellence Initiative – Research University"
framework. Acknowledgments We thank M. Mikołajczyk, R. Demkowicz-Dobrzański
and F. Sośnicki for insightful discussions.

Disclosures The authors declare no conflicts of interest.

## References

  * [1] C. H. Bennett and G. Brassard, “Quantum cryptography: Public key distribution and coin tossing,” in _Proceedings of IEEE International Conference on Computers, Systems and Signal Processing,_ (IEEE, 1984), pp. 175–179.
  * [2] A. Ekert, “Quantum cryptography based on Bell’s theorem,” Physical Review Letters 67, 661––663 (1991).
  * [3] H. Bechmann-Pasquinucci and W. Tittel, “Quantum cryptography using larger alphabets,” Physical Review A 61, 062308 (2000).
  * [4] N. J. Cerf, M. Bourennane, A. Karlsson, and N. Gisin, “Security of quantum key distribution using d-level systems,” Physical Review Letters 88, 127902 (2002).
  * [5] M. Paris and J. Rehacek, _Quantum state estimation_ , vol. 649 (Springer Science & Business Media, 2004).
  * [6] V. Ansari, J. M. Donohue, M. Allgaier, L. Sansoni, B. Brecht, J. Roslund, N. Treps, G. Harder, and C. Silberhorn, “Tomography and purification of the temporal-mode structure of quantum light,” Phys. Rev. Lett. 120, 213601 (2018).
  * [7] H.-H. Lu, K. V. Myilswamy, R. S. Bennink, S. Seshadri, M. S. Alshaykh, J. Liu, T. J. Kippenberg, D. E. Leaird, A. M. Weiner, and J. M. Lukens, “Bayesian tomography of high-dimensional on-chip biphoton frequency combs with randomized measurements,” Nature Communications 13, 4338 (2022).
  * [8] M. Kues, C. Reimer, P. Roztocki, L. R. Cortés, S. Sciara, B. Wetzel, Y. Zhang, A. Cino, S. T. Chu, B. E. Little, D. J. Moss, L. Caspani, J. Azaña, and R. Morandotti, “On-chip generation of high-dimensional entangled quantum states and their coherent control,” Nature 546, 622–626 (2017).
  * [9] L. Serino, J. Gil-Lopez, M. Stefszky, R. Ricken, C. Eigner, B. Brecht, and C. Silberhorn, “Realization of a multi-output quantum pulse gate for decoding high-dimensional temporal modes of single-photon states,” PRX Quantum 4, 020306 (2023).
  * [10] B. Jack, J. Leach, H. Ritsch, S. M. Barnett, M. J. Padgett, and S. Franke-Arnold, “Precise quantum tomography of photon pairs with entangled orbital angular momentum,” New Journal of Physics 11, 103024 (2009).
  * [11] M. Riebe, K. Kim, P. Schindler, T. Monz, P. O. Schmidt, T. K. Körber, W. Hänsel, H. Häffner, C. F. Roos, and R. Blatt, “Process tomography of ion trap quantum gates,” Phys. Rev. Lett. 97, 220407 (2006).
  * [12] M. Steffen, M. Ansmann, R. C. Bialczak, N. Katz, E. Lucero, R. McDermott, M. Neeley, E. M. Weig, A. N. Cleland, and J. M. Martinis, “Measurement of the entanglement of two superconducting qubits via state tomography,” Science 313, 1423–1425 (2006).
  * [13] J. Bavaresco, N. Herrera Valencia, C. Klöckl, M. Pivoluska, P. Erker, N. Friis, M. Malik, and M. Huber, “Measurements in two bases are sufficient for certifying high-dimensional entanglement,” Nature Phys. 14, 1032–1037 (2018).
  * [14] Y. Guo, X.-M. Hu, B.-H. Liu, Y.-F. Huang, C.-F. Li, and G.-C. Guo, “Experimental witness of genuine high-dimensional entanglement,” Phys. Rev. A 97, 062309 (2018).
  * [15] J. Schneeloch and G. A. Howland, “Quantifying high-dimensional entanglement with Einstein-Podolsky-Rosen correlations,” Phys. Rev. A 97, 042338 (2018).
  * [16] O. Gühne and G. Tóth, “Entanglement detection,” Physics Reports 474, 1–75 (2009).
  * [17] P. Hyllus, O. Gühne, D. Bruß, and M. Lewenstein, “Relations between entanglement witnesses and Bell inequalities,” Phys. Rev. A 72, 012321 (2005).
  * [18] M. A. Nielsen and I. L. Chuang, _Quantum computation and quantum information_ (Cambridge university press, 2010).
  * [19] B. Brecht, D. V. Reddy, C. Silberhorn, and M. G. Raymer, “Photon temporal modes: A complete framework for quantum information science,” Phys. Rev. X 5, 041017 (2015).
  * [20] M. Karpiński, A. O. Davis, F. Sośnicki, V. Thiel, and B. J. Smith, “Control and measurement of quantum light pulses for quantum information science and technology,” Advanced Quantum Technologies 4, 2000150 (2021).
  * [21] J. M. Lukens and P. Lougovski, “Frequency-encoded photonic qubits for scalable quantum information processing,” Optica 4, 8–16 (2017).
  * [22] P. Roztocki, S. Sciara, C. Reimer, L. R. Cortés, Y. Zhang, B. Wetzel, M. Islam, B. Fischer, A. Cino, S. T. Chu, B. E. Little, D. J. Moss, L. Caspani, J. A. na, M. Kues, and R. Morandotti, “Complex quantum state generation and coherent control based on integrated frequency combs,” J. Lightwave Technol. 37, 338–344 (2019).
  * [23] J. Roslund, R. M. De Araujo, S. Jiang, C. Fabre, and N. Treps, “Wavelength-multiplexed quantum networks with ultrafast frequency combs,” Nature Photonics 8, 109–112 (2014).
  * [24] D. Bruß and C. Macchiavello, “Optimal eavesdropping in cryptography with three-dimensional quantum states,” Phys. Rev. Lett. 88, 127901 (2002).
  * [25] L. Sheridan and V. Scarani, “Security proof for quantum key distribution using qudit systems,” Phys. Rev. A 82, 030301 (2010).
  * [26] A. Eckstein, B. Brecht, and C. Silberhorn, “A quantum pulse gate based on spectrally engineered sum frequency generation,” Opt. Express 19, 13770–13778 (2011).
  * [27] D. V. Reddy and M. G. Raymer, “High-selectivity quantum pulse gating of photonic temporal modes using all-optical Ramsey interferometry,” Optica 5, 423–428 (2018).
  * [28] V. Ansari, G. Harder, M. Allgaier, B. Brecht, and C. Silberhorn, “Temporal-mode measurement tomography of a quantum pulse gate,” Physical Review A 96, 063817 (2017).
  * [29] J. D. Franson, “Bell inequality for position and time,” Phys. Rev. Lett. 62, 2205 (1989).
  * [30] T. Zhong, H. Zhou, R. D. Horansky, C. Lee, V. B. Verma, A. E. Lita, A. Restelli, J. C. Bienfang, R. P. Mirin, T. Gerrits, S. W. Nam, F. Marsili, M. D. Shaw, Z. Zhang, L. Wang, D. Englund, G. W. Wornell, J. H. Shapiro, and F. N. C. Wong, “Photon-efficient quantum key distribution using time–energy entanglement with high-dimensional encoding,” New J. Phys. 17, 022002 (2015).
  * [31] N. T. Islam, C. Cahall, A. Aragoneses, A. Lezama, J. Kim, and D. J. Gauthier, “Robust and stable delay interferometers with application to $d$-dimensional time-frequency quantum key distribution,” Phys. Rev. Appl. 7, 044010 (2017).
  * [32] T. Brougham, S. M. Barnett, K. T. McCusker, P. G. Kwiat, and D. J. Gauthier, “Security of high-dimensional quantum key distribution protocols using Franson interferometers,” J. Phys. B 46, 104010 (2013).
  * [33] A. Widomski, S. Stopiński, K. Anders, R. Piramidowicz, and M. Karpiński, “Precise on-chip spectral and temporal control of single-photon-level optical pulses,” Journal of Lightwave Technology pp. 1–8 (2023).
  * [34] K. Goda and B. Jalali, “Dispersive Fourier transformation for fast continuous single-shot measurements,” Nature Photonics 7, 102–112 (2013).
  * [35] A. O. C. Davis, P. M. Saulnier, M. Karpiński, and B. J. Smith, “Pulsed single-photon spectrometer by frequency-to-time mapping using chirped fiber Bragg gratings,” Opt. Express 25, 12804–12811 (2017).
  * [36] M. Avenhaus, A. Eckstein, P. J. Mosley, and C. Silberhorn, “Fiber-assisted single-photon spectrograph,” Opt. Lett. 34, 2873–2875 (2009).
  * [37] M. Närhi, L. Salmela, J. Toivonen, C. Billet, J. M. Dudley, and G. Genty, “Machine learning analysis of extreme events in optical fibre modulation instability,” Nature Commun. 9, 4923 (2018).
  * [38] K. Sedziak-Kacprowicz, A. Czerwinski, and P. Kolenderski, “Tomography of time-bin quantum states using time-resolved detection,” Phys. Rev. A 102, 052420 (2020).
  * [39] T. Jannson and J. Jannson, “Temporal self-imaging effect in single-mode fibers,” J. Opt. Soc. Am. 71, 1373–1376 (1981).
  * [40] H. F. Talbot, “Facts relating to optical science. No. IV,” Philosophical Magazine 9, 401–407 (1836).
  * [41] B. H. Kolner, “Space-time duality and the theory of temporal imaging,” IEEE Journal of Quantum Electronics 30, 1951–1963 (1994).
  * [42] V. Torres-Company, J. Lancis, and P. Andres, “Space-time analogies in optics,” Prog. Opt. 56, 1–80 (2011).
  * [43] J. Azana and M. Muriel, “Temporal self-imaging effects: theory and application for multiplying pulse repetition rates,” IEEE Journal of Selected Topics in Quantum Electronics 7, 728–744 (2001).
  * [44] D. Bigourd, B. Chatel, W. P. Schleich, and B. Girard, “Factorization of numbers with the temporal Talbot effect: Optical implementation by a sequence of shaped ultrashort pulses,” Phys. Rev. Lett. 100, 030202 (2008).
  * [45] J. Azaña and X. Zhu, “Optical time-mapped spectrograms (I): from the time-lens Fourier transformer to the Talbot-based design,” Journal of Lightwave Technology (2023).
  * [46] S. Tainta, M. Erro, M. Garde, and M. Muriel, “Temporal self-imaging effect for periodically modulated trains of pulses,” Optics Express 22, 15251–15266 (2014).
  * [47] C. R. Fernández-Pousa, “A dispersion-balanced Discrete Fourier Transform of repetitive pulse sequences using temporal Talbot effect,” Optics Communications 402, 97–103 (2017).
  * [48] J. Zhou, “All-optical discrete fourier transform based on multimode interference couplers,” IEEE Photonics Technology Letters 22, 1093–1095 (2010).
  * [49] S. Barnett, _Quantum information_ , vol. 16 (Oxford University Press, 2009).
  * [50] R. O. Duda, P. E. Hart, and D. G. Stork, _Pattern classification_ (John Wiley & Sons, 2000).
  * [51] N. T. Islam, C. C. W. Lim, C. Cahall, J. Kim, and D. J. Gauthier, “Provably secure and high-rate quantum key distribution with time-bin qudits,” Sci. Adv. 3, e1701491 (2017).

sample


11institutetext: Université Paris-Saclay, CNRS, Institut d’Astrophysique
Spatiale, 91405, Orsay, France

# Image enhancement with wavelet-optimized whitening

F. Auchère Image enhancement with wavelet-optimized whiteningImage enhancement
with wavelet-optimized whitening E. Soubrié Image enhancement with wavelet-
optimized whiteningImage enhancement with wavelet-optimized whitening G.
Pelouze Image enhancement with wavelet-optimized whiteningImage enhancement
with wavelet-optimized whitening É. Buchlin Image enhancement with wavelet-
optimized whiteningImage enhancement with wavelet-optimized whitening

(Received 01 November 2022; accepted 02 December 2022)

###### Abstract

Context. Due to its physical nature, the solar corona exhibits large spatial
variations of intensity that make it difficult to simultaneously visualize the
features present at all levels and scales. Many general-purpose and
specialized filters have been proposed to enhance coronal images. However,
most of them require the ad hoc tweaking of parameters to produce subjectively
good results.

Aims. Our aim was to develop a general purpose image enhancement technique
that would produce equally good results, but based on an objective criterion.

Methods. The underlying principle of the method is the equalization, or
whitening, of power in the à trous wavelet spectrum of the input image at all
scales and locations. An edge-avoiding modification of the à trous transform
that uses bilateral weighting by the local variance in the wavelet planes is
used to suppress the undesirable halos otherwise produced by discontinuities
in the data.

Results. Results are presented for a variety of extreme ultraviolet (EUV) and
white light images of the solar corona. The proposed filter produces sharp and
contrasted output, without requiring the manual adjustment of parameters.
Furthermore, the built-in denoising scheme prevents the explosion of high-
frequency noise typical of other enhancement methods, without smoothing
statistically significant small-scale features. The standard version of the
algorithm is about two times faster than the widely used multiscale Gaussian
normalization (MGN). The bilateral version is slower, but provides
significantly better results in the presence of spikes or edges. Comparisons
with other methods suggest that the whitening principle may correspond to the
subjective criterion of most users when adjusting free parameters.

###### Key Words.:

Techniques: image processing – Methods: numerical – Sun: corona – Sun: UV
radiation – Sun: transition region

##  1 Introduction

Modern image sensors, like those used in the Extreme Ultraviolet Imager (EUI,
Rochus et al. 2020) on board the Solar Orbiter mission (Müller et al. 2020),
have dynamic ranges (ratio of the digital range to the read noise) of $2^{12}$
or more. As long as the observed scenes use the entire detector range, the
typical $2^{8}$ dynamic of monitors is insufficient to display the recorded
images unmodified without losing information (left panel of Fig. 1).
Therefore, most of the time the intensity values are rescaled to 8 bits prior
to display using a nonlinear function (or stretch), usually a power
$1/\gamma<1$ ($\gamma$-stretch) or a logarithm of the intensity (middle panel
of Fig. 1). It is worth noting that a square root transformation ($\gamma=2$)
is sometimes applied on board (e.g., on EUI) for it provides optimum intensity
sampling for images dominated by Poisson statistics (Nicula et al. 2005).
Nonlinear intensity transformations alone, however, are not sufficient to
bring out all the information present in the data. Filtering and enhancement
is a crucial step in the analysis of scientific images for it can reveal
features and phenomena that would otherwise remain undetected.

As illustrated in the right column of Fig. 1, the Fourier power spectrum of a
typical solar EUV image follows a power-law-like distribution, with more power
at low frequencies and less power at high frequencies. As per Parseval’s
theorem, the mean power of the Fourier spectrum is equal to the variance of
the signal. Thus, the larger (respectively lower) Fourier power at lower
(respectively higher) frequencies correspond to larger (respectively lower)
signal variance at larger (respectively smaller) spatial scales. The middle
column of Fig. 1, which is remapped using a logarithm intensity
transformation, indeed shows that the largest intensity variations come from
the global scales. However, while bringing out the faint large-scale features,
this simple intensity remapping produces a washed out image with poor local
contrast. Intensity variations in the active regions are mapped to a much
smaller range of output values than with the linear scaling of the left
column. This is why many enhancement methods are designed to reduce the signal
variance at low frequencies and to increase it at high frequencies.

Figure 1: Dynamic range and spatial scales in a solar EUV image. Left:
Fifteen-bit image taken by the Full Sun Imager (FSI) channel of EUI at
$17.4\text{\,}\mathrm{nm}$ on 2022 March 17 at 06:00:45 UT, displayed in
linear scale. Only the brightest features are visible. The profile is taken
along the row marked by the side ticks. Middle: Same image after logarithmic
scaling. Faint off-limb details are revealed, but the local contrast is
diminished. Right: Fourier power spectrum of the variance normalized image.
The profile is that of the bottom row. The larger the spatial scale, the more
power (and thus variability) in the image.

Specialized background subtraction methods (e.g., Morgan et al. 2006; Patel et
al. 2022) can be used to this end, but they are typically limited to the off-
disk emission. One of the earliest general-purpose methods, the unsharp mask,
is a simple high-pass filter that was developed to improve the quality of
photographic prints by analog composition of a negative with a blurred
positive (Yule 1944; Schreiber 1970; Levi 1974). Numerically, the convolution
kernel is generally taken to be Gaussian, but the frequency response of the
filter can be tailored to arbitrary shapes in the Fourier space. However,
Fourier filtering techniques are limited in that they ignore the local
variations of the image properties, unless the input is segmented into
independent tiles.

On the contrary, multiscale approaches are able to modify the spectral content
differently at different positions in the input image. The detail layer of the
multiscale Gaussian normalization (MGN, Morgan & Druckmüller 2014) is
essentially a weighted sum of unsharp-masked images of increasing Gaussian
kernel widths, normalized to their local standard deviations and remapped to
the [-1, 1] range. MGN is successful in equalizing the variances over a range
of scales, but also enhances the high-frequency noise in doing so (see bottom
right panel of Fig. 5), a common issue with image sharpening techniques.

The decomposition of images into scales using a wavelet transform offers the
possibility to simultaneously denoise (Starck & Murtagh 1994; Murtagh et al.
1995) and modulate the output power spectrum. The wavelet packet equalization
scheme of Stenborg & Cobelli (2003); Stenborg et al. (2008) consists in a two-
level à trous wavelet decomposition, followed by local noise reduction and
weighted synthesis.

Enhancement methods typically include free parameters that need to be adjusted
by the user: the kernel width for the unsharp mask; the kernel widths,
variance regularization function, and synthesis weights in MGN; and the noise
threshold and synthesis weights in Stenborg & Cobelli (2003). This provides
flexibility in tuning the output, with the drawback of subjectivity in
deciding what good results are.

In the present article we propose a transform that whitens the power spectrum,
resulting in equal variance at all spatial frequencies (section 2). While not
the only possible choice of normalization, a comparison of the results
(section 3) with other methods suggests that images considered to be well
enhanced tend to have a white spectrum. We summarize our work in section 4.

##  2 Method

###  2.1 À trous transform

Figure 2: A trous wavelet transform and coefficients whitening. Two top rows:
Chained convolutions $c_{s}(k)$ of the input image $c_{0}(k)$ with à trous
kernels dilated by a factor of two at each step. The left half of each panel
corresponds to regular convolutions; the right half to bilateral filtering.
Two middle rows: Corresponding wavelet coefficients
$w_{s+1}(k)=c_{s}(k)-c_{s+1}(k)$, that is the difference between successive
convolutions or bilateral filterings. The last coefficient is the result of
the last convolution or bilateral filtering. The original image ($c_{0}(k)$,
left) is the sum of all coefficients. The bilateral transform preserves the
edges in the coarse images (top rows). Two bottom rows: Whitened coefficients,
i.e., normalized to the local power. The filtered image (left) is the sum of
all whitened coefficients. For visualization purposes, the convolutions and
the coefficients are displayed using two different color scales.

We summarize below only those aspects of the à trous transform necessary to
understand the proposed method. We refer the reader to Holschneider et al.
(1989); Shensa (1992); Starck & Murtagh (2002) and references therein for
material on wavelet theory and the à trous transform in particular. See also
Stenborg & Cobelli (2003) for a concise introduction and application to images
of the solar corona.

The wavelet transform of a discrete signal $f(l)$ is defined by

$w_{s}(k)=\frac{1}{\sqrt{s}}\sum_{l}\psi^{*}\left(\frac{l-k}{s}\right)f(l),$
(1)

where $s$ is a scaling factor. In the case of the à trous transform, the real-
valued wavelet $\psi$ satisfies the relationship

$\frac{1}{2}\psi\left(\frac{x}{2}\right)=\phi(x)-\frac{1}{2}\phi\left(\frac{x}{2}\right)$
(2)

with $\phi(x)$ a scaling function that obeys the dilation equation

$\frac{1}{2}\phi\left(\frac{x}{2}\right)=\sum_{l}h(l)\,\phi(x-l),$ (3)

where $h$ is a low-pass filter kernel. We can then demonstrate (Shensa 1992)
that the à trous wavelet transform can be computed by the following iteration,
initialized with the input image $I(k)=c_{0}(k)$:

$\displaystyle c_{s+1}(k)$
$\displaystyle=\sum_{l\in\Omega}c_{s}(k+2^{s}l)\,h(l)=h_{s}(k)*c_{s}(k),$ (4)
$\displaystyle w_{s+1}(k)$ $\displaystyle=c_{s}(k)-c_{s+1}(k).$ (5)

Here $\Omega$ denotes the support of the compact kernel $h$ and $*$ the
convolution product. The set $\\{w_{1}(k)\ldots w_{S}(k),c_{S}(k)\\}$ forms
the wavelet transform of the input image, with $S$ the largest scale
compatible with the data. At scale $s$, the data values weighted by $h(l)$ are
separated by $2^{s}$, which is equivalent to a convolution with a sparse (with
holes, à trous in French) kernel $h_{s}(k)$. We note that this iteration lends
itself naturally to an implementation in which successive scales are
recursively computed by convolutions of the compact kernel with interleaved
subsamplings of one point out of two on each axis.

The successive convolutions form increasingly smooth representations
$c_{s}(k)$ of the data (left half images in the first two rows of Fig. 2). The
input image $c_{0}(k)$ used in this example was acquired with the
$17.4\text{\,}\mathrm{nm}$ high-resolution EUV channel of EUI
(HRI${}_{\text{EUV}}$, top left panel of Fig. 3) on 2022 March 17 at 04:02:57
UT, when Solar Orbiter was at 0.38 AU from the Sun. It is the one closest in
time to the FSI image of Fig. 1. Its wavelet transform with $S=8$ is
represented in the left half images of the third and fourth rows. The input
image can be synthesized with

$c_{0}(k)=c_{S}(k)+\sum_{s=1}^{S}w_{s}(k).$ (6)

In the remainder of this paper we use a basic cubic spline for the scaling
function $\phi$, which corresponds (Equation 3) to the one-dimensional kernel
$h=\left[1,4,6,4,1\right]/16$, its n-dimensional counterparts being obtained
by iterating $h^{T}h$. Other choices are possible, but the $B_{3}$-spline is
smooth, is continuously differentiable, and has minimum curvature and a
compact support, which makes it both robust to oscillations and
computationally efficient (see, e.g., Unser 1999). We use mirror boundary
conditions so that $c(k+N)=c(N-1-k)$ with $N$ the number of pixels on a given
axis.

The wavelet coefficients (third and fourth row of Fig. 2) can be manipulated
to achieve a variety of results after synthesis. Denoising will be treated in
section 2.3. In the next section, we discuss multiscale contrast enhancement.

###  2.2 Whitening

Figure 3: Three levels of WOW processing. Top left: Original
$\gamma$-stretched ($\gamma=2.4$) HRI${}_{\text{EUV}}$ image. The inset (top
left quarter) shows an enlargement of the area within the white square. Top
right: Output of the WOW algorithm. Bottom right: Same with added denoising.
Bottom right: Output of the edge-aware denoised version of the algorithm. For
all images the gray scale maps to the central 99.8% of values. An animated
version of this figure is available online.

It is possible to enhance features at a given scale, or more generally to
modify the relative importance of different scales, by multiplying the
corresponding coefficients $c_{k}^{s}$ with ad hoc weights $\beta_{s}$ at the
image synthesis stage, so that Equation 6 becomes

$\overline{c_{0}}(k)=c_{S}(k)+\sum_{s=1}^{S}\beta_{s}\,w_{s}(k).$ (7)

It is possible, for example, to suppress large-scale gradients by setting the
last weight to zero. This approach is used in the wavelet packet equalization
method of Stenborg & Cobelli (2003). While it yields interesting results, it
also produces undesirable ringing around sharp and bright features (see, e.g.,
the right panel of their Fig. 4). Nonetheless, building on this idea, we
wondered whether the weights, instead of being arbitrary, could be dictated by
the data content to produce an optimal output. Figure 1 suggests that a
possible scheme consists in equalizing the power at all scales (i.e.,
whitening the power spectrum). However, using one global weight per scale
ignores the information contained in the spatial variations of power in the
wavelet spectrum. In order to equalize the spectrum in both the spectral and
spatial dimensions, we normalize the coefficients using an estimate of the
local mean power $\overline{P_{s}}(k)$ obtained by convolution of the wavelet
power by the à trous kernel. The whitened coefficients $\overline{w_{s}}(k)$
are given by

$\overline{w_{s}}(k)=\frac{w_{s}(k)}{\sqrt{\overline{P_{s}}(k)}},\quad\mathrm{with}\quad\overline{P_{s}}(k)=w_{s}(k)^{2}*h_{s}(k).$
(8)

Unweighted synthesis of the whitened coefficients results in the wavelet-
optimized whitening (WOW) enhanced image:

$\overline{c_{0}}(k)=\overline{c_{S}}(k)+\sum_{s=1}^{S}\overline{w_{s}}(k).$
(9)

We note the similarity between Equation 8 and Equations 1 and 2 of Morgan &
Druckmüller (2014): both represent convolutions of the original image by
variable-sized kernels normalized to a quantity analogous to a local standard
deviation. However, the numerator in Equation 1 of Morgan & Druckmüller (2014)
is an unsharp mask (i.e., the difference between the original image and its
convolution by a Gaussian), while in Equation 8 it is an à trous wavelet
coefficient (i.e., the difference between two successive chained convolutions
by a scaling function). As we show in section 3.2, this difference between the
two methods results in significantly different output.

Figure 3 shows a $\gamma$-stretched high-resolution $17.4\text{\,}\mathrm{nm}$
EUI image (top left), and the corresponding output of the WOW process (top
right). Only the brightest regions of the input image would be visible without
$\gamma$-stretching, as in the left panel of Fig. 1. Multiscale whitening
tends to equalize the variance at all spatial frequencies and locations, which
attenuates the large-scale variations and comparatively reinforces the small-
scale structures. The large-scale contrast is reduced, while the small-scale
contrast is increased. The enlargement of the region in the white square shows
that the noise, which behaves like a small-scale feature, is amplified in the
process, especially in the originally faint regions where the noise
contributes the most to the signal variance. However, the à trous wavelet
transform offers a powerful framework for denoising.

###  2.3 Denoising

Table 1: Gain $g$ and read noise $r$ for the detectors used in this paper (Equation 14). Instrument | Gain | Read noise | Reference
---|---|---|---
| [DN/ph] | [DN] |
FSI304 | 3.88 | 1.50 | Gissot (2022)
HRI${}_{\text{EUV}}$ | 5.28 | 1.50 | Gissot (2022)
AIA174 | 1.12 | 1.15 | Boerner et al. (2012)
LASCO C2 | 0.07 | 0.3 | Brueckner et al. (1995)

The data can be denoised by attenuation of the wavelet coefficients that are
not statistically significant, followed by synthesis (Starck & Murtagh 1994;
Murtagh et al. 1995). This is achieved by comparing the $w_{s}(k)$ with the
amplitudes expected in the presence of noise only. In the case of Gaussian
white noise of unit standard deviation, the amplitudes at scale $s$ are
Gaussian-distributed with standard deviation $\sigma_{s}^{1}$, which can
easily be estimated numerically (Starck & Murtagh 2002). In the presence of
noise of standard deviation $\sigma$, the significant coefficients are those
for which

$\left|w_{s}(k)\right|>n_{s}\,\sigma_{s},\quad\mathrm{with}\quad\sigma_{s}=\sigma\,\sigma_{s}^{1},$
(10)

where $n_{s}$ sets the chosen significance level at scale $s$. The
coefficients can then be weighted by a function
$\alpha(n_{s}\,\sigma_{s},w_{s}(k))$ of their values relative to
$n_{s}\,\sigma_{s}$. Hard thresholding corresponds to

$\alpha(n_{s}\,\sigma_{s},w_{s}(k))=\cases{1}&\text{if}\quad\left|w_{s}(k)\right|>n_{s}\,\sigma_{s}\\\
0\text{otherwise}{},$ (11)

but tends to produce ringing in the synthesized image. There are several
possibilities for smooth soft-thresholding functions. We chose to weight the
coefficients by their probability of chance occurrence, so that

$\alpha(n_{s}\,\sigma_{s},w_{s}(k))=\mathrm{erf}\left(\frac{\left|w_{s}(k)\right|}{n_{s}\,\sigma_{s}}\right).$
(12)

In addition to the Gaussian component corresponding to the read and thermal
noises within the sensor, the data also contains the Poisson-distributed
component of the photon shot noise. The variance of the Poisson component is
proportional to the intensity and is thus variable across the image. In the
case of pure Poisson statistics, the transform

$T[I(k)]=\sqrt{I(k)+\frac{3}{8}}$ (13)

results in data with Gaussian noise of unit variance (Anscombe 1948) for
numbers of counts greater than about ten. It was generalized by Murtagh et al.
(1995) for the presence of additional Gaussian noise in the input.

Following Starck et al. (1997); Stenborg & Cobelli (2003), the alternative
approach adopted in this paper to account for the spatial dependence of the
noise consists in replacing $\sigma_{s}$ in Equations (10)–(12) by a spatially
variable $\sigma_{s}(k)=\sigma(k)\,\sigma_{s}^{1}$. The local standard
deviation of the noise $\sigma(k)$ can be estimated from the first wavelet
coefficients (Starck et al. 1997; Stenborg & Cobelli 2003) or, preferably, if
the detector gain $g$ and read noise standard deviation $r$ are known, with

$\sigma(k)=\sqrt{gI(k)+r^{2}}.$ (14)

In the end, for combined denoising and whitening, equation 8 becomes

$\overline{w_{s}}(k)=\frac{\alpha(n_{s}\,\sigma_{s}(k),w_{s}(k))\,w_{s}(k)}{\sqrt{\overline{P_{s}}(k)}};$
(15)

the expressions for $\overline{P_{s}(k)}$ and the synthesis (Equation 9)
remain unchanged.

The bottom left panel of Fig. 3 shows the denoised whitening of the top left
input image, to be compared with the regular whitening (Section 2.2) of the
top right panel. The detector gain and read noise (Equation 14) for
HRI${}_{\text{EUV}}$ are given in Table 1. Noise being inherently high-
frequency, it is usually significant only in the first few scales. In this
example, we used $n_{s}=\\{5,2,1\\}$ for the first three scales and $n_{s}=0$
for the others. This choice is arbitrary; larger values yield smoother output
at the expense of losing fainter details.

###  2.4 Edge-aware transform

In the à trous wavelet transform each convolution by the scaling function
(Equation 4) produces halos around sharp peaks and edges. These halos are
reinforced as genuine features by the whitening process, which tends to
produce glows and/or gradient reversals. As proposed by Hanika et al. (2011),
this unwanted side effect can be mitigated by replacing the convolution by
bilateral filtering in Equation (4). Range filtering is analogous to spatial
filtering, but instead of being a function of spatial distance, the kernel
weights are a function of the disparity of values. The bilateral filter
(Tomasi & Manduchi 1998) combines both spatial and range filtering, so that
Equation (4) becomes

$c_{s+1}(k)=\frac{\sum_{l\in\Omega}c_{s}(k+2^{s}l)\,h(l)\,b(k,k+2^{s}l))}{\sum_{l\in\Omega}{h(l)\,b(k,k+2^{s}l)}},$
(16)

where the denominator ensures unit normalization. A classical choice for the
range filter $b$ is a Gaussian

$b(k,k+2^{s}l)=\exp{\left(-\frac{1}{2}\left(\frac{c_{s}(k)-c_{s}(k+2^{s}l)}{\nu_{s}}\right)^{2}\right)}.$
(17)

Following the idea of equalization, we use for $\nu_{s}(k)^{2}$ the local
variance at scale $s$, so that the range weights do not depend upon the
absolute values of the wavelet coefficients. For computational efficiency,
this is approximated by

$\nu_{s}(k)^{2}=h_{s}(k)*c_{s}(k)^{2}-\left(h_{s}(k)*c_{s}(k)\right)^{2}.$
(18)

The resulting edge-aware transform is illustrated in the right halves of the
panels of Fig. 2. Compared to the regular transform, the edges are preserved
in the coarse images (top two rows) by the bilateral filtering, and the power
is correspondingly reduced in the coefficients (middle two rows). For
denoising with the edge-aware transform, the corresponding standard deviations
$\sigma_{s}^{1}$ of Gaussian white noise were estimated numerically. The
bottom right panel of Fig. 3 shows the output of the denoised edge-aware WOW
algorithm for the top left input image. With the same set of significance
levels $n_{s}$ as for the regular transform (bottom left panel), the output is
sharper, but denoising is less pronounced. An animated version of Fig. 3
containing 899 frames acquired on 2022 March 17 from 03:18:00 to 04:02:57 UT
is available online.

##  3 Results

###  3.1 Three examples

Figure 4: Examples of WOW processed images. For each row, the left image is
the original $\gamma$-stretched with $\gamma=2.4$, the right image is the
output of the edge-aware denoised version of the algorithm. Top: Crop of the
$30.4\text{\,}\mathrm{nm}$ FSI image taken on 2022 March 21 at 06:30:15 UT.
Middle: $17.1\text{\,}\mathrm{nm}$ AIA image taken on 2012 August 31 at
19:44:35 UT. Bottom: LASCO C2 image taken on 1998 June 2 at 13:31:06 UT.

Figure 4 shows three examples of images processed with the edge-aware version
of the WOW algorithm. The original $\gamma$-stretched ($\gamma=2.4$) images
are in the left column and the linearly stretched processed images, all
obtained using the same denoising significance levels $n_{s}=\\{5,3,1\\}$, are
in the right column. The gain and read noise for the different detectors are
given in Table 1. For all images, the black and white points of the linear
grayscale are mapped respectively to the 0.1 and 99.9 percentiles. In the FSI
$30.4\text{\,}\mathrm{nm}$ image (top row) the rapid off-disk intensity fall-
off masks the faint outer coronal structures (see also Fig. 1). In the top
right image, the whitening of the wavelet spectrum enhances the local contrast
and reveals the fainter extensions. The method also accentuates imperfections
in the calibration (vertical bands and pattern at the top and bottom edges)
because they are spatially correlated, and therefore are not considered as
noise by the algorithm. The Atmospheric Imaging Assembly (AIA, Lemen et al.
2012) image (middle row) can be compared to Fig. 2 of Druckmüller (2013)
processed with noise adaptive fuzzy equalization (NAFE, see next section). The
two images are visually similar, with WOW producing as sharp and contrasted an
output without the ad hoc adjustment of parameters nor the addition of a
$\gamma$-stretched component. In images from white-light coronagraphs (like
that in the bottom left of Fig. 4, from LASCO C2, Brueckner et al. 1995), the
dust corona and stray light form a diffuse background that dominates the
K-corona. By equalizing the power across spatial frequencies, this component
is attenuated with respect to smaller-scale features resulting, in the bottom
right image, in the enhancement of the twisted filamentary structure of the
erupting prominence. Comparison with the output of the wavelet packets
equalization method (Fig. 4 of Stenborg & Cobelli 2003) shows that WOW avoids
the ringing around the bright and sharp features, and also enhances the
disturbances in the surrounding K-corona. As already mentioned, the algorithm
treats coherent features irrespective of their solar or instrumental nature.
For optimal results, artifacts like the inner diffraction rings or the long
curved streaks should be either masked or calibrated out (Morgan 2015; Lamy et
al. 2020, 2022) prior to enhancement.

###  3.2 Comparison with NAFE and MGN

Figure 5: Comparison of the WOW, NAFE, and MGN algorithms. Top left: Original
$\gamma$-stretched ($\gamma=2.4$) HRI${}_{\text{EUV}}$ image. The inset (top
left quarter) shows an enlargement of the area within the white square. Top
right: Output of the edge-aware denoised version of the algorithm. Middle row:
NAFE (left) and corresponding detail layer (right). Bottom row: MGN (left),
and corresponding detail layer (right). In NAFE and MGN, the
$\gamma$-stretched original image (top left) is added to their respective
detail layers (right column). The left column thus corresponds to: no
enhancement, added NAFE detail layer, added MGN detail layer. The standard WOW
output does not include an arbitrary $\gamma$-stretched component. An animated
version of this figure is available online.

We chose to compare the output of WOW with that of NAFE (Druckmüller 2013) and
MGN (Morgan & Druckmüller 2014), two algorithms commonly used in the solar
physics community. The output $O(k)$ of both NAFE and MGN is the weighted sum
of the $\gamma$-stretched input $I(k)$ and of the filtered image $F(k)$
proper:

$O(k)=w\left(\frac{I(k)-a_{0}}{a_{1}-a_{0}}\right)^{1/\gamma}+(1-w)F(k).$ (19)

Here $a_{0}$ and $a_{1}$ are respectively the minimum and maximum input values
considered, and $w$ is a user-defined weight set by default to $0.2$ for NAFE
and $0.7$ for MGN. In order to have comparable output, we set $\gamma=2.4$ for
both NAFE and MGN. The detail layers $F(k)$ of NAFE and MGN depend on
additional arbitrary parameters. For NAFE we use the width $N$ of the
neighborhood for local histogram equalization and the denoising coefficient
$\sigma$ (set to 129 and 12 respectively in the following). For MGN we use the
scaling parameter $k$ of the $\arctan$ function, the synthesis weights (all
set to 1 by default), the number of scales (6 by default), and the widths of
the Gaussian kernels, which by default follow a geometric progression,
similarly to the à trous transform. Conversely, in WOW, the values of most
parameters are dictated by either the use of the à trous transform (the
progression of scales) or the whitening principle (the number of scales, the
absence of synthesis weights). The only free parameters are the denoising
significance levels.

As a test case, we picked an HRI${}_{\text{EUV}}$ image recorded on 2022 April
2 at 09:30:55 UT while Solar Orbiter was at 0.38 AU from the Sun. The
$\gamma$-stretched layer ($\gamma=2.4$) of NAFE and MGN is shown in the top
left panel of Fig. 5. An animated version containing 450 frames acquired from
09:19:15 to 10:34:05 UT is available online. For NAFE (middle row) and MGN
(bottom row), the output is on the left, and the corresponding detail layer on
the right. WOW does not have a $\gamma$-stretched layer, and thus its output
(top right panel) must be compared with the detail layers of NAFE and MGN
(right column). For all images the gray scale maps to the central 99.8% of
values.

NAFE and MGN produce a somewhat similar output, with NAFE being more
contrasted and MGN emphasizing the small scales to a greater degree. The
detail layer of MGN presents the typical dark fringes around bright features
produced by unsharp masking. Despite the similarities mentioned in Section 2,
WOW produces a significantly more contrasted output than MGN and no ringing,
with notably less noise. Judgment of the merits of the enhancement methods can
be subjective when free parameters are involved. The output of WOW is both as
sharp as MGN and as contrasted as NAFE, which suggests that the whitening
principle may correspond to what is qualitatively considered to be good
results.

###  3.3 Performance

Table 2: Comparison of execution times (in seconds) for NAFE (on 48 cores), MGN, and WOW, for the images presented in this paper. |  |  |  | WOW
---|---|---|---|---
Image | size | NAFE | MGN | standard | edge-aware | edge-aware (6 scales)
HRI${}_{\text{EUV}}$ 2022-03-17T04:02:57 | $2048\times 2048$ | 98.8 | 7.3 | 2.4 | 11.4 | 7.6
FSI${}_{\text{304}}$ 2022-03-21T06:30:15 | $3040\times 3072$ | 936.9 | 11.4 | 10.9 | 37.3 | 16.9
AIA${}_{\text{174}}$ 2012-08-31T19:44:35 | $4096\times 4096$ | 813.3 | 49.8 | 13.6 | 53.1 | 27.3
LASCO C2 1998-06-02T13:31:06 | $1024\times 1024$ | 190.7 | 1.3 | 0.5 | 2.3 | 1.9
HRI${}_{\text{EUV}}$ 2022-04-02T09:30:55 | $2048\times 2048$ | 100.4 | 5.8 | 2.4 | 11.3 | 7.4
FSI${}_{\text{174}}$ 2022-03-17T06:00:45 | $3072\times 3072$ | 1330.5 | 13.8 | 10.9 | 37.7 | 17.0

WOW was implemented in Python using `open-cv` bindings to compute the
convolutions; the code is available on GitHub111https://github.com/frederic-
auchere/wavelets/. We compare the performance to a parallelized NAFE
implementation in Python222https://git.ias.u-psud.fr/ebuchlin/aia-
movie/-/blob/master/medocimage/nafe.py validated by comparison with the
original implementation, and to the `sunkit-image`333sunkit-image
(https://docs.sunpy.org/projects/sunkit-image/) is a SunPy (The SunPy
Community et al. 2020) affiliated package. implementation of MGN. The results
of performance tests on a 2.40 GHz Intel Xeon Silver 4214R CPU are summarized
in Table 2. The execution time of WOW and MGN depends on the number of scales,
which is set to six by default in MGN. In WOW, the number of scales is set to
be $\ln(N)/\ln(2)$. For comparison with MGN, we therefore also include
execution times for the edge-aware version of WOW limited to six scales. The
operations required by the two algorithms are very similar, hence the similar
performance. However, the bilateral filter makes the edge-aware version of WOW
significantly slower than the standard one. Compared to NAFE, it is much
faster in all cases while providing equivalent results.

###  3.4 Back to subjectivity

WOW was developed to identify an objective criterion for image enhancement
that leads to a method free of arbitrary parameters. While the approach seems
successful, the whitening of the power spectrum is not necessarily the unique
viable choice. In addition, manual adjustment of the output may be desirable.
Naturally, a weighted $\gamma$-stretched input image can be added to the WOW
image, as in NAFE and MGN. However, any other transfer function can be used to
rescale the large dynamic range of the input onto a restricted range of output
values. A logarithmic scaling is commonly used for solar EUV images. Lupton et
al. (2004) proposed to stretch the data using
$\operatorname{arcsinh}(x/\beta)$, which is linear for $x\ll 1$ and
logarithmic for $x\gg 1$, with $\beta$ a free parameter. Synthesis weights can
also be applied at each scale (Equation 7) to provide enhanced control over
the result.

##  4 Summary

We introduced a new image enhancement algorithm that produces sharp and
contrasted outputs from a variety of images, without arbitrarily determined
coefficients nor the addition of a gamma-stretched image. While it is possible
to introduce ad hoc parameters at both the decomposition and the synthesis
stages, the power spectrum whitening principle at the heart of the method,
although arbitrary in itself, dictates that these parameters be set to one,
yielding an optimum parameter-free scheme. The optional denoising step still
requires significance levels to be set manually; however, based on the work of
Batson & Royer (2019), it may be possible to determine optimal values
automatically. The regular variant is significantly faster than MGN and
produces more contrasted results. The sensitivity of the nominal scheme to
spikes and edges is mitigated in the edge-aware variant, which is
significantly slower, but remains practical on a typical laptop computer.

###### Acknowledgements.

Solar Orbiter is a space mission of international collaboration between ESA
and NASA, operated by ESA. The EUI instrument was built by CSL, IAS, MPS,
MSSL/UCL, PMOD/WRC, ROB, LCF/IO with funding from the Belgian Federal Science
Policy Office (BELSPO/PRODEX PEA C4000134088); the Centre National d’Etudes
Spatiales (CNES); the UK Space Agency (UKSA); the Bundesministerium für
Wirtschaft und Energie (BMWi) through the Deutsches Zentrum für Luft- und
Raumfahrt (DLR); and the Swiss Space Office (SSO). AIA is an instrument on
board SDO, a mission for NASA’s Living With a Star program. SOHO/LASCO data
are produced by a consortium of the Naval Research Laboratory (USA), Max-
Planck-Institut für Aeronomie (Germany), Laboratoire d’Astronomie (France),
and the University of Birmingham (UK). SOHO is a project of international
cooperation between ESA and NASA. Gabriel Pelouze’s work was supported by a
CNES post-doctoral grant. This work used data provided by the MEDOC data and
operations centre (CNES / CNRS / Université Paris-Saclay),
http://medoc.ias.u-psud.fr/.

## References

  * Anscombe (1948) Anscombe, F. J. 1948, Biometrika, 35, 246
  * Batson & Royer (2019) Batson, J. & Royer, L. 2019, in Proceedings of Machine Learning Research, Vol. 97, Proceedings of the 36th International Conference on Machine Learning, ed. K. Chaudhuri & R. Salakhutdinov (PMLR), 524–533
  * Boerner et al. (2012) Boerner, P., Edwards, C., Lemen, J., et al. 2012, Sol. Phys., 275, 41
  * Brueckner et al. (1995) Brueckner, G. E., Howard, R. A., Koomen, M. J., et al. 1995, Sol. Phys., 162, 357
  * Druckmüller (2013) Druckmüller, M. 2013, ApJS, 207, 25
  * Gissot (2022) Gissot, S. 2022, personal communication
  * Hanika et al. (2011) Hanika, J., Dammertz, H., & Lensch, H. 2011, Computer Graphics Forum, 30, 1879
  * Holschneider et al. (1989) Holschneider, M., Kronland-Martinet, R., Morlet, J., & Tchamitchian, P. 1989, in Wavelets. Time-Frequency Methods and Phase Space, ed. J.-M. Combes, A. Grossmann, & P. Tchamitchian, 286
  * Lamy et al. (2020) Lamy, P., Llebaria, A., Boclet, B., et al. 2020, Sol. Phys., 295, 89
  * Lamy et al. (2022) Lamy, P. L., Gilardy, H., & Llebaria, A. 2022, Space Sci. Rev., 218, 53
  * Lemen et al. (2012) Lemen, J. R., Title, A. M., Akin, D. J., et al. 2012, Sol. Phys., 275, 17
  * Levi (1974) Levi, L. 1974, Computer Graphics and Image Processing, 3, 163
  * Lupton et al. (2004) Lupton, R., Blanton, M. R., Fekete, G., et al. 2004, PASP, 116, 133
  * Morgan (2015) Morgan, H. 2015, ApJS, 219, 23
  * Morgan & Druckmüller (2014) Morgan, H. & Druckmüller, M. 2014, Sol. Phys., 289, 2945
  * Morgan et al. (2006) Morgan, H., Habbal, S. R., & Woo, R. 2006, Sol. Phys., 236, 263
  * Müller et al. (2020) Müller, D., St. Cyr, O. C., Zouganelis, I., et al. 2020, A&A, 642, A1
  * Murtagh et al. (1995) Murtagh, F., Starck, J. L., & Bijaoui, A. 1995, A&AS, 112, 179
  * Nicula et al. (2005) Nicula, B., Berghmans, D., & Hochedez, J.-F. 2005, Sol. Phys., 228, 253
  * Patel et al. (2022) Patel, R., Majumdar, S., Pant, V., & Banerjee, D. 2022, Sol. Phys., 297, 27
  * Rochus et al. (2020) Rochus, P., Auchère, F., Berghmans, D., et al. 2020, A&A, 642, A8
  * Schreiber (1970) Schreiber, W. F. 1970, Pattern Recognition, 2, 117
  * Shensa (1992) Shensa, M. J. 1992, IEEE Transactions on Signal Processing, 40, 2464
  * Starck & Murtagh (1994) Starck, J.-L. & Murtagh, F. 1994, A&A, 288, 342
  * Starck & Murtagh (2002) Starck, J.-L. & Murtagh, F. 2002, Astronomical Image and Data Analysis (Springer-Verlag)
  * Starck et al. (1997) Starck, J.-L., Siebenmorgen, R., & Gredel, R. 1997, ApJ, 482, 1011
  * Stenborg & Cobelli (2003) Stenborg, G. & Cobelli, P. J. 2003, A&A, 398, 1185
  * Stenborg et al. (2008) Stenborg, G., Vourlidas, A., & Howard, R. A. 2008, ApJ, 674, 1201
  * The SunPy Community et al. (2020) The SunPy Community, Barnes, W. T., Bobra, M. G., et al. 2020, The Astrophysical Journal, 890, 68
  * Tomasi & Manduchi (1998) Tomasi, C. & Manduchi, R. 1998, in Sixth International Conference on Computer Vision (IEEE Cat. No.98CH36271), 839–846
  * Unser (1999) Unser, M. 1999, IEEE Signal Processing Magazine, 16, 22
  * Yule (1944) Yule, J. 1944, Photographic Journal, 84, 321