Optical coherence tomography (OCT) is a potentially noninvasive sensing optical imaging system for various clinical diagnosis applications, and has a lot of developments particularly in the past decade.1,2 As an emerging imaging system, OCT can provide high-resolution axial structure images inside tissue from two- (A-scan) to three-dimensional (A-scan and B-scan).3 It is an interferometry-based optical technique, collecting the backscattered light from the sample to further extract the optical characteristic details within the tissue. Currently in this field, precise analysis of optical properties from the inhomogeneous and dense tissue remains challenging work.4 Different models of OCT have been established to obtain the optical properties of samples with further consideration of multiple scattering effect.5,6 Among these works, Thrane et al. were the first to propose the use of modulated interferometer signals of an OCT system, which can be expressed analytically based on the extended Huygens-Fresnel (EHF) principle, and to further consider the effect of multiple scattering.7 In addition, EHF theory has commonly been utilized to numerically quantify optical properties from A-scan signals of OCT images. However, in the EHF theory, there are two main flexible parameters highly correlated with the slope of OCT signals, i.e., the scattering coefficient () and the anisotropy factor (). Koach et al.8 reported in 2011 that the scattering phase function in a backward direction could be determined as the ratio of the backscattering coefficient to the scattering coefficient (only single scattering) detected from the CCD, but without the consideration of multiple scatterings. Then in 2013, Nguyen et al.9 further improved the OCT system by introducing a new setup of the transmission OCT system, in order to detect transmitted signals, and they introduced the structure factor (Percus–Yevick theory) to show the concentration-dependent scattering coefficients. In 2015, Almasian et al.10,11 proposed that a concentration-dependent scattering coefficient of the sample could be obtained from the Mie theory by incorporating the structure factor of Percus–Yevick theory in the calculation procedure. In this study, we propose a new analysis flow by combining the OCT system with the finite-difference-time-domain (FDTD) simulation method, to quantify the scattering coefficient of mesoporous beads after introducing factors into EHF theory. The benefit of combining with the FDTD method is that the characteristics of objects can be flexibly modified, including the shape (not only sphere), the refractive index, and the arrangement for inhomogeneous morphology. For comparison to the Percus–Yevick theory mentioned above, the structure factor can only hypothesize for a set of two hard spheres with fixed diameter size at a specific condition.12 In addition, differing from the algorithm of calculating the radial correlation function (radial distribution function) between adjacent particles in P-Y theory, the near-field electrical field of the FDTD method can briefly demonstrate and show the light interaction when the object is scattered by incident light, and the far-field scattering pattern can be further determined.
First, the anisotropy factor represents the intrinsic scattering characteristic properties (angular scattering pattern). The multiple scattering phenomena depending on the interparticle spacing (IPS) will be discussed in Sec. 3. From this method, the scattering phase function  can be simulated to predict the corresponding anisotropy factor of mesoporous beads samples with different diameters varied from 20, 150, 300, to 500 nm, respectively. In Sec. 4, the A-scan profiles are extracted by signal processing procedures from the interference patterns and then fitted by the EHF theory with consideration for the corresponding anisotropy factor , which is introduced to the EHF model to fit the scattering coefficient of each mesoporous bead sample. Also, different calculations of the anisotropy factor between the Mie scattering model and the FDTD method are compared.
Methods and Experimental Setup
Since the development of the FDTD method,13 it has been utilized to solve Maxwell’s equation numerically in the time domain to describe the complex behavior of light interaction and widely applied in solving the scattering issue in biological fields.14–16 In this study, we have utilized the R-Soft software with fullwave mode based on the FDTD algorithm. The FDTD algorithm is the computation process for implementing time and space derivatives of the electric and magnetic fields in Yee’s cell.17 Solve the Maxwell’s curl equations by first discretizing the equations with central differences in time and space and these equations can be further numerically solved by the software. To estimate the optical scattering properties of beads with different diameters and embedded in thin-films, we align the scattering direction along with the line connecting the beads, and they are separated by the IPS, which varies from 0 to 1000 nm in an interval of 250 nm. The model used for FDTD is shown in Fig. 1(a): the wavelength of the incident plane wave is , while the refractive index of mesoporous beads is 18,19 and assumed to be homogeneous and isotropic in distribution. The polarization of the incident wave is set as TE polarized. The grid size of calculation is set as 0.005 nm, time step is 0.005 s, and the incident wave is a plane wave in the enclosed area, respectively. With a perfectly matched layer (PML), it is composed of several grid points as the edge of the domain for the enclosed area, which can absorb the incident lights and result in no reflection.
In order to obtain the far-field scattering patterns, the virtual box with appropriate PML boundary is set to enclose the scattering objects (mesoporous beads) and launch a plane wave inside an enclosed area, which is shown in Fig. 1(a). In the meantime, there are only scattering fields from objects propagating beyond the enclosed launch boundary, and unscattered light will be absorbed at the PML boundary. For obtaining the scattering phase function, the corresponding far field can be calculated using the near-field to far-field transformation from the simulation data in software.20
Preparation of Mesoporous TiO2 Beads Samples
Mesoporous bead samples with average diameters () of 20, 150, 300, and 500 nm were prepared by the sol–gel process. Initially, the powder weight percentage is 13 wt. % (US Research Nanomaterials, Inc.) mixed with different percentages of isopropyl alcohol and DI water. Adding polyethylene glycol to the mixture, the pastes were obtained after more than 1 week of stirring, and the pastes prepared on glass substrates were heated by an oven at 450°C for 30 min to evaporate the solvent. In the process, the were assumed to be homogeneously distributed. The final volume fraction of bead samples was fixed as 0.03 for beads with diameters varied from 20 to 500 nm, respectively. Afterward, the samples of mesoporous bead films on glass substrates were obtained.
For a suspension of colloidal system, the samples of beads with volume fraction of 0.03 for different diameters are assumed to be homogeneously distributed; the schematic diagram of beads in a unit cell () for different diameters from 20 to 500 nm is shown in Fig. 1(b). The IPS parameter is defined as the distance from the surface of one particle to the other [the line of connection between the two through the center of each bead, as shown in Fig. 4(a)], which changes according to the volume fraction in the colloidal system.21 The maximum volume fraction is with samples, while is the volume fraction of bead samples, and is the radius of beads. With the volume fraction of 0.03, the values of IPS are 25, 187, 373, and 623 nm for select diameters between 20 and 500 nm.
Extended Huygens-Fresnel Theory (Multiple Scattering) and Signal Analysis
First, as an assistant method, the FDTD simulation process was utilized to calculate the electromagnetic near field and then transformed to a far-field scattering pattern, which is also called the phase function , for different IPS between mesoporous beads with different diameters from 20, 150, 300, and 500 nm. Here, the far-field is indicated as the probability density function of photons for every scattering angle (),22 and the anisotropy factor is the integral of phase function over all scattering angles, which can be expressed as
Second, based on EHF theory, the amplitude of experimental OCT signals versus depth can be modulated as the function of two main flexible parameters for and . Once the multiple scattering effect appears, the right-hand side of point spread function (PSF) can be changed by slope variation at a larger axial depth as compared to that of a signal with only single scattering effect.23 In practice, the most common method to date to obtain the optical properties has been to analyze the axial A-scan signals based on EHF model at the right-hand side of PSF with increasing depth. Except for the fixed optical system parameters, the resulting heterodyne signals [Eqs. (3) and (4)] in depth can be mainly modulated by two competing factors, the anisotropy factor (, is average scattering angle) and the scattering coefficients ().23 In addition, the fixed parameters in the system include the power and measured at the reference and sample arms; the quantum conversion efficiency of the CCD; the backscattering cross-section , and the irradiance radii and for the absence and presence of scattering; the refractive index for medium ; the radius of the sample beam at the focal plane ; and the focal length of objective lens , as follows:
In the EHF theory, it should be noted that it is not easy to predict the optical coefficients of and at the same time due to the fact that they are competing parameters. Therefore, using the other method to obtain one of them is a general solution to approach. After obtaining the of mesoporous beads with different diameter sizes for various IPS from the FDTD method, introducing the value of for specific beads samples, the corresponding scattering coefficient can be further derived for the PSF of A-scan signal based on EHF theory.
Optical Coherence Tomography Setup
The schematic diagram of the SD-OCT system using an 853-nm super-luminescent diode as the source is shown in Fig. 2(a). A single mode fiber (Thorlabs 780HP) is used to direct the light into the reference and sample arms. This OCT system has a numerical aperture of 0.042, with an objective lens focal length of 17.9 mm. The interferometer signal strikes a transmission grating (1200 lines/mm) and the objective lens then collects and focuses the diffracted light onto the CCD (e2V, 2048 pixels, 28-kHz line rate). Each A-scan profile (in depth) was averaged over 1000 measurements and subtracted by the reference arm intensity, to obtain the interferometer signal without the DC term.
Scattering Phase Function of Single TiO2 Bead Simulated by FDTD
The optical behavior of a single TiO2 bead with different diameters has to be confirmed prior to the simulation of multiple beads, with the FDTD method analysis. Figure 3 shows the plots of angular-dependent scattering patterns of a single bead with diameters of 20, 150, 300, and 500 nm, obtained by using the FDTD method, and is compared with the results from the Mie scattering model.24–26 These comparison results show that both methods are consistent and verify that the settings in the FDTD method are appropriate. The resultant values of obtained for single bead with diameters from 20, 150, 300, and 500 nm are 0.003, 0.174, 0.525, and 0.327, respectively.
Interparticle Spacing Effects of Near-Field Scattering Patterns
To study the multiple scattering effects that occur within high density tissues, multiple beads (for smaller IPS in Sec. 2.2) are simulated, with near-field scattering patterns as shown in Fig. 4(b). In analyzing actual OCT signals, the presence of multiple scattering makes it difficult to acquire the corresponding anisotropy factors .27,28 We simplify the scatter geometry to in-line configuration of beads separated by different IPS in order to investigate the scattering phenomenon. In Fig. 4(a), the electrical field patterns for in-line beads with a diameter of 300 nm are shown; the IPS is varied from 0 to 1000 nm in intervals of 250 nm. The light propagates in the upward direction and is incident onto the beads from underneath, as shown in Fig. 4. According to the results of far-field scattering patterns, except for the main lobe that resulted from the Mie scattering effect, the presence of side-lobes was found by considering the effect of multiple scattering. This needs to be noted when the IPS is increased.
Interparticle Distance Effects-Dependent Phase Function
To approach the realistic [Eq. (1)], the far-field angular phase function  of in-line beads was further simulated using the FDTD method. Figure 5 shows the effect of IPS on the angular scattering patterns for mesoporous beads with diameters of 20, 150, 300, and 500 nm. In each case, we observe adjusting the value of (), where is the period for in-line beads, which can be defined as the sum of IPS and the diameter of bead. Changing of induces variation on the phase function  [Figs. 5(a)–5(d)], including the contribution ratio between backward and forward scatterings, and the appearance of side-lobes. Apart from the characteristic nature of Mie scattering , shown for each diameter size of bead, the incident wave also leads to scattering field interference through the propagation path.
For the 150-nm-diameter beads, the appearance of side-lobes can be found as the IPS is increased beyond 750 nm. When IPS is smaller than 500 nm, the scattering patterns are similar to Mie scattering ones; these can be maintained except for the main-lobe energy leakage with increasing IPS, without diffraction associated side-lobes. The side-lobes can be detected at scattering angles of 52 deg and 128 deg. The two lobes are symmetric with respect to the main lobe, and when reaching the condition of , the interference phenomenon induced by phase delay between adjacent beads [see Fig. 4(b)] produces a grating diffraction effect.29,30 The interference pattern can be likened to a multiple-slit interference pattern from the geometric arrangement of in-line beads [Fig. 4(a)], by analogizing the diameter of beads as the width of the slit, and the bead periodicity as the grating period. The detected intensity was simulated by the function which can be expressed31 as , where is the intensity of incident wave, is the signal of single slit diffraction, where , and is the signal of the wave interference by phase delay, where . The minor side-lobes can be approximated by using the above equations at the angles of 48.1 deg and 131.4 deg with some small, but acceptable within margin, amount of errors, as compared with the data from the FDTD method. As shown in Figs. 6(a)–6(d), the factors of beads with different diameters from 20 to 500 nm versus the IPS from 0 to 1000 nm were simulated and calculated by the FDTD method, which integrated the phase function overall the scattering angles. In addition, the curves as the function of IPS for these samples were fitted with the B-spline function.32 With the experimental volume fraction condition of 0.03, the corresponding values extracted from the curves in Figs. 6(a)–6(d) are 0.0025, 0.1685, 0.3181, and 0.3866 for bead diameters of 20, 150, 300, and 500 nm, respectively. Introducing the specific IPS values of 25, 187, 373, and 623 nm for diameters between 20 and 500 nm with the volume fraction of 0.03 into FDTD method, shows the slight difference within 5% error deviation as compared to the curve of as simulated in Figs. 6(a)–6(d). Therefore, the curves of versus IPS can be further used in scattering analysis, and the scattering coefficients () of different samples can be obtained by inserting the factors into EHF model. Therefore, the trend of IPS-dependent for different diameters of beads can be determined. As illustrated in Figs. 6(a) and 6(b), the with smaller diameters for 20 and 150 nm show factor results similar to those from the Mie scattering effect, which implies that the variation of IPS affects weakly for smaller beads. However, when increasing the size of up to 300 and 500 nm, the IPS-dependent factor deviates from the Mie scattering effect due to the diffraction and interference phenomena between adjacent beads. With this, an approximation of the anisotropy factor () of each bead sample can be made using the OCT system, and can be further introduced into the EHF theory to find the most fitted scattering coefficient ().
Extended Huygens-Fresnel Theory
After verifying the value of for each of the bead samples, the corresponding can be further fitted by the EHF theory. Figure 7(a) shows a plot of OCT A-scans for the experimental and fitted curves of the mesoporous sample with 300 nm diameter. With the fitted results, the slope of A-scan signals is dependent on the optical parameters () introduced by the EHF model, simulated by fixing the factor from the Mie scattering model and FDTD method. In Fig. 7(b), the scattering coefficient fitted for different bead diameters is shown. When increasing the diameter from 20 to 500 nm, the calculated scattering coefficient is , , , and , with the factors calculated from the FDTD method. In comparison, the calculated scattering coefficient is , , , and , respectively, when the factors are obtained by using the Mie scattering model. In Table 1, the calculated correlation factors (), between the experimental results and the fitted OCT results, show good agreement (), verifying the accuracy of the obtained IPS-dependent with factors for each bead sample. The can be calculated from the mean and stand deviation of the experiment, with fitted A-scan signals to estimate the linear dependence between these two signals.33 The highest can be calculated for different bead diameters of 20, 150, 300, and 500 nm, and the results are 0.9431, 0.9617, 0.9721, and 0.9241, respectively. This shows that, once the for specific IPS can be accurately determined, the corresponding scattering coefficient can be more accurately quantified with reduced error value from 0.6 to 0.3. Also from the results, the fitted can be observed to be overestimated, when the bead diameter is increased up to 300 nm, for which the multiple scattering effect becomes more pronounced along with the increased bead diameter. The oscillation behavior of scattering coefficient for different bead diameters from 100 to 300 nm is shown in Table 1 and Fig. 7(b) while the fitting is based on FDTD method. Due to the limitations of the experiment, the scattering coefficient trend can only be observed within the narrow range of variable bead sizes. However, similar results can also be found in other research.34 In addition, the fit can be obtained when factors calculated from the FDTD method show the smaller error values decreased from 0.6 to 0.3. The results of lower as compared with are consistent with the findings from other research,9,11 which had showed that the presence of the multiple scattering effect, due to larger diameter size and larger volume fraction, can induce the deviation from the calculated based on Mie scattering. Therefore, in these cases, the FDTD method is proved to be useful for accurately quantifying the optical properties for multiple scattering.
Correlation factor of scattering coefficient analysis for both FDTD method and Mie scattering model.
|TiO2 beads (D)||θrms (Mie)||μs (Mie)||Correlation factor|
|20 nm||1.5678||13.5 ± 0.6||0.9429|
|150 nm||1.3961||16.5 ± 0.6||0.9559|
|300 nm||1.0183||20.2 ± 0.6||0.9613|
|500 nm||1.2378||18.5 ± 0.6||0.9243|
|TiO2 beads (D)||(FDTD)||(FDTD)||Correlation factor|
|20 nm||1.5683||12.9 ± 0.3||0.9431|
|150 nm||1.4015||20.8 ± 0.3||0.9617|
|300 nm||1.2472||14.2 ± 0.3||0.9721|
|500 nm||1.1739||19.8 ± 0.3||0.9241|
Alternatively, one similar analysis was shown using the structure factor in Percus–Yevic theory (2015).11 It interpreted the spherical particles illuminated by scattering light that induce the nonlinearly scattering cross section, to show that the oscillatory behavior depends on the volume fraction. With scatter of larger diameter, despite taking into account the concentration dependence, the optical properties induced by multiple scattering cannot be completely described. Our study based on the FDTD method can be used to directly observe the far-field scattering pattern for different samples. In addition, with methods based on Monte Carlo simulation, it is difficult to perform proper algorithm modeling.35,36 Therefore, in this work, the OCT and FDTD combination method stands out for its flexibility in simulating complex structures (cellular morphology, cancer cells, etc.)37,38 to directly acquire the original phase function for further expansions in biological applications, and as opposed to making indirect predictions in the other methods mentioned above.
In this work, we successfully demonstrated that combining the OCT system with the FDTD method can be used to obtain the factor of high density tissues in the presence of multiple scatterings. Combining with the FDTD method, the light scattering interactions for near- and far-field patterns between adjacent beads can be determined to identify the anisotropy factor, which also can be utilized for other scatters with different shapes or arrangements. The can be further extracted by using a fitting process based on the EHF theory. The increase in the sample density is accompanied by smaller IPS, and the near-field scattering patterns show that the photons interact and interfere between adjacent scatters through propagating the sample. A new analysis flow to accurately quantify the optical properties of the samples through the use of OCT and FDTD combination method is flexible and has future potential in the biomedical field.
The authors have no relevant financial interests in this article and no conflicts of interest to disclose.
The authors gratefully acknowledge the support from the Ministry of Science and Technology of Taiwan under Contract Nos. MOST 106-2221-E-002-155-MY3, 105-2221-E-002-121, and MOST 107-3113-E-155-001-CC2 and National Taiwan University under the Aim for Top University Project 105R8908.
Ling-Hsuan Tsai is a PhD candidate at the Graduate Institute of Photonics and Optoelectronics, College of Electrical Engineering and Computer Science, National Taiwan University. She received her BS degree from Feng Chia University and her MS degree from the National Sun Yat-Sen University in 2008 and 2010, respectively. Her current research interests are utilizing optical coherence tomography system to analyze scattering optical properties information from tissue and further obtain accurate scattering parameters.
Po Nien Yang is a PhD candidate at the Graduate Institute of Photonics and Optoelectronics, College of Electrical Engineering and Computer Science, National Taiwan University. He received his BS degree from I-Shuo University and his MS degree from National Taiwan Normal University. His research interests include using optical coherence tomography systems to study the physics of photon scattering, theoretical high energy physics phenomenology, and quantum field theory.
Chung-Chih Wu received his PhD from Princeton University in 1997. Currently, he is a distinguished professor of Department of Electrical Engineering, Graduate Institute of Photonics and optoelectronics, and Graduate Institute of Electronics Engineering, National Taiwan University.
Hoang Yan Lin received his BS and PhD degrees from the Electrical Engineering Department and the Graduate Institute of Electrical Engineering at National Taiwan University in 1987 and 1993, respectively. His research interests are design of optical components and integration of optical systems. He is the representative for the Taipei section of the Institute of Electronics, Information, and Communication Engineers, and a member of the Optical Society of America, the Society for Information Display, and the 3D Interaction & Display Association.