The Redistribution Function

Scott Prahl

January 2021

version 2

import numpy as np
import matplotlib.pyplot as plt
import iadpython

Redistribution function.

The single scattering phase function \(p(\nu)\) for a tissue determines the amount of light scattered at an angle \(\nu=\cos\theta\) from the direction of incidence. The subtended angle \(\nu\) is the dot product of the unit vectors \(\hat{\bf s}_i\) and \(\hat{\bf s}_j\)

\[\nu=\hat{\bf s}_i\cdot\hat{\bf s}_j=\nu_i\nu_j+\sqrt{1-\nu_i^2}\sqrt{1-\nu_j^2}\cos\phi\]

where \(\hat{\bf s}_i\) is the incident and \(\hat{\bf s}_j\) is the scattered light directions

The redistribution function \({\bf h}_{ij}\) determines the fraction of light scattered from an incidence cone with angle \(\nu_i\) into a cone with angle \(\nu_j\). The redistribution function is calculated by averaging the phase function over all possible azimuthal angles for fixed angles \(\nu_i\) and \(\nu_j\),

\[h(\nu_i,\nu_j) = {1\over2\pi} \int_0^{2\pi} p(\nu_i\nu_j+\sqrt{1-\nu_i^2}\sqrt{1-\nu_j^2}\cos\phi)\,d\phi\]

Note that the angles \(\nu_i\) and \(\nu_j\) may also be negative (light travelling in the opposite direction). The full redistribution matrix may be expressed in terms a \(2\times2\) matrix of |n|\(\times\)|n| matrices

\[\mathbf{h}=\left[\matrix{\mathbf{h}^{--}&\mathbf{h}^{-+}\cr \mathbf{h}^{+-}&\mathbf{h}^{++}\cr} \right]\]

The first plus or minus sign indicates the sign in front of the incident angle and the second is the sign of the direction of the scattered light.

When the cosine of the angle of incidence or exitance is unity (\(\nu_i=1\) or \(\nu_j=1\)), then the redistribution function \(h(1,\nu_j)\) is equivalent to the phase function \(p(\nu_j)\). In the case of isotropic scattering, the redistribution function is a constant

\[h(\nu_i,\nu_j) = p(\nu) = {1\over4\pi}.\]

Other phase functions require numerical integration of the phase function. If the phase function is highly anisotropic, then the integration over the azimuthal angle is particularly difficult and care must be taken to ensure that the integration is accurate. This is important because errors in the redistribution function enter directly into the reflection and transmission matrices for thin layers. Any errors will be doubled with each successive addition of layers and small errors will rapidly increase.

Redistribution Matrix using Legendre Polynomials

One way to calculate the redistribution function is the \(\delta\)-\(M\) method of Wiscombe. This method works especially well for highly anisotropic phase functions. The number of quadrature points is specified by \(M\). The \(\delta\)-\(M\) method approximates the true phase function by a phase function consisting of a Dirac delta function and \(M-1\) Legendre polynomials

\[p^*(\nu)= 2 g^M\delta(1-\nu) + (1-g^M) \sum_{k=0}^{M-1} (2k+1)\chi_k^* P_k(\nu)\]

where

\[\chi_k^*={\chi_k-g^M\over 1-g^M} \qquad\mbox{and}\qquad \chi_k = {1\over2}\int_0^1 p(\nu) P_k(\nu) \,d\nu\]

When the \(\delta\)-\(M\) method substitutes \(p^*(\nu)\rightarrow p(\nu)\), then both the albedo and optical thickness must also be changed, \(a^*\rightarrow a\) and \(\tau^*\rightarrow\tau\). This approximation is analogous to the similarity transformation often used to improve the diffusion approximation by moving a part (\(g^M\)) of the scattered light into the unscattered component. The new optical thickness and albedo are

\[\tau^*=(1-ag^M)\tau \qquad\mbox{and}\qquad a^* = a {1-g^M\over1-ag^M}\]

This is equivalent transforming the scattering coefficient as \(\mu_s^* = \mu_s(1-g^M)\). The redistribution function can now be written as

\[h^*(\nu_i,\nu_j) = \sum_{k=0}^{M-1} (2k+1)\chi_k^* P_k(\nu_i)P_k(\nu_j)\]

For the special case of a Henyey-Greenstein phase function,

\[\chi_k^*={g^k-g^M\over1-g^M}.\]

The current implementation is somewhat inefficient, but it works.

s = iadpython.Sample(g=0.9, quad_pts=4)

hp, hm = iadpython.hg_legendre(s)

print(hp)
print()
print(hm)

[[ 1.57558022  1.43478008  0.6702228  -0.09855935]
 [ 1.43478008  1.78554993  1.42038122  0.65844447]
 [ 0.6702228   1.42038122  2.73731157  3.69901785]
 [-0.09855935  0.65844447  3.69901785  6.84908404]]

[[ 1.49114428  1.10817646  0.3889397  -0.08632955]
 [ 1.10817646  0.49081177  0.10073799  0.22945985]
 [ 0.3889397   0.10073799  0.09249519  0.22802648]
 [-0.08632955  0.22945985  0.22802648 -0.37394591]]

Redistribution Matrix using Elliptic Functions

For Henyey-Greenstein scattering, the redistribution function can be expressed in terms of the complete elliptic integral of the second kind \(E(x)\)

\[h(\nu_i,\nu_j) = {2\over\pi}{1-g^2\over (\alpha-\gamma)\sqrt{\alpha+\gamma} } \,E\left(\sqrt{2\gamma\over\alpha+\gamma}\,\right)\]

where \(g\) is the average cosine of the Henyey-Greenstein phase function and

\[\alpha=1+g^2-2 g \nu_i \nu_j \qquad\mbox{and}\qquad \gamma=2 g \sqrt{1-\nu_i^2}\sqrt{1-\nu_j^2}\]

The function \(E(x)\) may be calculated using scipy.special.ellipe().

The drawback to this approach is that the \(\delta-M\) method cannot be used and therefore it doesn’t work well for highly anisotropic phase functions.

s = iadpython.Sample(g=0.9, quad_pts=4)

hpe, hme = iadpython.hg_elliptic(s)

print(hp)
print()
print(hm)

[[ 1.57558022  1.43478008  0.6702228  -0.09855935]
 [ 1.43478008  1.78554993  1.42038122  0.65844447]
 [ 0.6702228   1.42038122  2.73731157  3.69901785]
 [-0.09855935  0.65844447  3.69901785  6.84908404]]

[[ 1.49114428  1.10817646  0.3889397  -0.08632955]
 [ 1.10817646  0.49081177  0.10073799  0.22945985]
 [ 0.3889397   0.10073799  0.09249519  0.22802648]
 [-0.08632955  0.22945985  0.22802648 -0.37394591]]