Optical coherence tomography (OCT) is an interferometric imaging modality that allows the detection of a complex signal. Its high acquisition rate, high axial resolution of several microns, and ability to create 3D images make it the ideal instrument for retinal imaging. OCT can typically scan at rates of around 100,000 A-lines per second and has a sensitivity of over 90 dB.

However, single OCT B-scans are often degraded by noise and additional processing is required to make the images useful for clinical interpretation. In addition, the OCT scattering signal acquired is not quantitative. For example in posterior eye imaging, the vitreous body appears with a weak random signal, where one expects a much lower signal because its purpose is to be transparent.

Recently there has been a growing interest in vitreo-macular pathologies [1]. Therefore we anticipate that accurate quantification of low-scattering regions will become important [2,3]. Although intensity averaging is widely applied to enhance the image contrast [4], this type of processing can result in erroneously bright values at low-intensity regions such as the vitreous. To compensate for these noise effects at low intensities, the threshold is often applied before averaging the signals. However, this causes an arbitrary loss of information.

Here we present a signal amplitude estimation method based on a statistical model of OCT signal and noise. This method is a maximum a posteriori (MAP) estimator of OCT signal intensity and is less biased at low-intensity regions than intensity averaging. Hence it provides a quantitative estimate of the true underlying backscatter intensity even at the low-intensity regions. This enables the generation of quantitative, high-contrast images and the estimation of intensities below the noise floor.

Intensity averaging is a widely used processing method to enhance the image contrast of OCT. While the averaging is simple and easy to compute, it is unbiased only if a signal under noise is symmetrically distributed about the true signal. This assumption would be applicable to each of the real or imaginary parts of complex OCT signals, as they are affected by additive white Gaussian noise (AWGN). However, due to the phase instability of OCT signals, averaging each of the real and imaginary parts separately cannot improve the image contrast in most practical cases. On the other hand, the amplitude or intensity of OCT signals is not affected by the phase instability. Hence its contrast can be improved by averaging. However, the amplitude distribution of OCT signals under complex AWGN is Rician. Similarly, the intensity (amplitude squared) of OCT is related to the noncentral chi-squared distribution. Neither are symmetric and both are positively skewed, especially if the signal amplitude is low. Hence, averaging of amplitude and intensity yields positively biased results, especially for low intensities.

A more theoretically sound approach would be to take into account the statistics of the OCT signal and noise. Each of the real and imaginary parts of an OCT signal are affected by AWGN. Therefore, the amplitude of the signal is more accurately modeled by the Rice distribution as it has been done for magnetic resonance imaging [5]. For instance, the probability density function of observed OCT signal amplitude, $a_n$, under a “true” signal amplitude of $\nu$ is given by [6]

Assume $a_n$ represents the signal amplitude of the $n$th B-scan when multiple, $N$, OCT B-scans were acquired at the same location of a sample. For measured data, $\{a_1,\cdots ,a_n,\cdots ,a_N\}$, the combined likelihood of the true OCT signal is given by the product of the likelihoods of each datum as $l(\nu )=\prod _{n=1}^{N}f(\nu ;a_n,{\sigma }^2)$. From the combined likelihood, the MAP estimation of the true signal under the measured data is obtained as

From the combined likelihood, one may also obtain the 95% credible interval based on the likelihood ratio statistic [8]. The reciprocal of this interval is used as a measure of estimation reliability. Quantifying the reliability of an estimation is important to avoid false interpretations of data. The likelihood ratio statistic is given by

In order to minimize the processing time, Eq. (2) was implemented as $\widehat{\nu }=\arg {\max }_{\nu }\underset{n=1}{\overset{N}{\mathrm{\Sigma }}}\log f(\nu ;a_n,{\sigma }^2)$ in our particular implementation. Each likelihood was computed prior to the MAP process for combinations of 500 true amplitude levels and 200 noise levels, and stored in a lookup table. The amplitude was equally spaced in log scale, while the noise levels were equally spaced in linear scale. While the support of a Rice distribution is from zero to $\infty$, we computed the likelihood only from $a_n=0$ (numerically ${10}^{-6}$) to the maximum amplitude in the measured OCT data volume. This is equivalent to setting the prior, $p(\nu )$, to one at these amplitudes. A specific likelihood is selected from the lookup table based on a measured signal intensity and a noise variance. The noise was obtained at each depth from an OCT B-scan image taken while obstructing the sample arm. The noise variance at each depth is calculated by summing the variances of each of the real and imaginary parts. In our implementation, the MAP estimation procedure took 140.3 s, using four parallel threads in Matlab 2014b on a Intel Core-i7-5930 K PC with 32 GB of memory for a data consisting of 2048 A-lines, each with 512 pixels, and the number of frames $N=24$.

To test the MAP estimation algorithm, a normal macula of the left eye of a 29-year-old man was scanned by using a custom-built 1.0 μm-spectral-domain OCT at 91,900 A-lines/s, with an axial resolution of 4.8 μm in tissue and with a system sensitivity of around 93 dB [9]. 2048 A-lines were taken per B-scan with a transverse range of 2.6 mm. Twenty-four B-scans obtained at the same location were co-registered by an A-line-wise, cross-correlation-based method. The registration shifted each A-line with respect to a central reference frame (twelfth B-scan). The co-registration was visually inspected in the form of movies and were considered satisfactory. These 24 B-scans were combined to form composite images using a MAP estimator as well as standard arithmetic averaging of the signal intensity. Composite images were also created for 1, 2, 4, 6, 8, 10, 13, 16, and 20 frames, respectively, to analyze estimation performance.

We also examined a pathologic case of a 62-year-old man diagnosed with polypoidal choroidal vasculopathy (PCV). This patient was scanned at an earlier date using a similar 1.0 μm-spectral-domain OCT system but with a scan rate of 47,000 A-lines/s. The axial resolution and axial pixel separation was the same as above, but the transverse scan range was 6.0 mm. Each scan had 1500 A-lines. Twenty-four frames at the same location were used for MAP estimation and averaging.

The data acquisition protocol adhered to the tenets of the Declaration of Helsinki and was approved by the institutional review board of Tokyo Medical University and the University of Tsukuba. Informed consent was obtained from the subjects.

Figure 1 shows the macular cross sections of the healthy subject. The MAP image [Fig. 1(a)] shows a darker appearance at the vitreous and deep choroidal (out-of-penetration) regions than the intensity-averaged image [Fig. 1(b)] or a single B-scan [Fig. 1(c)]. The intensity average at the low-signal region, where the probe power is lower than the noise floor, approaches the mean noise energy. On the other hand, MAP estimation approaches the true signal amplitude even at the low-signal region owing to the statistics of multiple pixels. This effect can be seen by comparing the vitreous intensity to the outer nuclear layer (ONL) intensity (Fig. 1). In the MAP image, the ONL (lower box) appears brighter than the vitreous (upper box), whereas in the intensity-averaged image, their intensities appear similar.

 

Fig. 1. (a) Macula of a healthy subject obtained by MAP estimation using 24 B-scans, (b) intensity averaging, (c) and a single representative B-scan. The SIR of the vitreous (delineated by the upper box) to ONL (lower box) was measured to be 2.09 dB for the MAP image and 0.36 dB for the intensity-averaged image. The lower black box indicates the region of RPE utilized for the analysis of Fig. 4.

Download Full Size | PDF

The representative A-line plots shown in Fig. 2(a), along the vertical dashed lines in Figs. 1(a) and 1(b), show around 6.9 dB suppression of the noise floor in the vitreous region with MAP estimation (blue line) compared with intensity averaging (red line). Here 0 dB indicates the ninety-ninth percentile intensity of the whole image. Figure 2(b) shows line profiles at the horizontal dashed lines depicted in Figs. 1(a) and 1(b). It demonstrates that MAP estimation also provides better contrast of the deep choroid region. Hence MAP estimation also provides better contrast of structure as observed around the ONL and the choroid.

 

Fig. 2. Representative line profiles of MAP estimation and intensity averaging along the vertical and horizontal dashed lines in Figs. 1(a) and 1(b). (a) Plot shows the vertical profile and (b) plot shows the horizontal profile.

Download Full Size | PDF

The histograms in Fig. 3 illustrate the advantage of using MAP filtering over averaging. The MAP-filtered image [Fig. 3(a)] maintains the same dynamic range of around 60 dB and the same general distribution of amplitudes with the unfiltered image [Fig. 3(c)], while also estimating very low (near zero) intensity values. This result is rational, as the vitreous is expected to have little to no scattering. On the other hand, the intensity-averaged image [Fig. 3(b)] appears to have a lower dynamic range as values below $-38\mathrm{dB}$ appear to be shifted upward. As a result, there appears to be a lower dynamic range and less contrast in Fig. 1(b). The cluster of pixels around $-130\mathrm{dB}$ in Fig. 3(a) represents estimations of zero intensity. It corresponds to the lowest possible estimation in our numerical implementation of the MAP estimator.

 

Fig. 3. Signal intensity histograms of images shown in Fig. 1. (a) is of the MAP estimation image, (b) is of intensity-averaged image, and (c) is of the single image. The signal intensity is in decibel scale, and the count is in logarithmic scale. Here 0 dB intensity indicates the ninety-ninth percentile intensity of the whole image. (d) MAP estimation reliability corresponding to Fig. 1(a). The decibel scale of (d) is defined as $-10{\log }_{10}$ (upper bound of the 95% credible interval in amplitude scale − its lower bound), which is not the same as the decibel scale of OCT images.

Download Full Size | PDF

Figure 3(d) shows the reliability of the MAP estimation, based on Eq. (3). Since the region with higher signal intensity or, technically, higher signal to noise ratio (SNR) shows higher reliability, the reliability distribution has a similar appearance to the MAP image [Fig. 1(a)]. Some pixels around the retinal pigment epithelium (RPE) complex appear to have low reliability. A possible reason is that even a small registration error causes high signal fluctuation among the frames at this region that has high intensity and fine structure.

As the RPE and the vitreous are expected to have very high and low scattering, respectively, the signal intensity ratio (SIR) between the RPE and the vitreous can be a measure of the effective SNR. The SIR of MAP ($\mathrm{\Delta }$) and averaged ($\square$) images are plotted in Fig. 4 over the number of frames used for the estimation. The SIRs were computed from the regions indicated in Fig. 1. MAP estimation shows higher SIR than averaging for any number of frames. It can be seen that additional frames beyond eight yield little improvement with MAP, so four to eight frames can be a good compromise between image quality and measurement speed. Furthermore, one-frame MAP estimation also shows an SIR improvement over a single-frame image, i.e., one-frame averaging.

 

Fig. 4. Plot of SIR between the RPE and the vitreous against number of frames for MAP ($\mathrm{\Delta }$) and intensity averaging ($\square$) for the healthy subject.

Download Full Size | PDF

In the pathologic (PCV) case, the MAP image provides high contrast of the hyaloid membrane (indicated by the arrow) as shown in Fig. 5. The SIR between the hyaloid membrane and the vitreous is 14.92 dB for the MAP image [Fig. 5(a)], while it is 7.11 dB for the intensity-averaged image [Fig. 5(b)]. This comparison was made using the same region of interest (that showed high contrast) for both images. The MAP estimation image also shows higher contrast in the subretinal pigment epithelial regions (circled regions). Since vitreous imaging has been attracting the interest of ophthalmologists recently [3], this property of MAP filtering is expected to have high clinical utility.

 

Fig. 5. (a) Macular cross section images obtained from a PCV case with MAP estimation using 24 B-scans, (b) intensity averaging, and (c) a single representative image. The arrows indicate the hyaloid membrane. The circled regions indicate areas of suspected subretinal fluid.

Download Full Size | PDF

MAP estimation provides more quantitative estimations of OCT signal amplitude and intensity with less bias than standard averaging, especially at the low-intensity regions. For example, averaging evidently gives positively biased estimation at the low-signal regions as shown in Figs. 1(b), 3(b), and 5(b). This is because averaging is only unbiased if the noise distribution is symmetric and centered at the true signal value, while the intensity-noise distribution becomes highly skewed at the low-signal regions. On the other hand, the MAP estimator takes into account the theoretical a priori knowledge of the signal distribution and can provide low-bias estimation, even for low intensities and low SNRs. Although histogram equalization can enhance contrast, unlike MAP estimation, it cannot improve statistical or quantitative analysis.

The MAP estimator we present can estimate signal values that are smaller than the shot-noise-limited noise floor of a single-image frame, although it cannot overcome the shot-noise-limit sensitivity defined by the total exposure time of all the frames used for MAP estimation. On the other hand, intensity averaging cannot even overcome the single-frame shot-noise limit sensitivity. The advantage of the MAP estimator is that it uses all the information contained in the probability density function, whereas averaging uses only the first-order statistics. Similar to the MAP estimator, complex averaging of frames also can overcome the single-frame, shot-noise-limit sensitivity, but it requires phase consistency among all frames, which is not likely to be obtained for in vivo measurement.

In conclusion, it was demonstrated that the MAP estimator provides less biased estimation of OCT signal amplitude and intensity. It also provides higher contrast images at the vitreous and choroid than standard signal intensity averaging and can enhance the clinical utility of OCT. The quantitative estimation of MAP can be used as preprocessing for further analysis of OCT. In particular, OCT intensity images obtained by MAP estimation can enhance the ability of other quantitative signal evaluation methods, such as attenuation coefficient imaging [10]. This will finally enable more accurate diagnosis of pathologies.