BIM-Sim: Interactive Simulation of Broadband Imaging Using Mie Theory

1. Introduction

Particle interactions within an incident electromagnetic (EM) field are difficult to simulate and visualize because the result is inherently high-dimensional. EM radiation is generally represented using a complex vector field in three spatial dimensions. When the incident light consists of broadband radiation, each wavelength λ interacts with materials in a unique way. Studying the behavior of EM fields generates a time-dependent component, as the user changes properties of the particles or incident light in order to study these effects on the scattered field. Accounting for all of these factors requires a five-dimensional (5D) simulation in (x, y, z, λ, t), which is both computationally expensive and requires a prohibitive amount of memory.

In this paper, we propose software for interactive simulation of broadband Mie scattering in EM fields. We describe a mathematical model based on Mie theory that allows efficient evaluation of the scattered field using hardware and algorithms that support highly parallel computation. The high-dimensional nature of the simulation is constrained by user-specified visualization parameters. By taking advantage of GPU-based computation, we allow the user to interactively adjust these parameters to perform the simulation in real-time. This has several advantages over traditional simulation of the complete 5D forward model.

1.1. Applications

Mie theory plays an important role in describing the propagation of EM radiation by providing a rigorous solution to Maxwell's Equations for a spherical scatterer. This theory has been used to characterize absorption in atmospheric simulations [1], and to approximate cellular structures in spectroscopic imaging of biological tissues [2]. There is also interest in using micro and nanospheres, such as quantum dots, to improve sensing in biomedicine [3]. Generalized Mie theory [4, 5] has been used in plasmonics applications to understand the near-filed response of nanoparticle aggregates to circularly polarized light [6, 7].

Mie theory is often used as a first-order approximation for scattering by general particles [8]. Of particular interest is the problem of inverse scattering in samples composed of particulate or spherical components. We have recently demonstrated a method for computing the refractive indices of spheres composed of polymers measured using mid-infrared point spectroscopy [9]. However, these nonlinear solutions require iterative computation of the forward model. Since Mie theory is computationally expensive, this solution becomes impractical for multidimensional spectroscopic images. In particular, a fast forward model is necessary for iterative calculations of inverse solutions [9], characterization of imaging systems with multiple interacting optical components [10], and iterative fitting of large numbers of particles [11].

1.2. Previous Work

Several applications have been developed for computing spectra based on the scattering properties of spheres. Algorithms like MiePlot [12] provide plotting functions for scattering amplitude as a function of wavelength, particle size, and material properties. In addition, many plotting algorithms are available for creating Mie scattering plots for specific cases such as coated [13] and multi-layered spheres in shaped beams [14]. Finite element methods are generally used for non-spherical particles [15], while general solutions are known for fibers [16, 17]. These computational methods are often based on the BHMIE method proposed by Bohren and Huffman [8]. Open-source implementations are also available [18].

However, the available tools focus on creating plots that characterize scattering in highly specific cases, which include incident plane waves or fields focused on the sphere. These cases generally allow for a closed-form analytical solution for the scattered field. We have found no interactive tools for exploring the multidimensional scattered field produced by spheres. This is primarily due to (a) the computational complexity of evaluating a scattered field, (b) the large amount of memory required to store a pre-computed field, and (c) a lack of visualization methods available for exploring simulation results.

Recent methods have been proposed for creating a general framework for simulating scattering in microspectroscopy [16, 19], with the goal of correcting artifacts. These methods evaluate cross-sections of the EM field near the sample. However, these techniques are currently very time-consuming and focus only on the evaluation of scattering through planar substrates.

Figure 1 shows an example of a polar scattering plot and a near-filed image produced by existing applications. While useful for understanding some Mie scattering properties in specific cases, 2-D scattering plots alone are not sufficient for understanding the multidimensional structure of a scattered EM field. Furthermore, the computational time required by MiePlot [12], for computing scattering plots based on Mie theory (Figure 1A), increases proportionally to the radius of the sphere. The T-matrix software [20] requires 15 s for a full CPU evaluation of scattered fields (600 × 600 slice resolution) at a single wavelength for a sphere positioned at the focal point (Figure 1B).

FIGURE 1

Figure 1. Examples of Mie scattering functionality provided from available tools. (A) Polar plot of scattered intensity vs. scattering angle (λ = 1 μm, sphere radius = 1 μm) created using MiePlot [12]. (B) Scattered field produced by a PMMA sphere (radius = 1 μm, λ = 1 μm) positioned at the focal point. Image was created using a Python-based GUI for the parallel Multi-Sphere T-Matrix software [20].

2. Materials and Methods

We simulate an optical system similar to those used in most microscopes. The incident light is emitted from a source with an intensity spectrum given by I₀(k), where $k=\frac{2\pi }{\lambda }$ is the wavenumber.

Our simulation focuses on the scattering effects of spheres specified by a radius a and a complex refractive index n(k) = η(k)+iκ(k). Here, the real part, η(k), denotes the refractive index and indicates the phase velocity, while the imaginary part, κ(k), quantifies absorption. The total field resulting from scattering by a single sphere is given by:

$\begin{array}{rc}E(r,\theta )=\{\begin{array}{ll}{E}_i(r,\theta ),\hfill & \text{if }r<a\hfill \\ E_f(r,\theta )+E_s(r,\theta ),\hfill & \text{otherwise}\hfill \end{array}& (1)\end{array}$

where E_f is the incident field produced by the light source and focusing optics, E_s is the scattered field, E_i is the field inside of the sphere, (r, θ) are the spherical coordinates of each point in the field, and a is the radius of the sphere. The internal and external fields are equal at the sphere surface: E_i(r, θ) = E_f(r, θ)+E_s(r, θ) for r = a.

2.1. Calculation of the Focused Incident Field (E_f)

The vector formulation for the electric field is given by

$\begin{array}{rc}{E}_f(k,r)=E_0{e}^{i{k}^{\top }r}& (2)\end{array}$

subject to the constraint that E₀^⊤ K = 0. Here, E₀ is a vector providing the amplitude and polarization direction of the field, k is a vector representing the direction of the incident plane wave, and r = (x, y, z)^⊤ is the position vector.

2.1.1. Scalar Model

The focused field is represented using a superposition of plane waves:

$\begin{array}{ll}{E}_f(r)=E_0{\displaystyle \sum _{j=1}^J{E}_f}(k_j,r)=E_0{\displaystyle \sum _{j=1}^J{e}^{i{k}_j{}^{\top }r}}& (3)\end{array}$

We implement Debye focusing [21, 22] by first applying the partial wave expansion of the plane wave to obtain

$\begin{array}{ll}{e}^{i{k}^{\top }r}={\displaystyle \sum _{l=0}^{\infty }(2l+1)}{\text{i}}^l{j}_l(kr)P_l(\cos \theta ),& (4)\end{array}$

where l is the order of the incident field, j_l(·) is an order-l spherical Bessel function of the first kind, P_l(·) is an order-l Legendre polynomial, and k = ||k|| is the wavenumber. We use the j subscript to indicate the field produced by a single propagating plane wave k_j. We also assume that the incident light is coherent, therefore ||k_j|| = k for all j ∈ J.

Using the addition theorem of the spherical harmonics [23] it is possible to write the sum of the plane waves as the following integral

$\begin{array}{l}{E}_{\text{f}}(p)=2\pi E_0{\displaystyle \sum _{l=0}^{\infty }(2l+1)}{\text{i}}^l{j}_l(kr)P_l(\cos \theta )\\ {\displaystyle {\int }_{{\alpha }_1}^{{\alpha }_2}{P}_l}(\cos {\theta }_k)\sin ({\theta }_k)\text{d}{\theta }_k& (5)\end{array}$

where α₁ and α₂ indicate the angles subtended by the objective (Figure 2), θ_k is the angle between the unit position vector and the propagation direction of the incident plane wave: $\cos {\theta }_k=\frac{k^{\top }r}{\left|\right|k\left|\right|\left|\right|r\left|\right|}$, and p = (r, θ). For a lens, α₁ = 0 and ${\alpha }_2={sin}^{-1}(NA)$. In the case of a cassegrain mirror, which is commonly used for mid-infrared spectroscopic measurements, the inner angle α₁ will be the angle subtended by the central obscuration. If α₁ and α₂ are known, the incident field can be computed using a closed-form solution:

$\begin{array}{ll}{E}_f(p)=2\pi E_0{\displaystyle \sum _{l=0}^{\infty }{i}^l}{j}_l(kr)P_l(\cos \theta )c_l,& (6)\end{array}$

where c_l = [P_l+1(cosα₁) − P_l+1(cosα₂) − P_l−1(cosα₁)+P_l−1(cosα₂)]. The solution to Equation (6) for various optical parameters produces the point spread function (PSF) shown in Figure 2.

FIGURE 2

Figure 2. Constructive interference at the focal point f results in a high-intensity point spread function (PSF). The magnitude of the field resulting from Equation (6) is shown for a transparent lens (top) and a Schwarzschild objective (bottom). Note that increasing the size of the center obscuration creates a PSF with the characteristic features of a Bessel beam.

2.2. Calculation of the Scattered and Internal Fields (E_s and E_i)

The scattered field for a single incident plane-wave k produced by a sphere with radius a positioned at the focal point p_f is given by:

$\begin{array}{ll}{E}_s(p)={\displaystyle \sum _{l=0}^{\infty }{B}_l(\lambda ,n,a)h_l^{(1)}(kr)P_l(\cos \theta )}& (7)\end{array}$

where $h_l^{(1)}(x)$ is the order-l spherical Hankel function [8]. The coefficients B that couple the internal and external scattered fields are given by:

$\begin{array}{ll}B(\lambda ,n,a)=(2l+1)i^l\frac{j_l(ka)j_l^{\prime }(kna)n-j_l(kna)j_l^{\prime }(k)}{j_l(kna)h_l^{{(1)}^{\prime }}(ka)-h_l^{(1)}(ka)j_l^{\prime }(kna)n}& (8)\end{array}$

where $j_l^{\prime }(x)$ is the first derivative of the spherical Bessel function of the first kind and $h_l^{(1)\prime }(x)$ is the first derivative of the spherical Hankel function. These are derived by enforcing the appropriate boundary conditions [23].

Computing the scattered field for a condenser of NA_c > 0 requires integrating the scattered field equation across the solid angle subtended by the condenser lens. Since the only term in Equation (7) dependent on the light direction is the Legendre polynomial P_l(cosθ), this solution can be found analytically. However, this is only valid for spheres centered at p_f. For more generally positioned objects, each plane wave must be shifted by an additional phase term in the form of a complex exponential. The need for this phase shift is demonstrated in Figure 3B. A change in sphere position parallel to the direction of propagation of the plane wave will result in a phase delay. Applying this phase shift results in the final scattering equation:

$\begin{array}{ll}{E}_s(p)=e^{i{k}^{\top }c}{\displaystyle \sum _{l=0}^{\infty }{B}_l(\lambda ,n,a)h_l^{(1)}(kr)P_l(\cos \theta )},& (9)\end{array}$

where c = p_s − p_f is the vector from the focal point p_f to the center of the sphere p_s.

FIGURE 3

Figure 3. Reference points for scattering through a sphere. (A) The terms used in the scattering definitions are graphically defined. (B) A phase shift is included in the scattering equations for spheres that are not located at the focal point p_f. The purple line gives the advancing wave-front for a plane wave. Note that the green sphere will not require a phase shift, but the red sphere will require a shift equal to e^ik⊤_x.

The internal field, specifying the field inside of a sphere, for an incident plane wave k is given by:

$\begin{array}{ll}{E}_i(p)=e^{i{k}^{\top }c}{\displaystyle \sum _{l=0}^{\infty }{A}_l(\lambda ,n,a)j_l(knr)P_l(\cos \theta )}& (10)\end{array}$

with the scattering coefficients

$\begin{array}{ll}A(\lambda ,n,a)=(2l+1)i^l\frac{j_l(ka)h_l^{{(1)}^{\prime }}(ka)-j_l^{\prime }(ka)h_l^{(1)}(ka)}{j_l(kna)h_l^{{(1)}^{\prime }}(ka)-h_l^{(1)}(ka)j_l^{\prime }(kna)n}& (11)\end{array}$

Note that these equations have no known closed-form integrals. Therefore, the scattered field produced by a focused beam must be computed numerically.

3. Discussion

Using the proposed formulation (Equations 6, 9, and 10), we demonstrate interactive performance by minimizing the amount of computation and taking advantage of GPU-based hardware. GPU-based methods provide an inexpensive means of high-performance computing by taking advantage of parallelism and efficient allocation of resources. We are able to achieve interactive performance by minimizing the computational load and limiting the simulation domain to regions visualized by the user. In addition, fast evaluation is useful in situations where iterative evaluation of the forward model is required to solve an inverse problem [9, 11].

3.1. Focused Field

Note that the focused field E_f (Equation 6) depends only on position, wavelength λ, and the condenser NA_c. We first describe a concise representation of the focused field that allows fast evaluation of E_f as a function of position. We then show that this representation can be constructed quickly, thereby allowing interactive selection of the wavelength and condenser NA.

3.1.1. Representation

We compute and store the field by taking advantage of the symmetry provided using a spherical condenser. In this case, the condenser aperture subtends a circular solid angle that results in a circularly-symmetric incident field centered on the focal point. Note that this formulation is equally valid for a cassegrain with a central obscuration.

Because of this symmetry, the incident field between the condenser and objective can be reconstructed from a single cross-section of the cylindrical volume (Figure 4). This cross-section is stored as a two-dimensional (R x R) 32-bit floating-point array mapped to a 2-channel GPU-based texture. This texture is referenced in terms of cylindrical coordinates p ∈ [u v] with the origin located at p_f. In addition, it is only necessary to store $\frac{1}{4}$ of this cross-section by taking advantage of the symmetry along the cylinder axis v. Note that direct mirroring of the slice results in a phase shift that is corrected by swapping the sign of the imaginary component for v < 0 (Figures 4B,C).

FIGURE 4

Figure 4. The use of spherical optics to focus the incident light results in a circularly-symmetric incident field E_f. (A) A single slice of the incident field is evaluated and stored as a texture map. Any slice through the incident field is computed by interpolating values in this texture using the cylindrical coordinates (u, v). (B,C) Note that extrapolation across the v = 0 plane requires switching the sign of the imaginary component of E_f.

This representation has two major advantages. First of all, hardware-accelerated texture units provide linear interpolation between samples. This reduces the array resolution required to build an accurate approximation of the incident field. Since the simulation domain and resolution are user-specified, the accuracy of the simulation is tightly controlled by defining the desired array resolution R.

Secondly, GPU-based texture maps provide spatially coherent caching, which allows a block of processors to acquire values of E_f using a single fetch. Computing E_f at any point in the field then requires at most a single texture fetch, and on average $\frac{1}{N}$, where N is the GPU warp size (32 on an nVidia GeForce GTX 970). In addition, these texture fetches are performed in parallel to additional computation on newer systems.

3.1.2. Computation

The focused field is evaluated using a GPU-based kernel developed in CUDA. The Legendre polynomials dependent on condenser angle α are constant for all points in E_f, and are therefore pre-computed. The second Legendre polynomial, P_l(cosθ), depends on position and is computed independently for each point using a recursive definition, requiring only 4 floating point operations for each order l = [0 N_l], where N_l is the maximum order of the field (Equation 12).

The final component of E_f is the spherical Bessel function j_l(kr). These functions are time-consuming to compute. However, the only dependence is distance from the focal point p_f. Since the simulation domain is specified by the user, the parameter-space is one-dimensional and constrained by the user-specified simulation domain and resolution. These values are readily pre-computed and stored in an N_l x R table accessed as a 1D 2-channel texture, where R is the domain resolution (specified above). This insures at least one sample per pixel.

In conclusion, we show that a high-order representation of the entire field can be pre-computed and stored efficiently when the user makes a change to the condenser NA or incident wavelength. Since the focused field is independent of downstream properties, such as the position and material of particles, the value of E_f is determined at any point using a single texture fetch.

3.2. Scattered Field

The implicit functions for the internal field and external scattered field produced by a sphere are given in Equations (9) and (10). The scattered and internal fields are only rotationally symmetric when the position of the sphere is at the focal point (p_s = p_f) or when the incident light is represented by a single plane wave (${\text{E}}_s={\text{E}}_s^0$). This is due to the phase shift introduced by moving a sphere relative to p_f (Figure 3). We utilize the symmetry in the plane-wave solution to represent a particle's scattering properties in 2D. This image is computed and re-sampled using Monte-Carlo integration to estimate the internal and scattered fields produced by a particle positioned at any point in the near-field.

3.2.1. Evaluation for a Single Plane Wave

The internal and scattered fields are separated into four components:

• the scattering coefficients given by B and A

• a propagation function given by $h_l^{(1)}$ in E_s (or j_l in E_i)

• the Legendre polynomial P_l(cosθ)

• the phase shift exponential dependent on the light direction and sphere position.

The scattering coefficients are independent of position and therefore constant for any wavelength λ and material m. These values are therefore pre-computed for each scatterer. This dramatically improves performance, since each scattering coefficient (Equations 8, 11) requires computing multiple complex-valued Bessel functions and their derivatives.

The attenuation functions $h_l^{(1)}(kr)$ and j_l(knr) are dependent on distance from the sphere center. Like the propagation function for the focused field (j_l in Equation 6), both of these parameters are 1D and bounded by R. Therefore, these functions are pre-computed for a range of distance values. The Hankel function in E_s is independent of any material properties, and is stored with j_l(kr) for the focused field E_f (Section 3.1). Since the Hankel function can be expressed as a linear combination of Bessel functions of the first (j_l(kr)) and second kind (y_l(kr)), both j_l(kr) and y_l(kr) are stored in separate channels of the same texture, allowing both attenuation functions for E_f and E_s to be evaluated using a single texture fetch.

The Bessel function used to compute the internal field, however, is dependent on the index of refraction n(k) of the particle. Therefore, this table is stored for each sphere. The parameter range is significantly smaller, requiring only values d < a, where a is the radius of the sphere. As with the focused field, the Legendre polynomials are computed recursively.

These methods are used to compute the scattering domain (Figure 5) for a sphere. This a 2D representation of ${\text{E}}_s^0$ and ${\text{E}}_i^0$ as a function of d ∈ $\left[0\frac{FOV}{2}\right]$ and cos(θ) ∈ [−1 1]. The accuracy of the simulation is controlled by R × R_θ, where R is the spatial resolution of the field slice and R_θ is the angular resolution of the scattered field emanating from the sphere surface both inward and outward. The required R_θ is directly dependent on the highest order Legendre polynomial, which has l−1 oscillations in the domain of cos(θ). Unlike E_f, the internal and scattered fields are scaled by the scattering coefficients B and A, which quickly converge to zero with increasing values of l. The maximum order N_l required for convergence is given by Bohren and Huffman [8]:

$\begin{array}{ll}{N}_l=\lceil \frac{2\pi a}{\lambda }+4\sqrt[3]{\frac{2\pi a}{\lambda }}+2\rceil & (12)\end{array}$

Application of Nyquist sampling implies that R_θ ≥ 2N_l is required to capture all of the oscillations, while greater values increase accuracy for points far from the sphere. Our simulations use R_θ = 1000.

FIGURE 5

Figure 5. (A) The pre-computed scattering domain ${\text{E}}_s^0$ for a sphere is a 2D function dependent on the position of a point p in terms of its distance from the sphere center, d, and orientation relative to an incoming plane wave. Color represents field magnitude. (B) Surfaces in the near field are shown mapped to the scattering domain. Values where d < r (orange) use the separate ${\text{E}}_i^0$ domain. Monte-Carlo samples are always collected along the y = cosθ direction near the curve, within a range of ±cosα as specified by the vertical bars.

3.2.2. Monte-Carlo Integration

The solutions for E_i and E_s are determined using Monte-Carlo integration based on a uniform distribution of sample points within the solid angle α defined by NA_c. The source points for plane waves are determined using a stratified uniform distribution projected onto the unit sphere based on Archimedes' principle [24]. A uniform distribution of points is computed in cylindrical coordinates in the range ϕ = [0 2π], z = [cosα 1] using jittered stratified sampling [25] and projected inward onto the unit sphere (Figure 6).

FIGURE 6

Figure 6. Monte-Carlo sampling of the solid angle α defined by the condenser NA_c. (A) The internal and scattered fields (E_i and E_s) are determined by integrating the results based on plane waves originating at points uniformly distributed on the unit sphere and bounded by the condenser aperture. (B) Uniform sampling is performed using Archimedes' principle by creating a uniform distribution on a cylinder in the range of ϕ = [0 2π], z = [1 cos(α)] and projecting inward onto the unit sphere.

The scattered and internal fields produced by a particle are approximated by:

$\begin{array}{l}{E}_s\approx \frac{2\pi E_0((1-\cos {\alpha }_2)-(1-\cos {\alpha }_1))}{M}{\displaystyle \sum _{j=1}^M{E}_s(p,k_j)}\\ E_i\approx \frac{2\pi E_0((1-\cos {\alpha }_2)-(1-\cos {\alpha }_1))}{M}{\displaystyle \sum _{j=1}^M{E}_i(p,k_j)},& (13)\end{array}$

where M is the number of Monte-Carlo samples, k_j is a single plane-wave direction given by a random sample, and E_s, E_i are the equations for a scattered and internal field produced by a single plane wave (Equations 9, 10).

3.3. Surface Fields

The scattered field at the sphere surface is an important characteristic, and often the subject of simulations using Mie theory. In addition, both the internal and external fields are identical at this interface. This is generally referred to as the scattering efficiency, and is often represented as a graph plot, given as a function of θ. However, it is difficult to characterize the scattering efficiency for a sphere when p_s ≠ p_f, since this introduces an additional ϕ-dependence that is time-consuming to simulate and difficult to visualize using a 1D plot.

We visualize the field using a geometric surface. This technique is similar to those applied to other spherical functions in bio-medical imaging [26], such as diffusion-tensor MRI [27], and provides a unique insight into the structure of the scattered field. This allows the user to explore the scattering efficiency of a particle in three-dimensions as a function of material and position (Figure 9) and validate the effectiveness of MC integration (Figure 7).

FIGURE 7

Figure 7. Monte-Carlo sampling of the scattered field at the surface of a 3 μm particle in a λ = 2.5 μm field focused by a 0.8 NA condenser. The error in Monte-Carlo integration falls off with $\frac{1}{\sqrt{N}}$, where N is the number of samples. Due to the constructive coherence in a focused field, the relative error is negligible near the focal point for N ≥ 400. The relative error will increase significantly with d = |p − p_f| since the signal decreases with the square of the distance. (A) Scattered field from a simulation using a single plane wave, (B) using 4 Monte-Carlo samples, (C) using 400 Monte-Carlo samples, (D) using 400 Monte-Carlo samples.

3.4. Simulated Imaging

The final step in the imaging process is to determine the field at the detector and produce an image. The field is focused onto a detector to produce an image of the objective focal plane. Note that the image plane contains the focal point of the objective lens, p_o.

3.4.1. Objective and Detector

According to the principles of Fourier optics, the objective aperture blocks high- and low-frequency components of the image plane given by the upper and lower cutoff frequencies: $f_u=\frac{N{A}_o}{\lambda },f_l=\frac{N{A}_{in}}{\lambda }$.

The resulting image is produced by evaluating a slice f(x, y) of the field at p_o. The user specifies the field-of-view (FOV) S, and field resolution R. We then perform a forward Fourier Transform (FT) of the field slice, resulting in a frequency-domain representation F(u, v) = ${F}$[f(x, y)]. Frequency components above the cutoff frequency f_u and below the cutoff frequency f_l are eliminated. The field at the detector is therefore $\widehat{f}(x,y)={{F}}^{-1}\left[F(u,v)A_p(u,v)\right]$, where the aperture function is given by

$\begin{array}{ll}{A}_p(u,v)=\{\begin{array}{ll}0,\hfill & \text{if}\sqrt{{(u\Delta u)}^2+{(v\Delta v)}^2}>f_u\text{or}<f_l\hfill \\ 1,\hfill & \text{otherwise}\hfill \end{array}& (14)\end{array}$

and $\Delta u=\Delta v=\frac{1}{S}$. The detector then transforms the complex field into the final intensity image $I=|\widehat{f}(x,y){|}^2$.

The forward and inverse Fourier transforms are performed using the GPU-based cuFFT Fast Fourier Transform (FFT) [28]. The resulting FFT and detector images for an array of particles with varying extinction coefficients κ are shown in Figure 8.

FIGURE 8

Figure 8. Objective affect on an image of 3 μm particles with n = 1.4+κ, where κ = 0.0 → 0.12. The incident field is λ = 2.5 μm and the field of view is 100 μm. (A) Fourier transform of the EM field at the image plane. (B,C) Intensity images produced at the detector for 1.0 and 0.6 NA objectives.

3.4.2. Absorbance

Common vibrational spectroscopic imaging techniques, such as mid-infrared spectroscopic imaging, rely on absorbance measurements of a tissue sample in order to estimate the extinction κ. The absorption is given by $A=-lo{g}_{10}\left(\frac{I}{I_0}\right)$, where I₀ is an image of the incident field, without any particles present.

A simulated image of multiple materials is shown in Figure 10, demonstrating the complex interactions between the incident field, particles, and optics in order to produce a final absorption image.

3.4.3. Extended Sources

We also allow the software to simulate extended (non-coherent) sources. Extended sources are created by simulating multiple point sources and integrating their intensity with the detector. An image of the extended source can be specified, allowing different sources (point, Gaussian, glowbar) to be simulated.

4. Results

We have developed a software package called the Broadband Interactive Mie Simulator (BIM-Sim), which interactively computes the Bjorn approximation for scattered fields produced by spherical particles. This includes the simulation of a focusing lens, objective lens, and imaging system. The complete field is computed both inside and outside of each particle. BIM-Sim is open-source and available online at http://stim.ee.uh.edu/resources/software/bimsim/.

We visualize the field using planar cross-sections of the volume placed anywhere between the objective and condenser (Figure 9). The complex components of the scattered field at the particle surface are computed and visualized using a 3D surface model (Figure 9). This facilitates the study of the scattering efficiency of a sphere in three dimensions. The shape of the scattered field is computed using Monte-Carlo sampling, which allows arbitrary positioning of the particle within the focal volume. The use of Monte-Carlo sampling also allows any apodization function to be used for the focusing optics, facilitating the study of beam quality on the resulting field [29, 30].

FIGURE 9

Figure 9. Scattered fields created by three 2 μm diameter spheres separated by 2 μm of vacuum in a focused EM beam produced by a 0.2 NA condenser. All three spheres have identical material properties. The magnitude and the real part of the scattered field at the sphere surface are shown (top). Note that the sphere positioned within the center of the focused beam produces radially symmetric scattering, while adjacent spheres show reduced intensity and an asymmetric scattered field. A cross section of the full field (Equation 1) through the spheres is shown (bottom).

The imaging process is simulated by band-limiting a cross-section of the field at the focal plane based on the objective aperture function. The intensity of the filtered cross-section is computed and re-sampled based on user-specified detector parameters to create a final image (Figures 10, 11). Computing both the total and incident field intensity also allows the simulation of absorption spectroscopy for any distribution of micro-spheres. Material properties for the spheres are specified at run-time or as a wavelength-dependent set of refractive indices. If the extinction spectrum is known, BIM-Sim uses the Kramers-Kronig relation to determine the phase speed [31] as a function of λ.

FIGURE 10

Figure 10. Absorption image and scattering efficiency of three spheres made of different materials with absorption spectra in the mid-infrared. (A) The absorption image (left) and scattering efficiency (right) is shown at λ = 3.1 μm. The real and imaginary parts of the index of refraction are shown for toluene. (B) The absorption image and scattering efficiency for λ = 6.6 μm is shown (top) with the index of refraction for water (H₂O, bottom).

FIGURE 11

Figure 11. Scattered fields and absorption image created by three 2 μm spheres made of noble metals: silver, gold, and platinum (from left to right). A cross section of the magnitude (|E|) and the real part [real(E)] of the full field through the spheres, generated by a single plane wave, are shown (top left). The focused EM beam (E_f using a field order of 100) and the absorption image at λ = 3.1 μm are shown (top right). The real and imaginary parts of the index of refraction are shown for silver, gold, and platinum (bottom).

4.1. Comparisons with Measured Data

We compared the predicted absorbance data (spectra and images) from the forward model described above with IR absorbance data of poly(methyl methacrylate) (PMMA) and polystyrene microspheres. PMMA microspheres were obtained from Bangs Laboratories, Fishers, IN, with diameters between 1 and 6.5 μm. Polystyrene microspheres were obtained from Bangs Laboratories, Fishers, IN, and Polysciences, Warrington, PA, with diameters between 4 and 6 μm. A small volume of spheres was dispersed onto a 3 mm thick calcium fluoride (CaF₂) substrate by gently tapping a loaded pipette tip, and regions with individual spheres on the substrate were found. Infrared imaging data were recorded using a Cary 670 Series FTIR Spectrometer coupled to a Cary 620 Series IR Microscope (Agilent Technologies, Santa Clara, CA). The system was equipped with a 128 × 128 pixel focal plane array (FPA) detector that was used for imaging. The NA of the detection Schwarzschild objective and condenser was 0.62 with an obscuration of NA = 0.34. The imaging was performed in high magnification mode (pixel size of 1.1 μm). Spectra were recorded at 16 cm⁻¹ resolution.

We estimated the complex refractive index as a function of wavelength from bulk measurements of PMMA and polystyrene absorption spectra. The effective refractive index, n, and the absorption coefficient, k, obey the Kramer's-Kronig relations, allowing one to be estimated using a measurement of the other. In this case, these values are computed using the following equations:

$\begin{array}{l}A(\lambda )\approx -\frac{2\pi }{n_0\lambda }{\chi }^{″}(\lambda )n(\lambda )\approx n_0+\frac{{\chi }^{\prime }(\lambda )}{2n_0},\end{array}$

where A is the absorption (per unit thickness), n₀ is the baseline refractive index (≈ 1.4 for most polymers). The real and imaginary components of the magnetic susceptibility

$\begin{array}{l}\chi (\lambda )={\chi }^{\prime }(\lambda )+i{\chi }^{″}(\lambda )\end{array}$

are related through the Hilbert transform, and the absorption coefficient

$\begin{array}{l}k=-\frac{A(\lambda )\lambda }{4\pi }\end{array}$

can be computed from the absorption. The complex refractive index for PMMA and polystyrene are shown in Figure 12.

FIGURE 12

Figure 12. Complex refractive index for PMMA (left) and polystyrene (right), computed from absorption measurements collected using FTIR spectroscopy.

Spectra and absorbance images from a number of PMMA and polystyrene spheres were recorded. The forward model is validated by comparing the predicted absorption spectra of any pixel from the simulated spheres to the actual measured spectra using the FTIR imaging system. The predicted spectra from the described forward model and the measured spectra for PMMA and polystyrene microspheres are shown in Figure 13. The results show a close match between the measured and the predicted spectra. We also compare the simulated absorption images at any wavelength with measured absorbance images (Figures 14, 15). It can be observed that the predicted and measured absorption images match closely up to the noise level.

FIGURE 13

Figure 13. Measured and predicted absorption spectra for PMMA (top) and polystyrene (bottom) microspheres. Spheres with a radius of 1.5, 2.2, and 2.5 μm are considered. The spectrum from the center pixel of each sphere is shown.

FIGURE 14

Figure 14. Measured and predicted absorption images for a PMMA microsphere of radius ≈2.2 μm at λ = 9.4, 8.2, 7.3 μm.

FIGURE 15

Figure 15. Measured and predicted absorption images for a polystyrene microsphere of radius ≈2.2 μm at λ = 9.4, 8.2, 7.3 μm.

There are a number of factors which can effect the measured spectrum of polymer microspheres, such as the noise level, position of the sphere in the FOV, the material properties of the sphere, the focal point, the spherical shape of the microsphere etc. We note that the signal-to-noise ratio for the measured data (recorded in high magnification mode) is reasonable but not exceptional while the forward model described here does not incorporate noise. The results shown here are based on running our forward model with the following assumptions about each imaged microsphere: it is positioned at the center of the FOV, it is composed of clear material (in this case PMMA and polystyrene), has a perfect spherical shape, and the microsphere is on focus. We tried to image the microspheres as close as possible to the center of the FOV. However, it was not possible to determine the exact position of the center pixel of each imaged sphere in the FOV using our imaging system (a small offset of the sphere from the center of the image can be observed in Figures 14, 15). The polystyrene microspheres where composed of additional ingredients, such as water and 2.6% latex. The commercially available microspheres are not guaranteed to have a perfect spherical shape and they are of varying diameter size within the same particle solution.

In order to validate our results for predicted spectra from different pixels in the sphere we compare a cluster of measured and predicted pixel spectra around the center of the spheres. Figure 16 shows that the predicted spectra from different pixel locations closely follow the measured spectra. Also, notice that there is more variation within measured spectra and more overlap within the predicted ones. Overall, our results show that the described theory for the forward model is capable of describing recorded data with commercially available microscopy systems.

FIGURE 16

Figure 16. Measured and predicted absorbance spectra for PMMA (top) and polystyrene (bottom) microspheres. Spheres with a radius of 1.5, 2.2, and 2.5 μm are considered. A cluster of spectra around the center pixel of each sphere is shown.

4.2. Accuracy and Timing Details

The most time-consuming operations are changes to the incident wavelength λ and the condenser NA, which require re-evaluation of the focused field and all scattered fields. Our technique allows the EM field to be computed to user-specified accuracy, subject to the limits of numerical precision, since the user can select the resolution of the simulation. We use 32-bit floating point precision.

Stratified sampling is used to limit variance across the condenser aperture [24, 25]. Sample points are reconstructed randomly when NA_c is changed. However, the same random seed value is used, making random samples deterministic in order to reduce flickering as these changes are reflected in the final image.

All of our results were tested on an nVIDIA Geforce GTX 970 with 4 GB of memory and 13 multiprocessors. For a 512 × 512 resolution slice, full evaluation of E_f requires 2–4 ms and evaluation of a single scattered field requires 400–500 ms. This includes both computation of the scatter domain and 400 Monte-Carlo samples of a single sphere. Evaluation time scales linearly with the number of particles. Our current simulations are bounded by the amount of time required to fetch multiple samples of the scattered domain for Monte-Carlo sampling (one fetch per sample).

5. Conclusions

In this paper, we address the need for interactive simulation and visualization of scattered EM fields. This is an important problem in the optics community, where there is significant interest in understanding how micro and nano-particles interact with incident radiation. Our two major contributions include a method for interactive simulation of scattered fields and a framework for visualizing these fields by coupling user interaction with sparse simulation.

The simulation methods that we describe are interactive, allowing exploration of scattered fields that would previously have been impractical. The resulting speedup allows the user to dynamically set all aspects of the simulation. This has many advantages over pre-computation. In particular, a pre-computed representation of a field at the demonstrated magnitude and dimension would require an impractical amount of storage space and would be difficult to query efficiently. In addition, a fast forward model of particle scattering provides a framework for solving inverse problems, such as spectral un-mixing [9] and localization of probes within an imaged field.

Author Contributions

SB and DM designed the GPU-based algorithms proposed in this paper. Tv and PC designed the analytical model and phase corrections. RB and PC designed the theoretical framework behind mid-infrared spectroscopic scattering models. All authors contributed to writing this paper.

Funding

This work has been supported in part by the Cancer Prevention and Research Institute of Texas (CPRIT) #RR140013, National Institutes of Health/National Library of Medicine #4 R00 LM011390-02, and Agilent Technologies University Relations #3938.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

The authors would like to thank Brad Deutsch and Rohith Reddy for their insights on analytical modeling of electromagnetic fields and Suzanne Leslie for her helpful tutorial on the relationship of measured absorbance spectra A(λ) to κ(λ).

References

1. Bond TC, Habib G, Bergstrom RW. Limitations in the enhancement of visible light absorption due to mixing state. J Geophys Res. (2006) 111:D20211. doi: 10.1029/2006JD007315

CrossRef Full Text | Google Scholar

2. Bassan P, Kohler A, Martens H, Lee J, Jackson E, Lockyer N, et al. RMieS-EMSC correction for infrared spectra of biological cells: extension using full Mie theory and GPU computing. J Biophotonics (2010) 3:609–20. doi: 10.1002/jbio.201000036

PubMed Abstract | CrossRef Full Text | Google Scholar

3. West JL, Halas NJ. Engineered nanomaterials for biophotonics applications: improving sensing, imaging, and therapeutics. Annu Rev Biomed Eng. (2003) 5:285–92. doi: 10.1146/annurev.bioeng.5.011303.120723

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Xu HX. A new method by extending Mie theory to calculate local field in outside/inside of aggregates of arbitrary spheres. Phys Lett A (2003) 312:411–9. doi: 10.1016/S0375-9601(03)00687-X

CrossRef Full Text | Google Scholar

5. Li Z, Xu H. Electromagnetic energy flow near metal nanoparticlesII: algorithms for the calculation of the light scattering of multi-spheres and photon energy transport via linear chains of Ag nanoparticles. J Quant Spectrosc Radiat Transf. (2007) 103:394–401. doi: 10.1016/j.jqsrt.2006.06.007

CrossRef Full Text | Google Scholar

6. Wang H, Li Z, Zhang H, Wang P, Wen S. Giant local circular dichroism within an asymmetric plasmonic nanoparticle trimer. Sci Rep. (2015) 5:8207. doi: 10.1038/srep08207

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Wang H, Miao L, Jiang Y, Lu S, Li Z, Li P, et al. Enhancing the saturable absorption and carrier dynamics of graphene with plasmonic nanowires. Phys Status Solidi B (2015) 252:2159–66. doi: 10.1002/pssb.201570362

CrossRef Full Text | Google Scholar

8. Bohren CF, Huffman DR. Absorption and Scattering of Light by Small Particles. Hoboken, NJ: John Wiley & Sons (2008).

Google Scholar

9. van Dijk T, Mayerich D, Carney PS, Bhargava R. Recovery of absorption spectra from Fourier transform infrared (FT-IR) microspectroscopic measurements of intact spheres. Appl Spectrosc. (2013) 67:546–52. doi: 10.1366/12-06847

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Mayerich D, Dijk Tv, Walsh MJ, Schulmerich MV, Carney PS, Bhargava R. On the importance of image formation optics in the design of infrared spectroscopic imaging systems. Analyst (2014) 139:4031–6. doi: 10.1039/C3AN01687K

PubMed Abstract | CrossRef Full Text

11. van Dijk T, Mayerich D, Bhargava R, Carney PS. Rapid spectral-domain localization. Opt Express (2013) 21:12822–30. doi: 10.1364/OE.21.012822

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Laven P. A Computer Program for Scattering of Light From a Sphere Using Mie Theory & the Debye Series (2011). Available online at: www.philiplaven.com (Accessed March 5, 2014).

13. Cai W, Kranendonk L, Lee T, Ma L. Characterization of composite nanoparticles using an improved light scattering program for coated spheres. Comput Phys Commun. (2010) 181:978–84. doi: 10.1016/j.cpc.2010.01.010

CrossRef Full Text | Google Scholar

14. Wu ZS, Guo LX, Ren KF, Gouesbet G, Grehan G. Improved algorithm for electromagnetic scattering of plane waves and shaped beams by multilayered spheres. Appl Opt. (1997) 36:5188–98.

PubMed Abstract | Google Scholar

15. Mishchenko MI, Hovenier JW, Travis LD. Light Scattering by Nonspherical Particles: Theory, Measurements, and Applications. Cambridge, MA: Academic Press (1999).

Google Scholar

16. Davis BJ, Carney PS, Bhargava R. Theory of mid-infrared absorption microspectroscopy: II. heterogeneous samples. Anal Chem. (2010) 82:3487–99. doi: 10.1021/ac902068e

PubMed Abstract | CrossRef Full Text | Google Scholar

17. van de Hulst HC. Light Scattering by Small Particles. Mineola, NY: Dover Publications (2012).

18. Matzler C. MATLAB Functions for Mie Scattering and Absorption, Version 2. IAP Research Report (2002). Available online at: http://www.atmo.arizona.edu/students/courselinks/spring08/atmo336s1/courses/spring09/atmo656b/maetzler_mie_v2.pdf

19. Davis BJ, Carney PS, Bhargava R. Theory of midinfrared absorption microspectroscopy: I. homogeneous samples. Anal Chem. (2010) 82:3474–86. doi: 10.1021/ac902067p

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Mackowski DW, Mishchenko MI. A multiple sphere T-matrix Fortran code for use on parallel computer clusters. J Quant Spectrosc Radiat Trans. (2011) 112:2182–92. doi: 10.1016/j.jqsrt.2011.02.019

CrossRef Full Text | Google Scholar

21. Debye P. Das verhalten von Lichtwellen in der nahe eines brennpunktes oder einer brennlinie. Annln Phys. (1909) 335:755–76.

Google Scholar

22. Wolf E, Li Y. Conditions for the validity of the Debye integral representation of focused fields. Opt Commun. (1981) 39:205–10.

Google Scholar

23. Joachain CJ. Quantum Collision Theory. North-Holland: North-Holland Publication (1975).

Google Scholar

24. Shao MZ, Badler N. Spherical Sampling by Archimedes' Theorem. University of Pennsylvania Department of Computer and Information Science Technical Report, Vol. MS-CIS-96-02, (1996). Available online at: http://repository.upenn.edu/cis_reports/184/?utm_source=repository.upenn.edu%2Fcis_reports%2F184&utm_medium=PDF&utm_campaign=PDFCoverPages

25. Green R. Spherical harmonic lighting: the gritty details. In: Archives of the Game Developers Conference, Vol. 56 (2003). Available online at: http://www.cs.columbia.edu/~cs4162/slides/spherical-harmonic-lighting.pdf

26. Ropinski T, Oeltze S, Preim B. Survey of glyph-based visualization techniques for spatial multivariate medical data. Comput Graph. (2011) 35:392–401. doi: 10.1016/j.cag.2011.01.011

CrossRef Full Text | Google Scholar

27. van Almsick M, Peeters THJM, Prckovska V, Vilanova A, and ter Haar Romeny B. GPU-based ray-casting of spherical functions applied to high angular resolution diffusion imaging. IEEE Trans Vis Comput Graph. (2011) 17:612–25. doi: 10.1109/TVCG.2010.61

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Lee D, Dinov I, Dong B, Gutman B, Yanovsky I, Toga AW. CUDA optimization strategies for compute- and memory-bound neuroimaging algorithms. Comput Methods Prog Biomed. (2012) 106:175–87. doi: 10.1016/j.cmpb.2010.10.013

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Carney PS, Gbur G. Optimal apodizations for finite apertures. J Opt Soc Am A (1999) 16:1638.

Google Scholar

30. Gbur G, Carney PS. Convergence criteria and optimization techniques for beam moments. Pure Appl Opt. (1998) 7:1221.

Google Scholar

31. Kuzmenko AB. KramersKronig constrained variational analysis of optical spectra. Rev Sci Instrum. (2005) 76:083108. doi: 10.1063/1.1979470

CrossRef Full Text | Google Scholar