1. Introduction

While developing optical microscope systems we experienced the difficulty to compare a new system with the previous setups in terms of speed, resolution, sensitivity and temporal behavior. To shortcut the trial and error phase of such a development we searched for a method which would allow to account for most linear optical effects in light propagation. Furthermore the method should include a simple representation of the factors contributing to the results and easily being implementable to computational simulations.

Therefore we chose to establish a flexible numerical simulation method accounting for all properties of the electromagnetic wave equation in non conducting and therefore non-absorbing media. Thereby this method is applicable to basically any far field microscope.

In an insulating, uncharged and homogeneous (in space and time) medium the electromagnetic wave equation is [1]:

Δ is the Laplace operator, E⃗ the electric field, n the refractive index of the media, c the speed of light and (r⃗, t) the space and time coordinates.

Any solution to Eq. (1) can be written as a superposition of an arbitrary amount of plane waves $\overrightarrow{\tilde{E}}(\overrightarrow{k},\omega )e^{i\left(\overrightarrow{k}\overrightarrow{r}-\omega t\right)}$ with frequency ω and wave-vector k⃗ [1]:

Equation (2) can be interpreted as the Fourier transformation of $\overrightarrow{\tilde{E}}(\overrightarrow{k},\omega )$ to $\overrightarrow{E}(\overrightarrow{r},t)$. The signs for plane waves and the Fourier transformation differ. By substituting $\overrightarrow{\tilde{E}}(\overrightarrow{k},\omega )$ with $-\overrightarrow{\tilde{E}}(-\overrightarrow{k},\omega )$ this problem can be omitted. For a better readability we chose to continue with $\overrightarrow{\tilde{E}}(\overrightarrow{k},\omega )$ in the following text. Since E⃗(r⃗, t) represents the electric field (and with the relation $I(\overrightarrow{r},t)\propto {\left|\overrightarrow{E}(\overrightarrow{r},t)\right|}^2$ also the intensity) in real space and time, e.g. in the focal volume, it also represents the central information we want to obtain with our simulation. Therefore $\overrightarrow{\tilde{E}}(\overrightarrow{k},\omega )$ describes the complete light field of the focal volume.

In Fourier optics the known back focal plane (BFP) of e.g. a microscope objective can be used to determine $\overrightarrow{\tilde{E}}(\overrightarrow{k},\omega )$ [2, 3]. To understand this step, it is convenient to start from a single monochromatic point source in the BFP. A point source in a focal plane of an ideal optical system, like a collecting lens, is collimated to a parallel beam, which is the equivalent of a single plane wave. The direction of the collimated beam is given by the ray projection function $g(\mathrm{\Theta })=\raisebox{1ex}{$\rho $}\!\left/ \!\raisebox{-1ex}{$f$}\right.$ where Θ is the angle between the propagation direction through the focus (in the focal volume around r⃗ = 0⃗) and the optical axis [3]. For common simulation situations the width of the collimated beam is small compared to the optics, which allows to neglect higher order aberrations in this assumption. $\rho =\left|\overrightarrow{\rho }\right|=\sqrt{x^2+y^2}$ describes the distance of the point source to the optical axis in the BFP (ρ_z = 0) with $\overrightarrow{\rho }=\left(\begin{array}{c}{\rho }_x\\ {\rho }_y\\ {\rho }_z=0\end{array}\right)$and f the focal length of the lens system. Most microscope objectives follow the sine condition g(Θ) = sinΘ which is assumed to hold for all following considerations [3]. With the known projection function it is possible to connect one point with known frequency of the back focal plane E⃗_BFP(ρ⃗, ω) to a plane wave and therefore one point of $\overrightarrow{\tilde{E}}(\overrightarrow{k}(\overrightarrow{\rho }),\omega )$. In spherical coordinates k⃗ is expressed as

With this relation all points E⃗_BFP(ρ⃗, ω) in the back focal plane (BFP) determine $\overrightarrow{\tilde{E}}(\overrightarrow{k},\omega )$.

S(ρ⃗) represents a scaling factor depending on the position (ρ⃗) in the back focal plane. We will return to the properties of S(ρ⃗) in the Gridding section.

This relation now allows to implement a numerical simulation, which by the knowledge of the electrical field properties in the back focal plane can derive the spatial and temporal behavior of the light field in the focal volume. To have an efficient implementation it is necessary to apply the Fourier transformation by using a fast Fourier transformation (FFT) algorithm. Therefore the simulated volume needs to be sampled in an equidistant manner.

Gridding

Due to sampling constrains E⃗_BFP(ρ⃗, ω) and $\overrightarrow{\tilde{E}}(\overrightarrow{k},\omega )$ are only known at discrete positions. The conversion with Eq. (3) makes it very unlikely, that the known positions of E⃗_BFP(ρ⃗, ω) coincide with the known grid positions of $\overrightarrow{\tilde{E}}(\overrightarrow{k},\omega )$. Therefore it is necessary to introduce an interpolation method, distributing contributions from single points (ρ⃗, ω) to their closest positions in (k⃗, ω). This approach is also known as gridding [4–7]. The interpolation combined with the FFT should deliver a result similar to the result obtained by solving the problem with a non equidistant matrix. This special combination has come to be known as non uniform FFT (NUFFT) and has been applied to various tomographic reconstruction techniques [8].

We chose four representative interpolation methods to be compared among each other: nearest-neighbor-(NN-) NN(s), cubic-spline-(CS-) [9] CS(s), scaled-windowed-sinc-(SWS-) SWS_H,ξ (s) [10–12] and power-of-hann-(POH-)interpolation ${\text{POH}}_{\text{H}}^{\text{L}}(s)$. In the context of E⃗_BFP(ρ⃗, ω) as source and $\overrightarrow{\tilde{E}}(\overrightarrow{k},\omega )$ as destination, the interpolation was always assumed to be destination oriented. This was achieved, by introducing a folding of the destination points $\overrightarrow{\tilde{E}}(\overrightarrow{k}(\overrightarrow{\rho }),\omega )$ with the interpolation function F_I(s) in k-space performing a subsequent grid sampling. Therefore the interpolation function F_I(s) will determine the values V_i at the grid positions x_i depending on the destination position x_des normed to the distance Δx between two grid points in the destination space.

The interpolation functions are defined as follows:

H is the width of the interpolation kernel, ξ the scaling factor for the sinc kernel, a_i the coefficients used for the window function [13] and L the power of the power-of-hann (POH) kernel.

To achieve consistent results it also is necessary to consider the non uniform distributions of the interpolated points. Thus the scaling function S(ρ⃗) has to be taken into account. The relative scaling from sample to sample can be derived by introducing a conservation of energy or, as an equivalent in the case of propagating light, power constrain [3].

We choose one grid-point ρ⃗ to represent one light beam, with a constant power P. This light beam will be transformed from one representation (e.g. BFP) to another (e.g. k-space (KS)). By this transformation the area δS_BFP being passed by the light beam in the BFP changes to another area δS_KS after the transformation to k-space. In a ray representation these areas are inversely proportional to the ray density $\propto \frac{1}{\delta S}$. To create a measure for the area, we use the absolute value of the cross product of the derivatives ∂_ρx, ∂_ρy along the perpendicular axes in the BFP. With $\overrightarrow{\rho }=\raisebox{1ex}{$\overrightarrow{\rho }$}\!\left/ \!\raisebox{-1ex}{$f$}\right.$ the normed position-vector in the BFP and ρ_z = constant both areas can be computed

Similar calculations can be found in [3]. Due to the type of interpolation, the sample density coincidences with the ray density and therefore implies the area scaling calculated above. A change in sample distances results in a scaled image after the FFT as expected by the Fourier scale property [14]. With the chosen interpolation approach, the number of samples does not change with the shrinking and enlarging of the mapped area. Due to this discrepancy to the scaling with a discrete Fourier transformation the sample values have to be prescaled by their local density. This can be approximated by the calculated areas δS_BFP and δS_KS. Therefore a numerical scaling error S_num(ρ⃗) of the E-field $\left|{\overrightarrow{E}}_{\text{i}}\right|=\frac{{\left|{\overrightarrow{E}}_{\text{i}}\right|}_{\mathit{num}}}{\delta S_i}$δS is introduced with i = KS or BFP and a scaling factor δS for physical consistency. Introducing these E-field equivalents, analog calculations to Eq. (9):

Power-of-hann-interpolation

The naming of the “power of hann interpolation” (following the concept of “power of two”) was chosen to describe the function ${\text{POH}}_{\text{H}}^{\text{L}}(s)$ which is a Hanning window: ${\text{cos}}^2\left(\pi \frac{s}{H}\right)$, to the power of L. Other publications choose Kaiser-Bessel and Gaussian windows for their kernel. The POH kernel has some advantageous properties, which can be derived from the publication by Nuttall [13].

Primary its first 2L − 1 derivatives at $\pi \frac{s}{H}=\frac{\pi }{2}$ are zero, which allows for a smooth blend out to the zero function outside the kernel area and therefore steep decline after their first dip in Fourier space as displayed in Fig. 1(a). Secondly the first dip can be placed on the edge of the destination space or Nyquist frequency by choosing H = 2L + 2 as shown in Fig. 1(b). This ensures the side-lobes outside the destination space of the Fourier-transformed kernel to drop-off steeply. Therefore $L=\frac{H-2}{2}$ is chosen for any following ${\text{POH}}_{\text{H}}^{\text{L}}={\text{POH}}_{\text{H}}^{\frac{\text{H}-2}{2}}$-kernel. Thirdly the computational cost of the kernel is reduced compared to a three term cosine series or higher order polynomial sums. Only one cos-function has to be evaluated followed by L or, in the case of clever usage of intermediate results (e.g. a = (x * x) ⇒ x⁴ = x * x * x * x = a * a, needs two instead of three multiplications.), less multiplications.

The kernel is freely scalable (see Fig. 1(a)), fast calculable and shows high quality drop-off behavior in Fourier space (see Fig. 1(b)). This is bought by attenuated values respectively false amplitudes after the Fourier transform as displayed in Fig. 1(b). These have to be accounted for by normalizing to the average amplitude error and subtracting the average phase error from the “final” result as explained in the Impulse response section.

Due to the fast drop to zero at the border a division by small numbers is provoked, which becomes worse for bigger kernels. Therefore the scalability is limited by the floating point accuracy (e.g. float vs. double) and in addition the border areas do not contain useful information.

 

Figure 1 Impulse response for multiple widths of the POH interpolation. a) Several POH kernels at the central grid position. b) Impulse response after Fourier transformation with identical color code. The resulting normed amplitude show the behavior of a ${\text{POH}}_{\text{H}}^{\frac{\text{H}-2}{2}}$ kernel for H = 4, 6 and 12. The first drop in frequency domain stays consistent on the Nyquist frequency $0.5\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$\text{grid}\hspace{0.17em}\text{points}$}\right.$. For bigger H the kernel broadens and drops at steeper rates for extreme frequencies above Nyquist.

Download Full Size | PDF

2. Methods

Impulse response

To investigate the behavior of different interpolation methods the numerical results of a one dimensional single point interpolation with a subsequent FFT can be compared to the expected, analytical Fourier transform of the corresponding delta peak δ(k_x − k):

A plane wave with a constant amplitude and phase defined by the given circular wave number k is expected. In other publications the amplitude of the analytical Fourier-transform has been at the center of investigations [11, 15, 16]. However the following results, especially for the nearest neighbor method, will show that for the case of an FFT with limited samples not only the amplitude but also the phase is useful to describe the whole quality of the interpolation.

The analysis was implemented using a one dimensional grid with 256 elements all set to zero. Then the value of one point, applying one of the described interpolation kernels, is set close to the zero position as displayed in Fig. 2(a). Thereafter the vector is transformed with a one dimensional FFT. This procedure is repeated 500 times at equally spaced positions between −0.5 and +0.5 grid points (distance between two elements) around the zero position. Some examples for the intermediate results are displayed in Fig. 2(b). By this all possible responses of the interpolation can be represented. Interpolation outside the range of −0.5 to +0.5 grid points will, besides of numerical accuracy limitations, lead to redundant results. Since the whole vector can be cyclic permuted, this would only add a steeper phase due to the larger k.

The difference of the results to the ideally expected e^−i(kr) are displayed (see Fig. 3) for further analysis. A worst case estimation (biggest deviation from the analytical value at any of the 500 interpolation positions) of the error can be extracted into one curve plotted over the positions in the Fourier-transformed grid. As post-FFT correction, to account for the needed normalization of the POH interpolation, every grid-position is normed to the average amplitude attenuation and corrected for the average phase deviation of the 500 translation samples.

 

Fig. 2 Impulse response for multiple positions p_i of the POH interpolation. a) ${\text{POH}}_6^2$ kernels at three different positions. The dots indicate the values used for the interpolation. b) Impulse response after Fourier transformation with identical color code. The resulting amplitude and phase distributions differ from position to position. For the amplitude distribution a constant (0 dB) would represent the ideal frequency response. The phase response should result in a linear phase with a slope proportional to the displacement Δ spanning from −90 to 90 deg for Δ = ±0.5. Differences from the ideal expected behavior have to be compensated by normalization. This approach is limited by the variations e.g. in the extreme frequencies of the normed amplitude.

Download Full Size | PDF

Similar normalizations also called grid correction or deconvolution were performed in [5, 15].

Oversampling

The second method to interpret the interpolation accuracy is to perform an intensity simulation of a fs-focus at different base sizes. By introducing oversampling it is assumed, that the result yields a higher accuracy. For the simulations at the standard sampling a grid size of 512 × 512 × 512 voxels (k_x × k_y × k_z) was chosen. The oversampled version was reduced in one dimension to achieve sufficient oversampling within memory restrictions and set to the x,z-plane with a grid of 4096 × 1 × 4096 voxels. The differences of the low sampled and high sampled version, both scaled to reach a maximum value of one, are calculated as absolute error (see Fig. 4(b)). In a second step the difference is divided by the high sampled version to obtain the relative error (see Fig. 4(c)). To omit differences due to the sampling in the back focal plane the sample dimensions of the BFP were fixed to 100 × 100 × 40 voxels (ρ_x × ρ_y × δω) dictated by the simulation parameters and the standard sampling grid. The physical properties of the simulation were chosen to base on a fs-laser light spectrum spanning λ = 700nm...900nm, with an amplitude behavior of a cos²(ω(λ))-window. Furthermore the simulation volume had the dimensions 50 μm × 50 μm × 50 μm and an NA of 0.2 was chosen. The ${\text{POH}}_6^2$-kernel was chosen for interpolation in all four dimensions.

Equidistant gridding

As a third method the focus simulation was compared by using an analytically defined $\left|\overrightarrow{\tilde{E}}(\overrightarrow{k},\omega )\right|$ not accounting for polarization effects. It was chosen to simulate the intensity of a fs-focus while minimizing aliasing effects introduced by the chosen field itself. Therefore cos² and sin² were introduced:

Spectral vortex simulation

To demonstrate the potential of this simulation approach the three dimensional time evolving field of a spectral optical vortex was simulated. The NA was chosen to be 0.6 and the spectrum of the mode locked laser-source spans λ = 779nm...850nm. This spectrum is distributed linearly with the polar angle (|k⃗_center| ∝ φ) in the BFP. At any point the local spectrum Δk_local is restricted to 10 % of the full spectrum Δk_local = 0.1 · Δk_laser of the fs-laser source. The local spectrum was chosen to show the amplitude behavior of a cos²(ω(λ))-window. The BFP is sampled with 100 × 100 × 40 voxels (ρ_x × ρ_y × δω) and the focal volume is sampled with 256 × 256 × 256 × 64 voxels (k_x × k_y × k_z × t), with the space and time dimensions of 50 μm × 50 μm × 50 μm × 2 ps. The ${\text{POH}}_6^2$-kernel was chosen for interpolation in all four dimensions.

Hard- & software

All simulations were performed on an Kubuntu based Linux system. The simulation was implemented in C++ with the FFT using the FFTW package. The system possessed two Intel Xeon E5603 1,6 GHz 4-Core CPUs and 24 GB of DDR3 ECC RAM.

3. Results

The analysis of the maximum error comparison, after compensating for the average amplitude attenuation and phase deviations as explained in Impulse response, is displayed in Fig. 3.

 

Fig. 3 Amplitude and phase error at all frequencies after a one dimensional FFT with normalization to constant errors. a) Maximum deviation from the ideal amplitude of 1 in any of the 500 tested interpolation positions after normalization to average amplitude after the FFT. b) Maximum deviation from the ideal phase −i(kr) in any of the 500 tested positions after subtraction of the average phase error after the FFT.

Download Full Size | PDF

The relative amplitude error of the nearest neighbor interpolation is zero (out of scaled region) at all grid positions since the value is transferred to the next useable grid point. However the relative phase error yields the largest error among all compared interpolation methods. The majority of the interpolation kernels deliver a similar behavior for both error types of amplitude and phase. Apart from the nearest neighbor method, the POH-interpolation with a kernel size of 12 shows inconsistent behavior as well.

 

Fig. 4 Top: ${\text{POH}}_6^2$ interpolation method compared to oversampled version. Bottom: ${\text{POH}}_6^2$ interpolation method compared to equidistant sampling of an analytical $\left|\overrightarrow{\tilde{E}}(\overrightarrow{k},\omega )\right|$. a) & d): Results of the ${\text{POH}}_6^2$ interpolation method. b) & e) The absolute error calculated as the difference between the standard simulation and the reference simulation, both with their maximum intensity normed to one. c) & f) Relative error calculated by dividing the absolute errors (b & e) by their reference simulations.

Download Full Size | PDF

By changing the kernel position s from float to double floating point accuracy and leaving the value representation ${\text{POH}}_{12}^5(s)$ at float accuracy, the difference is avoided and again a similar behavior is achieved.

Figs. 4(a)–4(c) display the results from the comparison between simulation results at two different sample rates. Samples from a border with 1% width of the simulation are not displayed due to undefined behavior after normalization (compare Power-of-hann-interpolation section). The simulation result displays intensities normed to the peak intensity. The absolute difference, for normed central maxima, as explained in the 2 section, we will denote as error. The error does not rise above 10⁻⁵. The error relative to the simulation values of the oversampled version, as described in the 2 section, is only at 10⁻⁵ for the center values (see Fig. 4(c)). Towards the borders there is a continuous rise and relative errors above 100% are not displayed. Figs. 4(d)– 4(f) displays the results from the comparison between simulation results at two different sample rates while the higher sample rate uses an analytical $\left|\overrightarrow{\tilde{E}}(\overrightarrow{k},\omega )\right|$ function to create an equidistantly sampled E⃗-field. Samples from a border with 1% width of the simulation are not displayed due to undefined behavior after normalization. The results are quite similar to the comparison by simple oversampling. The absolute difference/error, for normed central maxima, does not rise above 10⁻⁵ while the error relative to the simulation values of the oversampled version is at 10⁻⁵ only for the center values. To the borders there is a continuous rise. Relative errors above 100% are not displayed.

 

Fig. 5 Simulation of a spectral vortex. Side by side view of the x,z-plane and the x,y-plane (indicated by the red dashed lines) with the time evolving from top left to bottom right. The diameter of the vortex reduces over time. After a complete collapse it starts expanding again. Scale bar represents 10 μm. (see Media 1

Download Full Size | PDF

To demonstrate a more complex setting, an optical vortex has been simulated. In this case the vortex is not rising with a rotary phase [17] rather than with a rotary spectrum. This leads to a focal distribution which expands and collapses in time with the vortex-typical shape of a doughnut mode (see Fig. 5). At close examination the doughnut shape is interrupted at several occasions. These interruptions mainly correlate with the direction in which the discontinuity of the rotary spectrum takes place. Furthermore the highest intensity is reached for the collapsed focus, while the expanding focus distributes the field power over the whole doughnut.

4. Discussion

To interpret the results shown in Fig. 3 one at first may notice that the utmost interest in simulating a focal distribution lies in the correct knowledge of the amplitude-point spread function of the focal volume. Therefore the phase error might be dismissed as non contributing. This is not true for this type of analysis scheme. If for example two independent points are interpolated, their corresponding plane waves will interfere after the Fourier transformation. For equally high amplitudes and destructive interference the relative phase error may produce a remaining term similar in amplitude to the relative amplitude error. Therefore both error types can be estimated to contribute equally to the final simulation error. This as well explains the strong coincidence of both error types for almost all interpolation kernels.

Calculating the sum of both error types and counting the number of grid points at which one kernel results in the smaller error sum, the kernels can be arranged in their quality (the higher the error the less the quality) in ascending order:

In other words the further the interpolation kernel is arranged to the “right” the more grid points have a lower error.

The nearest-neighbor as well as the cubic-spline-interpolation methods are only suitable at the center of the simulated region. Therefore the step to interpolation methods with a higher accuracy is the logical consequence for the pursued task. This typically implies the use of sinc-type kernels which preserve the higher frequencies and therefore reduce the error at the borders of the spectrum. By rescaling after Fourier transform, as applied in Fig. 3, the ideal preservation of high frequencies is no longer necessary. Rescaling allows to use kernels with low responses at the higher frequencies, but with a steeper drop-off behavior above their sampling limits which leads to comparable or, in the case of the POH interpolation, even less erroneous results than with the sinc interpolation. Nevertheless the mayor benefit is the reduction of the kernel sizes, in this example 6 and 12, reach comparable results to an 11 and 31 elements wide sinc kernel.

The practical suitability of this interpolation method has been shown in the comparisons of simulations between normal and oversampled results. Figure 4 displays the absolute as well as the relative errors produced for a fs-pulsed laser focus. The primary comparison of the absolute error shows that the simulation combined with the POH-interpolation is capable of reaching an accuracy well below 10⁻⁵ of the simulated peak intensities up close to the border of the simulated region. This is surprising since the single point analysis suggested a maximum accuracy of 10⁻³ for amplitude or 10⁻⁶ for the intensity relative to the normed wave amplitude. Therefore, the use of several interpolation points to create a $\overrightarrow{\tilde{E}}(\overrightarrow{k},\omega )$ increases the dynamic range, which can be seen in the simulated data allowing the intensity values to drop to 10⁻¹¹. Nevertheless, this steep drop reduces the relative accuracy of the simulation significantly. By normalizing the absolute error to the simulation data the relative error can be displayed. Again the accuracy drops to 10⁻⁵ and rises continuously to the borders. As simulation values drop below 10⁻⁹ the relative error rises above 1. If relative accuracy is crucial this may lead to high inaccuracies and should always be taken in to account. By this the simulation is limited to values above 10⁻¹⁰ of the peak amplitude of the focus. Close to the simulation borders (< 10% of width) this accuracy degrades rapidly for all attempted analysation methods and should be avoided by choosing a slightly bigger simulation volume.

As an example a spectral vortex has been introduced. The behavior in the focal volume can be explained by examine the center wavelength of the spectrum at a given angle φ in the back focal plane. The spectrum was designed to give a rise to the wave-vector proportional to the angle φ. For a mode-locked fs-spectrum this can be interpreted as all modes having the same phase at the start time t₀. This phase ΔΦ rises proportional to the frequency separation Δω and the time t, ΔΦ = tΔω, which creates a typical phase vortex for any given time t. Therefore the slope of the phase rises in time and blends through several vortices until the highest slope representable by the used modes is reached and the process reverses. This behavior is reproduced in the simulation (see Fig. 5). As expected the more the doughnut mode expands the dimmer it becomes, due to the light distributing over a bigger area. Nevertheless, more complex $\overrightarrow{\tilde{E}}(\overrightarrow{k},\omega )$ might achieve a compensation for this effect and will be at interest for future investigations.

5. Conclusion

We were able to demonstrate a method to efficiently simulate solutions of the wave equation concerning the focus formation of a fs-laser scanning microscope. Furthermore by interpolation analysis and oversampling we were able to demonstrate absolute error rates (difference between the simulation and a reference normed to one) below 10⁻⁵ with a fast POH based interpolation. Finally we demonstrated the possibility of an “expanding optical vortex” with a time dependent radius.