Research Papers: General

Open source software for electric field Monte Carlo simulation of coherent backscattering in biological media containing birefringence

[+] Author Affiliations
Andrew J. Radosevich, Jeremy D. Rogers, İlker R. Çapoğlu, Nikhil N. Mutyal, Prabhakar Pradhan, Vadim Backman

Northwestern University, Biomedical Engineering Department, Tech E310, 2145 Sheridan Road, Evanston, Illinois 60208

J. Biomed. Opt. 17(11), 115001 (Nov 02, 2012). doi:10.1117/1.JBO.17.11.115001
History: Received July 13, 2012; Revised October 1, 2012; Accepted October 2, 2012
Text Size: A A A

Open Access Open Access

Abstract.  We present an open source electric field tracking Monte Carlo program to model backscattering in biological media containing birefringence, with computation of the coherent backscattering phenomenon as an example. These simulations enable the modeling of tissue scattering as a statistically homogeneous continuous random media under the Whittle-Matérn model, which includes the Henyey-Greenstein phase function as a special case, or as a composition of discrete spherical scatterers under Mie theory. The calculation of the amplitude scattering matrix for the above two cases as well as the implementation of birefringence using the Jones N-matrix formalism is presented. For ease of operator use and data processing, our simulation incorporates a graphical user interface written in MATLAB to interact with the underlying C code. Additionally, an increase in computational speed is achieved through implementation of message passing interface and the semi-analytical approach. Finally, we provide demonstrations of the results of our simulation for purely scattering media and scattering media containing linear birefringence.

Figures in this Article

Diffuse reflectance measurements provide a method to noninvasively characterize the physical composition of biological tissue with known sensitivities to the concentration of chromophores (e.g., hemoglobin and melanin) as well as scattering structures as small as tens of nanometers in size.1 Typically, models of diffuse reflectance neglect the presence of birefringent materials due to the assumption that their contribution to the measured signal should be small. Yet, biological tissue contains a large number of structures which exhibit either linear birefringence due to structural alignment (e.g., lipid bilayers, collagen fibers, and muscle fibers) or circular birefringence, also known as optical activity, due to the presence of chiral molecules (e.g., glucose and certain amino acids). Because of the common presence of such substances in biological tissue, it is plausible that in general their effects should not be neglected. Wang et al. were the first to model the effect of linear birefringence on the shape of the spatial reflectance profile using Monte Carlo simulation.2 Their results demonstrate alterations occuring at subdiffusion length-scales (i.e., source-detector separations less than a transport mean free path ls*) due to the presence of linear birefringence. One easy-to-implement experimental technique that enables the measurement of such changes is coherent backscattering (CBS).3

CBS, also known as enhanced backscattering (EBS), is a coherence phenomenon in which rays traveling time-reversed paths constructively interfere to form an angular intensity peak centered in the backscattering direction.47 The shape of the CBS peak is related to the spatial reflectance profile through Fourier transformation.3,8 Because of this relationship, CBS is sensitive to the optical scattering, absorption, and polarization properties that alter the spatial reflectance profile. Employing these sensitivities, CBS has been used to study such objects as fractal aggregates,9,10 amplifying random media,11 cold atoms,12 liquid crystals,13,14 and biological tissue1517 to name a few.

For use in biomedical applications, CBS offers a noninvasive tool to interrogate the optical properties of biological materials at subdiffusion length-scales where information about the scattering phase function is preserved.17 As such, CBS can be used to quantify the absorption coefficient μa, the scattering coefficient μs, the anisotropy factor g, as well as a second shape parameter of the phase function D using a single spectral measurement.17,18 Combining these sensitivities with the ability to selectively interrogate different layers of tissue through implementation of a partial spatial coherence source, CBS has become a promising technique for the characterization and detection of colorectal and pancreatic cancers.19,20

In order to accurately characterize a tissue sample using CBS, it is necessary to understand the dependencies of the peak shape on tissue structural composition. While various analytical formalisms have been developed to describe the CBS peak shape for different sample properties, each of these calculations rely on simplifying assumptions (e.g., scalar approximation21 or double scattering22) which cannot fully describe the complex sensitivities of the CBS peak. As a more rigorous but time-consuming alternative, polarized light Monte Carlo simulations give results that are in experimental agreement provided that a sufficient number of photon realizations are computed and the underlying model accurately describes the medium under observation.3,8,17 Numerous Monte Carlo codes have been developed to simulate CBS in a wide variety of different materials. Of particular importance for this paper are CBS codes that have been developed to model and calculate tissue optical properties,18,23,24 electric field tracking codes,25,26 and codes that implement the semi-analytical approach (also known as the partial photon technique).8,27,28

In this paper, we have developed a polarized light Monte Carlo algorithm written in the C programming language that tracks the progression of the electric field in scattering and absorbing media containing birefringent materials. This code enables the modeling of tissue as a statistically homogeneous continuous random media under the Whittle-Matérn model or as a composition of discrete spherical scatterers under Mie theory. For ease of operator use and data processing, a graphical user interface (GUI) written in MATLAB is used to interact with the underlying C code. In addition, speed improving techniques using message passing interface (MPI) for parallel computing and the semi-analytical approach are employed.

The paper is organized as follows: In Sec. 2, we first provide a summary of the theoretical origin of the CBS peak. We then discuss the methodologies used to compute the coherent spatial reflectance profile, including the computation of the amplitude scattering matrix and Jones N-matrix for implementation of birefringence. In Sec. 3, we describe the methods used in our software to provide improved computational speed, accuracy, and usability. Finally, in Sec. 4 we provide demonstrations of results from our simulations first for purely scattering media and second, for media containing both scattering and linear birefringence.

Coherent Backscattering

Detailed discussion of the nature of the CBS phenomenon can be found in a number of different publications.37,17 Briefly, the experimentally measurable CBS peak shape ICBS(θx,θy) is the Fourier transform of the product of four functions: Display Formula

ICBS(θx,θy)=p(xs,ys)·pc(xs,ys)·s(xs,ys)·c(xs,ys)eik(xssinθx+yssinθy)dxsdys,(1)
where function p(xs,ys) is the spatial impulse-response of multiply scattered light in the exact backscattering direction (i.e., antiparallel to the incident direction), function pc(xs,ys) is the degree of phase correlation between the forward and reverse path of all rays exiting at a particular (xs,ys) separation, function s(xs,ys) is a modulation due to a finite illumination spot size, and function c(xs,ys) is a modulation due to finite spatial coherence of the illumination. Functions p(xs,ys) and s(xs,ys) assume values between 0 and 1, while functions pc(xs,ys) and c(xs,ys) assume values between 1 and 1. Note that in the above notations, the subscript s indicates a relative separation between any two points (i.e., xs=x2x1).

Within the Fourier integral of Eq. (1), functions p(xs,ys) and pc(xs,ys) represent sample dependent properties that we will simulate numerically using electric field Monte Carlo while functions s(xs,ys) and c(xs,ys) are instrumental properties that can be found analytically as follows: Function s(xs,ys) can be calculated as the normalized autocorrelation of the spatial illumination intensity distribution A(x,y) incident on the scattering sample:3,17Display Formula

s(xs,ys)=A(x,y)A(xxs,yys)dxdyA2(x,y)dxdy.(2)
Under the assumptions of the van Cittert-Zernike theorem, function c(xs,ys) can be calculated as the normalized Fourier transform of the angular intensity distribution I(θx,θy) of light incident on the scattering sample29: Display Formula
c(xs,ys)=I(θx,θy)eik(xsθx+ysθy)dθxdθyI(θx,θy)dθxdθy.(3)
Equations for functions s(xs,ys) and c(xs,ys) are presented here for completeness. However, in the remainder of this paper we will focus on the shapes of functions p(xs,ys) and pc(xs,ys) which are calculated using our Monte Carlo code. Numerical MATLAB calculation of functions s(xs,ys) and c(xs,ys) using top-hat distributions for functions A(x,y) and I(θx,θy), respectively, are posted on our laboratory website.30

Numerical Calculation of Functions p(xs,ys) and pc(xs,ys) with Monte Carlo

Electric field Monte Carlo simulations provide a numerical solution to the vector radiative transport equation in situations where an analytical solution is either difficult or impossible to derive. The basic structure of a light Monte Carlo code is well known. In brief, photons are first injected into a material and allowed to propagate according the material properties until the photon is either absorbed or exits the medium. Along the way, the polarization state of each photon is tracked and any number of parameters characterizing light propagation (e.g., reflectance, absorbance, time of flight, etc.) can be calculated.

For simulation of CBS, the interference between time-reversed photon path-pairs must be considered. One way to rigorously treat this coherence property is using the Jones calculus formalism to track the evolution of the electric field as it encounters optical components (e.g., polarizers), refractive index contrasts that induce scattering, and birefringent materials.31 In our simulations, which use a heavily modified version of the Stokes vector meridian plane Monte Carlo code written by Ramella-Roman et al.,32 the electric field undergoes four linear transformations for each scattering event: Display Formula

[EE]=[cosγsinγsinγcosγ][m1m4m3m2][S2(θ,ϕ)S3(θ,ϕ)S4(θ,ϕ)S1(θ,ϕ)][cosϕsinϕsinϕcosϕ][EE],(4)
Display Formula
E˜=R(γ)MS(θ,ϕ)R(ϕ)E˜,(5)
where R(ϕ) is a rotation from the meridian plane into the scattering plane, S(θ,ϕ) is the amplitude scattering matrix, M is the transformation due to propagation through birefringent media, and R(γ) is a rotation from the scattering plane back into the meridian plane. Note that the notation for the elements of matrices S(θ,ϕ) and M follow the conventions of van de Hulst33 and Jones,34 respectively. Computation of the matrix elements of S(θ,ϕ) and M will be discussed in the following two subsections.

For a multiple scattering sample, the transformations in Eq. (5) accumulate for each scattering event until the light escapes from the medium: Display Formula

E˜=Rn(γ)MnSn(θ,ϕ)Rn(ϕ)R1(γ)M1S1(θ,ϕ)R1(ϕ)M0E˜,=M¯E˜,(6)
where the subscript in Eq. (6) indicates the numbers of scattering events and M¯ is the effective complex scattering matrix for a ray scattered n times. Each matrix element of M¯ can assume an infinite number of different complex values depending on the sample composition, geometry, and photon visitation history.

Generalizing the matrix M¯ for the forward propagating path (denoted with subscript ) we can write: Display Formula

M¯=[ar+iaibr+ibicr+icidr+idi],(7)
where a, b, c, and d are arbitrary variables with subscripts indicating the real (r) and imaginary parts (i) of each term. In simulations of CBS, it is necessary to find both the forward propagating and reverse propagating (denoted with subscript ) matrices in order to accurately calculate the degree of interference between the two rays. Fortunately, once we have calculated M¯, M¯ can be found trivially according to the reciprocity theorem as:25Display Formula
M¯=[ar+iaicricibribidr+idi].(8)
After M¯ and M¯ have been calculated, the electric fields exiting the medium for the forward and reverse paths can be found by multiplying the matrix M¯ by the Jones matrix for the incident polarizer, Pincident. Display Formula
E˜=M¯PincidentE˜,E˜=M¯PincidentE˜.(9)
In order to calculate function p(xs,ys), we first convert E˜ into observable intensities specified by the Stokes parameters:35Display Formula
I=EE*+EE*,Q=EE*EE*,U=EE*+EE*,V=i(EE*EE*).(10)
The unnormalized function p(xs,ys) can then be found for various polarization channels by incoherently summing the appropriate combination of Stokes parameters for the forward propagating path over all multiply scattered photon realizations exiting in the exact backscattering direction. In our simulations, we calculate function p(xs,ys) for four different polarization channels: linear co-polarized xx, linear cross-polarized xy, helicity preserving ++, and orthogonal helicity +. These can be found as: Display Formula
pxx(xs,ys)=nW·[1+Q(xs,ys)/I(xs,ys)],pxy(xs,ys)=nW·[1Q(xs,ys)/I(xs,ys)],p++(xs,ys)=nW·[1+V(xs,ys)/I(xs,ys)],p+(xs,ys)=nW·[1V(xs,ys)/I(xs,ys)],(11)
where W is the photon weight and n is the number of photons. In addition to scoring the intensities in Cartesian coordinates, we also score the intensities as p(rs,z) where rs is the exit radius and z is the maximum penetration depth. Mathematically, these grids are related by p(rs,z)=02πp(x=rscosϕ,y=rssinϕ,z)dϕ.

In order to maintain conservation of energy, we score all photons that exit outside of our grid or are single scattered in the peripheral pixels of each grid. In this paper, we normalize function p(xs,ys) such that the integral over all spatial coordinates plus the single scattered intensity for unpolarized light equals 1: Display Formula

poo(xs,ys)dxsdys+SSoo=1,(12)
where the subscript oo indicates unpolarized illumination and collection and SS is the single scattered intensity. In terms of the component polarizations, we can calculate the unpolarized case by summing the four components: poo=pxx+pyy+pxy+pyx=p+++p−−+p++p+.

Since CBS is a coherence phenomenon, it is required that both the polarization and phase of the time-reversed photon path-pairs are the same in order for interference to occur. As a result, function p(xs,ys) is never independently measurable using CBS and instead can only be measured as the product of p(xs,ys)·pc(xs,ys). One caveat of this statement is that for the polarization preserving channels (e.g., xx and ++), the reciprocity theorem requires that the forward and reverse propagating paths exit with the same accumulated phase and pc(xs,ys)=1. For the polarization nonpreserving channels (e.g., xy and +) the degree of coherence DOC for a single time-reversed path-pair is found as:8Display Formula

DOC=2R[E(xs,ys)E*(xs,ys)]|E(xs,ys)|2+|E(xs,ys)|2.(13)
The product of functions p(xs,ys)·pc(xs,ys) for the orthogonal polarization channels is then found as: Display Formula
pxy(xs,ys)·pcxy(xs,ys)=nW·DOCxy·[1Q(xs,ys)/I(xs,ys)],p+(xs,ys)·pc+(xs,ys)=nW·DOC+·[1V(xs,ys)/I(xs,ys)].(14)

Computation of the Amplitude Scattering Matrix S(θ,ϕ) and Phase Function P(θ,ϕ)

The amplitude scattering matrix S(θ,ϕ) in Eq. (5) succinctly summarizes the effects that a single scattering event has on the transformation of the incident electric field. For an arbitrary scattering media, each of the four matrix elements is independent of the others and are a function of both the polar angle θ and the azimuthal angle ϕ. However, a number of simplifications to the elements of S(θ) can be made through knowledge about the geometry and composition of the scattering material. In media composed of spherically symmetrical scatterers (e.g., spheres or statistically homogeneous random media), the matrix elements S3 and S4 are identically equal to zero and the matrix elements S1 and S2 are solely functions33 of θ. In this paper, we implement two spherically symmetrical models for the composition of the scattering material: first, discrete spheres and second, statistically homogeneous continuous random media with refractive index distributions specified by the Whittle-Matérn family of correlation functions.1,17,24,36,37

Discrete spheres—Mie theory

For a scattering media composed of discrete spherical particles, the matrix elements S1(θ) and S2(θ) can easily be calculated numerically according to Mie theory.33 In our simulation, we implemented a vectorized version of the Mie code written by Mätzler following the formalism put forth by Bohren and Huffman.35,38 Once the scattering amplitudes are calculated, the shape of the phase function P(θ,ϕ) can be determined as: Display Formula

P(θ,ϕ)=S11(θ)·Io+S12(θ)·(Qocos2ϕ+Uosin2ϕ)02π0π[S11(θ)·Io+S12(θ)·(Qocos2ϕ+Uosin2ϕ)]sinθdθdϕ,(15)
where [Io, Qo, Uo, Vo] are the Stokes parameters for incident illumination and S11(θ) and S12(θ) are elements of the Mueller scattering matrix:35Display Formula
S11(θ)=|S2(θ)|2+|S1(θ)|22,S12(θ)=|S2(θ)|2|S1(θ)|22.(16)

Continuous random media—Born approximation

The Born approximation, also known as Rayleigh-Gans-Debye theory, provides a simplification to scattering theory which enables solutions for otherwise intractable problems (e.g., scattering in continuous random media). Provided that the fluctuations in refractive index are sufficiently weak such that the phase shift of the incident wave is small, we can approximate the total field within the scattering medium as the incident field. The scattered field can then be computed as the coherent summation of the dipole scattering pattern from each position within the medium. In the following paragraphs, we begin by describing the calculation of the scattering amplitude matrix under the Born approximation. We then apply these calculations to computation of scattering for a continuous random media defined under the Whittle-Matérn model.

The scattering amplitude matrix for an infinitesimal dipole can be found as:33,35Display Formula

dS(θ)=ik3·dα·[cosθ001],(17)
where dα is the differential polarizability and k is the wavenumber. For weak refractive index contrast dα, can be approximated as: Display Formula
dαnΔ2πdV,(18)
where nΔ is the excess refractive index within the differential volume dV relative to the mean medium refractive index [(ndV/nmedium)1]. To find the total scattered field for the ensemble scattering medium, we coherently sum the contribution from each dipole over the entire volume: Display Formula
S(θ)=ik3·VnΔ(r)2πeiks·rdV·[cosθ001],(19)
where |ks|2|k|sinθ2. The term eiks·r accounts for the relative phase difference of the field originating from different positions within the medium and observed in the direction ks. Upon inspection, we see that the integral in Eq. (19) is essentially the three dimensional Fourier transform of nΔ(r).

For continuous random media, there is no elementary scattering particle with decipherable boundaries. As such, scattering must be defined in terms of the amplitude scattering matrix per unit volume polarizability s(θ): Display Formula

s(θ)=ik3·f(ks)·[cosθ001],(20)
where function f is the scattering form factor of the particular scatterer under investigation. Assuming spherical symmetry for the scattering medium, function f can be calculated through reduction of the three dimensional Fourier transform in Eq. (19): Display Formula
f(ks)=2α0Mn(r)rsin(ksr)ksdr,(21)
Display Formula
α=20r2Mn(r)dr,(22)
where the statistical function Mn(r) is the particulate equivalent medium and replaces nΔ(r) in Eq. 19. Conceptually, Mn(r) can be thought of as the effective particle which gives the same scattered intensity as the random process described by nΔ(r) [see Fig. 1(a)]. Note that Mn(r) uses the scalar r (implying spherical symmetry), while nΔ(r) uses the vectorial r (implying lack of symmetry).

Graphic Jump LocationF1 :

Numerical validation of our treatment of scattering using function Mn(r). (a) Randomly generated 3-D medium under the Whittle-Matérn model for D=3. (b) Comparison of function |f| calculated numerically and analytically for the medium shown in panel a.

One attractive model for Mn(r) in a continuous random scattering medium originates from the Whittle-Matérn family of correlation functions.37 Under this model, the distribution of refractive index fluctuations is defined through the medium’s refractive index correlation function Bn(rs): Display Formula

Bn(rs)=An(rslc)D32KD32(rslc),(23)
where Kν is the modified Bessel function of the second kind with order ν, lc describes the length-scale of tissue heterogeneity, An is the fluctuation strength, and D is a parameter which determines the shape of the distribution (e.g., Gaussian, stretched exponential, exponential, and power law distributions).37 Implementation of the Whittle-Matérn model provides the flexibility to mimic the distributions of refractive index expected for a wide range of different biological tissue types.18,24,37 It should be noted that when D=3, this model predicts a scattering phase function that is identical to the commonly used Henyey-Greenstein case. Figure 1(a) shows a single realization of the three dimensional excess refractive index distribution for D=3.

According to the convolution theorem, the random process nΔ(r) is related to Bn(rs) and Mn(r) by: Display Formula

Bn(rs)=F1[|F[nΔ(r)]|2]=F1[|F[Mn(r)]|2],(24)
where the symbol F indicates the Fourier transform operation. Inverting Eq. (24), we can derive Mn(r): Display Formula
Mn(r)=224π3/4AnΓ(D/2)Γ(D/4)lc3/2(rlc)D64KD64(rlc).(25)
It should be noted that Mn(r) is a statistical function that provides the same information as Bn(rs) but allows for the calculation of the scattered electric field needed for electric field Monte Carlo. One caveat of the previous statement is that the phase information is not represented and so we can only correctly calculate the modulus of function f. Still, the shape of the scattering phase function which is proportional to the square of function |f| is unaffected by this distinction. Performing the integrations in Eqs. (21) and (22), under the Whittle-Matérn model function |f| can be found as: Display Formula
|f(ks)|=(1+ks2lc2)D/4.(26)
Equation (26) is an accurate estimate of scattering provided σn2(klc)21.39

Figure 1 provides a numerical validation of our treatment of scattering using function Mn(r). Figure 1(a) shows a three dimensional random medium generated using a publicly available code for D=30.40 Using this medium, we calculated function |f| by numerically computing the modulus of the three-dimensional (3-D) Fourier transform in Eq. (19), rotationally averaging the resulting function, and finally normalizing by the first point. Figure 1(b) shows a close match between this numerically calculated function and the analytical result from Eq. (26).

Computation of the Jones M-Matrix for Implementation of Birefringence

A number of biological tissues exhibit some combination of linear birefringence due to structural alignment and optical activity due to the presence of chiral molecules. Since the matrix transformations of the electric field are noncommutative, the order in which these effects are applied can potentially alter the observable signal.41 In this paper, we assume a stochastic medium in which there is no preferential order for the effects of linear birefringence and optical activity. In order to combine these two effects into a single matrix operation, the Jones N-matrix formalism can be utilized.34

The Jones N-matrix represents the transformation of the electric field for an infinitesimally small propagation path length ds. Following the derivation by Jones, this enables the combination of multiple polarization-altering effects (e.g., birefringence and dichroism) into a single M-matrix which represents a transformation over the entire path length. The N-matrix for a specimen containing both linear birefringence and optical activity can be written as41: Display Formula

N=(dMds)M1=[n1n4n3n2]=[i·g0ωωi·g0],(27)
where the elements ±i·g0 represent changes in the phase between the two orthogonal linear polarization states that occurs due to linear birefringence and the elements ±ω represent rotations of polarization state due to optical activity.

The parameter g0 can be found as:34Display Formula

g0=πλΔnb(θb),(28)
where λ is the wavelength of light within the medium and Δnb is dependent on the ordinary refractive index no, the extraordinary refractive index ne, and the angle between the photon trajectory and the optic axis θb: Display Formula
Δnb(θb)=nonene2cos2θb+no2sin2θbno.(29)
The photon trajectory is defined by the direction cosine e^prop and the optic axis41 by the “birefringence unit vector” b. In order to describe a particular media’s material property (which must be independent of θb), we will refer to its value of Δnb,max=Δnb(π/2)=neno.

The parameter ω can be found as the product of the chiral molecule’s specific rotation [α]λT at temperature T (degrees Celsius) and its concentration ρ: Display Formula

ω=[α]λT·ρ.(30)
Both g0 and ω are calculated in units of radians per centimeter.

The M-matrix elements corresponding to the convention in Eq. (4) can be calculated for a given path length s as: Display Formula

m1=i·g0sinh(QN·s)QN+cosh(QN·s),m2=i·g0sinh(QN·s)QN+cosh(QN·s),m3=ωsinh(QN·s)QN,m4=ωsinh(QN·s)QN,(31)
where QN is found as: Display Formula
QN=g02ω2.(32)
As defined in Eq. (31), the M-matrix elements must first be rotated by an angle β to a reference frame in which the parallel component of the electric field e^ is aligned with the maximum refractive index difference b before application to Eq. (4). Wood et al. provide a nice discussion and schematic of the steps involved in this rotation.41 Mathematically, the vector b can be found as b=e^prop×(b×e^prop). The angle β can then be found by vector multiplication between e^ and b. The properly rotated M-matrix is then found as: Display Formula
M=[cos(β)sin(β)sin(β)cos(β)][m1m4m3m2][cos(β)sin(β)sin(β)cos(β)].(33)

The general details of meridian plane Monte Carlo simulations as well as an open source code are discussed in the publication by Ramella-Roman et al.32 Here, we implement a heavily modified version of this code to model a pencil beam normally incident on a nonabsorbing semi-infinite medium with refractive index matching at the boundary. The rationale for assuming index matching at the boundary is the excellent agreement between such a code and experimental measurements.3,17 The collection geometry to model p(xs,ys) is a square grid with x and y resolution of 0.01·μs* and extent ranging from 5·μs* to 5·μs* in both x and y. Additionally, p(rs,z) is recorded as a square grid with rs and z resolution of 0.01·μs* and extent ranging from 0 to 5·μs* in both rs and z. As a reminder, both p(xs,ys) and p(rs,z) record only photons that exit the medium exactly in the direction of the surface unit normal vector (i.e., antiparallel to the incident pencil beam with numerical aperture=0).

In addition to the theoretical formalisms for scattering in continuous random media and propagation in birefringent materials discussed in Sec. 2, we have implemented three major speed and usability improvements to our Monte Carlo code. These include the implementation of, first, the semi-analytical approach, second, a message passing interface (MPI), and third, a graphical user interface (GUI) written in MATLAB.

Semi-Analytical Approach

In order to accurately model CBS, function p(xs,ys) must be calculated for reflected light that is antiparallel to the incident direction. However, under the traditional Monte Carlo approach, only an infinitesimally small number of photons will exit the scattering medium exactly in this direction. As a result, in order to achieve any computational efficiency, it is necessary to collect all photons that exit the medium within some finite collection angle around the backscattering direction. Unfortunately, this generates a trade-off between computational accuracy and efficiency since the spatial distribution of light reflected at different angles is not constant. As a compromise, previous publications have found that collection of photons exiting the medium within an angle of 10 degrees around the backscattering direction produce no noticeable deviations from the theoretical shape of function p(xs,ys).18,23 Still, under the traditional approach only a relatively small number of the total photons contribute to the final result.

In order to improve the computational efficiency of the traditional approach to Monte Carlo simulations, we implement the semi-analytical approach (also known as the partial photon technique).8,27,28,42,43 Using this method, a portion of scattered intensity is collected after a photon reaches each new position within the medium. This “partial photon” intensity Ipartial is calculated by multiplying the probability that the photon is scattered in the direction of the surface F(θs,0) by the probability that it will reach surface according to the Beer-Lambert law e(μs+μa)z: Display Formula

Ipartial=F(θs,0)·e(μs+μa)z,(34)
where θs is the angle between e^prop and the surface normal vector and z is the distance to the surface. Modifying Eq. (11) under the semi-analytical approach, we find function p(xs,ys) for different the different polarization channels as: Display Formula
pxx(xs,ys)=nIpartial·W·[1+Q(xs,ys)/I(xs,ys)],pxy(xs,ys)=nIpartial·W·[1Q(xs,ys)/I(xs,ys)],p++(xs,ys)=nIpartial·W·[1+V(xs,ys)/I(xs,ys)],p+(xs,ys)=nIpartial·W·[1V(xs,ys)/I(xs,ys)],(35)
where the index of summation n now represents the number of scattering events. The product of functions p(xs,ys)·pc(xs,ys) for the orthogonal polarization channels can be found in analogy with Eq. (36) as: Display Formula
pxy(xs,ys)·pcxy(xs,ys)=nIpartial·W·DOCxy·[1Q(xs,ys)/I(xs,ys)],p+(xs,ys)·pc+(xs,ys)=nIpartial·W·DOC+·[1V(xs,ys)/I(xs,ys)].(36)
Scoring photons in this way enables both computational efficiency through collection of information at each scattering event and accuracy through collection of intensity in the exact backscattering direction.

Figure 2 shows a visual comparison of the noise level between the semi-analytical (panel a) and traditional (panel b) technique in a Rayleigh scattering sample simulated for the xx polarization channel. Figure 2(c) compares the pxx(rs·μs*) and pxx(z·μs*) distributions achieved by summing pxx(rs·μs*,z·μs*) over rows and columns, respectively. Note that since p(rs,z) scales linearly with the reduced scattering coefficient μs*=1/ls*, we show each distribution as a function of the dimensionless parameters rs·μs* and z·μs*. Quantitative comparison of the computational efficiency between the two techniques is shown in Fig. 2(d) and calculated by taking the ratio of the number of traditional photons over the number of semi-analytical photons needed to achieve the same simulation noise level. For isotropic scattering (i.e., g=0), the semi-analytical approach is an exceptional 200 times faster than the traditional approach. The efficiency of the semi-analytical approach decreases for increasing anisotropy factor, but remains superior for values of g<0.96 (which encompasses the majority of tissue types44). In order to ensure optimal efficiency for each simulation, our code scores photons using both the traditional and semi-analytical technique.

Graphic Jump LocationF2 :

Comparison between the semi-analytical and traditional Monte Carlo method in a sample of Rayleigh scatterers simulated for the xx polarization channel. Function pxx(rs·μs*,z·μs*) for (a) the semi-analytical method and (b) the traditional method. (c) Functions pxx(rs·μs*) and pxx(zs·μs*) achieved by summing pxx(rs·μs*,z·μs*) over rows and columns, respectively. The symbols indicate the traditional technique while the solid line indicates the semi-analytical approach. (d) Computational efficiency measured by taking the ratio of the number of traditional photons over the number of semi-analytical photons needed to achieve the same simulation noise level. These simulations utilize the Whittle-Matérn model with D=3.

MPI Implementation

One of the computational benefits of performing Monte Carlo simulations is that the noise level scales inversely with the number of recorded photons. Due to the stochastic nature of the Monte Carlo method, information from each photon is independent and calculations from an indeterminate number of processors can be linearly combined to reduce the noise variance. In order to take advantage of this capability, we have implemented MPI into our simulations. This allows multiple processors to independently calculate a predetermined number of photon histories, and subsequently combine the results after each processor has finished. A simulation of Rayleigh scattering for 108 photons on a dual quadcore workstation (2.1 GHz AMD Opteron Processor 2352) can be run in less than 1.5 h.

MATLAB GUI

Within the engineering community, MATLAB is the preferred platform on which to analyze data. However, to perform highly repetitive calculations in the minimum amount of time, the C programming language is advantageous. In order to combine the usability of MATLAB and speed of C code, we have implemented a software program that integrates these two environments. User interaction with the simulation is carried out through the MATLAB GUI shown in Fig. 3. After specifying the desired parameters, a simulation can be imported into a C code environment for rapid calculation of functions p(xs,ys) and pc(xs,ys) on multiple processors with a single button click. After the simulation is completed, all relevant data is automatically compiled and saved into a single MATLAB format file under a user specified file name.

Graphic Jump LocationF3 :

MATLAB GUI for performing Monte Carlo simulations of CBS. (a) Specification of the general Monte Carlo parameters. (b) Specification of scattering model to be implemented. (c) Specification of the birefringence properties of the sample.

The MATLAB GUI consists of three main panels used to specify the relevant simulation parameters. The panel in Fig. 3(a) allows the user to specify the general Monte Carlo parameters such as illumination polarization, photon and processor numbers, optical coefficients μs* and μa, as well as grid size and discretization. The panel in Fig. 3(b) enables simulations to be carried out using either a discrete sphere model under Mie Theory or a continuous distribution of refractive index under the Whittle-Matérn model as described in Sec. 2.3. The panel in Fig. 3(c) allows the user to set birefringence parameters such as the ordinary and extraordinary refractive index, the optical rotation, and the birefringence unit vector. Additionally, the GUI provides the capability to plot a summary of the simulation data as well as export data to and compile data from a remote computer cluster. Further specific details on the operation of the GUI can be found in a user manual posted along with our code.30

In this section, we provide demonstrations of the results of our simulation first for purely scattering samples and then for scattering samples containing linear birefringence. These results are presented in spatial coordinates rather than the conventional way of displaying CBS data in angular coordinates. The rationale for presentation in this way is that we believe understanding light transport in terms of spatial coordinates is more intuitive than in angular coordinates. Although CBS measurements must be acquired in angular coordinates, according to Eq. (1) a simple inverse Fourier transform of these results provides easily interpretable information about how light is transported through biological tissue.3 Additionally, we would like to stress the comparison between CBS and conventional diffuse reflectance measured in the exact backscattering direction.

Trends for Purely Scattering Samples

For any particle form factor (spheres, continuous random media, etc.), the shape of the differential scattering cross-section converges to that of Rayleigh scattering when the characteristic length-scale of the particle is much smaller than the wavelength of illumination. Because of this common thread between all scattering form factors, we begin by analyzing the shape of the CBS peak for Rayleigh scattering.

One way to quantify the shape of the CBS peak is through the enhancement factor E or the relative height of the peak at θ=0 divided by the total unpolarized incoherent intensity. Mathematically, this can be found as: Display Formula

Eν=0pν(rs)·pcν(rs)d(rs),(37)
where the subscript ν indicates the specific polarization channel under analysis (e.g., xx, xy, etc.). Note that the value calculated in Eq. (37) is for plane wave illumination with infinite spatial coherence length (i.e., c(rs)=s(rs)=1).

Table 1 contains a summary of the values of E in the various polarization channels for nonabsorbing Rayleigh scatterers with index matching at the boundary and results rounded to the fourth decimal place. Comparison with the benchmark values calculated by Mishchenko as well as Amic et al. using two separate numerical techniques show errors below 2% in each case.45,46 Additionally, the value of Exx is in agreement with a Monte Carlo code which takes an alternative semi-analytic approach.47 We note that the values given in Refs. 45 and 46 are converted into the normalization used in this paper before display in Table 1.

Table Grahic Jump Location
Table 1Values of E in various polarization channels for nonabsorbing Rayleigh scatterers with index matching at the boundary. The values given for Refs. 45 and 46 are converted into the normalization used in this paper before display.

In the idealized scalar case E would be exactly equal to 1, meaning that the coherent intensity is the same magnitude as the incoherent intensity. However, for real electromagnetic waves Eoo is less than 1 due to first, single scattering and second, partially reversible photon paths (i.e., paths in which DOC1). For the xx and ++ polarization preserving channels, function pc(rs) is identically equal to 1. Because of this, Exx and E++ are nearly one quarter of the total unpolarized incoherent intensity. For the xy and + orthogonal polarization channels, Exx and E++ are greatly reduced due to the decay of function pc at large values of rs.

As the characteristic length-scale of the particle approaches the wavelength of illumination, the specific scattering form factor begins to dominate the shape of the differential scattering cross-section. Under Mie theory, the characteristic length-scale is the radius of the spherical particle, a. Figure 4 shows E as a function of the dimensionless size parameter ka. Figure 4(a) shows these trends for a scattering medium with a relative refractive index m=nsphere/nmedium corresponding to polystyrene microspheres suspended in water. For very small ka, the results converge to the values of Rayleigh scattering given in Table 1. As ka increases, the trends in each polarization channel exhibit an oscillatory pattern due to the spherical form factor. Interestingly, the trends shown in Fig. 4(a) remain essentially the same with the reduced refractive index contrast shown in Fig. 4(b) and 4(c). From these results it can be concluded that refraction of the light wave within the scattering particle has a smaller effect on the shape of the CBS peak than the underlying spherical particle scattering form factor.

Graphic Jump LocationF4 :

Enhancement factor E versus size parameter ka for various polarization channels under Mie theory. (a) Trends for relative refractive index m=nsphere/nmedium corresponding to aqueous suspension of polystyrene microspheres. (b) Trends for m=1.1. (c) Trends for m=1.001.

Under the Whittle-Matérn model, the characteristic length-scale is the parameter lc. Figure 5 shows E as a function of the dimensionless parameter klc for various values of D. Similar to the Mie theory results above, for very small klc the values converge to those of Rayleigh scattering. For larger values of klc, E exhibits a smooth change trend due to the smoothly decaying form factor underlying scattering in the Whittle-Matérn model.

Graphic Jump LocationF5 :

Enhancement factor E versus klc for various polarization channels under the Whittle-Matérn model. (a) Trends for D=2. (b) Trends for D=3. (c) Trends for D=4.

Effects of Birefringence

To demonstrate the general effects of biological birefringence on the shape of the CBS peak, we once again turn to the case of Rayleigh scattering. To study these effects within the biologically relevant regime,48 we performed simulations for values of Δnb,max ranging from 0 to 1×103 with a birefringence unit vector oriented along the x-axis (i.e., b=[1,0,0]).

Figure 6 demonstrates the effects of linear birefringence on the spatial reflectance profiles for linear polarized illumination and collection with the arrows indicating increasing values of Δnb,max. Noting that pcxx=1 for all length-scales, Fig. 6(a) demonstrates that the pxx and pxy distributions which would be measured using conventional incoherent techniques exhibit minimal sensitivity to the presence of birefringence. However, a noteworthy change in the coherent pxy·pcxy distribution (measurable using CBS) occurs even for small values of Δnb,max. The explanation for this observation is that the presence of birefringence reduces the reversibility of the path travelled by each multiply scattered photon, resulting in a more sharply decaying function pcxy shown in Fig. 6(b). Mathematically, the additional polarization rotations imparted by the presence of birefringence reduce the symmetry of the effective complex scattering matrix M¯ and therefore cause the forward and reverse propagating rays to be less correlated with each other on average.

Graphic Jump LocationF6 :

Effects of linear birefringence on the spatial reflectance profiles under linear polarized illumination and collection for a sample of Rayleigh scatterers. (a) The normalized intensity distributions pxx·pcxx and pxy·pcxy which are measurable using CBS. The dashed line indicates the shape of the incoherent pxy distribution which is not directly measurable using CBS. (b) Function pcxy. Due to the full reversibility of the xx polarization channel, pcxx=1 for all length-scales. In both panels, the arrow indicates increasing Δnb,max=0, 4×104, 7×104, and 1×103.

Figure 7 shows the effects of birefringence on the spatial reflectance profiles for circularly polarized illumination and collection with increasing Δnb,max. Similar to the effect described for linear polarization, function pc+ shown in Fig. 7(b) is more strongly decaying for large values of Δnb,max. This results in a strong attenuation of function p+·pc+ relative to function p+ at large length-scales. Interestingly, for circular polarization both functions p++ and p+ are altered at short length scales due to the presence of birefringence. The reason for this observation can be understood by considering the rotation of polarization that occurs due to birefringence prior to the first scattering event. As the circularly polarized photon enters the medium and propagates to the first scattering event, it encounters the birefringence material and becomes elliptically polarized. Thus, when the photon encounters the first scattering event, the shape of the scattering phase function is altered relative to the absence of birefringence case (note: although we use the first scattering event as conceptual description, the described effect will occur at each scattering event). The degree of alteration in the shape of the phase function is directly related to the strength of the birefringence; the stronger the birefringence, the larger the alteration in the phase function. Because of this change in the phase function, the shapes of p++ and p+ are altered at short length scales due to the presence of birefringence. Further demonstration of this effect is seen in Fig. 8 with discussion found below.

Graphic Jump LocationF7 :

Effects of linear birefringence on the spatial reflectance profiles under circularly polarized illumination and collection for a sample of Rayleigh scatterers. (a) The normalized intensity distributions p++·pc++ and p+·pc+ which are measurable using CBS. The dashed line indicates the shape of the incoherent p+ distribution which is not directly measurable using CBS. (b) Function pc+. Due to the full reversibility of the ++ polarization channel, pc++=1 for all length-scales. In both panels, the arrow indicates increasing Δnb,max=0, 4×104, 7×104, and 1×103.

Graphic Jump LocationF8 :

Alterations in the shape of p(xs,ys)·pc(xs,ys) due to birefringence for different polarization channels. To emphasize the shape of the various distributions, each image shows log10[p(xs,ys)·pc(xs,ys)] in the same intensity scale. (Rows) Each row corresponds to the polarization channel specified on the far right of the figure: top row shows xx polarization, second row shows ++ polarization, third row shows xy polarization, and the bottom row shows + polarization. (Columns) The left column shows the log10[p(xs,ys)·pc(xs,ys)] distributions for Δnb,max=0 and the middle column shows the distributions for Δnb,max=1×103. The right column shows the angular distribution found by converting to polar coordinates, summing over radius, and normalizing by the mean.

Trends of E as a function of physiological levels of Δnb,max are shown in Fig. 9(a). Beginning with Δnb,max=0, the values for E are the same as in Table 1. As Δnb,max increases, the value of E weakly increases for the polarization preserving channels and more strongly decreases for the orthogonal polarization channel. The percent change in E for the various polarization channels relative to the absence of birefringence is shown in Fig. 9(b). For Δnb,max=1×103 (on the order of the highest value found in biological tissue), there is an approximately 10% increase in E for the polarization preserving channels attributable mainly to the change in the shape of the phase function. On the other hand, in the orthogonal polarization channels E decreases by 37% for the + channel and 50% for the xy channel due the destruction of reversibility attributable to the presence of birefringence.

Graphic Jump LocationF9 :

Changes in E as a function of Δnb,max for various polarization channels. (a) Trends of versus E. (b) Percent change in E from the absence of birefringence case versus Δnb,max.

Figure 8 shows the comparison in the shape of p(xs,ys)·pc(xs,ys) between Δnb,max=0 (left column) and Δnb,max=1×103 (center column). The angular distribution shown in the right column is found by converting the Cartesian coordinates into polar coordinates, summing over radius, and normalizing by the mean. Each of the four rows corresponds to the polarization channel shown at the far right of the figure. Starting from the top row we show the xx, ++, xy and + polarization channels. In the top row, we observe a minimal change in the shape of pxx(xs,ys), while in the second row the shape of p++ is noticeably elongated along the 135°/315° direction. The explanation for ++ channel is that the incident right circular illumination becomes elliptically polarized along the 45-deg/225-deg direction as it propagates through the birefringent crystal. This causes a decreased probability of scattering in the direction of the ellipticity (known as the dipole factor)17,37 which in turn results in more light intensity being scattered orthogonal to this direction. The direction of the elongation depends on the magnitude and sign of Δnb,max as well as the helicity of incident light. In the third row, the increased value of function Δnb,max shrinks the spatial extent of pxy(xs,ys)·pcxy(xs,ys) due to a more strongly decaying function pcxy(xs,ys) as described above. Finally, in the last row the shape of p+(xs,ys)·pc+(xs,ys) is altered from a symmetric shape to an oblong “X” shape with increasing birefringence. The reason for this shape is the conversion of the incident circular light to an elliptical state that mimics the “X” pattern of the xy polarization channel.

In this paper, we presented the methodologies needed for performing rigorous electric field tracking Monte Carlo simulations of CBS in biological media containing birefringence. We began by reviewing the dependence of the angular CBS peak on the spatial reflectance profile. Next, we detailed the calculation of the scattering amplitude matrix for continuous random media under the Whittle-Matérn model and the calculation of the Jones N-matrix for light propagation in birefringent media. We then described the particular computational methods used to improve the speed and usability of our code. Using a dual quadcore workstation (2.1 GHz AMD Opteron Processor 2352), a simulation of Rayleigh scattering for 108 photons can be run in less than 1.5 h. Future developments to our code can implement GPU acceleration and peer-to-peer networking to further improve the speed of our simulations.49 Finally, we provided demonstrations of the results from simulations for purely scattering samples and samples containing scattering and linear birefringence. These simulation results demonstrate a strong dependence of the shape of the coherent spatial reflectance profile on polarization, illumination wavelength, scattering form factor, and degree of linear birefringence.

For samples containing linear birefringence, the shape of the coherent spatial reflectance profile was altered due to, first, changes in the shape of the phase function attributable to changes in polarization state and, second, a more quickly decaying function pc(xs,ys) for the orthogonal polarization channels caused by loss of reversibility between the forward and reverse propagating rays. The change in the shape of the phase function causes the two circular polarization channels to exhibit large differences in azimuthal distribution, while for the two linear polarization channels the azimuthal distribution remains in large part the same. The loss of reversibility in the orthogonal polarization channels causes the enhancement factor to decrease by 37% for the + channel and 50% for the xy channel for a sample of Rayleigh scatterers with Δnb,max=1×103 and birefringence vector oriented along the x-axis.

The speed and usability of this software makes it ideal for studying the alteration in the shape of the CBS peak for different biological sample compositions. One application of this software will be to study early alterations in biological tissue caused by carcinogenesis. Recent work suggests that one harbinger of early carcinogenesis is a profound reorganization of the extracellular matrix which causes an alignment of the collagen fibers.50 This collagen fiber alignment results in an increase in the degree of birefringence which may be detectable using CBS. Our software will help elucidate the sensitivities of the CBS peak to nanoscale mass density fluctuations as well as levels of birefringence caused by extracellular matrix reorganization. While the results presented in this paper primarily focus on the CBS phenomenon, we note that these simulations are also accurate for conventional diffuse reflectance measurements in the backscattering direction.

This study was supported by National Institute of Health Grant Nos. RO1CA128641 and R01EB003682. A.J. Radosevich is supported by a National Science Foundation Graduate Research Fellowship under Grant No. DGE-0824162.

Gomes  A. J. et al., “In vivo measurement of the shape of the tissue-refractive-index correlation function and its applicationto detection of colorectal field carcinogenesis,” J. Biomed. Opt.. 17, (4 ), 047005  (2012). 1083-3668 CrossRef
Wang  X., Wang  L. V., “Propagation of polarized light in birefringent turbid media: a Monte Carlo study,” J. Biomed. Opt.. 7, (3 ), 279 –290 (2002). 1083-3668 CrossRef
Radosevich  A. J. et al., “Measurement of the spatial backscattering impulse-response at short length scales with polarized enhanced backscattering,” Opt. Lett.. 36, (24 ), 4737 –4739 (2011). 0146-9592 CrossRef
Kuga  Y., Ishimaru  A., “Retroreflectance from a dense distribution of spherical particles,” J. Opt. Soc. Am. A. 1, (8 ), 831 –835 (1984). 0740-3232 CrossRef
Tsang  L., Ishimaru  A., “Backscattering enhancement of random discrete scatterers,” J. Opt. Soc. Am. A. 1, (8 ), 836 –839 (1984). 0740-3232 CrossRef
Wolf  P. E., Maret  G., “Weak localization and coherent backscattering of photons in disordered media,” Phys. Rev. Lett.. 55, (24 ), 2696 –2699 (1985). 0031-9007 CrossRef
Albada  M. P. V., Lagendijk  A., “Observation of weak localization of light in a random medium,” Phys. Rev. Lett.. 55, (24 ), 2692 –2695 (1985). 0031-9007 CrossRef
Lenke  R., Maret  G., Multiple Scattering of Light: Coherent Backscattering and Transmission, Scattering in Polymeric and Colloidal Systems. , Brown  W., Mortensen  K., Eds., pp. 1 –73,  Dover ,  London  (2000).
Dogariu  A., Uozumi  J., Asakura  T., “Enhancement factor in the light backscattered by fractal aggregated media,” Opt. Rev.. 3, (2 ), 71 –82 (1996). 1340-6000 CrossRef
Ishii  K., Iwai  T., Asakura  T., “Polarization properties of the enhanced backscattering of light from the fractal aggregate of particles,” Opt. Rev.. 4, (6 ), 643 –647 (1997). 1340-6000 CrossRef
Wiersma  D. S., van Albada  M. P., Lagendijk  A., “Coherent backscattering of light from amplifying random media,” Phys. Rev. Lett.. 75, (9 ), 1739 –1742 (1995). 0031-9007 CrossRef
Labeyrie  G. et al., “Coherent backscattering of light by cold atoms,” Phys. Rev. Lett.. 83, (25 ), 5266 –5269 (1999). 0031-9007 CrossRef
Vithana  H. K. M., Asfaw  L., Johnson  D. L., “Coherent backscattering of light in a nematic liquid crystal,” Phys. Rev. Lett.. 70, (23 ), 3561 –3564 (1993). 0031-9007 CrossRef
Sapienza  R. et al., “Anisotropic weak localization of light,” Phys. Rev. Lett.. 92, (3 ), 033903  (2004). 0031-9007 CrossRef
Yoo  K. M., Tang  G. C., Alfano  R. R., “Coherent backscattering of light from biological tissues,” Appl. Opt.. 29, (22 ), 3237 –3239 (1990). 0003-6935 CrossRef
Yoon  G., Roy  D. N. G., Straight  R. C., “Coherent backscattering in biological media: measurement and estimation of optical properties,” Appl. Opt.. 32, (4 ), 580 –585 (1993). 0003-6935 CrossRef
Radosevich  A. et al., “Polarized enhanced backscattering spectroscopy for characterization of biological tissues at subdiffusion length-scales,” IEEE J. Sel. Top. Quant. Electron.. 18, (4 ), 1313 –1325 (2012). 1077-260X CrossRef
Turzhitsky  V. et al., “Measurement of optical scattering properties with low-coherence enhanced backscattering spectroscopy,” J. Biomed. Opt.. 16, (6 ), 067007  (2011). 1083-3668 CrossRef
Roy  H. K. et al., “Association between rectal optical signatures and colonic neoplasia: potential applications for screening,” Cancer Res.. 69, , 4476 –83 (2009). 0008-5472 CrossRef
Turzhitsky  V. et al., “Investigating population risk factors of pancreatic cancer by evaluation of optical markers in the duodenal mucosa,” Dis. Markers. 25, (6 ), 313 –321 (2008). 0278-0240 
Akkermans  E., Wolf  P. E., Maynard  R., “Coherent backscattering of light by disordered media: Analysis of the peak line shape,” Phys. Rev. Lett.. 56, (14 ), 1471 –1474 (1986). 0031-9007 CrossRef
Ishii  K., Iwai  T., “Double-scattering approximation theory of enhanced backscatterings of light produced from aggregated particles,” Jpn. J. Appl. Phys. 1. 41, (Part 1, No. 8 ), 5155 –5159 (2002). 0021-4922 CrossRef
Eddowes  M. H., Mills  T. N., Delpy  D. T., “Monte Carlo simulations of coherent backscatter for identification of the optical coefficients of biological tissues in vivo,” Appl. Opt.. 34, (13 ), 2261 –2267 (1995). 0003-6935 CrossRef
Turzhitsky  V. et al., “A predictive model of backscattering at subdiffusion length scales,” Biomed. Opt. Express. 1, (3 ), 1034 –1046 (2010). 2156-7085 CrossRef
Martinez  A. S., Maynard  R., “Faraday effect and multiple scattering of light,” Phys. Rev. B. 50, (6 ), 3714 –3732 (1994). 1098-0121 CrossRef
Sawicki  J., Kastor  N., Xu  M., “Electric field Monte Carlo simulation of coherent backscattering of polarized light by a turbid medium containing mie scatterers,” Opt. Express. 16, (8 ), 5728 –5738 (2008). 1094-4087 CrossRef
Lenke  R., Maret  G., “Magnetic field effects on coherent backscattering of light,” Eur. Phys. J. B. 17, (1 ), 171 –185 (2000). 1434-6028 CrossRef
Lenke  R., Tweer  R., Maret  G., “Coherent backscattering of turbid samples containing large mie spheres,” J. Opt. A—Pure Appl. Opt.. 4, (3 ), 293  (2002). 1464-4258 CrossRef
Born  M., Wolf  E., Principles of Optics: Electromagnetic Theory of Propagation, Interference and Diffraction of Light. , 7th ed.,  Cambridge University Press ,  Cambridge  (1999).
Radosevich  Andrew J., Rogers  Jeremy D., Matlab functions, Biophotonics Laboratory at Northwestern University, http://biophotonics.bme.northwestern.edu/resources/index.html (30 October 2012).
Jones  R. C., “A new calculus for the treatment of optical systems,” J. Opt. Soc. Am.. 31, (7 ), 488 –493 (1941). 0030-3941 CrossRef
Ramella-Roman  J., Prahl  S., Jacques  S., “Three Monte Carlo programs of polarized light transport into scattering media: part I,” Opt. Exp.. 13, (12 ), 4420 –4438 ( Jun 2005). 1094-4087 CrossRef
van de Hulst  H., Light Scattering by Small Particles. ,  Dover ,  New York  (1981).
Jones  R. C., “A new calculus for the treatment of optical systems. VII. properties of the n-matrices,” J. Opt. Soc. Am.. 38, (8 ), 671 –683 (1948). 0030-3941 CrossRef
Bohren  C. F., Huffman  D. R., Absorption and Scattering of Light by Small Particles. ,  John Wiley and Sons ,  New York  (1983).
Guttorp  P., Gneiting  T., “Studies in the history of probability and statistics XLIX on the Matérn correlation family,” Biometrika. 93, (4 ), 989 –995 (2006). 0006-3444 CrossRef
Rogers  J. D., Çapoğlu  R., Backman  V., “Nonscalar elastic light scattering from continuous random media in the born approximation,” Opt. Lett.. 34, (12 ), 1891 –1893 (2009). 0146-9592 CrossRef
Mätzler  C., “Matlab functions for mie scattering and absorption,” Research Report No. 2002-08, Institut für Angewandte Physik, Hamburg (2002).
Çapoğlu  R. et al., “Accuracy of the born approximation in calculating the scattering coefficient of biological continuous random media,” Opt. Lett.. 34, (17 ), 2679 –2681 (2009). 0146-9592 CrossRef
Capoglu  I. R., “MATLAB script for generating Whittle-Matern correlated 3-D random media,” http://goo.gl/hU6En ( July 2012).
Wood  M. F. G., Guo  X., Vitkin  I. A., “Polarized light propagation in multiply scattering media exhibiting both linear birefringence and optical activity: Monte Carlo model and experimental methodology,” J. Biomed. Opt.. 12, (1 ), 014029  (2007). 1083-3668 CrossRef
Poole  L. R., Venable  D. D., Campbell  J. W., “Semianalytic Monte Carlo radiative transfer model for oceanographic lidar systems,” Appl. Opt. 20, (20 ), 3653 –3656 (1981). 0003-6935 CrossRef
Tinet  E., Avrillier  S., Tualle  J., “Fast semianalytical Monte Carlo simulation for time-resolved light propagation in turbid media,” J. Opt. Soc. Am. A. 13, (9 ), 1903 –1915 (1996). 0740-3232 CrossRef
Cheong  W., Prahl  S., Welch  A., “A review of the optical properties of biological tissues,” IEEE J. Quant. Electron.. 26, (12 ), 2166 –2185 (1990). 0018-9197 CrossRef
Mishchenko  M. I., Travis  L. D., Lacis  A. A., Multiple Scattering of Light by Particles: Radiative Transfer and Coherent Backscattering. ,  Cambridge University Press ,  Cambridge  (2006).
Amic  E., Luck  J. M., Nieuwenhuizen  T. M., “Multiple Rayleigh scattering of electromagnetic waves,” J. Phys. I. 7, (3 ), 445 –483 (1997). 1155-4304 CrossRef
Kuzmin  V. L., Meglinski  I. V., “Numerical simulation of coherent backscattering and temporal intensity correlations in random media,” Quant. Electron.. 36, (11 ), 990 –1002 (2006). 1063-7818 CrossRef
Everett  M. J. et al., “Birefringence characterization of biological tissue by use of optical coherence tomography,” Opt. Lett.. 23, (3 ), 228 –230 (1998). 0146-9592 CrossRef
Doronin  A., Meglinski  M., “Peer-to-peer Monte Carlo simulation of photon migration in topical applications of biomedical optics,” J. Biomed. Opt.. 17, (9 ), 090504  (2012). 1083-3668 CrossRef
Nadiarnykh  O. et al., “Alterations of the extracellular matrix in ovarian cancer studied by second harmonic generation imaging microscopy,” BMC Cancer. 10, (1 ), 94  (2010).CrossRef
© 2012 Society of Photo-Optical Instrumentation Engineers

Citation

Andrew J. Radosevich ; Jeremy D. Rogers ; İlker R. Çapoğlu ; Nikhil N. Mutyal ; Prabhakar Pradhan, et al.
"Open source software for electric field Monte Carlo simulation of coherent backscattering in biological media containing birefringence", J. Biomed. Opt. 17(11), 115001 (Nov 02, 2012). ; http://dx.doi.org/10.1117/1.JBO.17.11.115001


Figures

Graphic Jump LocationF5 :

Enhancement factor E versus klc for various polarization channels under the Whittle-Matérn model. (a) Trends for D=2. (b) Trends for D=3. (c) Trends for D=4.

Graphic Jump LocationF2 :

Comparison between the semi-analytical and traditional Monte Carlo method in a sample of Rayleigh scatterers simulated for the xx polarization channel. Function pxx(rs·μs*,z·μs*) for (a) the semi-analytical method and (b) the traditional method. (c) Functions pxx(rs·μs*) and pxx(zs·μs*) achieved by summing pxx(rs·μs*,z·μs*) over rows and columns, respectively. The symbols indicate the traditional technique while the solid line indicates the semi-analytical approach. (d) Computational efficiency measured by taking the ratio of the number of traditional photons over the number of semi-analytical photons needed to achieve the same simulation noise level. These simulations utilize the Whittle-Matérn model with D=3.

Graphic Jump LocationF3 :

MATLAB GUI for performing Monte Carlo simulations of CBS. (a) Specification of the general Monte Carlo parameters. (b) Specification of scattering model to be implemented. (c) Specification of the birefringence properties of the sample.

Graphic Jump LocationF1 :

Numerical validation of our treatment of scattering using function Mn(r). (a) Randomly generated 3-D medium under the Whittle-Matérn model for D=3. (b) Comparison of function |f| calculated numerically and analytically for the medium shown in panel a.

Graphic Jump LocationF4 :

Enhancement factor E versus size parameter ka for various polarization channels under Mie theory. (a) Trends for relative refractive index m=nsphere/nmedium corresponding to aqueous suspension of polystyrene microspheres. (b) Trends for m=1.1. (c) Trends for m=1.001.

Graphic Jump LocationF6 :

Effects of linear birefringence on the spatial reflectance profiles under linear polarized illumination and collection for a sample of Rayleigh scatterers. (a) The normalized intensity distributions pxx·pcxx and pxy·pcxy which are measurable using CBS. The dashed line indicates the shape of the incoherent pxy distribution which is not directly measurable using CBS. (b) Function pcxy. Due to the full reversibility of the xx polarization channel, pcxx=1 for all length-scales. In both panels, the arrow indicates increasing Δnb,max=0, 4×104, 7×104, and 1×103.

Graphic Jump LocationF7 :

Effects of linear birefringence on the spatial reflectance profiles under circularly polarized illumination and collection for a sample of Rayleigh scatterers. (a) The normalized intensity distributions p++·pc++ and p+·pc+ which are measurable using CBS. The dashed line indicates the shape of the incoherent p+ distribution which is not directly measurable using CBS. (b) Function pc+. Due to the full reversibility of the ++ polarization channel, pc++=1 for all length-scales. In both panels, the arrow indicates increasing Δnb,max=0, 4×104, 7×104, and 1×103.

Graphic Jump LocationF9 :

Changes in E as a function of Δnb,max for various polarization channels. (a) Trends of versus E. (b) Percent change in E from the absence of birefringence case versus Δnb,max.

Graphic Jump LocationF8 :

Alterations in the shape of p(xs,ys)·pc(xs,ys) due to birefringence for different polarization channels. To emphasize the shape of the various distributions, each image shows log10[p(xs,ys)·pc(xs,ys)] in the same intensity scale. (Rows) Each row corresponds to the polarization channel specified on the far right of the figure: top row shows xx polarization, second row shows ++ polarization, third row shows xy polarization, and the bottom row shows + polarization. (Columns) The left column shows the log10[p(xs,ys)·pc(xs,ys)] distributions for Δnb,max=0 and the middle column shows the distributions for Δnb,max=1×103. The right column shows the angular distribution found by converting to polar coordinates, summing over radius, and normalizing by the mean.

Tables

Table Grahic Jump Location
Table 1Values of E in various polarization channels for nonabsorbing Rayleigh scatterers with index matching at the boundary. The values given for Refs. 45 and 46 are converted into the normalization used in this paper before display.

References

Gomes  A. J. et al., “In vivo measurement of the shape of the tissue-refractive-index correlation function and its applicationto detection of colorectal field carcinogenesis,” J. Biomed. Opt.. 17, (4 ), 047005  (2012). 1083-3668 CrossRef
Wang  X., Wang  L. V., “Propagation of polarized light in birefringent turbid media: a Monte Carlo study,” J. Biomed. Opt.. 7, (3 ), 279 –290 (2002). 1083-3668 CrossRef
Radosevich  A. J. et al., “Measurement of the spatial backscattering impulse-response at short length scales with polarized enhanced backscattering,” Opt. Lett.. 36, (24 ), 4737 –4739 (2011). 0146-9592 CrossRef
Kuga  Y., Ishimaru  A., “Retroreflectance from a dense distribution of spherical particles,” J. Opt. Soc. Am. A. 1, (8 ), 831 –835 (1984). 0740-3232 CrossRef
Tsang  L., Ishimaru  A., “Backscattering enhancement of random discrete scatterers,” J. Opt. Soc. Am. A. 1, (8 ), 836 –839 (1984). 0740-3232 CrossRef
Wolf  P. E., Maret  G., “Weak localization and coherent backscattering of photons in disordered media,” Phys. Rev. Lett.. 55, (24 ), 2696 –2699 (1985). 0031-9007 CrossRef
Albada  M. P. V., Lagendijk  A., “Observation of weak localization of light in a random medium,” Phys. Rev. Lett.. 55, (24 ), 2692 –2695 (1985). 0031-9007 CrossRef
Lenke  R., Maret  G., Multiple Scattering of Light: Coherent Backscattering and Transmission, Scattering in Polymeric and Colloidal Systems. , Brown  W., Mortensen  K., Eds., pp. 1 –73,  Dover ,  London  (2000).
Dogariu  A., Uozumi  J., Asakura  T., “Enhancement factor in the light backscattered by fractal aggregated media,” Opt. Rev.. 3, (2 ), 71 –82 (1996). 1340-6000 CrossRef
Ishii  K., Iwai  T., Asakura  T., “Polarization properties of the enhanced backscattering of light from the fractal aggregate of particles,” Opt. Rev.. 4, (6 ), 643 –647 (1997). 1340-6000 CrossRef
Wiersma  D. S., van Albada  M. P., Lagendijk  A., “Coherent backscattering of light from amplifying random media,” Phys. Rev. Lett.. 75, (9 ), 1739 –1742 (1995). 0031-9007 CrossRef
Labeyrie  G. et al., “Coherent backscattering of light by cold atoms,” Phys. Rev. Lett.. 83, (25 ), 5266 –5269 (1999). 0031-9007 CrossRef
Vithana  H. K. M., Asfaw  L., Johnson  D. L., “Coherent backscattering of light in a nematic liquid crystal,” Phys. Rev. Lett.. 70, (23 ), 3561 –3564 (1993). 0031-9007 CrossRef
Sapienza  R. et al., “Anisotropic weak localization of light,” Phys. Rev. Lett.. 92, (3 ), 033903  (2004). 0031-9007 CrossRef
Yoo  K. M., Tang  G. C., Alfano  R. R., “Coherent backscattering of light from biological tissues,” Appl. Opt.. 29, (22 ), 3237 –3239 (1990). 0003-6935 CrossRef
Yoon  G., Roy  D. N. G., Straight  R. C., “Coherent backscattering in biological media: measurement and estimation of optical properties,” Appl. Opt.. 32, (4 ), 580 –585 (1993). 0003-6935 CrossRef
Radosevich  A. et al., “Polarized enhanced backscattering spectroscopy for characterization of biological tissues at subdiffusion length-scales,” IEEE J. Sel. Top. Quant. Electron.. 18, (4 ), 1313 –1325 (2012). 1077-260X CrossRef
Turzhitsky  V. et al., “Measurement of optical scattering properties with low-coherence enhanced backscattering spectroscopy,” J. Biomed. Opt.. 16, (6 ), 067007  (2011). 1083-3668 CrossRef
Roy  H. K. et al., “Association between rectal optical signatures and colonic neoplasia: potential applications for screening,” Cancer Res.. 69, , 4476 –83 (2009). 0008-5472 CrossRef
Turzhitsky  V. et al., “Investigating population risk factors of pancreatic cancer by evaluation of optical markers in the duodenal mucosa,” Dis. Markers. 25, (6 ), 313 –321 (2008). 0278-0240 
Akkermans  E., Wolf  P. E., Maynard  R., “Coherent backscattering of light by disordered media: Analysis of the peak line shape,” Phys. Rev. Lett.. 56, (14 ), 1471 –1474 (1986). 0031-9007 CrossRef
Ishii  K., Iwai  T., “Double-scattering approximation theory of enhanced backscatterings of light produced from aggregated particles,” Jpn. J. Appl. Phys. 1. 41, (Part 1, No. 8 ), 5155 –5159 (2002). 0021-4922 CrossRef
Eddowes  M. H., Mills  T. N., Delpy  D. T., “Monte Carlo simulations of coherent backscatter for identification of the optical coefficients of biological tissues in vivo,” Appl. Opt.. 34, (13 ), 2261 –2267 (1995). 0003-6935 CrossRef
Turzhitsky  V. et al., “A predictive model of backscattering at subdiffusion length scales,” Biomed. Opt. Express. 1, (3 ), 1034 –1046 (2010). 2156-7085 CrossRef
Martinez  A. S., Maynard  R., “Faraday effect and multiple scattering of light,” Phys. Rev. B. 50, (6 ), 3714 –3732 (1994). 1098-0121 CrossRef
Sawicki  J., Kastor  N., Xu  M., “Electric field Monte Carlo simulation of coherent backscattering of polarized light by a turbid medium containing mie scatterers,” Opt. Express. 16, (8 ), 5728 –5738 (2008). 1094-4087 CrossRef
Lenke  R., Maret  G., “Magnetic field effects on coherent backscattering of light,” Eur. Phys. J. B. 17, (1 ), 171 –185 (2000). 1434-6028 CrossRef
Lenke  R., Tweer  R., Maret  G., “Coherent backscattering of turbid samples containing large mie spheres,” J. Opt. A—Pure Appl. Opt.. 4, (3 ), 293  (2002). 1464-4258 CrossRef
Born  M., Wolf  E., Principles of Optics: Electromagnetic Theory of Propagation, Interference and Diffraction of Light. , 7th ed.,  Cambridge University Press ,  Cambridge  (1999).
Radosevich  Andrew J., Rogers  Jeremy D., Matlab functions, Biophotonics Laboratory at Northwestern University, http://biophotonics.bme.northwestern.edu/resources/index.html (30 October 2012).
Jones  R. C., “A new calculus for the treatment of optical systems,” J. Opt. Soc. Am.. 31, (7 ), 488 –493 (1941). 0030-3941 CrossRef
Ramella-Roman  J., Prahl  S., Jacques  S., “Three Monte Carlo programs of polarized light transport into scattering media: part I,” Opt. Exp.. 13, (12 ), 4420 –4438 ( Jun 2005). 1094-4087 CrossRef
van de Hulst  H., Light Scattering by Small Particles. ,  Dover ,  New York  (1981).
Jones  R. C., “A new calculus for the treatment of optical systems. VII. properties of the n-matrices,” J. Opt. Soc. Am.. 38, (8 ), 671 –683 (1948). 0030-3941 CrossRef
Bohren  C. F., Huffman  D. R., Absorption and Scattering of Light by Small Particles. ,  John Wiley and Sons ,  New York  (1983).
Guttorp  P., Gneiting  T., “Studies in the history of probability and statistics XLIX on the Matérn correlation family,” Biometrika. 93, (4 ), 989 –995 (2006). 0006-3444 CrossRef
Rogers  J. D., Çapoğlu  R., Backman  V., “Nonscalar elastic light scattering from continuous random media in the born approximation,” Opt. Lett.. 34, (12 ), 1891 –1893 (2009). 0146-9592 CrossRef
Mätzler  C., “Matlab functions for mie scattering and absorption,” Research Report No. 2002-08, Institut für Angewandte Physik, Hamburg (2002).
Çapoğlu  R. et al., “Accuracy of the born approximation in calculating the scattering coefficient of biological continuous random media,” Opt. Lett.. 34, (17 ), 2679 –2681 (2009). 0146-9592 CrossRef
Capoglu  I. R., “MATLAB script for generating Whittle-Matern correlated 3-D random media,” http://goo.gl/hU6En ( July 2012).
Wood  M. F. G., Guo  X., Vitkin  I. A., “Polarized light propagation in multiply scattering media exhibiting both linear birefringence and optical activity: Monte Carlo model and experimental methodology,” J. Biomed. Opt.. 12, (1 ), 014029  (2007). 1083-3668 CrossRef
Poole  L. R., Venable  D. D., Campbell  J. W., “Semianalytic Monte Carlo radiative transfer model for oceanographic lidar systems,” Appl. Opt. 20, (20 ), 3653 –3656 (1981). 0003-6935 CrossRef
Tinet  E., Avrillier  S., Tualle  J., “Fast semianalytical Monte Carlo simulation for time-resolved light propagation in turbid media,” J. Opt. Soc. Am. A. 13, (9 ), 1903 –1915 (1996). 0740-3232 CrossRef
Cheong  W., Prahl  S., Welch  A., “A review of the optical properties of biological tissues,” IEEE J. Quant. Electron.. 26, (12 ), 2166 –2185 (1990). 0018-9197 CrossRef
Mishchenko  M. I., Travis  L. D., Lacis  A. A., Multiple Scattering of Light by Particles: Radiative Transfer and Coherent Backscattering. ,  Cambridge University Press ,  Cambridge  (2006).
Amic  E., Luck  J. M., Nieuwenhuizen  T. M., “Multiple Rayleigh scattering of electromagnetic waves,” J. Phys. I. 7, (3 ), 445 –483 (1997). 1155-4304 CrossRef
Kuzmin  V. L., Meglinski  I. V., “Numerical simulation of coherent backscattering and temporal intensity correlations in random media,” Quant. Electron.. 36, (11 ), 990 –1002 (2006). 1063-7818 CrossRef
Everett  M. J. et al., “Birefringence characterization of biological tissue by use of optical coherence tomography,” Opt. Lett.. 23, (3 ), 228 –230 (1998). 0146-9592 CrossRef
Doronin  A., Meglinski  M., “Peer-to-peer Monte Carlo simulation of photon migration in topical applications of biomedical optics,” J. Biomed. Opt.. 17, (9 ), 090504  (2012). 1083-3668 CrossRef
Nadiarnykh  O. et al., “Alterations of the extracellular matrix in ovarian cancer studied by second harmonic generation imaging microscopy,” BMC Cancer. 10, (1 ), 94  (2010).CrossRef

Some tools below are only available to our subscribers or users with an online account.

Related Content

Customize your page view by dragging & repositioning the boxes below.