1. Introduction

Optical coherence tomography (OCT) is now used routinely as part of ophthalmic clinical practice [1] and is approaching this status in a variety of other areas of clinical practice such as cardiology [2] and dermatology [3]. Alongside these developments, there is an immense amount of work being applied in earlier stage research into OCT technical and application development [4]. Consistent with other medical imaging techniques such as X-ray, magnetic resonance and ultrasonographic imaging, the development and clinical application of OCT can be assisted by accurate models of image formation. This aids image interpretation, particularly where features of interest are at, or near, the resolution limit of the OCT system. Models of image formation are also crucial in the design and optimisation of imaging systems as well as the development of image processing and quantification techniques.

Models of image formation have, until recently, fallen broadly into one of two categories: those based on a point spread function (PSF) formalism and those based on the Monte Carlo method. Image formation, in the case of a single scatterer, is readily calculated using theory established for modelling coherent microscopes [5]. Modelling of image formation in turbid media, as is usually the case in imaging of biological tissue, has been achieved using different approximations. Perhaps the simplest of these is the first-order Born approximation, where the light incident upon a particular refractive index perturbation is assumed to be unperturbed by the sample. This approximation thus allows image formation to be described by the linear superposition of unperturbed PSFs and has been employed in the context of inverse scattering [6]. Coupland and Lobera [7] developed a general formalism of image formation in OCT, using rigorous scalar diffraction theory, which treats image formation as a linear filtering operation. Villiger and Lasser [8] demonstrated an equally rigorous formalism based on the notion of coherent transfer functions. The extended Huygens-Fresnel formalism [9, 10] extends the first-order Born approach by including the scattering induced, depth-dependent, OCT signal attenuation and the contribution due to multiple scattering. A key aspect of this approach is that these extensions are based upon scatterer ensemble-average contributions, rather than deterministic scatterer distributions. Furthermore, neither the change in coherence of light due to propagation in tissue, nor the interference of scattered and reference light are explicitly modelled by PSF based approaches.

There are several models of OCT image formation (for example [11–16]) based upon the Monte Carlo method for modelling light propagation in biological tissue [17]. These models have revealed much about OCT image formation and Monte Carlo modelling is considered to be the gold standard technique in some branches of biomedical optics. Despite this, Monte Carlo based models possess some limitations when used to model image formation in OCT. Monte Carlo methods represent tissue by its spatially resolved, statistically averaged, properties, which are assumed to “extend uniformly over small units of tissue volume” [17]. Furthermore, although some effort has been made along these lines, wave properties such as polarisation, coherence and interference are not naturally treated using the particle formalism intrinsic to the Monte Carlo method.

The limitations of both PSF and Monte Carlo models can be overcome by full wave models based upon Maxwell’s equations. To our knowledge, the first such comprehensive model was recently introduced [18], which is limited to two-dimensional optical systems and samples. A full wave analytic model for simulating image formation for single cylindrical scatterers was also recently introduced [19]. A similar model for spheres has also been introduced [20], however this model is limited by approximations such as treating the focused illumination, incident upon the sphere, as plane wave illumination. Full wave modelling of optical imaging finds its origins in modelling image formation in confocal microscopes, as was necessary due to the high numerical aperture lenses employed. Image formation models based on Maxwell’s equations were introduced for point scatterers [21], finite sized spheres [22] and for general samples [23], the latter of which is the inspiration for this current work. There have been steps toward full wave modelling of OCT image formation in two-dimensions which have fallen short of calculating image formation. For example, the finite-difference time-domain method has been used to probe the electromagnetic field within the sample [24–28], but aside from other approximations, these models do not consider the detection of the interference between scattered and reference light, only the distribution of scattered light within the sample.

The three-dimensional full wave model presented in this paper is inspired by one developed for modelling image formation in coherent optical microscopes [23] and extends one developed for OCT that was limited to two dimensional optical systems and samples [18]. The two-dimensional model could not be trivially extended to three dimensions, as the addition of a spatial dimension would result in a computationally infeasible model. This is because a three-dimensional model must consider volumes of tissue which are large compared with other applications of full wave electromagnetic modelling. Furthermore, OCT employs broadband light, thus adding a dimension to the full wave simulation, over and above the three spatial dimensions. The model presented in this paper has been adapted from the two-dimensional model through two principal developments. The first is an extension to the memory-efficient pseudospectral time-domain (PSTD) method enabling arbitrary focused beams to be employed efficiently and accurately [29]. The second is a method for efficiently calculating the detection of scattered and reference light which is embedded in the PSTD simulation [30].

In the remainder of this paper we will describe the full wave model’s key components: the illumination, scattering and detection sub-models. We will also describe how these components are integrated together to form a complete model. Some examples of image formation are given, and we evaluate the accuracy of the model. Finally, the paper is concluded by considering the model’s weaknesses, areas for improvement and outlook for application by the field.

2. The model

2.1. Outline

We model the system depicted schematically in Fig. 1, which is representative of spectral domain and swept source OCT systems. In the case of spectral domain systems, an ultra-broadband light source is employed and the detector is actually a spectrometer. In the case of swept source systems, a swept source laser is employed along with a detector that has no spectral resolution. As is explained in further detail below, owing to second-order coherence phenomena, these two systems can be treated identically within the framework of the mathematical model that we present.

 

Fig. 1 Schematic diagram of the modelled OCT system. ϕ_s is the electric field produced by focusing the fiber mode into the sample space and E_sc is the field scattered by the sample back towards the objective lens. n_s is the refractive index distribution of the sample.

Download Full Size | PDF

The model is composed of a sequence of algorithms as depicted schematically in Fig. 2, each of which is closely associated with a physical aspect of an OCT system. The model is vectorial in nature and so polarisation phenomena are implicitly modelled. The coordinate systems used in the following discussion are illustrated in Fig. 3, which is a schematic diagram of the focusing optics. Note that although this optical system is depicted as a 4f system, it is understood that this is an approximation, since scanning mirrors will generally be present in the focal plane common to the collimator and objective lenses in Fig. 3. Despite this, study of the OCT system used to generate the experimental results in this paper reveals that this is a sound approximation. It is, however, possible to modify this component of the model in cases where this approximation is not valid.

 

Fig. 2 Flow diagram illustrating the principal components of the imaging model.

Download Full Size | PDF

 

Fig. 3 Schematic diagram and notation of the optical focusing system. Each position vector r_i = (x_i, y_i) denotes a transverse position in the space denoted in the diagram. R_a is the physical radius of the aperture of the system, z_obs is the axial location of plane where the scattered field is sampled and n_b is the refractive index of the material in which the sample is embedded, which has an interface with air at z₃ = −h.

Download Full Size | PDF

We now briefly discuss how the model treats coherence. The model is most analogous to a system in which a pulsed laser illuminates the sample and the light scattered by the sample is sensed by a spectrometer, after being coupled into the optical fiber. Using only the theory of linear time-invariant systems, neglecting the finite size of the spectrometer’s pixel, this is equivalent to calculating the light coupled into the fiber for a set of coherent wavefields, one for each wave number within the system’s spectral range, which is in turn equivalent to a swept source system. Spectral domain systems differ in that they employ broadband sources, which emit a partially coherent wavefield, that is described by a stationary random process. Since, however, we are ultimately interested in the correlation between light from the sample and reference arms, the coherent-mode representation [31] allows us to model the illumination as a superposition of coherent modes, as in the swept source case. Note that it is only because we ultimately measure spectral density that these two systems are mathematically equivalent.

We now outline the principal components of the model as depicted in Fig. 2. As just explained, it is thus valid, for the purposes of this model, to consider that the optical fiber emits a coherent mode with wave number k which we denote φ_f (r₁). This field is assumed in this paper, without loss of generality, to be polarised in the x-direction. The model calculates the vectorial electric field due to this mode being focused into the sample space using vector diffraction theory [32], denoted ϕ_s (r₃, z₃, k). We allow for the case where the sample may be embedded within a bulk material with refractive index n_b which has an interface with air at z₃ = −h as shown in Fig. 3. Having calculated the light incident upon the sample, defined by a spatially resolved refractive index distribution, n_s (r₃, z₃), the light scattered by the sample is calculated using the PSTD method [33]. The final step is to calculate how the scattered field is coupled back in to the fiber. This is achieved using a specially designed algorithm embedded within the PSTD algorithm [30] which alleviates, without approximation, the need to explicitly propagate the scattered field through the optical system. This step results in a modal amplitude for both the scattered (α_sc (k)) and reference (α_ref (k)) light, which can be used to calculate an OCT A-scan. We now describe each of these steps in further detail.

2.2. Calculation of incident illumination

Low numerical aperture (NA) objectives, such as those used in OCT, do not significantly alter the polarisation state of light. A vectorial treatment of the incident light is thus unnecessary from this perspective. However, there is interest in how samples alter the polarisation state of light. Furthermore, the PSTD method solves Maxwell’s equations, thus requiring the source field to satisfy Maxwell’s equations. We thus employ the Debye-Wolf formalism [34, 35] which calculates the field, satisfying Maxwell’s equations, in the vicinity of the focus of an objective lens, which may be of arbitrary NA.

By employing the Gaussian approximation for weakly guiding optical fibers it is reasonable to assume that the mode of the single mode optical fiber can be represented as a polarised wavefield. We assume, without loss of generality, that the field is linearly polarised in the x-direction. The x-component of the field at the exit face of the fiber is thus assumed to take the form

2.3. Calculation of light scattered by the sample

Having calculate the mathematical form of the light incident upon the sample, the light scattered by the sample, E_sc (r₃, z₃, k), must then be calculated. Formally, E_sc (r₃, z₃, k) + ϕ_s (r₃, k) (E_sc is the scattered field and E_sc + ϕ_s is the total field) is the field which satisfies Maxwell’s equations and the constitutive relations for a sample refractive index distribution n_s (r₃, z₃) and source field ϕ_s (r₃, k). Since we are interested in simulating image formation for general samples, we employ a numerical method to solve for E_sc (r₃, z₃, k). Our previous two-dimensional model employed the finite-difference time-domain (FDTD) method [36, 37] to calculate the scattered field. The FDTD method discretizes space on a grid with typical spacing of approximately λ/20 (λ is the minimum wavelength of interest) and Maxwell’s equations are expressed as difference equations. This fine spatial sampling requirement means that the FDTD method is only tractable in the two-dimensional case when considering typical OCT imaging depths. Consider, for example, simulating the formation of an A-scan for a sample 1mm deep. The transverse physical size of the FDTD simulation will typically be at least 50μm depending on the resolution of the system being modelled. Supposing a source with center wavelength λ = 800nm is employed, with a grid spacing of λ/20, approximately 1.7 terabytes of computer memory and infeasible computation time would be required to simulate a single A-scan.

In order to overcome this computational limitation we have employed the PSTD method [33, 38], which allows space to be discretized on a grid approaching a spacing of λ/2. The spatial derivatives inherent to Maxwell’s equations may then be expressed using derivatives evaluated in the Fourier domain. In particular, the derivative of a band limited function, g(x), sampled on a grid of spacing Δ, according to the Nyquist criterion, resulting in a sequence g_j = g(jΔx) of length N, may be evaluated according to $g_j^{\prime }={\mathcal{D}}^{-1}\left\{i2\pi a_j/N\mathcal{D}\left\{g_j\right\}\right\}$ where N is assumed to be even, $\mathcal{D}$ and ${\mathcal{D}}^{-1}$ are the discrete Fourier and inverse Fourier transform operators and:

Thus, using a grid spacing of λ/2 would require three orders of magnitude less computer memory (a tractable 1.7 gigabytes) to simulate the above A-scan. Like the FDTD method, the PSTD method, being a time domain simulation, is able to evaluate the scattered field for the full range of wavelengths required to simulate A-scan formation. This is vital since this model would be computationally intractable if a separate simulation were required to be performed for all wavelengths. We employed the FFTW implementation of the fast Fourier transform in this work [39].

Irrespective of whether the FDTD or PSTD methods are employed, the incident illumination, ϕ_s (r₃, z₃, k), calculated in Eq. (2), is introduced into the simulation at a plane z₃ = z_s, by converting the time-harmonic field ϕ_s (r₃, z_s, k) into a time-variant magnetic current density of the form [29]:

The PSTD method imposes some restrictions upon the model. In particular, the spatial resolution of the sample’s refractive distribution is dictated by the simulation’s grid spacing. Thus, if a grid of spacing λ/2 is employed, the sample can only be represented by cubes of side λ/2, which must have constant refractive index. An anisoptropic grid spacing may, however, be employed. A perfectly matched layer, which absorbs outgoing radiation with very low reflection, must be placed around the simulation space in order to approximate unbounded scattering. We also note that the plane z₃ = z_s, at which the incident illumination is introduced, is limited by the finite size of the PSTD grid, thus leading to some truncation of the incident beam. The latter two limitations can be made insignificant [23] whilst the first limitation is the subject of ongoing work.

2.4. Calculation of scattered light coupled back in to the optical fiber

As explained in Sec. 2.3, the PSTD method directly yields the scattered field in the time-domain. For clarity, however, we initially consider E_sc (r₃, z₃, k), the time-harmonic form of the scattered field, before explaining how the light coupled into the fiber is calculated directly from the time-domain field data. The wave number specific modal coefficient of scattered light coupled into the optical fiber could be evaluated by storing E_sc (r₃, z_obs, k) for each PSTD simulation and evaluating E_sc (r₁, k) using imaging theory. The modal coefficient for the light coupled in to the fiber would then be given by:

We now consider how Eq. (8) may be evaluated directly from time-domain data as the PSTD simulation progresses. In particular, assume that we have access to U_sc (r₃, z_obs, nΔt), and thus ${\tilde{U}}_{sc}\left({\mathit{q}}_3,z_{obs},n\mathrm{\Delta }t\right)$, through the application of a spatial Fourier transform, at time step n for a PSTD time step of Δt. The time-harmonic field would normally be evaluated using a discrete Fourier transform as $U_{sc}\left({\mathit{r}}_3,z_{obs},k\right)=\frac{1}{N_t}{\displaystyle {\sum }_{n=0}^{N_t-1}{U}_{sc}\left({\mathit{r}}_3,z_{obs},n\mathrm{\Delta }t\right)\exp \left(ikcn\mathrm{\Delta }t\right)}$, where the summation is updated at each of a total of N_t time steps of the PSTD simulation. It is then easy to see that if this definition of $U_{sc}\left({\mathit{r}}_3,z_{obs},k\right)\left(\text{or}{\tilde{U}}_{sc}\left({\mathit{r}}_3,z_{obs},k\right)\right)$ is used in Eq. (8), the summation over the time step (n) can be moved outside the integral over r₃ (or q₃). We note that it is numerically more efficient to use the first integral in Eq. (8) since the k dependence of P is mathematically simpler. We also note that the term $ℱ\left\{P\left({\mathit{q}}_3,z_{obs},k\right){\tilde{\varphi }}_f\left(\frac{f_2}{f_1}{\mathit{q}}_3\right)\right\}$ in Eq. (8) is equivalent to the appropriate component of the incident field discussed in Sec. 2.2

2.5. Calculation of A-scan

Having now calculated α_sc (k) and α_ref (k), it now remains to evaluate the detector current as:

2.6. Mitigation of numerical dispersion

Close attention must be paid to numerical dispersion when using the PSTD method. For completeness, we consider both the PSTD and FDTD methods. As in physical wave propagation, numerical dispersion leads to waves having wavelength-dependent group and phase velocities. The analytic expression of the FDTD numerical dispersion relationship for a plane wave propagating in a homogeneous medium is [37]:

Numerical dispersion can be reduced by employing higher order time stepping procedures [37], which, however, require additional computer storage and computation. Since this application requires the conservation of computing resources, we employ an approach not requiring additional computational resources. This is achieved by choosing the angular frequency, ω, at which the quantities ${\alpha }_{sc}(\tilde{k})$ and ${\alpha }_{ref}(\tilde{k})$ are sampled, such that the numerical wave number, $\tilde{k}$, matches the desired wave number, k. This is able to be achieved since the temporal discrete Fourier transform discussed in Sec. 2.4 can be performed at values of ω not given by kc, but by the appropriate dispersion relationships, Eqs. (11) and (12). This guarantees that, for a homogeneous medium, the correct wave numbers are sampled. This is suitable for OCT in biological tissue, since we are generally interested in modelling scatterers embedded in a largely homogeneous medium. This approach remains valid for the PSTD case even for inhomogeneous volumes since $\tilde{k}c$ is independent of refractive index. A small error results when using FDTD, however, as is shown in Sec. 3, this error is negligible when refractive index contrasts compatible with biological tissue are employed. We note however, that FDTD is rarely employed due to its dense spatial sampling requirements.

3. Results

All of the simulations presented here share common parameters based upon the experimental system, in particular, a Thorlabs Telesto-II spectral domain OCT system with an LSM03 objective. On the basis of discussions with Thorlabs, the Telesto-II was approximated by the optical system shown in Fig. 3 using parameters detailed in Table 1. A spectrum centered on wavelength λ₀ = 1300nm with a 170nm bandwidth was employed throughout the simulations. In the interests of accuracy, we note that after completing the simulations presented in this paper, it was revealed that the Telesto-II in fact has an effective spectrum, combining the influences of the source spectrum, spectrometer bandwidth and spectral shaping, of approximately 240nm wide. This 240nm spectrum will be used for future simulations when comparison with experimental images will be performed. The PSTD simulations employed a Yee cell of size λ₀/4 in the transverse directions and λ₀/6 in the axial direction. The simulation employed 180 Yee cells in each of the lateral dimensions and 1660 cells in the axial direction. The transverse simulation size was chosen to ensure that the incident beam and detection algorithm operate correctly. These criteria are in fact equivalent [30]. The main restriction upon the transverse simulation size is that the incident beam should enter the PSTD grid without being significantly truncated by the transverse limits of the PSTD grid. At the same time, the transverse size should be kept as small as possible to minimize the computational complexity of the simulation. It is also desirable for the number of cells in each dimension to be equal to a product of small prime numbers as this influences the computational efficiency of the fast Fourier transform employed [39]. This was able to be achieved for the transverse dimension of 180, but not for the axial dimension. Each A-scan evaluated using the PSTD method took approximately 23.5 hours to compute using a 12-core, Intel Xeon E5-2690 V3 “Haswell” processor.

Table 1. Base parameters of the numerical simulations. Note that MFD stands for mode field diameter and symbols are defined in Fig. 3.

View Table

3.1. Ideal point spread function

This first simulation acts as both a verification and exemplar of the model. It is well understood that the lateral PSF of an OCT system deteriorates the further a scatterer is from the focus of the objective. In this example we have calculated 15 independent PSFs at different axial locations as demonstrated in Figs. 4(a)–4(c). In particular, Fig. 4(a) shows a maximum intensity projection depiction of the OCT system’s beam, Fig. 4(b) shows the scatterers which were embedded within the PSTD simulation to generate the 15 PSFs shown in the maximum intensity projection of the simulated OCT image shown in Fig. 4(c). The OCT system was focused on the first scatterer along the optical axis. The scatterers had the size of a single Yee-cell with refractive index 1.425, chosen so as to not perturb the OCT beam such that each of the 15 PSFs can be considered to have been simulated in isolation from the others. The C-scan is built up by executing one PSTD simulation per A-scan, each of which considers a different transverse displacement of the scatterers illustrated in Fig. 4.

 

Fig. 4 Volume, surface, line and error plots of unaberrated, depth dependent, PSFs.

Download Full Size | PDF

Figure 4(d) is a surface plot where the surface height represents the magnitude of the OCT signal as a function of scatterer location in the xz plane. The red plotted lines have been calculated using focusing theory [32] at the center wavelength. The red line linking each of the PSFs is the so-called confocal function of the system and is proportional to the irradiance of the focused beam. It has been plotted to show that the PSTD based model implicitly includes the confocal effect. The other lines plotted along various values of constant z are also proportional to the incident beam’s irradiance, which is equivalent to the transverse PSF of a coherent scanning microscope [40]. This result was calculated as a check that the PSTD model was implemented correctly and could have been calculated using, for example, the extended Huygens-Fresnel formalism. An error was calculated for each of the 15 transverse PSFs displayed in Fig. 4. In particular, if we denote by z_j, the axial location of the jth scatterer, there is one red line plot on the lower axis of Fig. 4 at each z = z_j. An error value was then evaluated as:

3.2. OCT B-scan of a sphere

We calculated the OCT B-scan of a finite sized sphere since the OCT image may be calculated analytically. The general computational model was thus validated against this analytic model of OCT image formation for a single sphere, based upon Mie’s solution for scattering by dielectric spheres, and is a modified version of a previously reported analytic model applicable to image formation in confocal microscopes [41, 42]. The model used in this study has been modified in two principal ways. Firstly, detection of the light scattered by the sphere is performed in the sample space as discussed in Sec. 2.3. Secondly, a numerically efficient scheme for scanning the sphere, or equivalently the beam, has been implemented in a manner similar to Brenner et al. [19]. We note that a model of OCT image formation for single spheres has already been reported [20], however, as discussed in Sec. 1, this model is limited by approximations such as treating the focused illumination, incident upon the sphere, as plane wave illumination.

We briefly describe the analytic model. Consider Eq. (2) which describes how the incident illumination is calculated. Incorporation of a phase term, $\exp \left(-ik{\tilde{\mathit{r}}}_2\cdot {\mathit{r}}_s\right)$, in the kernel of Eq. (2) allows for the incident beam to be translated in the sample space. Each point within the aperture, r₂, gives rise to a plane wave, in the sample space, which takes the form:

The strategy for calculating the field scattered by a sphere for focussed illumination, is to integrate the scattered fields due to each plane wave expressed in Eq. (15). The sphere’s rotational symmetry allows the scattered field for a plane wave, with arbitrary angle of incidence, to be calculated using an implementation of Mie’s series for a single angle of incidence, combined with coordinate system transformations [42]. Thus, if we denote by ${\mathit{e}}_{\text{Mie}}\left({\tilde{\mathit{r}}}_2;{\mathit{r}}_3\right)$, the field scattered by a sphere at the origin, evaluated at r₃, for an incident plane wave of unit amplitude, with propagation and polarisation vectors in common with ${\mathit{e}}_p\left({\tilde{\mathit{r}}}_{\mathbf{2}},{\mathit{r}}_s,{\mathit{r}}_3\right)$, the scattered field due to the focussed beam in Eq. (15) can be evaluated as ${\mathit{E}}_{sc}({\mathit{r}}_3)={\displaystyle {\int }_{|{\tilde{\mathit{r}}}_2|\le {\text{NA}}_2}A({\tilde{\mathit{r}}}_2)\exp (ik{\varphi }_s({\tilde{\mathit{r}}}_2,{\mathit{r}}_s)){\mathit{e}}_{\text{Mie}}}({\tilde{\mathit{r}}}_2;{\mathit{r}}_3){\mathrm{d}}^2{\tilde{\mathit{r}}}_2$, where the dependence of both E_sc(r₃) and ϕ(r₃) on k is not explicitly noted for brevity. It now remains to evaluate α_sc (k) using Eq. (5) projected into the sample space as:

Once α_sc (k) is evaluated, it only remains to evaluate α_ref (k) before the OCT A-scan can be evaluated according to Eq. (10). α_ref (k) may be evaluated using an approach similar to Eq. (16) but with the term $(\widehat{\mathit{i}}\cdot {\mathit{e}}_{\text{Mie}}({\tilde{\mathit{r}}}_2^{\prime };{\mathit{r}}_3))$ in Eq. (17) replaced by a reflection of the incident wave by a plane mirror.

We compared image formation for a sphere of radius 10µm and refractive index 1.42 × 1.1 within a medium of refractive index 1.42 (i.e., 10% refractive index contrast) for the analytic and computational models. The imaging system specifications were otherwise identical to that used in Sec. 3.1. The sphere was located with its centre 10µm beyond the focus of the objective. For this test, however, a Yee cell size of λ₀/20 was used to minimise the stair-casing error that exists at the boundary of the sphere. We used the FDTD method instead of the PSTD method in this case, since the small Yee cell permits it, and it is only necessary to simulate just over 20µm in depth. Furthermore, our implementation of the FDTD algorithm is significantly more computationally efficient than the PSTD algorithm for equal numbers of cells in the computational grid. Each of the 21 FDTD simulations making up Fig. 5(b) took approximately 17 hours to execute using dual Intel Xeon E5-2650 v2 processors, each with 8 cores.

 

Fig. 5 (a) Shows a high resolution image of the magnitude of the OCT B-scan calculated analytically. The transparent circle indicates the location of the sphere, and has been stretched in the z direction by 10% to account for the refractive index of the sphere. (b) and (c) show the magnitude of the OCT B-scan calculated using the FDTD and analytic methods, respectively, with equal pixel size. (c) and (d) show plots of the OCT magnitude and real part for an A-scan through the sphere center, for the FDTD and analytic methods, respectively.

Download Full Size | PDF

The results of the comparison are shown in Fig. 5. In order to provide some quantitative meaning to the results, the OCT data has been normalised by the complex OCT signal corresponding to reflection by a dielectric slab of refractive index matching the sphere, placed at the focus. No further normalisation has been applied. The z-axes in Fig. 5 have been scaled by the refractive index 1.42 and not by a numerical group refractive index, when numerical dispersion mitigation is not employed [18]. Figure 5(a) is a high resolution image calculated using the analytic model. It shows that the OCT B-scan of such a sphere has a non-trivial structure. Discussion on this structure’s origin, for cylinders, is given by Brenner at al. [19]. Good agreement is found between the B-scan magnitudes as shown in Figs. 5(b) and 5(c). The most noticeable difference occurs in the region denoted by the grey box in Fig. 5(b). We performed a quantitative comparison between the FDTD and analytic based models by calculating the following error metrics:

Figure 5(d) shows the magnitude of an A-scan through the centre of the sphere, demonstrating a very good correspondence between the two models. It is interesting to note that the signals appearing in the central A-scan, in the vicinity of z = 60µm agree well in terms of amplitude and less well when the real part is considered. This error is due to numerical dispersion and has been exacerbated by increased propagation within the sphere.

3.3. Discretized scatterer design

Tissue mimicking phantoms used in OCT often make use of scatterers such as silica, PMMA and titanium dioxide microspheres to create regions of known scattering coefficient. In Sec. 3.4 we calculate a PSF which has been aberrated by a layer containing scatterers. In order to control the scattering coefficient of the scattering layer, it is necessary to consider how scattering by discretized spheres, which must conform to the computational grid, differ from that of ideal spheres. In this work we consider titanium dioxide microspheres of diameter 1µm as have been used to construct structured phantoms for OCT [44]. Spherical scatterers which are on the scale of the wavelength are not, in general, well represented in numerical solvers of Maxwell’s equations which must all approximate scatterers using finite elements or cells. In our case we employ a Yell cell of size (λ₀/4) × (λ₀/4) × (λ₀/6), where λ₀ = 1300nm, which make approximating a sphere of diameter 1µm difficult. This is illustrated in Fig. 6(b) which depicts an ideal sphere of diameter 1µm and the scale of the PSTD method’s computational grid. The axes in Fig. 6(b) are in terms of PSTD grid index and, thus, such a sphere can be represented only approximately in the PSTD grid. The discretized scatterers in Fig. 6(c) and Fig. 6(d) are composed of cells with indices satisfying {(i, j, k)|(iλ₀/4)² + (jλ₀/4)² + (kλ₀/6)² < R} where R = 0.5µm in case c) and R = 0.66µm in case d). Since we wish to control the scattering coefficient in our numerical simulations, we studied the scattering cross-sections of the discretized scatterers and compared them with that of a sphere using Mie theory. The scattering cross-sections of the discretized scatterer were calculated from first principles using the PSTD method with plane wave illumination. We present the results in Fig. 6(a) in terms of scattering coefficient by assuming the scatterers have a concentration of 4.5×10⁻³ scatterers per cubic micron, as is the design concentration of the physical phantom which inspired this simulation [44].

 

Fig. 6 (a) Plots showing the scattering coefficient (µ_s) for a sphere of diameter 1µm (illustrated for reference in (b)) and the discretized scatterers shown in (c) and (d) as a function of refractive index. The shape corresponding to “Discretized - 1” is shown in (c) and “Discretized - 2” in d). Axis labels in (b)-(d) are PSTD grid indices where i, j and k correspond to indices in the x, y and z-directions, respectively.

Download Full Size | PDF

The plots in Fig. 6a) show that discretized scatterer 1 under approximates the scattering cross-section of the sphere whilst discretized scatterer 2 over approximates the scattering cross-section. Two values of scatterer refractive index have been calculated as n_s = 1.56 and 1.62, as shown by the dots in Fig. 6(a), resulting in scattering coefficients of 2 and 4mm⁻¹, respectively. These values are employed in Sec. 3.4 to create layers of known scattering coefficient. Light scattering is, however, described by more than just the scattering cross-section, which ignores the direction of scattering. To demonstrate this, a component of the time averaged Poynting vector, S, for the two discretized scatterers and a perfect sphere has been plotted in Fig. 7. To achieve a scattering coefficient of 2 or 4mm⁻¹, discretized scatterer 1 must have refractive index n_s = 1.70 or 1.82, respectively, and a sphere must have a refractive index of 1.65 or 1.74, respectively. In each plot in Fig. 7, the scatterer was placed at the origin. We define $\widehat{\mathit{n}}$ as the vector normal to each plane, directed away from the scatterer. The quantity plotted in Fig. 7 is $\mathit{S}·\widehat{\mathit{n}}/I_i$, where I_i is the irradiance of the incident plane wave propagating in the positive z-direction. Integration of this quantity on a surface enclosing the scatterer yields the scattering cross-section. These plots show that, for a given scatterer, the distribution of scattered light does not vary appreciably between the two scattering coefficients, only the magnitude. Perhaps most importantly, scattering by scatterer 2 is more tightly focussed around the forward scattering direction and scatters more light in the direction opposite to the incident waves direction of propagation.

 

Fig. 7 Plots of the component of the time averaged Poynting vector of scattered light directed away from the scatterer, which is centered on the origin, and normal to each plane, normalised by the irradiance of the incident plane wave. The depicted planes thus correspond to forward and back scattered light.

Download Full Size | PDF

3.4. Aberrated point spread function

Having calculated the unaberrated PSF of the simulated OCT system in Sec. 3.1, we now calculate a PSF aberrated by a scattering overlayer. The scattering overlayers had a thickness of 300µm with discretized scatterer 2 randomly arranged throughout the layer at concentration 4.5 × 10⁻³ scatterers per cubic micron. The refractive index of the scatterers was set to 1.42, 1.56 and 1.62, resulting in scattering coefficients of µ_s =0 (i.e. scatterer free), 2 and 4mm⁻¹, respectively. An aberrated PSF was calculated by placing an isolated scatterer just below the scattering overlayer. A C-scan containing 289 A-scans (i.e. 17 A-scans in each lateral direction) was calculated for each scattering coefficient as shown in the top row of images in Fig. 8. Each A-scan corresponds to a particular amount by which every scatterer within the simulated scattering geometry was shifted in a transverse direction. These plots are, to our knowledge, the first rigorously simulated C-scans ever published. The en-face PSFs are plotted on a linear scale in the lower row of Fig. 8. Note that each of these en-face PSFs is plotted on its own linear color axis and the peaks of the aberrated PSFs are 2.2 and 6.5dB below that of the unaberrated peak for the 2 and 4mm⁻¹ cases, respectively. Aside from the reduction in signal, the shape of the PSFs is seen to depart from the rotationally symmetric shape of the unaberrated PSF. This simulation is not possible using any of the existing models of OCT image formation.

 

Fig. 8 Calculation of aberrated PSFs. Each image in the top row shows slices through a C-scan on the same log scale corresponding to scattering coefficients of 0 (i.e. no scatterers), 2 and 4 mm⁻¹ in the aberrating layer. The lower row of images shows en-face planes through the PSF target where each image is displayed on its own linear scale.

Download Full Size | PDF

4. Discussion and conclusions

We have presented a highly realistic OCT image formation model which implicitly includes phenomena such as multiple-scattering and polarisation. We have published what we believe to be the first rigorously simulated OCT C-scans. This model has been made possible by exploiting the PSTD method which allows for light scattering to be calculated for samples which are several hundreds of microns thick. A new method of launching a focussed beam into the PSTD simulation had to be developed [29], along with a numerically efficient method for calculating how scattered light is imaged onto, and coupled into, an optical fiber [30]. The presented model is very flexible in that, as shown in Fig. 2, it is partitioned into several sub-models, each of which may be independently modified to incorporate different phenomena. For example, sensitivity roll-off could be modelled by extending Eq. (10) to consider the spectrometer’s pixel width. Anisotropic materials can be modelled by extending the PSTD method to consider anisotropic materials. Different beam types or abberations can be modelled by incorporating them into Eq. (2), and so on.

We anticipate that this model will enable a variety of pressing questions, regarding OCT image formation, to be addressed. It will be particularly useful for investigating instances where OCT is used to probe tissue possessing features which lie at, or just below, the resolution of the imaging system, thus assisting image interpretation. The model will also be very useful in validating emerging quantitative techniques such as, for example, parametric imaging [45] and cell type differentiation [46].

Since this is the first demonstration of the model, there are several ways in which the model can be improved. One key area is computational speed. The PSTD method is already parallelised for shared memory platforms using openMP, however it may be desirable to implement distributed memory parallelisation using MPI. Implementation on graphical processing units (GPUs) and even GPU clusters promises dramatically reduced computation times and is thus of high priority.

The primary weaknesses in the model are currently the source condition and perfectly matched layer (PML). In particular, the source condition is exact only at the centre wavelength, with only a small error introduced away from the centre wavelength as demonstrated in Sec. 3. This error remains small, however, only at low numerical apertures, and so if the model is to be extended to high numerical apertures, the source condition will need to be improved. This is a computational, rather than fundamental problem, since a solution to this problem has already been developed which increases the computational cost of the PSTD method [47]. The PML was discovered to perform poorly when the Yell cell approaches λ₀/2, which is one reason why the Yee cell in Sec. 3 did not exceed λ₀/6 in the axial direction. Overcoming this will allow for larger samples to be modelled on fixed resource computer hardware. Finally, we note that the refractive index distribution is generally known for phantoms composed of discrete scatterers. This is, however, not the case for biological tissue, where the refractive index distribution is generally not known. It is anticipated that this model can contribute to overcoming this problem by simulating image formation for hypothesized refractive index structures, derived from higher resolution imaging modalities, for comparison with experimental OCT data.