Exploiting data redundancy in computational optical imaging

We present an algorithm which exploits data redundancy to make computational, coherent, optical imaging more computationally efficient. This algorithm specifically addresses the computation of how light scattered by a sample is collected and coherently detected. It is of greatest benefit in the simulation of broadband optical systems employing coherent detection, such as optical coherence tomography. Although also amenable to time-harmonic data, the algorithm is designed to be embedded within time-domain electromagnetic scattering simulators such as the psuedo-spectral and finite-difference time domain methods. We derive the algorithm in detail as well as criteria which ensure accurate execution of the algorithm. We present simulations that verify the developed algorithm and demonstrate its utility. We expect this algorithm to be important to future developments in computational imaging. © 2015 Optical Society of America OCIS codes: (050.1755) Computational electromagnetic methods; (180.0180) Microscopy; (110.4500) Optical coherence tomography; (050.1960) Diffraction theory; (110.1758) Computational imaging. References and links 1. P. Elbau, L. Mindrinos, and O. Scherzer, “Mathematical methods of optical coherence tomography,” in Handbook of Mathematical Methods in Imaging, O. Scherzer, ed. (Springer, 2015). 2. T. Brenner, D. Reitzle, and A. Kienle, “An algorithm for simulating image formation in optical coherence tomography for cylinder scattering,” in “European Conferences on Biomedical Optics,” Proc. SPIE 9541, 95411F (2015) 3. P. R. T. Munro, A. Curatolo, and D. D. Sampson, “Full wave model of image formation in optical coherence tomography applicable to general samples,” Opt. Express 23, 2541–2556 (2015). 4. 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). 5. 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). 6. I. R. Capoglu, J. D. Rogers, A. Taflove, and V. Backman, “Computational optical imaging using the finitedifference time-domain method,” in Advances in FDTD Computational Electrodynamics: Photonics and Nanotechnology, A. Taflove, A. Oskooi and S. G. Johnson, ed. (Artech House, 2013). 7. R. L. Coe and E. J. Seibel, “Computational modeling of optical projection tomographic microscopy using the finite difference time domain method,” J. Opt. Soc. Am. A 29, 2696–2707 (2012). 8. 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). 9. V. S. Ignatowsky, “Diffraction by a lens of arbitrary aperture,” T. Opt. Inst. Petrograd 1, 1–36 (1919). 10. R. Luneburg, Mathematical Theory of Optics (University of California Press, 1966). 11. K. Yee, “Numerical solution of inital boundary value problems involving maxwell’s equations in isotropic media,” IEEE T. Antenn. Propag. 14, 302–307 (1966). 12. A. Taflove and S. Hagness, Computational Electrodynamics, Third Edition (Artech House, 2005). 13. P. Török, P. R. T. Munro, and E. E. Kriezis, “Rigorous nearto far-field transformation for vectorial diffraction calculations and its numerical implementation,” J. Opt. Soc. Am. A 23, 713–722 (2006). 14. 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). 15. J. Goodman, Introduction to Fourier Optics (McGraw-Hill, 1988). 16. 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). 17. K. Matsushima and T. Shimobaba, “Band-limited angular spectrum method for numerical simulation of freespace propagation in far and near fields,” Opt. Express 17, 19662–19673 (2009). 18. 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, 55995613 (2014). 19. 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). 20. M. Born and E. Wolf, Principles of Optics (Cambridge University Press, 1999), Seventh Edition. 21. M. Frigo and S. Johnson, “The design and implementation of FFTW3,” Proc. IEEE 93, 216–231 (2005). 22. V. Čı́žek, Discrete Fourier transforms and their applications (Adam Hilger, 1986).


Introduction
There is currently a growing interest in rigorous simulation of broadband, coherent optical imaging systems, particularly those aimed at resolution in the tens of microns, such as optical coherence tomography (OCT) [1][2][3].This is being enabled by advancements in both algorithms and computer hardware which, in combination, make the simulation of image formation for practically significant tissue volumes feasible.Broadband simulation adds an entire dimension to the simulation parameter space, when compared with monochromatic simulation.This additional parameter, wavelength, poses a new computational bottleneck for already computationally burdensome models.We have developed a method to remove this bottleneck by exploiting data redundancy.This development leads to a reduction in both the computational complexity and computer memory needs of our model.This reduced computer memory requirement is particularly important as execution of models on memory-limited graphical processing units (GPU) is necessary for computational optical imaging (COI) to realize its potential.
COI of three-dimensional samples is well established [4][5][6][7].These models are amenable to broadband simulation yet, with the exception of [6], exhibit monochromatic examples, principally due to the computational bottleneck of broadband simulation.Rigorous COI is performed by decoupling the four principal components of a coherent optical microscope [5].The first component of the model is that of the focussed illumination, which uses the established Debye-Wolf vectorial focussing formalism [8][9][10].The second component of the model uses the finite-difference time-domain (FDTD) method [11,12] to calculate how the focussed beam is scattered by an arbitrary sample.The third component rigorously re-samples the scattered field, calculated within a relatively small, densely sampled grid, onto a larger grid [13].The reason for this is that the fourth component, a numerical method for propagating the scattered field through the optical system [14], requires the field to be sampled on a plane large enough to have a negligible field at its boundary.
The greatest computational burden of monochromatic COI is the FDTD method which calculates light scattered by the sample.The scattered field is recorded for a single wavelength, only on the surface of the FDTD method's computational grid, thus contributing negligibly to the simulation's total computer memory requirement.The second greatest computational burden of the simulation is re-sampling the scattered field on a plane larger than the transverse extent of the FDTD's computational grid.Each sample point on this larger plane results from an integration over the entire surface of the FDTD grid, a considerable computational burden.This burden becomes even more significant if it must be performed for a range of wavelengths, potentially exceeding the computational burden of the FDTD method.Furthermore, storing the scattered field for a range of wavelengths on the surface of the FDTD grid leads to a significant increase in its computer memory consumption, which may be prohibitive on some computer architectures such as GPUs.
The computational bottleneck of broadband simulation can, however, be circumvented when performing low NA, coherent, COI.The first step towards this comes by recognising that depolarisation due to focussing can be neglected at low NA (say NA< 0.5).The second step comes by noting that there is a large degree of redundancy in the data contained within a COI simulation.In particular, we do not need to know the spatial distribution of scattered light in the space of the detector, only the integral of a coherent detector sensivitiy function multiplied by the scattered field.We have developed an algorithm which exploits this data redundancy to eliminate one component of COI models.This algorithm leads to gains in both computational and computer memory efficiency in broadband COI.
The paper begins by introducing the optical system that this algorithm applies to.We then derive the new algorithm and derive criteria which must be observed for this algorithm to produce accurate results.We conclude the paper with examples which demonstrate how this algorithm can be applied to the simulation of image formation.

Background to algorithm
The optical system under study is depicted in Fig. 1 and is representative of a typical optical system found in a scanning coherent optical imaging system such an optical coherence tomography (OCT) system.In particular, incident light is delivered by an optical fiber and focussed onto the sample by a system well approximated by a 4f system.Incident light is scattered within the sample space, which is subsequently focussed and coupled in to the optical fiber.We note that the system in Fig. 1 is an idealisation as practical systems require space for elements such as scanning mirrors, however, such systems are still represented to a good approximation by a 4f system.
As discussed in the introduction, the propagation of scattered light to the optical fiber has previously been simulated using a two-stage computation [5] when high numerical aperture focussing is employed.The use of low numerical aperture focussing, however, alleviates the need for a two step computation as depolarisation due to focussing can be neglected.Indeed, we have previously demonstrated how image formation for the system shown in Fig. 1 can be rigorously simulated for arbitrary scattering configurations using a single stage of computation, but within a two-dimensional approximation [3].In particular, the lens and sample were assumed to extend uniformly in one of the transverse directions.The main reason for the twodimensional restriction was to make the model computationally feasible.The model works by using the finite-difference time-domain (FDTD) method to model light propagation in the sample.The numerically calculated scattered light is recorded at the surface of the FDTD grid and, at the conclusion of the FDTD simulation, the detection of this light is calculated by a further Fourier-optics based algorithm.As we are now developing a three-dimensional implementation of this model it is necessary to develop methods which reduce the computational burden of the two-dimensional model.An important development is the use of the psuedospectral time-domain (PSTD) method in place of the FDTD method, allowing for a significantly reduced spatial sampling requirement in the sample space.However, when modelling broadband imaging systems such as OCT, the scattered field must be recorded on a transverse plane of the PSTD grid for a range of wavelengths.In the case of OCT, this means that the computer memory required to store the scattered field can equal, or even exceed, the amount of memory required to perform the PSTD simulation.As previously mentioned, there is redundancy in this data since we require only the integral of a coherent detector sensitivity function multiplied by the scattered field.We have thus developed an algorithm which is embedded in the PSTD simulation which calculates how scattered light is coupled back into an optical fiber, or indeed any coherent detector, as a PSTD simulation progresses, thus significantly reducing the computer storage required to perform rigorous coherent COI.We would like to point out that the presented algorithm works equally well embedded in a FDTD simulation.

Derivation of algorithm
Most of the derivations and results in this paper are for scalar fields.We note however, that vectorial simulations can be performed by considering the orthogonal field components separately.This is possible because this work is limited to low numerical aperture (NA), say less than 0.5, objectives where focussing can be accurately represented within the scalar approximation.Here, NA, refers to the right hand lens in Fig. 1 as this is the lens which faces the sample.Due to the sine condition, the NA of the left hand lens must be equal to ( f 2 / f 1 )NA.We have defined three distinct spaces in the optical system of Fig. 1, in particular, the optical fiber space, the focal plane common to both lenses and the sample space, to which we apply indices of 1, 2 and 3, respectively.We use the symbols r, q, U and u to represent spatial coordinates, spatial-frequency coordinates, complex amplitudes and time-domain fields, respectively, within a particular space.We associate each parameter and complex amplitude with a particular space by applying the index of the applicable space as a subscript.Finally, we denote the Fourier transform of a field by applying a tilde (e.g., U).
We develop our algorithm using a time-harmonic formalism.The conventional approach would be to apply such an algorithm to the final result of a PSTD simulation, thus on timeharmonic data.Time-harmonic fields are extracted from the intrinsically time-domain fields simulated in the PSTD simulation by performing temporal discrete Fourier transforms of the time-domain fields at the required spatial and temporal-frequency sample points.In particular, time-harmonic fields U 3 (r 3 , f ) can be extracted at observation points r 3 and temporal frequencies f , from the time-domain field u 3 (r 3 , n∆t) according to where the index n refers the time step index of the PSTD simulation, N t is the total number of time steps in the PSTD simulation and ∆t is the PSTD time step.We note that, in general, ∆t is smaller than that required to correctly sample the fields, and so in practice, Eq. (1) may be evaluated at every Mth time steps so long as One way [3] of calculating how the scattered field is coupled in to the optical fiber is to take the field U 3 (r 3 , f ) resulting from the PSTD simulation and use Fourier optics to propagate the field through the optical system.This amounts to performing discrete Fourier transforms with respect to the spatial variable r 3 , applying a spatial-frequency domain filter, followed by a further discrete Fourier transform, giving the field in the fiber space as: where F represents that discrete Fourier transform operator with respect to spatial coordinates r 3 and r 2 = q 3 f 2 λ , where q 3 is the spatial-frequency parameter corresponding to r 3 .This shows that we may perform the Fourier optics transformation of the field from the sample to fiber space concurrently with the PSTD simulation.
Having explained that the transformation from sample space to detector space is performed concurrently with the PSTD simulation, we now derive this transformation.We will do this by supposing that we indeed have access to the time-harmonic field as is permitted by the equivalence of Eq. ( 3) and Eq. ( 4).The field scattered by the sample is observed at the plane z 3 = z obs and is denoted U 3 (r 3 , z obs ) where r 3 is a two-dimensional vector denoting a position in the plane of observation.We allow for the observation plane to be located within a medium of refractive index n s with an air interface at z 3 = −h.Our objective is to evaluate where φ 1 (r 1 ) represents a detector sensitivity function such as a fiber mode and U 1 (r 1 ) is the complex amplitude of the scattered field in the detector space.As explained previously, we use the notation of U(q) and U(r) to denote Fourier transform pairs.Using an angular spectrum approach, it is straightforward to show that [15]: where and P(q) is the pupil function defined as We can then find where ⊗ represents convolution.Substituting Eq. ( 9) into Eq.( 5) then gives where the second line follows from the definition of Fourier transforms and the third line from application of the convolution theorem.F {U(r)} denotes Fourier transform which we define as: Then, making use of Plancherel's theorem which states that for square integrable functions g(x) and and neglecting the negative sign, we then can write Eq. ( 10) as We note that Eqs. ( 13) and ( 14) are the main results of this paper.In particular, Eq. ( 13) shows how the field coupled in to the fiber may be evaluated within the central focal plane of the 4f system and Eq. ( 13) shows how this quantity can be calculated in the sample space.We also note that a similar result to Eq. ( 14) has been derived by Ralston et.al [16], however their result was applicable only to Gaussian detector sensitivity functions.As these equations will be applied to sampled data, the integrals must be replaced by summations and the Fourier transforms by discrete Fourier transforms.We address some of the limitations imposed by this in Sec. 4. sample space detector space optical fiber Fig. 1.Schematic diagram of the optical system under study.Each position vector r 1 , r 2 and r 3 is a two dimensional vector representing transverse displacement relative to a location on the optical axis.The field is observed on plane z 3 = z obs and there is an infinite half-space of material with refractive index n s beginning at z 3 = −h.z 3 is measured relative to the nominal focus of the right hand lens, in air.Angles θ 1 and θ 2 represent the semi-convergence angles of typical rays.R a is the radius of the aperture.

Translation of coherent detector
We note that in some instances it is useful to be able to simulate the translation of the coherent detector, from the on-axis position to an arbitrary in-focus position r d .This means modifying Eq. ( 5) to: This is easily included in Eqs. ( 13) and ( 14) by multiplying φ 1 (( . This is useful for two principal reasons.The first reason is that it allows for efficient simulation of image formation in wide-field coherent optical microscopes.
The second reason is that it suggests a formalism by which a scattered field can be sensed in a compressed manor.
3.2.Implementation of Eqs. ( 13) and ( 14) We now need to show how Eq. ( 13) and Eq. ( 14) are incorporated into the PSTD simulation.As previously explained, the PSTD simulation gives access to time-domain fields u 3 (r 3 , n∆t) at iteration n of the PSTD simulation.∆t is required to satisfy the Courant stability criterion [12] and will generally be smaller than the temporal-sampling period required to correctly temporally sample the time-domain fields.As a result, a temporal-sample time step is usually adopted as M∆t, where M satisfies Eq. ( 2).For the sake of clarity, we will however assume a temporalsampling period of ∆t and N t sampling periods.We note also that u 3 (r 3 , n∆t) is defined on a discrete spatial grid with points r 3 = (x 3 , y 3 ) = (i, j)∆ 3 , for integers i and j, as employed by the PSTD simulation.We shall assume that there are a total of I = J grid points in both the transverse directions.The discrete Fourier transform of the field on the observation plane will be sampled on a spatial-frequency grid defined by q 3 = (i, j)/(I∆ 3 ).Note that the spatialsampling period, ∆ 3 , is required by the PSTD method to be smaller than λ 0 /(2n s ).In typical PSTD implementations, a time-harmonic field is evaluated concurrently with the PSTD simulation according to: We now substitute Eq. ( 16) into Eqs.( 13) and ( 14), replace integrals with summations and convert Fourier transforms to discrete Fourier transforms to yield where r 3 and q 3 both take on their discrete forms.Either of Eq. ( 17) or Eq. ( 18) may be used depending upon the particular application.Consider first Eq.( 17), at each iteration of the PSTD algorithm, represented by the sum on n, three terms are required: φ 1 , K and u 3 .We now discuss whether these terms should be computed once per PSTD iteration, or computed once and then stored, with reference to computational efficiency.φ 1 is the Fourier transform of the detector sensitivity function which may be calculated analytically or numerically by a discrete Fourier transform.Whichever is used, it is likely that this term will be stored in computer memory for use at each iteration since it does not depend on the optical frequency f .The second term, K, is known analytically and can be computed relatively cheaply.Furthermore, storage of this term is likely to be inefficient in the case of broadband simulations since it depends upon f .When this is the case it is evaluated using Eq. ( 7), for each optical-frequency of interest, at each iteration of the PSTD algorithm.The final term, u 3 must be evaluated at each iteration of the PSTD algorithm by a two-dimensional discrete Fourier transform since it changes each iteration.
Consider now the alternative formulation, Eq. ( 18), which is principally useful in timeharmonic simulations since it requires knowledge of the optical-frequency dependent term F φ 1 f 2 f 1 q 3 K (q 3 ) .The Fourier transform in this term could be performed either analytically of numerically.If performed analytically, it is likely to result in a Fresnel integral, requiring numerical computation.Either way, it is only efficient to perform such computation for a small number of optical-frequencies, or, to store this term for use in each iteration of the PSTD algorithm.The advantage of using this formulation is that the sampling requirements of K are relaxed, as discussed in the following section.This is because F φ 1 can be calculated accurately, either analytically, or using discrete Fourier transforms.If discrete Fourier transforms are used, it may be evaluated on a grid of q 3 with sampling denser than what may be allowed by the transverse size of the PSTD simulation.This formulation is most efficiently implemented by computing F φ 1 f 2 f 1 q 3 K (q 3 ) once and storing the result for subsequent use in the PSTD simulation.This means that the transverse size of the PSTD simulation can be chosen independently of the need to sample K correctly.

Sampling requirements of K(q)
As foreshadowed in Sec.3.2, K(q) must be sampled at the Nyquist rate, or better, for Eqs. ( 17) and ( 18) to be evaluated correctly.This requirement comes directly from the use of a discrete Fourier transform in the case of Eq. ( 18).In the case of Eq. ( 17), however, it is necessary to use a result derived in Appendix 1, which shows that integration over the entire domain of a band limited function can be calculated exactly by summing all sampled values of that function.It is thus necessary that K(q) be sampled at, or better than, the Nyquist rate for Eq. ( 17) to be evaluated correctly.We note, however, that in some cases it may be necessary to employ a smoothed pupil function in order to ensure that K(q) is band limited, however, this was not found to be necessary in this work.We now derive a constraint upon the transverse size of the PSTD grid which ensures that K(q) is correctly sampled.
We follow Matsushima and Shimobaba [17] and write K(q 3 ) within the pupil as and thus Φ 1 (q 3 ) is given by If we now define q 3 = (q 3x , q 3y ), we can find the local frequencies of K(q 3 ) as [17] f q 3x = 1 2π it is straight forward to show that where q i can be either q 3x or q 3y .We define ∆q 3 = ∆q 3x = ∆q 3 as the sampling period of the spatial-frequency vectors in the focal plane common to both lenses in Fig. 1.In order to satisfy the Nyquist sampling criterion we require which, may be evaluated by first noting that 0 ≤ λ |q 3 | ≤ NA holds due to the aperture.We then note that | f q i |, as defined in Eq. ( 22), will in general have its maximum value at λ |q 3 | = NA, apart from when z obs is very close to the focus, in which case | f q i | has a maximum value within the domain 0 < λ |q 3 | < NA.It is easily verifiable, however, that the near-focus regime can be neglected since it results in only a negligible departure from the result predicted by assuming the maximum value of | f q i | occurs at λ |q 3 | = NA, and, within this regime, the sampling requirement of K is less restrictive than that of ensuring the sensitivity function is contained within the computational domain, as presented in Sec. 5.
In order to proceed, we make use of the rotational symmetry of Φ(q 3 ) and consider only a single component, q i , of q 3 .We now continue to substitute Eq. ( 22) into Eq.( 23), assuming that the maximum value of | f q i | occurs at λ q i = NA, giving the result: This, in turn, places a restriction on the transverse size of the PSTD simulation space since discrete Fourier transforms are to be used to transform from u(r 3 ) to u(q 3 ).In particular, if the PSTD grid is assumed to be square in a plane normal to the z-axis, and the square has a total width of max(x 3 ) − min(x 3 ), since ∆q 3 = 1/(max(x 3 ) − min(x 3 )) , we can say that max

Reduction of frequency sampling period by interpolation
In many cases of practical importance, the satisfaction of Eq. ( 25) may entail overly burdensome computational sizes.This may, however, be overcome by interpolating u 3 (q 3 ) resulting in a finer sampling and thus smaller sampling periods ∆q 3x and ∆q 3y .Note that so long as the field on the observation plane is sampled at, or better than, the Nyquist rate, this can be done without any approximation.Note that this requires that the observed field decays to zero at the edges of the observation plane, or an appropriate window function is employed.Interpolation can be implemented by multiplying the observed field by a linear phase function as described in [18], or the observed field data can be zero padded, prior to discrete Fourier transformation.We employed zero padding in our implementation.

Support of image of sensitivity function
In order for Eqs. ( 13) and ( 14) (and thus Eqs. ( 17) and ( 18)) to accurately simulate image formation, the transverse extent of the computational domain must be sufficiently large such that the detector sensitivity function, imaged onto the observation plane, is negligibly small at the boundary of the computational domain.In general, this means that φ 3 (r 3 ) = F φ 1 f 2 f 1 q 3 K (q 3 ) in Eq. ( 14) must be negligible at the boundary of the computational domain.The analysis may be simplified significantly in the common case in which a single mode optical fibre is employed.Then, noting that the projected sensitivity function must be identical to the focussed beam resulting from light emergent from the fiber, the projected sensitivity function may be evaluated as [19]: where θ 2 is uniquely determined by the expression n s sin θ 2 = sin θ 1 , J 0 (.) is the Bessel function of the first kind with order 0, F = ( f 2 / f 1 )W πNA/(2λ ) and W is the mode field diameter of the fiber.Note also that a scalar constant has been omitted from Eq. (26) for simplicity.One way in which an appropriate lower bound on the width of the computational domain can be found is to determine the minimum width, D, of a disc in the plane of observation such that which is analogous to saying that the computational domain is wide enough such that a fraction, 1 − ε, of the beam's energy is contained within it.It is important to note that this is only an approximation to energy since irradiance is not a conserved quantity [20], yet at low numerical apertures it is a good approximation.It is straightforward to evaluate R 2 |φ 3 (r 3 )| 2 d 2 r 3 as: This result is useful as it allows a constraint upon the minimum simulation size to be evaluated as: max(x 3 ) − min(x 3 ) ≥ 2D, where where I(D) is numerically evaluated according to and an appropriate value for ε is 0.01, such that 99% of the reciprocal focussed beam's energy enters the simulation space.

Evaluation and analysis
We present simulations which demonstrate the utility and properties of the developed method.
The first example verifies that the method correctly calculates the lateral detection point spread function (PSF) of a low numerical aperture (NA) scanning coherent optical microscope.We use this simulation to plot the limitations upon the transverse extent of the simulation required to sample the angular spectrum propagation term, K(q), correctly and also to ensure that the image of the sensitivity function is well contained within the simulation space.The second example demonstrates how the method may be used to calculate the axial PSF of the same optical microscope.We conclude with an example of how this method may be used to efficiently simulate image formation in a scanning coherent microscope for general samples.
All three examples are based upon the same primary simulation parameters.In particular, we model a detection system as depicted in Fig. 1 with the initial system parameters outlined in Table 6, which we understand to be representative of a typical optical coherence tomography system such as the Thorlabs Telesto-II with an LSM03 objective.In all three cases, normally incident plane wave illumination was assumed since our interest is in the detection, rather than illumination, PSF.For the remainder of this paper, PSF is thus taken to implicitly mean the detection PSF.

Lateral point spread function simulation
The lateral detection PSF was simulated by performing a series containing 23 psuedospectral time-domain (PSTD) simulations in which a single PSTD scatterer of refractive index 1.45 was scanned within an in-focus plane, through lateral positions r 3 = (iλ 0 /2, 0) for i = 0, 1, . . ., 22.
The PSTD simulation had a transverse size of 1112 × 1112 grid points, with an isotropic grid spacing of λ 0 /6, yielding a transverse width of just less than 241µm.An abnormally large transverse width was chosen so that interpolation of the field within the pupil was not required, as the method described by Eq. ( 13) was employed.A perfectly matched layer 10 cells thick, which absorbs, with low reflection, radiation incident upon it was employed in the PSTD simulation.The PML attenuates the field at the edges of the simulation space thus removing the need to employ a window function prior to performing discrete Fourier transformation, which we performed using the FFTW implementation of the fast Fourier transform [21].Although the mode field diameter (MFD) of the fiber we modelled is 9.2µm, for the purpose of illustration we additionally modelled MFDs of 1 and 18.4µm.We thus define three detector sensitivity functions φ 1 (r 1 ) = ϕ 1 (r 1 ), ϕ 2 (r 1 ) and ϕ 3 (r 1 ), corresponding to MFDs of 1µm, 9.2µm and 18.4µm respectively.We first present the quantities which form part of Eq. ( 13).In order to make interpretation clearer, since the non field-based quantities are rotationally symmetric in the pupil, we have plotted them only along the axis q y = 0 and only for positive values of q x .We have plotted the real part of K(q x ) in both the continuous and sampled cases.Note that for this simulation Eq. ( 25) dictates that the NA should not exceed approximately 0.35, corresponding with one of the vertical lines shown in Fig. 2. The meaning of Eq. ( 25) can be understood by noting that line λ 0 q x = 0.35 marks the point where it appears as though there is no longer one sample point per half-period of the chirped sinusoid.Note however, that for all MFDs except for that represented by ϕ 3 , this lack of sampling will not produce any inaccuracy in the detected field as K is apodized by φ1 and φ2 , respectively.Fig. 2. Plots of the real part of K(λ 0 q x ) in the continuous and discrete cases.K was sampled at the points denoted with dots and the dashed line added to enhance readability.The sensitivity functions transformed into the pupil have also been plotted for each simulated mode field diameter.The vertical black lines denote two NAs that were considered: the design NA of R a / f 2 and 0.35.
We now plot in Fig. 3(a) the lateral PSFs for each MFD which resulted from the PSTD simulations for an NA of R a / f 2 ≈ 0.097.To verify the accuracy the simulation we have plotted the analytically calculated PSF by noting that, in the low NA case, the detection PSF must match the reciprocal illumination PSF as given by the irradiance of the field emitted by the optical fiber and focussed into the sample space.This was calculated using a rigorous model of focussing theory developed previously [19].Very good agreement is achieved for all three mode field diameters resulting in a normalized mean square error, between the numerical and analytic results, of 4.5×10 −5 , 4.4×10 −6 and 1.4×10 −8 for sensitivity functions ϕ 1 , ϕ 2 and ϕ 3 , respectively.We define the normalized mean square error as where I an and I nu are the analytically and numerically evaluated PSFs and the x i are the coordinates of where I nu was evaluated.The PSFs resulting from NA = 0.35 are plotted in Fig. 3(b) resulting in normalized mean square errors of 3.1×10 −4 , 2.8×10 −7 and 1.3×10 −8 for sensitivity functions ϕ 1 , ϕ 2 and ϕ 3 , respectively.As expected, due to the apodization of K, there are no perceivable differences between the two NAs considered for the ϕ 2 and ϕ 3 cases.There is, however, a difference between the PSFs, for the ϕ 1 case, calculated analytically and using the PSTD method due to insufficient sampling of K in the pupil.The increase in normalized mean square error from 4.5×10 −5 to 3.1×10 −4 is small as the NA was set to the limit where sampling becomes inadequate.It was noted, however, that the error in the numerically calculated PSF was greatest for x i ≥ 5.5µm and when the normalized mean square error is recalculated for this range only it is observed to increase from 2.8×10 −3 to 1.1×10 −2 .
x (7m) The simulations employed a value of z obs which was only 10 PSTD cells above the focus in the simulation, thus why K is well sampled up until NA = 0.35.The required simulation size, however, increases as z obs departs from the focus.This is demonstrated in Fig. 4(a) which shows how the simulation size required to correctly sample K changes with z obs and such that the image of φ 1 in the sample space is well contained within the simulation space as measured by evaluating Eq. ( 29) for ε = 0.01.As previously discussed, the criterion described in Eq. ( 25) is inaccurate very close to the focus and this is the reason why the two plots in this figure have minima at slightly different values of z obs .We note, however, that this is of no practical significance since the minimum simulation size is determined by the greatest of the two plotted lower bounds, and the φ 1 constraint dominates for values of z obs in the vicinity of the focus whilst the K constraint dominates elsewhere.
z obs (7m) We now show in Fig. 4(b) how the two simulation size criteria vary with NA whilst z obs is fixed at the value used to evaluate the transverse PSFs presented in this section, for a fiber of MFD 9.2µm (φ 1 = ϕ 2 ).The plot shows that as NA increases, the K sampling criterion becomes dominant as is expected, since this is not constrained by the fiber mode ϕ 2 .The image of φ 1 support criterion is more restrictive at low NA since the lower the NA, the broader the image of φ 1 in the sample space.We note that at NA 0.09, the image of φ 1 becomes increasingly large as the Gaussian beam in the pupil becomes increasingly truncated by the aperture.At NA 0.09, there is negligible truncation by the aperture and so image of NA 0.09 is predicted by the theory of ideal Gaussian beams.

Axial point spread function simulation
We now consider the simulation of the axial PSF for the system considered in Sec.6.1.The simulation parameters were largely the same as in Sec.6.1 with the exception that the PSTD simulation had a transverse size of 350 × 350 PSTD cells and a depth of 1600 PSTD cells.A PSTD cell size of λ 0 /4 was employed to make the simulation of significant axial depths feasible.A total of 13 simulations were performed in which a single PSTD cell scatterer was scanned axially in increments of 130 PSTD cells (42.25µm), beginning at the in-focus position.The same MFDs as in Sec.6.1 were considered, leading to three distinct axial PSFs as plotted in Fig. 5.The analytic PSFs were evaluated using rigorous focussing theory developed previously [19], thus taking advantage of the reciprocity between illumination and detection point spread functions.We emphasize however, that plane wave illumination was employed in order to isolate the detection PSF.The normalized mean square errors between the numerical and analytic PSFs were 4.1×10 −4 , 2.1×10 −4 and 1.1 ×10 −5 for the three cases, respectively.Some deviation between the numerical and analytically evaluated PSFs is evident.One of the most evident discrepancies is in the ϕ 3 case, where the numerical PSF appears to oscillate slightly around the analytic PSF.This is due to the poor performance of the perfectly matched layer, which is supposed to absorb all radiation upon it.In practice, some of this radiation is reflected back to the detector interfering with the light from the scatterer.Another reason why the numerical PSFs underestimate the analytic PSFs is numerical dispersion.Numerical dispersion corresponds to the phenomenon by which, within the computational domain, a plane wave of given optical frequency ν, propagates in the PSTD grid with a dispersion relationship ñs (ν) = c/(ν λ ) where c is the speed of light in a vacuum.If we we instead calculate the analytic PSFs for the numerical values of wavelength and refractive index ( λ0 = 1307.2nmand ñs = 1.3925) a better match is obtained.Fig. 5. Plots of the axial PSFs for each MFD calculated analytically (solid line) and numerically using a PSTD simulation (dot markers) for NA = R a / f 2 .

Image formation example
We conclude this section by showing an example of how an image may be simulated using a single PSTD simulation.The image simulated was that of a shape stamped out of a sheet of dielectric material of refractive index 1.45 and thickness λ 0 /6.The resulting scatterer is shown in Fig. 6(a) where the gray corresponds to the material of refractive index 1.45.The white regions take on the background refractive index (1.42).A plane wave, linearly polarised at π/4 to the x (vertical) direction, was incident upon the scatterer.The magnitude of the incident plane wave components in the x and y directions were set to unity.17).This is less than than the 146 minutes required by the previously developed rigorous model [5].Although we were able to ommit the re-sampling step of the rigorous model [13] it could be made more efficient if it were optimized for low NA systems.The presented algorithm is, thus, numerically efficient relative to the rigorous model, particularly when multiple detector positions and optical wavelengths are considered for a single detector sensitivity function.This is because a significant cost of the entire simulation (i.e.including the PSTD simulation) is attributed to fetching data from the random access memory (RAM) into the central processing unit (CPU) cache memory.Once u 3 and φ 1 have been evaluated and stored in CPU cache memory, it is relatively fast to evaluate Eq. ( 17) for different values of λ and/or different detector positions.Computation would be slowed more significantly if different detector sensitivity functions were also considered.

Discussion and conclusions
In developing the numerical algorithm to be embedded in the PSTD and FDTD methods, we have also derived, in Eqs. ( 13) and ( 14), general formula which describe image formation in coherent optical imaging systems.Although these equations will be intuitive to many in the field, we are not aware of them having been previously derived in such a general manner.This new algorithm has reduced by one the number of components in a rigorous COI model.This has also enabled the amount of computer memory required to perform such simulations to be significantly reduced when considering broadband sources.This saving will enable the simulation of an OCT A-scan to be performed within a single graphical processing unit, which was previously unfeasible due to memory limitations.We thus anticipate that this algorithm will enable a significant increase in the uptake of rigorous simulation in OCT and related fields.
Future plans for this algorithm revolve around developing a variety of implementations of Eqs. ( 17) and (18).Our current implementation employs Eq. ( 17) in a manner whereby φ 1 is calculated using a discrete Fourier transform at the beginning of the PSTD simulation and K is calculated at each iteration of the PSTD algorith as required.We would like to allow for analytic forms of φ 1 to be available in the PSTD simulation.We would also like the alternative formulation, Eq. ( 18), to be available in the PSTD algorithm whereby F φ 1 K may be calculated using discrete Fourier transform or analytically.
An interesting direction is to consider a set of detector sensitivity functions, {ϕ i }, such that {ϕ i K} forms an orthonormal set.The set of detected {a i } values would allow the field in the PSTD simulation to be reconstructed to an an accuracy dependent on the construction of {ϕ i }.This would provide a way of storing scattered fields with a reduced sampling requirement.Finally, we would like to extend this algorithm to the high numerical aperture case where vectorial effects such as depolarisation by lenses cannot be neglected [14].

Fig. 3 .
Fig. 3. Plots of the lateral PSFs for each MFD calculated analytically (solid line) and numerically using a PSTD simulation (dot markers) for NA = R a / f 2 (a) and NA = 0.35 (b).

? 1 Fig. 4 .
Fig. 4. Plot of minimum simulation size as a function of z obs (a) and of NA (b) required to sample K correctly and to contain the image of φ 1 = ϕ 2 within the sample space.
Figures 6(b) and (c) show the magnitude of the x and y scattered field components on a plane 2λ 0 above the top of the scatterer.Figures 6(d)-6(f) show images simulated using the parameters listed in Table 6, but for different values of z obs in order to simulate defocussing.In particular, for each of the images Figs.6(d)-6(f), an additional 0, 200µm and 400µm, respectively was added to z obs .Each image was simulated by implementing the scanning technique outlined in Sec.3.1 using an unnecessarily high number of values, r d , in order to obtain a high resolution image.Note, however, that although the PSTD simulation had a transverse dimension of 1400×1400, the resulting image contains only 175×175 pixels, a 64 fold reduction in the computer memory requirement.The simulations were performed on a workstation with two Intel ® (Santa Clara, California, USA) Xeon ® E5-2650V2 processors.Each image in Figs.6(d)-6(f) took approximately 125 minutes to compute using Eq. (

Fig. 6 .
Fig. 6. a) through the simulated scattering object which was λ 0 /6 thick and embedded in a material of refractive index 1.4.b) and c) show the magnitudes of the x and y components of the scattered field just above the scatterer.d) -e) show images when the scatterer was 0, 200µm and 400µm, respectively, from the focus.

Table 1 .
Base parameters of the numerical simulations, λ 0 is the wavelength in air and MFD stands for mode field diameter.