The instrument development and design of a prototype frequency-domain optical imaging device for breast cancer detection is described in detail. This device employs radio-frequency intensity modulated near-infrared light to image quantitatively both the scattering and absorption coefficients of tissue. The functioning components of the system include a laser diode and a photomultiplier tube, which are multiplexed automatically through 32 large core fiber optic bundles using high precision linear translation stages. Image reconstruction is based on a finite element solution of the diffusion equation. This tool for solving the forward problem of photon migration is coupled to an iterative optical property estimation algorithm, which uses a Levenberg-Marquardt routine with total variation minimization. The result of this development is an automated frequency-domain optical imager for computed tomography which produces quantitatively accurate images of the test phantoms used to date. This paper is a description and characterization of an automated frequency-domain computed tomography scanner, which is more quantitative than earlier systems used in diaphanography because of the combination of intensity modulated signal detection and iterative image reconstruction.
©1997 Optical Society of America
In recent years, there has been a significant amount of research on the development of improved methods for optical imaging of human tissues. One of the most promising applications of imaging with diffuse light is cancer-detection in the human breast. Transillumination, or diaphanography, have been used in various forms for several decades [1–4], in an attempt to use the transmitted intensity of red or near-infrared light to image the interior of the breast. These early attempts ended with little success mainly because of poor image resolution and an inability to detect lesions less than 2 cm in diameter [4–6]. Since the design of these early imaging systems, methods to improve image resolution and detectability have been realized through the use of computed tomography reconstruction [7–14]. This paper describes the essential features of an optical breast cancer imaging instrument which is being developed and will be evaluated in early clinical trials.
In general, there have been two main theoretical developments in diffuse optical imaging, namely improved methods for solving both the forward problem to predict light propagation in tissue [15–19], and the inverse problem for estimating the tissue optical properties based upon projection measurements along the surface [8–10, 20, 21]. The forward prediction of light propagation in a highly scattering medium such as tissue has been modeled accurately using several radiation transport methods. Perhaps the most simple and yet generally applicable for near-infrared light propagation in tissue is the diffusion model of light transport. In the present prototype imaging method, a finite element solution of the diffusion equation has been demonstrated to accurately predict the detected light remitted from tissue phantoms, and serves as the basis for the computational reconstruction of images [10–13]. Efforts to understand and predict light propagation in tissue, using diffusion theory solutions have also been undertaken here. Independent confirmation of how light propagates in tissue based on a Monte Carlo statistical approach has been included to confirm the quality of the measurements obtained with the new data acquisition system.
The instrument design has followed many of the recent developments in frequency-domain tissue spectroscopy [22–27]. The use of intensity modulated light allows a measurement of two parameters, namely intensity and phase, which are influenced uniquely by the absorption and diffusion coefficients. These two sets of measurements then permit quantitative optical property reconstruction. The major technological objective is to design a source-detector array which is low in cost while also being able to acquire multiple measurements rapidly. A single source-single detector scheme has been used here which incorporates fast multiplexing to minimize the data acquisition time, yet is able to record 256 projection measurements through the tissue. The system has been specifically designed to maximize light detection, thereby maximizing the thickness of tissue which can be imaged.
This paper outlines the design and development of a fully functional breast cancer screening device which should be able to provide quantitative images of absorption and diffusion coefficient maps, and the use of multiple projections with a reconstruction algorithm should provide increased spatial resolution over diaphanography [10-14]. While this type of imaging has been studied extensively both theoretically and in test phantoms, there have been few descriptions of automated systems for computational image reconstruction for breast cancer detection presented in the published literature [28,29]. This type of system is described here along with some representative images obtained from test phantoms. Plans for future clinical use are discussed.
2. Instrument design
The light source for the optical imaging system is a diode laser (SDL Inc., CA) at 800 nm wavelength and output power of 100 mW. The current to this laser is modulated by mixing continuous and radio frequency (RF) current at a bias Tee (Picosecond Pulse Labs model 5545, Boulder CO). The RF current is supplied by a signal generator (Marconi Instruments model 2023, Ft. Worth TX), typically run at 100 MHz and power output of 13 dBm. Light is delivered to and from the tissue phantom via fiber optic bundles with 1/4 inch (approx. 6 mm) inner diameter (Cuda Products, FL).
The light detection is accomplished with a single photomulitiplier tube (PMT) (Hamamatsu Inc model R928, Japan), using a specially designed RF modulated housing (ISS Instruments, Champaign, IL). The detector is modulated at a slightly offset frequency, so that detection can be achieved at the heterodyne beat frequency, which in this case is an offset frequency of 100.001 MHz which yields a beat frequency of 1 kHz.
The source and detector are multiplexed rapidly through the use of 16 source optical bundles and 16 detector optical bundles (see Fig. 1 and Fig. 2). The fibers were translated laterally in front of both the source and the detector via two linear translation stages, which can be moved at 6 inches per second (approx. 15.3 cm/sec) and with a positioning accuracy of 5 microns (PCLM series SP-4, Anorad Corp., Hauppage NY). For the detection multiplexing, the system has to be able to measure intensity changes over 6 orders of magnitude. To increase the dynamic range of intensities that the PMT can detect, an automated filter wheel fitted with pre-calibrated neutral density filters was added to the system (Oriel Instruments model 77371, Stratford CT). By automatically adjusting the filter before the detector to achieve maximal signal, the dynamic range of intensities has been increased from 3 to 8 orders of magnitude. An additional advantage of the configuration used here, is that high voltage components are spatially separated from the imaging array, which contacts the patient.
The current signal from the PMT is sampled with a 10 kHz data acquisition board using an alternating sampling of the reference and phantom detectors. Before detection, the signal is passed through a low pass filter with a cut off of 4 kHz to reduce the high frequency noise.
2.2 Data Aquisition
All data acquisition, signal processing and instrument control subroutines were written within the LABVIEW language (National Instruments, Austin TX). For each measurement, a 100 ms sample of both the source and detector signals are taken at 10 kHz, and these data samples are processed with software for lock-in detection at 1 kHz to detect the real and imaginary components of the AC signal. Preliminary measurements are shown in Fig. 3 (discussed in Section 2.2) for typical light signals using this short acquisition time. Currently the shifting time between samples is 0.5 sec, predominantly due to the linear translation stages.
Approximately 30 seconds are required additionally for all filter changes during an image acquisition. The total data acquisition time for 256 measurements is 7 minutes, at the present time. Once all the real and imaginary values are obtained for the 256 source-detector pairs, the data is written to a file for image reconstruction. The data aquistion time is longer than desired, but should be sufficient to start initial clinical trials for feasability measurements. In the future the data aquisition time can be reduced considerably by the addition of more photomultiplier tubes for parallel detection. At present all fibers are held in a plastic block holder (see Fig. 2), and can be moved manually along the radial direction to adjust to different object sizes. The fibers need to be indirect contact with the object being imaged, however feasability studies are being done to see if Intralipid filling the space between fiber and tissue can be used to help coupling to breast tissue in the clinical tests.
2.3 System performance
The performance of the system is dependent upon the characteristics of the modulated light source and detector. In particular, a high DC intensity and a high modulation ratio are important in order to maintain high signal to noise ratio . The DC intensity can be maintained at a high level by simply collecting enough light to achieve a good relative signal, however, the modulation ratio is dependent upon the system characteristics.
Measurements of AC amplitude and phase error (phase error is in degrees) are plotted in Fig. 3(a), as a function of AC amplitude, which demonstrates that in this situation an AC intensity of 0.1V and an AC error of 0.001V, results in a S/N ratio of 100, while decreasing the AC intensity to 0.01V results in S/N=1. Note here that 0.1 V corresponds to 10% modulation of the signal since the DC intensity was fixed at 1.0 V, which is easily achieved in imaging through 8 cm tissue phantoms. From these calibration measurements, it is clear that maintaining a high AC amplitude at the detector is crucial for reliable data since the S/N ratio scales with the square of the AC amplitude. The S/N response is a point of ongoing study to better understand the limitations of the system.
The frequency response of the source-detector system is plotted in Fig. 3(b), demonstrating that the optimal frequency range of operation is between 50 and 250 MHz where the modulation ratio is approximately 40%. It is interesting to note that if a S/N ratio of 1 is taken as the useful limit, then reliable measurements can be made with this system up to modulation frequencies of 800 MHz, and a S/N ratio of 100 can be maintained up to 500 MHz. Although it is important to appreciate that at these high frequencies the AC light signal is highly attenuated in human tissue, and may not be that useful for imaging.
Another important parameter to examine is the reproducibility for successive measurements of an test object. The two main factors which determine this error are the drift of the detectors, and the reproducability of the fiber optic multiplexer. Both of these add to the systematic error of the system, however in practice we have found that the drift of the photomultiplier tubes is the limiting error presently. The translation stages which are used to multiplex the fiber optics have a repeatabilty of approximately 10 microns, allowing very small variation in the measured signals. Alternatively the drift in the detectors are difficult to quantify exactly but we have found that during the course of repeated measurements on the same phantom, that the repeatability is approximately 0.5 degrees in phase and 0.5% in AC intensity. This measurement is complicated by the fact that not all measurements have the same AC or DC intensities, and that the data aquisition time is flexible. It is important to note that these errors are systematic and will manifest themselves in a small shifts of the reconstructed optical properties, rather than added noise in the image.
3. Theoretical tools to analyze photon migration in tissue
We have developed several numerical tools based on theoretical descriptions of light propagation in tissue which have been used to further examine the nature of the data measured by our hardware system. In addition, the ability to reconstruct quantitatively accurate images depends upon being able to precisely match the measured data with theoretical calculations.
where instantaneous scattering is assumed. The radiance I(r⃗, s⃗, t) is a function of both the position in space as well as the direction of the radiated energy within a solid angle dΩ. The radiance depends on the interaction coefficient μt and the scattering phase function p(s⃗, s⃗'). The light source is described by ϵ(r⃗, s⃗, t) and c is the velocity of light in the medium. Though Eq. (1) provides a valid description for most problems related to optical tomography, there are only a few simple cases for which exact solutions can be found. Therefore, the diffusion approximation of transport theory is commonly used instead .
While transport theory serves as a reference, the diffusion approximation is more attractive from an engineering point of view, because numerical calculations can be performed with better efficiency. This becomes particularly important for tomography, where iterative methods are used and light propagation has to be simulated multiple times to solve the inverse problem. In addition, the diffusion equation can be solved analytically for various geometries . Here we briefly summarize their properties, because they complement the experimental facilities described in section 2.
3.1 Computation of radiative transfer by Monte Carlo simulation
A very flexible way to calculate photon propagation in turbid media is Monte Carlo simulation. The probability laws of photon migration can be chosen freely which allows all assumptions of transport theory to be incorporated. Photon packages are traced within an absorbing and scattering medium and at every interaction the photon weight diminishes according to the ratio of absorption and scattering. The scattering characteristic is described by the Henyey - Greenstein equation , where we typically assume a value of 0.9 for the mean cosine of the scattering angle in our simulations.
The simulation algorithm we use for tracing the photon path is based on a program released previously . The original code was extended to allow the simulation of amplitude modulated light sources in the frequency domain. To this end we incorporated a similar approach to the one outlined in reference : the propagation distance along the photon path is accumulated. At the detector site the path lengths can be associated with a phase delay. Phase and photon weight define a complex value. The interference of all contributing photon packages then provides an estimate of the detected intensity.
The simulation program also allows the definition of different kinds of geometries, including a cylinder with absorbing boundaries (see Fig. 4) which corresponds to the phantom we use in our imaging experiments.
3.2 Solving the diffusion equation
The diffusion approximation turns the problem of photon migration into a simple boundary value problem. For a harmonically modulated light source we seek a solution of 
with boundary conditions of type
where n⃗ is the normal of the boundary surface. The solution is given in terms of the fluence ϕ(r⃗, ω) which depends on the modulation frequency ω, the absorption coefficient μa, and the diffusion coefficient D. The light source is given by S(r⃗, ω). The coefficient α describes the structure of the boundary.
To solve Eq. (2) numerically we use a finite element (FEM) algorithm. This algorithm calculates the forward problem of photon migration, within the diffusion approximation. It is also incorporated into the iterative reconstruction method which retrieves the inhomogeneous distributions of the absorption and scattering coefficient by minimizing the difference between simulated and measured intensity at the detector locations as well as the total-variation over the entire image. This minimization is performed by a Levenberg-Marquardt algorithm. Details of the finite element algorithm as well as of the reconstruction method can be found in references [10–13].
Fig. 5 shows calculations for a homogenous cylindrical phantom. On the one hand, these calculation prove, that our theoretical assumptions, incorporated into our Monte Carlo simulation meet the experimental conditions. On the other hand, they demonstrate, that the diffusion approximation provided by the FEM calculation is indeed a very good approximation for the problem with which we are concerned. This is even true at the boundary, where the detectors are located, although in general the diffusion approximation is only valid far from sources and boundaries . The measured data here fits the FEM calculation better than the Monte Carlo simulation, however the differences in both cases are small and is likely attributable to systematic problems with matching the boundary conditions of our tissue-simulating phantom.
4. Image reconstruction
The FEM routine for image reconstruction is called from within LABVIEW and processes the data into images of absorption and scattering coefficients. The reconstruction program works to minimize both the least-squares difference between the measured data and the calculated data from a 2-dimensional FEM of the diffusion equation, as well as minimizing the total variation of the image [10–13]. In practice the data fitted is the real and imaginary components of the AC signal, which is equivalent to fitting to the AC magnitude and phase. In these measurements, there was an inter-fiber intensity variation which could corrupt the data. For this reason, all of the following reconstructions have used data measurements which have been normalized to the homogenous sample. In practice, normalized data is created by precalibration measurements on a homogenous phantom which closely mimics breast tissue. In the future, this may also be achieved by imaging the same object at two different wavelengths , and normalizing the measurements. The images shown below were taken with an early version of the system which employed manual multiplexing of the source and detector. The computational time per iteration was approximately 2 minutes, but this can be decreased to approximately 1 minute per iteration through the use of a dual mesh reconstruction algorithm . A total of 10 iterations were obtained for each of the images reported here. A minimization method described previously  was used with total-variation and regularization parameters of 10-7 and 10, respectively.
Fig. 6 shows the reconstruction from a cylindrical 25 mm diameter high contrast absorbing inhomogeneity placed in the center of a 86 mm diameter scattering phantom. The background optical properties were D = 0.67 mm and μa = 0.0062 mm-1, with an inhomogeneity having D = 0.67 mm and μa = 0.056 mm-1. In this case the reconstruction is quantitatively correct for both absorption and diffusion coefficient, and the full width at half maximum for the inhomogeneity is 22 mm. Fig. 7 shows the reconstruction from the same size object as in Fig. 6, but with optical properties of D=0.2 mm and μa = 0.0062 mm-1. In this case the lower diffusion coefficient object is correctly reconstructed, but without quantitatively accurate results.
In the past 7 or 8 years, there have been many reports published discussing the possibilities of near infrared tissue imaging, yet very few papers have addressed the technical challenges of developing an integrated system for this purpose [28, 29]. One of the obvious applications for near infrared computed tomography is breast cancer detection, mainly because recent developments in diffuse computed tomography may afford better means of tissue discrimination over older clinical devices. Unfortunately, better methods for image reconstruction can also be more sensitive to measurement error, and it is crucial to design an optical system in which both the random errors of the measurements and the systematic mismatch between measurement and theory are all minimized. In this paper, these errors are characterized for our prototype system, and preliminary images of tissue phantoms have been presented.
Systematic errors in optical signal detection need to be eliminated by careful normalization of the data or by some other means (i.e. multiple wavelengths, frequencies, etc.). In this multiple source-detector arrangement variations in transmitted intensity between fibers ranged up to 30%. This level of systematic intensity variation must be corrected in the measurements before the data can be used in the reconstruction algorithm. In the cases presented here the heterogeneous phantom image data was normalized by data collected from a homogenous phantom. At this point, it is still not clear whether the limiting noise in the system will be systematic or statistical, although this depends upon the signal from each source-detector pair. At detector locations near the source, systematic error tends to impose the limitation on precision, while at detectors far from the source, where there is not much light signal, statistical errors from the detector dominate. In spite of these limitations, accurate absorption coefficient images can be created with the system for imaging phantoms which are similar to breast tissue [36–38].
The reconstruction of the absorbing object is quantitatively correct and is encouraging for images based upon changes in absorption (see Fig. 7). The second image of a diffusing object does not reconstruct properly, because there is a slight artifactual reconstruction of the diffusing effect into the absorption coefficient image (i.e. the presence of a pure absorbing and diffusing inhomogeneity cannot be exactly separated given the measurements obtained with this instrument, see Fig. 7). These results suggest that measurements are required which better discriminate between absorption and scattering, such as the use of higher frequency modulation . Another possibility might be to increase the total number of independent measurements through the use of more source-detector pairs. Alternatively, biasing the regularization method to better fit the phase shift can preferentially enhance the contrast of a diffusing object over an absorbing object , and this type of approach may turn out to be the most practical. While the inability to quantitatively reconstruct diffusion coefficient changes is problematic, most changes between normal and tumor tissue are thought to be predominately from blood volume [26, 36], and hence changes in absorption would dominate in this case. It is important to note that while the images shown in this paper were taken with manual multiplexing of the source and detector fibers, that identical images can be obtained presently with the automated system. Also there has not been a conclusive test as to how well the reconstruction algorithm will work in highly heterogenous media such as real breast tissue, however this will be one of the major goals of the early clinical testing.
In summary, the system presented here is a complete description of a computed tomography system which will be used in clinical trials to evaluate the potential for breast cancer screening. The architecture of this system provides an automated and accurate method for data collection and image generation.
This work was supported by NIH grant RO1CA69544 awarded by the National Cancer Institute. The authors would like to thank David Rinehart and Jennifer Beane for help in the development of the system.
References and links
1. M. Cutler, “Transillumination as an aid in the diagnosis of breast lesions,” Surg. Gyn. Obst. 48, 721–729 (1929).
2. D. J. Watmough, “Transillumination of breast tissues: factors governing optimal imaging of lesions,” Radiol. 147, 89–92 (1982).
3. R. J. Bartrum and H. C. Crow, “Transillumination light scanning to diagnose breast cancer: a feasibility study,” Am. J. Roentgenol. 142, 409–414 (1984).
4. A. Alveryd, I. Andersson, K. Aspegren, G. Balldin, N. Bjurstam, G. Edstrom, G. Fagerberg, U. Glas, O. Jarlman, and S. A. Larsson, et al., “Lightscanning versus mammography for the detection of breast cancer in screening and clinical practice. A Swedish multicenter study,” Cancer; Diag. Treat. Res. 65, 1671–1677 (1990).
7. P. C. Jackson, P. H. Stevens, J. H. Smith, D. Kear, H. Key, and P. N. T. Wells, “The development of a system for transillumination computed tomography,” Br. J. Radiol. 60, 375–380 (1987). [CrossRef] [PubMed]
8. J. R. Singer, F. A. Grunbaum, P. D. Kohn, and J. P. Zubelli, “Image reconstruction of the interior of bodies that diffuse radiation,” Science 24,: 990–993 (1990). [CrossRef]
9. S. R. Arridge, P. van der Zee, M. Cope, and D. T. Delpy “Reconstruction methods for infrared absorption imaging,” Proc. SPIE 1431, 204–215 (1991). [CrossRef]
10. 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). [CrossRef]
11. H. Jiang, K. D. Paulsen, U. L. Osterberg, and M. S. Patterson, “Frequency-domain optical image reconstruction for breast imaging: initial evaluation in multi-target tissue-like phantoms,” Med. Phys. (in press):, 1997. [PubMed]
14. S.R. Arridge and M. Schweiger, “Image reconstruction in optical tomography,” Philos. Trans. R. Soc. London Ser. B. 352, 717–726 (1997). [CrossRef]
16. S. T. Flock, M. S. Patterson, B.C. Wilson, and D. R. Wyman “Monte Carlo modeling of light propagation in highly scattering tissues -I: Model predictions and comparison with diffusion theory,” IEEE trans. Biomed. Eng. 36, 1162–1168 (1989). [CrossRef] [PubMed]
17. R. F. Bonner, R. Nossal, and S. Havlin, “Model for photon migration in turbid biological media,” J. Opt. Soc. Am. A. 4, 423–432 (1988). [CrossRef]
19. M. S. Patterson, B. C. Wilson, and D. R. Wyman, “The propagation of optical radiation in tissue: I. Models of radiation transport and their application,” Lasers Med. Sci. 6, 155–168 (1992). [CrossRef]
20. M. A. O'Leary, D. A. Boas, B. Chance, and A. G. Yodh, “Experimental images of heterogeneous turbid media by frequency-domain diffusing-photon tomography,” Opt. Lett. 20, 426–428 (1995). [CrossRef] [PubMed]
21. S. Walker, S. Fantini, and E. Gratton, “Image reconstruction by backprojection from frequency-domain optical measurements in highly scattering media,” Appl. Opt. 36, 170–179 (1997). [CrossRef] [PubMed]
22. M. S. Patterson, J. D. Moulton, B. C. Wilson, K. W. Berndt, and J. R. Lakowicz, “Frequency-domain reflectance for the determination of the scattering and absorption properties of tissue,” Appl. Opt. 30, 4474–4476 (1991). [CrossRef] [PubMed]
23. J. Fishkin, E. Gratton, M. J. van de Ven, and W. W. Mantulin, “Diffusion of intensity modulated near infrared light in turbid media,” Proc. SPIE 1431, 122–135 (1991). [CrossRef]
26. B. J. Tromberg, O. Coquoz, J. B. Fishkin, T. Pham, E. R. Anderson, J. Butler, M. Cahn, J. D. Gross, V. Venugopalan, and D. Pham “Non-invasive measurements of breast tissue optical properties using frequency-domain photon migration,” Philos. Trans. R. Soc. London Ser. B. 352, 661–668 (1997). [CrossRef]
28. J. Hoogenraad, J. M. van der Mark, S. B. Colak, G. W. t'Hooft, and E. S. van der Linden, “First results of the Phillips optical mammoscope.“ Proc. SPIE 3194 (in press).
29. 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 techniques enhance optical mammography: initial clinical results.“ Proc. Nat. Acad. Sci. USA 94(12), 6468–73 (1997). [CrossRef] [PubMed]
30. V. V. Sobolev, A Treatise on Radiative Transfer, (Van Nostrand-Reinhold, Princeton, 1963), pp. 240–244.
31. A. Ishimaru, Wave Propagation and Scattering in Random Media, (Academic Press, New York, 1978), Vol. 1, pp. 147–167 and pp.175–190.
32. S. R. Arridge, M. Cope, and D. T. Delpy, “The theoretical basis for the determination of optical pathlengths in tissue: temporal and frequency analysis,” Phys. Med. Biol. 37, 1531–1560 (1992). [CrossRef] [PubMed]
33. L. Wang and S. L. Jacques, Monte Carlo Modeling of Light Transport in Multi-layered Tissues in Standard C, (University of Texas M. D. Anderson Cancer Center, Houston, 1992-93).
34. I. V. Yaroslavsky, A. N. Yaroslavsky, V. V. Tuchin, and H.-J. Schwarzmaier, “Effect of scattering delay on time-dependent photon migration in turbid media,” Appl. Opt. 36, 6529–6538 (1997). [CrossRef]
35. B. W. Pogue, M. S. Patterson, and T. Farrell, “Forward and inverse calculations for 3-D frequency-domain diffuse optical tomography,” Proc. SPIE 2389, 328–339 (1995) [CrossRef]
36. V. G. Peters, D. R. Wyman, M. S. Patterson, and G. L. Frank, “Optical properties of normal and diseased breast tissue in the visible and near infrared,” Phys Med. Biol. 35, 1317–1334, (1990). [CrossRef] [PubMed]
37. T. L. Troy, D. L. Page, and E. M. Sevick-Muraca, “Optical properties of normal and diseased breast tissues: prognosis for optical mammography,” J. Biomed. Opt. 1, 342–355 (1996). [CrossRef] [PubMed]
38. K. Suzuki, Y. Yamashita, K. Ohta, M. Kaneko, M. Yoshida, and B. Chance, “Quantitative measurement of optical parameters in normal breasts using time-resolved spectroscopy: In vivo results of 30 Japanese women,” J. Biomed. Opt. 1(3), 330–334 (1996). [CrossRef] [PubMed]