Convolution for optical broad-beam responses in scattering media
Encyclopedia
Photon transport theories, such as the Monte Carlo method, are commonly used to model light propagation in tissue
. The responses to a pencil beam
incident on a scattering medium are referred to as Green’s functions
or impulse response
s. Photon transport methods can be directly used to compute broad-beam responses by distributing photons over the cross section of the beam. However, convolution
can be used in certain cases to improve computational efficiency.
, linear
, and translation invariant
. Time invariance implies that a photon beam delayed by a given time produces a response shifted by the same delay. Linearity indicates that a given response will increase by the same amount if the input is scaled and obeys the property of superposition
. Translational invariance means that if a beam is shifted to a new location on the tissue surface, its response is also shifted in the same direction by the same distance. Here, only spatial convolution is considered.
Responses from photon transport methods can be physical quantities such as absorption
, fluence
, reflectance, or transmittance
. Given a specific physical quantity, G(x,y,z), from a pencil beam in Cartesian space and a collimated light source with beam profile S(x,y), a broad-beam response can be calculated using the following 2-D convolution formula:
Similar to 1-D convolution, 2-D convolution is commutative between G and S with a change of variables
and 
:
Because the broad-beam response
has cylindrical symmetry, its convolution integrals can be rewritten as:
where
. Because the inner integration of Equation 4 is independent of z, it only needs to be calculated once for all depths. Thus this form of the broad-beam response is more computationally advantageous.
, the intensity profile is given by
Here, R denotes the
radius of the beam, and S0 denotes the intensity at the center of the beam. S0 is related to the total power P0 by
Substituting Eq. 5 into Eq. 4, we obtain
where I0 is the zeroth-order modified Bessel function.
of radius R, the source function becomes
where S0 denotes the intensity inside the beam. S0 is related to the total beam power P0 by
Substituting Eq. 8 into Eq. 4, we obtain
where
where the first term results from the first interactions and the second, from subsequent interactions.
For a Gaussian beam, we have
For a top-hat beam, we have
For a Gaussian beam, no simple upper integration limits exist because it theoretically extends to infinity. At r >> R, a Gaussian beam and a top-hat beam of the same R and S0 have comparable convolution results. Therefore, r ≤ rmax − R can be used approximately for Gaussian beams as well.
ation (FFT and IFFT) according to the convolution theorem
. To calculate the optical broad-beam response, the impulse response of a pencil beam is convolved with the beam function. As shown by Equation 4, this is a 2-D convolution. To calculate the response of a light beam on a plane perpendicular to the z axis, the beam function (represented by a b × b matrix) is convolved with the impulse response on that plane (represented by an a × a matrix). Normally a is greater than b. The calculation efficiency of these two methods depends largely on b, the size of the light beam.
In direct convolution, the solution matrix is of the size (a + b − 1) × (a + b − 1). The calculation of each of these elements (except those near boundaries) includes b × b multiplications and b × b − 1 additions, so the time complexity
is O
[(a + b)2b2]. Using the FFT method, the major steps are the FFT and IFFT of (a + b − 1) × (a + b − 1) matrices, so the time complexity is O[(a + b)2 log(a + b)]. Comparing O[(a + b)2b2] and O[(a + b)2 log(a + b)], it is apparent that direct convolution will be faster if b is much smaller than a, but the FFT method will be faster if b is relatively large.
Monte Carlo method for photon transport
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-tissue interaction...
. The responses to a pencil beam
Pencil beam
In optics, a pencil or pencil of rays is a geometric construct used to describe a beam or portion of a beam of electromagnetic radiation or charged particles, typically in the form of a narrow cone or cylinder....
incident on a scattering medium are referred to as Green’s functions
Green's function
In mathematics, a Green's function is a type of function used to solve inhomogeneous differential equations subject to specific initial conditions or boundary conditions...
or impulse response
Impulse response
In signal processing, the impulse response, or impulse response function , of a dynamic system is its output when presented with a brief input signal, called an impulse. More generally, an impulse response refers to the reaction of any dynamic system in response to some external change...
s. Photon transport methods can be directly used to compute broad-beam responses by distributing photons over the cross section of the beam. However, convolution
Convolution
In mathematics and, in particular, functional analysis, convolution is a mathematical operation on two functions f and g, producing a third function that is typically viewed as a modified version of one of the original functions. Convolution is similar to cross-correlation...
can be used in certain cases to improve computational efficiency.
General convolution formulas
In order for convolution to be used to calculate a broad-beam response, a system must be time invariantTime-invariant system
A time-invariant system is one whose output does not depend explicitly on time.This property can be satisfied if the transfer function of the system is not a function of time except expressed by the input and output....
, linear
Linear
In mathematics, a linear map or function f is a function which satisfies the following two properties:* Additivity : f = f + f...
, and translation invariant
Translational symmetry
In geometry, a translation "slides" an object by a a: Ta = p + a.In physics and mathematics, continuous translational symmetry is the invariance of a system of equations under any translation...
. Time invariance implies that a photon beam delayed by a given time produces a response shifted by the same delay. Linearity indicates that a given response will increase by the same amount if the input is scaled and obeys the property of superposition
Superposition principle
In physics and systems theory, the superposition principle , also known as superposition property, states that, for all linear systems, the net response at a given place and time caused by two or more stimuli is the sum of the responses which would have been caused by each stimulus individually...
. Translational invariance means that if a beam is shifted to a new location on the tissue surface, its response is also shifted in the same direction by the same distance. Here, only spatial convolution is considered.
Responses from photon transport methods can be physical quantities such as absorption
Absorption (electromagnetic radiation)
In physics, absorption of electromagnetic radiation is the way by which the energy of a photon is taken up by matter, typically the electrons of an atom. Thus, the electromagnetic energy is transformed to other forms of energy for example, to heat. The absorption of light during wave propagation is...
, fluence
Fluence
In physics, fluence is the flux integrated over time. For particles, it is defined as the total number of particles that intersect a unit area in a specific time interval of interest, and has units of m–2...
, reflectance, or transmittance
Transmittance
In optics and spectroscopy, transmittance is the fraction of incident light at a specified wavelength that passes through a sample. A related term is absorptance, or absorption factor, which is the fraction of radiation absorbed by a sample at a specified wavelength...
. Given a specific physical quantity, G(x,y,z), from a pencil beam in Cartesian space and a collimated light source with beam profile S(x,y), a broad-beam response can be calculated using the following 2-D convolution formula:
Similar to 1-D convolution, 2-D convolution is commutative between G and S with a change of variables
and :Because the broad-beam response
has cylindrical symmetry, its convolution integrals can be rewritten as:where
. Because the inner integration of Equation 4 is independent of z, it only needs to be calculated once for all depths. Thus this form of the broad-beam response is more computationally advantageous.Gaussian beam
For a Gaussian beamGaussian beam
In optics, a Gaussian beam is a beam of electromagnetic radiation whose transverse electric field and intensity distributions are well approximated by Gaussian functions. Many lasers emit beams that approximate a Gaussian profile, in which case the laser is said to be operating on the fundamental...
, the intensity profile is given by
Here, R denotes the
radius of the beam, and S0 denotes the intensity at the center of the beam. S0 is related to the total power P0 bySubstituting Eq. 5 into Eq. 4, we obtain
where I0 is the zeroth-order modified Bessel function.
Top-hat beam
For a top-hat beamTophat beam
In optics, a tophat beam is a laser beam with a near-uniform fluence within a circular disk. It is typically formed by diffractive optical elements from a Gaussian beam. Tophat beams are often used in industry, for example for laser drilling of holes in printed circuit boards...
of radius R, the source function becomes
where S0 denotes the intensity inside the beam. S0 is related to the total beam power P0 by
Substituting Eq. 8 into Eq. 4, we obtain
where
First interactions
First photon-tissue interactions always occur on the z axis and hence contribute to the specific absorption or related physical quantities as a delta function. Errors will result if absorption due to the first interactions is not recorded separately from absorption due to subsequent interactions. The total impulse response can be expressed in two parts:where the first term results from the first interactions and the second, from subsequent interactions.
For a Gaussian beam, we have
For a top-hat beam, we have
Truncation error
For a top-hat beam, the upper integration limits may be bounded by rmax, such that r ≤ rmax − R. Thus, the limited grid coverage in the r direction does not affect the convolution. To convolve reliably for physical quantities at r in response to a top-hat beam, we must ensure that rmax in photon transport methods is large enough that r ≤ rmax − R holds.For a Gaussian beam, no simple upper integration limits exist because it theoretically extends to infinity. At r >> R, a Gaussian beam and a top-hat beam of the same R and S0 have comparable convolution results. Therefore, r ≤ rmax − R can be used approximately for Gaussian beams as well.
Implementation of convolution
There are two common methods used to implement discrete convolution: the definition of convolution and fast Fourier transformFast Fourier transform
A fast Fourier transform is an efficient algorithm to compute the discrete Fourier transform and its inverse. "The FFT has been called the most important numerical algorithm of our lifetime ." There are many distinct FFT algorithms involving a wide range of mathematics, from simple...
ation (FFT and IFFT) according to the convolution theorem
Convolution theorem
In mathematics, the convolution theorem states that under suitableconditions the Fourier transform of a convolution is the pointwise product of Fourier transforms. In other words, convolution in one domain equals point-wise multiplication in the other domain...
. To calculate the optical broad-beam response, the impulse response of a pencil beam is convolved with the beam function. As shown by Equation 4, this is a 2-D convolution. To calculate the response of a light beam on a plane perpendicular to the z axis, the beam function (represented by a b × b matrix) is convolved with the impulse response on that plane (represented by an a × a matrix). Normally a is greater than b. The calculation efficiency of these two methods depends largely on b, the size of the light beam.
In direct convolution, the solution matrix is of the size (a + b − 1) × (a + b − 1). The calculation of each of these elements (except those near boundaries) includes b × b multiplications and b × b − 1 additions, so the time complexity
Analysis of algorithms
To analyze an algorithm is to determine the amount of resources necessary to execute it. Most algorithms are designed to work with inputs of arbitrary length...
is O
Big O notation
In mathematics, big O notation is used to describe the limiting behavior of a function when the argument tends towards a particular value or infinity, usually in terms of simpler functions. It is a member of a larger family of notations that is called Landau notation, Bachmann-Landau notation, or...
[(a + b)2b2]. Using the FFT method, the major steps are the FFT and IFFT of (a + b − 1) × (a + b − 1) matrices, so the time complexity is O[(a + b)2 log(a + b)]. Comparing O[(a + b)2b2] and O[(a + b)2 log(a + b)], it is apparent that direct convolution will be faster if b is much smaller than a, but the FFT method will be faster if b is relatively large.
Computational examples
The fate of photons can be modeled using a Matlab implementation of the Monte Carlo method (nrel = 1, μa = 0.1, μs=100, g = 0.9, 100,000 photons). Using this Matlab model, the fluence of a 3 × 3 × 3 cm3 region is recorded and the fluence distribution of a broad-beam response is plotted. Figure 1 and Figure 2 show the responses to a pencil beam and a 1-cm top-hat broad-beam, respectively. Direct convolution was used to calculate the broad-beam response in Figure 2. Figure 3 shows the broad-beam response calculated using the FFT method. When the diameter of the light beam is 0.2 cm, direct convolution costs 1.93 seconds, and the FFT method costs 7.35 seconds. When the diameter of the light beam is 2 cm, direct convolution costs 90.1 seconds, and FFT method costs 16.8 seconds. Of course, the absolute computation time depends on the processing speed of the computer being used. These two comparisons were made on the same computer. Although the computation times differ, the plots in Figures 2 and 3 are indistinguishable.See also
- Radiative transfer equation and diffusion theory for photon transport in biological tissueRadiative transfer equation and diffusion theory for photon transport in biological tissuePhoton transport in biological tissue can be equivalently modeled numerically with Monte Carlo simulations or analytically by the radiative transfer equation . However, the RTE is difficult to solve without introducing approximations. A common approximation summarized here is the diffusion...
- Monte Carlo methodMonte Carlo methodMonte Carlo methods are a class of computational algorithms that rely on repeated random sampling to compute their results. Monte Carlo methods are often used in computer simulations of physical and mathematical systems...
- Monte Carlo method for photon transportMonte Carlo method for photon transportModeling 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-tissue interaction...