1. Introduction

Optical coherence tomography (OCT) is a technique that is widely used to diagnose pathologies through morphological changes in tissue structure. The standard OCT imaging of tissue morphology based on the back-scattered intensity provides qualitative images that are not quantitatively related to the tissue’s local optical properties. Since certain pathologies affect the optical parameters of tissues, such as the attenuation coefficient, quantification of the attenuation coefficient could be useful for medical diagnosis [1]. Accurate and robust depth resolved models for quantifying tissue optical parameters from OCT measurements are therefore needed. Recently Girard et al. and Vermeer et al. developed a method to quantify the depth resolved attenuation coefficient [2,3]. Such models should not only take into account tissue specific properties but also the effects of detection geometry [3]. Many OCT systems are single mode fibre (SMF) based and therefore experience a so-called confocal effect [4,5]. Earlier work by van Leeuwen et al. [5] derives an expression for the confocal function, $\left (1 + \left ( \frac {z - z_f}{\alpha z_R}\right )^2\right )^{-1}$, that depends on the Rayleigh range ($z_R$), focus position ($z_f$) and the type of reflection that happens, with $\alpha = 1$ for (specular) mirror reflections and $\alpha = 2$ for diffusely reflecting surfaces. The extended Huygens-Fresnel (EHF) model for OCT measurements by Thrane et al. [6,7] contains a confocal function implicitly in the scattered and unscattered beam waist parameters which change with depth and depend on $z_f$ and $z_R$. The scattered beam waist parameter additionally depends on the scattering properties of the probed medium.

In this work we present a single scatter model based on Huygens-Fresnel theory where the scattering by small particles is described by spherical waves. This model, in contrast to [5], predicts $\alpha = 1$ for both specular and diffusely reflecting surfaces. We verified our model by extensive focus series measurements on a mirror, diffusely scattering Spectralon surfaces as well as on the surfaces of different dilutions of Intralipid samples. The effect of scattering on the confocal function was further investigated by measuring the confocal function and its fitted Rayleigh range as a function of depth with respect to sample surface for a series of Intralipid dilutions. Finally, the confocal function with its parameters extracted from the focus series measurements was used in combination with the exponential decay from Lambert-Beer’s law to extract the attenuation coefficient from OCT depth profiles (A-lines) of different Intralipid dilutions.

2. Theory

In fibre based OCT systems, the collection efficiency of the backscattered light depends on the overlap integral of the incident E-field and the modes of the single mode fibre (SMF). Given that the intensity distribution and phase front of a backreflected or backscattered E-field changes with the depth position of the reflection/scatter plane, the collection efficiency of the SMF will be depth dependent. In a free space setup the overlap integral between backscattered and SMF mode fields would be replaced by an overlap integral between the reference arm and backscattered E-fields at the detector. A useful measure of the collection efficiency of a SMF is the power transmission coefficient $T$. $T$ can be calculated with an overlap integral according to [8,9]

Now we will consider a setup where the same SMF acts both as the source of the illumination beam and as the detector of the backreflected beam. The SMF mode can be approximated by a scalar Gaussian beam which is focused by a lens resulting in a Gaussian beam in the sample arm with focus position $z_f$, Rayleigh range $z_R$ and a beam waist radius in the focus of $w_0$, see Appendix Eq. (15). In the sample arm a mirror is placed perpendicular to the beam optical axis at plane $S$ and using energy conservation (no limiting apertures) the power transfer coefficient is evaluated after the lens at plane $D'$ instead of at the SMF surface (plane $D$), see Fig. 1. The $E_{beam}$ field in Eq. (1) can then be determined by propagating the Gaussian beam from the SMF to the mirror at $S$, reflecting it and then propagating it to $D'$. Instead of using the known propagation of Gaussian beams to determine the beam at $D'$, we take the field at $S$ and use Huygens-Fresnel to propagate it to $D'$ to make the comparison for incoherent backreflection later on more clear. Defining the positive $z$-direction from $S$ to $D'$, the beam just before the reflection would be $E_{SMF}(\boldsymbol r_{\perp S}, -z)$ and just after reflection $E_{SMF}(\boldsymbol r_{\perp S}, z)$. So the power transfer coefficient in the coherent case becomes

 

Fig. 1. Sketch of SMF detection system. Beam from the SMF in solid lines, the backreflected (incident) beam indicated by the dashed lines. If there are no apertures the energy conservation principle can be used to evaluate the field overlap integral at any $D'$ along the optical axis instead of at $D$. Note an axial offset, distance between beam foci, of $2(z - z_f)$.

Download Full Size | PDF

To calculate $T$ for a scattering or diffusely reflecting object instead of a mirror we assume that all fields on the backscattering surface $S$ have uncorrelated phases. Therefore each point on surface $S$ contributes an E-field uncorrelated with the other contributions from the surface and the coupled intensity becomes the intensity of each individual point (incoherent sum), integrated over the backscattering surface. As in the coherent case we use the Huygens-Fresnel principle to propagate the E-field of a point at the scattering surface $S$ to the detector surface $D'$. However, instead of propagating the field from scatter surface $S$ to $D'$ where the power transfer coefficient is calculated, each point from scatter surface $S$ has its power transfer coefficient calculated at $D'$ individually before the integration over $S$, resulting in an incoherent sum over the scattering surface $S$. Therefore this case of an incoherent sum over $S$ results in the following equation for the power transfer function

$T_{coh}$ and $T_{inc}$ are identical up to a prefactor, which means we do not find a factor of 2 in front of $z_R$ in the fraction, contrary to what was previously described by van Leeuwen [5] et al. who used a different model to describe the coupling coefficient for a scattering or diffusely reflecting object.

In van Leeuwen’s model a diffusely reflecting object is treated similarly to the mirror by using the same formula from Yuan et al. [9], and therefore assuming a perfectly coherent backscattered Gaussian beam, with the difference that this backscattered Gaussian beam has its focus at the scatter plane with a waist equal to that of the sample beam at the scatter plane. This implies that not only do the scatterers reflect in a perfectly coherent mirror-like manner, also the scatterers transform any sample beam wavefront into a perfectly flat wavefront upon backscattering.

3. Methods and setup

3.1 Overview of the confocal function determination principle

We performed measurements to determine a) the confocal function for coherent and diffuse reflecting surfaces and b) the confocal function inside scattering media as a function of the depth location inside the media. For the first set of measurements, the reflective surface was translated along the optical axis and the backscattered intensity from the surface was measured as a function of the focus position with respect to the sample surface. We fitted Eq. (3) to the measurements to determine the confocal parameter $z_R$. For the second set of measurements, solutions with different concentrations of Intralipid in water were prepared. Depth profiles were measured for each sample as a function of the focus position by translating the sample along the optical axis. The Rayleigh range $z_R$ for a particular depth inside the sample could then be determined by fitting Eq. (3) to the measured backscattered intensity at a particular depth with respect to the surface as a function of the translation of the focus. The advantage of this approach is that since the backscattered intensity is determined at a fixed depth, the attenuation of the beam propagating back and forth to this depth due to scattering and absorption is the same for all measured points in the fit, and therefore will cancel in the determination of $z_R$. A second advantage is that $z_R$ can be determined for different depths inside the sample, and the effect of scattering on $z_R$ as a function of depth can be determined.

3.2 OCT system

Measurements were performed with a SMF based swept source (SL132120, Thorlabs) OCT (SS OCT) device with a central wavelength of 1300 nm, bandwidth of 90 nm, a sweep rate of 200 kHz and a depth scan range of 8 mm in air. In the sample arm a beam is outcoupled from the SMF (SMF-28-J9, Thorlabs), collimated with a reflective collimator (RC04APC-P01, Thorlabs), directed by a Galvo mirror scanning system (GVS002, Thorlabs) and finally focused (SL50-3P - Scan Lens, Thorlabs) resulting in a beam with $z_R = 560$ $\mu \text {m}$ and focus spot radius $w_0 = 15$ $\mu \text {m}$. Without the lens the collimated beam directly after the Galvo mirror scanner unit, see Fig. 2, was calculated to have a radius of 1360 $\mu \text {m}$ which was close to the measured radius of 1386 $\mu \text {m}$ with a beam profiler (BP209-IR/M, Thorlabs).

 

Fig. 2. Swept source PS-OCT system setup adapted from [10]. The Swept source laser providing a "k-clock" signal for A-line triggering and a laser light output that gets split in two arms of a Mach-Zehnder interferometer split in a reference arm with a transmission delay line and a sample arm with a polarisation dispersion unit which splits the sample arm light in two polarisation states, P and S, and introduces a delay between the two states. A circulator (C) redirects the sample arm light to a reflective collimator after which the collimated light enters a Galvo mirror scanning unit that scans the collimated beam over a scan lens which leads focused beam scanning on the sample. Sample arm light back reflected from the sample is directed via the circulator to the polarisation diversity detection module where it interferes with the reference arm light and gets split in the orthogonal P and S polarisation components which get detected on separate balanced detectors. PC are polarisation controllers. One polarisation state in the sample arm is blocked for OCT measurements without polarisation sensitivity.

Download Full Size | PDF

3.3 Angle scanning Galvo mirrors

Accurately measuring the backscattered intensity into a SMF from a mirror while translating the mirror over a few millimeters is difficult. Slight angular misalignment between the sample arm mirror and optical axis of the system will strongly affect the backcoupled intensity. This will not only reduce the power transferred to the SMF but will also vary with the sample arm mirror axial position and therefore affect the confocal function measurement. To ensure a perpendicular beam on the sample arm mirror, the Galvo mirror scanning system was adjusted. Conventionally the scanning mirrors are positioned in the back focal plane of the scan lens so that when the scan mirror angle is changed, the beam scans laterally after the scan lens with the beam parallel to the optical axis of the setup. By moving the scan mirrors slightly out of the scan lens focus, the angle of the beam after the scan lens becomes a function of the deflection angle of the scanning mirrors. For a deflection angle of $\phi$, a focal length of $f$ and an axial offset of the scan mirrors of $\Delta z$, the beam angle after the scan lens becomes $\phi _s = -\frac {\Delta z}{f} \phi$ in the small angle approximation. Then, even with a small sample mirror angular misalignment the maximum coupling coefficient can be found from an area scan (C-scan) of the mirror. The SMF power transfer is highest when no angular sample misalignment is present [9], thus the peak of the Galvo angle C-scan can be used to extract the correct value. The mirror reflection intensities for a Galvo angle C-scan are shown in Fig. 3. For each mirror axial position the peak of a paraboloid fit to the angle C-scan intensities is used as the highest backreflected intensity.

 

Fig. 3. Back reflected Intensities mirror surface for Galvo mirrors angle scan. Peak of plotted intensity surface corresponds to confocal function intensity for mirror axial position of this angle scan measurement.

Download Full Size | PDF

3.4 OCT system roll-off

To accurately measure the confocal function, the system roll-off, i.e., the decrease in sensitivity as function of depth, needs to be compensated for. A common source of roll-off in swept source OCT systems is the finite instantaneous coherence length $l_{ci}$ of the wavelength swept laser. Due to the narrow instantaneous linewidth of our source, the system has an $l_{ci}$ in the order of 30 cm leading to no roll-off over our measurement range of 8 mm. While the system has no roll-off in the optical signal, there is roll-off introduced by the system electronics due to a low-pass filter (Mini-Circuits SBPL-200+) in the signal path from the detector to the digitizer board to suppress aliasing of signals above the Nyquist frequency. Taking into account the aliasing effect, the measured signal $I(z)$ over a depth range $0 \leq z \leq Nq$, with $Nq$ the Nyquist frequency (or depth) can be described by true sample signal $A(z)$, low-pass filter function $f(z)$ and a Nyquist frequency (or depth) $Nq$ as

Since $f(z)$ decays rapidly around the Nyquist depth $Nq$ the higher order alias terms will be small and are neglected. Dividing signal $I(z)$ by filter function $f(z)$ will overestimate the true sample signal $A(z)$ by a factor of 2 at $z = Nq$. Noting that for $z = 0$ we have $f(z) >> f(2Nq - z)$ and therefore $I(z) \approx f(z)A(z)$ and that for $z=Nq$ we get $I(z) = 2 f(Nq)A(Nq)$. We therefore approximate the measured signal as

The filter function $f(z)$ was determined by measuring the background signal with a signal analyser (PXA Signal Analyzer N9030 3 Hz-26.5 GHz) with the filter in place and dividing this by the background signal from the signal analyser without filter, which was then fitted with the response function for a second order Butterworth filter,

The signal, filter function fit and aliased filter function ($f(z) + f(2Nq - z)$) are shown in Fig. 4, with fit parameter $b = 0.001$ mm$^{-4}$ in the depth domain and $b = 6 \cdot 10^{-10}$ MHz$^{-4}$ in the frequency domain.

 

Fig. 4. OCT system roll-off due to electronic filter, in blue the filter response measured with a signal analyser, in orange the filter function (Eq. (9) ) fit to the signal analyser curve and in green the sum of the fitted filter function with its alias. Response given both as function of signal frequency (bottom) and penetration depth (top) of the OCT system

Download Full Size | PDF

3.5 Data acquisition and processing

Three types of samples were used to investigate the confocal function: a mirror for coherent specular reflections, 2% and 99% Spectralon samples for surface diffusive reflections and Intralipid 20% dilutions for back-scattering at different depths. To measure the confocal function, focus series measurements were performed by first positioning the sample surface in the focus of the beam, then increasing the sample distance by 2 mm resulting in the beam focus being positioned 2 mm above the sample surface. Then the sample was translated up in steps of 0.2 mm with a translation stage, moving the focus towards and into the sample and measuring depth profiles at each step over a range of 4 mm. The measurement was repeated by moving the sample in the reverse direction over 4 mm, thus measuring twice at each sample position. At each sample position the mirror was scanned with the angle scanning method described in Section 3.3. The Spectralon samples were scanned in both lateral directions by parallel beams, 500 A-lines per B-scans and 500 B-scans per C-scan, resulting in 250000 A-lines which were averaged into a single A-line to reduce speckle effects and also resulting in an ensemble averaged A-line. The scanner assembly was rotated by an angle of $4^{\circ }$ with respect to the vertical to suppress specular reflections from the Intralipid sample surface. Scans were performed along a line perpendicular to the sample inclination, 500 A-lines per B-scan. Since the Intralipid samples have Brownian motion, speckle could be averaged by repeating the B-scans at the same location 500 times to again calculate an (ensemble) average A-line from 250000 A-lines. The averaged A-lines from both detectors, each measuring a polarisation orthogonal to the other, were further processed according to the following steps. First, a Hann window was applied to suppress sidelobes that would be introduced by a Fourier transformation of a finite time interval signal. Second, numerical dispersion compensation was performed. Third, zero padding was performed to increase the A-line pixel density by a factor of 2. Fourth, the signals were Fourier transformed and the mirror images resulting from the Fourier transform were removed. Fifth, the roll-off was compensated for as explained in Section 3.4. And finally sixth, the signals of the 2 detectors/polarisation states were squared and then summed to get the total intensity signal.

3.6 Properties of Intralipid samples

To asses the effect of different scattering strengths on the confocal function, a dilution series of Intralipid 20% with distilled water was made. Using the wavelength dependent expression for the scattering coefficient ($\mu _s$) of stock solution Intralipid 20% from Aernouts et al. [11]

Table 1. Properties Intralipid dilutions compared to stock Intralipid 20%. Independent scattering values based on curve from [11,14]. Dependent scattering values based on curve from [14]. For attenuation coefficient ($\mu _t = \mu _a + \mu _s$) add 0.11 mm$^{-1}$ to the $\mu _s$ values.

View Table | View all tables in this article

3.7 Confocal function determination

To determine the confocal function from focus series measurements of backreflected intensities from sample surfaces, Eq. (11) is fitted, with $z_s$ the position of the surface as the independent variable and with the normalisation factor (peak value) $A$, focus position $z_f$ and Rayleigh range $z_R$ as free parameters.

To determine the confocal function inside the scattering media at fixed distances $\Delta z$ from the sample surface, the procedure described above needs to be revised to take into account the shift of the focus position to $z_f'$ by refraction of the incident beam at the air-sample interface due to a change in refractive index. The shift in focus position depends on the distance between the air-sample interface $z_s$ and the unshifted focus position $z_f$ and the sample refractive index $n$. The shifted focus position is given by $z_f' = z_f + \Delta z_f$, where $\Delta z_f = (n-1)(z_f - z_s)$. See Fig. 5 and note that the beam focus position in the confocal function is with respect to the intensity being backscattered inside the sample, so when a negative focus shift ( $z_s > z_f$) occurs, the (now virtual) focus position $z_f'$ of the beam propagating in the sample should still be used, and not the physical and unshifted focus position ($z_f$) located in air. Given a focus series measurement at a fixed distance inside the sample ($z' = z_s + \Delta z$) and how the (possibly virtual) $z_f'$ shifts depending on $z_s$, the confocal function can be rewritten as

So for a fixed depth position $\Delta z$ with respect to the surface the independent variable will be $n(z_s - z_f)$. Finally the Rayleigh range of a Gaussian beam is proportional to the refractive index so the confocal fits inside the sample are expected to produce $z_{R, \text {intralipid}} = nz_{R, \text {air}}$.

 

Fig. 5. Shift in focus due to jump in refractive index. The solid line represents actual beam and dashed line the beam if it did not experience refraction. The focus position in air indicated by $z_f$, the position of air sample interface by $z_s$ and the focus position of the refracted beam by $z_f'$

Download Full Size | PDF

3.8 Attenuation coefficient fit OCT depth profiles

Having determined the confocal function parameters $z_f$ and $z_R$, the attenuation coefficient ($\mu _t = \mu _a + \mu _s$) can now be extracted based on the singly scatter OCT model for a homogeneous medium given by

Using this as fit function for each averaged depth profile for each focus position and sample, with given beam focus position in the sample $z_f$ and Rayleigh range of the beam inside the sample $z_R$ and distance from the surface $z$ as independent variable. The normalisation constant $A$, attenuation coefficient $\mu _t$ and noise floor $C$ can be fitted for and extracted. To determine the attenuation coefficient from an OCT depth profile we need to transform the OCT depth profiles as a function of optical pathlength (OPL) to depth profiles as a function of depth. This is done by first determining the positions of all the interfaces where a jump in $n$ occurs; in our case this is only at the air-Intralipid interface. Then based on the interface positions we assign each depth profile pixel with a refractive index $n_i$. The physical depth size of pixel $i$ can then be calculated using the known system OPL axial resolution $\Delta l$ via $\Delta z_{i} = \frac {\Delta l}{n_i}$. Finally the depth of the depth profile is then simply $z = \sum _{i} \Delta z_{i}$.

4. Results

4.1 Confocal function of mirror and Spectralon samples

The depth profiles from a focus series measurement for both the Spectralon 2% and 99% with the confocal function, Eq. (11), fitted on the surface reflection signals are shown in Fig. 6 A and B. Note that at the surfaces there are no specular reflection peaks visible. The penetration depth for Spectralon 2% is much lower compared to the 99% due to the strong absorption of Spectralon 2%. The normalised confocal function from surface measurements for a mirror, Spectralon 2% and 99% measurement are shown in Fig. 6(C). It clearly shows the mirror and Spectralon 2% confocal functions nearly overlapping and the Spectralon 99% confocal function being significantly wider compared to the confocal function of the other 2 samples. For each of the three samples three focus series measurements were performed and the weighted mean of the $z_R$ extracted from these fits are shown with standard deviation in Table 2. The fitted $z_R$ of the mirror and Spectralon 2% are close to one another as well as the theoretical beam $z_R$ (564 $\mu \text {m}$), while the fitted $z_R$ for Spectralon 99% sample is about 1.4 times the theoretical beam $z_R$.

 

Fig. 6. A: The (averaged) depth profiles of a focus series measurements of a Spectralon 2% sample with the dashed line the fitted confocal function. B: The same as A for a Spectralon 99% sample. Figure C shows the extracted surface reflections for the Spectralon 2%, Spectralon 99% and mirror sample with blue circles, orange triangles and green crosses, respectively. For each sample the confocal is fit drawn in the same colour as their respective marker, the intensities have been normalised based on the peak value from the confocal function fit.

Download Full Size | PDF

Table 2. Beam $z_R$ based on theoretical calculations and weighted mean of the confocal fit $z_R$ for 3 focus series measurements with standard deviation.

View Table | View all tables in this article

4.2 Confocal function of Intralipid dilution samples

For the Intralipid dilution measurements the confocal function was determined for the signal from the sample surfaces and as a function of depth inside the the samples, so at each position along the depth profiles. For the confocal function measurements inside the sample the change in refractive index needed to be accounted for as described in Section 3.7. In Fig. 7 the averaged depth profiles of the focus series measurements are shown for the different Intralipid dilutions with the confocal function fit, Eq. (11), on the surface reflections shown as black dashed lines. Although specular reflections where suppressed by scanning the samples under an angle of $4^{\circ }$, some samples still show specular reflection as can be seen by the additional peaks in the depth profiles at the sample surfaces. This is especially visible for the samples with lower scattering coefficient, such as the I01 sample ($\mu _s = 0.1$ mm$^{-1}$), due to a relatively lower signal contribution from light scattered from within the sample. As mentioned in Section 3.5 all sample positions are measured twice in a focus series measurement, a forward and a backward axial sweep. The strength of the specular reflection changed between the axial sweeps for I025 and I2, necessitating separate confocal fits on the forward and backward sweeps for these measurements. In Fig. 7 only a single sweep is shown for each sample for clarity.

 

Fig. 7. Averaged depth profiles of a focus series measurement for each Intralipid sample with a black dashed line indicating the confocal function fit on the backscatter intensities from the sample surface. Depth in mm from the zero optical path length difference with OCT reference arm.

Download Full Size | PDF

To determine the confocal function inside the Intralipid samples, Eq. (12) was used with as fit parameters the peak value $A$, $z_R$, and $nz_f$ with $z_f$ the focus position in air. Each fit was done at fixed depth $\Delta z$ with respect to the sample surface, while sample surface position $z_s$ of the averaged depth profiles of the focus series measurement was the independent variable. The results of the confocal function measurements for 3 different depths for each Intralipid sample are shown in Fig. 8 with the sample scatterer concentration increasing from top to bottom. The left plots of Fig. 8 show the averaged depth profiles for each sample with their surfaces aligned, visualising fixed $z_s$ and variable $z_f$. The vertical dashed lines indicate the depth positions at which the confocal functions were evaluated. The right hand plots of Fig. 8 show normalised intensities of the averaged depth profiles at the 3 depths as a function of distance with respect to the beam focus position in Intralipid. The beam focus positions in Intralipid were calculated from the sample surface positions and beam focus position in air, see Section 3.7, and normalised based on the peak of the fitted confocal functions. The shallowest depth is indicated by blue circles, the middle depth by orange squares and the deepest depth by green crosses. The dashed lines indicate the fitted confocal function, each with colour corresponding to the marker colour of that depth. For the top plots, which have lower scatter coefficient, the normalised confocal functions for the different depths overlap. For the bottom 2 plots, I4 and I10, which had the highest scatter coefficients, the normalised confocal functions become wider with depth and for the deepest depth of I10 the confocal function becomes almost a flat line. Our single scatter based theory for the confocal function predicts no change in the confocal function with depth, however a broadening of the normalised confocal function with depth is observed for samples with these high scatterer concentrations.

 

Fig. 8. Left Intralipid dilution (averaged) depth profiles of focus series measurements with the surfaces of the depth profiles within a focus series aligned, dashed lines indicate depths from surface at which the confocal functions were evaluated, sample scatterer concentration increases from top to bottom. On the right the normalised intensities at the different depths as a function of distance to the beam focus position in the Intralipid dilutions, the blue circles are from the shallowest depth, the orange triangles from a deeper depth and the green crosses from the deepest depth. The dashed lines indicate confocal function fits, their colours are identical to the colour of the data point they are fitted to. The distances from surface are also indicated in the depth profile plots on the right following the same colour coding.

Download Full Size | PDF

4.3 Fitted $z_R$ confocal functions

While the confocal function derived in this work predicts a half width half maximum (HWHM) equal to the beam Rayleigh range for both coherent (specular) and diffuse reflections, work by van Leeuwen [5] predicts a factor of 2 increase in confocal function width in case of diffuse reflections. To determine which is correct the fitted $z_R$ of all the different samples are compared in Fig. 9. For the fitted $z_R$ measured inside the Intralipid a weighted average of the $z_R$ extracted over a depth range from 100 to 200 $\mu \text {m}$ was calculated. This averaging region for the depth fitted $z_R$ was chosen to avoid the regions affected by specular reflection and depth dependence as best as possible. Figure 9 shows that for the surface reflection measurements (crosses), the coherent (mirror) and all the diffuse reflection measurement (Spectralon and Intralipid), the fitted $z_R$ are close to the expected beam $z_R$ of $564 \ \mu \text {m}$ (dashed line). The only exception is the Spectralon 99% measurement with a fitted $z_R$ that is about 1.4 times the expected beam $z_R$. This could be due to significant contribution of multiple scattered light whose effect on the confocal function (coupling efficiency) is not taken into account by the model. The circles in Fig. 9 show the $z_R$ resulting from fits of back scattered light from within the medium for the Intralipid samples. Due to a change of refractive index the expected beam $z_R$ inside Intralipid is $750 \ \mu \text {m}$ as indicated by the second dashed line. The Intralipid samples with lower scatterer concentration (I01, I025, and I05) have a fitted $z_R$ inside the sample that are close to the expected beam $z_R$, however the higher scatterer concentration samples have higher fitted $z_R$ that increases with the scatterer concentration of the Intralipid dilutions.

 

Fig. 9. The fitted confocal function $z_R$ for different samples, crosses are measurements on surface reflections while circles indicate confocal function measurements inside samples. The dashed lines indicate the beam Rayleigh range in air and Intralipid, $564 \ \mu \text {m}$ and $750 \ \mu \text {m}$ respectively. The fitted $z_R$ for inside the samples are weighted averages from a depth range of $100$ to $200 \ \mu \text {m}$.

Download Full Size | PDF

This could be attributed to changes in the beam shape due to attenuation [15,16], or multiple scattered light [6,7], which both are not accounted for in the model. From these observations we conclude that the width of the confocal function in the single scatter approximation depends in the same manner on the Gaussian beam parameter $z_R$ for both coherent and diffuse reflecting surfaces and Intralipid dilutions. Additionally the width of the confocal function increases as multiple scattered light becomes more significant.

4.4 Depth dependence of confocal function $z_R$

To better visualise the broadening of the confocal function as function of depth, the width (fitted $z_R$) of the confocal function at each depth for each Intralipid dilution is shown in Fig. 10(A). This figure shows that at the surface all Intralipid samples have a fitted $z_R$ close to the expected beam Rayleigh range in Intralipid of $750 \ \mu \text {m}$ indicated by the black dashed line. For Intralipid dilutions with $\mu _s < 1$ mm$^{-1}$ the fitted $z_R$ remains constant with depth, while for dilutions with $\mu _s \geq 1$ mm$^{-1}$ the fitted $z_R$ shows an increase as the depth increases.

 

Fig. 10. In A) the Rayleigh range extracted from the confocal function fits as function of depth, the black dashed line indicates the Rayleigh range of the beam in Intralipid $z_R = 750 \ \mu \text {m}$. In B) The same as in A) but as a function of the (dimensionless) number of scatter mean free paths given by $2z \mu _s$.

Download Full Size | PDF

Since the confocal function derived in this work is based on single scattering and the fraction of detected multiple scattered light increases with sample scatterer concentration and depth, it would be useful to have the fitted $z_R$ as function of a parameter that includes both these effects. To that end we introduce a dimensionless quantity called the scatter optical depth or the depth expressed as a number of mean free paths $\lambda _{mfs} = 1/\mu _s$, which is a measure for the amount of (multiple) scattered light and defined as $\int \mu _s(l) dl = 2z \mu _s$, with $l$ the light path and $\mu _s(l)$ the scatter coefficient along that path. Since each sample has a constant $\mu _s$ it simplifies to depth of the signal multiplied with the scatter coefficient with a factor of 2 due to the round trip. The fitted $z_R$ as function of this scatter optical depth is shown in Fig. 10(B). The increase of $z_R$ as a function of scatter optical depth shows a universal behaviour with a steady increase up to 6 mean free paths. The curves for samples with measurements past 6 mean free paths (I2, I4, I10) each show a different increase in slope causing the curves to diverge with respect to each other.

The universal behaviour of the fitted $z_R$ curves for the different Intralipid dilutions when plotted as a function of the scatter optical depth indicates that the fitted $z_R$, and therefore the confocal function, are dependent on the amount of (multiple) scattering. For lower scattering optical depths, and therefore lower fractions of multiple scattered light, the change in the confocal function due to multiple scattered light can be approximated as a linear broadening of the confocal function for scatter optical depth (SOD) $\leq 6$. Since the Rayleigh range is a measure for the amount of (de)focusing of a Gaussian beam, the broadening of the confocal function (increase in fitted $z_R$) could be ascribed to broadening of the beam radius due to scattering inside the sample as predicted by Andersen [6,7] in case of multiple forward scattered light.

At higher scatter optical depths (SOD $\geq 6$) however the divergence of the fitted $z_R$ curves for the different Intralipid dilutions suggest a break down of the model confocal function, leading to erratic behaviour in the fitting.

4.5 Attenuation coefficient fits OCT depth profile Intralipid dilutions

The aim is to be able to extract attenuation coefficients from single focus position measurements with known $z_f$ and $z_R$. Therefore in this section each averaged Intralipid dilution depth profile is fitted with Eq. (13) with a fixed focus position $z_f$ and a fixed $z_R$ of the beam inside the Intralipid dilutions. The $z_f$ were calculated from the surface confocal fits and the depth profile surface positions, and $z_R$ from the depth averaged fitted $z_R$ from Fig. 9. The resulting fits for a selection of averaged depth profiles for all of the Intralipid samples are shown in Fig. 11, with the dashed lines indicating the fit curves. For the curve fitting procedure the first few pixels from the surface were ignored to remove the effect of the surface specular reflections on the fit, as can be most clearly seen for I05 and I1. Ignoring these specular reflections the fitted curves closely overlap the averaged depth profiles for all samples except for I10. For I10 the fitted curve has trouble capturing both initial slope right after the surface and the change in slope at around 1 mm in depth from the surface. This was to be expected as the fits assume a constant confocal function which was shown to be depth dependent for more strongly scattering media in Section 4.2. In fact it is surprising how well the I2 and I4 curves fit given that its fitted $z_R$ also increased with depth. It is therefore expected that, while the I2 and I4 fits do seem to correspond well with the measured curves, the extracted $\mu _t$ for I2 and I4 are unreliable as the fitting procedure allows for adjusting $\mu _t$ to compensate for an inaccurate confocal function model. The same issues are at play for I10 but there it is easily visible that the fits are inaccurate, i.e. the OCT depth profile model is a poor description of the actual OCT depth profile signal for high scattering.

 

Fig. 11. Averaged depth profiles Intralipid samples with black dashed lines indicating OCT depth profile fits based on the singly scattered model with confocal function, where $\mu _t$ was fitted for. On the x-axis the depth positions in OPL with respect to the the OCT reference arm.

Download Full Size | PDF

The attenuation coefficients ($\mu _t = \mu _a + \mu _s$) from the fitted curves are shown in Fig. 12(A) with the nominal $\mu _t$ for each sample indicated by a dashed line with a colour corresponding to the samples marker colour. The nominal $\mu _t$ for I10 ($\mu _t = 10.11$ mm$^{-1}$ ) is not shown in the figure due to it being outside the range of the vertical axis. The Intralipid dilutions with low scatter coefficients ($\mu _s \leq 1$ mm$^{-1}$) have fitted $\mu _t$ that are close to the nominal values of those dilutions and their fitted $\mu _t$ are constant as a function of the focus position of the beam. However, Intralipid dilutions with higher scatter coefficients ($\mu _s \geq 2$ mm$^{-1}$) have fitted $\mu _t$ that underestimate the nominal $\mu _t$ values. Additionally the fitted $\mu _t$ for focus positions above the sample surface are lower compared to those with a focus position inside the medium for these samples. To get a single $\mu _t$ value for each sample a weighted average of the $\mu _t$ of all the different focus position depth profiles was calculated with weights being equal to one over the variance as returned by the Scipy function scipy.optimize.curve_fit. The results are shown in Fig. 12(B) with a black dashed line indicating equal averaged fitted $\mu _t$ and nominal $\mu _t$ for the different Intralipid dilutions. This figure shows even more clearly what could be seen from the previous plot, averaged fitted $\mu _t$ for Intralipid dilutions with $\mu _s \leq 1$ mm$^{-1}$ are in agreement with the nominal value, but for dilutions with $\mu _s \geq 2$ mm$^{-1}$ the averaged fitted $\mu _t$ underestimates the nominal $\mu _t$ and this underestimation gets larger with dilution scatter coefficient.

 

Fig. 12. In A) fitted $\mu _t$ of the averaged depth profiles from the Intralipid dilutions, on the horizontal axis the distance from the sample surface to the focus position of the beam inside the sample for each averaged depth profile. The dashed lines indicate the nominal $\mu _t$ for the samples with the colour of the lines corresponding to the colour of the markers. Note that the nominal $\mu _t$ for I10 is not shown as it is outside the plot range. In B) For each sample the $\mu _t$ averaged from the different focus position depth profiles, on the horizontal line the nominal $\mu _t$ and the vertical axis shows the averaged $\mu _t$ from the fits. The black dashed line indicates line of equal nominal and averaged fitted $\mu _t$.

Download Full Size | PDF

Based on these results, Eq. (13) is a good description for OCT depth profiles of homogeneous media in case of low scatter coefficients (single scatter approximation) but breaks down as the scatter coefficient, and therefore multiple scattered light, increases. This failure of the OCT depth profile model is partly due to the inaccurate confocal function for these samples with higher $\mu _s$ (I2, I4 and I10), as in Section 4.4 a depth dependent fitted $z_R$ is shown while the model assumes a $z_R$ constant with depth. Finally, it is also expected that the Lambert-Beer decay in the model becomes inaccurate as $\mu _s$ and consequently the multiple scattering contribution increases.

5. Discussion

We have derived a confocal function model for reflections from specular and diffusely reflecting surfaces. The model predicts (except for a prefactor) the same confocal function with identical dependence on $z_R$ for both specular and diffuse reflections. This contradicts earlier work by van Leeuwen et al. [5] which predicted a $2z_R$ dependence for diffusely reflecting surfaces. Extensive focus series measurements in this work on a mirror, Spectralon 2% and 99% reflectivity and on Intralipid dilution solutions as a function of depth confirmed the prediction of our model for weakly scattering media and low number of mean-free scatter paths. The fitted $z_R$ increased with $\mu _s$ and number of mean-free scatter paths, also observed by Kübler et al. [4], until the model failed when $2 z \mu _s > 6$. This could be explained by changes in beam shape due to attenuation [15,16] or multiple scattered light [6,7], which are unaccounted for in the model.

To investigate the accuracy of OCT attenuation coefficient extraction using a single scatter Lambert-Beer model including our confocal function, averaged depth profiles of Intralipid dilutions for different focus positions and different sample scatterer concentration were measured. While there was a good overlap between the measured depth profiles and the fits for all samples except for I10, the fitted $\mu _t$ for the I2 and I4 samples are unreliable due to the confocal function model being inaccurate for these more strongly scattering samples. A discrepancy between the nominal $\mu _t$ and fitted $\mu _t$ for samples I2, I4 and I10 is therefore expected. For the less strongly scattering samples the nominal $\mu _t$ and fitted $\mu _t$ were in good agreement.

Nominal $\mu _s$ values of the Intralipid dilutions were based on the assumption of independent scattering, i.e. a linear relation between scatterer concentration and $\mu _s$, using Aernouts et al. [11,14] curve for independent scattering. For high intralipid concentrations, dependent scattering provides a correction to the linear relation between concentration and $\mu _s$ as described by Aernouts [14], changing the nominal $\mu _s$ of I2, I4 and I10 from 2, 4 and 10 mm$^{-1}$ to 1.8, 3.2 and 5.5 mm$^{-1}$, respectively. The nominal $\mu _s$ values for dependent scattering would bring the I2, I4, and I10 points in Fig. 12(B) closer to the fitted values. However the confocal function model was demonstrated to be inaccurate for these samples as is evident from the increase of the fitted $z_R$, affecting the fitted $\mu _t$. In addition, we observed that using the independent scattering curve over the dependent scattering curve for the dimensionless x-axis in Fig. 10(B) shows a much clearer universal behaviour. We hypothesise that OCT-based determination of the effect of dependent scattering can be partially contributed to the failure of the confocal model for high scattering.

Robust and accurate quantification of the attenuation coefficient of tissue using OCT, enabled by appropriate modelling of the confocal function, could be used to quantify pathological changes in tissue and help to earlier diagnose these changes. For example, determination of the tissue attenuation coefficient could allow for detection of atherosclerotic plaques [17,18] and distinguish cancerous from healthy tissue [19–21]. In ophthalmology the attenuation coefficients of the retinal tissue layers could help distinguishing between healthy and glaucomatous eyes [22,23] as well as help with the diagnosis of diabetes [24] and pituitary adenoma [25]. The model for the confocal function in SMF OCT systems developed and demonstrated in this paper should be of use to further investigate the clinical value of quantification of the tissue attenuation coefficient for these and other applications.

The presented confocal function model is limited to weakly scattering media and lower number of mean-free scattering paths due to beam profile changes resulting from attenuation [15,16] and multiple scattering [6,7] unaccounted for within the model. Additionally, since both the beam profile changes due to attenuation and multiple scattering are beam-NA dependent the $\mu _s$ and the number of mean-free scattering paths at which the model becomes inaccurate is expected to be NA dependent. Finally, relative increases in forward scattering, or a scattering anisotropy $g$ closer to 1, increases the fraction of multiple scattered light within the OCT signal, reducing the accuracy of the single scatter approximation.