1. Introduction

A variety of methods for simulating optical coherence tomography (OCT) image formation have been developed [1–11]. These methods have been motivated by the need to interpret images, to understand and manipulate phenomena which influence image formation in OCT, and to provide a means of testing new system designs without the need to build them. Such models have also been important in the development of quantitative techniques, such as parametric imaging [12] and solving the OCT inverse problem [6].

All such mathematical models make approximations and the scale of these approximations varies according to a model’s objectives and constraints. Rigorous simulation of image formation in OCT requires the calculation of how a stochastic electromagnetic beam is focused by an optical system, how the focused beam interacts with a sample and, finally, how that light scattered back by the sample is focused onto a detector, where it interferes with the reference light. Further processing in the time or spectral domains is then required to extract axial information from the recorded signals. Such a model cannot, in general, be evaluated analytically for samples described by an arbitrary refractive index distribution, even using mathematically convenient representations for the incident beam and detection optics.

Here, we are interested in models based on a wave theory of light. Most such models developed to date employ a scalar approximation. Some models have been extended to allow for the simulation of vectorial phenomena (e.g., [1–3]); however, the scalar approach is still predominant. The coherence properties of the incident and scattered fields are generally also not explicitly modelled, in the sense that the effect that propagation has on the coherence of the beam is not modelled. Coherence is incorporated either directly through a point spread function (PSF) [4] or by assuming that the light is delta-correlated [5]. Combining the scalar approximation with the first-order Born approximation allows the construction of a model which has been used as a means of solving the linear inverse problem [6]. The first-order Born approximation, however, implies that light scattered back to the detector undergoes only a single interaction with the sample, i.e., the contribution of each scatterer is calculated as though it was located in homogeneous space. Whilst OCT’s coherence gating aids in the satisfaction of this approximation, it does not hold, in general, as is demonstrated by the observation of attenuation and multiple scattering in OCT images.

The extended Huygens-Fresnel formalism has been used to construct a model which partially relaxes the first-order Born approximation [4]. In this model, light which travels to and from a scatterer is attenuated by the traversed part of the sample according to a statistical model. The remainder of the model is essentially based upon linear systems theory, whereby the coherence gate is embedded in a PSF, which is used to calculate the image formed by a deterministic arrangement of scattering potentials. It is important to emphasize that, in this model, samples are represented only in a quasi-deterministic manner since beam attenuation and multiple scattering are based on an ensemble average representation. This useful model has influenced many subsequent models, including those based upon the Monte Carlo method [7, 8], which derive their image-formation model from it.

Mesoscale models have also been developed which bypass wave optics. For example, the radiative transport equation (RTE), which is based upon energy conservation in a particle transport regime, may be used to model how light propagates in a sample. Image formation can be modelled using a result which relates the complex amplitude of light in the sample to the OCT signal at the corresponding depth [9]. An example of such a model makes use of a solution to the RTE for the case of layered media [5]. However, there are very few samples for which analytic solutions are available. One way of overcoming this lack of solutions is to use the Monte Carlo method [10], combined with an image formation model based upon the extended Huygens-Fresnel formalism [7]. The advantage of this approach is that the Monte Carlo method can be used to calculate the propagation of light in an arbitrary sample, so long as that sample is well described by mesoscopic optical properties, which are expected to remain constant over “small units of tissue volume” [10]. The disadvantage of this approach is that the wave nature of light, including polarization and coherence, is not naturally included in the model.

In other branches of optical imaging, an unavoidable need has been the catalyst for the development of rigorous models of image formation. For example, vectorial effects, which are inescapable in high-numerical aperture confocal microscopy, prompted the development of the first rigorous model of such a microscope valid for general samples [13]. The model was made possible by the use of rigorous numerical techniques, such as the finite-difference time-domain method, to calculate how light propagates in general samples. Capoglu et al. [14] published a good review of rigorous computational models of optical microscopes. Few OCT models have used rigorous numerical techniques to model light propagation in the sample and those that do, do not rigorously model image formation. In one of only two examples we are aware of that employ rigorous numerical techniques, Reed et al. [15] report a model of OCT which uses the time-resolved Poynting vector in the sample to approximate the resulting OCT image. In the second example, Hung et al. [16] developed a model in which the impact of surface roughness upon refractive index recovery, for simple scattering geometries, was considered, but examples of OCT image formation were not given.

The purpose of this paper is to introduce a vectorial model, an extension of a previously developed model of coherent optical microscopes [13], by which image formation in OCT can be rigorously simulated. By this, we mean that the simulation is derived directly from the physics of each system component and the sample, applying as few approximations as possible. In this way, physical phenomena, such as different beam types, deterministic heterogeneous samples, polarization, and anisotropic scatterers, are implicitly included in the model. For example, using a vectorial description of light allows polarization-sensitive OCT to be modelled rigorously, and, as foreshadowed by Schmitt and Knüttel [4], the coherence properties of an OCT system require a vectorial model in order to be treated “more exactly”. Further, a vectorial model is required in order to treat scattering as a continuum of interactions rather than a heuristic partitioning between scattered and unscattered light.

The model is intrinsically three-dimensional. The examples in this paper, however, correspond to two-dimensional scenarios due to computational limitations. Removal of these limitations will be the subject of future work, which will enable the model to be applied to three-dimensional samples. In principle, the model presented here is limited only by the phenomena which are able to be modelled by the FDTD method, the set of which is extremely broad [17]. Therefore, we provide a foundation to study, in principle, all phenomena observed in OCT images.

In the remainder of this paper, we give an outline of the model that we propose before giving detailed explanations of its three principal components: 1) Calculation of light incident on the sample; 2) Calculation of how light interacts with the sample; and 3) Calculation of how light scattered by the sample and light from the reference arm interfere at the detector. A number of examples of the application of the model are then given in order to demonstrate its versatility. We describe image formation for a stratified medium, for a numerical point-spread function phantom and for a numerical phantom based on the physical three-dimensional structured phantom used in our laboratory, for which we perform a comparison with experiment.

2. Details of the model

2.1. Overview

We make reference to the diagrams in Fig. 1 to explain how our model is constructed. Although the illumination and detection optical systems are depicted as 4f systems in air, the mathematical formalism applies equally as well to fiber-based systems which may employ optical systems more complex than depicted in Fig. 1. We do not consider noise or system components with non-ideal characteristics, even though it is possible to do so. The model applies equally well to spectral- and time-domain OCT systems due to both the Wiener-Khintchine theorem and the coherent-mode decomposition of fields [18]. We have, however, described the model in the spectral-domain, which seems to be most common in the literature. One of the model’s defining characteristics is that the principal functional components of the OCT system are mathematically decoupled and treated independently, in a similar manner to models of high-NA coherent microscopes [13]. Consider Fig. 1(a): the model’s ultimate objective is to spatially integrate the irradiance, at the detector D, as a function of wavenumber, k, of the field which results from interfering light from the reference and sample arms. This is achieved by the sequential application of three steps:

 

Fig. 1 a) A schematic diagram of a spectral-domain OCT system containing a light source (LS), lenses (L1 – L4), beam splitter (BS), sample (S) and detector (D); b) A typical ray of the geometrical optics field, ε₀(s, k). Phase and amplitude modulation are described by A_i(s) whilst refraction by L2 is described by [L_ij]; c) A 4f system which images the light scattered by the sample onto a detector or into the guided mode of an optical fiber with any aberrations being represented by the phase modulation A_d(s).

Download Full Size | PDF

As shown in Fig. 1, r_s and r_d are positions in the sample and detector space, respectively. We note also that bold face symbols are used to denote vector quantities. The light incident on the detector from the reference arm, E_ref(r_d,k), is calculated in the same way as E_sc(r_d,k), except that the sample is replaced with a mirror modelled as either a dielectric or perfect electrical conductor. When detection is coherent, such as by a single mode fiber, the spectrally resolved detector signal may then be evaluated according to:

2.2. Illumination

Any incident beam with a mathematical description satisfying Maxwell’s equations may be employed. The Debye-Wolf integral is a flexible formulation that has been widely used to calculate the form of coherent beams focused by high-NA lenses [20–22], which remains valid for low NA lenses. It has been applied in a wide variety of applications, for example, in focusing exotic beams such as cylindrical-vector beams [23] and Gauss-Laguerre beams [24], as well as focussing these and more conventional beams through microscope cover slips [25]. A variety of beams such as Gaussian, Bessel and Airy, may be modelled and arbitrary states of polarization are also able to be represented without any additional complexity. With reference to Fig. 1(a) and Fig. 1(b), the Debye-Wolf integral is formally valid for calculating the time-harmonic field in the vicinity of the focus of L2 if L2 and L4 form a 4f system and the aperture stop (not shown in Fig. 1) is placed at their common focal plane. The 4f system will then be afocal and telecentric from the spaces of both the sample and light source. This is important because it means the back focal plane of L2 appears at infinity from the sample space, which in turn means that the field at the back focal plane may be specified using geometrical optics [26]. This enables the incident field to be calculated according to:

2.3. Scattered light

We employ the finite-difference time-domain (FDTD) method to calculate how light is scattered by the sample. The advantage of using this method is that scattering, at all wavelengths of interest, may be calculated in a single simulation. It is beyond the scope of this paper to give a detailed account of these well established methods; instead, we refer readers to the comprehensive book by Taflove and Hagness [17].

The sample is modelled by associating with each cell in the FDTD grid, a relative permittivity and permeability, conductivity and, possibly, a dispersion relationship. When employing the FDTD method, a grid cell size not exceeding λ₀/10 is required; however, a smaller size may be employed subject to the amount of computer memory available. When computer memory is limited, the pseudospectral time-domain (PSTD) method may be used, which requires significantly less computer memory than the FDTD method, as a minimum grid spacing of λ₀/2 may be approached; however, our recent work [29] shows that a discretization not exceeding λ₀/2.5 is required in practice.

A pulse with a Gaussian temporal profile is introduced into the FDTD simulation space at a plane in the vicinity of the the sample surface, such as S_inc in Fig. 2. This plane must be large enough such that the majority of the beam’s energy flows through it to prevent the non-physical diffraction of the focussed beam. The incident field is introduced to the simulation using equivalent electromagnetic sources [30], which, in the case of FDTD, is equivalent to the well known total-field scattered-field source condition [31]. The properties of the source condition employed in our model are most easily understood by considering the case in which only equivalent magnetic current sources are employed. In particular, we assume that a sheet of magnetic current density is located at plane S_inc in Fig. 2. This magnetic current density has a time-dependent form [29]:

 

Fig. 2 Schematic diagram of where the source field is introduced, S_inc, and where the scattered field is recorded, S_sc, in FDTD simulations. PML stands for perfectly matched layer.

Download Full Size | PDF

The complex amplitudes of the scattered field are recorded, at all wavelengths of interest, on a plane also in the vicinity of the sample surface, such as S_sc in Fig. 2. This is done by applying a temporal discrete Fourier transform at each wavenumber of interest as the FDTD simulation proceeds. The resulting complex amplitudes are equivalent to those that would be obtained by performing individual steady state simulations for each wavelength of interest [17]. Scattering in an infinite space is simulated by surrounding the sample with a perfectly matched layer, as proposed by Berenger [27], which absorbs, without reflection, radiation incident upon it. We have arrived at this treatment by tacitly assuming that the light source is well modelled as a point source. If this is not the case, the coherent mode representation may be used to treat beams with arbitrary states of coherence [32]. We employ an FDTD code developed in house, originally for modelling high-NA confocal microscopes [13, 28].

2.4. Image formation

After calculating the light scattered by the sample at plane S_sc in the sample space, it must be focussed onto the detector or coupled into a fiber. In the general case, where the NA of the objective may be high, a vectorial technique such as equivalent dipole decomposition, based on the Debye-Wolf integral [33], must be employed. This formalism allows aberrations in the optical detection system to be modelled. However, when a low NA objective is used, the low aperture approximations of Richards and Wolf [20] can be used to show that equivalent dipole decomposition reduces to scalar diffraction, where the two transverse field components may be considered separately. This allows Fourier optics [34] to be employed, whereby the angular spectrum of plane waves of the scattered field is propagated through the 4f system, shown in Fig. 1(c), to the detector. This results in the scattered field, E_sc(r_d,k), being known in the vicinity of the detector. The reference field, E_ref (r_d,k), may be similarly evaluated in the vicinity of the detector. Substituting the calculated values of E_sc(r_d,k) and E_ref (r_d,k) in to Eq. (1), and expanding, allows I_tot to be expressed as

3. Evaluation and analysis

This model has been previously verified for coherent optical microscopes [13]. This section, thus, focusses on how low-coherence interferometry is integrated into the model. We consider, first, the example of a stratified medium to verify the depth-sectioning performance of the model. We follow this with the simulation of OCT point spread functions over a range of axial positions. We conclude the section with a comparison of simulated and experimental OCT images of a three-dimensional structured phantom. The simulation parameters, such as source spectral width, FDTD cell size, etc., remain constant for each simulation.

3.1. Numerical dispersion

The phenomenon of numerical dispersion [17], well known to the FDTD community, can impact significantly upon the simulation of OCT image formation if it is not compensated. This form of dispersion has a different origin from that arising in real systems, which is not considered further here. Numerical dispersion refers to the phenomenon by which waves propagated in a homogeneous FDTD grid have wavelength-dependent phase and group velocities. Numerical dispersion occurs as a result of the solution of the difference equations derived from Maxwell’s equations diverging from the solution of the associated, continuous-case differential equation. In particular, in the continuous case, the phase refractive index n_p relates the optical frequency ν to the wavelength according to n_p = νλ/c, where c is the speed of light in vacuo. In the discrete case, it is necessary to define a numerical phase refractive index ñ_p and a numerical wavelength $\tilde{\lambda }$ which are related by ${\tilde{n}}_p=\nu \tilde{\lambda }/c$. Analogously to the continuous-case group refractive index n_g which satisfies n_g = n_p −λ dn_p/dλ [35], we define the discrete-case group refractive index ñ_g which satisfies

In practice, ñ_p is calculated by extracting $\tilde{\lambda }$ for a plane wave propagating along the z-axis in a homogeneous FDTD simulation and ñ_g is then calculated using Eq. (5). As an example, we calculated ñ_p and ñ_g for separate cases of media having refractive indices n₀ = n_p = n_g equal to 1 and 1.42. The results of these simulations are shown in Fig. 3 for a free space wavelength band centered on λ₀ = 1325nm. The FDTD simulation employed an isotropic cell size of λ₀/15. The group refractive index is seen to vary more than the phase refractive index in both cases and the deviation is worse in the case of n₀ = 1.42. Once calculated, these parameters can be used in more complex simulations. These results are important in the example following in Sec. 3.2, in which the OCT image of a stratified medium is simulated.

 

Fig. 3 Plots of numerical phase refractive index ñ_p and numerical group refractive index ñ_g for air (n₀ = 1) and dielectric material (n₀ = 1.42). The horizontal axis represents the continuous case wavelength.

Download Full Size | PDF

Numerical dispersion may be reduced by employing a smaller FDTD cell size, however, some additional compensation will generally be required. In the first instance, where simulations are composed of scatterers embedded in a background medium, the background medium ñ_p should be calculated to allow the field to be sampled at numerical wavenumbers matching the required continuous wavenumbers. This step needs to be made explicit only because complex amplitudes are extracted at specified frequencies, not wavenumbers, by application of a temporal discrete Fourier transform. As a second step, numerical dispersion within isolated scatterers can be compensated for by selecting a refractive index, n₀, which results in a mean value of ñ_p matching the desired refractive index. Finally, we note that FDTD simulations are performed for both the sample and reference arms so that both have the same numerical dispersion relationship.

3.2. Simulation of the OCT image of a stratified medium

The imaging model presented in this paper has been rigorously tested previously for the case of image formation using monochromatic light [13]. The objective of the current test is to verify the acquisition of axial information using low-coherence interferometry. We, thus, compared image formation for a stratified medium, using our model, with the analytic case. We employed a two-dimensional model, in order to reduce computational complexity, whereby the illumination and sample were considered to be uniform in the y-direction. We employed the transverse magnetic (TM) case of a previously developed two-dimensional formulation of the Debye-Wolf integral [36], although the model applies equally well to the transverse electric (TE) case. The FDTD simulation employed an isotropic cell size of λ₀/15 with center wavelength λ₀ = 1325nm. A wavelength band 200nm wide, centered on λ₀ was modelled. S(k) took the form of a Hanning window over the corresponding range of k. The FDTD computational grid had 6,000 cells in the transverse direction and 3,000 cells in the axial direction. The modelled stratified medium was composed of a background of refractive index 1 and sheets, one FDTD cell size thick, with refractive index 1.42. The light from the reference arm was also calculated using an FDTD simulation with a perfect electric conductor positioned at the focal plane of the focussed illumination. This ensures that the same numerical dispersion compensation outlined in Sec. 3.1 can be used for both the sample and reference arms, since the thin sheets contribute negligibly to the sample arm numerical dispersion relationship.

A cylindrical objective lens, of NA 0.056, and over-filled aperture was employed, with the additional lenses in Fig. 1(a) having the same NA for simplicity. The detector sensitivity function, D, was set equal to unity for each of the field components over the region of the detector in which the total field was non-negligible. Note that the longitudinal component of the electric field incident upon the detector was negligible due to the low NA of the detector lens (lens L3 in Fig. 1(a).

The A-scan calculated using the FDTD-based imaging model was compared with that calculated analytically using the characteristic matrix treatment of stratified media [37], using an angular spectrum of plane waves approach, as has previously been used as a platform for solving the inverse problem of OCT [38]. The reference beam is also calculated in this manner. Figure 4 shows two simulated A-scans, for 5 and 20 reflecting sheets in the stratified medium. The positions of the sheets are indicated by vertical lines. In the example with only 5 sheets, all sheets are resolvable and the axial resolution of the system is evident from the A-scan peak appearing at the lowest value of z. The simulation of 20 sheets shows examples of where not all sheets are resolvable, resulting in an axial speckle pattern. Both examples show how the proposed model predicts the OCT signal due to multiply reflected light, as is indicated by the signal appearing after the final sheet. The normalized mean squared error, calculated as ∑|OCT_n − OCT_a|²/∑|OCT_a|², where OCT_n and OCT_a are the complex A-scans calculated numerically and analytically, respectively, is 1.96×10⁻³ and 1.51 × 10⁻³ for the left and right plots, respectively. The level of agreement between the two models is, thus, sufficiently good to conclude that the low-coherence interferometry aspect of the model is correct.

 

Fig. 4 OCT A-scans for 5 (left) and 20 (right) thin sheets of material with refractive index 1.42, embedded in a background of refractive index 1. The locations of the sheets are indicated by vertical lines.

Download Full Size | PDF

3.3. Calculation of an OCT point spread function

Before demonstrating the simulation of practical image formation, we consider an uncomplicated example of system PSF simulation. The PSFs, as shown in Fig. 5, demonstrate how the imaging properties of a confocal microscope combine with low-coherence interferometry to yield an OCT PSF [39]. In particular, transverse resolution is obtained in the same way as a confocal microscope employing monochromatic light. Axial resolution is influenced by two phenomena. One is the depth sectioning possessed by microscopes employing monochromatic light, both confocal and conventional, and the other is coherence gating. Coherence gating dominates when a low numerical aperture objective is employed.

 

Fig. 5 Images and plots which demonstrate the depth-dependent PSFs employed in OCT by simulating images of 24 scatterers arranged equidistantly along the optical axis. (left) Each image corresponds to a different detector sensitivity function: D = (D_pt,0,0)^T, D = (D_int,0,0)^T and D = (D_ext,0,0)^T. (upper right) Transverse line plots corresponding to PSFs at three axial depths for each detector sensitivity function. (lower right) Axial line scans corresponding to the three detector sensitivity functions. The envelope function of the PSFs has been plotted to illustrate the confocal function for each case.

Download Full Size | PDF

A model identical to that of the stratified medium (Sec. 3.2) is employed, except with the FDTD grid extended to 6,000 cells deep and without the stratified medium. A series of 24 single FDTD cell scatterers, with refractive index 1.01, were arranged equidistantly along the optical axis, allowing depth-dependent PSFs to be calculated, as shown in Fig. 5. The PSFs were calculated by performing 120 separate FDTD simulations, for a series of transverse scatterer positions each differing by 0.52μm to ensure that the signals from different scatterers do not overlap. The TM and TE polarizations yielded practically identical results. A low refractive index contrast between scatterer and background ensured that each of the 24 PSFs shown in Fig. 5 were practically identical to that obtained by calculating each PSF separately.

Each of the three images displayed in the left portion of Fig. 5 correspond to different detector sensitivity functions: D = (D_pt,0,0), D = (D_int,0,0) and D = (D_ext,0,0). D_pt corresponds to a true point detector, D_int is unity over a width of 114μm and D_ext corresponds to an extended detector which detects nearly all light imaged into the detector space. The width of the support of D_int was chosen because it results in imaging properties approximately in between the point and extended cases. OCT systems generally employ a single-mode optical fibre to sample the reference and scattered light. In this case, the detector sensitivity function will take the form of the fibre mode. Although not shown, the simulated PSF for a single-mode optical fibre will be similar to the case of the point detector. The OCT system was focused by over-filling the aperture on the upper-most scatterer, which appears at the top of the PSFs shown in Fig. 5. The PSFs calculated using each of the three detector sensitivity functions have been normalized by the value of the corresponding PSF at the in-focus position. The images demonstrate how the PSF deteriorates with increased axial imaging distance from the focus. The transverse line scans demonstrate how, as expected, the PSF is narrower in the transverse direction when a pinhole detector is employed and increases as the detector radius increases. The non-Gaussian nature of the transverse line scans, which becomes more pronounced with axial distance from the focus, is due to an over-filled aperture being employed, although we note that, in practice, an under-filled aperture is generally employed.

The axial line scans have an envelope function fitted to the peak of each PSF. This envelope is the so-called confocal function, the plots of which show how the depth of focus can be increased by increasing the extent of the detector, albeit to the detriment of lateral resolution. We note, however, that in the case of the extended detector, the depth of focus is determined by geometrical effects rather than by a confocal pinhole. As is expected, aside from the amplitude modulation due to the confocal function, each of the individual axial PSFs has a profile matching the Fourier transform of S(k). Finally, as noted previously, the transverse profile of the PSFs has been verified previously in a model of confocal microscopes employing monochromatic light [13].

3.4. Simulation of OCT image formation for a phantom

We simulated image formation by using a numerical phantom based upon a physical three-dimensional structured phantom which has been developed in our laboratory [40]. The objective of this example is to demonstrate that the model is able to simulate the salient features of experimental OCT images, rather than to match a particular experimental result. OCT image formation for this phantom has also been recently modelled using a linear systems model of OCT, accounting only for single scattering phenomena [11]. The phantom is depicted schematically in Fig. 6(a). The sample feature consists of the letters “OBEL”, constructed by dispersing TiO₂ (diameter = 1μm, n=2.488 at λ₀ = 1300nm) particles into a silicone elastomer casting (n=1.42 at λ₀ = 1300nm) at a concentration of 10mg/mL. The scattering coefficient within the letters was calculated to be μ_s = 12.6mm⁻¹ at λ₀ = 1300nm using the theory of extinction by a slab of particles. The particle scattering cross section was calculated using Mie theory [42]. Experimental images of the physical phantom, presented in Fig. 6(b) and Fig. 8(c), were acquired using a Thorlabs Telesto II, spectral-domain OCT system possessing a source with central wavelength of 1300nm and axial resolution of 5.5μm, implying an effective 3db-bandwidth of 135nm. It has a sample arm objective with NA of 0.056, which is under-filled resulting in a lateral resolution (full-width at half-maximum) of 13μm. The experimentally acquired auto-correlation image is shown in Fig. 8(d) for comparison with the simulated images. The auto-correlation image results from the component of Eq. (1) due to purely scattered light, I_sc(k), and is usually separated from the meaningful OCT image by adjusting the position of the reference mirror relative to the sample, as is shown in Fig. 6(b). Although not generally of importance in dual-arm OCT systems, the autocorrelation image demonstrates the comprehensive nature of our model.

 

Fig. 6 a) Schematic diagram of the physical phantom reported in [40]. Imaging is performed from above the sample as demonstrated by the objective lens (not to scale). b) OCT image of the phantom acquired with a Telesto II system, showing the lettering feature at the centre and the auto-correlation image at the top. The nearly horizontal line depicts the top of the phantom’s embedding medium, from which light reflects, leading to the faintest “OBEL” image.

Download Full Size | PDF

A number of approximations are required to model, in two dimensions, a three-dimensional imaging system and sample. It is natural to represent spherical scatterers by cylinders in two dimensions as shown in Fig. 7. However, the radius of such cylinders will not, in general, be equal to that of the spheres they represent. This is because the geometrical cross section, perunit length, of a cylinder exceeds the geometrical cross section of a sphere with the same radius, at the length scale of interest. A further approximation arises when using the FDTD method to calculate light scattering, as cylindrical scatterers are represented approximately due to space being discretized with a resolution determined by the FDTD cell size. Thus, when approximating cylinders with radii approaching only a few FDTD cells, it is necessary to calculate the scattering cross section, per-unit length, of such scatterers using the FDTD method itself. One may, thus, vary the shape, refractive index and concentration of two-dimensional scatterers and illumination polarization to achieve the desired optical properties, such as scattering coefficient and anisotropy.

 

Fig. 7 Diagram showing how a collection of spherical scatterers in a three-dimensional beam (left) can be approximated in a two-dimensional system (right) where fields and scatterers extend uniformly in the y direction. The scatterers are, thus, cylinders, depicted as circles in the xz plane.

Download Full Size | PDF

The FDTD simulation was identical to that of Sec. 3.3 except for the scatterers. The illumination was focussed on the top of the letters in the phantom. Direct simulation of image formation, i.e., three-dimensional simulation, is currently not feasible due to computational limitations. The current simulation, thus, amounts to approximating image formation in an intrinsically three-dimensional system in two-dimensions. The task of matching the optical properties between two- and three-dimensional situations can be achieved through parametric studies of the effect of shape, refractive index and concentration of two-dimensional scatterers as described above. Our current interest is in providing a proof of principle demonstration of our model. To this end, we chose to represent the 1μm nominal diameter TiO₂ scatterers in the physical three-dimensional structured phantom with cylindrical scatterers of diameter 250nm and a refractive index matching that of TiO₂ (n=2.488 at λ = 1325nm). This enabled Mie theory to be used to estimate the scattering coefficient of the resulting arrangement of scatterers. A numerical phantom was then constructed by calculating a random arrangement of non-overlapping circles, in the xz plane, with diameter 250nm having centers r_n = (x_n,z_n) and density 5 × 10⁴mm⁻², where n is a circle index. Circles falling outside of the stencil of the letters “OBEL” were eliminated and the indices of scattering cells, within the FDTD simulation, were determined according to the list $\{(i,j)\in {\mathbb{Z}}^2:\sqrt{{(i\mathrm{\Delta }-x_n)}^2+(j\mathrm{\Delta }-z_n)}\le 125\mathrm{nm},\text{for any}n\}$ where Δ is the FDTD cell size.

The scattering coefficients, within the letters of “OBEL”, were 4.3 and 22 mm⁻¹ for the TM and TE cases, respectively. A good explanation of how to calculate scattering coefficients is given by Saidi et al [41] and Bohren and Huffman [42]. The scattering coefficient of an ensemble of particles may be calculated from the product of the volume density of particles and the scattering cross-section defined as C_sc = (1/I_i) ∫_S S_s · ndA, where I_i is the irradiance of the incident plane wave, S_s is the Poynting vector of the scattered field and S is a surface enclosing the particle with outward normal n. In the case of cylinders, the scattering cross-section, and thus scattering coefficient, are calculated on the per unit length basis as the cylinders are infinitely long. The scattering cross-section of a circle represented by FDTD grid cells was calculated for both the TE and TM illuminations using the FDTD method. In particular, S_s was calculated on a contour surrounding the scatterer, enabling the scattering cross-section to be evaluated. For both the case of discrete approximations to cylinders and ideal cylinders, theory shows that cylinders with a sub-wavelength diameter scatter light much more significantly when the electric field vector is parallel to the cylinder (i.e., the TE case).

The results of image formation using the numerical phantom are shown in Fig. 8. The scatterers in Fig. 8(a) are depicted with a larger diameter than was actually modelled for clarity. Figure 8b) shows the magnitude of the principal field components for one particular beam position, at the central wavelength, for the TM and TE cases, respectively. The difference in scattering cross sections, for the TM and TE cases, of each cylinder is evident from this plot. The field magnitude has been normalized by the magnitude of the field at the focus, in air, and in the absence of scatterers.

 

Fig. 8 a) Schematic diagram of the numerical phantom for which OCT image formation was simulated. The scatterers are not to scale; b) Magnitudes of the principal field component for the TM and TE cases, respectively, normalized by the in-focus free space field magnitude; c) Experimental OCT image of the physical phantom; d) Experimentally acquired auto-correlation image of the physical phantom; e) Simulated OCT image of the numerical phantom for the TM case; f) Simulated auto-correlation image of the numerical phantom for the TM case; g) Simulated OCT image of the numerical phantom for the TE case; h) Simulated auto-correlation image of the numerical phantom for the TE case. The scale bars in each image indicate 50μm. The axial scale represents physical distance, scaled using the refractive index of silicone (1.42) for the experimental images and the mean group numerical refractive index for the simulated cases (1.48, see Fig. 3).

Download Full Size | PDF

The experimentally acquired OCT and auto-correlation images are shown in Fig. 8(c) and Fig. 8(d), respectively. They are displayed on a scale covering 0 to 37dB of OCT signal-to-noise ratio. Figure 8(e) and Fig. 8(g) show the simulated OCT images for the TM and TE cases, respectively. The simulated images are displayed on a logarithmic scale such that the peak OCT signal in the simulated image corresponds to the peak value of OCT signal-to-noise ratio in the experimental images.

The scattering coefficient within the letters for the experimental image (μ_s = 12.6mm⁻¹) lies between that of the two simulated images (μ_s=4.3 and 22 mm⁻¹, respectively). As a result, multiple scattering is just visible below the letters in Fig. 8(c), absent in Fig. 8(e) and prominent in Fig. 8(g). Furthermore, attenuation of the OCT signal is observable in Fig. 8(c) in regions near the base of each letter, but is much more pronounced in Fig. 8(g). The simulated images, thus, represent extreme cases of low and high scattering and the experimental image represents an intermediate case.

We have presented the auto-correlation images to show the model’s inherent ability to model contributions to image formation from sources other than those which are usually assumed to make the principal contribution to OCT images. Furthermore, in other experiments we employ a common-path configuration in our laboratory which uses the interface between a glass slide and the sample to achieve high phase sensitivity [43]. In this case, the auto-correlation image is included in the acquired image of the superficial region of the sample. Simulation of the auto-correlation image will, thus, provide the opportunity to explore the impact that the auto-correlation image has on phase sensitivity.

In contrast to the auto-correlation image, the magnitude of the OCT signal depends on the reflectivity of the reference mirror. We have, thus, displayed the simulated OCT and auto-correlation images on the same scale, noting that the effect of different reference mirror reflectivities can be modelled simply by offsetting the logarithmic scale of the OCT image, in accordance with the reflectivity of the reference mirror. Finally, we note that it is also possible to model the degradation of sensitivity due to parasitic reflections in the optical system simply by building these into the simulation of the optical system. It is also possible to model shot noise and sensitivity roll-off, since we calculate the irradiance at the detector surface.

We conclude this section by considering the computational burden involved with simulating the images in Fig. 8. The OCT B-scans in Fig. 8 are all derived from the same set of 220 FDTD simulations, representing adjacent A-scans separated by 2.2μm, requiring 38 minutes per A-scan to complete using all 6 cores of an Intel Xeon E5645 CPU (2.4 GHz clock speed, 12 Mb cache) and required just less than 4 Gb of RAM.

4. Discussion and conclusions

We have presented a framework, and an example implementation, of a rigorous, full wave model of OCT. At the core of the model are rigorous, time-domain electromagnetic scattering solvers, such as the FDTD method, which allow the light scattered by the deterministic sample to be calculated in a way which automatically includes phenomena such as multiple scattering and vectorial effects.

A range of phenomena arising from OCT imaging can now be examined. For example, in the future, we plan to consider applications such as displacement measurement using phase sensitive detection [11], parametric imaging [46], the use of non-Gaussian beams, such as Bessel beams [47, 48], and to test hypotheses regarding unresolved features observed in a variety of medical and biomedical OCT images [15].

The model itself has scope for improvement. Most notably, if applications require it, the modelling of dispersive media in the FDTD method can be implemented as can anisotropic media. Sources of noise can also be incorporated. In fact, one of the prime strengths of our model is that its modular nature allows phenomena arising in, for example, a sample, optical elements and detector, to be modelled by modifying only one component of the model.

The capability of modern desktop (and superior) computers, along with the emergence of open source FDTD implementations [44, 45] mean that this kind of simulation will eventually be accessible to non-specialists. Access to institutional computer clusters will enable volume scans to be evaluated in on the order of a day, which is a time short enough to be of practical use.