1. 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–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–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–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.

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

 

Fig. 1 Schematic diagram of the optical system under study. Each position vector r₁, r₂ and r₃ is a two dimensional vector representing transverse displacement relative to a location on the optical axis. The field is observed on plane z₃ = z_obs and there is an infinite half-space of material with refractive index n_s beginning at z₃ = −h. z₃ is measured relative to the nominal focus of the right hand lens, in air. Angles θ₁ and θ₂ represent the semi-convergence angles of typical rays. R_a is the radius of the aperture.

Download Full Size | PDF

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 two-dimensional 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.

3. 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₂/f₁)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., Ũ).

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 time-harmonic 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₃(r₃, f) can be extracted at observation points r₃ and temporal frequencies f, from the time-domain field u₃(r₃, nΔt) according to

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₃ = z_obs and is denoted U₃(r₃, z_obs) where r₃ 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₃ = −h. Our objective is to evaluate

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.

3.1. 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:

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₃(r₃, 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 temporal-sampling period of Δt and N_t sampling periods. We note also that u₃(r₃, nΔt) is defined on a discrete spatial grid with points r₃ = (x₃, y₃) = (i, j)Δ₃, 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₃ = (i, j)/(IΔ₃). Note that the spatial-sampling period, Δ₃, is required by the PSTD method to be smaller than λ₀/(2n_s). In typical PSTD implementations, a time-harmonic field is evaluated concurrently with the PSTD simulation according to:

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: ϕ̃₁, K and ũ₃. We now discuss whether these terms should be computed once per PSTD iteration, or computed once and then stored, with reference to computational efficiency. ϕ̃₁ 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, ũ₃ 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 time-harmonic simulations since it requires knowledge of the optical-frequency dependent term $\mathit{ℱ}\left\{{\tilde{\varphi }}_1\left(\frac{f_2}{f_1}{\mathit{q}}_3\right)K({\mathit{q}}_3)\right\}$. 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 $\mathit{ℱ}\left\{{\tilde{\varphi }}_1\left(\frac{f_2}{f_1}{\mathit{q}}_3\right)K({\mathit{q}}_3)\right\}$ 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₃ with sampling denser than what may be allowed by the transverse size of the PSTD simulation. This formulation is most efficiently implemented by computing $\mathit{ℱ}\left\{{\tilde{\varphi }}_1\left(\frac{f_2}{f_1}{\mathit{q}}_3\right)K({\mathit{q}}_3)\right\}$ 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.

4. 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₃) within the pupil as

In order to proceed, we make use of the rotational symmetry of Φ(q₃) and consider only a single component, q_i, of q₃. We now continue to substitute Eq. (22) into Eq. (23), assuming that the maximum value of |f_qi| occurs at λq_i = NA, giving the result:

4.1. 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 ũ₃ (q₃) 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.

5. 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 ${\varphi }_3({\mathit{r}}_3)=\mathit{ℱ}\left\{{\tilde{\varphi }}_1\left(\frac{f_2}{f_1}{\mathit{q}}_3\right)K({\mathit{q}}_3)\right\}$ 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]:

6. 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 1, 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.

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

View Table

6.1. 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₃ = (iλ₀/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 λ₀/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 ϕ₁(r₁) = φ₁(r₁), φ₂(r₁) and φ₃(r₁), 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 λ₀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 φ₃, this lack of sampling will not produce any inaccuracy in the detected field as K is apodized by φ̃₁ and φ̃₂, respectively.

 

Fig. 2 Plots of the real part of K(λ₀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₂ and 0.35.

Download Full Size | PDF

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₂ ≈ 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⁻⁵, 4.4×10⁻⁶ and 1.4×10⁻⁸ for sensitivity functions φ₁, φ₂ and φ₃, respectively. We define the normalized mean square error as $\epsilon ={\sum }_i{\left|I_{\mathit{an}}(x_i)-I_{\mathit{nu}}(x_i)\right|}^2/{\sum }_i{\left|I_{\mathit{an}}(x_i)\right|}^2$, 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⁻⁴, 2.8×10⁻⁷ and 1.3×10⁻⁸ for sensitivity functions φ₁, φ₂ and φ₃, respectively. As expected, due to the apodization of K, there are no perceivable differences between the two NAs considered for the φ₂ and φ₃ cases. There is, however, a difference between the PSFs, for the φ₁ 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⁻⁵ to 3.1×10⁻⁴ 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⁻³ to 1.1×10⁻².

 

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₂ (a) and NA = 0.35 (b).

Download Full Size | PDF

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 ϕ₁ 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 ϕ₁ constraint dominates for values of z_obs in the vicinity of the focus whilst the K constraint dominates elsewhere.

 

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 ϕ₁ = φ₂ within the sample space.

Download Full Size | PDF

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 (ϕ₁ = φ₂). 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 φ₂. The image of ϕ₁ support criterion is more restrictive at low NA since the lower the NA, the broader the image of ϕ₁ in the sample space. We note that at NA ≲ 0.09, the image of ϕ₁ 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.

6.2. 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 λ₀/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⁻⁴, 2.1×10⁻⁴ and 1.1 × 10⁻⁵ 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 φ₃ 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 (λ̃₀ = 1307.2nm and ñ_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₂.

Download Full Size | PDF

6.3. 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 λ₀/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. Figures 6(b) and 6(c) show the magnitude of the x and y scattered field components on a plane 2λ₀ above the top of the scatterer. Figures 6(d)–6(f) show images simulated using the parameters listed in Table 1, 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.

 

Fig. 6 a) Cross-section through the simulated scattering object which was λ₀/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.

Download Full Size | PDF

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. (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 ũ₃ and ϕ̃₁ 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.

7. 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 ϕ̃₁ 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 ϕ̃₁ to be available in the PSTD simulation. We would also like the alternative formulation, Eq. (18), to be available in the PSTD algorithm whereby ℱ̂{ϕ̃₁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 {φ_iK} 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].