Signal averaging improves signal-to-noise in OCT images: But which approach works best, and when?

The high acquisition speed of state-of-the-art optical coherence tomography (OCT) enables massive signal-to-noise ratio (SNR) improvements by signal averaging. Here, we investigate the performance of two commonly used approaches for OCT signal averaging. We present the theoretical SNR performance of (a) computing the average of OCT magnitude data and (b) averaging the complex phasors, and substantiate our findings with simulations and experimentally acquired OCT data. We show that the achieved SNR performance strongly depends on both the SNR of the input signals and the number of averaged signals when the signal bias caused by the noise floor is not accounted for. Therefore we also explore the SNR for the two averaging approaches after correcting for the noise bias and, provided that the phases of the phasors are accurately aligned prior to averaging, then find that complex phasor averaging always leads to higher SNR than magnitude averaging. Published by The Optical Society under the terms of the Creative Commons Attribution 4.0 License. Further distribution of this work must maintain attribution to the author(s) and the published article’s title, journal citation, and DOI.

Introduction

Optical coherence tomography (OCT) [1] performs biomedical imaging at high speed and with high sensitivity. State-of-the art OCT systems provide data rates of ∼100 million pixels per second and frame rates of ∼100-200 frames per second and beyond [2,3]. OCT provides a distinctively high sensitivity of typically ∼100 dB, meaning that backscatter signals as weak as of a mirror reflection can still be detected [4-6]. At the same time, OCT enables imaging with a very high dynamic range spanning several tens of decibels between the strongest and the weakest signal in the image.

The strong performance of OCT in detecting weak signals can be improved even more by image processing. Noise in OCT images limits the detection capabilities for weakly scattering structures, and thus a variety of different approaches for reducing OCT image noise have been proposed [7-15]. The high acquisition speeds of modern OCT systems enable the improvement of the signal-to-noise ratio by averaging multiple signals, for instance by fusing OCT frames or even volumes quickly repeated at the same sample position. In recent years, several methods were proposed for improving the detection sensitivity of OCT by averaging the complex-valued OCT signals rather than just averaging their magnitudes [16-20]. In 2013, Szkulmowski and Wojtkowski published a thorough analysis of signals and noise subject to different averaging approaches [21]. In their analysis, the authors found a much stronger reduction of the noise floor by complex averaging as compared to magnitude averaging but also observed a heterogeneous outcome in terms of signal-to-noise performance for different imaging scenarios.

In this article, we set out to answer the question: Which signal averaging approach works best for improving signal-to-noise in OCT images, and when? We first introduce signals and noise in OCT (section 2.1) based on the analysis in Ref. [21]. We then analyze the signals and noise after magnitude and complex averaging, respectively, for a given pixel in an OCT image (section 2.2) and present simple expressions for the resulting signal-to-noise ratios. Next, we compare the somewhat surprising theoretical performance of the averaging schemes for different input signal levels and for different numbers of averaged signals and introduce a noise bias corrected signal-to-noise ratio (SNR) analysis (section 2.3). After substantiating our theoretical analysis with data from simulations (section 3.1) and experimentally acquired OCT data (section 3.2), we conclude the paper with a brief summary of the findings and their implications on actual OCT image processing (section 4). An overview of the terminology as well as a brief section on phase correction required for complex phasor averaging is provided in the appendix.

Analysis of signals and noise in OCT

OCT signals and noise

OCT signals can be described by complex phasors of the form where denotes the amplitude and the phase of the OCT signal at location and time . Here the amplitude represents the length of the phasor and the phase corresponds to the polar angle in the complex plane. The noise in OCT images can be specified by phasors too, namely by random phasors. In the complex plane, the real and imaginary parts of random phasors, and , can be described by normal distributions with zero mean and variance [22]. Hence, the complex noise signal is characterized by a binormal (or Beckmann) distibution in the complex plane [22,23]: This 2D Gaussian probability density function (PDF) describes the probability for a noise phasor of an image pixel to have the real part and the imaginary part . The probability density function of the noise amplitude is represented by the Rayleigh distribution (which can be derived by transforming Eq. (2) to polar coordinates and integrating over all angles ) [22] Analogous to , the PDF of the noise amplitude represents the probability of the noise amplitude to take specific values . Figure 1(a-d) shows an example of the complex representation of noise phasors in an OCT image and the corresponding Beckmann and Rayleigh distributions. Note that the Beckmann distribution is symmetric and centered at the origin while the Rayleigh distribution is skewed and only takes positive values. The mean amplitude and variance of a general PDF is given by and , respectively [22]. The mean amplitude and the variance of can be computed as

Fig. 1.

Complex phasor representation of noise and signals in OCT and their probability density functions. (a) Cartoon of Beckmann distribution of noise phasors around the origin of the complex plane. A representative phasor with real part and imaginary part is shown in green. (b) Complex OCT signals of 100 repeated noise measurements in the same pixel from real-world OCT data. (c) Histogram of noise amplitudes (1000 repeats, gray line) and Rayleigh PDF (red line) computed from the standard deviation of the Beckmann distribution in (b) by Eq. (3). (d) Histogram of the noise intensity (gray line) and PDF (red line) computed from in (b) by Eq. (9). (e) Cartoon of an OCT signal affected by noise. The green arrow represents the signal phasor. (f) Complex OCT signals of 100 repeated measurements of a weak reflection in the same pixel. (g) Histogram of signal amplitudes (1000 repeats, gray line) and Rice distribution (blue line) computed from the mean signal amplitude and using Eq. (6). (h) Histogram of the signal intensity (gray line) and PDF (blue line) computed from the mean intensity and by Eq. (10). The + in (b) and (f) indicates the origin of the complex plane. The PDFs in (c,d,g,h) were scaled to match the count levels of the respective histograms.

In every OCT image pixel, the detected signal is a combination of the actual OCT signal and the noise . By choosing the OCT signal to have and thus to point into the direction of the positive real axis (Fig. 1(e)) similar to Goodman [22], the PDF of the measured amplitudes incorporating both signal and noise is represented by a Rice distribution where B0 is a modified Bessel function of the first kind and zero order. An example of the PDF of the signal amplitudes in presence of noise seen in an OCT image is shown in Fig. 1(g). Note that the signal gives rise to a somewhat spread probability accumulation at a high amplitude level. The mean amplitude and variance in an image pixel are observed as where denotes the Laguerre polynomial of degree . yields steadily increasing, positive values for increasingly negative arguments .

OCT amplitude data are usually squared in order to display a quantity proportional to sample reflectivity, , henceforth called intensity. Therefore it makes sense to investigate the behavior of the PDFs of the squared signal and noise amplitude data. For this purpose, a variable transform of is performed such that the PDFs are rendered into . The resulting PDFs for noise () and signal intensities () are where . Note that is a negative exponential probability function whereas is a combination of a negative exponential and a monotonically increasing Bessel function. Examples of the above PDFs for noise and signal intensities are shown in Figs. 1(d) and (h), respectively. The mean values and variances for noise and signal intensities, respectively, are given by By comparing Eqs. (11) and (12), one can observe that the mean noise intensity equals the standard deviation . Therefore, for the case of both non-zero signal and noise, the observed mean intensity corresponds to and the intensity variance to . Next, we will investigate the effect of averaging the magnitudes of the phasors as well as the effect of averaging the complex phasors on OCT intensity data.

Averaging magnitudes and phasors

For averaging OCT data, usually the absolute values (magnitudes) of the signals are used. For some applications, the complex signals have been exploited for adding or increasing image contrast in one way or another [31-36]. Here, we systematically analyze the effect of averaging magnitude and complex signals in OCT images.

Magnitude averaging

Magnitude averaging uses the absolute values of spatially and/or temporally separated signals (e.g., pixels in a 2D or 3D kernel or repeated measurements at the same spatial pixel location but at different time points): The PDF representing the average of data points can be computed numerically by convolving the PDF of a single data point (N − 1)-times [21]. Unfortunately, no closed form expressions are available for the distribution describing an N-fold average of the above PDFs and only some approximations and bounds have been derived [37]. Nevertheless, the mean intensity values and intensity variances can be calculated for Rayleigh distributed noise and Rician signals affected by random noise, respectively, as The above expressions show that the mean noise level as well as the observed mean signal level remain unchanged upon averaging and retain the mean intensity values calculated for single data points. However, the variances of the observed signal and noise, and , are reduced by and thus reducing the standard deviations and by a factor .

Complex averaging

In contrast to magnitude averaging, complex averaging also includes the phase information into the averaging process. This inclusion exploits the fact that the noise phasors (in absence of a signal) randomly fluctuate about the origin according to a 2D Gaussian PDF with a standard deviation of , see Eq. (2). By averaging spatially and/or temporally independent complex noise phasors, the standard deviation of the Beckmann distribution in Eq. (2) is reduced by a factor [22]. This reduction of the noise amplitude variance translates into a reduction of both the mean noise intensity as well as of its variance: Compared to the noise performance after magnitude averaging (Eqs. (16) and (17)), which maintains the mean noise level and reduces the noise intensity variance by , complex averaging leads to a decrease of the noise floor intensity and to a smaller noise intensity variance.

Recalling that we initially assumed that the signal vector was pointing in the direction of the positive real axis (i.e. Fig. 1(e)), the precondition for averaging signals in a complex fashion is that the phase of remains constant at . (Likewise, a different signal phase could be chosen, as long as the phases of all signals to be averaged are aligned at the same angle . The signal phase () is chosen here out of convenience to ensure a purely real signal which facilitates the calculations described above.) For such perfectly aligned signals in the presence of random noise, the following mean signal intensity and intensity variance will be observed: Unlike magnitude averaging, which maintains the original mean signal intensity , the mean signal intensity does not stay constant but is reduced by after complex averaging. This signal reduction converges to an amount of for large . At the same time, also the noise intensity variance is reduced by an amount compared to magnitude averaging, which converges to a reduction by for large .

Signal-to-noise performance

The signal-to-noise ratio is the measure of choice to describe the influence of noise on signal measurements. An overview of methods for determining the SNR and in particular the sensitivity of an OCT system has recently been published by Agrawal et al. [24]. The SNR in OCT is typically defined as the ratio of the average signal intensity to the standard deviation of the noise intensity, [24-28] (or alternatively for averaged signals). As shown in the previous sections and visualized in Fig. 2, the measured signal intensity is the sum of the pure signal intensity and the noise floor . Historically, the bias caused by the noise floor is not specifically factored out from the SNR analysis, however some investigations of OCT image data did exclude the background [29,30]. Because of the multiple approaches for OCT data processing, it makes sense to ask, "how important is noise bias correction when measuring SNR?". In the following section, we will investigate the SNR for the total measured signal (i.e. including the contribution of the noise floor) with and without averaging. Subsequently, we will account for the signal offset caused by the noise floor and evaluate the SNR with noise bias correction.

Fig. 2.

Intensity signal and noise background in a schematic OCT depth profile. The noise floor is characterized by its average intensity and its variance . The measured OCT intensity signal consists of the pure signal intensity biased by the average noise level .

Signal-to-noise performance without noise bias correction

Using the signal and noise pairs in Eqs. (13) and (12), (18) and (17), and (22) and (21), and the relation in Eqs. (11) and (12), we find the following SNRs for single signals, magnitude averaged signals, and complex averaged signals, respectively: In short, averaging the magnitudes of signals improves the SNR by a factor whereas complex averaging of signals changes the SNR by a factor less .

The relative SNR improvements SNR/SNR1 and SNR/SNR are plotted for N up to 100 and for typical OCT image dynamics with SNR1 between 5 dB and 20 dB in Fig. 3. In particular for small , complex averaging provides a considerable SNR improvement over magnitude averaging. At the same time, the SNR1 dependence of complex averaging has a considerable impact on its performance: While for strong signals with SNR1 ≫ 1, complex averaging outperforms magnitude averaging by a factor , the signal-to-noise improvement becomes much less for weaker input signals.

Fig. 3.

Relative SNR improvement by signal averaging for strong input signals (without noise bias correction). (a) The ratios SNR/SNR1 and SNR/SNR1 are shown for from 1 to 100. SNR/SNR1 is plotted as a dash-dotted line, whereas the SNR1-dependent ratio SNR/SNR1 is plotted in rainbow colors for several SNR1 values between 5 dB and 50 dB. Note that SNR/SNR1 converges to an N-fold improvement. (b) The ratio is plotted for the spectrum of SNR1 values used in (a). Note that in particular for strong input signals with high SNR1, complex averaging outperforms magnitude averaging and converges to a -fold better SNR performance. As the SNR profiles converge for large SNRs, the curves in (a) and (b) start to overlap for values greater than 15 dB.

In the next section, we are going to explore the SNR enhancement of the two averaging approaches particularly for very weak signals on the order of the noise background, however still without correcting for the signal bias caused by the noise floor. Will magnitude averaging or complex averaging be superior in terms of recovering such small signals?

Averaging domains and borderline SNR

When the signal bias caused by the noise is not accounted for, the SNR performance averaging depends on the approach taken. Equations (25) and (26) revealed a dependence on the number of averaged signals, , and – for complex averaging – on the SNR of the input signals, SNR1. In this section, we are investigating this dependence for different numbers of averaged signals, , and for a great range of signals – from much smaller than the noise variance to several orders of magnitude greater. Recalling the identities for the noise variance and the observed signal intensity in Eqs. (12) and (13), we choose the quantity , which is identical to SNR1 − 1, as the benchmark for the input signal.

Figure 4(a) shows three plots of the signal-to-noise ratios calculated for magnitude averaging (red) and complex averaging (blue) where relatively small signals , 0.5 and 0.1 were chosen as the respective inputs. Note that for (left plot), complex averaging always yields greater SNR than magnitude averaging for . However, for and a small number of signals , magnitude averaging provides a more effective SNR improvement. Figure 4(b) plots three pairs of SNR profiles for a fixed number of averaged signals (, 10 and 100), this time as a function of the relative signal intensity . Again, for weak input signals, the superior performance of magnitude averaging can be observed, while complex averaging excels beyond a crossover point of . This relative borderline signal is plotted for up to 100 in Fig. 4(c) alongside the borderline input SNR1 which is the input SNR1 for which magnitude averaging and complex averaging perform equally well.

Fig. 4.

Influence of input signal level and number of averaged signals on the SNR (without noise bias correction). (a) SNR and SNR after magnitude and complex averaging of signals plotted for relative signal levels of (left), (middle), and (right), respectively. (b) SNR and SNR after magnitude and complex averaging of signals ranging from 0.01 through 10, plotted for averages of (left), (middle), and signals (right), respectively. Green arrows in (a) and (b) indicate the intercepts of the SNR profiles, i.e. the borderline SNR where magnitude and complex averaging perform equally. (c) Borderline plots of as well as SNR1 as described in Eq. (27).

An overview of the relative SNR performance SNR/SNR at various inputs and is shown as a heat map in Fig. 5. At first glance, two domains can be differentiated: For stronger signals and larger , complex averaging provides a greater SNR improvement than magnitude averaging (see the area in blue in Fig. 5). For signals comparable to or smaller than the noise level, magnitude averaging yields a better SNR enhancement (see the area in red). In between these two domains, the borderline with equal SNR performance of the two approaches is indicated in white, SNR/SNR = 1.

Fig. 5.

Relative SNR performance for magnitude and complex averaging (without noise bias correction). The heat map displays the ratio SNR/SNR in decibels for relative input signals ranging from -20 dB to +40 dB and up to a number of averaged signals . For small input signals below the noise level, magnitude averaging yields a better SNR improvement (red range), while complex averaging performs better for greater and stronger input signals (blue range). The borderline SNR where SNR = SNR separates these two domains (white plot).

SNR in absence of a signal (noise floor only)

An interesting scenario is posed by the calculation of the SNR when the signal intensity is zero. Then, using the relation from Eqs. (11) and (12), the SNRs for a single signal, after magnitude averaging and after complex averaging, respectively, become: An obvious disunity of the SNRs calculated for pure noise signal can be observed. This calls for a revised definition of the SNR, this time accounting for the signal bias imposed by the noise floor.

SNR analysis with noise bias correction

In the limit of a very weak OCT signal, the small meaningful signal contribution sits on top of a comparatively large noise floor. This noise bias is evident as the term related to in Eqs. (24–26) and is also visualized in Fig. 2. Thus, it makes sense to actually consider the impact of this noise bias in the SNR analysis – in particular for weak signals. In this section, we introduce and evaluate an SNR analysis with noise bias correction. A similar approach has for instance also been used by Makita et al. [39] and recently been modified by our group [40] to correct for noise contributions in PS-OCT images and improve the computation of the degree of polarization uniformity for weak signals. To some extent, this SNR with noise bias correction resembles the SNR definitions in Refs. [29,30] and has formal similarities with previous definitions of the contrast-to-noise ratios in Refs. [21,27].

The noise bias of the average OCT signals can be removed by subtracting the average noise intensity (, , and in Eqs. (11), (16) and (20), respectively) from the average signal intensity (, , and in Eqs. (13), (18) and (22), respectively). The noise bias corrected average signal intensities then read Note that the three corrected average signal intensities now match the pure signal intensity . Further, using the noise variances in Eqs. (12), (17) and (21), the respective noise bias corrected SNRs can be calculated as equates the quantity already known from the analyses in the previous sections. Even more strikingly, the noise bias corrected SNRs after magnitude and complex averaging now exhibit a simple proportionality of and with (Fig. 6), which also means that SNR′ will be zero if there is no signal for Eqs. (34–36). Using the noise bias correction, complex averaging thus always yields a -times better SNR than magnitude averaging, even for weak input signals and small .

Fig. 6.

Theoretical improvement of the signal-to-noise ratio SNR′ after noise bias correction plotted on (a) linear scales and (b) log scales. A - and N-fold improvement of the SNR′ of a single signal can be observed for magnitude and complex averaging, respectively. Note that, unlike for the noise-afflicted SNR calculations in Figs. 3 through 5, neither of the averaging approaches depend on the input signal strength.

Experimental validation

Numerical simulation

OCT signals and noise were simulated according to the probability density functions described in Eqs. (3) and (6), respectively. For the noise contribution, binormally distributed complex signals with similar standard deviations along both the real axis and the imaginary axis were generated. For the signal contribution, the binormal distribution was generated around a real-valued signal such that the phasor distribution was shifted from the origin of the complex plane to (see also Fig. 1(e)). A total of averaged signals were computed for every simulation run. For complex averaging, in every run to complex-valued phasors were averaged. For magnitude averaging, the absolute values were computed for every single noise and signal phasor prior to averaging. Finally, the respective SNRs and SNR′s were calculated from the average signal intensity and the variance of the noise signals for every .

Results of the simulations performed in MATLAB (R2014a, MathWorks) are shown in Fig. 7. To showcase the effect of averaging for strong and weak signals, SNR simulations for two distinct relative signal levels of and are presented in Fig. 7(a) and (b), respectively. Unlike magnitude averaging, complex averaging markedly decreased the signal intensity but at the same time also reduced noise much more by averaging. For the standard deviation of the noise fluctuations, a decrease by and was observed for magnitude and complex averaging, respectively. The resulting SNRs are shown in the four panels in Fig. 7. Here, for averaging up to 100-fold, SNR reveals a better SNR improvement for the strong input signal (Fig. 7(a)) when the noise bias is not corrected for, while magnitude averaging appears more effective for the weaker input signal in Fig. 7(b). After noise bias correction, however, no more dependence upon the input signal level can be observed, and complex averaging provides an N-fold increase in SNR′ whereas the SNR′ is only enhanced by a factor after magnitude averaging.

Fig. 7.

Simulation of the effect of averaging signals with relative strength in (a) and in (b). Shown are the averaged signal-to-noise ratios calculated from simulated phasors () alongside the corresponding theoretical plots (−). The SNRs without and with noise bias correction are plotted in the left and right panels, respectively. Without noise bias correction, the SNR shows a better performance for the strong signal in (a), whereas SNR dominates for to 100 both for theoretical calculation and simulation. With noise bias correction (rightmost column), complex averaging similarly provides an N-fold improvement of the respective SNR′ = I/2σ2 whereas an -fold improvement can be observed for magnitude averaging.

Averaging experimental OCT data

For demonstrating the effect of the different averaging approaches on OCT images, we used a spectral domain (SD) OCT system in our lab [38]. This polarization sensitive SD-OCT system was used for imaging a stationary eye phantom. This system is based on a multiplexed superluminescent diode (λ = 840 nm, Δλ = 100 nm) as a light source, a free-space Michelson-type interferometer, and a polarization-sensitive detection unit including two spectrometers. For the investigations presented here, only the co-polarized detection channel was analyzed; thus the PS-OCT system was reduced to what would be considered a standard SD-OCT with high-resolution imaging capabilities - 3.6 μm axial resolution (assuming a refractive index of 1.35), 83 kHz A-scan rate - similar to state-of-the-art commercial SD-OCT scanners for retinal imaging.

The eye phantom, composed of several layers of transparent nail polish on a glass bead to produce a laminar reflectance pattern [38], was imaged using B-scan and M-scan protocols at different levels of attenuation (from 0 dB to -40 dB) to observe the effects of averaging in low-signal conditions. The B-scan protocol scanned a 1 mm lateral range with 100 repeats to generate cross-sectional images of the phantom (Fig. 8(a-d)). Representative depth profiles are shown in Fig. 8(e).

Fig. 8.

Experimental verification of SNR improvement by the different averaging approaches in a layered retina phantom. (a) Single B-scan image of the phantom (no attenuation). (b) Single B-scan image after attenuating the sample beam by 30 dB. (c) B-scan image after averaging the magnitudes of 100 repeated frames (with 30 dB attenuation). (d) B-scan image after averaging the phasors of 100 repeated frames (with 30 dB attenuation). Note that all B-scan images in (a-d) are displayed with identical dynamic ranges of 40 dB where 0 dB refers to the maximum signal intensity in the frame. (e) Depth profiles of a single A-scan before and after attenuation, 100 magnitude averaged, and 100 complex averaged A-scans at the locations indicated by the dotted lines in (a-d). Due to a beam offset caused by the ND filter, the scattering profile of the unattenuated case has a slightly different structure. Dynamic range as in (a-d). (f) SNR improvement without noise bias correction for three pixels with weak (left), borderline (middle) and strong signal strength (right), respectively. Pixel locations are indicated by orange boxes numbered with 1-3 in panel (d). SNR curves are shown for magnitude and complex averaging of 1-100 repeated M-scan signals for the experimental data () alongside the corresponding theoretical plots (−). (g) SNR′ improvement after noise bias correction for the data sets shown in panel (f). Note that the experimental data () slightly fluctuates around the theoretical profiles (−). SNR′ data fluctuating below SNR′ = 0 is not shown. The axes are scaled as in the respective plots in (f) in order to enable a direct comparison between the two SNR analyses.

The M-scan protocol was used to acquire 1000 repeated A-lines from a fixed position within the phantom for a precise comparison of averaging methods. The M-scan protocol was chosen to minimize phase differences between consecutive A-lines and ensure that averaging represents a best-case real-world scenario. The 1000 A-lines were split into 10 temporally independent bins of 100 consecutive A-lines each. Complex and magnitude averaging for up to 100 was performed on each of the 10 bins and the resulting signal curves were averaged together before calculating SNR and SNR′ in order to yield more stable SNR curves, particularly for lower numbers of averages. Here, the noise variances were computed from the noise pixel data of 100 consecutive A-lines, similar to the simulation above. Three characteristic SNR and SNR′ curves from the magnitude-averaging dominant, borderline, and complex-averaging dominant domains are shown in Fig. 8(f) and (g). Theoretical SNR profiles based on SNR1 as described in Eqs. (25) and (26) demonstrate good agreement with the experimental data in Fig. 8(f). These results show that the theoretical magnitude-averaging dominant domain is reproducible with real world OCT data when the noise floor induced bias is not accounted for. Similarly, the SNR′ plots from the experimental data with noise bias correction follow the expected slopes of and (Fig. 8(g)).

Discussion and conclusion

Signal averaging is one of the most commonly used procedures for OCT image processing. Here, we investigated the SNR performance of magnitude and complex averaging in theory and backed our findings with experimental results from simulations and actual OCT image data. We also studied the effect of the background noise bias on the measured SNR and calculated a noise bias corrected SNR. In the following paragraphs, the main observations are summarized and discussed in order to provide a better understanding of the strengths and weaknesses of the two approaches.

Impact of averaging on noise and signals

Complex averaging reduces the noise variance by , and thus provides a much stronger noise reduction than magnitude averaging, which only reduces the noise by . At the same time, complex averaging also reduces the signal – in contrast to magnitude averaging, which maintains the signal level and only reduces noise. The main cause of the average signal decrease by complex averaging is the reduction of the average noise level. Both signal and noise characteristics have to be considered to understand the impact of averaging on OCT data.

Leaving or removing the noise floor bias for the SNR calculation

The signal-to-noise ratio is a commonly used measure to assess the image quality and to describe the system performance of OCT machines. Usually, the ratio of the average intensity of a signal peak and the standard deviation of the noise intensity is used to calculate the SNR [24,28]. While it is rather clear how to estimate the noise variance (e.g., by considering the noise intensity in a signal-free image area), the definition of the average intensity is not that obvious. As visualized in Fig. 2, taking the entire measured signal intensity from the zero line to the (average) peak intensity includes two portions, namely the actual signal intensity term and the background noise level .

For relatively strong signals, the measured signal is dominated by the intensity term and . In contrast, when the actual signal contribution to the measured signal is small (i.e. ), the noise floor bias prevails. We hence investigated the signal-to-noise performance of magnitude and complex averaging for (a) leaving the noise floor bias as part of the measured signal and (b) correcting for the noise bias. As discussed in detail in the following sections, similar results were observed for both SNR analyses when strong signals were investigated. However, for rather weak signals, the impact of manifested in the observation of dissimilar SNR characteristics.

Noise-afflicted signal-to-noise ratios after averaging strong and weak signals

For sufficiently large signals , the term in Eq. (26) dominates such that complex averaging converges on an N-fold SNR improvement, while magnitude averaging only enhances the SNR -fold. In contrast, for recovering weak signals comparable to or smaller than the noise, magnitude averaging appeared to outperform complex averaging – in particular for small . When the input SNR1 was just slightly larger than unity (i.e., just slightly greater than zero), the right-hand side of Eq. (26) was essentially neutralized such that the SNR improvement was far less than the -fold enhancement yielded by magnitude averaging, especially in the limit of small . Hence, magnitude averaging could be considered more effective for boosting small signals as they may be found in the outer nuclear layer in the retina when the noise floor induced signal bias is not accounted for. However, by taking a look at the SNR in absence of an actual signal (), odd SNR characteristics were observed (see section 2.3.3) which underscored the importance of a noise bias correction.

SNR calculations with noise bias correction

When SNRs are calculated in the classical way described in section 2.3.1, the average noise floor level contributes to the intensity measured as "signal". This noise bias impacts the measured SNR, in particular for weak signals. Akin to noise offset removal in PS-OCT processing [39,40], we performed SNR calculations with noise bias correction in sections 2.3.4 and 3. By using this modified approach, a superior SNR performance was always observed for complex averaging, regardless of the input signal strength. Compared to the noise bias corrected SNR of a single signal prior to averaging, magnitude and complex averaging improved the SNR (after noise bias correction) by a factor of and , respectively. The theoretically predicted performance agreed well with the performance observed in simulated and experimental OCT data (Figs. 7 and 8). This finding suggests that noise bias correction definitely is an important processing step in SNR analyses – especially when it comes to averaging weak signals.

Implications for practical OCT image averaging

State-of-the-art swept source and spectral domain OCT devices deliver complex-valued signals (see Eq. (1)) right out of the box. Hence both magnitude and complex averaging can be easily implemented using the image data. For effective complex averaging, it is imperative that the phases of the signals are accurately aligned before averaging as described in more detail in Appendix B. While this means additional computational steps, it may be well worth the effort when the SNR is to be increased in images of scattering structures such as the retina where OCT has been frequently applied. Retinal OCT images may easily span a dynamic range of 40 dB with hyperscattering structures such as the nerve fiber layer and the pigment epithelium that can serve as landmarks for phasor alignment [18]. Other retinal tissues such as the ganglion cell layer and the outer nuclear layer [41], and also structures in the vitreous, which are usually much less reflective [42], may then be visualized by averaging multiple OCT images together. Complex averaging has also been found promising to detect signals in settings with multiple scattering [20].

Magnitude averaging on the other hand can be implemented at less computational expense and may be the averaging approach of choice for images containing mostly weak signals of interest, which would render phase alignment for complex averaging difficult or even impossible. Also in scenarios where images include a lot of varying motion, e.g. when visualizing weakly scattering structures such as flow in lymphatic vessels [43], albeit theoretically less effective by in terms of SNR improvement, magnitude averaging may likely trump complex averaging.

Future perspectives

Finally, we would like to point out that, while the analysis presented here particularly focused on the SNR improvement for different averaging approaches, it may be interesting to also investigate other image metrics such as their contrast-to-noise characteristics and/or their efficiency in terms of speckle reduction. For this purpose and to explore scenarios with tissue-like scattering, recently described OCT signal models could be particularly interesting candidates [28,44-46]. Additionally, OCT images are often displayed on a logarithmic scale in order to compress signals covering a large dynamic range. The logarithm pulls up weak signals but also the noise floor [21], such that the SNR performance for averaged logarithmic amplitude data will deviate from that for uncompressed OCT data discussed here. An in-depth analysis similar to the one performed for linear data in section 2.2 may provide more insight in the advantages and disadvantages of magnitude- and phasor-based averaging approaches for enhancing OCT image quality.

Table 1.

Overview of variables and abbreviations

Symbol	Explanation	Equation
A(x,t)	Amplitude of phasor representing the OCT signal	(1)
A¯	Mean signal amplitude	(7)
A_noise	Amplitude of complex phasor representing the OCT noise signal	(3)
A_noise¯	Mean noise amplitude	(4)
Δϕ	Phase difference	(39,42)
I	Signal intensity, I=A^2
I¯	Mean signal intensity	(13)
I_x¯^′	I_x¯ with noise bias correction, (x = void,MAG,CPX)	(31–33)
I_noise	Noise intensity, I_noise=Anoise2
I_noise¯	Mean noise intensity	(11)
i_noise	Imaginary part of complex phasor representing the OCT noise signal	(2)
N	Number of averaged signals
p_A(A_noise)	PDF describing amplitudes of random phasors (Rayleigh distribution)	(3)
p_Δϕ(Δϕ,A)	PDF describing distribution of phase differences	(43)
p_I(I_noise)	PDF describing intensities of random phasors	(9)
p_I(I,I_noise)	PDF describing intensities of signals in presence of noise	(10)
p_A(A,A_noise)	PDF describing amplitudes of signal phasors in presence of noise (Rice distribution)	(6)
p_ri(r_noise,i_noise)	Binormal PDF describing random phasors	(2)
P(N)	Penalty factor after averaging N signals	(37)
PDF	Probability density function	(2,3,6,9,10)
ϕ(x,t)	Phase of phasor representing the OCT signal	(1)
r_noise	Real part of complex phasor representing the OCT noise signal	(2)
S_OCT(x,t)	Complex-valued OCT signal	(1)
SOCT(moco)	OCT signal after motion correction (moco)	(40,41,45)
σ^2	Variance of real/imaginary part of p_ri	(1)
σA2	Variance of signal amplitudes for p_A(A,A_noise)	(8)
σI2	Variance of signal intensity for p_I(I,I_noise)	(14)
σI_noise2	Variance of noise intensity for p_I(I_noise)	(12)
σnoise2	Variance of noise amplitudes for p_A(A_noise)	(5)
SNR_1	SNR of non-averaged signal	(24)
SNR_1,borderline	SNR1 of non-averaged signal for SNR_MAG=SNR_CPX	(27)
SNR_MAG	SNR after magnitude averaging	(25)
SNR_CPX	SNR after complex phasor averaging	(26)
SNRx′	SNR_x with noise bias correction, (x=1,MAG,CPX)	(34–36)
⟨⋅⟩_MAG	Quantity ⋅ after magnitude averaging	(16–19)
⟨⋅⟩_CPX	Quantity ⋅ after complex averaging	(20–23)