Spectroscopic Optical Coherence Tomography by Using Multiple Multipole Expansion

This paper presents a pre-processing method to remove multiple scattering artifacts in spectroscopic optical coherence tomography (SOCT) using time–frequency analysis approaches. The method uses a multiple multipole expansion approach to model the light fields in SOCT. It is shown that the multiple scattered fields can be characterized by higher order terms of the multiple multipole expansion. Hence, the multiple scattering artifact can thus be eliminated by applying the time–frequency transform on the SOCT measurements characterized by the lower order terms. Simulation and experimental results are presented to show the effectiveness of the proposed pre-processing method.


Introduction
Spectroscopic optical coherence tomography (SOCT) [1] extends the function of optical coherence tomography (OCT) to obtain the depth resolved spectra of samples [2].The spectra could provide invaluable information such as the health condition of the cells or tissues [3,4], authenticity of Chinese freshwater pearls [5], and mapping of microvascular hemoglobin [6].The heart of SOCT is the time-frequency analysis of the SOCT measurements.Various time-frequency analysis approaches such as the short-time Fourier transform method [4], the Wigner transform method [7], and the Dual window Fourier transform method [8] have been developed successfully and demonstrated capable of obtaining depth resolved spectra.In application of the time-frequency analysis, it is assumed that the SOCT measurement consists of fields scattered directly by the sample features of interest.However, SOCT measurement also consists of those fields that are scattered multiple times from the features [9].The multiple scattered fields reduce image contrast and make it more difficult to localize features or estimate scattered spectra, particularly at increased imaging depth [10].Hence, SOCT could gain benefits from combining the time-frequency analysis with some pre-or post-processing methods to mitigate the multiple scattering artifacts.
Post-processing methods such as spatial averaging or angular averaging have been proposed to mitigate the multiple scattering effect.Each of these methods have advantages and disadvantages.For example, the method proposed in [9] is able to distinguish multiple scattering from single scattering, but it is limited to the samples that are accessible from both sides.More recently, pre-processing methods based on signal correlation [11] or on a mixture of adaptive and computational optics [12] have also been developed, but these methods requires additional apparatus and they are not scalable to existing OCT setups.In this paper, we propose to model the light fields in the sample by a multiple multipole expansion (MMP) approach [13].The SOCT measurements are characterized by the positions and parameters of the multipoles.It is found that the MMP terms of higher orders are sensitive to the multiple scatterings.Thus, the multiple scattering artifacts can be removed if time-frequency analysis is applied only on those SOCT measurements that are fitted by the MMP terms of lower orders.Simulation and experimental studies show that the pre-processing is effective for any time-frequency based SOCT methods.
The rest of this paper is organized into the following sections.Section 2 formulates the proposed pre-processing method of using a multiple multipole expansion.Section 3 illustrates a simulation study where multiple scatterings are simulated and their effects on high order MMP terms are disclosed.In Section 4, the performance of the proposed method is demonstrated experimentally on a home-made phantom.

Spectroscopic OCT Using MMP Model
In the proposed method, each scatterer also representing a sub-surface feature of interest in SOCT is modeled as a multipole with boundaries formed by circles, semi-circles, connected circular arcs, or combinations of them, as shown in Figure 1.Specifically, we consider using M multipoles, each individual position is marked by, respectively, the radial distance r m and the azimuthal angle α m with reference to the center of the focusing lens in an OCT system.The light field detected at the focusing center can be expressed as the sum of scattered fields from the M multipoles as follows [13]: where J n is the Bessel function of order n, A n and B n are unknown expansion coefficients, and κ is 2π times the reciprocal of the wavelength.
Ignoring the self-correlation terms and DC terms due to the sample scattered spectrum and the mirror scattered spectrum, the depth resolved spectrum at the depth given by r m cos(α m ) can be obtained by where g(κ) is the spectrum of the broadband light source used in the SOCT system.The SOCT measurement can be expressed in the following compact format where M m (κ, r, α) is a vector with elements of 2g(κ)J 0 (κr) cos(j.α)and 2g(κ)J 0 (κr) sin(j.α)for j = 0, 1, . . ., N, and It is noted that the model has factored into multiple scatterings.As shown in Section 3, multiple scatterings are contained only in the MMP terms of higher order.Hence, by removing the fields characterized by the high order terms from the SOCT measurement, the multiple scattering artifacts is removed from the SOCT results.Now, the problem of the proposed pre-processing method can be stated as the one of solving unknown parameters r m , α m , A n , and B n from the available I(κ).Once these unknowns are obtained, the pre-processed SOCT measurements are obtained by the summation in Equation ( 3) over those multipoles of pre-set orders.
For easy calculation, we take spectral measurements at different κ i for i = 1, 2, . . ., L and (L > 2N + 2).L equations of 2N + 2 unknowns are obtained as where 4) can be solved by the following optimization: Considering noises in measurement, Equation ( 4) can be solved by the following mixed where denotes the noise level.The pre-processed SOCT measurements can be computed using the solution of Equation ( 5) with MMP terms of higher order removed.

Effect of Multiple Scattering on MMP Expansion
In this section, simulation studies are carried out to evaluate the contributions of multiple scatterings to the MMP expansion of SOCT measurements.We consider a sample with two spherical scatterers with diameters of 1 µm and 4.5 µm, respectively.The two scatterers are 10 µm apart in space, as shown in Figure 2. The multiple scattering process is simulated in the following way.Assume a single light pulse of 6 fs starts from L and travels through the scatterers sequentially, we calculate the amplitudes of light fields at various positions and time using finite-difference time-domain (FDTD) method.This method is used to simulate the field distribution as it neither needs to assume the geometry of the scatterers nor the extent of multiple scattering.Thus, the method is able to consider the contributions of multiple scattering effectively in the calculated field distribution.We used commercial FDTD software (Lumerical) working in temporal mode to simulate the light fields over a time duration of 300 fs.Other simulation parameters included setting the maximum mesh size to 50 nm to ensure that the mesh size is smaller than the smallest wavelength used, and using perfectly matched layer (PML) boundary condition to terminate the field.Then, the simulation was run and the amplitudes at locations marked by D1, D2, and D3 are recorded at different time, as shown in Figure 3a.
First, we considered the scattering events that contribute to the single scattered electric field received by D1.The light pulse first passes D2 and gives the response marked by C in Figure 3a, and then it encounters S1.One portion of the light is backscattered by S1 towards D1 and gives the response marked by A in Figure 3a.The other portion continues travelling towards D3 and gives a response marked by E in Figure 3a.After the light encounters S2, a portion of it carries on its forward propagation and does not return, whereas a backscattered portion travels towards D2.The backscattered portion propagates towards detector D2 and D1 and produces responses marked by D and B in Figure 3a, respectively.Beyond this instance, the scatterings that take place are multiple scattering fields and when the light pulse returns to D1, it generates the multiple scattering fields in SOCT measurement.Earlier when S1 is first illuminated by the light pulse, S1 scatters a portion of the field towards D3 to give the response F in Figure 3a.Likewise, the field described above that generates the response marked by D in Figure 3a is backscattered by S1 towards S2 to generate the response G in Figure 3a.This process of scattering to and fro between spheres S1 and S2 goes on several times and gives the total field received at D1 containing both single and multiple scatterings.Thus, by filtering the response at D1 at different time durations, we can calculate, respectively, the spectrum due to electric field from single scattering and the spectrum due to the total received field including the multiple scattering portion.Furthermore, we can obtain the corresponding spectra.The difference between the two spectra is plotted in Figure 3b.For comparison, we repeated the simulation for the case when two scatterers are separated further by S = 35 µm.The corresponding spectrum difference is included in the figure .The comparison clearly shows that the spectrum difference due to multiple scatterings becomes more significant when the two scatterers are closer.It is the spectrum difference that contributes the multiple scattering artifacts in SOCT.
Next, we show that MMP terms of higher order fit well the multiple scattering fields in SOCT measurements.We solve the optimization problem in Equation ( 5) on the SOCT measurements consisting of multiple scatterings using N = 3 and N = 5, respectively.We take r m = [0 µm, 10 µm, 20 µm ... 1000 µm] and α m = [−0.05,−0.03...0.05].The resultant SOCT measurements are shown in Figure 3c.It is evident that the MMP terms of higher order can be used to characterize effectively the multiple scattering artifacts in the SOCT measurements.Hence, it is possible to filter the multiple scattering fields in SOCT measurements by limiting the order of the multiple multipole expansion.The filtered SOCT measurements from Equation ( 5) can then be used to calculate the depth resolved spectroscopy of SOCT.

Setup
The performance of the proposed method was evaluated in an experiment by imaging a homemade phantom with an OCT setup, as shown in Figure 4.The OCT setup is a typical Michelson interferometer.Low coherence light comes from a superluminescent diode (QSDM-810-2, QPhotonics, Ann Arbor, USA) with center wavelength of 813 nm and FWHM linewidth of 19 nm.The light enters the interferometer via a collimator (NA = 0.26, Thorlabs, New Jersey, USA) and it is split into two paths by a 50:50 beam splitter (25 mm near-infrared non-polarizing cube beam-splitter, Edmund Optics, Singapore).The split light is then focused, respectively, on the reference mirror and the sample through near infrared achromatic lens (40 mm focal length, Edmund Optics, Singapore).The light reflected from the mirror and from the sample interferes again at the beam splitter.The interference light goes through a collimation module consisting of two identical achromatic lens with a 25 µm pinhole placed in the middle at the common focal point of two lenses.The collimated light is dispersed and projected to a charged coupled device (CCD, Dalsa 1M60, Waterloo, Canada) through a grating (800 lines per mm, Edmund Optics, Singapore), and an achromatic lens.A separate calibration shows the axial and lateral resolutions of the OCT setup are 20 µm and 10 µm, respectively.

Phantom
The homemade phantom was made by coating a film on a microscope glass slide.The film contains a mixture of microspheres with diameters of 2.0 µm and 4.5 µm microspheres.The film was made in the following steps.A thin layer of glue is first applied onto the surface of the microscope glass slide.Then, a 5 µL droplet of solution consisting of a mixture of polystyrene microspheres of 2.0 µm and 4.5 µm diameters is dropped onto the microscope glass slide.The droplet of solution is spread out by gently pressing a cover glass with a thin layer of glue facing down onto the surface of the microscope glass slide.A microscopy picture of the phantom is shown in Figure 5.

Imaging Procedures
In the experiment, we imaged the film from its side rather than its surface, as shown in Figure 5.We aligned the phantom along the Z-direction (not the X-direction) to the light beam of a typical SOCT setup.The spectra of A-scans in a cross-section of the sample were collected in a process of scanning the sample with a motorized translational stage (25 mm translational stage, Zaber, Vancouver, Canada).The scanning was done in steps of 1 µm along the cross-section of interest.At every scan position, an A-scan spectrum was captured.A total of 500 A-scans were obtained for three cross-sections: the first cross-section contains only 2 µm diameter microspheres, the second one contains only 4.5 µm diameter microspheres, and the third one contains both 2 µm and 4.5 µm diameter microspheres.
Then, the proposed method was applied to process the A-scan spectra through solving the optimization problem in Equation ( 6) on the SOCT measurements by using N = 3, = 0.056, r m = [0 µm , 10 µm, 20 µm ... 1000 µm] and α m = [−0.05,−0.03...0.05].This solution provided us the estimates of the parameters for the scatterers, namely A n , B n , α m and r m .Then, we used these parameters to reconstruct the spectra reflected from each scatterer by using Equation (2).After this pre-processing, a short-time Fourier transform was performed on the processed A-scan spectra and we obtained a SOCT image by calculating the contrast measure described in the section below.

Visualization
For better visualization of the SOCT image, a new contrast measure is introduced here, and it is defined by A 2 1 + B 2 1 .The measure reflects the scatterer size because both A 1 and B 1 in the MMP expansion are directly proportional to the radius of curvature of scatterer boundary.To display the SOCT image of the phantom, the contrast measure is rounded about two values that, respectively, represent the 2 µm and 4.5 µm diameter microspheres.The contrast is coded into the green and red channels.As a result, the image obtained by the proposed method are obtained.We also process the SOCT measurement by short-time Fourier transform method.The SOCT images for the cross-section containing only 2 µm diameter microspheres and the one containing only 4.5 µm diameter microspheres are shown in Figures 6 and 7, respectively.The corresponding microscope images are also given for reference.Finally, we process in the same manner the SOCT measurement obtained from the cross-section containing a mixture of 2 µm and 4.5 µm diameter microspheres.We also process the SOCT measurements by the standard FD-OCT method and the short-time Fourier transform method.The obtained images are included in Figure 8 for comparison.Comparing the SOCT images obtained by short-time Fourier transform method and the proposed method in Figures 6 and 7, we can see that the images obtained by short-time Fourier transform displayed artifacts due to multiple scatterings when two scatterers that are located close together, whereas the SOCT image from the proposed method does not have this artifact.
Comparing the images in Figure 8 draws two observations.First, the proposed method can locate the scatterers of 4.5 µm and 2 µm, as shown in Figure 5. Second, compared with the one obtained by short-time Fourier transform, the proposed method is able to eliminate the artifacts due to multiple scatterings.This can be attributed to the ability of the multiple multipole expansion to differentiate the multiple scatterings from single scatterings.

Conclusions
In summary, this paper uses a multiple multipole expansion method to model the light fields in a sample under SOCT imaging.Since OCT and SOCT setup uses the same set of measurement, the pre-processing method is also applicable to OCT.Using the proposed method, multiple scatterings from features in the sample can be effectively modeled by a parameterized optimization problem.With the parameters obtained, multiple scattering fields can be filtered out from the SOCT measurements before applying routine time-frequency analysis of SOCT.The multiple scattering artifacts are removed in the SOCT imaging results.Experiments on a calibrated phantom has exemplified that the proposed method can improve the image quality.The challenges to integrate this method to existing time-frequency analysis methods are with obtaining a suitable scattering model of the features and minimizing the spectral distortion that could arise from dispersion.In the case of industrial samples, suitable scattering models would be determined by characterizing the scattered spectra and dispersion of the different materials used to manufacture the samples.

Figure 1 .
Figure 1.Schematic illustration of an A-scan.

Figure 2 .Figure 3 .
Figure 2. Schematic illustration of an A-scan showing the locations of optical axis and reference plane.The reference plane has zero optical path difference (OPD) with the reference mirror in OCT setup.

Figure 4 .
Figure 4. Schematic illustration of the OCT setup and the components used.

Figure 5 .
Figure 5. Cross-sectional view of the phantom sample imaged by Olympus upright microscope with 100× magnification (the while bar denotes 10 µm).

Figure 6 .
Figure 6.Imaging of the cross-section with 2 µm microspheres only: (Left) SOCT image of the cross-section obtained by short-time Fourier transform method; (Middle) SOCT image of the cross-section obtained by the proposed method; and (Right) cross-sectional view of the area under investigation.SOCT images are normalized and displayed in logarithmic scale.

Figure 7 .Figure 8 .
Figure 7. Imaging of the cross-section with 4.5 µm microspheres only: (Left) SOCT image of the cross-section obtained by short-time Fourier transform method; (Middle) SOCT image of the cross-section obtained by the proposed method; and (Right) cross-sectional view of the area under investigation.SOCT images are normalized and displayed in logarithmic scale.In the microscope image, there are small air bubbles that are not in the imaging XZ plane and they do not contribute our SOCT measurement.