1.

Introduction

Most biological tissues are optically anisotropic turbid media, containing fibrous microstructures and birefringence optical polarization effects.^1–3 It is important to mimic and interpret the structural and optical anisotropies using polarized scattering models or phantoms for tissue characterization.^4–7 There have been many polarized scattering models for biological tissues.⁸ Wang et al. presented a sphere birefringence model, which contains spherical particles randomly suspended in linearly birefringent media.⁹ He and Yun et al. proposed a sphere-cylinder scattering model (SCSM) to characterize the anisotropic scattering property of fibrous tissues, such as skeletal muscles.^10,11 Furthermore, Du et al. extended the SCSM to the sphere-cylinder birefringence model (SCBM) by introducing birefringence to the medium between scatterers.¹² In a previous study, we examined the tissue anisotropy contributed by both the scattering of cylindrical scatterers and the birefringent media using Mueller matrix polar decomposition.¹³ Lately, it has been noticed that there exists intrinsic birefringence on intracellular microtubules, microfilaments, and other fibrous scatterers in tissues.¹⁴ Recently, in a study on the microstructures of textiles, Peng et al. proposed a birefringent cylinder scattering model, in which the textile fibers have different refractive indices in the axial and radial directions; a corresponding Monte Carlo simulation was also reported for textile fibers randomly distributed in the $xy$ plane.¹⁵

In this article, we extend further the SCBM to include the intrinsic birefringence on the cylindrical scatterers of arbitrary spatial distributions. We also developed a new Monte Carlo simulation program to calculate the effects of well-ordered birefringent fibers to mimic the structure of biological tissues.¹⁶ After verifying the validity of the new anisotropic tissue model and the Monte Carlo simulation program, we then analyzed the contributions of the well-ordered birefringent cylinders to the retardance obtained from Mueller matrix polar decomposition. We present the comparison of polarization characterization for the birefringent cylinders model and the model of nonbirefringent cylinders with birefringence medium between them. The study will help us to understand and mimic the polarization behavior of photons in complicated tissues and explain the optical anisotropy and characteristic features of tissues based on polarized photon scattering measurements.

2.

Theory

2.1.

Polarized Photon Scattering at a Inhomogeneous Birefringent Cylinder

We have developed a SCSM and a SCBM to describe tissue anisotropy.^11,12 We assumed that the small spheres and infinite cylinders are embedded in the medium of linear birefringence. When photons transport through the media, they will be scattered by the spherical and cylindrical scatterers. We extended the SCBM by replacing the nonbirefringent cylinders with birefringent ones. Peng’s calculation of single scattering for birefringent cylinders is referenced in this paper, but we corrected a mistake in Eq. (9), and the dielectric coefficient and permeability are replaced with refractive index, consistent with our previous notations.¹¹

We can establish a coordinate system shown below in Fig. 1. The refractive index of birefringent cylinders can be written as axisymmetric tensors

Eq. (1)

$$n̿=n_{\perp }^2\widehat{x}\widehat{x}+n_{\perp }^2\widehat{y}\widehat{y}+n_{||}^2\widehat{z}\widehat{z},$$

where $\widehat{x}$, $\widehat{y}$, and $\widehat{z}$ are the transverse and axial unit vectors.

Fig. 1

Schema of the single light scattering by infinite cylinder.

Then, we can write the internal electromagnetic field quantities in the form

Eq. (2)

$$\begin{gathered} \overrightarrow{E}={\overrightarrow{E}}_t+\widehat{z}{E}_z \\ \overrightarrow{H}={\overrightarrow{H}}_t+\widehat{z}{H}_z\mathrm{.} \end{gathered}$$

Subscript $t$ expresses the transverse plane, which is perpendicular to the $z$ axis. According to these, Maxwell’s equation can be expressed as^17,18

Eq. (3)

$$D_1({\nabla }_t)·H_z+D_3({\nabla }_t)·E_z=0D_2({\nabla }_t)·E_z-D_4({\nabla }_t)·H_z=0,$$

where $D_1({\nabla }_t)$, $D_2({\nabla }_t)$, $D_3({\nabla }_t)$, and $D_4({\nabla }_t)$ are quadratic functions of ${\nabla }_t$, and ${\nabla }_t$ is a differentiating operator. By using plane wave expansion method, we can solve Eq. (3) to get internal electromagnetic field quantities $E_z$ and $H_z$, and get ${\overrightarrow{E}}_t$ and ${\overrightarrow{H}}_t$. Using the continuity conditions on the boundary, we can get scattered field ${\overrightarrow{E}}_s$ (parallel component ${\overrightarrow{E}}_{||s}$ and perpendicular component ${\overrightarrow{E}}_{\perp s}$)

Eq. (4)

$$\left(\begin{array}{l}{E}_{||s}\\ E_{\perp s}\end{array}\right)=\frac{1}{\sin \zeta }\sqrt{\frac{2}{\pi k\rho \sin \zeta }}{e}^{j\frac{\pi }{4}}{e}^{-jk(\rho \sin \zeta +z\cos \zeta )}\left(\begin{array}{cc}{J}_1& J_4\\ J_3& J_2\end{array}\right)\left(\begin{array}{l}{E}_{||i}\\ E_{\perp i}\end{array}\right),$$

Eq. (5)

$$\{\begin{array}{l}{J}_1=\sum _{-\infty }^{\infty }{b}_{nI}{e}^{-\mathrm{in}\mathrm{\Theta }}=b_{0I}+2\sum _{n=1}^{\infty }{b}_{nI}\cos (n\mathrm{\Theta })\\ J_2=\sum _{-\infty }^{\infty }{a}_{nII}{e}^{-\mathrm{in}\mathrm{\Theta }}=a_{0II}+2\sum _{n=1}^{\infty }{a}_{nII}\cos (n\mathrm{\Theta })\\ J_3=\sum _{-\infty }^{\infty }{a}_{nI}{e}^{-\mathrm{in}\mathrm{\Theta }}=-2i\sum _{n=1}^{\infty }{a}_{nI}\sin (n\mathrm{\Theta })\\ J_4=\sum _{-\infty }^{\infty }{b}_{nII}{e}^{-\mathrm{in}\mathrm{\Theta }}=-2i\sum _{n=1}^{\infty }{b}_{nII}\sin (n\mathrm{\Theta })\end{array},$$

Eq. (6)

$$\{\begin{array}{l}{a}_{nI}=-\frac{2}{\pi \tau }{F}_n^{(2)}{J}_n({\xi }_{+})/\{\delta ·H_n^{(2)}(\tau )\}\\ b_{nI}=\{\frac{2}{\pi \tau \eta }{F}_n^{(1)}{J}_n({\xi }_{-})/\delta -\sin \zeta J_n(\tau )\}/H_n^{(2)}(\tau )\\ a_{nII}=\{-\sin \zeta J_n(\tau )-\frac{2}{\pi \tau }\eta F_n^{(4)}{J}_n({\xi }_{+})/\delta \}/H_n^{(2)}(\tau )\\ b_{nII}=\frac{2}{\pi \tau }{F}_n^{(3)}{J}_n({\xi }_{-})/\{\delta ·H_n^{(2)}(\tau )\}\\ \delta =F_n^{(1)}·F_n^{(4)}-F_n^{(2)}·F_n^{(3)}\end{array},$$

Eq. (7)

$$\{\begin{array}{l}{F}_n^{(1)}=-j\frac{\eta }{\sin \zeta }{H}_n^{(2)\prime }(\tau )J_n({\xi }_{+})+g_n^{(1)}{J}_n^{\prime }({\xi }_{+})\\ F_n^{(2)}=\frac{\cos \zeta (n_0^2-n_{\perp }^2)}{ak{\sin }^2\zeta (n_{\perp }^2-n_0^2{\cos }^2\zeta )}n{H}_n^{(2)}(\tau )J_n({\xi }_{-})\\ F_n^{(3)}=\frac{\cos \zeta (n_0^2-n_{\perp }^2)}{ak{\sin }^2\zeta (n_{\perp }^2-n_0^2{\cos }^2\zeta )}n{H}_n^{(2)}(\tau )J_n({\xi }_{+})\\ F_n^{(4)}=j\frac{1}{\eta \sin \zeta }{H}_n^{(2)\prime }(\tau )J_n({\xi }_{-})-g_n^{(1)}{J}_n^{\prime }({\xi }_{-})\end{array},$$

Eq. (8)

$$\{\begin{array}{l}{g}_n^{(1)}=j\frac{\eta n_0^2{\chi }_{+}}{k(n_{\perp }^2-n_0^2{\cos }^2\zeta )}·H_n^{(2)}(\tau )\\ g_n^{(2)}=j\frac{n_{\perp }^2{\chi }_{-}}{\eta k(n_{\perp }^2-n_0^2{\cos }^2\zeta )}·H_n^{(2)}(\tau )\end{array},$$

Eq. (9)

$$\begin{gathered} \eta =\frac{1}{n_0}\sqrt{\frac{{\mu }_0}{{\epsilon }_0}} \\ \left\{\begin{array}{l}\tau =ak\sin \zeta \\ {\xi }_{\pm }=a·{\chi }_{\pm }\end{array}\right\{\begin{array}{l}{\chi }_{+}^2=k^2(\frac{n_{\perp }^2}{n_0^2}-{\cos }^2\zeta )\\ {\chi }_{-}^2=k^2(\frac{n_{||}^2}{n_0^2}-\frac{n_{||}^2}{n_{\perp }^2}{\cos }^2\zeta )\end{array}(n_{||}\ge n_{\perp })\{\begin{array}{l}{\chi }_{+}^2=k^2(\frac{n_{||}^2}{n_0^2}-\frac{n_{||}^2}{n_{\perp }^2}{\cos }^2\zeta )\\ {\chi }_{-}^2=k^2(\frac{n_{\perp }^2}{n_0^2}-{\cos }^2\zeta )\end{array}(n_{||}<n_{\perp }), \end{gathered}$$

where $\rho$ and $z$ represent the radial and axial distances in the cylindrical coordinates, $a$ is the radius of the cylinder, $n_0$ is the refractive index and $\eta$ is the wave impedance of the ambient medium, $k=1/\lambda$ is the wave number, $\zeta$ is the angle between the direction of incident light and the direction of cylinder, and $\mathrm{\Theta }$ is the azimuth angle of the scattered light on the scattering cone. $J_n$ is the Bessel function, ${H_n}^{(1)}=J_n+i{Y}_n$ and ${H_n}^{(2)}=J_n-i{Y}_n$, respectively, expressing the first and second kind of Hankel function.

From the above formulas, we can get Mueller matrix $\mathbf{M}(\zeta ,\mathrm{\Theta })$ and its elements as follows¹⁹

Eq. (10)

$$M(\zeta ,\mathrm{\Theta })=\frac{2}{\pi k\rho {\sin }^3\zeta }\left(\begin{array}{cccc}{m}_{11}& m_{12}& m_{13}& m_{41}\\ m_{21}& m_{22}& m_{23}& m_{42}\\ m_{31}& m_{32}& m_{33}& m_{43}\\ m_{41}& m_{42}& m_{43}& m_{44}\end{array}\right),$$

Eq. (11)

$$\{\begin{array}{l}{m}_{11}=(|J_1{|}^2+|J_2{|}^2+|J_3{|}^2+|J_4{|}^2)/2\\ m_{12}=m_{21}=({|J_1|}^2-{|J_2|}^2)/2\\ m_{13}=-m_{31}=\mathrm{Re}\{J_1{J}_4^{*}+J_2{J}_3^{*}\}\\ m_{14}=m_{41}=\mathrm{Im}\{J_1{J}_4^{*}-J_2{J}_3^{*}\}\\ m_{22}=({|J_1|}^2+{|J_2|}^2-{|J_3|}^2-{|J_4|}^2)/2\\ m_{23}=-m_{32}=\mathrm{Re}\{J_1{J}_4^{*}-J_2{J}_3^{*}\}\\ m_{24}=m_{42}=\mathrm{Im}\{J_1{J}_4^{*}+J_2{J}_3^{*}\}\\ m_{33}=\mathrm{Re}\{J_1^{*}{J}_2+J_3^{*}{J}_4\}\\ m_{34}=-m_{43}=\mathrm{Im}\{J_1{J}_2^{*}+J_3{J}_4^{*}\}\\ m_{44}=\mathrm{Re}\{J_1^{*}{J}_2-J_3^{*}{J}_4\}\end{array}.$$

2.2.

Validity for Single Scattering Process

In the present model, scattering by a birefringent cylinders with equal axial and radial refractive indices should be equivalent to a nonbirefringent cylinders studied in our previous work.¹¹ Hence, we test the validity of the new model and the new simulation program by checking if the scattering behavior of birefringent cylinders with equal axial and radial refractive indices coincides with the situation that the cylinder scatters without the birefringence.

Figure 2 shows the Mueller matrix elements as functions of the azimuth angle $\mathrm{\Theta }$ for a single scattering of a birefringent cylinder with the different axial and radial refractive indices marked by solid black line, and a birefringent cylinder with same axial and radial refractive indices marked by dotted black line. All are calculated using our newly developed program. The Mueller matrix elements of nonbirefringent cylinder calculated with our previous program are shown as the red dotted line completely overlapping with the dotted black line,¹¹ which confirms the validity of the developed program for birefringent cylinders. In Fig. 2, Mueller matrix $m_{ij}$ except for $m_{11}$ ($i$, $j=1,\dots ,4$) are all normalized by $m_{11}[m_{ij}(\mathrm{\Phi })/m_{11}(\mathrm{\Phi })]$. And $m_{11}$ is normalized by $m_{11}(\mathrm{\Phi }=0°)$, $m_{11}(\mathrm{\Phi })/m_{11}(\mathrm{\Phi }=0°)$.

Fig. 2

Single scattering Mueller matrix of a birefringent cylinder with the different axial and radial refractive indices of 1.57 and 1.53 (solid black lines), and a birefringent cylinder (dotted black lines) with the same axial and radial refractive indices of 1.57. Calculations for a nonbirefringent cylinder with the refractive index of 1.57 are shown as the red dotted lines, which overlap completely with the dotted black lines. The incident angle of light is $\zeta =45\mathrm{deg}$ and the diameter of the cylinder is 1.5 μm. The abscissa is azimuth 0–180 deg.

Figure 2 shows that the Mueller matrix elements for single scattering by a birefringent cylinder are clearly different from those by a nonbirefringent cylinder. To examine the polarization features of fibrous tissue microstructures, contributions by the intrinsic birefringence on the cylinders need to be considered.

3.

Samples and Experiments

For further verification of our Monte Carlo program for birefringent cylinders and demonstration of the influence of birefringent cylinders on polarization measurements, we present the simulation and experimental results of forward scattering Mueller matrix of nonbirefringent and birefringent cylinders, respectively. In the experiments, we prepare two samples by winding either well-aligned nonbirefringent glass fibers or birefringent silk fibers around small frames.²⁰ The experimental setup has been used in a previous study for forward scattering Mueller matrix measurements.⁵ Before the experiments, we measured air and a polarizer as the standard samples, and the calibration errors can be estimated at about 5%.¹³

Figure 3(a) represents the experimental result for the silk sample, and Fig. 3(c) is for the glass fiber sample. Earlier studies have shown that the silk fibers contain fibrous substructures of 1.5-μm diameter,⁵ and the refractive indices along radial and axial directions are 1.53 and 1.57.²⁰ The glass fibers have a refractive index of 1.547 and about 10 μm in diameter (observed by SEM).^13,20 The corresponding simulation results of birefringent and nonbirefringent cylinders are shown in Figs. 3(b) and 3(d) with the scattering coefficient ${\mu }_s=30\mathrm{c}{\mathrm{m}}^{-1}$ and the other parameters the same as the experimental conditions; the wavelength in the simulation is 633 nm. Figures 3(a)–3(d) show very good agreement between experiments and simulations and confirm the validity of the scattering model and the corresponding Monte Carlo simulation program for scattering medium containing orderly arranged birefringent cylinders.

Fig. 3

Experimental and simulation results of forward scattering Mueller matrix. (a) and (b) are experimental and simulated results of birefringent cylinders (well-aligned silk fiber). (c) and (d) are experimental and simulated results of nonbirefringent cylinders (well-aligned glass fiber).

It can be seen from Fig. 3, for both the well-aligned nonbirefringent or birefringent cylinders, $m_{12}$ and $m_{21}$ are nearly equal, and $m_{22}$ is significantly bigger than the other polarization dependent elements. For the Mueller matrix elements $m_{33}$ and $m_{44}$, however, they are different for the birefringent and nonbirefringent cylinders, probably due to the intrinsic birefringence on cylindrical scatters. Considering the possible small fluctuation of arrangement and size for cylinders during the preparation of experimental phantoms, the slight differences between experiments and simulations are possibly due to the uncertainty of the simulation approximation. Further simulations with different parameters in the model show that the effects on the Mueller matrix elements due to the birefringence in the cylindrical scatterers are complicated and sensitive to the diameter and refractive indices of the cylinders.

4.

Results and Discussion

The Mueller matrix can characterize polarization properties of a sample. There are several other methods to decompose Mueller matrix; in a previous work,¹³ we used the Lu-Chipman Mueller matrix polar decomposition method to “decompose” a Mueller matrix $M$ into three constituent “basic” matrices, representing depolarization ($M_{\mathrm{\Delta }}$), retardance ($M_R$), and diattenuation ($M_D$).²¹ From Eq. (10), we can see that the value of Mueller matrix elements depends on various factors, but we focus on the relationship between Mueller matrix and intrinsic birefringence of cylindrical scatterers. By Monte Carlo simulations, we examine how the birefringent cylinders affect the retardance ($\delta$) for an anisotropic sample. We simulate Mueller matrix images using different birefringence values of the cylinders, evaluate the retardance ($\delta$), depolarization ($\mathrm{\Delta }$), and diattenuation ($D$) pixel by pixel using Mueller matrix polar decomposition, and then calculate their mean values.

Figure 4 shows the retardance ($\delta$) from the simulated data as functions of the radial and axial refractive indices difference $\mathrm{\Delta }{n}_c$ and diameters of the cylinders. The diameters of the cylinders are from 0.6 to 2 μm, and the values of birefringence $\mathrm{\Delta }{n}_c$ are from 0 to $9\times {10}^{-3}$; the thickness of the medium is 1 cm. Figures 4(a)–4(c), respectively, represent three cases of different scattering coefficients, ${\mu }_s=10$, ${\mu }_s=20$, ${\mu }_s=30\mathrm{c}{\mathrm{m}}^{-1}$. It can be seen that the contributions by the birefringent cylinders to the retardance strongly depend on their diameter and density. For $\mathrm{\Delta }{n}_c=0$, the retardance decreases as the diameter of the cylinders increases and increases as the scattering coefficient increases, as found in the previous work.¹³ There is a good positive linear relationship between $\delta$ and the intrinsic birefringence of the cylinders, which appear to be very similar to the effects due to the birefringent medium between scatterers.¹³ The slope of the retardance-birefringence curves is shown in Fig. 4, which represents how sensitively the retardance varies with the birefringence on the cylinder, and also increases with the diameter of the birefringent cylinders and the scattering coefficient. We extract the values of slopes of the curves by linear function fitting. We plot in Figs. 5(a) and 5(c) these slopes as functions of the diameter and scattering coefficient of the birefringent cylinders. After a normalized transformation of abscissa axis by $d·{\mu }_s$, Figs. 5(a) and 5(c) change into Figs. 5(b) and 5(d); there is a clear positive and approximately linear correlation between the slope and the product of cylinder diameter and density.

Fig. 4

The relationship between the birefringence of the cylinders and the retardance $\delta$. Different curves represent different diameters of birefringent cylinders, (a)–(c) represent different scattering coefficients of the sample.

Fig. 5

Slopes extracted from Figs 4(a)–4(c), (a) and (c) represent the relationship between the diameter and scattering coefficient ${\mu }_s$ and $\delta$, (b) and (d) represent the relationship between the product of diameter $d$ and scattering coefficient ${\mu }_s$ and $\delta$.

In addition, we did some investigation about depolarization ($\mathrm{\Delta }$) and diattenuation ($D$) in the birefringent cylinders model. We can find that, compared with the retardance, $\mathrm{\Delta }$ and $D$ just have a small fluctuation with the birefringence on the birefringent cylinders.

In the previous study,¹³ we have analyzed the contributions of retardance by the nonbirefringent cylinders with the birefringent medium between scatterers. According to our previous research work, Fig. 6 shows $\delta$ as functions of the medium birefringence in the cylinder-birefringence sample. The retardance $\delta$ increases linearly with the birefringence of the ambient medium $\mathrm{\Delta }{n}_m$ and decreases with the diameter of the cylinders. Here, we compare two types of birefringence, respectively, in the birefringent cylinders and the cylinder-birefringence sample. For the birefringent-cylinder scattering samples (Figs. 4 and 5), the diameter of the cylinders plays an important role on the retardance. However, for the cylinder-birefringence model, the birefringent medium surrounding the nonbirefringent cylinder is the major source of the total retardance.

Fig. 6

Change in the birefringence of the medium surrounding the nonbirefringent cylinder. The diameters of the cylinder were 1.0, 1.4, and 2.0 μm, and the scatter coefficient is ${\mu }_s=20\mathrm{c}{\mathrm{m}}^{-1}$; the birefringence of the ambient medium $\mathrm{\Delta }{n}_m$ varies from 0 to $5\times {10}^{-6}$.

In Fig. 6, the well-aligned nonbirefringent cylinders are in parallel to the birefringence axis of the birefringent medium, where the $\mathrm{\Delta }{n}_m$ range is set to a smaller scale to ensure $\delta$ stays in the range of 0 to $\pi$. The scattering coefficient of cylinders affects only the initial value of the retardance, but not the slope of the curve.

For the birefringent medium, the retardance can be written as $\delta =2\pi n_{m}s\mathrm{\Delta }{n}_m/\lambda$, where $s$ is the projection along the $x$ direction of the moving distance between two successive scattering events, $\mathrm{\Delta }{n}_m$ is the birefringence of the medium, $n_m$ is the average refractive index, and $\lambda$ is the wavelength of light. If we set the same birefringence parameter as Fig. 6, the coefficient $2\pi n_{m}s/\lambda$ is nearly equal to the slope of curves in Fig. 6, which is around ${10}^4–{10}^5$. By contrast, in the birefringent cylinders samples, the slope of the curves shown in Fig. 5 is only on the order of magnitude of ${10}^2$, indicating that the intrinsic birefringence of cylindrical scatterers has weaker influence on the total retardance than the birefringence of the ambient medium.

5.

Conclusion

We established an anisotropic scattering model including birefringent cylinders with different refractive indices along the radial and axial directions and developed the corresponding Monte Carlo simulation program. The validity of simulation program is tested by experiments using well-aligned nonbirefringent and birefringent cylinders. The present study shows that the retardance due to scattering of well-ordered birefringent cylinders comes from two sources: scattering by nonbirefringent cylinders and an additional part associated with the intrinsic birefringence on cylinders. Moreover, the retardance increases linearly with $\mathrm{\Delta }{n}_c$ and almost linearly with the diameter of the birefringent cylinder. We also compare the influence of the birefringent cylinders on the retardance with our previous cylinder-birefringence model. For the same birefringence, the birefringent cylinders usually contribute much less to the retardance $\delta$ than the birefringent medium.