A full wave model of image formation in optical coherence tomography applicable to general samples

We demonstrate a highly realistic model of optical coherence tomography, based on an existing model of coherent optical microscopes, which employs a full wave description of light. A defining feature of the model is the decoupling of the key functions of an optical coherence tomography system: sample illumination, light-sample interaction and the collection of light scattered by the sample. We show how such a model can be implemented using the finite-difference time-domain method to model light propagation in general samples. The model employs vectorial focussing theory to represent the optical system and, thus, incorporates general illumination beam types and detection optics. To demonstrate its versatility, we model image formation of a stratified medium, a numerical point-spread function phantom and a numerical phantom, based upon a physical three-dimensional structured phantom employed in our laboratory. We show that simulated images compare well with experimental images of a three-dimensional structured phantom. Such a model provides a powerful means to advance all aspects of optical coherence tomography imaging. © 2014 Optical Society of America OCIS codes: (170.3660) Light propagation in tissues; (050.1755) Computational electromagnetic methods; (180.0180) Microscopy; (110.4500) Optical coherence tomography. References and links 1. B. J. Davis, S. C. Schlachter, D. L. Marks, T. S. Ralston, S. A. Boppart, and P. S. Carney, “Nonparaxial vectorfield modeling of optical coherence tomography and interferometric synthetic aperture microscopy,” J. Opt. Soc. Am. A 24, 2527–2542 (2007). 2. M. Kirillin, I. Meglinski, V. Kuzmin, E. Sergeeva, and R. Myllylä, “Simulation of optical coherence tomography images by Monte Carlo modeling based on polarization vector approach,” Opt. Lett. 18, 21714–21724 (2010). 3. I. Meglinski, M. Kirillin, V. Kuzmin, and R. Myllyla, “Simulation of polarization-sensitive optical coherence tomography images by a Monte Carlo method,” Opt. Lett. 33, 1581–1583 (2008). 4. J. M. Schmitt and A. Knüttel, “Model of optical coherence tomography of heterogeneous tissue,” J. Opt. Soc. Am. A 14, 1231–1242 (1997). 5. I. V. Turchin, E. A. Sergeeva, L. S. Dolin, V. A. Kamensky, N. M. Shakhova, and R. Richards-Kortum, “Novel algorithm of processing optical coherence tomography images for differentiation of biological tissue pathologies,” J. Biomed. Opt. 10, 064024 (2005). 6. T. S. Ralston, D. L. Marks, P. S. Carney, and S. A. Boppart, “Inverse scattering for optical coherence tomography,” J. Opt. Soc. Am. A 23, 1027–1037 (2006). 7. A. Tycho, T. M. Jorgensen, H. T. Yura, and P. E. Andersen, “Derivation of a Monte Carlo method for modeling heterodyne detection in optical coherence tomography systems,” Appl. Opt. 41, 6676–6691 (2002). 8. Q. Lu, X. Gan, M. Gu, and Q. Luo, “Monte Carlo modeling of optical coherence tomography imaging through turbid media,” Appl. Opt. 43, 1628–1637 (2004). 9. L. S. Dolin, “A theory of optical coherence tomography,” Radiophys. Quantum El. 41, 850–873 (1998). 10. L. H. Wang, S. L. Jacques, and L. Q. Zheng, “MCML Monte-Carlo modeling of light transport in multilayered tissues,” Comput. Meth. Prog. Bio. 47, 131–146 (1995). 11. L. Chin, A. Curatolo, B. F. Kennedy, B. J. Doyle, P. R. T. Munro, R. A. McLaughlin, and D. D. Sampson, “Analysis of image formation in optical coherence elastography using a multiphysics approach,” Biomed. Opt. Express 5, 2913–2930 (2014). 12. J. A. Izatt, E. A. Swanson, J. G. Fujimoto, M. R. Hee, and G. M. Owen, “Optical coherence microscopy in scattering media,” Opt. Lett. 19, 590–592 (1994). 13. P. Török, P. R. T. Munro, and E. E. Kriezis, “High numerical aperture vectorial imaging in coherent optical microscopes,” Opt. Express 16, 507–523 (2008). 14. I. R. Capoglu, J. D. Rogers, A. Taflove, and V. Backman, “Chapter – 1: The microscope in a computer: image synthesis from three-dimensional full-vector solutions of Maxwell’s equations at the nanometer scale,” in Progress in Optics, E. Wolf, ed. (Elsevier, 2012), 57, pp. 1–91. 15. D. C. Reed and C. A. DiMarzio, “Computational model of OCT in lung tissue,” Proc. SPIE 7570, 75700I (2010). 16. Y.-T. Hung, S.-L. Huang, and S. H. Tseng, “Full EM wave simulation on optical coherence tomography: impact of surface roughness,” Proc. SPIE 8592, 859216 (2013). 17. A. Taflove and S. Hagness, Computational Electrodynamics, Third Edition (Artech House, 2005). 18. L. Mandel and E. Wolf, Optical Coherence and Quantum Optics (Cambridge University, 1995). 19. P. D. Higdon, P. Török, and T. Wilson, “Imaging properties of high aperture multiphoton fluorescence scanning microscopes,” J. Mod. Opt. 193, 127–141 (1998). 20. B. Richards and E. Wolf, “Electromagnetic diffraction in optical systems ii. structure of the image field in an aplanatic system,” Proc. Roy. Soc. A 253, 358–379 (1959). 21. V. S. Ignatowsky, “Diffraction by a lens of arbitrary aperture,” T. Opt. Inst. Petrograd 1, 1–36 (1919). 22. R. Luneburg, Mathematical Theory of Optics (University of California, Berkeley and Los Angeles, 1966). 23. K. S. Youngworth and T. G. Brown, “Focusing of high numerical aperture cylindrical-vector beams,” Opt. Express 7, 77–87 (2000). 24. P. Török and P. R. T. Munro, “The use of Gauss-Laguerre vector beams in STED microscopy,” Opt. Express 12, 3605–3617 (2004). 25. P. Török and P. Varga, “Electromagnetic diffraction of light focused through a stratified medium,” Appl. Opt. 36, 2305–2312 (1997). 26. P. Török, “An imaging theory for advanced, high numerical aperture optical microscopes,” D.Sc. Thesis (2004). 27. J.-P. Berenger, “A perfectly matched layer for the absorption of electromagnetic waves,” J. Comput. Phys. 114, 185–200 (1994). 28. P. R. T. Munro and P. Török, “Vectorial, high numerical aperture study of Nomarski’s differential interference contrast microscope,” Opt. Express 13, 6833–6847 (2005). 29. P. R. T. Munro, D. Engelke, and D. D. Sampson, “A compact source condition for modelling focused fields using the pseudospectral time-domain method,” Opt. Express 22, 5599–5613 (2014). 30. A. Oskooi and S. G. Johnson, “Electromagnetic wave source conditions,” in Advances in FDTD Computational Electrodynamics: Photonics and Nanotechnology, A. Taflove, A. Oskooi and S. G. Johnson, eds. (Artech House, 2013), pp. 65–96. 31. D. E. Merewether, R. Fisher, and F. W. Smith, “On implementing a numeric Huygen’s source scheme in a finite difference program to illuminate scattering bodies,” IEEE T. Nucl. Sci. 27, 1829–1833 (1980). 32. I. R. Capoglu, A. Taflove, and V. Backman, “FDTD simulation of a partially-coherent Gaussian Schell-model beam,” in IEEE International Symposium on Antennas and Propagation (IEEE, 2011), pp. 2286–2288. 33. P. R. T. Munro and P. Török, “Calculation of the image of an arbitrary vectorial electromagnetic field,” Opt. Express 15, 9293–9307 (2007). 34. J. W. Goodman, Introduction to Fourier Optics (McGraw-Hill, New York, 1988). 35. J. R. Rogers and M. D. Hopler, “Conversion of group refractive index to phase refractive index,” J. Opt. Soc. Am. A 5, 1595–1600 (1988). 36. C. J. R. Sheppard, “Cylindrical lenses–focusing and imaging: a review [invited],” Appl. Opt. 52, 538–545 (2013). 37. M. Born and E. Wolf, Principles of Optics, Seventh Edition (Cambridge University, 1999). 38. O. Bruno and J. Chaubell, “One-dimensional inverse scattering problem for optical coherence tomography,” Inverse Probl. 21, 499524 (2005). 39. J. A. Izatt and M. A. Choma, “Theory of optical coherence tomography,” in Optical Coherence Tomography, Theory and Applications, W. Drexler and J. G. Fujimoto, eds. (Springer, 2008), pp. 47–72. 40. A. Curatolo, B. F. Kennedy, and D. D. Sampson, “Structured three-dimensional optical phantom for optical coherence tomography,” Opt. Express 19, 19480–19485 (2011). 41. I. S. Saidi, S. L. Jacques and F. K. Tittel, “Mie and rayleigh modeling of visible-light scattering in neonatal skin”, Appl. Opt. 34, 7410–7418 (1995). 42. C. F. Bohren and D. R. Huffman, Absorption and Scattering of Light by Small Particles (Wiley Interscience, 1983). 43. B. F. Kennedy, R. A. McLaughlin, K. M. Kennedy, L. Chin, A. Curatolo, A. Tien, B. Latham, C. M. Saunders and D. D. Sampson, “Optical coherence micro-elastography: mechanical-contrast imaging of tissue microstructure,” Biomed. Opt. Express 5, 2113–2124 (2014). 44. A. F. Oskooi, D. Roundy, M. Ibanescu, P. Bermel, J. D. Joannopoulos, and S. G. Johnson, “MEEP: A flexible free-software package for electromagnetic simulations by the FDTD method,” Comput. Phys. Commun. 181, 687–702 (2010). 45. P. Wahl, D.-S. Ly-Gagnon, C. Debaes, D. Miller, and H. Thienpont, “B-calm: An open-source GPU-based 3DFDTD with multi-pole dispersion for plasmonics,” Opt. Quant. Electron. 44, 285–290 (2012). 46. P. Gong, R. A. McLaughlin, Y. Liew, P. R. T. Munro, F. M. Wood, and D. D. Sampson, “Assessment of human burn scars with optical coherence tomography by imaging the attenuation coefficient of tissue after vascular masking,” J. Biomed. Opt. 19, 021111 (2013). 47. Z. Ding, H. Ren, Y. Zhao, J. S. Nelson, and Z. Chen, “High-resolution optical coherence tomography over a large depth range with an axicon lens,” Opt. Lett. 27, 243–245 (2002). 48. D. Lorenser, C. C. Singe, A. Curatolo and D. D. Sampson, “Energy-efficient low-Fresnel-number Bessel beams and their application in optical coherence tomography,” Opt. Lett. 39, 548–551 (2014).


Introduction
A variety of methods for simulating optical coherence tomography (OCT) image formation have been developed [1][2][3][4][5][6][7][8][9][10][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][2][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 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 imageformation 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 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, are 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.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].This will 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 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.

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: 1. Analytic calculation of the light emitted from the source, propagated through the optical system and incident upon the sample, E inc (r s , k); 2. Numerical calculation of the light scattered back by the sample to the objective lens, E sc (r s , k); 3. Imaging of the light scattered back by the sample onto the detector, E sc (r d , k); 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 re f (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.The spectrally resolved detector signal may then be evaluated according to: where is the effective spectral intensity distribution of the source, combining the spectrally dependent power of the source, any source spectral shaping and detector characteristics, diag(D(r d )) is the diagonal matrix with diagonal entries taken from the vector D(r d ), where D(r d ) is the detector sensitivity function, such as the mode of a single-mode fiber, D is the spatial extent of the detector and the norm is defined by |a| 2 = a * a, where * represents complex conjugate transpose.Note that although we represent the effective spectral intensity distribution as a scalar function, there is no reason why it cannot be represented by a tensor which operates on field quantities.Furthermore, different spectral intensity distributions can be used for the reference and sample arms if required, particularly as the two are calculated using independent FDTD simulations.The OCT A-scan may then be found by inverse Fourier transform of I tot (k).We explain how Steps 1-3 above are performed in the remainder of this section.

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][21][22], which remains valid for low NA lenses.It has been applied in a wide variety of applications, such as 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: where, as shown in Fig. 1(b), s = (s x , s y , s z ) is the unit vector of a typical ray which passes through the nominal focus of the objective, NA is the numerical aperture of the lens in air, A i (s) represents possible phase (including aberrations) and amplitude modulation, [L i j ] accounts for refraction by the lens and is evaluated using the generalized Jones matrix formalism [19], E 0 (s, k) is the geometrical optics approximation to the electric field vector incident upon the back focal plane of the lens and f is the lens focal length.A similar expression exists for the magnetic field.

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 λ 0 /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 λ 0 /2 may be approached; however, our recent work [29] shows that a discretization not exceeding λ 0 /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 nonphysical 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]: where k is the unit vector aligned with the z-axis, E inc (r s , k 0 ) is the incident field evaluated using Eq. ( 2) at mean wavenumber, k 0 , ω 0 = c/k 0 , where c is the speed of light, t 0 is set such that the components of J * s (t) are vanishingly small at t = 0 and ℜ denotes the real part.W controls the temporal width of the incident pulse and, thus, the spectral width of the source, allowing the scattered field at wavelengths throughout the spectrum of interest to be calculated in a single simulation, making this model computationally tractable.
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].

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 re f (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 re f (r d , k) in to Eq. ( 1), and expanding, allows I tot to be expressed as where 2 is due solely to the reference arm light, 2 is due solely to the scattered light and )dS} is due to the interference of the scattered and reference arm light.Again we note that the use of S(k) does not imply that the sample and reference arms have the same spectral dependence.Any differing spectral dependence can be accounted for in the independent calculations of each field quantity.Inverse Fourier transformation of I tot with respect to k results in an OCT A-scan.After performing this transformation, I re f results in a "DC" contribution, I sc leads to the auto-correlation terms and I int contributes the cross-correlation terms, which are generally considered to constitute the OCT image.In practice, E re f (r d , k) and E sc (r d , k) may be calculated at only a finite, but arbitrary, number of values of k and so a discrete inverse Fourier transform is employed.The main task in constructing a full wave model of OCT is, thus, the calculation of E re f (r d , k) and E sc (r d , k) over a sufficient range of wavenumbers, using the FDTD method as outlined in Sec.2.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 at 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.

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, λ which are related by ñp = ν λ /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 λ 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 0 = 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 λ 0 = 1325nm.The FDTD simulation employed an isotropic cell size of λ 0 /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 0 = 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.
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 0 , 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.

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 ñ n ñp /n 0 , n 0 = 1 ñg /n 0 , n 0 = 1 ñp /n 0 , n 0 = 1.42 ñg /n 0 , n 0 = 1.42 Fig. 3. Plots of numerical phase refractive index ñp and numerical group refractive index ñg for air (n 0 = 1) and dielectric material (n 0 = 1.42).The horizontal axis represents the continuous case wavelength.
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 λ 0 /15 with center wavelength λ 0 = 1325nm.A wavelength band 200nm wide, centered on λ 0 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 where OCT n and OCT a are the complex A-scans calculated numerically and analytically, respectively, are 1.96×10 −3 and 1.51 × 10 −3 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.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.

Calculation of an OCT point spread function
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.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].

Simulation of OCT image formation for a phantom
We simulated image formation using a numerical phantom based upon a physical threedimensional 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 schemat-ically in Fig. 6(a).The sample feature consists of the letters "OBEL", constructed by dispersing TiO 2 (diameter = 1μm, n=2.488 at λ 0 = 1300nm) particles into a silicone elastomer casting (n=1.42 at λ 0 = 1300nm) at a concentration of 10mg/mL.The scattering coefficient within the letters was calculated to be μ s = 12.6mm −1 at λ 0 = 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  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.The FDTD simulation was identical to that of Sec.3.3 except for the scatterers.One point to note is that the illumination was focussed on the top of the letters in the phantom.Direct sim-ulation 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 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 2 scatterers in the physical three-dimensional structured phantom with cylindrical scatterers of diameter 250nm and a refractive index matching that of TiO 2 (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 nonoverlapping circles, in the xz plane, with diameter 250nm having centers r n = (x n , z n ) and density 5 × 10 4 mm −2 , 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) ∈ Z 2 : (iΔ − x n ) 2 + ( jΔ − z n ) ≤ 125nm, for any n} where Δ is the FDTD cell size.
The scattering coefficients, within the letters of "OBEL", were 4.3 and 22 mm −1 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.
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-tonoise 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 −1 ) lies between that of the two simulated images (μ s =4.3 and 22 mm −1 , 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 autocorrelation 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 autocorrelation 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 Ascan 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.

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 FDTD can be implemented as can anisotropic media.Sources of noise can also be incorporated.In fact, one of the prime strengths of this method is that the modular nature of the model 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.

Fig. 1 .
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, E 0 (s, k).Phase and amplitude modulation are described by A i (s) whilst refraction by L2 is described by L i j ; 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).

Fig. 2 .
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.

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

Fig. 5 .
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 functions 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.
(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 3dbbandwidth 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. 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 .Fig. 7 .
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) and b) an 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.

Fig. 8 .
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).

3
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).
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).