 Convolution for optical broadbeam responses in scattering media

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 responses. Photon transport methods can be directly used to compute broadbeam responses by distributing photons over the cross section of the beam. However, convolution can be used in certain cases to improve computational efficiency.
Contents
General convolution formulas
In order for convolution to be used to calculate a broadbeam response, a system must be time invariant, 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 broadbeam response can be calculated using the following 2D convolution formula:
Similar to 1D convolution, 2D convolution is commutative between G and S with a change of variables and :
Because the broadbeam 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 broadbeam response is more computationally advantageous.
Common beam profiles
Gaussian beam
For a Gaussian beam, the intensity profile is given by
Here, R denotes the radius of the beam, and S_{0} denotes the intensity at the center of the beam. S_{0} is related to the total power P_{0} by
Substituting Eq. 5 into Eq. 4, we obtain
where I_{0} is the zerothorder modified Bessel function.
Tophat beam
For a tophat beam of radius R, the source function becomes
where S_{0} denotes the intensity inside the beam. S_{0} is related to the total beam power P_{0} by
Substituting Eq. 8 into Eq. 4, we obtain
where
Errors in numerical evaluation
First interactions
First photontissue 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 tophat beam, we have
Truncation error
For a tophat beam, the upper integration limits may be bounded by r_{max}, such that r ≤ r_{max} − 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 tophat beam, we must ensure that r_{max} in photon transport methods is large enough that r ≤ r_{max} − 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 tophat beam of the same R and S_{0} have comparable convolution results. Therefore, r ≤ r_{max} − 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 transformation (FFT and IFFT) according to the convolution theorem. To calculate the optical broadbeam response, the impulse response of a pencil beam is convolved with the beam function. As shown by Equation 4, this is a 2D 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)^{2}b^{2}]. 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)^{2}b^{2}] 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 (n_{rel} = 1, μ_{a} = 0.1, μ_{s}=100, g = 0.9, 100,000 photons). Using this Matlab model, the fluence of a 3 × 3 × 3 cm^{3} region is recorded and the fluence distribution of a broadbeam response is plotted. Figure 1 and Figure 2 show the responses to a pencil beam and a 1cm tophat broadbeam, respectively. Direct convolution was used to calculate the broadbeam response in Figure 2. Figure 3 shows the broadbeam 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 tissue
 Monte Carlo method
 Monte Carlo method for photon transport
Links to other Monte Carlo resources
References
 L.H. Wang and H.I. Wu. Biomedical Optics: Principles and Imaging. Wiley 2007.
 L.H. Wang, S. L. Jacques, and L.Q. Zheng, "Monte Carlo modeling of photon transport in multilayered tissues," Computer Methods and Programs in Biomedicine 47, 131–146 (1995).
 L.H. Wang, S. L. Jacques, and L.Q. Zheng, "Convolution for responses to a finite diameter photon beam incident on multilayered tissues," Computer Methods and Programs in Biomedicine 54, 141–150 (1997). Download article.
Categories:
Wikimedia Foundation. 2010.