Diffuse optical 3D-slice imaging of bounded turbid media using a new integro-differential equation

A new integro-differential equation for diffuse photon density waves (DPDW) is derived within the diffusion approximation. The new equation applies to inhomogeneous bounded turbid media. Interestingly, it does not contain any terms involving gradients of the light diffusion coefficient. The integro-differential equation for diffusive waves is used to develop a 3D-slice imaging algorithm based on the angular spectrum representation in the parallel plate geometry. The algorithm may be useful for near infrared optical imaging of breast tissue, and is applicable to other diagnostics such as ultrasound and microwave imaging ©1999 Optical Society of America OCIS codes: (170.3010) Image reconstruction techniques; (170.5280) photon migration; (170.3830) mammography; (170.5270) Photon density waves; (170.3660) Light propagation in tissues References and links 1. B. Chance, Q. M. Luo, S. Nioka, D. C. Alsop and J. A. Detre, “Optical investigations of physiology: a study of intrinsic and extrinsic biomedical contrast," Phil. Trans. Royy. Soc. London B. 352, 707 (1997). 2. A. Villringer and B. Chance, “Non-invasive optical spectroscopy and imaging of human brain functions," Trends. Neurosci. 20, 435 (1997) . 3. Y. Hoshi and M. Tamura, “Near-Infrared Optical Detection of Sequential Brain Activation in The Prefrontal cortex during mental tasks," Neuroimage. 5, 292 (1997). 4. J. H. Hoogenraad, M. B.van der Mark, S. B.Colak, G.W.'t Hooft, E.S.van der Linden, ``First Results from the Philips Optical Mammoscope'', Proc.SPIE / BiOS-97 (SanRemo, 1997 ). 5. J. B. Fishkin, O. Coquoz, E. R. Anderson, M. Brenner and B. J. Tromberg, “Frequency-domain photon migration measurements of normal and malignant tissue optical properties in a human subject," Appl. Opt. 36, 10 (1997). 6. M. A. Franceschini, K. T. Moesta, S. Fantini, G. Gaida, E. Gratton, H. Jess, W. W. Mantulin, M. Seeber, P. M. Schlag and M. Kaschke. “Frequency-domain instrumentation techniques enhance optical mammography: Initial clinical results” Proc. Natl. Ac ad. Sci. USA, 94, 6468-6473 (1997). 7. S. K. Gayen and M. E.Zevallos, B. B. Das, R. R. Alfano, “Time-sliced transillumination imaging of normal and cancerous breast tissues," in Trends in Opt. And Photonics , ed. J. G. Fujimoto, M. S. Patterson. 8. B. W. Pogue, M. Testorf, T. McBride, U. Osterberg and K. Paulsen ,” Instrumentation and design of a frequency-domain diffuse optical tomography imager for breast cancer detection ," Opt. Express 1, 391 (December 1997). http://epubs.osa.org/oearchive/source/2827.htm 9. D. Grosenick, H. Wabnitz, H. H. Rinneberg, K. T. Moesta and P. M. Schlag, “Imaging and Characterization of Breast tumors using a laser-pulse mammograph," SPIE 3597 (1999). 10. A. G. Yodh, B. Chance, “Spectroscopy and imaging with diffusing light," Phys. Today 48, 34-40 (March 1995). #9231 $15.00 US Received March 05, 1999; Revised April 09, 1999 (C) 1999 OSA 12 April 1999 / Vol. 4, No. 8 / OPTICS EXPRESS 231 ______________________________________________________________________________________________ 11. H. Jiang, K. D. Paulsen, U. L. Osterberg, B. W. Pogue and M.S. Patterson, “Optical Image reconstruction using frequency-domain data: Simulations and experiments," J. Opt. Soc. Am. A 13, 253-266 (1996). 12. Y. Yao, Y. Wang, Y. Pei, W. Zhu and R. L. Barbour, “Frequency-domain optical imaging of absorption and scattering distributions using a born iterative mehod," J. Opt. Soc. Am. A 14, 325-342 (1997). 13. S. R. Arridge, “Forward and inverse problems in time-resolved infrared imaging," in Medical Optical Tomography: Functional Imaging and Monitoring , ed. G. Muller, B. Chance, Rl. Alfano, S. Arridge, J. Beuthan, E. Gratton, M Kaschke, B. Masters, S. Svanberg, P. van der Zee, Proc SPIE IS11, 35-64 (1993). 14. An integro-differential equation valid for infinite medium,but containing no gradient of Diffusion constant was used for image reconstruction by M. A. O’Leary, D. A . Boas, B. Chance, A.G. Yodh, “Experimental Images of heterogeneous turbid media by frequency-domain diffusing-photon tomotography," Opt. Lett. 20, 426-428 (1995). 15. An anonymous reviewer has brought to our attention a recent paper by S.R.Arridge and W.R.B. Lionheart, "Non-uniqueness in diffusion based optical tomography," Opt. Lett. 23, 882-884 (1998), which deals with the second order derivative of the diffusion coefficient. 16. D. N.Pattanayak and E.Wolf , " Resonance States as Solutions of the Schrodinger Equation with a Non-Local Boundary Condition, " Phys. Rev. D13, 2287(1976). 17. J. L.Ye,R.P.Millane, K. J.Webb and T. J.Downar, "Importance of the gradD term in frequency resolved optical diffusion imaging,"Opt. Letters,23,1998. 18. K. M.Case and P. F.Zweifel, in Linear Transport Theory (Addison -Wesley,MA,1967). 19. A.Ishimaru, Wave Propagation and Scattering in Random Media (Academic Press, New York, 1978) Vol.. 19. 20. P. M.Morse and H. Feshbach, Methods of Theoretical Physics (McGraw-Hill,1953). 21. Claus Muller, Foundations of the Mathematical theory of Electromagnetic Waves (Springer-Verlag,New York,1969). 22. H. Weyl, Ann. Phys. 60, 481, (1919). 23. See, also , A.J. Banos, in Dipole Radiation In the Presence of a Conducting Half-Space (Pergamon Press, New York, 1966). 24. E. Wolf, "Principles and Development of Diffraction Tomography” in Trends in Optics , ed. A. Consortini (Academic Press, San Diego, 1996). 25. A. J.Devaney, “Reconstructive tomography with diffracting Wavefields," Inv. Probl.,2,161-18391986). 26. For an application of diffraction tomography to near field diffusion wave imaging and an analytical expression relating theoretical resolution and tissue thickness see, GE Class I Technical report (publicly available on request) by D. N. Pattanayak, “Resolution of Optical Images Formed by Diffusion Waves in Highly Scattering Media," GE Tech. Info. Series 91CRD241 (1991). 27. C. L. Matson, N. Clark, L. McMackin and J. S. Fender, “Three-dimensional Tumor Localization in Thick Tissue with The Use of Diffuse Photon-Density Waves," Appl. Opt. 36, 214-219 (1997). 28. X.D.Li, T. Durduran, A.G. Yodh, B. Chance and D.N. Pattanayak, “Diffraction Tomography for Biomedical Imaging With Diffuse Photon Density Waves," Opt. Lett. 22, 573-575 (1997). 29. Chen BQ, Stamnes JJ, Stamnes K, “Reconstruction algorithm for diffraction tomography of diffuse photon density waves in a random medium”. Pure Appl. Opt. 7, 1161-1180 (1998). 30. X. Cheng and D. Boas, “Diffuse Optical Reflection Tomography Using Continous Wave Illumination," Opt. Express 3, 118-123 (1998); http://epubs.osa.org/oearchive/source/5663.htm 31. S. J. Norton, T. Vo-Dinh, “Diffraction Tomographic Imaging With Photon Density Waves: an Explicit Solution, J. Opt. Soc. Am. A 15, 2670-2677 (1998). 32. J. C. Schotland, “Continuos Wave Diffusion Imaging," J. Opt. Soc. Am. A 14, 275-279 (1997). 33. J. Ripoll, M. Nieto-Vesperinas, “Reflection and Transmission Coefficients of Diffuse Photon Density Waves," (to be published). 34. J. Ripoll and M. Nieto-Vesperinas, “Spatial Resolution of Diffuse Photon Density Waves, " J. Opt. Soc. Am.A (to be published). 35. ]C.L.Matson and H. Liu, "Analysis of the forward problem with diffuse photon density waves in turbid media by use of a diffraction tomography model," J. Opt. Soc. Am. A 16, 455-466 (1999). 36. For the angular spectrum representation in a slab field and for an excellent account of the angular spectrum representation see, L.Mandel and E.Wolf, Optical Coherence and Quantum Optics (Cambridge University Press, New York,1995).


Introduction
Tomographic imaging of deep tissues with diffusive waves makes possible functional imaging based on novel optical contrasts derived from tissue chromophores, structure and metabolism [1][2][3]. The near-infrared spectroscopic properties of breast tumors for example, have been found to differ from adjacent normal tissues. Such spectroscopic signatures hold promise for increased tumor detection sensitivity and specificity [4][5][6][7][8][9][10]. In order to take full advantage of these new contrasts it is critical to develop high fidelity three-dimensional optical imaging algorithms for diffusing light in turbid media such as the breast. To this end a number of image reconstruction schemes have been developed [11][12][13]. Many of these schemes use integral equations based on Born or Rytov approximations to generate a set of linear equations which are then solved to update the absorption and scattering coefficients associated with each voxel in the reconstructed volume. The integral formulation is attractive because of its speed. However, to our knowledge most of these formulations ignore terms involving gradients of the light diffusion coefficient [14,15] Although this approximation is often reasonably accurate, a recent paper [17] has established the importance of this gradient term, showing that its neglect is responsible for cross coupling between absorption and scattering images.
In this paper, we derive a new integro-differential equation for diffuse photon density waves (DPDW) within the diffusion approximation. The new equation, which is developed for bounded turbid media does not explicitly contain any terms involving gradients of the light diffusion coefficient. We then use this integro-differential equation and develop a theoretical inverse scattering algorithm for three-dimensional image reconstruction. The theoretical framework follows the principles of near-field diffraction tomography based on the angular-spectrum representation; in particular we show how to develop three-dimensional slice images of a breast compressed between two parallel plates. The technique employs a series of plane diffuse photon density waves with different modulation frequencies. For this set of plane incident diffusive waves, the algorithm requires two-dimensional FFT operations, and a one-dimensional matrix inversion. Thus the method is non-iterative in two space dimensions and is therefore computationally fast.

Integro-differential equation for diffuse photon density waves
We begin with the equations satisfied by the diffuse photon density, Φ and the diffuse photon flux density, J Here, ω denotes the modulation frequency of the source intensity, µ α and µ s ′ are the absorption and transport scattering coefficients of the bounded turbid media, c is the velocity of light in the medium and S 0 (r,ω ), S 1 (r,ω) are the monopole and dipole moments of the intensity modulated optical source respectively. From these equations we eliminate the photon flux density and obtain, after a lengthy calculation, the following wave equation involving only the photon density function is the square of the complex wavenumber k . The bar over a physical parameter denotes d spatial average of that parameter. Notice that for unmodulated source intensity (the continuous dc domain in contrast to the frequency domain), the photon density is a nonpropagating damped wave. The wave number in this case is purely imaginary. For frequency modulated source intensity, the wave number is complex and the photon density acquires the characteristics of a damped or attenuated wave. Because of its differential nature Eq. (3) will admit many solutions and boundary conditions are necessary to obtain the physical solution. The problem of finding the proper boundary condition is complex and is generally further exacerbated because the diffusion approximation is not valid near the boundary. Inspite of this basic difficulty, diffusion theory is used widely due to its simplicity and the applicability of numerous analytical and numerical techniques for its solution in finite domains of arbitrary geometry. Approximate boundary conditions are generally adopted which have been shown to be fairly accurate. Here, we adopt the general boundary condition [19] & & n ⋅ J & (r ,ω ) = αΦ(r ,ω ).

˘
The value of α can be defined based on boundary considerations or taken as a fitting parameter for a given interface. In Eq. (5) n is the unit outward normal to the boundary surface. Using Eq. (5), it is now possible to transform Eq. (3) into an integral equation, the solution of which will provide one with a physical picture of the interaction of the diffuse photon density wave with the turbid media. The basic interaction of light with the molecules of the turbid media is already included in the absorption and scattering parameters. The determination of these material parameters is at the heart of the optical modality for cancer diagnostics. The integral representation forms a basis for extraction of these parameters from measurements of the diffusion waves at the boundary. We next derive the integro-differential equation using Eqs. (3) and (5). We follow the standard Green's function technique for the derivation of an integral equation from a partial differential wave equation [20]. However, we deliberately choose a simple Green's function corresponding to the operator on the left hand side (L.H.S.) of Eq.
(3) that satisfies the equation and is given by Here Im(k ) > 0. Notice that our choice of an infinite space Green's function will not limit In deriving Eq. (9), we also made use of Eqs. (1)-(2) and (5). V denotes the volume of the turbid medium and the surface bounding the volume is denoted by S. The first two terms on the right hand side (R.H.S.) of Eq.(9) are related directly to the source and the second term on the R.H.S. of Eq.(9) arises because of the presence of the boundary. The solution to this integro-differential equation is unique as long as there exists no non-trivial solutions to the corresponding homogeneous integral equation. It is possible that for certain modulation frequencies the resulting homogeneous integral equation may possess non-trivial solutions. These frequencies will generally be complex and will correspond to resonance behavior of the photon density waves. We shall not dwell with this singular case, but refer the reader to a related paper on resonance theory by one of the authors [15]. Our Eq.(9) differs from the corresponding equation used in reference [14], in which the boundary effect was ignored. The treatment of the source contribution to the integral equation is also different.
For a homogeneous diffusing media, the integro-differential equation becomes  (9) and (10) are the basic equations of our paper. They are derived using the diffusion approximation to photon transport and without any other approximations. Eq. (9) is an intro-differential equation for the photon density wave and involves only the absorption and transport scattering coefficients. Notice that Eq. (9) does not explicitly contain gradients of light diffusion coefficient. Our derivation is exact and does not assume any condition on the gradient of the light diffusion coefficient. As has been pointed out recently [17] neglect of the gradient term leads to cross coupling between scattering and absorption images. Eq. to be on the surface. This limiting operation must be done carefully because the integrals involving the gradient of the Green's function are not continuous when the interior point approaches a point on the surface. The other integrals involving the Green's function however are continuous [21]. From an experimental point of view, if the photon density on the surface surrounding the turbid medium is measured, then one can use Eq. (10) to determine the diffuse photon density wave anywhere inside the medium by simply carrying out the surface integrals in Eq. (10). If the value of this estimated photon density matches that of the measured value at all points on the surface, then an absence of tissue heterogeneity is indicated. In practice there will never be perfect cancellation, but a larger difference indicates the presence of stronger heterogeneities inside the tissue.
Another approximate way to evaluate the background field is to evaluate the surface integrals in Eq. (10) iteratively. First ignore the third surface term involving the boundary values of the photon density in Eq. (10) and use only the first two terms involving the source terms to calculate the first order background field everywhere. Then update the value of the background field by using the first order values in the last integral in Eq. (10). It is important to devise a way to estimate a reasonable background field, because when this field is subtracted from the total measured diffuse photon density field, then the remainder is the field directly attributable to the presence of the heterogeneity.
We have already noted that unlike other integral formulations our integro-differential equation does not contain an explicit term involving the gradient of the transport scattering coefficient. This is important in the inverse problem where the absorption and transport scattering heterogeneity are to be estimated from the measured values of the field. The formalisms that include the gradient term will of necessity have to use a numerical estimate of the gradient value. Such procedures may lead to numerical instabilities and numerical artifacts. We next develop a 3D-slice image reconstruction technique based on the integrodifferential equation for the photon density, i.e. the Eqs. (9) & (10).

3.Inversion algorithm for 3D-slice imaging
We make the assumption that the heterogeneities are weak and do not perturb the background field substantially. Then we can approximate the photon density within the integral operator with the background field. With this first order Born approximation, Eq. (9) becomes The superscript M implies measured field and the superscript P denotes the processed (scattered) field values. In Eq.(12), we have used the identity The quantity Φ (r,ω) often identified as the background field can be determined from the knowledge of the source and the measured field at the boundary from Eq. (12). Although our integro-differential equation is valid for arbitrary geometry, from now on we shall restrict ourselves to the case of a plane parallel plate configuration. The plane parallel plate is an accepted clinical diagnostic configuration for breast imaging. It offers the possibility of compressing the breast to decrease its thickness for increased signal strength. We also take the transverse dimension of the plate to be sufficiently larger than the plate separation. Because of the negligibly small value of the field at the periphery of sufficiently large plates, in practice the transverse plate extent can be taken to be infinite. In essence, Eq.
& & (11) is an integral equation for µ (r ) , and µ ′(r ) with all the other variables such as the a s background field, the Green's function and the background optical parameters either given or estimated.
In the new integro-differential equation, we use the Weyl representation [22][23] of The direction z is chosen perpendicular to the surface normal and the x and y directions are chosen along the transverse directions. There is a singularity in the z direction as can be seen from the absolute value operator in the exponent of the Weyl representation. This representation is also known as the angular spectrum representation and is fundamental to the theory and algorithm now commonly known as Diffraction tomography [24][25]. Diffraction tomography was applied to far field coherent optical imaging systems, but recently this technique has been extended to the near field diffuse photon density wave based imaging [26][27][28][29][30][31][32][33][34]. With this representation of Green's function the integrals in x and y become convolution integrals which turn into algebraic products in p, q space. Taking the two-dimensional  Fourier transform of the so-called tumor functions, T a and T which are products of the s scattering parameters and the background field in real space, will be convolution integrals in the Fourier space. This convolution in the spatial frequency space between the tumor parameters and the background field couples temporal frequency, ω with spatial frequency (p,q). This implies that the use of a set of modulation frequencies does not in any easy manner let us solve for the spatial Fourier amplitudes of the tumor functions. This coupling is eliminated if the background field is spatially in the form of a plane wave rather than a spherical wave. The usefulness of plane wave sources is well known in the field of diffraction tomography [24]. We will use a series of plane waves whose central modulation frequencies are stepped to span some frequency bandwidth. For the plane wave case the tumor functions, Tˆ, Tˆ , can a s be written as a product of two functions. One function, t or t D , contains information about a the heterogeneities and depends upon (p, q, z) and the other function depends upon the background plane waves and modulation frequency (see Eq. (18) below). One can then utilize this temporal frequency bandwidth to obtain resolution in the longitudinal (z) direction. The key advantage of using the series of stepped central modulation frequencies is that the detection system is still narrow band and one can use maximum permissible signal strength at all discrete frequencies. We plan to publish in a later paper the details of this technique with simulated data, but here we point out the salient features of the method. For the plane wave case we find   10)) to be a sum of both forward and backward propagating waves [33] Φ ( x, y, z,ω) Here, Φ + , and Φ − are complex constants obtained from solution or measurements as discussed earlier in the section on the background field. We approximate the integration in Eq. (14) in the z variable by means of a sum over N slices between the parallel plates (the slice thickness can be different and is denoted by Δz j with j ranging from 1: N), and use Eqs. (15)- (19) to obtain the following result matrix inversion as shown below or iteratively by following methods such as SIRT, ART [13] . t consists of the sum of the two N dimensional vectors t and t D . For the Breast imaging case, Eq. (24) provides one with 3D-slice images of the absorption and scattering coefficients at the selected "virtual" slice positions in the interior of the breast compressed between the two plates. These images may be reconstructed with light of different wavelengths within the near-infrared window. Such measurements improve the usefulness of nearfield optical imaging modality employing diffuse photon density waves.

4.Conclusion
We presented a rigorous derivation of a new integro-differential equation for the diffuse photon density wave in an inhomogeneous bounded turbid medium within the diffusion approximation. The resulting equation contains only the absorption and scattering coefficients and does not include gradients of the light diffusion coefficient. We also derived a novel surface integral equation for the diffuse photon density waves for the case of a bounded homogeneous turbid medium. We used the theory to develop a detailed imaging algorithm that provides one with 3D-slice images of the absorption and scattering coefficients of tissue (breast) that is compressed between two parallel plates. We believe this algorithm will be useful in increasing the specificity of the near infrared optical imaging modality. Our algorithm being based on inhomogeneous wave equation of the Helmholtz type may also be applicable to other imaging modalities such as ultrasound.