Photon path distributions in turbid media: applications for imaging

Near-IR optical tomography is thwarted by the highly scattering nature of light propagation in tissue. We propose a weighted back-projection method to produce a spatial map of an optical parameter which characterized the investigated medium. We have studied the problem of the choice of the back-projection weight function for the absorption coefficient ((mu) a) and for the reduced scattering coefficient ((mu) s') of tissuelike phantoms. Working in frequency-domain optical imaging, we have initially approached the problem of quantifying the effect caused by a small absorbing defect embedded in the medium on the measured DC intensity, AC amplitude, and phase. The collection of DC, AC, and phase changes during a 1 mm step raster scan of the absorbing defect provides information on the photon path distributions and, in general, on the probed spatial region when DC, AC, and phase are, respectively, the measured parameters. We report experimentally determined weight functions for (mu) a and (mu) s'. They indicate that absorption and scattering maps can significantly differ in terms of resolution.


ETRODUCTION
Optical tomography in the near-infrared is an attractive tool for non-invasive, real time imaging of biological tissue. However, since most tissues act as strongly scattering media for near-infrared light. photons propagate diffusively from light source to optical detector. For this reason, instead of a single photon path, there is a distribution of photon paths inside the medium. This constitutes a complication of the imaging problem treated in the framework of x-ray Computerized Tomography (CT). A possible approach to tackle the inverse problem in optical tomography is to initially find a way to solve the forward problem. Then one can proceed iteratively to find the distribution of values of the optical properties inside the reconstruction region which is consistent with the results of the experimental measurement.13 This approach, which is rigorous and formally correct, presents several difficulties in the critical choice ofthe computer algorithm, in the possibility ofleading to non-unique solutions, and in the long time required to get the result. A simplified method for tomographic reconstruction is a weighted back-projection scheme, which has been applied in optical imaging by Barbour e. at. to obtain intensity images of scattering media.4 '5 In order to obtain spatial maps of a physical property f of the investigated medium, we propose an alternative back-projection scheme in which the measured value of f from a particular source-detector configuration is assigned to a broad spatial region. Each point in this region is also assigned a weight, specified by an appropriate weight function. The tomographic image of 340/SPIE Vol. 2389 O-8194-1736-X/95/$6.OO the measured physical parameter is finally built by superimposing, point by point, the weighted contributions from all source-detector configurations. This paper treats the problem of determining the weight function and relates such weight function to the photon path distributions in the medium.
2.1 X-ray tomography 2. BACK-PROJECTION OPERATOR In x-ray tomography, where the paths of the probing particles (x-photons) from source to detector can be assumed to be straight lines, it is possible to define a back-projection operator in terms of onedimensional integrals. Figure 1 defines the rectangular (x,y) and the polar (r,4) coordinates used in modeling the imaging problem. L is the straight line joining source and detector, K is the line through the origin (0) and perpendicular to L, and P is the intersection point between L and K. Let f(r,) be a physical parameter which is function of the position in the reconstruction region (1), and let us assume ftr,) = 0 for (r,) T. A projection operator, which projects all the values of J(r,4) in L on the line K, can be defined in terms of the Radon transform R (Ref. [6]): [f](,a) = ,a + tan1(z / In Eq. (2) the integral sums up the contributions from all the lines L passing through (r,).

Near-infrared tomography
In near-infrared optical tomography, the photon paths in tissue are distributed over a 3 -D region, which is sometimes called a light bundle, so that the line integral in Eq. (1) is replaced by a volume integral. Furthermore, the photon path density is not uniformly distributed over the light bundle, and for this reason a position dependent weight function J4'(r) must be introduced. In optical imaging, a normalized operator (Z) which generalizes Eq. (1) can thereby be defined as a weighted average of the values ofj(r) over the volume ofthe light bundle (i.e. over the domain of Wt (r)): where cx is the six-component vector which defines the source-detector geometrical configuration (by defining the positions of source and detector). Similar to the x-ray case, a back-projection operator can be defined. For a discrete number N of source-detector configurations a.., it can be written as follows: The approach to the problem of imaging presented in this Section requires that the contributions to the measured value off from each voxel is independent of that from all the other voxels. In optical tomography this is not true, but as long as small perturbations are considered, it constitutes a reasonable approximation.

Continuous wave (CW) light
Light propagation in turbid media is modeled by the diffusion approximation to the Boltzmann transport equation. Analytical solutions for the photon flux distribution in macroscopically homogeneous media and for several boundary conditions have been determined. 7 The presence of macroscopic inhomogeneities which have small influence on the measured quantities have been treated within the framework of perturbation theory.8 '9 In the case that the inhomogeneity is a small (linear dimensions much less than source-detector distance) totally absorbing defect, one can assume that the decrease in intensity (defined as l-IDC(r)/IDCO where IDC(r) is the detected intensity when the defect is in r and 'DCO is the detected intensity in the absence of the defect) caused by this defect is proportional to the photon path density at the position of the defect. In this way, the photon path distributions in turbid media have been studied both experimentally10 and theoretically.8 In this perspective, the probability PDC(r,rS,rd) that a photon emitted at r = r,, (source position) and detected at r = rd (detector position) passes through r is analytically given, in the infinite geometry and for CW light, by Ref. [9] as: where cDC(r,rS,r) is defined as IDC(r)/IDCO, a is the absorption coefficient of the medium, D is the diffusion coefficient defined as 11(3 +3 ') with ji' reduced scattering coefficient, Ajt is the difference between the absorption coefficient of the defect and that of the background, and zV is the volume occupied by the defect. The probability pDC(r,S,rd) does not directly define the light bundle. We can define the light bundle by giving its sectional surface S; at every point i of the line joining source and detector, which we will call the x-axis. In the infinite geometry, we define S as the circle in the plane x = i and with the center in the x-axis which obeys: JfpDc(i,y,z,r, rd)dydz =0.8 (6) JJPDC(X,Y,Z, r, rd)dydz where the value 0.8 is arbitrarily chosen and corresponds to the selection of 80% of the photon paths through x = SPIE Vol. 2389 / 343

Intensity modulated light
In frequencydomain imaging, the intensity of the light source is modulated at a frequency of the order of 100 MHz, and, besides average intensity (DC), there are two additional measurable parameters: the amplitude of the intensity oscillations (AC) and the phase () of the intensity wave at the modulation frequency. The effect of a totally absorbing defect on the measured value of AC is in general different than that on the measured DC. The difference is a function of the modulation frequency. It is possible to express the variations in measured AC and phase (6Ac and c, respectively) due to the presence of a totally absorbing point-like defect in r. cAC(r,rS,r) is defined as IAC(r)/IACO, while c(r,rS,rd) is defined as the difference '1(r)-c10 with c1(r) phase measured with the defect in r and phase measured in the absence ofthe defect. By following the approach presented in Ref. [9] we have: cAc(r,rS,rd)=Abs) We note the different physical meaning ofthe three quantities c(r,r,r), 6AC(r,rS,r) and c(r,rS,rd). cDC(r,S,rd), as already discussed, is related to the photon path distribution in the medium in CW conditions; CAC(r,rS,rd) has a similar interpretation, but it refers to the AC signal detected in the frequency-domain. The physical reason for the difference between cDC(r,r,r) and cAC(r,rS,r) is that the attenuation of the DC (or CW) component oflight intensity is smaller than that of the AC component of light intensity; c(r,r,r) gives information on the spatial region effectively probed when the phase is the measured parameter. Larger values of e(r,rS,rd) mean a larger effect on the measured phase due to the deletion of photon paths through r. In Fig. 2, panels (a)  and c(r,r,r) plotted in Fig. 2 (a), (b), and (c)). In Fig. 2 (d), (e), and (1)

DETERMINATION OF THE WEIGHT FUNCTION W$(r)
The DC, AC, and phase defect-induced variations discussed in the last Section give information on the sensitivity of each directly measured quantity to the deletion of photon paths through a specific point of space. Ifflr) is a physical parameter of the medium and fm S its measured value with a particular source-detector pair a as a function of DC, AC, and phase, we can write f = DC ,AC, t) . We define the weight function forJ(r) as: where F0 is the measured value offin the absence ofthe defect. By taking into account the fact that we are considering small perturbations we can also write the weight function in terms of the changes in DC, AC, and phase as: where DC0, AC0, and are the measured DC, AC, and J? in the absence ofthe defect. We have applied Eq. (9) to the results ofthe experiment in which we scanned the small absorbing defect in the plane containing source and detector. For each position of the defect we calculated and t' with the pre-calibrated measurement protocol described 1 1 The results for the weight functions for j.x and j.j', which quantify the effect of the deletion of photon paths through a specific point of space on the recovered values of ji and are plotted in

DISCUSSION AND CONCLUSION
We stress that Fig. 2 does not provide a representation of the light bundle, which is rigorously defined by Eq. (bundle) in terms ofp(r,r,r). This light bundle, which contains the majority of the photon paths from source to detector, specifies the geometrical properties of diffusive light propagation from source to detector. We are currently exploring the connection among these geometrical properties, the wavelength of the photon-density wave, and the resolution attainable in optical imaging. In the case in which DC, AC, and phase are combined to yield a physical quantity, the weight function for that quantity provides resolution information. In this perspective, Fig. 3 suggests that, at least in our experimental conditions, a scattering image should be better resolved than an absorption image. In fact, this is what we found in a series of experiments on phantoms containing scattering and absorbing nh1 A problem in the determination of the weight function for the quantity Jis that the weight function depends on the value off of the background medium. This difficulty could be overcome by performing a measurement of the average value off over the reconstruction region and by calculating the weight function using that average value. We are currently evaluating the sensitivity of W and W on changes in the optical properties ofthe background medium.
The back-projection method presented in this paper constitutes a simplification of the mathematically rigorous reconstruction algorithm based on the inversion of the integral operator SPIE Vol. 2389 I 347 defined in Eq. (3). The major advantage of the back-projection method is that it is fast, in terms of computation time, and always yields a unique solution, For this reason, in addition to fast acquisition times for optical data, can be used to provide real time optical tomography, which would constitute a considerable step forward with respect to current tomographic techniques. We estimate that an optical image of a region with linear dimensions of about 5 cm could be updated every few seconds. On the other hand, the back-projection method is likely to produce blurred images. The reason for this is that we assign a single reading value, which could be mainly due to a localized inhomogeneity, to a broad spatial region. An appropriate weight function will partially correct for this problem, but it still remains as a characteristic feature of back-projection. Also, a highly resolved and contrasted image requires a high number N of source detector configurations, i.e. the summation in Eq. (4) should contain many significant (J1) 0) terms.
In conclusion, we have described a weighted back-projection scheme designed to yield a spatial map of a physical property ofthe investigated medium. Working in frequency-domain imaging, we have studied the effect ofthe deletion ofthe photon paths through a specific point of space on the measured DC, AC, and phase. In particular, we have related the effect on the DC intensity to the photon path distributions between light source and optical detector. Finally, we have defined and experimentally derived the weight functions for absorption (W,(r)) and for reduced scattering coefficient (W'(r)). They should be interpreted as the relative change in the measured values of ha and respectively, caused by the deletion of the photon paths through r. The shape of these weight functions resembles that of EDC(r) (i.e. IDC(r)/IDCO) but their width is determined by the particular functional dependence on the phase.

ACKNOWLEDGMENTS
This work was performed at the Laboratory for Fluorescence Dynamics at the University of Illinois at Urbana-Champaign (UJUC), which is supported by the National Institutes of Health (NIH), grant RRO3 155 and by TJITJC. This research is also supported by grant CA57032 from the NIH.