Originally published as Biophys J. BioFAST on February 24, 2005.
doi:10.1529/biophysj.104.058305
Biophysical Journal 88:3205-3211 (2005)
© 2005 The Biophysical Society
Estimation of Anisotropic Optical Parameters of Tissue in a Slab Geometry
Olga K. Dudko and
George H. Weiss
Mathematical and Statistical Computing Laboratory, Division of Computational Bioscience, Center for Information Technology, National Institutes of Health, Bethesda, Maryland
Correspondence: Address reprint requests to Dr. Olga K. Dudko, Tel.: 301-402-8702; E-mail: dudko{at}mail.nih.gov.
 |
ABSTRACT
|
|---|
The scattering and absorption coefficients of many homogeneous biological tissues such as muscle, skin, white matter in the brain, and dentin are often anisotropically oriented with respect to their bounding interface. In consequence the curves of equal intensity of re-emitted light on the surface of the slab will no longer be circular. We here consider the problem of determining the parameters allowing one to estimate the angles defining anisotropy, directional bias of diffusive spreading, and scattering and absorbing coefficients from data obtained from time-gated measurements of light intensity transmitted through a slab of the tissue. Our model can be solved exactly and leads to accurate approximations in which measured values of the surface intensity are shown to be elliptical. The parameters of the ellipses suffice to estimate the anisotropy of the tissue interior. A summary of the parameter estimates with the observables from which they are found is given in a table. Our analysis is based on a diffusion model.
|
INTRODUCTION
|
|---|
The use of optical techniques in the red and near-infrared spectrum to probe physical properties of human tissue is attractive because it is non-invasive, inexpensive, real-time, and non-carcinogenic (Benaron et al., 1997
). These techniques, in contrast to other imaging modalities based on ionizing radiation (e.g., x-rays, PET scans), are also sensitive to absorptive properties of tissue, which suggests the possibility of measuring the concentration of endogenous molecules such as hemoglobin (Cope and Delpy, 1988). An obvious drawback in using optical imaging techniques is the blurring effect of photons randomly scattered as they migrate through the target tissue. However, the photons that are scattered and ultimately detected are the ones carrying information about optical properties of the tissue interior. One technique that partially overcomes this problem is the transillumination measurement in which a slab of material is scanned by a point source of light, the emerging time-gated light then being detected at points either on the source plane (the reflectance mode) or the plane opposite (the transmission mode) (Andersson-Engels et al., 1990; Hebden et al., 1997). The first transillumination experiment appears to have been one carried out by Cutler (1929), who reported the ability to differentiate between normal and pathological breast tissue using white light and looking for pathology in light transmitted through the tissue.
Scattering and absorption are the two most significant mechanisms that affect photon transport through tissue. The former is caused by the different refractive indices of tissue components and is generally a dominant mechanism in the absence of tissue inhomogeneities. In many tissues, e.g., muscle, skin, white matter in the brain, and dentin, light is scattered primarily in the direction of the fibers which is one form of anisotropy. Diffuse reflectance measurements have been developed to allow a mapping of fiber anisotropy (Nickell et al., 2000; Kienle et al., 2003).
To date there are several experimental methodologies that permit the mapping of tissue directionality (Tower and Tranquillo, 2001a,b; Napadow et al., 2001). Recently we have enumerated some effects of anisotropy on time-gated and continuous-wave measurements in semi-infinite turbid media to estimate the angles describing the anisotropy of optical parameters in terms of experimental parameters (Dagdug et al., 2003; Dudko et al., 2004). Those analyses dealt with a semi-infinite tissue. An experiment making use of a phantom applied to a semi-infinite medium, described by Sviridov et al. (2004), gave results in good agreement with the theory developed by Dagdug et al. (2003). Similar results were reported for an experiment for a slab geometry in an article by Hebden et al. (2004). The results in that article compared favorably with the theory developed by Dagdug et al. (2003), based on the theory of random walks as applied to optical problems (Gandjbakhche and Weiss, 1995; Weiss et al., 1998). In this article we extend the analysis to treat diffusion in a slab of width L, whose optical properties lie at arbitrary, but fixed, angles with respect to the planar interfaces. Typical slab thicknesses in these experiments are from 40 to 100 mm. Such media can be treated by diffusion theory.
The object of this article is to derive relations between possible measurements that can be utilized to estimate optical parameters of the tissue such as the angles that determine fiber alignment, directional bias, and the scattering and absorption coefficients. The mathematical development makes use of diffusion theory, which is now defined in terms of a diffusion matrix rather than by a single diffusion constant. Our analysis requires a generalization of the method of images to take this feature into account.
|
MODEL FORMULATION
|
|---|
The configuration to be analyzed consists of a slab of width L. The position of an arbitrary point in the slab will be denoted by r = (x, y, z) where the slab faces are at z = 0 and z = L so that z measures the perpendicular distance from the upper face of the slab into the slab interior. The two transverse coordinates satisfy –
x, y . Photons generated by a laser beam initially enter the tissue at the surface point (0,0,0) and a photon is assumed to experience the first scattering at depth z0, so that we can identify z0 as being the scattering length. In a transmission experiment the usable information is based on time-gated measurements made of the intensity profile, Itrans(R;
), measured at the detecting plane at R = (X, Y, L) at a dimensionless time = kt, where k is the scattering rate, which is related to the physical parameters by
 | (1) |
where c is the speed of light in the tissue. Two optical parameters of specific interest are µ's, the transport-corrected scattering coefficient, and µa, the absorption coefficient. We will assume that these are constant along axes skewed with respect to the planar interfaces, as shown in Fig. 1. In the most general case the directionality of the optical parameters is defined by three Euler angles (Goldstein, 1950). However, when the z axis is fixed perpendicular to the boundary planes, two angles suffice to describe the directionality. We define the angles by requiring the x axis be first rotated by an angle 
around the z axis. The next step is a rotation of angle 
around the y axis.
In effect, we need to keep track of two coordinate systems, one defined by the surfaces of the slab, and the second pointing in the direction along which the optical parameters are constant. The first set will be referred to as the laboratory coordinates, i.e., r; and the second as the skewed coordinates. A point in this set of coordinates will be denoted by r' = (x', y', z'). Since r and r' are related by rotations, the transformation between them is affine. Before measurements are made, the experimenter will be assumed to have no information relating to r' since the angles defining the anisotropy are generally unknown. The problem, as so far defined, requires that the diffusion problem be mapped onto a 3 x 3 diffusion matrix rather than being characterized by a single diffusion constant and a single bias parameter. The diffusion matrix will be expressed in terms of the skewed coordinates as
 | (2) |
where B will be termed the bias and D is equal to (3µ's)–1. Notice that, in phrasing the diffusion equation in terms of dimensionless time, the dimension of the diffusion constant is changed to (length)2. We emphasize this point by writing D0 = BD/k as found from Eq. 1, where D0 is the z' component of D0. The dimensionless parameter B is a measure of diffusive spreading along the z' axis as opposed to the spreading along the remaining two axes. It will prove convenient to introduce a dimensionless parameter
 | (3) |
so that isotropy corresponds to B = 1 or 
= 0. When B > 1, diffusive motion along the z' axis proceeds more quickly than that along the remaining two axes.
To derive a formula for the transmitted light at R (on either of the surfaces), boundary conditions must be imposed at the two planar interfaces. We take these to correspond to perfect absorption, so that any photon reaching either surface contributes to the light intensity appearing at the surface. In addition, we assume that photons can be internally absorbed. The absorption kinetics will be assumed to follow the Beer-Lambert law. The survival probability will therefore be written as exp(–
), where
 | (4) |
Any calculation of the light intensity at the surface depends on the propagator, which is the probability density for finding a photon at r at time . Let the propagator at an arbitrary point at dimensionless time be denoted by p(r, ). The flux of photons through the plane opposite to the source will be denoted by Itrans(R;). It is related to the propagator p(r, ) in terms of the flux by Fick's law Itrans(r;) = –D0
p/z, which is to be evaluated either at z = 0 in the reflectance mode, or at z = L in the transmission mode . When B differs from 1 in the transmission mode, the point at which the maximum intensity appears on the detecting surface will shift from (0, 0, L) to another point on the surface that will be determined in the course of the calculation. In the case of the transmission measurement, we will base our analysis on the logarithmic ratio
 | (5) |
which has the advantage that several parameters common to numerator and denominator cancel, thereby simplifying the form of the results. Values of J(X, Y, L;) can be calculated from data available to the experimenter.
There are a potentially large number of candidate relations to be considered as a basis for estimating parameters in a slab geometry, since both transmitted and reflected light can be used for that purpose. In the following exposition we begin by considering the transmission mode, which provides enough information to estimate the required optical parameters. Later we give comparable expressions for reflectance measurements.
|
DERIVATION OF THE DETECTED FLUX
|
|---|
The propagator
Since the transformation between the two sets of coordinates simply involves two successive rotations, we can, by using standard formulae for rotation matrices as in Goldstein (1950), express the resulting set of transformations in terms of the angles and . To abbreviate the notation we use the symbolism c = cos , s = sin , and t = tan . In this notation the forward and backward transformations are
 | (6) |
and
 | (7) |
The derivation of the intensity is based on the propagator in free space,
 | (8) |
as in Dudko et al. (2004). To convert the free space propagator to one satisfying the absorbing boundary conditions at z = 0 and z = L, one needs to successively add and subtract the appropriate set of image terms, so that in the limit of an infinite set of terms the resulting series vanishes at both boundaries. This series generalizes the standard method of images (Weiss, 1994) by incorporating effects of anisotropy. The final expression for the propagator in the slab is
 | (9) |
where z0Cn, z0En, z0Fn, and z0Gn are image points. These are chosen recursively so as to cancel the propagator at 2nL ± z0 (n = ···. –2, –1, 0, 1, 2, ···). Formulae for the image points in Eq. 9 are given in Appendix A.
|
THE TRANSMITTED INTENSITY
|
|---|
An extensive series of transformations, based on the propagator in Eq. 9, and the relation
 | (10) |
leads to an involved expression for the intensity on the detecting surface. Rewriting the expressions for the intensity in terms of laboratory coordinates using the transformations in Eq. 6 we arrive at an exact expression for the transmitted intensity in laboratory coordinates given in Appendix B. Later we will use the exact expression in Eq. B2 for the transmitted intensity by evaluating it numerically as a test of the approximation to be developed below. Now, rather than work with the full expression in Eq. 10, we observe that p(r;) must vanish when z0 = 0. In general, z0 will be small in comparison with the order of magnitude of lengths of interest. Our simplification is formalized by requiring the relation L >> z0 to be satisfied. This validates expanding the intensity around z0 = 0 and retaining only the linear term:
 | (11) |
When the ratio of intensities is taken in Eq. 5 defining the function J(X, Y, L;), the scattering length, z0, cancels out of the result.
One further approximation will be made which considerably simplifies the analysis. Equations 9 and 10 indicate that the expression for the intensity takes the form of an infinite series. Working with such a series, although correct, requires very complicated calculations with results that are hard to interpret. To get around this difficulty, we note, from the expressions for the image points in Eqs. A1 and A2 in Appendix A, that when n differs from zero, each of the terms contains a factor L/z0. For the diffusion model to be valid, the inequality L/z0 >> 1 must be satisfied. As an approximation, we therefore drop all terms with n 
0, leaving us with only a single term for the intensity. This will be referred to as the n = 0 approximation. Later we show, by numerically evaluating the exact expression, that the n = 0 approximation is sufficiently accurate for practical purposes.
The model described to this point has five parameters to be estimated from the data, the angles and , the bias, B, and µ's and µa. These last two provide estimates of D0 and . We suggest a number of relations from which the physically interesting parameters can be derived from transmission mode data. A point to be emphasized is that the set of these relations to be given is not exhaustive.
The approximation to the function J(X, Y, L;) defined in Eq. 5 reduces to a quadratic form in X and Y, which is
 | (12) |
The function 
can be expressed in terms of the linear combination
 | (13) |
as
 | (14) |
The curves of equal intensity (i.e., J = constant) in Eq. 12 are seen to be ellipses, which follows from the n = 0 approximation. This was tested for several sets of parameters by fitting the curves of equal intensity in Eq. 12, with the exact elliptic curves. Fig. 2 shows the effects of changing the parameters , , and L. A change in rotates the ellipse around its center (Fig. 2 a), a change in changes the eccentricity of the ellipse (Fig. 2 b), and a change in L shifts the center of the ellipse and reduces the size of the ellipse as shown in Fig. 2 c. In the case of a completely isotropic medium, B = 1, the curves reduce to circles centered at (0, 0, L).
Three functions describe the effect of anisotropy on the ellipses of equal intensity. The first is the angle 
between the major axis of the ellipse and the x axis. The second is the distance 
between the point (0, 0, L) and the point on the z = L plane at which the intensity is a maximum. The third is the ratio of the minor to the major axis, (R1/R2), of the ellipse.
To obtain expressions in closed form for these functions we apply a rotation to the quadratic form, Eq. 12. That is, we introduce two new coordinates u and 
by writing
 | (15) |
and require that the term proportional to u should vanish. Equation 12 can then be written in terms of the new coordinates u and in the form of the equation for an ellipse,
 | (16) |
where, as seen from Eq. 12, –4D0J is positive. The parameters appearing in Eq. 16 are found in the course of diagonalizing the quadratic form in Eq. 12. They are
 | (17) |
The angle of rotation is found to be
 | (18) |
Two further parameters, measurable from the ellipse, yield relations between B and and between D0 and . The coordinates of the point at which the intensity is a maximum, located at the ellipse centers, are easily found in the n = 0 approximation. These coordinates can be written in terms of B and the angular parameters as
 | (19) |
so that, in this approximation, the distance from origin in laboratory coordinates and (Xm, Ym) is
 | (20) |
which is clearly a measurable quantity. When B = 1 or = 0, which defines the isotropic medium, Eq. 20 indicates that the maximum remains at the origin so that the ellipses of equal intensity reduce to circles whose centers remain at the origin. The parameters Xm and Ym are both proportional to slab thickness. In the limit B 
0 both Xm and Ym go to infinity when = 0. This is to be expected, since this limit corresponds to transport confined to a plane parallel to the slab surfaces. The same limiting behavior is observed in the limit B when is set equal to 90°, since this limit also corresponds to photons constrained to move only in the plane z = z0.
A relation between B and follows directly from Eq. 20 and is
 | (21) |
where is an observable. Equation 21 can be supplemented by a relation obtained by measuring the ratio between minor and major axes of an equi-intensity ellipse. This ratio can be obtained directly from Eq. 16 with 
and ß defined in Eq. 17:
 | (22) |
The combination of this equation and Eq. 21 allows us to find B and separately since R1/R2 can be measured rather easily.
We next observe that, at the point on the detector plane diametrically opposite to the position of the source, we have
 | (23) |
From this it follows that the parameter can be estimated from the long-time behavior of the slope of 
when plotted as a function of . Since = kt, the observed slope, as a function of the physical time t, is
 | (24) |
which can be used to provide an estimate of µa.
Another relation between observable and theoretical parameters is provided in terms of the time at which the maximum value of Itrans(0, 0, L;) is attained. Typical times for such measurements for L = 50 mm are of the order of 600 ps, which are experimentally measurable. Fig. 3 indicates the evolution in time at the set of points (0, Y, L). The figure demonstrates that the intensity has a single maximum as a function of time. The time at which the maximum is attained is found to be
 | (25) |
where
 | (26) |
When absorption is completely negligible so that, in effect, = 0, an expansion of Eq. 25 yields the result
 | (27) |
Equation 25 provides a relation for 
and in terms of the time at which the maximum occurs. This relation allows us to express D0 in terms of B and s2.
Another possibility for relating B, , and D0 (or its equivalent, µ's) can be formulated in terms of J(Xm, Ym;). This has the more complicated form
 | (28) |
The bracketed terms can be regarded as a function of only if the expression for B in Eq. 21 is substituted into the right-hand side of this equation. Therefore, since J(Xm, Ym;) is a measurable quantity, it provides an expression for D0. A rough idea of the accuracy obtainable from the suggested estimators is summarized in Table 1, where Xm(ex) is the value of Xm calculated from the exact solution and Xm(est) is the estimate furnished by the n = 0 approximation. The estimated results agree with the exact ones to within a relative error of 6%.
View this table:
[in this window]
[in a new window]
|
TABLE 1 Values of the coordinates of the center of the ellipses formed by the curves of equal intensity of the detected light
|
|
Reflectance measurements
For completeness we summarize results for reflectance measurements. Except for a few small changes the analysis for reflectance is quite similar to that which has been used to solve the transmission problem, since both are based on the free-space propagator in Eq. 8. The small z0 approximation in Eq. 11 is also invoked in finding the solution for the case of reflectance. In writing the reflectance intensity one replaces the function 
trans in Eq. 13 by
 | (29) |
and the function by 
where, in the n = 0 approximation,
 | (30) |
The later development follows the theory for the case of transmission, with the exception that, since the intensity at the origin is infinite, it can no longer provide usable results for parameter estimates.
|
CONCLUDING REMARKS
|
|---|
The earlier theory of Hebden et al. (2004) allowed for a restricted class of anisotropy consisting of two angles only. In this work we allow for a continuum of values for the two angles that define the anisotropy. This generalization can be analyzed by extending the method of images to encompass arbitrary anisotropy. Although it is possible to solve the resulting equations exactly, the resulting expressions are extremely complicated and it would be hard to generate qualitative information without more extensive numerical calculations. However, a number of approximations are available that produce accurate results when checked against more exact numerical calculations.
One of the approximations allows us to write J(R;) in Eq. 5 as being inversely proportional to D0/L2. In doing so, a term proportional to the square of this quantity has been neglected, which is equivalent to the assumption that L2 >> D0. This, however, is a necessary condition for a diffusion model to be physically realistic (i.e., there must be many scattering events to validate the diffusion picture). A further approximation is based on the neglect of an infinite set of terms in the series defining the propagator, equivalent to the assumption that L >> z0, which again is implicit in the use of a diffusion model. In this approximation the set of curves of equal intensity are ellipses, which is a common feature of data from tissue having anisotropic parameters.
The theory developed can be utilized in a large number of ways to estimate the angles defining our model, which is based on time-gated measurements. A summary of parameter estimates with the observables from which these are found is given in Table 2.
It would be of further interest to explore parameter space to elucidate the set of measurements that are most sensitive to variations in the parameters. Another extension that suggests itself is that of the possibilities available through the use of continuous-wave rather than time-gated measurements.
|
APPENDIX A: EXPRESSIONS FOR THE SERIES OF IMAGE POINTS IN EQ. 9
|
|---|
The image points in Eq. 9 are
 | (A1) |
and
 | (A2) |
It is readily verified, from these expressions, that in the absence of both bias (B = 1) and anisotropy ( = 0) the image vectors reduce to z0Cn = z0Fn = 0, z0En = z0 + 2Ln, and z0Gn = –z0 + 2Ln, as can otherwise be found by the standard method of images. Another significant feature of these results is that the variables here are each proportional to z0, which will allow us to express the results in the much simplified form given in the text.
|
APPENDIX B: DETAILED EXPRESSION FOR THE TRANSMITTED INTENSITY IN LABORATORY COORDINATES
|
|---|
Define the function
 | (B1) |
where trans is defined in Eq. 13. The expression for the intensity then becomes
 | (B2) |
where expressions for the factors Cn, En, Fn, and Gn are specified in Eqs. A1 and A2.
|
ACKNOWLEDGEMENTS
|
|---|
We are deeply indebted to Drs. Victor Chernomordik and Sinisa Pajevic for several enlightening discussions.
Submitted on December 17, 2004;
accepted for publication February 14, 2005.
|
REFERENCES
|
|---|
Andersson-Engels, S., R. Berg, S. Svanberg, and O. Jarlman. 1990. Time-resolved illumination transilluminations for medical diagnosis. Opt. Lett. 18:1179–1181.
Benaron, D. A., C. Wai-Fung, and D. K. Stevenson. 1997. Imaging enhanced tissue optics. Science. 5321:2002–2003.
Cope, M., and D. T. Delpy. 1988. System for long-time measurement of cerebral blood and tissue oxygenation on newborn-infants by near-infrared transillumination. Med. Biol. Eng. Comp. 26:289–294.[CrossRef][Medline]
Cutler, M. 1929. Transillumination of the breast. Surg. Gyn. Obstet. 48:721–727.
Dagdug, L., G. H. Weiss, and A. H. Gandjbakhche. 2003. Effects of anisotropic optical properties on photon migration in structured tissues. Phys. Med. Biol. 48:1361–1370.[CrossRef][Medline]
Dudko, O. K., G. H. Weiss, V. Chernomordik, and A. H. Gandjbakhche. 2004. Photon migration in turbid media with anisotropic optical properties. Phys. Med. Biol. 49:3979–3989.[CrossRef][Medline]
Gandjbakhche, A. H., and G. H. Weiss. 1995. Random walk and diffusion-like models of photon migration in turbid media. Prog. Opt. 34:333–402.
Goldstein, H. 1950. Classical Mechanics. Addison-Wesley, Cambridge, MA.
Hebden, J. C., S. R. Arridge, and D. T. Delpy. 1997. Optical imaging in medicine. Phys. Med. Biol. 42:835–840.
Hebden, J. C., J. J. G. Guerrero, V. Chernomordik, and A. H. Gandjbakhche. 2004. Experimental evaluation of an anisotropic scattering model for a slab geometry. Opt. Lett. 9:2518–2520.[CrossRef]
Kienle, A., F. K. Forster, R. Diebolder, and R. Hibst. 2003. Light propagation in dentin: influence of microstructure on anisotropy. Phys. Med. Biol. 48:N7–N14.[CrossRef][Medline]
Napadow, V. J., Q. Chen, V. Mai, P. T. C. So, and R. J. Gilbert. 2001. Quantitative analysis of three-dimensional-resolved tissue using NMR and optical imaging methods. Biophys. J. 80:2968–2975.[Abstract/Free Full Text]
Nickell, S., M. Hermann, M. Essenpreis, T. J. Farrell, U. Kramer, and M. S. Patterson. 2000. Anisotropy of light propagation in human skin. Phys. Med. Biol. 45:2873–2886.[CrossRef][Medline]
Sviridov, A., V. Chernomordik, M. Russo, A. Eidsath, P. Smith, and A. H. Gandjbakhche. 2005. Intensity profiles of linearly polarized light backscattered from skin and tissue-like phantoms. J. Biomed. Opt. 10:014012.[CrossRef]
Tower, T. T., and R. T. Tranquillo. 2001a. Alignment maps of tissues. I. Microscopic elliptical polarimetry. Biophys. J. 81:2954–2963.[Abstract/Free Full Text]
Tower, T. T., and R. T. Tranquillo. 2001b. Alignment maps of tissues. II. Fast harmonic analysis for imaging. Biophys. J. 81:2964–2971.[Abstract/Free Full Text]
Weiss, G. H., J. M. Porrà, and J. Masoliver. 1998. The continuous-time random walk description of photon motion in an isotropic medium. Opt. Comm. 146:268–276.[CrossRef]
Weiss, G. H. 1994. Aspects and Applications of the Random Walk. North-Holland, Amsterdam, The Netherlands.
Copyright © 2005 by the Biophysical Society.
TOP
APPENDIX B: DETAILED EXPRESSION...
ACKNOWLEDGEMENTS
/3, D0 = 1,