Open Access
Add to BibSonomy
Abstract
The attenuation coefficient provides a quantitative parameter for tissue characterization and can be calculated from optical coherence tomography (OCT) data, but accurate determination requires compensation for the confocal function. We present extensive measurement series for extraction of the focal plane and the apparent Rayleigh length from the ratios of OCT images acquired with different focus depths and compare these results with two alternative approaches. By acquiring OCT images for a range of different focus depths the optimal focus plane difference is determined for intralipid and titanium oxide (TiO2) phantoms with different scatterer concentrations, which allows for calculation of the attenuation coefficient corrected for the confocal function. The attenuation coefficient is determined for homogeneous intralipid and TiO2 samples over a wide range of concentrations. We further demonstrate very good reproducibility of the determined attenuation coefficient of layers with identical scatter concentrations in a multi-layered phantom. Finally, this method is applied to in vivo retinal data.
© 2021 Optical Society of America under the terms of the OSA Open Access Publishing Agreement
1. Introduction
Optical coherence tomography (OCT) is primarily used to image the morphology of tissue by visualizing different structures in arbitrary units to diagnose pathologies in medical specialties such as ophthalmology, cardiology, dermatology, and urology. However, the state of a tissue can also be described in absolute units by its inherent optical properties, such as the attenuation coefficient, providing a quantitative tissue parameter. Vermeer et al. [1] introduced and experimentally demonstrated a method that allows estimation of the attenuation coefficient with a pixel level resolution from OCT depth profiles. If the focal plane is located close to or within the sample the confocal function for the OCT measurement geometry [2] has to be taken into account to precisely determine the attenuation coefficient. Smith et al. [3] suggested a method for automated, depth-resolved estimation of the attenuation coefficient with the focal plane located within the sample; however, this method requires a priori knowledge of the confocal function. Determining the confocal function constitutes estimating the position of focus and the apparent Rayleigh length.
The confocal function also depends on the sample that is examined. In many applications, especially ophthalmology, the confocal function parameters, i.e., focus location and apparent Rayleigh length vary. It would be beneficial to extract these parameters from data acquired with the matching setup and sample or subject. Ghafaryasl et al. [4] have suggested a method to estimate the confocal function and the attenuation coefficient directly from a single scan. However, their method is currently limited to homogeneous samples. Dwork et al. [5,6] have demonstrated that the ratio of two B-scans can be used to estimate the focus location and the apparent Rayleigh length and that this can be applied to inhomogeneous samples. Stefan et al. [7] extended this technique by accounting for the curvature and tilt of the focal plane.
In this paper, we analyze the applicability and limitations of extracting the focus depth and/or the apparent Rayleigh length from OCT images using three different methods (Section 3.7 focus series, Section 3.8 ratio fit, and Section 3.10 direct fit). We find that the most robust and versatile method to extract the confocal function parameters in the sample is to use the ratio fit of two B-scans acquired with different focal depths, which is based on the work of Dwork et al. [5,6]. We investigate this method by imaging homogeneous intralipid samples and titanium oxide (TiO2) phantoms of different scatterer concentrations with an OCT system through a series of different focal depths from the sample surface. The data series are analyzed for the optimal focal depth offset for confocal function determination. The obtained focus depths and Rayleigh lengths are used to correct for the confocal function and subsequently calculate the attenuation coefficients. For comparison, we also test a technique where we directly fit the confocal function and attenuation coefficient in a single A-scan. Finally, the ratio fit method is demonstrated for a layered phantom and in vivo retinal data by determining the curved focal planes and calculating attenuation coefficient maps.
2. Theory
In this section, we summarize the main aspects of the Gaussian beam model needed to understand how raw OCT data can be corrected for the confocal function. Single backscattered light detected within an OCT A-scan can be described by [1]
where z is the depth within the A-scan, r(z) is the sensitivity roll-off, h(z) is the confocal function, α is the backscattering fraction, and µt(z) is the attenuation coefficient. The attenuation coefficient µt is the sum of the scattering coefficient µs and the absorption coefficient µa, µt = µs + µa. For a Gaussian beam, the confocal function h(z) can be described by [2] with the focus position zf and the apparent Rayleigh length zR. For a scattering sample with refractive index n the apparent Rayleigh length can be calculated as [2] with the central wavelength $\lambda $ and the $1/{e^2}$ radius of the lateral intensity distribution of the focused Gaussian beam ${\omega _0}$ which can be estimated as with the focal length f and the beam diameter at the pupil d.To determine zf and zR one can exploit the common characteristics of two scans of the same sample location [5]. If two A-scans of the same structure are acquired with different focus depth settings, all the terms in Eq. (1) are the same, except for h(z). Therefore, the ratio of the intensity of the numerator A-scan I1(z) and the denominator A-scan I2(z) can be described by the ratio of the two confocal functions as
3. Methods
3.1 OCT system
OCT measurements were performed with a high-resolution spectral domain OCT device based on the SPECTRALIS platform (Heidelberg Engineering GmbH, Heidelberg, Germany) with a central wavelength of 853 nm and a spectrometer bandwidth of 137 nm. Using the active tracking capability of the device many B-scans of the same location can be averaged together. The focus of the device can be changed by moving the objective relative to the rest of the device using a focusing knob. The focus setting in diopters (1/m) is displayed in the graphical user interface of the acquisition software and saved with each data set.
3.2 Data acquisition and B-scan processing
For each sample, a focus series of B-Scans was acquired at various depths of the focal plane referenced to the sample surface. All the data processing was performed on the raw spectral data stored by the device using custom-written software in Python 3.7 (Python Software Foundation) and MATLAB 2020b (The MathWorks, Inc.).
First, to reduce fixed pattern noise [8], the raw spectra were divided by the mean spectrum of all A-scans and centered around zero by subtracting one. Second, the nonlinear relation between wavenumber k and wavelength λ was compensated for by cubic interpolation to linearize the signal in wavenumber [9]. Third, to suppress sidelobes caused by the Fourier transformation of a finite interval, the spectrum was windowed using the familiar Hamming window (a0 = 0.53836 and a1 = 0.46164). Fourth, numerical dispersion compensation was performed [10,11]. Fifth, the Fourier transformation of the signal was squared. Sixth, the depth dependent background was subtracted, which was measured by acquiring a B-scan without a sample. Seventh, to account for sensitivity changes with depth, the obtained B-scans were corrected for roll-off which was measured by moving a mirror through the sampling depth of the OCT scan [12]. Finally, to reduce speckle noise and average over imaging artifacts, 120 B-scans were averaged (768 A-scans/B-scan) to produce one B-scan in which each A-scan was averaged 120 times. These A-scans were used for further processing.
3.3 Model eye
For all phantom measurements, a model eye was used featuring an 18 mm focal length planoconvex lens (01LPX015, Melles Griot, anti-reflection coating on the convex side) and a water filled chamber that can hold different types of samples.
3.4 Intralipid samples
To assess the effect of different scattering strengths, a serial dilution ranging from 6.22% volume to 0.062% volume intralipid (Table 1) was produced by diluting an intralipid 20% emulsion (I141-100ML, Sigma-Aldrich, Saint Louis, MO, USA) with distilled water. Cuvettes (Macro cell 100-5-20, Hellma GmbH & Co. KG, Müllheim, Germany) were filled with the dilutions and imaged inside the model eye with a scan angle of 9 degrees. The nominal scattering coefficients for the OCT wavelength of the different dilutions were calculated by linear interpolation of the scattering coefficient for intralipid 20% described by Michels et al. [13]. The refractive index of the samples was assumed to be equal to the refractive index of water and the absorption coefficient of the intralipid samples was assumed to be negligible (µs ≫ µa, and µt = µs) [13].
Table 1. Intralipid samples
3.5 Titanium oxide phantoms
Silicone-based phantoms were fabricated similar to the method outlined by de Bruin et al. [14], where the scattering coefficient is controlled by varying the % weight of TiO2 in the final elastomer. A TiO2 (634662, Sigma Aldrich, Saint Louis, MO, USA) and PDMS (481955, Sigma Aldrich) stock dispersion was thoroughly incorporated in an elastomer base (SYLGARD 184 Silicone Elastomer, Dow, Inc., Midland, MI, USA) then cast into the final phantom geometry and cured. Seven bulk phantoms and one layered phantom were fabricated with TiO2% weight concentrations ranging from 0.011 to 0.730% weight, outlined in Table 2. The layered phantom (Fig. 1) was created by casting three thin films of the substrate with different scatterer concentrations on top of each other and afterwards stacking five of these layers. This procedure was chosen to ensure that films with the same nominal scattering coefficient have the exact same physical properties. In our experience, even multiple pours of the same formulation can be unreliable in producing the same physical properties. Although the samples are flat, the apparent curvature of the sample, that can be observed in Fig. 1, is a typical OCT imaging artifact originating from the optical path length differences of the beam depending on the scan angle. The curvature depends on the position of the scan pupil relative to the focusing lens and on the lens itself. All phantoms were imaged inside the model eye with a scan angle of 9 degrees. The refractive index of the phantoms was assumed to be 1.42 [14], and the absorption coefficient was neglected, µt = µs, as no absorbers were added.
Fig. 1. Upper part of a B-scan of the layered TiO2 phantom. The scale bars are 200 µm in each dimension.
Table 2. TiO2 Phantoms
3.6 In vivo data
To demonstrate the applicability of our analysis to in vivo data, scans of the retinas of two subjects were obtained at the locations shown in Fig. 12(C) and (F) using the tracking capability of the device to ensure acquisition at the same location for all focus settings. To avoid accommodation during the measurement an external fixation target was used. No drugs were applied to dilate the pupil or paralyze accommodation. This paper represents a feasibility study with the focus on the evaluation and correction algorithms and not a systematic analysis of clinical data. A written declaration of informed consent was obtained from both subjects.
3.7 Expected focus depth and Rayleigh length estimation using the focus series
To estimate the theoretical expected focus depth and the actual apparent Rayleigh length, versus the theoretical apparent Rayleigh length calculated using Eq. (3), a focus series fit was performed. To suppress high frequency noise, a two-dimensional Gaussian filter with a standard deviation of 3 pixels was applied to the B-scans and to further reduce noise a mean A-scan was calculated for each B-scan by averaging the 61 central A-scans (location shown in Fig. 2(A)). Next, the intensity just below the surface of the sample was determined for every focus position (Fig. 2(B), red dashed line). This location was chosen to avoid the surface reflection and minimize multiple scattering contributions expected to increase with depth. Each determined intensity value corresponds to a focus setting in diopters (1/m) which was saved by the system together with the B-scan data. The focal length f can be calculated from the diopter value D in 1/m and the focal length in air of the lens used for imaging flens in meters according to
with the refractive index of the medium n. Because the medium above the sample has a different refractive index than the sample itself, the theoretical focus distance from the sample surface Δzfsample was determined to fit for the Rayleigh length. Therefore, the focus setting at which the sample surface would be perfectly in focus Dsample was estimated by fitting A([(D - Dsample)/R]2+1)−1 to the intensity values, with Dsample, R and A as fitting parameters. Using Eq. (6) fsample was calculated from Dsample. The theoretical focus distance from the sample surface was calculated as Δzfsample= f - fsample, where f and fsample were calculated with the respective refractive index depending on if the focus was above the sample (D > Dsample) or within the sample (D < Dsample). To estimate the apparent Rayleigh length the Gaussian beam confocal function (Eq. (2)) multiplied by a constant A was fitted to the intensity values over Δzfsample as described by Faber et al. [15] with A, zf and zR as fitting parameters (Fig. 2(C)).Fig. 2. (A) Examples of B-scans from intralipid sample I3 (nominal µs of 2 mm-1) acquired with different focus depths (127 µm, 398 µm, and 987 µm below surface) showing the location of the mean A-scans plotted in (B). The scale bars are 200 µm in each dimension. (B) Mean A-scans calculated from the 61 central A-scans for different focus depths, the black arrow indicates increasing focus depth. The sample depth used to measure the intensities for the confocal function fit is marked with a red, dashed line. (C) Measured confocal PSF (blue, asterisk) and best fit (red, solid).
3.8 Determination of confocal function using the ratio fit
To extract the confocal function parameters, Eq. (5) describing the intensity ratio of two A-scans (numerator I1(z) and denominator I2(z)) acquired with two different focus depths was used in the decibel form as
Fig. 3. Example of a ratio fit (intralipid sample I3, nominal µs of 2 mm-1). The intensity ratio data of two mean A-scans in dB acquired with a focus depth location of 127 µm and 398 µm below the surface is shown in green, the resulting ratio fit with a fitted focus depth of 127 µm (standard error 0.4 µm) below the surface and a fitted apparent Rayleigh length of 290 µm (standard error 0.7 µm) is depicted with the blue dashed line.
Fig. 4. Matrices of fitted focus depth (A) and Rayleigh length (B) for all possible focus setting ratios of intralipid sample I3 (nominal µs of 2 mm-1). The values in the Matrices represent the ratio fit result. On the horizontal and vertical axes, the focus depth number of the corresponding focus setting of the denominator and numerator used for the ratio are given, respectively. Additionally, the vertical axes are labeled with the expected focus depths of the corresponding focus numbers with respect to the sample surface in square brackets. The color coding represents the standard error of the corresponding ratio fit and does not correspond to the value of the numbers in the matrix.
3.9 Calculation of attenuation coefficient
To correct for the axial point spread function (PSF), the A-scans were divided by the estimated confocal function described by the mean focus depth and Rayleigh length obtained for each respective focus setting. Subsequently, the depth resolved attenuation coefficient profile µt(i) was calculated using the method as described by Girard et al [16] and Vermeer et al. [1] as
For the homogeneous phantoms, a reference value for the attenuation coefficient was estimated using the commonly applied slope fitting method [15,17] where a single scattering model in the form I(z) ∝ exp(-2µz) is fitted to the A-scans.
3.10 Direct fit for confocal function and attenuation coefficient
Additionally, for the homogeneous phantoms, the method introduced by Ghafaryasl et al. [4] was tested. It has the advantage that focal position and Rayleigh length do not have to be determined beforehand and only a single scan is needed, removing the necessity to acquire two scans from the exact same position. The single scattering model
4. Results
Sections 4.1 to 4.3 describe the results obtained from the homogeneous samples, which were analyzed by calculating a mean A-scan from the 61 central A-scans of an averaged B-scan. This simplified our analysis because within these A-scans we assumed a uniform focus depth. The results for the layered phantom (Section 4.4) and the in vivo data (Section 4.5) were obtained by analyzing the whole averaged B-scans which were filtered as described in Section 4.4.
4.1 Ratio fit data processing
As an example, Fig. 4 illustrates the ratio fit result matrices for the intralipid sample I3 (nominal µs of 2 mm-1), showing the results for all possible combinations of different focus depths fitted in Eq. (7) (color-coded for the standard errors of each fit). The fit matrices for the other intralipid concentrations are provided in the Supplement 1.
Ratio combinations with both focus depths (numerator and denominator) above the sample show a very large standard error and mostly fitting the model (Eq. (7)) fails which is reflected in fitted apparent Rayleigh lengths of zero. The lowest fit errors are obtained for ratios where at least one focus location is located well inside the sample with a focus difference between the two tested A-scans of more than 100 µm. For focus locations very deep in the sample an increasing fitting error can be observed.
Calculating the error-weighted mean fit result for each horizontal line of the ratio fit result matrices (Fig. 4) results in the values given in Fig. 5. The mean fitted focus depths agree well with the expected focus depths obtained from the focus settings of the device, and the fitted mean Rayleigh lengths are located around the expected Rayleigh length calculated based on the beam parameters (Eq. (3)).
Fig. 5. Weighted-mean result of the ratio fit method, focus depth zf (blue, circles) and Rayleigh length zR (red, asterisks) for the different focus settings in the intralipid sample I3 (nominal µs 2 mm-1). The error bars indicate the error-weighted standard deviations. The theoretical focus depth (blue) based on the diopter settings of the system and the theoretical Rayleigh length (320 µm, red) calculated using Eq. (3) are indicated with dashed lines.
Using the confocal function parameters recorded in Fig. 5, the respective mean A-scans shown in Fig. 6(A) were corrected for the confocal function. The uncorrected mean A-scans (Fig. 6(A)) acquired with different focus depths clearly show the effect of the confocal function, which is eliminated in the corrected mean A-scans (Fig. 6(B)).
Fig. 6. Mean A-scans before confocal correction (A) and after (B) for the intralipid sample I3 (nominal µs 2 mm-1). The legend indicates the theoretical focus offsets from the sample surface in µm.
4.2 Fitted Rayleigh length
We determined the Rayleigh length using three different methods (as described in Section 3.7 focus series fit, Section 3.8 ratio fit, and Section 3.10 direct fit) for the homogeneous intralipid and TiO2 phantoms with several different scatterer concentrations. For each method and scatterer concentration, the average Rayleigh length was determined (Fig. 7) by calculating the weighted mean of the Rayleigh lengths calculated for focus depths from 0 µm to 1200 µm, weighted with the corresponding standard error. For nominal scattering coefficients greater than 2 mm−1 the standard deviation of the direct fit exceeded 50 µm, which we considered as a failed fit and values are therefore not shown. Figure 7 shows that the fitted apparent Rayleigh length increases with increasing scatterer concentration for the ratio fit method. The same trend can be seen for the direct fit results. However, for the focus series fit the trend of increasing apparent Rayleigh length with increasing scatterer concentration is not as pronounced. These values were calculated from the intensities just below the surface minimizing the contribution of multiple scattering. For the intralipid samples (Fig. 7(A)) the data series are more consistent and show less fluctuation than for the TiO2 phantoms (Fig. 7(B)). Here, it is worth mentioning that the averaged OCT B-scans of the liquid intralipid samples (Fig. 2(A)) appear much more uniform and show less speckle noise compared to the B-scans of the solid TiO2 phantoms (Fig. 1) because Brownian motion averages out the speckle.
Fig. 7. Apparent Rayleigh length calculated with different fitting procedures for different intralipid (A) and TiO2 (B) concentrations. The horizontal axis indicates the nominal scattering coefficients of the phantoms, the corresponding scatterer concentration is referenced in Table 1 and Table 2. For scattering coefficients greater 2 mm-1 the direct fits failed and therefore are not shown. The theoretical Rayleigh lengths for the intralipid phantom (320 µm) and TiO2 phantom (343 µm) calculated with Eq. (3) are indicated with a dashed line.
4.3 Attenuation coefficient
To demonstrate the validity of the ratio fit method for determining the focus position and Rayleigh length in the sample of a particular mean A-scan, the attenuation coefficients calculated from the mean A-scans acquired with different focus settings are depicted in Fig. 8 (corrected for the fitted Rayleigh lengths and focus depths (Fig. 5)). Both, the depth resolved method and the slope fit method result in a calculated attenuation coefficient that is more or less constant as a function of the focus depth. If the confocal function is not corrected for, a reliable estimation of the attenuation coefficient is not possible which can also be seen in Fig. 8, where the values for the uncorrected mean A-scans vary considerably.
Fig. 8. Calculated attenuation coefficient from mean A-scans with different focus settings for intralipid sample I3 (nominal µs 2 mm-1). Results for the uncorrected mean A-scans are indicated by the dots and crosses for the slope fit and depth resolved method, respectively. Results corrected for the confocal function are indicated by the circles and asterisks for the slope fit and depth resolved method, respectively.
The weighted mean attenuation coefficients for focus locations from 0 µm to 1200 µm below the surface, weighted with the standard error of the corresponding fit, are shown in Fig. 9 for intralipid (A) and TiO2 (B) samples with different scatterer concentrations. The direct fit results were obtained directly from uncorrected mean A-scans, whereas slope and depth resolved values were calculated from the confocal function corrected mean A-scans using the parameters from the ratio fit. For scattering coefficients lower than 6 mm-1 the calculated attenuation coefficients agree well with the theory independent of the algorithm used. For strong scattering, the attenuation coefficients are consistently underestimated by both the depth resolved and the slope method (the direct fit failed in this range).
Fig. 9. Attenuation coefficients calculated with the different procedures for intralipid samples (A) and TiO2 phantoms (B). The theoretical scattering coefficients obtained from [13] and [14] by linear extrapolation are indicated with a dashed line. The direct fit results (yellow triangles) were obtained directly from uncorrected mean A-scans. For nominal scattering coefficients greater than 2 mm−1 the standard deviation of the direct fit exceeded 50 µm, which we considered as a failed fit and values are therefore not shown. Slope (red circles) and depth resolved (blue asterisks) values were calculated from the confocal function corrected mean A-scans using the parameters from the ratio fit.
4.4 Layered phantom
In general, the focal plane cannot always be assumed to be flat and is rather curved and/or tilted within the sample [7]. For the layered TiO2 phantom, the focal planes of the different focus settings were determined by applying the ratio fit method to every A-scan extracting any curvature or tilt of the focal plane. Prior to the confocal function determination, a two-dimensional Gaussian filter with a standard deviation of 5 pixels and a moving mean filter in the x direction with a window size of 11 pixels was applied to suppress high frequency noise. The parameters for the two filters were chosen as a compromise between noise reduction and conserving spatial resolution. For the ratio fit we were interested in the confocal function which varies slowly across the B-scan, therefore preserving hard edges was not a priority and some blurring was accepted. To reduce noise in the raw two-dimensional map of the extracted focus plane and Rayleigh lengths, we fitted a second order polynomial to the extracted raw zf and zR values and used the smoothed values for the further analysis. The resulting smooth curved focal planes can be seen in Fig. 10(A). The resulting confocal functions were used to correct the unfiltered B-scans, from which attenuation coefficient maps were calculated using the depth resolved method. The attenuation coefficient map obtained for the B-scan with the focal plane position indicated by the white arrow in Fig. 10(A) is shown in Fig. 10(B). Layers with the same scatterer concentration exhibit similar absolute attenuation coefficient values. The mean attenuation coefficients were calculated for each layer marked in Fig. 10(C) and every focal plane location. The resulting values are given in Fig. 10(D) demonstrating the reproducibility of the attenuation coefficient calculation and the effectiveness of the confocal function correction. Independent of the depth within the phantom, layers produced from the same pour (A1 = A2, B1 = B2, C1 = C2) exhibit similar attenuation coefficients. The obtained µt are in the same range as the nominal µs (Table 2) indicated with the dashed lines but have an offset.
Fig. 10. (A) Upper part of a B-scan with the focal planes in the layered TiO2 phantom calculated using the ratio fit method. The different line stiles indicate the focus planes that are mainly above the sample surface (dotted), less than 500 µm within the sample (solid), and further than 500 µm below the sample surface (dashed). (B) Attenuation coefficient map of the confocal function corrected B-scan acquired with the focal plane at the location marked with an arrow in (A). (C) Layer regions that are used to calculate the mean attenuation coefficient for each layer. The scale bars are 200 µm in each dimension. (D) Confocal function corrected focus dependent attenuation coefficient of the different layers indicated in (C) calculated using the depth resolved method and the nominal scattering coefficients (Table 2) indicated with the dashed lines.
4.5 In vivo data
The in vivo OCT images from two normal human retinas were calculated by averaging 120 B-scans for every focus setting using the registration information from the SPECTRALIS system. The data was processed in the same way as the layered TiO2 phantom images, extracting the curved focal planes and the Rayleigh length, correcting for the confocal function, and calculating attenuation coefficient maps. Figure 11 shows an example of the resulting attenuation coefficient maps for subject 1 and 2 for a focus position in the lateral center of the scans 200 µm and 147 µm below the inner limiting membrane (ILM), respectively. The different retinal layers visually show comparable attenuation coefficients for the two subjects, however, the map of subject 1 shows lower values in the center of the B-scan, where the retinal nerve fiber layer (RNFL) is thinnest. Note, that the scan for subject 1 was acquired closer to the fovea than for subject 2 (see Fig. 12(C) and (F)). A-scans with apparent blood flow present (Fig. 12(A) and (D), marked in red) were excluded from further calculations because the strong absorption of blood clearly violates the assumption of a constant backscattering fraction for the depth resolved method.
Fig. 11. Confocal function corrected attenuation coefficient maps for subject 1 (A) and subject 2 (B) acquired with the focus depth in the center of the B-scan 200 µm and 147 µm below the surface, respectively. Scale bars are 200 µm in each dimension.
Fig. 12. Upper part of B-scan of subject 1 (A) and subject 2 (D) where locations with blood flow are marked in red because they are not used for the calculations. Locations of B-scans (A) and (D) are indicated in the SLO images (C) and (F), respectively. Segmented retinal nerve fiber layer (RNFL), ganglion cell and inner plexiform layer (GCL&IPL), inner nuclear layer (INL), outer plexiform layer (OPL), outer nuclear layer (ONL), and photo receptor and retinal pigment epithelium (PR&RPE). Scale bars are 200 µm in each dimension. Confocal function corrected mean attenuation coefficients of the different retinal layers for subject 1 (B) and subject 2 (E) are plotted against the focus setting. The white part indicates the range where the ratio fit exhibits the smallest error and is used to calculate the mean attenuation coefficients in Table 3.
Table 3. Mean attenuation coefficients µt and corresponding standard deviations (std) of retinal layers calculated from the values within the white part of Fig. 12(B) and (E) for subject 1 and 2, respectively.
The mean attenuation coefficients calculated for the different layers marked in Fig. 12(A) and (D) for different focus locations are shown in Fig. 12(B) and (E) for subject 1 and 2, respectively. The segmentation was computed using open-source software [18,19]. Table 3 lists the mean and standard deviation of the attenuation coefficients acquired with central focus depths of 0 µm to 500 µm below the ILM. The µt for subject 1 are mostly slightly lower and show a higher standard deviation than for subject 2.
5. Discussion
We investigated a technique to extract the Rayleigh length and focus depth from a set of OCT B-scans acquired with different focus settings. This technique is based on the method introduced by Dwork et al. [5,6] and extended by Stefan et al. [7]. We applied this ratio fit method to a range of different foci and showed experimentally that the ratio fit works best in a focus depth range from about 0 µm to 500 µm and a focus offset of more than 100 µm by analyzing scattering phantoms of different concentrations. The ratio fit tends to fail if both B-scans used for the ratio calculation are acquired with a focus location above the sample. The obtained focus depth and Rayleigh length were used to correct for the confocal function and subsequently calculate the attenuation coefficients, which agreed well with the nominal scattering coefficients expected from linear extrapolation. The observed deviation of the calculated attenuation coefficients from the nominal scattering coefficients greater than 6 mm−1 can be explained with the increasing presence of multiple and dependent scattering [20–22]. As a reference we also used a method where we directly fitted the confocal function and attenuation coefficient in a single A-scan [4], this has the advantage that no change in focus is needed during the measurement. However, the applicability is limited to homogeneous samples with a scattering coefficient lower than 4 mm−1.
We observed a concentration dependent increase of the apparent Rayleigh length (Fig. 7) when using the ratio fit method. This method is based on the single scattering model and might be increasingly inaccurate with a higher frequency of multiple scattering events. This is expected to happen with both an increased scatterer concentration and at larger imaging depths. The experimental determination of the Rayleigh length using the focus series fit also showed a trend of longer Rayleigh lengths with increasing scatterer concentration, however, this trend is less pronounced than with the ratio fit method. In the focus series fit, multiple scattering effects should be less pronounced because only one location at the top of the sample was analyzed.
The applicability of the ratio fit method was also demonstrated for a layered phantom, where the two-dimensional focal planes could consistently be determined. The reproducibility of the determined attenuation coefficients was excellent, as demonstrated by the equal attenuation coefficients found for the same scattering layers at different depths from the sample surface.
Finally, the ratio fit method was used to extract the confocal function from in vivo retinal scans and calculate attenuation coefficient maps. For the two subjects analyzed we observed a small mismatch between the calculated attenuation coefficients of the retinal layers, however, this might be explained by the different locations of the B-scans relative to the fovea and natural intra- and inter-subject variability. The results on just two subjects seem to indicate that the attenuation coefficient of the RNFL becomes smaller closer to the fovea (Fig. 11(A) and Table 3). Additionally, the less consistent focus dependent attenuation coefficients for subject 1 (Fig. 12(B)) might be explained by accidental accommodation during the measurement making the focus reading of the device inaccurate and possibly also leading to averaging B-scans with different focus locations.
Reliably estimating the attenuation coefficient would add an additional parameter to quantify pathological changes in tissue and help to diagnose those changes earlier. In the field of ophthalmology, it has been reported that the attenuation coefficient might help to distinguish glaucomatous eyes from healthy eyes [23,24] and that diagnostic information could be derived for diabetes [25] and pituitary adenoma [26]. Additionally, extracting the confocal function parameters and subsequently correcting for the axial PSF could be of great value for other quantitative OCT analysis.
For the theory of the ratio fit and the depth resolved attenuation coefficient calculation to apply (as needed for the calculations presented here), we are limited by prior conditions that need to be fulfilled. Explicitly, Eq. (1-6) describing the intensity distribution and the confocal function are only valid if Gaussian optics can be applied, i.e., low numerical aperture and undistorted probe beam are given [15]. Both Eq. (5) (ratio of two A-scans) and Eq. (8) (depth resolved attenuation coefficient profile) are for single scattering only and become more inaccurate with increasing multiple scattering. Furthermore, the depth resolved method can only be applied reliably if the detected light is a fixed fraction of the attenuated light, and if the signal is completely attenuated within the depth of the scan [1]. Some of the mentioned requirements are only partly met for retinal imaging and further refinement is needed.
Funding
Heidelberg Engineering GmbH; Health~Holland, Topsector Life Sciences & Health; Holland High Tech, Topsector High Tech Systems and Materials.
Acknowledgments
We acknowledge funding by Heidelberg Engineering GmbH. The collaboration project is co-funded by the PPP Allowance made available by Health∼Holland, Topsector Life Sciences & Health, and by Holland High Tech, Topsector High Tech Systems and Materials, to stimulate public-private partnerships.
Disclosures
JK: Heidelberg Engineering GmbH (F,E), JFdB: Heidelberg Engineering GmbH (F), JF: Heidelberg Engineering GmbH (E).
Data availability
Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.
Supplemental document
See Supplement 1 for supporting content.
References
1. K. A. Vermeer, J. Mo, J. J. Weda, H. G. Lemij, and J. F. de Boer, “Depth-resolved model-based reconstruction of attenuation coefficients in optical coherence tomography,” Biomed. Opt. Express 5(1), 322–337 (2014). [CrossRef]
2. T. G. van Leeuwen, D. J. Faber, and M. C. Aalders, “Measurement of the axial point spread function in scattering media using single-mode fiber-based optical coherence tomography,” IEEE J. Sel. Top. Quantum Electron. 9(2), 227–233 (2003). [CrossRef]
3. G. T. Smith, N. Dwork, D. O’Connor, U. Sikora, K. L. Lurie, J. M. Pauly, and A. K. Ellerbee, “Automated, depth-resolved estimation of the attenuation coefficient from optical coherence tomography data,” IEEE Trans. Med. Imaging 34(12), 2592–2602 (2015). [CrossRef]
4. B. Ghafaryasl, K. A. Vermeer, J. Kalkman, T. Callewaert, J. F. de Boer, and L. J. Van Vliet, “Analysis of attenuation coefficient estimation in Fourier-domain OCT of semi-infinite media,” Biomed. Opt. Express 11(11), 6093–6107 (2020). [CrossRef]
5. N. Dwork, G. T. Smith, J. M. Pauly, and A. K. Ellerbee Bowden, “Automated estimation of OCT confocal function parameters from two B-scans,” in Conference on Lasers and Electro-Optics, OSA Technical Digest (online) (Optical Society of America, 2016), AW1O.4.
6. N. Dwork, G. T. Smith, T. Leng, J. M. Pauly, and A. K. Bowden, “Automatically determining the confocal parameters from OCT B-scans for quantification of the attenuation coefficients,” IEEE Trans. Med. Imaging 38(1), 261–268 (2019). [CrossRef]
7. S. Stefan, K. S. Jeong, C. Polucha, N. Tapinos, S. A. Toms, and J. Lee, “Determination of confocal profile and curved focal plane for OCT mapping of the attenuation coefficient,” Biomed. Opt. Express 9(10), 5084–5099 (2018). [CrossRef]
8. N. Nassif, B. Cense, B. H. Park, S. H. Yun, T. C. Chen, B. E. Bouma, G. J. Tearney, and J. F. de Boer, “In vivo human retinal imaging by ultrahigh-speed spectral domain optical coherence tomography,” Opt. Lett. 29(5), 480–482 (2004). [CrossRef]
9. M. Wojtkowski, R. Leitgeb, A. Kowalczyk, T. Bajraszewski, and A. F. Fercher, “In vivo human retinal imaging by Fourier domain optical coherence tomography,” J. Biomed. Opt. 7(3), 457–463 (2002). [CrossRef]
10. J. F. de Boer, C. E. Saxer, and J. S. Nelson, “Stable carrier generation and phase-resolved digital data processing in optical coherence tomography,” Appl. Opt. 40(31), 5787–5790 (2001). [CrossRef]
11. D. Hillmann, T. Bonin, C. Luhrs, G. Franke, M. Hagen-Eggert, P. Koch, and G. Huttmann, “Common approach for compensation of axial motion artifacts in swept-source OCT and dispersion in Fourier-domain OCT,” Opt. Express 20(6), 6761–6776 (2012). [CrossRef]
12. S. Yun, G. Tearney, B. Bouma, B. Park, and J. de Boer, “High-speed spectral-domain optical coherence tomography at 1.3 mum wavelength,” Opt. Express 11(26), 3598–3604 (2003). [CrossRef]
13. R. Michels, F. Foschum, and A. Kienle, “Optical properties of fat emulsions,” Opt. Express 16(8), 5907 (2008). [CrossRef]
14. D. M. de Bruin, R. H. Bremmer, V. M. Kodach, R. de Kinkelder, J. van Marle, T. G. van Leeuwen, and D. J. Faber, “Optical phantoms of varying geometry based on thin building blocks with controlled optical properties,” J. Biomed. Opt. 15(2), 025001 (2010). [CrossRef]
15. D. Faber, F. van der Meer, M. Aalders, and T. van Leeuwen, “Quantitative measurement of attenuation coefficients of weakly scattering media using optical coherence tomography,” Opt. Express 12(19), 4353–4365 (2004). [CrossRef]
16. M. J. Girard, N. G. Strouthidis, C. R. Ethier, and J. M. Mari, “Shadow removal and contrast enhancement in optical coherence tomography images of the human optic nerve head,” Invest. Ophthalmol. Vis. Sci. 52(10), 7738–7748 (2011). [CrossRef]
17. P. H. Tomlins, O. Adegun, E. Hagi-Pavli, K. Piper, D. Bader, and F. Fortune, “Scattering attenuation microscopy of oral epithelial dysplasia,” J. Biomed. Opt. 15(6), 066003 (2010). [CrossRef]
18. S. J. Chiu, X. T. Li, P. Nicholas, C. A. Toth, J. A. Izatt, and S. Farsiu, “Automatic segmentation of seven retinal layers in SDOCT images congruent with expert manual segmentation,” Opt. Express 18(18), 19413–19428 (2010). [CrossRef]
19. P.-y. Teng, “Caserel-an open source software for computer-aided segmentation of retinal layers in optical coherence tomography images,” Zenodo. 10.5281/zenodo.17893 (2013).
20. J. Kalkman, A. V. Bykov, D. J. Faber, and T. G. van Leeuwen, “Multiple and dependent scattering effects in Doppler optical coherence tomography,” Opt. Express 18(4), 3883–3892 (2010). [CrossRef]
21. V. D. Nguyen, D. J. Faber, E. van der Pol, T. G. van Leeuwen, and J. Kalkman, “Dependent and multiple scattering in transmission and backscattering optical coherence tomography,” Opt. Express 21(24), 29145–29156 (2013). [CrossRef]
22. M. Almasian, N. Bosschaart, T. G. van Leeuwen, and D. J. Faber, “Validation of quantitative attenuation and backscattering coefficient measurements by optical coherence tomography in the concentration-dependent and multiple scattering regime,” J. Biomed. Opt. 20(12), 121314 (2015). [CrossRef]
23. J. van der Schoot, K. A. Vermeer, J. F. de Boer, and H. G. Lemij, “The effect of glaucoma on the optical attenuation coefficient of the retinal nerve fiber layer in spectral domain optical coherence tomography images,” Invest. Ophthalmol. Vis. Sci. 53(4), 2424–2430 (2012). [CrossRef]
24. K. A. Vermeer, J. van der Schoot, H. G. Lemij, and J. F. de Boer, “RPE-normalized RNFL attenuation coefficient maps derived from volumetric OCT imaging for glaucoma assessment,” Invest. Ophthalmol. Vis. Sci. 53(10), 6102–6108 (2012). [CrossRef]
25. D. C. DeBuc, W. Gao, E. Tátrai, L. Laurik, B. Varga, V. Ölvedy, W. E. Smiddy, R. Tchitnga, A. Somogyib, and G. M. Somfái, “Extracting diagnostic information from optical coherence tomography images of diabetic retinal tissues using depth-dependent attenuation rate and fractal analysis,” in Biomedical Optics and 3-D Imaging, OSA Technical Digest (Optical Society of America, 2012), BTu3A.74.
26. M. Sun, Z. Zhang, C. Ma, S. Chen, and X. Chen, “Quantitative analysis of retinal layers on three-dimensional spectral-domain optical coherence tomography for pituitary adenoma,” PLoS One 12, e0179532 (2017).
