Optical coherence tomography images simulated with an analytical solution of Maxwell’s equations for cylinder scattering

Thomas Brenner; Dominik Reitzle; Alwin Kienle

doi:10.1117/1.JBO.21.4.045001

1 April 2016 Optical coherence tomography images simulated with an analytical solution of Maxwell’s equations for cylinder scattering

Thomas Brenner, Dominik Reitzle, Alwin Kienle

Author Affiliations +

Journal of Biomedical Optics, Vol. 21, Issue 4, 045001 (April 2016). https://doi.org/10.1117/1.JBO.21.4.045001

Abstract

An algorithm for the simulation of image formation in Fourier domain optical coherence tomography (OCT) for an infinitely long cylinder is presented. The analytical solution of Maxwell’s equations for light scattering by a single cylinder is employed for the case of perpendicular incidence to calculate OCT images. The A-scans and the time-resolved scattered intensities are compared to geometrical optics results calculated with a ray tracing approach. The reflection peaks, including the whispering gallery modes, are identified. Additionally, the Debye series expansion is employed to identify single peaks in the OCT A-scans. Furthermore, a Gaussian beam is implemented in order to simulate lateral scanning over the cylinder for two-dimensional B-scans. The fields are integrated over a certain angular range to simulate a detection aperture. In addition, the solution for light scattering by layered cylinders is employed and the various layers are identified in the resulting OCT image. Overall, the simulations in this work show that OCT images do not always display the real surface of investigated samples.

1. Introduction

Optical coherence tomography (OCT) is an emerging tool in clinical diagnostics, which enables high-resolution imaging of the internal structure of biological tissue.¹^–³ It is, for example, used to detect glaucoma in ophthalmology and coronary intervention can be facilitated in cardiology. In a clinical environment, the processed interferograms are interpreted based on experience because the link between the microscopic properties of the scatterers influencing light propagation and the generated images is not fully understood.⁴^–⁷ Therefore, the formation of images in OCT is investigated with computational models. Depending on the accuracy needed for simulated images, different OCT models are employed in literature. Many OCT simulations use Monte Carlo approaches for the radiative transfer equation⁸^–¹³ or they use a numerical solution of Maxwell’s equations without accounting for the interferometer setup.¹⁴^,¹⁵ Several proceedings dealing with Maxwell’s equations and OCT exist as well.¹⁶^–¹⁸ In general, the use of the radiative transfer equation in a Monte Carlo approach is an approximation that neglects interference phenomena. Additionally, omitting the interferometer setup changes the properties of the resulting images. To the knowledge of the authors, only recently a full-wave approach including the interferometer setup has been implemented for simulation of image formation in OCT.¹⁹^,²⁰

In this work, the analytical solutions of Maxwell’s equations for scattering of a plane wave and of a two-dimensional (2-D) Gaussian beam by an infinitely long dielectric homogeneous cylinder have been used to implement an interferometer setup with broadband illumination. The plane wave solution generates one-dimensional A-scans while the Gaussian beam can be scanned laterally over the cylinder to generate B-scans. To the knowledge of the authors, the analytical Maxwell solution for cylinder scattering has not yet been used for the simulation of OCT images. Additionally, the time-resolved fields scattered by a cylinder are calculated so that the intensities correspond to measurements with an ultrafast detector without an interferometer setup. For comparison of the time-resolved scattered intensity with geometrical optics, a 2-D ray tracing algorithm based on Snell’s law has been implemented for plane wave illumination. The simulated images based on Maxwell’s equations show the whispering gallery modes (WGMs) that appear in scattering by a cylinder. There is an approach for simulating OCT images with a solution of Maxwell’s equations for spherical scatterers,²¹ but the appearing WGMs have not been investigated further in their work and the single signals have not been associated with specific pathways. In this work, the single pathways are labeled and the behavior of the surface waves is explained based on the simulated OCT images. Beyond the geometrical optics simulation, we employ the Debye series to identify the OCT peaks in the Maxwell solution. To the knowledge of the authors, the Debye series has not been used for the simulation of OCT scans before. The basic understanding of the effects of a single scatterer with a smooth surface is the first step toward understanding more complicated scattering processes that influence the features of OCT images. Based on this knowledge, it will be possible to investigate more complicated scatterers in the future.

1.1.

Theory

OCT uses a broadband interferometer setup to reconstruct the internal structure of semitransparent scatterers. The changes in refractive index and polarization affect the backscattered light so that the correlation between the reference beam and the probing beam contains depth information about the scatterer. Either the length of the reference arm is shifted so that an axial profile is recorded (time domain OCT) or the whole information is recorded in the Fourier domain since the depth profile of the scatterer is also contained in the frequency spectrum according to the Wiener–Khintchine theorem. As shown in Fig. 1, a Fourier-domain OCT setup is considered in this work. It is assumed that the incident pulse in time domain has the form $E (t) = E_{0} \sqrt{2} e^{- \frac{t^{2}}{b^{2}}} \cos (ω_{0} t)$ so that the beam splitter with a splitting ratio of $\frac{1}{2} ∶ \frac{1}{2}$ yields

Eq. (1)

E_{r} (t) = E_{i} (t) = E_{0} e^{- \frac{t^{2}}{b^{2}}} \cos (ω_{0} t)

in the reference arm and the sample arm, respectively.

E_{r}

is the field in the reference arm and

E_{i}

is the field in the sample arm incident on the scatterer.

E_{0}

is the field amplitude,

ω_{0}

is the central frequency, and

2 \sqrt{\log 2} b

is the amplitude-related FWHM of the envelope of the pulse in time domain. The Fourier transform yields

Eq. (2)

E_{r} (ω) = E_{i} (ω) = E_{0} \frac{b}{2 \sqrt{2}} (e^{{- (\frac{b ω + b ω_{0}}{2})}^{2}} + e^{{- (\frac{b ω - b ω_{0}}{2})}^{2}}) .

Fig. 1

Simplified scheme of a Fourier-domain OCT setup. A broadband light source sends a light signal that is split at the beam splitter into a reference and a probe beam so that the scattered wave from the sample arm interferes with the signal from the reference arm at the detector. The interferogram is measured with a spectrometer and processed with the inverse Fourier transform to obtain temporal data. In the simulation, $r_{1} = d_{1} + d_{2}$ denotes the distance of the light source to the cylinder surface. $r_{2} = d_{2} + d_{4}$ is the distance the scattered light travels from the cylinder surface to the detector and $r_{3} = d_{1} + 2 d_{3} + d_{4}$ is the distance the reference pulse travels in total.

For the spectrum needed, we only consider positive frequencies and in order to keep the original peak height, the half-sided spectrum is multiplied by a factor of 2 and by the phase corresponding to the distance the pulse travels in each case (compare with the caption in Fig. 1):

Eq. (3)

E_{i} (ω) = 2 E_{0} \frac{b}{2 \sqrt{2}} (e^{- {(\frac{b ω - b ω_{0}}{2})}^{2}}) e^{i k r_{1}},

Eq. (4)

E_{r} (ω) = 2 E_{0} \frac{b}{2 \sqrt{2}} (e^{- {(\frac{b ω - b ω_{0}}{2})}^{2}}) e^{i k r_{3}},

where

r_{1}

is the distance from the light source to the cylinder surface and

r_{3}

is the total distance the reference pulse travels. Often, the central wavelength

λ

and the full-width half-maximum (FWHM)

Δ λ

of the source pulse are the given parameters for OCT devices. The central frequency is related to the central wavelength over

ω_{0} = (2 π c / λ)

with

c

being the speed of light in vacuum. The FWHM in frequency is related to the FWHM in

λ

as

Δ ω \approx (2 π c Δ λ / λ^{2})

. From the spectral width follows the temporal FWHM as

b = (4 \sqrt{\ln 2} / Δ ω)

for the field amplitude. The scattered field is calculated for

N

spectral values of the pulse ranging from

ω_{\min}

to

ω_{\max}

. The first frequency is set to

ω_{\min} = 0

. The last frequency value is the frequency belonging to the cutoff amplitude, where

E_{i} (ω_{\max}) = ε (b / 2 \sqrt{2}) E_{0}

. This yields

Eq. (5)

ω_{\max} = ω_{0} - \frac{2 \ln ε}{\sqrt{- b^{2} \ln ε}} .

The spectral values used for further calculation are, therefore, equally spaced

ω = (k / c) n_{1}

.

ε

is chosen to be

10^{- 30}

. Given the spectrum, the next step is to calculate the scattered field and the reference field. The scattered field is the solution of the Helmholtz equation for a homogeneous, dielectric cylinder. The detected light in the far field for a plane wave with perpendicular incidence

E_{i} (ω)

scattered by the cylinder in the sample arm is given by Bohren and Huffman²² as

Eq. (6)

(\begin{matrix} E_{s ∥} (r, θ, k) \\ E_{s ⊥} (r, θ, k) \end{matrix}) = e^{i \frac{3}{4} π} \sqrt{\frac{2}{π k (r_{2} + a)}} e^{i k (r_{2} + a)} (\begin{matrix} T_{1} (θ) & T_{4} (θ) \\ T_{3} (θ) & T_{2} (θ) \end{matrix}) (\begin{matrix} E_{i ∥} (k) \\ E_{i ⊥} (k) \end{matrix}),

∥

is the TM mode with polarization parallel to the cylinder axis and

⊥

is the TE mode with polarization perpendicular to the cylinder axis. The distance between the point detector and the origin at the cylinder center is

r_{2} + a

, where

a

is the cylinder radius.

r_{2}

denotes the distance the scattered wave travels from the cylinder surface to the detector. The

T

-matrix elements can be calculated according to²²

Eq. (7)

T_{1} (θ) = \sum_{n = - \infty}^{\infty} b_{n ∥} e^{- i n θ} = b_{0 ∥} + 2 \sum_{n = 1}^{\infty} b_{n ∥} \cos (n θ),

Eq. (8)

T_{2} (θ) = \sum_{n = - \infty}^{\infty} a_{n ⊥} e^{- i n θ} = a_{0 ⊥} + 2 \sum_{n = 1}^{\infty} a_{n ⊥} \cos (n θ),

Eq. (9)

T_{3} (θ) = \sum_{n = - \infty}^{\infty} a_{n ∥} e^{- i n θ} = - 2 i \sum_{n = 1}^{\infty} a_{n ∥} \sin (n θ),

Eq. (10)

T_{4} (θ) = \sum_{n = - \infty}^{\infty} b_{n ⊥} e^{- i n θ} = - 2 i \sum_{n = 1}^{\infty} b_{n ⊥} \sin (n θ) = - T_{3} (θ) .

The truncation criterion for the infinite sum and the shape of the expansion coefficients

a_{n}

and

b_{n}

as used in this work can be found in Ref. 22. The orientation of the coordinate system for the scattering problem is shown in Fig. 2.

z

is the axial coordinate and

y

is the lateral coordinate. We refer to the cylinder surface facing the illumination source (e.g., a plane wave coming from negative

z

, as shown in the figure) as the cylinder front side and the cylinder surface on the opposite side at positive

z

-values as the cylinder back side. The scatterer is symmetric with respect to the

z

-axis in Fig. 2. This symmetry axis is referred to as the optical axis.

Fig. 2

Orientation of the coordinate system for the scattering problem. The origin of the coordinate system is at the center of the cylinder cross section. $n_{1}$ is the refractive index of the surrounding medium and $n_{2}$ is the refractive index of the cylinder. The optical axis is the line between the cylinder front side and the cylinder back side at $y = 0 μ m$ .

In order to decompose the scattered field into different interactions, the Debye series based on Refs. 23 24.25.–26 is employed. Again, only perpendicular incidence is considered. The general Debye series for all incident angles can be found in Ref. 26. The Debye series introduces a geometrical series with transmission coefficients $T_{n}^{21}$ and $T_{n}^{12}$ and reflection coefficients $R_{n}^{11}$ and $R_{n}^{22}$ in order to identify single interactions. The superscript 2 stands for the surrounding medium while the superscript 1 stands for the region filled with the scatterer. $T_{n}^{21}$ , for example, is therefore the coefficient for transmission of light from region 2 to region 1 at the cylinder front side. The first reflection at the cylinder front side can be described with $R_{n}^{22}$ . If the light is transmitted into the cylinder through the front side ( $T_{n}^{21}$ ), any number of $p - 1$ reflections inside the scatterer ( $\frac{1}{2} \sum_{p = 1}^{\infty} T_{n}^{21} {(R_{n}^{11})}^{p - 1} T_{n}^{12}$ ) may occur before the light is transmitted ( $T_{n}^{12}$ ) again into the surrounding medium. Therefore, the expansion coefficients can be described as a geometrical series with $p - 1$ reflections inside the scatterer as²³^,²⁴

Eq. (11)

b_{n ∥} = \frac{1}{2} (1 - R_{n ∥}^{22} - \sum_{p = 1}^{\infty} (T_{n ∥}^{21} {(R_{n ∥}^{11})}^{p - 1} T_{n ∥}^{12})),

Eq. (12)

a_{n ⊥} = \frac{1}{2} (1 - R_{n ⊥}^{22} - \sum_{p = 1}^{\infty} (T_{n ⊥}^{21} {(R_{n ⊥}^{11})}^{p - 1} T_{n ⊥}^{12})) .

The factor

\frac{1}{2}

stands for diffraction.²⁶^–²⁸ Thus, a geometrical picture is used to understand the structure of the scattered fields, but the equations contain the exact Maxwell solution.

For simulations with a Gaussian beam in two dimensions in the paraxial approximation, the expansion coefficients given in Ref. 22 are modified, as described in Ref. 29. For simplicity, only the case of TM polarization and perpendicular incidence is considered here. The expansion coefficient can then be modified as

Eq. (13)

b_{n ∥}^{'} = C_{n 1} b_{n ∥},

where

Eq. (14)

C_{n 1} = (1 + \frac{n}{6} {(s \sqrt{2 i Q_{0}})}^{3} (Z^{3} - 3 Z)) C_{n 0},

and

Eq. (15)

C_{n 0} = e^{i k z_{0}} \sqrt{i Q_{0}} e^{- i Q_{0} \frac{{(k y_{0} + n)}^{2}}{w_{0}^{2} k^{2}}} .

As shown in Fig. 3,

z_{0}

is the focus position where the beam waist has the minimum value

w_{0}

.

y_{0}

is the lateral position of the beam, used for scanning over the sample. Furthermore, the abbreviations

Q_{0} = {[i - (2 z_{0} / k w_{0}^{2})]}^{- 1}

,

Z = s \sqrt{2 i Q_{0}} (n + k y_{0})

, and

s = (1 / k w_{0})

have been used (compare with Ref. 29). The validity of the paraxial approximation has been successfully tested for the used values for the beam waist

w_{0}

by comparing the beam to the beam constructed by the angular spectrum of plane waves method and by evaluating the integral given in Ref. 29 numerically.

Fig. 3

Gaussian beam in the paraxial approximation. The analytical solution allows scanning in $y$ -direction and shift of the focus in $z$ -direction (where the beam waist radius has the minimum value $w_{0}$ when the amplitude has dropped by $1 / e$ ). The focus is at $z = - z_{0}$ and $y = y_{0}$ .

The last modification of the expansion coefficients used in this work is for the case of layered cylinders. Light scattering by a layered cylinder can be calculated with a matrix formulation, as described in Refs. 30 31.–32. It is mentioned in passing that a part of the scattering algorithm for layered cylinders has been adapted from Ref. 31. For a layered cylinder with $N_{L}$ layers, radii $r_{j}$ ( $j = 1, \dots, N_{L}$ counted from the inner layers to the outer layers), refractive indices $n_{j}$ and perpendicular incidence of the illumination beam, the expansion coefficients $a_{n ⊥}$ and $b_{n ∥}$ are calculated as

Eq. (16)

a_{n ⊥} = \frac{\det A}{\det B},

Eq. (17)

b_{n ∥} = \frac{\det C}{\det D},

where

A_{i j}

,

B_{i j}

,

C_{i j}

, and

D_{i j}

are

2 N_{L} \times 2 N_{L}

matrices.³⁰

2. Calculation of an Optical Coherence Tomography Signal

The described equations are used to calculate the scattered field $E_{s} (ω)$ (compare Fig. 1). In Fourier-domain OCT, the detector measures the interferogram of the scattered field and the reference field and the measured data are transformed with the inverse Fourier transform to obtain a depth profile. For the scattering angle $θ = 180 \deg$ , the OCT intensity is given by

Eq. (18)

I_{OCT} (t, θ = 180 \deg) = \frac{1}{2} c ε_{0} n_{1} F^{- 1} ({| E_{s} (ω, 180 \deg) + E_{r} (ω) |}^{2}) .

Since the detector collects the spectral data coherently over a certain angular range, the integral describing a 2-D detection aperture can be written as

Eq. (19)

I_{OCT} (t) = \frac{1}{2} c ε_{0} n_{1} F^{- 1} ({| \int_{θ_{1}}^{θ_{2}} E_{s} (ω, θ) + E_{r} (ω) d θ |}^{2}) .

This is valid in the paraxial approximation under the assumption that the light coming from the cylinder is imaged onto the detector through a perfect imaging system without effects from finite lenses or aberrations. The integral can be evaluated analytically as

Eq. (20)

I_{OCT} (t) = \frac{1}{2} c ε_{0} n_{1} F^{- 1} (| e^{i \frac{3}{4} π} \sqrt{\frac{2}{π k (r_{2} + a)}} e^{i k (r_{2} + a)} (b_{0 ∥} (θ_{2} - θ_{1}) + \underset{\begin{matrix} n \neq 0 \end{matrix}}{\sum_{\begin{matrix} n = - \infty \end{matrix}}^{\infty}} b_{n ∥} (\frac{1}{- i n} e^{- i n θ_{2}} + \frac{1}{i n} e^{- i n θ_{1}})) + (θ_{2} - θ_{1}) E_{r} (ω) |^{2}) .

The single terms of the OCT signal can also be decomposed as¹

Eq. (21)

F^{- 1} ({| E_{s} (ω) + E_{r} (ω) |}^{2}) = F^{- 1} (E_{s} (ω) E_{s}^{*} (ω) + E_{r} (ω) E_{r}^{*} (ω) + E_{r} (ω) E_{s}^{*} (ω) + E_{s} (ω) E_{r}^{*} (ω)) .

These terms represent the autocorrelation term, the offset term at the origin $z = 0$ , and the cross-correlation signals, respectively. Instead of increasing the distance between the autocorrelation and the cross-correlation in the image by increasing $r_{3}$ relative to $r_{1} + r_{2}$ as shown in Fig. 1, (which requires higher $N$ because of the numerical Fourier transform), the term

Eq. (22)

I_{cc} (t) = \frac{1}{2} c ε_{0} n_{1} F^{- 1} (E_{s} {(ω)}^{*} E_{r} (ω))

can be directly evaluated for the cross-correlation intensity. For clarification, a detailed calculation is given for a simplified model for a Fourier-domain OCT A-scan: assuming for simplicity, a reference pulse of the form

E_{1} e^{- \frac{{(t - τ_{1})}^{2}}{b_{1}^{2}}}

with field amplitude

E_{1}

, temporal width

b_{1}

and delay

τ_{1}

because of the length of the reference arm, the field in the Fourier domain can be expressed as

E_{1} \frac{b_{1}}{\sqrt{2}} e^{- \frac{b_{1}^{2}}{4} ω^{2}} e^{i ω τ_{1}}

. Now, it is assumed for simplicity that a scatterer causes exactly two pulses in backward direction and it is assumed that the other signals can be neglected. The OCT signal is then the inverse Fourier transform of the sum of the reference pulse and the two signal pulses:

Eq. (23)

I_{OCT} (t, θ = 180 \deg) = \frac{1}{2} c ε_{0} n_{1} F^{- 1} ({| E_{1} \frac{b_{1}}{\sqrt{2}} e^{- \frac{b_{1}^{2}}{4} ω^{2}} e^{i ω τ_{1}} + E_{2} \frac{b_{2}}{\sqrt{2}} e^{- \frac{b_{2}^{2}}{4} ω^{2}} e^{i ω τ_{2}} + E_{3} \frac{b_{3}}{\sqrt{2}} e^{- \frac{b_{3}^{2}}{4} ω^{2}} e^{i ω τ_{3}} |}^{2}) .

Equation (23) can be simplified by using

{| α + β + γ |}^{2} = (α + β + γ) {(α + β + γ)}^{*}

and by executing the inverse Fourier transform. The constant

(1 / 2) c ε_{0} n_{1}

is omitted in the following explanation. The calculation gives nine terms: Three are of the form

E_{i}^{2} \frac{b_{i}}{2} e^{- \frac{t^{2}}{2 b_{1}^{2}}}, i \in {1,2, 3}

. These terms form the offset at the origin because they do not depend on

τ

. Furthermore, there are six correlation terms of the form:

Eq. (24)

E_{i} E_{j} \frac{b_{i} b_{j}}{\sqrt{2} \sqrt{b_{i}^{2} + b_{j}^{2}}} e^{- \frac{{(τ_{j} - τ_{i} + t)}^{2}}{b_{i}^{2} + b_{j}^{2}}}, i \in {1,2, 3}, j \in {1,2, 3}, i \neq j .

With

τ_{1}

being the time delay of the reference arm, the terms containing

τ_{1} - τ_{2}

,

τ_{1} - τ_{3}

,

τ_{2} - τ_{1}

, and

τ_{3} - τ_{1}

are the cross-correlation terms of the OCT A-scan. The A-scan is symmetric and the reference arm is chosen to be either much longer (

τ_{1} > τ_{2}, τ_{3}

) or much shorter (

τ_{1} < τ_{2}, τ_{3}

) than the sample arm so that the cross-correlations and autocorrelations do not overlap. In the first case, the pulse with

τ_{1} - τ_{2}

is on the negative side of the

z

-axis; in the second case, it is on the positive side. It can be seen from Eq. (24) that the OCT signal is wider than the width of the illumination pulse in the sample arm when comparing the denominators in the exponential function of the OCT cross-correlation signal to the original pulse from Eq. (1). The height of the OCT cross-correlation signal depends not only on

E_{1} E_{j}

but also on

b_{1}

and

b_{j}

. Based on the descriptions in Refs. 1, 33, and 34 for the detected OCT signal, the maxima of the OCT cross-correlation signal

I_{cc}

are compared to the square root of the time-resolved backscattered intensity

\sqrt{I_{s}} = \sqrt{{| E_{s} |}^{2}}

(compare Fig. 1). We scale

\sqrt{I_{s}}

to the maxima of

I_{cc}

by a factor of

C_{sca}

, which is defined by the relation [compare Eq. (24)]:

Eq. (25)

E_{j} E_{1} b_{1} b_{j} \frac{1}{\sqrt{2}} \frac{1}{\sqrt{b_{1}^{2} + b_{j}^{2}}} = C_{sca} E_{j} .

In Eq. (25),

C_{sca}

depends on

b_{1}

,

b_{j}

, and

E_{1} = \sqrt{I_{1}}

.

b_{1}

and

E_{1}

are known constants from the incident pulse, whereas

b_{j}

is the width of the signal after scattering. If it is assumed that the width of the pulses after the scattering process is not changed compared to

b_{1}

, which the reference and the incident pulse

E_{r}

and

E_{i}

had at the beam splitter (

b_{j} = b_{1}

), the scaling factor is calculated to be

Eq. (26)

C_{sca} = E_{1} \frac{b_{1}}{2} .

Therefore, the maxima of the OCT cross-correlation signal intensity are predicted to equal the square root of the scattered intensity $\sqrt{I_{s}}$ as long as the width of the pulses is not changed in the scattering process.

It is mentioned in passing that the simulations for geometrical optics have been conducted with a ray tracing algorithm based on Snell’s law and Fresnel’s equations and that a random number is drawn for the starting position of the light package between $z = - a$ , $y = - a$ and $z = - a$ , $y = a$ to simulate a plane wave incident on the cylinder cross section. The detected counts are converted to intensities as

Eq. (27)

I_{m c} = \frac{1}{2} c ε_{0} {| E_{0} |}^{2} \frac{n_{counts}}{\sum n_{counts}} \frac{A_{ill}}{A_{bin}},

where

n_{counts}

is the number of counts,

E_{0}

is the amplitude from Eq. (1),

A_{ill} = 2 a

is the geometrical illumination line from

- a

to

a

, and

A_{bin} = [2 π (r_{2} + a)] / n_{θ}

is the angular bin size when

n_{θ}

angles are calculated. Further literature concerning ray tracing and its modifications can be found in Govaerts et al.³⁵ and in Lugovtsov et al.³⁶

3. Results

3.1.

Time-Resolved Intensity for Plane Wave Scattering

Before the OCT results are shown, the scattered intensities for plane wave illumination are discussed from the point of view of geometrical optics. This allows the attribution of pathways to every signal. Unless otherwise noted, only scattering has been considered in the simulations so that the imaginary part of the refractive index is zero throughout this work. Figure 4 shows the graphical output of the implemented ray tracing algorithm for a cylinder cross section based on Fresnel’s equations and Snell’s law for three selected pathways for plane wave illumination of a cylinder with radius $a = 20 μ m$ and $n_{2} = 1.4$ in air $n_{1} = 1$ . The pathways are chosen so that the rays exit at $θ \approx 180 \deg$ . Figure 4 shows that the light packages either travel along the optical axis ( $y \approx 0 μ m$ ) so that the path lengths $d$ counted from $z = - a$ are multiples of the radius $a$ times the refractive index $n_{2}$ so that $d = 4 m a n_{2}$ with $m \in N_{0}$ or that the light packages follow pathways that form triangles or starlike shapes. For further reference, the rays in the geometrical optics picture are numbered according to how often the light package interacts with the boundary of the cross section. If the interaction is linked to a starlike pathway that can exit with $θ = 180 \deg$ for the given refractive indices $n_{1}$ and $n_{2}$ , the interaction number is denoted with a prime (′). Thus, the labels allow to distinguish between multiple reflections that occur on the optical axis and signals from more complex pathways. This number system is used throughout the rest of this work. The interactions in Fig. 4 are 1, 6′ and 7′, accordingly. In Fig. 5, the calculation is extended to all scattering angles $θ$ and the Maxwell solution is compared to the geometrical optics picture. $E_{0}$ is set to $1 \frac{V}{m}$ throughout this work [compare Eq. (1)]. Figure 5 shows the scattered 2-D intensity that an ultrafast detector would register for each scattering angle $θ$ computed with (a) the implemented Maxwell solution from Eq. (6) and (b) the ray tracing algorithm for geometrical optics from Eq. (27) for a cylinder with radius $a = 20 μ m$ and a refractive index of $n_{2} = 1.4$ in air ( $n_{1} = 1$ ). The light is polarized parallel to the cylinder axis (TM case) and it is assumed that the refractive index does not depend on the wavelength. The Maxwell algorithm calculates Eq. (6) for $λ = 1305 \pm 400 nm$ . The first seven interactions are shown on the right-hand side of Fig. 5. The signals 1 and 3 labeled on the right side of Fig. 5 indicate the cylinder cross section. The dashed black line in Fig. 5 shows the geometrically calculated travel distance of the light from the source to the cylinder front side and to the detector [compare Fig. 5(c)]. The intensity $I_{m c}$ has been calculated, as described in Eq. (27). In both simulations, the front and partially back side of the cylinder are visible as well as multiple signals behind the scatterer. Both the Maxwell solution and the ray tracing solution show roughly the same pattern including the cylinder front side 1 (labeled on the right side of Fig. 5), the back side 3, the multiple reflections from the optical axis and the more complex geometrical pathways at higher depths in $z$ ( $y \neq 0 μ m$ ). The measured distance between signal 1 and 3 in the Maxwell simulation in Fig. 5 at $θ = 180 \deg$ is $d_{sim} = 55.9 μ m$ while the calculated distance is $d = 2 n_{2} a = 56 μ m$ . The pathways in the geometrical optics simulation in Fig. 5(b) are labeled as described in Fig. 4. Signal 2 in Fig. 5(b) reaches only the right half space behind the cylinder back side ( $z > a$ , $θ \in [0, \dots, 90 \deg)$ , $θ \in (270, \dots, 360 \deg]$ ), as is expected from a ray that is transmitted through the scatterer. Signal 3 in Fig. 5(b) appears right behind the cylinder side and consists of two parts: The reflection along the optical axis, marking the backside of the cylinder ( $z = 2 n_{2} a$ ), and the starlike pathway that is longer. It is not possible for the latter pathway to exit at $θ = 180 \deg$ because of the chosen refractive index. In the Maxwell solution in Fig. 5(a), however, this signal reaches all the way to $θ = 180 \deg$ and beyond. The reason why light is able to exit at this particular angle is the effect of surface waves. These modes that travel along the circumference of the cylinder are the WGMs, which are well-described in literature³⁷^–⁴⁰ with useful applications.⁴¹^,⁴² The patterns in the Maxwell solution are in general broader and reach positions that are not reached in geometrical optics, which shows the wave behavior of light. The Maxwell solution gives an image that is fairly similar to the geometrical optics picture when the periodically appearing structures that are the WGMs radiating into the far field are ignored. The ray tracing solution and the Maxwell solution agree better when the cylinder diameter is very large compared to the wavelengths of the incident pulse.

Fig. 4

Graphical illustration of a 2-D ray tracing simulation for three selected light packages for plane wave illumination of a cylinder with a radius $a = 20 μ m$ and $n_{2} = 1.4$ in air. The incident geometrical rays are marked in red, the internal and the outgoing ones have different colors. The number of interactions with the boundary is 1 for the middle ray and 6′ and 7′ for the others, respectively. The light is polarized parallel to the cylinder axis.

Fig. 5

The Maxwell solution is compared to the geometrical optics picture: (a) Maxwell and (b) ray tracing simulation for a cylinder with $20 μ m$ radius. (c) The simulated setup. The noise from the inverse Fourier transform has been removed so that only intensities greater than or equal to $e^{- 22}$ have been plotted in the Maxwell case. The 2-D intensities coded in the color values are on a logarithmic scale (to base $e$ ) in both plots. The distances are calculated from $z = t c$ .

3.2.

Optical Coherence Tomography A-Scan for Plane Wave Scattering

With the basic effects of cylinder scattering explained for a plane wave at perpendicular incidence, the computational results for a plane wave OCT are discussed. With Eqs. (6) and (18), the A-scan shown in Fig. 6 is computed. The offset, the autocorrelation, and the cross-correlation terms are visible. The graph is symmetric as expected from Eq. (21). The $z$ -axis is shifted so that the origin is at the cylinder front side just like in Fig. 5. The autocorrelation signals overlap slightly with the cross-correlation signals. This is due to the fact that $r_{3}$ is not different enough from $r_{1}$ and $r_{2}$ , but increasing $r_{1}$ , $r_{2}$ , or $r_{3}$ increases the calculation time because of the required $N$ (compare with Sec. 1.1). The cross-correlations are investigated in Fig. 7 with Eq. (22).

Fig. 6

A-scan for an OCT system with $λ = 845 \pm 300 nm$ with $8 \times 10^{4}$ calculated spectral data points in backscattering direction ( $θ = 180 \deg$ ). The offset in the middle, the autocorrelations and the cross-correlations are visible. The distances used are $r_{1} = 6 a$ , $r_{2} = 4 a$ , and $r_{3} = 180 a$ . The cylinder has a refractive index of $n_{2} = 1.4$ in air and its radius is $a = 20 μ m$ .

Fig. 7

$\sqrt{I_{m c}}$ , $\sqrt{I_{s}}$ and the OCT cross-correlation signal $I_{c c}$ [compare Eqs. (6) and (22)] versus $z$ . The ray tracing simulation has 361 angular bins, 5001 temporal bins and $10^{9}$ light packages. The OCT parameters are the same as in Fig. 6 except for $N = 8 \cdot 10^{5}$ , $r_{1} = a$ , $r_{2} = a$ , and $r_{3} = 3 a$ . For the OCT signal, only the cross-correlation part according to Eq. (22) has been used.

Figure 7 shows the OCT cross-correlation signal from Fig. 6 calculated with Eq. (22), the square root of the backscattered intensity, and the results of the ray tracing algorithm for $θ = 180 \deg$ . The counts of the ray tracing algorithm are converted with Eq. (27). The ray tracing plot is drawn with a broader line for clarity. The square root of the backscattered field $\sqrt{I_{s}}$ (that an ultrafast detector for time-resolved measurements would detect) is scaled to the maxima of the OCT peaks with $\frac{1}{2} c ε_{0} n_{1} C_{sca}$ with Eq. (26) just like $I_{m c}$ . The different insets show the signals 1 and 3 in detail as well as the signals 5, 6′, 7′, and 11′ belonging to longer pathways and the WGMs. The first six WGMs are marked with “WGM.” It can be seen that the ray tracing algorithm for geometrical optics predicts the cylinder front side (interaction 1), the back side (interaction 3) and signals at $z = 112 μ m$ (interaction 5, compare Fig. 4), $z = 130.4 μ m$ (the ray that touches the cylinder boundary 6′ times, compare Fig. 4), $z = 138.6 μ m$ (7′ interactions), and $z = 217.8 μ m$ (11′). The intensity of the OCT signals that can be associated with a pathway from geometrical optics decreases faster than the intensity of the gallery modes. Furthermore, the intensity of the OCT cross-correlation signal $I_{cc}$ is proportional to the square root of the backscattered intensity $\sqrt{I_{s}}$ , as calculated in Eqs. (24) and (26) as long as the width of the scattered pulses does not change. As expected from Eq. (24), the width of the OCT peaks is in general broader than the width of $\sqrt{I_{s}}$ .

3.3.

Near Fields and Whispering Gallery Modes

Figure 8 shows snapshots of the cylinder near fields. The incident plane wave in (a) hits the cylinder boundary from the left and part of it is reflected and a part of it is transmitted into the cylinder. The wave outside where $n_{1} = 1$ travels faster than the part of the wave that has to travel through the cylinder with $n_{2} = 1.4$ . At the cylinder back side, a part of the wave is transmitted, a part of it is reflected in 180-deg direction again and the WGMs appear. The part that is reflected backward is either reflected again or transmitted, but in both cases, the intensity has decreased a lot compared to the intensity of the gallery modes, which form at the back side of the cylinder in Fig. 8(b) and which start to travel along the circumference (c). Because of dispersion, their width changes the longer they travel (d). While the light travels a distance $d = n 4 a n_{2}, n \in N_{0}$ along the optical axis when detected in backscattering direction $θ = 180 \deg$ at $y = 0 μ m$ (compare Fig. 7), the first gallery modes WGM 1 travel a distance $d = 2 a n_{2} + π a {\tilde{n}}_{2}$ . $2 a$ appears because they first form at the boundary after the light has passed through the cylinder. ${\tilde{n}}_{2}$ indicates that the refractive index for the surface waves corresponds neither exactly to the refractive index $n_{2}$ of the cylinder nor $n_{1}$ of the surrounding medium. This is due to the fact that the WGMs do not exactly travel a circular pathway around the cylinder.⁴⁴

Fig. 8

Temporal snapshots of the near fields around a cylinder with radius $a = 3.5 μ m$ and a refractive index of $n_{2} = 1.4$ in air at $λ = 845 nm$ . The WGM 1 (compare Fig. 7) appear first at the back side of the cylinder after the light has passed through its diameter. The part of the wave that is inside the cylinder with $n_{2} = 1.4$ is distorted relative to the part of the wave in the surrounding medium with $n_{1} = 1$ depending on the refractive index ratio. The intensity of the gallery modes decreases only slowly while the other signals decay faster. The width of the gallery modes changes the longer they travel along the circumference. The colormap used to display the results in Figs. 8, 10, 12, and 13 is taken from Ref. 43.

3.4.

Plane Wave Optical Coherence Tomography with the Debye Series

Another method to understand the origin of the cross-correlation signals is to decompose the scattered field into single interaction terms. This is done in the Debye series for plane wave scattering (compare Refs. 23 24.25.–26). To the knowledge of the authors, the Debye series has not been implemented for an OCT algorithm before. Figure 9 shows the OCT cross-correlation signal from Eq. (22) for different $p$ [compare Eqs. (11) and (12)] for $θ = 180 \deg$ . A chosen $p$ includes all $p + 1$ interactions. It can be seen in Figs. 9(a)–9(f) that again the contributions from the cylinder front side, the back side, and the WGMs are recovered. When $p = 2$ is calculated in (b), not only the signal at $z = 2 a n_{2}$ appears from the optical axis but also the peak belonging to the WGM 1. This is in agreement with the observation in the near field FDTD simulation (Fig. 8) that the gallery modes WGM 1 are formed at the cylinder back side and that they travel along the cylinder circumference. Lock et al. show for calculations with the Debye series that scattered light reaches regions that are not accessible in geometrical optics, but their calculations are not within the context of OCT imaging.²⁵^,²⁷^,⁴⁵

Fig. 9

Simulation of an OCT cross-correlation signal with the Debye series. The geometrical series that constructs the expansion coefficients reveals the single interaction terms. $p - 1$ denotes the number of internal reflections, the peaks are numbered as interactions. The incident light is a plane wave with $λ = 845 \pm 300 nm$ and the cylinder has a radius of $3.5 μ m$ and a refractive index of $n_{2} = 1.4$ in air $n_{1} = 1.0$ . Only $θ = 180 \deg$ is considered and $N = 8000$ frequencies have been calculated. The interactions show (a) the signal from the cylinder front side, (b) the signal from the front side, the back side and the first whispering gallery mode and (c)–(f) more signals from internal reflections and the gallery modes. The intensity is calculated with Eq. (22).

3.5.

Gaussian Beam Optical Coherence Tomography

In reality, most OCT devices employ a focused beam that scans over the sample to improve the lateral resolution of the images. To include the lateral scanning, a Gaussian beam has been implemented as shown in Fig. 3. The 2-D Gaussian beam has been implemented according to Eqs. (13) and (14) with an axial focus point at $z_{0} = 0 μ m$ at the cylinder center. The used positions for the Gaussian beam in Fig. 10 in lateral direction range from $y_{0} = - 26 μ m$ to $y_{0} = 26 μ m$ with a distance of $Δ y = 2 μ m$ between the channels. Only the cross-correlation part of the 2-D intensity in the units ${Wm}^{- 1}$ has been considered according to Eq. (22). The zero-position $y_{0} = 0 μ m$ is scanned as well. Thus, in total, 27 $y_{0}$ -values have been used. The cylinder is outlined in black. The distances in the OCT system are $r_{1} = a$ , $r_{2} = a$ , and $r_{3} = 3 a$ and the detected signals are labeled up to interaction 11 and WGM 5. The figure shows the resulting image with the Gaussian beam for the scattering angle $θ = 180 \deg$ . The cylinder has the same specifications as in Fig. 5 and the beam waist of the Gaussian beam has been chosen to be $w_{0} = 5 μ m$ . Since the beam waist is smaller than the cylinder diameter, the surface waves, labeled WGM 1, WGM 2, WGM 3, and so on, are only excited when $y_{0}$ is in the proximity of the cylinder radius. Illumination in the vicinity of the optical axis $y_{0} \approx 0 μ m$ does not excite surface waves, but the other signals that have an associated pathway in the geometrical optics picture can be seen. Since only the scattering angle $θ = 180 \deg$ is considered, no curvature of the cylinder surfaces is visible. Signals from other angles are not detected. The distance between the cylinder front and back side (interactions 1 and 3) at $y_{0} = 0 μ m$ is the expected value $d = 2 a n_{2}$ . At the upper and lower lateral sides of Fig. 3 (when $y_{0} \neq 0 μ m$ ), patterns that are not labeled can be seen in the background. It has been shown in previous work⁴⁶ that they stem from the limited length of the numerical inverse Fourier transform and that they decrease for higher $N$ . At large axial depths $z$ , the gallery modes show again the expected change in width.

Fig. 10

B-scan for a cylinder with a refractive index of $n_{2} = 1.4$ and radius $a = 20 μ m$ in air ( $n_{1} = 1$ ). The wavelengths used are $λ = 845 \pm 300 nm$ with $N {= 24 \times 10}^{4}$ data points. The plot shows the logarithm of the 2-D intensity in the units ${Wm}^{- 1}$ to base $e$ .

Figure 11 shows the pulse that is reflected from the first surface (interaction 1) during the lateral scan over $y_{0}$ . The values are the same as those used for Fig. 10 except for $z_{0} = a$ . Thus, we can compare the lateral width of the reflected signal to the minimum width $w_{0}$ of the incident Gaussian beam. The width that is found in Fig. 11 has the width $w_{0}$ of the incident Gaussian pulse in $y_{0}$ direction. This is due to the fact that only the part of the Gaussian that is backscattered in $θ = 180 ‐ \deg$ direction is detected. Since the focus with the minimum beam waist $w_{0}$ is located at the cylinder surface $z_{0} = a$ , the beam waist $w_{0}$ is the width that is found by the fitting algorithm with the fitted function $f (x) = a_{1} e^{{(\frac{x - b_{1}}{c_{1}})}^{2}}$ .

Fig. 11

The fitting algorithm is applied to the data at $z = 0 μ m$ along the $y_{0}$ -values. The fitting parameter $c_{1} = 5.000 \pm 0.017 μ m$ is found for the width.

3.5.1.

Optical coherence tomography with an aperture

In Fig. 12, the OCT cross-correlation signal simulated with an aperture is shown for a Gaussian beam with a beam waist of $w_{0} = 3.9 μ m$ . The angular range included in Eq. (19) is $180 \pm 80 \deg$ . The integral from Eq. (19) is solved numerically for a discrete number of angles in the first two images (from left to right) and the third image shows the analytical solution from Eq. (20). It can be seen that the numerical solution with 1001 angles on the left and the analytical solution on the right agree well with each other. The numerical solution with only 11 angles has not converged to the correct result yet. For 1001 angles and for the analytical solution, the signals at the lateral cylinder sides do not lie exactly on the surface and additionally, fewer signals are observed compared to 11 angles. The interference in the coherent detection mode prevents imaging of the curvature of the smooth cylinder surface. The curvature of the cylinder surface is visible for the range of $θ = 180 \pm 80 \deg$ in Fig. 5. This is due to the fact that in Fig. 5, no interferometer setup and no aperture integral for coherent detection is calculated. Similar simulated and experimental OCT images are presented in Yi et al.²¹ for calculations with spheres, but the WGMs and the individual pathways according to geometrical optics were not discussed in depth. For rough surfaces, the curvatures of the scatterers are visible both in simulated and experimental OCT images.¹^,⁸ We expect to see the WGMs and the missing curvature in the tomograms when the scatterer is similar to a cylinder or a sphere with a smooth surface (when absorption is negligible). For example, we predict that measuring a glass fiber with a smooth surface will produce OCT images showing the gallery modes.

Fig. 12

OCT cross-correlation signal for calculation with an aperture according to Eq. (19). The two images (a) and (b) are numerical solutions and image (c) is the analytical solution with Eq. (20). The used values are $N = 8000$ , $λ = 1300 \pm 300 nm$ , $y_{0}$ from $- 76 μ m$ to $76 μ m$ with $Δ y_{0} = 2 μ m$ , $z_{0} = 0 μ m$ , $r_{1} = a$ , $r_{2} = a$ and $r_{3} = 3 a$ . The cylinder radius is $a = 70 μ m$ . The logarithm of the 2-D intensity is plotted to base $e$ .

3.5.2.

Optical coherence tomography simulation for a layered cylinder

Figure 13 shows the simulated cross-correlation OCT image for a layered cylinder with three layers for TM polarization according to Eq. (17). Equation (20) has been used to simulate the collection of the signal over a certain angular range in the paraxial approximation. Only $I_{cc}$ has been calculated according to Eq. (22). The angular range has been chosen to be $θ = 180 \pm 10 \deg$ while the beam waist of the Gaussian beam is $w_{0} = 5 μ m$ . The layered cylinder consists of a core with a diameter of $4 μ m$ and a refractive index 1.2, a layer with $4 μ m$ and 1.5 and an outer layer with a thickness of $2 μ m$ and a refractive index of 1.4. The cylinder layers scaled by the refractive indices are drawn in black. The image shows that OCT can identify the single layers, but at the same time additional signals appear between and behind the real layers of the cylinder. The intensity of the signals depends a lot on the choice of the refractive indices. Another important observation is that with the chosen angular detection range, no curvature of the surface is visible. A comparison with Fig. 5 shows that the angular range is too small for the curvature of the scattered field to be visible.

Fig. 13

Cross-correlation signal of an OCT B-scan for a layered cylinder. $z$ has again been shifted and divided by a factor of 2. The incident Gaussian beam has a beam waist of $w_{0} = 5 μ m$ and $λ = 845 \pm 300 nm$ has been used with $N = 8000$ values. The used positions for the Gaussian beam in lateral direction range from $- 9 μ m$ to $9 μ m$ with a distance of $1 μ m$ between the channels. In total, 19 $y_{0}$ -values have been used. The focus is at the cylinder center. The logarithm to base $e$ of the 2-D intensity is shown. The distances in the OCT system are $r_{1} = a$ , $r_{2} = a$ , and $r_{3} = 3 a$ . The cylinder layers are drawn in black.

Similar to the way the simulated tomograms in this work do not show the whole cross section of the circular structures investigated, OCT images of phantoms mimicking lung tissue only show signals from the front and back side of the scatterer in literature.⁴⁷^,⁴⁸ Furthermore, in clinical OCT images of swine lung tissue, artifacts are described that show extra layers⁴⁹ similar to the extra layers due to multiple reflections in this work. We expect WGMs to contribute to clinical OCT tomograms, especially since the resonance condition for surface waves is not only fulfilled for the case of circular cross sections.⁵⁰^–⁵² However, many texts in OCT literature do not distinguish gallery modes from other types of artifacts or they do not mention these artifacts at all. Therefore, the investigation of WGMs continues to be an interesting topic for future work.

4. Conclusions

A 2-D algorithm for the simulation of image formation in OCT for scattering of a plane wave and of a Gaussian beam by an infinitely long homogeneous dielectric cylinder has been implemented. The basic effects of cylinder scattering have been investigated and the pathways from the geometrical optics picture have been compared to the full-wave Maxwell solution. The relation between the OCT signals and the backscattered light has been calculated and the scaling factor is used to compare the square root of the backscattered light and the OCT signal directly. The single contributions from the scattering interactions in geometrical optics have been labeled and linked to photon pathways. With the Debye series and with the ray pathways calculated from geometrical optics, the single signals have been explained and the differences between geometrical optics and the wave nature of light have been linked to the observed OCT signals. Both for a plane wave and for a scanning Gaussian beam, the WGMs, which travel along the circumference of the cylinder, have been identified in the Maxwell solution. The simulated OCT image of a layered cylinder has shown the different layers scaled with their respective refractive indices. The integral over a certain angular range for the 2-D aperture has shown that the curvature of the smooth cylinder in lateral direction is not visible in the computed simulations. It is shown that OCT images do not always display the real surface of the investigated sample.

Acknowledgments

We acknowledge the support of Deutsche Forschungs-gemeinschaft DFG. Thomas Brenner acknowledges the support by the Landesgraduiertenförderungsgesetz Baden-Württemberg.

References

1.

W. Drexler and J. G. Fujimoto, Optical Coherence Tomography: Technology and Applications, Springer, Berlin (2008). Google Scholar

2.

A. Knüttel et al., “Spatially confined and temporally resolved refractive index and scattering evaluation in human skin performed with optical coherence tomography,” J. Biomed. Opt., 5 (1), 83 –92 (2000). http://dx.doi.org/10.1117/1.429972 JBOPFO 1083-3668 Google Scholar

3.

A. Knüttel et al., “New method for evaluation of in vivo scattering and refractive index properties obtained with optical coherence tomography,” J. Biomed. Opt., 9 (2), 265 –273 (2004). http://dx.doi.org/10.1117/1.1647544 JBOPFO 1083-3668 Google Scholar

4.

A. Dunn and R. Richards-Kortum, “Three-dimensional computation of light scattering from cells,” IEEE J. Sel. Top. Quantum Electron., 2 (4), 898 –905 (1996). http://dx.doi.org/10.1109/2944.577313 Google Scholar

5.

R. Drezek, A. Dunn and R. Richards-Kortum, “Light scattering from cells: finite-difference time-domain simulations and goniometric measurements,” Appl. Opt., 38 (16), 3651 –3661 (1999). http://dx.doi.org/10.1364/AO.38.003651 APOPAI 0003-6935 Google Scholar

6.

A. A. Lindenmaier et al., “Texture analysis of optical coherence tomography speckle for characterizing biological tissues in vivo,” Opt. Lett., 38 (8), 1280 –1282 (2013). http://dx.doi.org/10.1364/OL.38.001280 OPLEDP 0146-9592 Google Scholar

7.

M. Kirillin et al., “Speckle statistics in OCT images: Monte Carlo simulations and experimental studies,” Opt. Lett., 39 (12), 3472 –3475 (2014). http://dx.doi.org/10.1364/OL.39.003472 OPLEDP 0146-9592 Google Scholar

8.

M. Kirillin et al., “Simulation of optical coherence tomography images by Monte Carlo modeling based on polarization vector approach,” Opt. Express, 18 (21), 21714 –21724 (2010). http://dx.doi.org/10.1364/OE.18.021714 OPEXFF 1094-4087 Google Scholar

9.

I. Meglinski et al., “Simulation of polarization-sensitive optical coherence tomography images by a Monte Carlo method,” Opt. Lett., 33 (14), 1581 –1583 (2008). http://dx.doi.org/10.1364/OL.33.001581 OPLEDP 0146-9592 Google Scholar

10.

G. Yao and L. V. Wang, “Monte Carlo simulation of an optical coherence tomography signal in homogeneous turbid media,” Phys. Med. Biol., 44 (9), 2307 (1999). http://dx.doi.org/10.1088/0031-9155/44/9/316 PHMBA7 0031-9155 Google Scholar

11.

A. Kienle and J. Schäfer, “Simulated and measured optical coherence tomography images of human enamel,” Opt. Lett., 37 (15), 3246 –3248 (2012). http://dx.doi.org/10.1364/OL.37.003246 OPLEDP 0146-9592 Google Scholar

12.

A. E. Hartinger et al., “Monte Carlo modeling of angiographic optical coherence tomography,” Biomed. Opt. Express, 5 (12), 4338 –4349 (2014). http://dx.doi.org/10.1364/BOE.5.004338 BOEICL 2156-7085 Google Scholar

13.

R. Su et al., “Optical coherence tomography for quality assessment of embedded microchannels in alumina ceramic,” Opt. Express, 20 (4), 4603 –4618 (2012). http://dx.doi.org/10.1364/OE.20.004603 OPEXFF 1094-4087 Google Scholar

14.

T. B. Swedish et al., “Computational model of optical scattering by elastin in lung,” Proc. SPIE, 7904 79040H (2011). http://dx.doi.org/10.1117/12.875707 PSISDG 0277-786X Google Scholar

15.

D. C. Reed and C. A. DiMarzio, “Computational model of OCT in lung tissue,” Proc. SPIE, 7570 75700I (2010). http://dx.doi.org/10.1117/12.842461 PSISDG 0277-786X Google Scholar

16.

A. S. F. Silva and A. L. Correia, “From optical coherence tomography to Maxwell’s equations,” in IEEE 3rd Portuguese Meeting in Bioengineering, 1 –4 (2013). http://dx.doi.org/10.1109/ENBENG.2013.6518419 Google Scholar

17.

Y. T. Hung, S. L. Huang and S. H. Tseng, “Full EM wave simulation on optical coherence tomography: impact of surface roughness,” Proc. SPIE, 8592 859216 (2013). http://dx.doi.org/10.1117/12.2006107 PSISDG 0277-786X Google Scholar

18.

S. H. Huang, S. J. Wang and S. H. Tseng, “Tomographic reconstruction of melanin structures of optical coherence tomography via the finite-difference time-domain simulation,” Proc. SPIE, 9328 93281T (2015). http://dx.doi.org/10.1117/12.2079309 PSISDG 0277-786X Google Scholar

19.

P. R. Munro, A. Curatolo and D. D. Sampson, “Full wave model of image formation in optical coherence tomography applicable to general samples,” Opt. Express, 23 (3), 2541 –2556 (2015). http://dx.doi.org/10.1364/OE.23.002541 OPEXFF 1094-4087 Google Scholar

20.

P. R. Munro, “Exploiting data redundancy in computational optical imaging,” Opt. Express, 23 (24), 30603 –30617 (2015). http://dx.doi.org/10.1364/OE.23.030603 OPEXFF 1094-4087 Google Scholar

21.

J. Yi, J. Gong and X. Li, “Analyzing absorption and scattering spectra of micro-scale structures with spectroscopic optical coherence tomography,” Opt. Express, 17 (15), 13157 –13167 (2009). http://dx.doi.org/10.1364/OE.17.013157 OPEXFF 1094-4087 Google Scholar

22.

C. F. Bohren and D. R. Huffman, Absorption and Scattering of Light by Small Particles, John Wiley & Sons, New York (2008). Google Scholar

23.

R. Li et al., “Debye series of normally incident plane-wave scattering by an infinite multilayered cylinder,” Appl. Opt., 45 (24), 6255 –6262 (2006). http://dx.doi.org/10.1364/AO.45.006255 APOPAI 0003-6935 Google Scholar

24.

R. Li, X. Han and K. F. Ren, “Generalized Debye series expansion of electromagnetic plane wave scattering by an infinite multilayered cylinder at oblique incidence,” Phys. Rev. E, 79 (3), 036602 (2009). http://dx.doi.org/10.1103/PhysRevE.79.036602 PLEEE8 1539-3755 Google Scholar

25.

E. A. Hovenac and J. A. Lock, “Assessing the contributions of surface waves and complex rays to far-field Mie scattering by use of the Debye series,” J. Opt. Soc. Am. A, 9 (5), 781 –795 (1992). http://dx.doi.org/10.1364/JOSAA.9.000781 JOAOD6 0740-3232 Google Scholar

26.

J. A. Lock and C. L. Adler, “Debye-series analysis of the first-order rainbow produced in scattering of a diagonally incident plane wave by a circular cylinder,” J. Opt. Soc. Am. A, 14 (6), 1316 –1328 (1997). http://dx.doi.org/10.1364/JOSAA.14.001316 JOAOD6 0740-3232 Google Scholar

27.

J. A. Lock and P. Laven, “Mie scattering in the time domain. Part II. The role of diffraction,” J. Opt. Soc. Am. A, 28 (6), 1096 –1106 (2011). http://dx.doi.org/10.1364/JOSAA.28.001096 JOAOD6 0740-3232 Google Scholar

28.

H. C. Hulst, Light Scattering by Small Particles, Courier Corporation, New York (1957). Google Scholar

29.

Z. W. L. Guo, “Electromagnetic scattering from a multilayered cylinder arbitrarily located in a Gaussian beam, a new recursive algorithms,” Prog. Electromagn Res., 18 317 –333 (1998). http://dx.doi.org/10.2528/PIER97071100 PELREX 1043-626X Google Scholar

30.

J.-P. Schäfer, “Implementierung und Anwendung analytischer und numerischer Verfahren zur Lösung der Maxwellgleichungen für die Untersuchung der Lichtausbreitung in biologischem Gewebe,” Ulm (2011). Google Scholar

31.

J. P. Schäfer, “Matscat,” (2014) http://www.mathworks.com/matlabcentral/fileexchange/36831-matscat October ). 2014). Google Scholar

32.

M. Kerker, Scattering of Light and Other Electromagnetic Radiation, Academic Press, New York (1969). Google Scholar

33.

Y. Pan et al., “Low-coherence optical tomography in turbid tissue: theoretical analysis,” Appl. Opt., 34 (28), 6564 –6574 (1995). http://dx.doi.org/10.1364/AO.34.006564 APOPAI 0003-6935 Google Scholar

34.

V. M. Kodach et al., “Quantitative comparison of the OCT imaging depth at 1300 nm and 1600 nm,” Biomed. Opt. Express, 1 176 –185 (2010). http://dx.doi.org/10.1364/BOE.1.000176 BOEICL 2156-7085 Google Scholar

35.

Y. M. Govaerts and M. M. Verstraete, “Raytran: a Monte Carlo ray-tracing model to compute light scattering in three-dimensional heterogeneous media,” IEEE Trans. Geosci. Remote Sens., 36 (2), 493 –505 (1998). http://dx.doi.org/10.1109/36.662732 Google Scholar

36.

A. E. Lugovtsov, S. Y. Nikitin and A. V. Priezzhev, “Ray-wave approximation for calculating laser radiation scattering by a transparent dielectric spheroidal particle,” Quantum Electron., 38 (6), 606 –611 (2008). http://dx.doi.org/10.1070/QE2008v038n06ABEH013941 QUELEZ 1063-7818 Google Scholar

37.

A. N. Oraevsky, “Whispering-gallery waves,” Quantum Electron., 32 (5), 377 –400 (2002). http://dx.doi.org/10.1070/QE2002v032n05ABEH002205 QUELEZ 1063-7818 Google Scholar

38.

O. Sbanski et al., “Elastic light scattering from single microparticles on a femtosecond time scale,” J. Opt. Soc. Am. A, 17 (2), 313 –319 (2000). http://dx.doi.org/10.1364/JOSAA.17.000313 JOAOD6 0740-3232 Google Scholar

39.

M. Gomilšek, Whispering Gallery Modes, University of Ljubljana, Ljubljana Google Scholar

40.

H. Bech and A. Leder, “Particle sizing by time-resolved Mie calculations—a numerical study,” Optik—Int. J. Light Electron Opt., 117 (1), 40 –47 (2006). http://dx.doi.org/10.1016/j.ijleo.2005.06.008 Google Scholar

41.

V. Braginsky, M. Gorodetsky and V. Ilchenko, “Quality-factor and nonlinear properties of optical whispering-gallery modes,” Phys. Lett. A, 137 (7), 393 –397 (1989). http://dx.doi.org/10.1016/0375-9601(89)90912-2 PYLAAG 0375-9601 Google Scholar

42.

S. Arnold et al., “Shift of whispering-gallery modes in microspheres by protein adsorption,” Opt. Lett., 28 (4), 272 –274 (2003). http://dx.doi.org/10.1364/OL.28.000272 OPLEDP 0146-9592 Google Scholar

43.

K. D. Moreland, “Diverging color maps for scientific visualization,” 2009). Google Scholar

44.

A. Boleininger et al., “Whispering gallery modes in standard optical fibres for fibre profiling measurements and sensing of unlabelled chemical species,” Sensors, 10 (3), 1765 –1781 (2010). http://dx.doi.org/10.3390/s100301765 SNSRES 0746-9462 Google Scholar

45.

J. A. Lock and P. Laven, “Mie scattering in the time domain. Part I. The role of surface waves,” J. Opt. Soc. Am. A, 28 (6), 1086 –1095 (2011). http://dx.doi.org/10.1364/JOSAA.28.001086 JOAOD6 0740-3232 Google Scholar

46.

T. Brenner, D. Reitzle and A. Kienle, “An algorithm for simulating image formation in optical coherence tomography for cylinder scattering,” Proc. SPIE, 9541 95411F (2015). http://dx.doi.org/10.1117/12.2183748 PSISDG 0277-786X Google Scholar

47.

H. Tang, A. Gouldstone and C. DiMarzio, “Bubble raft-an optical phantom for analyzing artifacts in OCT imaging of lung,” in 41st Annual Northeast Biomedical Engineering Conf., 1 –2 (2015). http://dx.doi.org/10.1109/NEBEC.2015.7117200 Google Scholar

48.

A. Golabchi, “Corrections and improvements of lung imaging under optical coherence tomography (OCT),” (2012). Google Scholar

49.

C. I. Unglert et al., “Validation of two-dimensional and three-dimensional measurements of subpleural alveolar size parameters by optical coherence tomography,” J. Biomed. Opt., 17 (12), 126015 (2012). http://dx.doi.org/10.1117/1.JBO.17.12.126015 JBOPFO 1083-3668 Google Scholar

50.

S. Lacey and H. Wang, “Free space excitation of chaotic whispering gallery modes in a deformed glass microsphere,” in Technical Digest. Summaries of Papers Presented at the Quantum Electronics and Laser Science Conf., 26 (2002). http://dx.doi.org/10.1109/QELS.2002.1031038 Google Scholar

51.

J. Gao et al., “Observations of interior whispering gallery modes in asymmetric optical resonators with rational caustics,” Appl. Phys. Lett., 91 (18), 181101 (2007). http://dx.doi.org/10.1063/1.2800308 APPLAB 0003-6951 Google Scholar

52.

Y. Baryshnikov et al., “Whispering gallery modes inside asymmetric resonant cavities,” Phys. Rev. Lett., 93 (13), 133902 (2004). http://dx.doi.org/10.1103/PhysRevLett.93.133902 PRLTAO 0031-9007 Google Scholar

Biography

Thomas Brenner is a PhD student at the Institut für Lasertechnologien in der Medizin und Meßtechnik (ILM). He works in the Materials Optics and Imaging Group under the supervision of Prof. Dr. Alwin Kienle. He received his master’s degree from the University of Ulm and studied for one semester in Copenhagen, Denmark. His research interests include analytical and numerical solutions of Maxwell’s equations, optical coherence tomography, and light propagation in biological tissues.

Dominik Reitzle is a PhD student at the Institut für Lasertechnologien in der Medizin and Meßtechnik (ILM), Ulm, Germany, and part of the Materials Optics and Imaging Group at ILM. He studies physics at the University of Ulm.

Alwin Kienle is the vice director (science) of the Institut für Lasertechnologien in der Medizin und Meßtechnik (ILM), Ulm, Germany, and head of the Materials Optics and Imaging Group at ILM. In addition, he is a professor at the University of Ulm. He studied physics and received his doctoral and habilitation degrees from the University of Ulm. As postdoc, he worked with research groups in Hamilton, Canada, and Lausanne, Switzerland.

Citation Download Citation

Thomas Brenner, Dominik Reitzle, and Alwin Kienle "Optical coherence tomography images simulated with an analytical solution of Maxwell’s equations for cylinder scattering," Journal of Biomedical Optics 21(4), 045001 (1 April 2016). https://doi.org/10.1117/1.JBO.21.4.045001

Published: 1 April 2016

Access the abstract

JOURNAL ARTICLE
12 PAGES

DOWNLOAD PAPER SAVE TO MY LIBRARY

GET CITATION

CITATIONS

Cited by 16 scholarly publications.

Explore citations on Lens.org

RIGHTS & PERMISSIONS

Get copyright permission Get copyright permission on Copyright Marketplace

KEYWORDS

Optical coherence tomography

Scattering

Geometrical optics

Refractive index

Gaussian beams

Ray tracing

Biological research

1.

Introduction

1.1.

Theory

Eq. (1)

Eq. (2)

Fig. 1

Eq. (3)

Eq. (4)

Eq. (5)

Eq. (6)

Eq. (7)

Eq. (8)

Eq. (9)

Eq. (10)

Fig. 2

Eq. (11)

Eq. (12)

Eq. (13)

Eq. (14)

Eq. (15)

Fig. 3

Eq. (16)

Eq. (17)

2.

Calculation of an Optical Coherence Tomography Signal

Eq. (18)

Eq. (19)

Eq. (20)

Eq. (21)

Eq. (22)

Eq. (23)

Eq. (24)

Eq. (25)

Eq. (26)

Eq. (27)

3.

Results

3.1.

Time-Resolved Intensity for Plane Wave Scattering

Fig. 4

Fig. 5

3.2.

Optical Coherence Tomography A-Scan for Plane Wave Scattering

Fig. 6

Fig. 7

3.3.

Near Fields and Whispering Gallery Modes

Fig. 8

3.4.

Plane Wave Optical Coherence Tomography with the Debye Series

Fig. 9

3.5.

Gaussian Beam Optical Coherence Tomography

Fig. 10

Fig. 11

3.5.1.

Optical coherence tomography with an aperture

Fig. 12

3.5.2.

Optical coherence tomography simulation for a layered cylinder

Fig. 13

4.

Conclusions

Acknowledgments

References

Biography

Show All Keywords

Keywords/Phrases

Search In:

Publication Years