Monte Carlo method for photon transport

Last updated

Modeling photon propagation with Monte Carlo methods is a flexible yet rigorous approach to simulate photon transport. In the method, local rules of photon transport are expressed as probability distributions which describe the step size of photon movement between sites of photon-matter interaction and the angles of deflection in a photon's trajectory when a scattering event occurs. This is equivalent to modeling photon transport analytically by the radiative transfer equation (RTE), which describes the motion of photons using a differential equation. However, closed-form solutions of the RTE are often not possible; for some geometries, the diffusion approximation can be used to simplify the RTE, although this, in turn, introduces many inaccuracies, especially near sources and boundaries. In contrast, Monte Carlo simulations can be made arbitrarily accurate by increasing the number of photons traced. For example, see the movie, where a Monte Carlo simulation of a pencil beam incident on a semi-infinite medium models both the initial ballistic photon flow and the later diffuse propagation.

Contents

The Monte Carlo method is necessarily statistical and therefore requires significant computation time to achieve precision. In addition Monte Carlo simulations can keep track of multiple physical quantities simultaneously, with any desired spatial and temporal resolution. This flexibility makes Monte Carlo modeling a powerful tool. Thus, while computationally inefficient, Monte Carlo methods are often considered the standard for simulated measurements of photon transport for many biomedical applications.

Monte Carlo simulation of a pencil beam incident on a semi-infinite scattering medium. MonteCarloSemiInf.gif
Monte Carlo simulation of a pencil beam incident on a semi-infinite scattering medium.

Biomedical applications of Monte Carlo methods

Biomedical imaging

The optical properties of biological tissue offer an approach to biomedical imaging. There are many endogenous contrasts, including absorption from blood and melanin and scattering from nerve cells and cancer cell nuclei. In addition, fluorescent probes can be targeted to many different tissues. Microscopy techniques (including confocal, two-photon, and optical coherence tomography) have the ability to image these properties with high spatial resolution, but, since they rely on ballistic photons, their depth penetration is limited to a few millimeters. Imaging deeper into tissues, where photons have been multiply scattered, requires a deeper understanding of the statistical behavior of large numbers of photons in such an environment. Monte Carlo methods provide a flexible framework that has been used by different techniques to reconstruct optical properties deep within tissue. A brief introduction to a few of these techniques is presented here.

Radiation therapy

The goal of radiation therapy is to deliver energy, generally in the form of ionizing radiation, to cancerous tissue while sparing the surrounding normal tissue. Monte Carlo modeling is commonly employed in radiation therapy to determine the peripheral dose the patient will experience due to scattering, both from the patient tissue as well as scattering from collimation upstream in the linear accelerator.

Photodynamic therapy

In Photodynamic therapy (PDT) light is used to activate chemotherapy agents. Due to the nature of PDT, it is useful to use Monte Carlo methods for modeling scattering and absorption in the tissue in order to ensure appropriate levels of light are delivered to activate chemotherapy agents.

Implementation of photon transport in a scattering medium

Presented here is a model of a photon Monte Carlo method in a homogeneous infinite medium. The model is easily extended for multi-layered media, however. For an inhomogeneous medium, boundaries must be considered. In addition for a semi-infinite medium (in which photons are considered lost if they exit the top boundary), special consideration must be taken. For more information, please visit the links at the bottom of the page. We will solve the problem using an infinitely small point source (represented analytically as a Dirac delta function in space and time). Responses to arbitrary source geometries can be constructed using the method of Green's functions (or convolution, if enough spatial symmetry exists). The required parameters are the absorption coefficient, the scattering coefficient, and the scattering phase function. (If boundaries are considered the index of refraction for each medium must also be provided.) Time-resolved responses are found by keeping track of the total elapsed time of the photon's flight using the optical path length. Responses to sources with arbitrary time profiles can then be modeled through convolution in time.

In our simplified model we use the following variance reduction technique to reduce computational time. Instead of propagating photons individually, we create a photon packet with a specific weight (generally initialized as unity). As the photon interacts in the turbid medium, it will deposit weight due to absorption and the remaining weight will be scattered to other parts of the medium. Any number of variables can be logged along the way, depending on the interest of a particular application. Each photon packet will repeatedly undergo the following numbered steps until it is either terminated, reflected, or transmitted. The process is diagrammed in the schematic to the right. Any number of photon packets can be launched and modeled, until the resulting simulated measurements have the desired signal-to-noise ratio. Note that as Monte Carlo modeling is a statistical process involving random numbers, we will be using the variable ξ throughout as a pseudo-random number for many calculations.

Schematic for modeling photon flow in an infinite scattering and absorbing medium with Monte Carlo simulations. MonteCarlo.png
Schematic for modeling photon flow in an infinite scattering and absorbing medium with Monte Carlo simulations.

Step 1: Launching a photon packet

In our model, we are ignoring initial specular reflectance associated with entering a medium that is not refractive index matched. With this in mind, we simply need to set the initial position of the photon packet as well as the initial direction. It is convenient to use a global coordinate system. We will use three Cartesian coordinates to determine position, along with three direction cosines to determine the direction of propagation. The initial start conditions will vary based on application, however for a pencil beam initialized at the origin, we can set the initial position and direction cosines as follows (isotropic sources can easily be modeled by randomizing the initial direction of each packet):

Step 2: Step size selection and photon packet movement

The step size, s, is the distance the photon packet travels between interaction sites. There are a variety of methods for step size selection. Below is a basic form of photon step size selection (derived using the inverse distribution method and the Beer–Lambert law) from which we use for our homogeneous model:

where is a random number and is the total interaction coefficient (i.e., the sum of the absorption and scattering coefficients).

Once a step size is selected, the photon packet is propagated by a distance s in a direction defined by the direction cosines. This is easily accomplished by simply updating the coordinates as follows:

Step 3: Absorption and scattering

A portion of the photon weight is absorbed at each interaction site. This fraction of the weight is determined as follows:

where is the absorption coefficient.

The weight fraction can then be recorded in an array if an absorption distribution is of interest for the particular study. The weight of the photon packet must then be updated as follows:

Following absorption, the photon packet is scattered. The weighted average of the cosine of the photon scattering angle is known as scattering anisotropy (g), which has a value between 1 and 1. If the optical anisotropy is 0, this generally indicates that the scattering is isotropic. If g approaches a value of 1 this indicates that the scattering is primarily in the forward direction. In order to determine the new direction of the photon packet (and hence the photon direction cosines), we need to know the scattering phase function. Often the Henyey-Greenstein phase function is used. Then the scattering angle, θ, is determined using the following formula.

And, the polar angle φ is generally assumed to be uniformly distributed between 0 and . Based on this assumption, we can set:

Based on these angles and the original direction cosines, we can find a new set of direction cosines. The new propagation direction can be represented in the global coordinate system as follows:

For a special case

use

or

use

C-code:

/*********************** Indicatrix ********************* *New direction cosines after scattering by angle theta, fi. * mux new=(sin(theta)*(mux*muz*cos(fi)-muy*sin(fi)))/sqrt(1-muz^2)+mux*cos(theta) * muy new=(sin(theta)*(muy*muz*cos(fi)+mux*sin(fi)))/sqrt(1-muz^2)+muy*cos(theta) * muz new= - sqrt(1-muz^2)*sin(theta)*cos(fi)+muz*cos(theta) *--------------------------------------------------------- *Input: * muxs,muys,muzs - direction cosine before collision * mutheta, fi - cosine of polar angle and the azimuthal angle *--------------------------------------------------------- *Output: *  muxd,muyd,muzd - direction cosine after collision *--------------------------------------------------------- */ void Indicatrix (double muxs, double muys, double muzs, double mutheta, double fi, double *muxd, double *muyd, double *muzd) {  double costheta = mutheta;  double sintheta = sqrt(1.0-costheta*costheta); // sin(theta)  double sinfi = sin(fi);  double cosfi = cos(fi);  if (muzs == 1.0) {    *muxd = sintheta*cosfi;    *muyd = sintheta*sinfi;    *muzd = costheta;  } elseif (muzs == -1.0) {    *muxd = sintheta*cosfi;    *muyd = -sintheta*sinfi;    *muzd = -costheta;  } else {    double denom = sqrt(1.0-muzs*muzs);    double muzcosfi = muzs*cosfi;    *muxd = sintheta*(muxs*muzcosfi-muys*sinfi)/denom + muxs*costheta;    *muyd = sintheta*(muys*muzcosfi+muxs*sinfi)/denom + muys*costheta;    *muzd = -denom*sintheta*cosfi + muzs*costheta;  } }

Step 4: Photon termination

If a photon packet has experienced many interactions, for most applications the weight left in the packet is of little consequence. As a result, it is necessary to determine a means for terminating photon packets of sufficiently small weight. A simple method would use a threshold, and if the weight of the photon packet is below the threshold, the packet is considered dead. The aforementioned method is limited as it does not conserve energy. To keep total energy constant, a Russian roulette technique is often employed for photons below a certain weight threshold. This technique uses a roulette constant m to determine whether or not the photon will survive. The photon packet has one chance in m to survive, in which case it will be given a new weight of mW where W is the initial weight (this new weight, on average, conserves energy). All other times, the photon weight is set to 0 and the photon is terminated. This is expressed mathematically below:

Graphics Processing Units (GPU) and fast Monte Carlo simulations of photon transport

Monte Carlo simulation of photon migration in turbid media is a highly parallelizable problem, where a large number of photons are propagated independently, but according to identical rules and different random number sequences. The parallel nature of this special type of Monte Carlo simulation renders it highly suitable for execution on a graphics processing unit (GPU). The release of programmable GPUs started such a development, and since 2008 there have been a few reports on the use of GPU for high-speed Monte Carlo simulation of photon migration. [1] [2] [3] [4]

This basic approach can itself be parallelized by using multiple GPUs linked together. One example is the "GPU Cluster MCML," which can be downloaded from the authors' website (Monte Carlo Simulation of Light Transport in Multi-layered Turbid Media Based on GPU Clusters): http://bmp.hust.edu.cn/GPU_Cluster/GPU_Cluster_MCML.HTM

See also

Related Research Articles

In physics, the cross section is a measure of the probability that a specific process will take place in a collision of two particles. For example, the Rutherford cross-section is a measure of probability that an alpha particle will be deflected by a given angle during an interaction with an atomic nucleus. Cross section is typically denoted σ (sigma) and is expressed in units of area, more specifically in barns. In a way, it can be thought of as the size of the object that the excitation must hit in order for the process to occur, but more exactly, it is a parameter of a stochastic process.

<span class="mw-page-title-main">Euler's formula</span> Complex exponential in terms of sine and cosine

Euler's formula, named after Leonhard Euler, is a mathematical formula in complex analysis that establishes the fundamental relationship between the trigonometric functions and the complex exponential function. Euler's formula states that, for any real number x, one has where e is the base of the natural logarithm, i is the imaginary unit, and cos and sin are the trigonometric functions cosine and sine respectively. This complex exponential function is sometimes denoted cis x. The formula is still valid if x is a complex number, and is also called Euler's formula in this more general case.

<span class="mw-page-title-main">Electroweak interaction</span> Unified description of electromagnetism and the weak interaction

In particle physics, the electroweak interaction or electroweak force is the unified description of two of the four known fundamental interactions of nature: electromagnetism (electromagnetic interaction) and the weak interaction. Although these two forces appear very different at everyday low energies, the theory models them as two different aspects of the same force. Above the unification energy, on the order of 246 GeV, they would merge into a single force. Thus, if the temperature is high enough – approximately 1015 K – then the electromagnetic force and weak force merge into a combined electroweak force. During the quark epoch (shortly after the Big Bang), the electroweak force split into the electromagnetic and weak force. It is thought that the required temperature of 1015 K has not been seen widely throughout the universe since before the quark epoch, and currently the highest human-made temperature in thermal equilibrium is around 5.5×1012 K (from the Large Hadron Collider).

<span class="mw-page-title-main">Spherical coordinate system</span> Coordinates comprising a distance and two angles

In mathematics, a spherical coordinate system is a coordinate system for three-dimensional space where the position of a given point in space is specified by three real numbers: the radial distancer along the radial line connecting the point to the fixed point of origin; the polar angleθ between the radial line and a polar axis; and the azimuthal angleφ as the angle of rotation of the radial line around the polar axis. (See graphic re the "physics convention".) Once the radius is fixed, the three coordinates (r, θ, φ), known as a 3-tuple, provide a coordinate system on a sphere, typically called the spherical polar coordinates.

<span class="mw-page-title-main">Navier–Stokes equations</span> Equations describing the motion of viscous fluid substances

The Navier–Stokes equations are partial differential equations which describe the motion of viscous fluid substances. They were named after French engineer and physicist Claude-Louis Navier and the Irish physicist and mathematician George Gabriel Stokes. They were developed over several decades of progressively building the theories, from 1822 (Navier) to 1842–1850 (Stokes).

<span class="mw-page-title-main">Ellipsoid</span> Quadric surface that looks like a deformed sphere

An ellipsoid is a surface that can be obtained from a sphere by deforming it by means of directional scalings, or more generally, of an affine transformation.

<span class="mw-page-title-main">Unit vector</span> Vector of length one

In mathematics, a unit vector in a normed vector space is a vector of length 1. A unit vector is often denoted by a lowercase letter with a circumflex, or "hat", as in .

In probability theory, the Borel–Kolmogorov paradox is a paradox relating to conditional probability with respect to an event of probability zero. It is named after Émile Borel and Andrey Kolmogorov.

In rotordynamics, the rigid rotor is a mechanical model of rotating systems. An arbitrary rigid rotor is a 3-dimensional rigid object, such as a top. To orient such an object in space requires three angles, known as Euler angles. A special rigid rotor is the linear rotor requiring only two angles to describe, for example of a diatomic molecule. More general molecules are 3-dimensional, such as water, ammonia, or methane.

<span class="mw-page-title-main">Superquadrics</span> Family of geometric shapes

In mathematics, the superquadrics or super-quadrics are a family of geometric shapes defined by formulas that resemble those of ellipsoids and other quadrics, except that the squaring operations are replaced by arbitrary powers. They can be seen as the three-dimensional relatives of the superellipses. The term may refer to the solid object or to its surface, depending on the context. The equations below specify the surface; the solid is specified by replacing the equality signs by less-than-or-equal signs.

In Bayesian statistics, the Jeffreys prior is a non-informative prior distribution for a parameter space. Named after Sir Harold Jeffreys, its density function is proportional to the square root of the determinant of the Fisher information matrix:

<span class="mw-page-title-main">Viviani's curve</span> Figure-eight shaped curve on a sphere

In mathematics, Viviani's curve, also known as Viviani's window, is a figure eight shaped space curve named after the Italian mathematician Vincenzo Viviani. It is the intersection of a sphere with a cylinder that is tangent to the sphere and passes through two poles of the sphere. Before Viviani this curve was studied by Simon de La Loubère and Gilles de Roberval.

<span class="mw-page-title-main">Hopf bifurcation</span> Critical point where a periodic solution arises

In the mathematical theory of bifurcations, a Hopfbifurcation is a critical point where, as a parameter changes, a system's stability switches and a periodic solution arises. More accurately, it is a local bifurcation in which a fixed point of a dynamical system loses stability, as a pair of complex conjugate eigenvalues—of the linearization around the fixed point—crosses the complex plane imaginary axis as a parameter crosses a threshold value. Under reasonably generic assumptions about the dynamical system, the fixed point becomes a small-amplitude limit cycle as the parameter changes.

<span class="mw-page-title-main">Prolate spheroidal coordinates</span>

Prolate spheroidal coordinates are a three-dimensional orthogonal coordinate system that results from rotating the two-dimensional elliptic coordinate system about the focal axis of the ellipse, i.e., the symmetry axis on which the foci are located. Rotation about the other axis produces oblate spheroidal coordinates. Prolate spheroidal coordinates can also be considered as a limiting case of ellipsoidal coordinates in which the two smallest principal axes are equal in length.

In physics and mathematics, the solid harmonics are solutions of the Laplace equation in spherical polar coordinates, assumed to be (smooth) functions . There are two kinds: the regular solid harmonics, which are well-defined at the origin and the irregular solid harmonics, which are singular at the origin. Both sets of functions play an important role in potential theory, and are obtained by rescaling spherical harmonics appropriately:

<span class="mw-page-title-main">Radiative transfer equation and diffusion theory for photon transport in biological tissue</span>

Photon transport in biological tissue can be equivalently modeled numerically with Monte Carlo simulations or analytically by the radiative transfer equation (RTE). However, the RTE is difficult to solve without introducing approximations. A common approximation summarized here is the diffusion approximation. Overall, solutions to the diffusion equation for photon transport are more computationally efficient, but less accurate than Monte Carlo simulations.

In mathematics, vector spherical harmonics (VSH) are an extension of the scalar spherical harmonics for use with vector fields. The components of the VSH are complex-valued functions expressed in the spherical coordinate basis vectors.

In fluid dynamics, the Oseen equations describe the flow of a viscous and incompressible fluid at small Reynolds numbers, as formulated by Carl Wilhelm Oseen in 1910. Oseen flow is an improved description of these flows, as compared to Stokes flow, with the (partial) inclusion of convective acceleration.

References

Inline references