1.

Introduction

The last decade has seen growing interest to the propagation of coherent polarized light in turbid tissue-like scattering media.^1,2 The majority of recent studies associated with the polarized light are focused on the development of new diagnostic modalities for noninvasive characterization of biological tissues with the special attention to cancer screening.^3,4 The incident polarized light is multiply scattered along its propagation within the biological tissue and becomes depolarized. The depolarization ratio (DR) depends strongly on the size and shape of scattering particles,^5,6 as well as on the number of scattering events,^7,8 and is independent of the state of polarization of incident light.⁹ Thus, scattering properties of a turbid medium can be evaluated quantitatively by monitoring the evolution of the polarization state of scattered light in comparison with the polarization of incident light. The Mueller matrix approach is typically utilized to assess the properties of the medium based on polarization measurements.² Potentially, this approach can be used for noninvasive cancer diagnosis, e.g., for colon cancer detection.¹⁰ In fact, typically utilized in diagnostic practice, the backscattered polarized light contains not only the light scattered within the medium but also the light reflected by the surface of the medium. Presently, there is an escalating interest in the surface roughness that potentially can be used in cancer diagnosis, e.g., for an assessment of grade of skin neoplasia.¹¹ To introduce a roughness score, the malignant features have been observed with the reflectance confocal microscopy.¹² Relief patterns of the skin surface of benign and malignant lesions of the skin have been studied by microtopography.¹³ The polarization imaging with high angles of incidence has been extensively used in routine clinical studies of skin roughness.¹⁴ Tchvialeva et al.¹⁵ suggested a methodology of quantifying skin surface roughness by analyzing the laser speckle contrast of the backscattered linearly polarized light. In spite of numerous studies in the field, a clear understanding of influence of surface roughness on the formation of polarized laser speckles is still required.

In this paper, in the framework of further development of a unified computational model for the needs of biomedical optics,^16,17 we present the electric field Monte Carlo (MC) approach specially developed for simulation of backscattering of coherent polarized light from turbid tissue-like scattering media with rough surfaces. The developed computational model is based on the Jones formalism and takes into account the wave properties of light, including temporal coherence, polarization, phase change on the reflection and/or refraction at the rough medium boundary, and interference. Using the developed electric field MC approach, we explore the influence of the surface roughness on the resulted depolarization of coherent linearly polarized light backscattered from various turbid media of known optical properties, which are in order of human skin lesions. In particular, we examine the laser speckle patterns formation and dependence of DR on the surface roughness of a phantom medium.

2.

Theory

2.1.

Basic Elements of Light Wave Propagation in Turbid Medium

The concept of light as a wave is fundamental for the phenomena of coherence and polarization.¹⁸ Coherence is a quantitative measure of the degree of phase correlation of the light wave, whereas the polarization is defined by the electric field vector $\overrightarrow{\mathbf{E}}$ lying in the plane perpendicular to the direction of the light wave propagation. The light wave with the electric field $\overrightarrow{\mathbf{E}}$ is incident normally to the surface of the medium and is propagating along the $z$-direction in frame of the standard Cartesian coordinate system. The complex field components with amplitudes $E_{0_x}$, $E_{0_y}$ and phases ${\varphi }_x$, ${\varphi }_y$ at time $t$ are presented in terms of the two-element matrix, known as a “Jones vector”¹⁹

Eq. (1)

$$\overrightarrow{\mathbf{E}}=\left[\begin{array}{c}{E}_x\\ E_y\end{array}\right]=\left[\begin{array}{c}{E}_{0_x}{e}^{i(kz-\omega t+{\varphi }_x)}\\ E_{0_y}{e}^{i(kz-\omega t+{\varphi }_y)}\end{array}\right],$$

where $i$ is the imaginary unit, $k$ is the wave number ($k=2\pi /\lambda$, $\lambda$ is the wavelength of incident light), $\omega$ is the angular frequency, and $\mathrm{\Delta }\varphi ={\varphi }_x-{\varphi }_y=0$ (i.e., linearly polarized light is considered). Absorption of light in the medium changes the amplitude of the components of the field presented in Eq. (1). The complex refractive index $\tilde{n}=n_{\mathrm{R}}+{\mathrm{in}}_{\mathrm{I}}$ is used to define the complex phase shift induced by the medium

Eq. (2)

$${\overrightarrow{\mathbf{E}}}_{\text{out}}={\overrightarrow{\mathbf{E}}}_{\text{in}}{e}^{i\frac{2\pi }{\lambda }{n}_{\mathrm{R}}S}{e}^{-\frac{2\pi }{\lambda }{n}_{\mathrm{I}}S},$$

where $S$ is the trajectory of a light wave in the medium, defined as a sum of pathlengths $l_i$ between successive scattering events, counted from the point of entrance into the medium to the point of output ($S=\sum l_i$), $\mathrm{Re}(\tilde{n})=n_{\mathrm{R}}$ is defined by the ordinary refractive index, and the term $e^{i(2\pi /\lambda )n_{\mathrm{R}}S}$ describes the wave retarded by $(2\pi /\lambda )n_{\mathrm{R}}S$; $\mathrm{Im}(\tilde{n})=n_{\mathrm{I}}$ is known as the extinction coefficient and determines the rate of absorption of light in the medium. The term $e^{-(2\pi /\lambda )n_{\mathrm{I}}S}$ defines the amplitude attenuation of the Jones vector components [see Eq. (1)] and is related to the well-known Beer–Lambert–Bouguer law

Eq. (3)

$$I=I_0{e}^{-{\mu }_{t}S},$$

where ${\mu }_t$ is the attenuation coefficient, ${\mu }_t={\mu }_{\mathrm{s}}+{\mu }_{\mathrm{a}}$, ${\mu }_{\mathrm{a}}$ is the absorption coefficient $({\mu }_{\mathrm{a}}=4\pi n_{I}S/\lambda )$, and ${\mu }_{\mathrm{s}}$ is the scattering coefficient.

The superposition of $N_{\mathrm{ph}}$ light waves at the detecting area for the absorption-free (${\mu }_{\mathrm{a}}=0$) medium is defined as

Eq. (4)

$$\begin{gathered} E_x^2=\sum _{m=1}^{N_{\mathrm{ph}}}{E}_{0_{xm}}^2+2\sum _{m>v}^{N_{\mathrm{ph}}}\sum _{v=1}^{N_{\mathrm{ph}}}{E}_{0_{xm}}{E}_{0_{xv}}\cos ({\alpha }_v-{\alpha }_m)\exp [-{\left(\frac{\mathrm{\Delta }{L}_{m,v}}{l_{\mathrm{c}}/\sqrt{2}}\right)}^2], \\ {E}_y^2=\sum _{m=1}^{N_{\mathrm{ph}}}{E}_{0_{ym}}^2+2\sum _{m>v}^{N_{\mathrm{ph}}}\sum _{v=1}^{N_{\mathrm{ph}}}{E}_{0_{ym}}{E}_{0_{yv}}\cos ({\alpha }_v-{\alpha }_m)\exp [-{\left(\frac{\mathrm{\Delta }{L}_{m,v}}{l_{\mathrm{c}}/\sqrt{2}}\right)}^2], \end{gathered}$$

where ${\alpha }_{m,v}=k{S}_{m,v}+{\phi }_{m,v}$ shows the phases for $x$ and $y$ components, respectively, $\mathrm{\Delta }{L}_{m,v}$ is the pathlength difference ($\mathrm{\Delta }{L}_{m,v}=S_m-S_v$), $l_{\mathrm{c}}$ is the temporal coherence length of the incident light, and $\phi$ is the phase shift of individual waves indicated by indices $m$ and $v$.

Reflection and refraction of light waves at the medium boundary are described by the Fresnel’s reflection coefficients and defined for the transverse magnetic (TM) and transverse electric (TE) components (or, respectively, for $x$ and $y$ components) as¹⁹

Eq. (5)

$$\begin{gathered} R_{\mathrm{TM}}={\left|\frac{n_1\cos {\theta }_{\mathrm{t}}-n_2\cos {\theta }_{\mathrm{i}}}{n_1\cos {\theta }_{\mathrm{t}}+n_2\cos {\theta }_{\mathrm{i}}}\right|}^2, \\ {R}_{\mathrm{TE}}={\left|\frac{n_1\cos {\theta }_{\mathrm{i}}-n_2\cos {\theta }_{\mathrm{t}}}{n_1\cos {\theta }_{\mathrm{i}}+n_2\cos {\theta }_{\mathrm{t}}}\right|}^2, \end{gathered}$$

where $n_1$ and $n_2$ are the refractive indices of the external and internal media, respectively, ${\theta }_{\mathrm{i}}$ is the angles of light incidence, and ${\theta }_{\mathrm{t}}$ is the angle refraction on the medium boundary. Figure 1 shows the Fresnel’s reflection coefficients counted by Eq. (5) for external $(n_1<n_2)$ and internal $(n_1>n_2)$ reflection.

Fig. 1

Fresnel’s reflection coefficients $R_{\mathrm{TM}}$ and $R_{\mathrm{TE}}$ for the external reflection, when $n_1=1.0$, $n_2=1.5$ (solid line), and for the internal reflection, when $n_1=1.5$, $n_2=1.0$ (dashed line), plotted as the functions of incidence angle ${\theta }_{\mathrm{i}}$.

In the case of external reflection ($n_1<n_2$), a $\pi$-phase shift occurs at any angle of incidence for the TE field component and for the TM field component when ${\theta }_{\mathrm{i}}<{\theta }_{\mathrm{p}}$. Here, ${\theta }_{\mathrm{p}}$ is the polarizing angle at which $R_{\mathrm{TM}}=0$, also known as Brewster’s angle. For the internal reflection ($n_1>n_2$) the phase shift is defined as¹⁹

Eq. (6)

$${\phi }_{\mathrm{TM}}=\{\begin{array}{ll}0,& {\theta }_{\mathrm{i}}<{\theta }_{\mathrm{p}},\\ \pi ,& {\theta }_{\mathrm{p}}<{\theta }_{\mathrm{i}}<{\theta }_{\mathrm{c}},\\ -2\arctan \left(\frac{\sqrt{{\sin }^2{\theta }_{\mathrm{i}}}-n^2}{n^2\cos {\theta }_{\mathrm{i}}}\right)+\pi ,& {\theta }_{\mathrm{i}}>{\theta }_{\mathrm{c}},\end{array}$$

Eq. (7)

$${\phi }_{\mathrm{TE}}=\{\begin{array}{ll}0,& {\theta }_{\mathrm{i}}<{\theta }_{\mathrm{c}},\\ -2\arctan \left(\frac{\sqrt{{\sin }^2{\theta }_{\mathrm{i}}}-n^2}{\cos {\theta }_{\mathrm{i}}}\right),& {\theta }_{\mathrm{i}}>{\theta }_{\mathrm{c}},\end{array}$$

where ${\theta }_{\mathrm{c}}$ is the critical angle $[{\theta }_{\mathrm{c}}=\arcsin (n_2/n_1)]$ when total internal reflection occurs. The phase shifts $\phi$ on the reflection for TE and TM modes depending on the angle of incidence are shown in Fig. 2.

Fig. 2

Phase shifts on the reflection ${\phi }_{\mathrm{TE}}$, ${\phi }_{\mathrm{TM}}$ depending on the angle of light incidence ${\theta }_{\mathrm{i}}$: (a) external and (b) internal TE reflection modes; (c) and (d) are, respectively, external and internal TM reflection modes, counted for the refractive indices: $n_1=1.0$, $n_2=1.5$.

Thus, considering superposition of light waves at the detector, the relative phase shift between $m$ and $v$ light waves due to their pathlength difference within the medium $(k{S}_{m,v})$, the phase shift due to reflection/refraction on the medium boundary ${({\phi }_{\mathrm{TE}},{\phi }_{\mathrm{TM}})}_{m,v}$, as well as the mutual phase shift due to variations of roughness heights from point to point

Eq. (8)

$$\mathrm{\Delta }\phi {(x,y)}_{m,v}=\left(\frac{4\pi }{\lambda }\right)\mathrm{\Delta }{h}_{\mathrm{rms}}{(x,y)}_{m,v}\sin {\theta }_i,$$

should be taken into account. Here, $\mathrm{\Delta }{h}_{\mathrm{rms}}=h_m-h_v$ is the difference between $m$ and $\nu$ root–mean-square heights, as shown in Fig. 3.

Fig. 3

Schematic presentation of randomly generated rough surface outlining and its key features, such as $h_{\mathrm{rms}}$ and correlation length $L_{\mathrm{corr}}$, the parameters used for quantitative characterization of rough surfaces.²⁰

In the framework of this study, with the final aim to simulate the speckle patterns formation and to understand the role of surface roughness in the depolarization of linearly polarized light backscattered from turbid tissue-like scattering media, both the mutual phase shifts due to the pathlength difference within the medium and reflection/refraction on the rough surface of the medium were implemented into the electric field MC model.

3.

Materials and Methods

Figure 4 schematically presents the experimental setup used for the MC model validation. A phantom with a rough surface is illuminated normally by a laser diode (5 mW, $\lambda =663\mathrm{nm}$, $l_{\mathrm{c}}=30\mu \mathrm{m}$, Flex, B&WTek Inc.) with a $200\mu \mathrm{m}$-diameter of the beam. The laser speckle patterns are observed by a CCD camera (Matrix Vision GmbH, mvBlueFOX-M124G—8 bit and mvBlueFOX-M224G—8 bit) utilizing a closed-circuit television (CCTV) lens (50 mm, $F/1.4$, Pentax Inc.). The polarizers (Thorlabs Inc., WP25L-UB) are used to control the input–output states of polarization. The Thorlabs pellicle beam splitter equally splits the co- and cross-polarized components without retardation. The processing of detected co- and cross-polarized images includes dark signal subtraction and accurate pixel-to-pixel correspondence adjustment. Spatial distribution of the DR is calculated by Eq. (16) utilizing spatial intensity distributions of co- and cross-polarized [$I_{\parallel }(x,y)$ and $I_{\perp }(x,y)$] speckle patterns.

Fig. 4

Schematic presentation of the laser speckle imaging system used in the experimental studies.

The phantoms used in the experimental studies have been prepared by utilizing silicone resin (the refraction index $n_{\mathrm{R}}\approx 1.41$) with the embedded silicone pigments.⁴⁶ Optical properties of the phantoms were assessed by the standard integrated sphere-based approach,⁴⁷ utilizing a red diode laser (B&WTek Flex Inc.) presented in Table 1, where ${\mu }_{\mathrm{s}}^{\prime }$ is the reduced scattering coefficient,⁴⁰ defined as: ${\mu }_{\mathrm{s}}^{\prime }={\mu }_{\mathrm{s}}(1-g)$. The surface of phantoms used in the study has been replicated using metal standard (Microshurf #334, Rubert Co. Ltd.) with $h_{\mathrm{rms}}$ in the range from 10 to $66\mu \mathrm{m}$ and $L_{\mathrm{corr}}=100\mu \mathrm{m}$. This range of surface roughness parameters corresponds to the typical roughness of human skin.⁴⁷ The actual parameters of surface roughness of phantoms have been confirmed by measurements with the WYKO NT1100 optical profilometer (Veeco) providing vertical resolution $0.05\mu \mathrm{m}$.

Table 1

Optical properties of phantoms used in the experiment.

The standard algorithm⁴⁸ was adopted to generate profiles of the surfaces of the modeling media with the known roughness. In this algorithm, an uncorrelated Gaussian distribution of random numbers is generated. When convolved with the Gaussian filter, a series of points representing spatial distribution of points along the surface corresponding to the randomly selected values of $h_{\mathrm{rms}}$ are produced. Then generated values of $h_{\mathrm{rms}}$ are implemented into the surface as a mesh (see Fig. 3). Thus, the surface roughness profiles of modeling media presented in Fig. 5 were generated in accordance with the surface profiles of phantoms used in the experiment.

Fig. 5

The surface roughness profiles generated by the algorithm adopted from Tsang et al.,⁴⁸ from the top: $h_{\mathrm{rms}}$ are 2.5, 34.4, and $65.8\mu \mathrm{m}$, and $L_{\mathrm{corr}}=100\mu \mathrm{m}$ for all profiles.

4.

Results and Discussion

Spatial intensity distributions of co- and cross-polarized light [$I_{\parallel }(x,y)$ and $I_{\perp }(x,y)$] and the DR counted with the electric field MC approach described above are presented in Fig. 6. The optical properties of the medium are summarized in Table 1. The parameters of the rough surface are varied $h_{\mathrm{rms}}=0$, 2.5, 34.4, $65.8\mu \mathrm{m}$, and $L_{\mathrm{corr}}=100\mu \mathrm{m}$.

Fig. 6

Spatial intensity distributions of co- and cross- polarized light $I_{\parallel ,\perp }(x,y)$ and depolarization ratio $\mathrm{DR}(x,y)$ counted for the modeling media with optical properties corresponding to the phantom I (see Table 1). The roughness from top row to bottom correspond, respectively, to $h_{\mathrm{rms}}=0$, 2.5, 34.4, and $65.8\mu \mathrm{m}$, and $L_{\mathrm{corr}}=100\mu \mathrm{m}$.

The intensity of backreflected copolarized light [$I_{\parallel }(x,y)$] in the area of illumination is dropped with the increase of $h_{\mathrm{rms}}$, whereas the contrast of DR is increased. The reduction of $I_{\parallel }(x,y)$ is explained by the growing variability of incidence angle ${\theta }_{\mathrm{i}}$ with an increase of $h_{\mathrm{rms}}$. As a result in the far-zone within the fixed narrow detection angle ($<7\mathrm{deg}$), the Fresnel’s reflection coefficients (see Fig. 1) and the phase shifts (due to reflection, see Fig. 2) become highly variable. Therefore, the intensities of co- and cross-polarized light $I_{\parallel ,\perp }(x,y)$ are dropped as prescribed by Eqs. (14) and (15). Reflection from the medium surface is not contributing to the speckles’ formation outside the area of illumination, therefore, the intensity distribution of $I_{\parallel ,\perp }(x,y)$ is formed there only by the photon packets multiply scattered within the medium and experiencing internal reflection and refraction.

To illustrate the influence of surface roughness on the detected signal formation, Fig. 7 shows the spatial distributions of the detector depth sensitivity⁴⁹ for co- and cross-polarized components ($I_{\parallel }$ and $I_{\perp }$) and their difference, counted as an analogy of the DR [Eq. (16)]. The effective optical pathlengths $S$ within the scattering medium counted by MC defines the paths that photon packets traveled within the medium from the point of incidence to the point of exit. Thereby, spatial distribution of the effective optical pathlengths within the medium can be considered in terms of the temporal point-spread function and this defines the detection sensitivity along the depth,⁴⁹ known also as a sampling volume.⁵⁰

Fig. 7

Spatial intensity distributions of the detecting sensitivity along the depth for $I_{\parallel ,\perp }(x,y)$ and their difference counted as DR for the media of the same properties, as presented in Fig. 6.

The results presented in Fig. 7 clearly show the influence of surface roughness on the detector depth sensitivity profile for co- and cross-polarized light ($I_{\parallel }$ and $I_{\perp }$), as well as their difference, defined by Eq. (16). The detector depth sensitivity is significantly reduced with the increase of $h_{\mathrm{rms}}$ (see Fig. 7), whereas the speckle patterns’ contrast for $I_{\parallel }(x,y)$ and $I_{\perp }(x,y)$ with the higher values of $h_{\mathrm{rms}}$ becomes less comparable with the smoother surfaces (see Fig. 6).

In a similar manner, Figs. 8 and 9 show, respectively, the spatial intensity distributions and the detecting depth sensitivities for co-/cross-polarized light ($I_{\parallel }$, $I_{\perp }$) and DR counted for modeling media that are similar to the phantom II, with the same values of $h_{\mathrm{rms}}$ and $L_{\mathrm{corr}}$ (i.e., $h_{\mathrm{rms}}=0$, 2.5, 34.4, $65.8\mu \mathrm{m}$ and $L_{\mathrm{corr}}=100\mu \mathrm{m}$), and the optical properties presented in Table 1. As one can see with the same roughness spatial intensity distributions and detector depth sensitivity for $I_{\parallel ,\perp }(x,y)$ and DRs are mainly influenced by optical properties of the scattering medium. Again, the bright circular spots seen at the copolarized $I_{\parallel }$ intensity distribution (see Figs. 8 and/or 9) correspond to the contribution of both specular reflection and diffuse reflection from the medium, although, for the cross-polarized light ($I_{\perp }$), the speckle patterns are formed mainly by the diffuse photons experienced internal reflection and refraction on the surface of the medium.

Fig. 8

Spatial intensity distributions of co- and cross-polarized light $I_{\parallel ,\perp }(x,y)$ and depolarization ratio $\mathrm{DR}(x,y)$ counted for the modeling media with optical properties corresponding to the phantom II (see Table 1). The roughness from top row to bottom correspond, respectively, to $h_{\mathrm{rms}}=0$, 2.5, 34.4, and $65.8\mu \mathrm{m}$, and $L_{\mathrm{corr}}=100\mu \mathrm{m}$.

Fig. 9

Spatial intensity distributions of the detecting sensitivity along the depth for $I_{\parallel ,\perp }(x,y)$ and their difference counted as DR for the media with the same properties, as presented in Fig. 8.

Bearing in mind that the results presented in Figs. 7 and 9 cannot be measured experimentally, to quantitatively compare the results of MC modeling and the results of experimental measurements of the spatial intensity distribution, we treated the speckle patterns utilizing second-order statistics.⁴⁶ In this approach, the spatial distributions of intensities $I_{\parallel ,\perp }(x,y)$ are averaged along the selected radial distance $d$ around the area of illumination as shown in Fig. 10.

Fig. 10

The results of experimental measurements of $I_{\parallel }(x,y)$, raw data (left) and schematic presentation of the procedure of averaging along distances $d$ around the area of illumination (right). The dimensions of the presented image are similar to those used in the simulation.

Thus, by collecting the average intensities of co- and cross-polarized light ($⟨I_{\parallel }(d)⟩$ and $lang;I_{\perp }(d)⟩$) the DR is counted as

Radial distribution of the averaged DR [Eq. (17)] in terms of optical thickness $d/l*$ is presented in Fig. 11 in comparison with the results of experimental measurements for the phantoms I and II (presented in Table 1); $l*$ is the transport length ($l^{*}=1/{\mu }_{\mathrm{s}}^{\prime }$) defining the length over which the direction of photon packet’s propagation is randomized. The parameters of roughness used in the simulation for both phantoms are $h_{\mathrm{rms}}=34.4\mu \mathrm{m}$ and $L_{\mathrm{corr}}=100\mu \mathrm{m}$.

Fig. 11

Radial distribution of DR measured experimentally for the phantom I (□) and II (▿) in comparison with the results obtained by the electric field MC (▪ and ▾). The solid and dashed lines represent the linear fits. The standard deviation for the results of experimental studies is in a range of 5–7%.

As one can see, the results of the developed computational model agree well with the results of the experiment (see Fig. 11). Similar agreements were found between the results of modeling and the experiments for the phantoms with other roughness parameters, these are not presented here for brevity.

Figure 12 shows relative slopes of DR [Eq. (17)], averaged along radial distance in a similar manner as has been done in the experiment (see Fig. 10) and plotted as the result of normalized linear fit for the samples with different values of $h_{\mathrm{rms}}$. The relative DR is an alternative representation of the results presented in Fig. 11 and apparently can be used as a quantitative measure of the relative changes of surface roughnesses.

Fig. 12

Radial distribution of normalized DR plotted for the phantom I with $h_{\mathrm{rms}}=0$ (▪), $h_{\mathrm{rms}}=2.5\mu \mathrm{m}$ (▴), $h_{\mathrm{rms}}=34.4\mu \mathrm{m}$ (•), and $h_{\mathrm{rms}}=65.8\mu \mathrm{m}$ (♦).

5.

Summary and Conclusions

Finally, with the further development of the unified computational model for the needs of biomedical optics and biophotonics,^16,17 we present the extension of the electric field MC for simulation of backreflectance and backscattering of coherent polarized light from turbid tissue-like scattering media with rough surfaces. The mutual phase shifts due to the pathlength difference within the medium and due to the reflection/refraction on the rough surface of the medium are taken into account. The validation of the model includes the creation the phantoms of various roughness and optical properties, measurements of co- and cross-polarized components of the backscattered/reflected light, and its analysis and extensive computer modeling accelerated by the parallel computing on the NVIDIA GPUs using CUDA. To characterize the speckle patterns of the backscattered light, the spatial distribution of co- and cross-polarized light were measured/simulated and DR was analyzed taking into account radial symmetry. The DR was analyzed versus surface roughness for the phantoms for which the optical properties are of the same order as human skin lesions. The results of the computer modeling agree reasonably well with the results of experimental measurements. The developed computational approach allows straight modeling of speckle patterns’ formation, their texture analysis, and, arguably, a quantitative assessment of surface roughness from the spatial distribution of the intensity of co- and cross-polarized scattering light. The analysis of the spatial intensity distribution is based on the second-order statistics that shows a strong correlation with the surface roughness, both with the results of modeling and experiment. The developed MC tool can be used for the direct simulation of light scattered by the tissue-like phantoms of various surface roughness. Further development of the model might include the periodic, anisotropic, and chiral surfaces as well as birefringent properties of the medium.