1. Introduction
High speed near-infrared spectral tomography of tissue can be used to measure functional properties of tissues such as metabolism, hemodynamics and contrast agents to provide useful diagnostic information. Imaging at a speed much faster than the tissue activity rate of interest is important to provide quantitative information with high signal to noise ratio (SNR), and potentially to allow signal processing methods to be used which can amplify low SNR. But practically building a fast scanning tomography system for biomedical application remains a challenging task, as the diffuse pattern of light propagation in biological tissue demands that each source-detector (S-D) pair be encoded properly through delivery and decoded afterwards. Accumulating enough photons at the detectors in a limited exposure time is difficult, especially when the tissue is thicker, e.g. more than 5–6 cm for imaging most critical organs or tissue volumes. Methods for mechanical switching of S-D pairs have been routinely used, with imaging speeds reported up to 3 Hz [
1
C. H. Schmitz, M. Löcker, J. M. Lasker, A. H. Hielscher, and R. L. Barbour, “Instrumentation for fast functional optical tomography,” Rev. Sci. Instrum.
73(2), 429 (2002).
[CrossRef]
]. Simultaneous source delivery with frequency encoding [
2
S. Fantini, M. A. Franceschini, G. Gaida, E. Gratton, H. Jess, W. W. Mantulin, K. T. Moesta, P. M. Schlag, and M. Kaschke, “Frequency-domain optical mammography: edge effect corrections,” Med. Phys.
23(1), 149–157 (1996).
[CrossRef]
] can be done considerably faster, but is limited by the overwhelming effect of closer strong signals dominating over the farther weaker signals that propagate through a few cm of tissue. This approach, while tremendously successful in sub-surface imaging [
3A. M. Siegel, J. J. A. Marota, and D. A. Boas, “Design and evaluation of a continuous-wave diffuse optical tomography system,” Opt. Express 4(8), 287–298 (1999).
], cannot be realized in deep tissue tomography, because of the severe dynamic range variation seen in signals that propagate different distances in tissue.
Our earlier work reported a configuration with spectrally encoded sources around a circular coil, where the spectral encoding was introduced as a way to eliminate signal interference at the detector [
4
D. Piao, H. Dehghani, S. Jiang, S. Srinivasan, and B. W. Pogue, “Instrumentation for video-rate near-infrared diffuse optical tomography,” Rev. Sci. Instrum.
76(12), 124301 (2005).
[CrossRef]
]. The sources were identical LDs with different temperature-controlled wavelength shifts, and using a single imaging spectrometer detection at 30 frames per second was achieved with this setup. The design demonstrated the concept, with limited dynamic range capability because of the single CCD detector and the circular geometry of the coil. It worked for small circular-shaped tissues with a radius of less than 3 cm, but a revised design was needed for larger tissue volumes. In this study, a parallel detector design is implemented which demonstrates the ability to image through much thicker tissues and therefore has applications in breast tissue or brain tissue tomography.
Fig. 1. A diagram of the video-rate near infrared tomography system is shown. Eight spectrally encoded LD systems were integrated onto one cart [5]. Eight fiber-coupled spectrometers with high-resolution CCDs were integrated onto another cart and set to respond to an external TTL trigger signal. The signal was generated by a data acquisition board and then split into eight channels by a customized splitter circuit, and delivered to the EXT SYNC port on each CCD simultaneously. Imaging data are transferred back to the computer for post processing. A flexible fiber mount was customized to hold the phantom or tissue in between.
The new system design has similar spectral encoded LD sources but is used with a slab-shaped interface and multiple spectrometers for detection. Eight spectrometers that are synchronized through external TTL trigger signals were used, providing much better signal to noise ratio than the previous effort [
4
D. Piao, H. Dehghani, S. Jiang, S. Srinivasan, and B. W. Pogue, “Instrumentation for video-rate near-infrared diffuse optical tomography,” Rev. Sci. Instrum.
76(12), 124301 (2005).
[CrossRef]
]. The use of parallel plate slab geometry is also important, because the signals intensities that propagate from source to detector locations are then all closer in value, allowing a more regular dynamic range. The system was constructed and tested in tissue phantoms, and the performance is characterized here. The implications and planned use is discussed. The system has the capacity to image through perhaps up to 8–10 cm of tissue, which would make it possible to be used to track hemoydnamics in peripheral limbs, breast tissue or part of the cranium.
2. System Configuration
The video-rate tomography system [
5
Z. Li, V. Krishnaswamy, K. D. Paulsen, and B. W. Pogue, “Video-rate near infrared tomography for imaging thick tissue with dynamically varying absorption properties,” Proc. SPIE
7174, 71742G (2009).
[CrossRef]
], depicted in
Fig. 1, consists of several sub systems. Eight spectrally encoded laser diodes (LDs) each with a temperature control module and a current control module were integrated on to one cart. Eight fiber-coupled spectrometers, each equipped with a high-resolution charge-coupled device (CCD), were integrated on to another cart and set to respond to an external TTL trigger signal. The signal was generated by a data acquisition (DAQ) board from National Instrument (NI) and programmed through LabVIEW software. A customized splitter circuit was built to convey the signal into eight channels and finally to the EXT SYNC port on each CCD simultaneously. Imaging data are transferred back to the console for post processing through USB connection with the CCDs [
5
Z. Li, V. Krishnaswamy, K. D. Paulsen, and B. W. Pogue, “Video-rate near infrared tomography for imaging thick tissue with dynamically varying absorption properties,” Proc. SPIE
7174, 71742G (2009).
[CrossRef]
]. A customized flexible fiber mount and phantom/tissue interface was built to hold the phantom in the parallel plate slab geometry with the fibers in contact with it. Each component is described in detail below.
2.1 Spectrally encoded sources
In conventional NIR tomography systems, light is delivered to one plane at the same wavelength. Therefore sources are launched in sequence to avoid interference at the detectors. So the time needed for one image is proportional to the number of sources in use, unless temporal frequency encoding methods are used. The new video-rate tomography system uses 8 identical single mode continuous wave (CW) LDs (Hitachi HL-7851G) in the 785 nm band. They were spectrally encoded to have approximately 1.5 nm interval in between by shifting the temperature of each LD differently. This interval is wide enough for the 1200 l/mm grating on the CCD to separate each source signal clearly, but still very narrow so that the tissue optical properties remains almost the same. The difference of the absorption coefficient in this range for whole blood is within 5%. Therefore, all source LDs can be launched simultaneously to improve imaging speed without interference with each other. The stability of the LDs under wavelength shift has been verified by sampling over long periods of time [
5
Z. Li, V. Krishnaswamy, K. D. Paulsen, and B. W. Pogue, “Video-rate near infrared tomography for imaging thick tissue with dynamically varying absorption properties,” Proc. SPIE
7174, 71742G (2009).
[CrossRef]
].
2.2 Optical Detectors
The optical detection system consisted of 8 Princeton/Acton Insight: 400F Integrated Spectroscopy Systems (Acton, MA) residing in one custom designed wheeled carts (8020, Columbia City, IN) [
6
S. C. Davis, B. W. Pogue, R. Springett, C. Leussler, P. Mazurkewitz, S. B. Tuttle, S. L. Gibbs-Strauss, S. S. Jiang, H. Dehghani, and K. D. Paulsen, “Magnetic resonance-coupled fluorescence tomography scanner for molecular imaging of tissue,” Rev. Sci. Instrum.
79(6), 064302 (2008).
[CrossRef]
]. Insight 400F consists of a 0.3 m length, F3.9 imaging spectrograph and a low noise, front illuminated charge coupled device (CCD) (Pixis 400F) cooled to −70 °C. The 16- bit CCD provided a large dynamic range of 4.8 orders and was vertically binned to maximize detector area/wavelength. A 1200 l/mm grating in each spectrometer provided a sufficiently fine spectral resolution, near 0.045 nm, which is much smaller than the 1.5 nm laser bandwidth interval, and a wide enough spectral range of 60 nm to allow simultaneous imaging of all wavelengths.
2.3 Tissue/Phantom Interface
Previous attempts at a NIR tomography system used a conventional circular coil design with source fibers and detector fibers interspersed around the coil [
4
D. Piao, H. Dehghani, S. Jiang, S. Srinivasan, and B. W. Pogue, “Instrumentation for video-rate near-infrared diffuse optical tomography,” Rev. Sci. Instrum.
76(12), 124301 (2005).
[CrossRef]
]. This worked well for fast small object imaging where the optical dynamic range was not larger than 3.0 OD between the different paths from each source to each detector. But for thick tissue imaging, e.g. clinical breast cancer imaging in vivo with diameters of up to 10 cm, it is technically impossible to manufacture reliable optical detectors with up to 8 orders of magnitude dynamic range. The problem is that even if the detectors are 16 bit CCDs, a signal that is 8 orders of magnitude smaller than another could not be detected. Meanwhile it is clear that the optical path length difference of the nearest and farthest detectors from the same source is much smaller in the slab geometry. Thus the slab configuration substantially lowers the requirement of dynamic range to closer to 3–4 OD range, and well fits for high speed thick tissue imaging. A flexible fiber holder, as shown on the right side of
Fig. 5, consisting of two parallel panels with two slide bars across below was built to hold the source and detector fibers against the two sides of the tissue/phantom. On each panel, 8 holes were drilled with 8 mm intervals across to hold source or detector fibers, and set screws from the top fixed the fibers in position. This configuration is flexible for different phantom thicknesses and stable enough for the experiment, as well as for future use in tissue imaging.
2.4 Exposure Synchronization
To image with fast and continuous CCD exposures was the goal of this work. The Pixis 400F provides a good signal to noise ratio, but uses a slow mechanical shutter in front of the CCD to control input light, so the shutter was kept open during these experiments. Moreover, all 8 CCDs were synchronized to start and stop exposing at exactly the same time to get meaningful timeline data set. To achieve this, an external TTL trigger signal was generated by a NI DAQ board, and split into 8 channels by a customized splitter circuit, and finally reached the external sync port of each CCD to accurately synchronize all CCDs. The shape and pulse of the TTL trigger signal can be easily specified through LabVIEW. But for the Pixis 400F, there is no standard mode to work with multiple CCDs both fast enough and synchronically though the external trigger port is provided. Using the Scientific Imaging Tool Kit (SITK) for LabVIEW, all 8 Pixis 400F CCDs were configured to work in a “non-standard” high speed acquisition mode to solve this problem. The shutter status on each CCD is set to open before an acquisition sequence and keep open until the sequence is finished. The exposure mode was set to Strobed Mode, meaning each exposure requires an external trigger signal to start. Each CCD was set to Synchronous Focus Mode. This mode alone would make one CCD keep exposing at a preset rate and only transfer the latest frame when the console was ready to receive, so no synchronization could be realized on multiple CCDs. But it was discovered that together with the Strobed Mode exposure setting, the CCD actually kept exposing until the next trigger signal arrived, and then it transferred the data and cleaned the built up charges and started a new exposure, which perfectly fit the design requirement. Through repeated measurements, it was proven that the whole transfer-clean-prepare process for all 8 CCDs took around 13 milliseconds. Since the shutter remained open, the actual exposure time was the trigger interval, which was not less than the nominal preset CCD exposure time, plus 13 milliseconds. With 100 µm slit-width on all spectrographs and 10 mW input laser beams, 20 data sets per second could be reached easily on a 64 mm thick resin phantom with satisfactory signal intensity.
3.3 Pulsatile Phantom Experiment
The clinical goal of this video-rate NIR tomography system is to image tissue hemodynamics. Thus a pulsatile phantom experiment was carried out to test the ability of this system to image an abnormity in real-time whose absorption property varied at a preset frequency with moderate effect upon the transmitted signal. The experimental setup is illustrated in
Fig. 5. A solution made of 1% intralipid and low concentration of Indian ink was filled into a 72 mm thick slab-shaped high-density polypropylene bottle (32 oz HDPE Rectangular Bottle, Nalgene) to provide a homogeneous background with low µ
a. The plastic bottle was colorless with a cloudy appearance that prevented straight light channeling laterally, so it was a fairly good container for this experiment. Several sides of the bottle were sealed with black tape to minimize the impact of boundary mismatch. Another tank held solution made of the same 1% intralipid and 3 times higher concentration of Indian ink than that in the Nalgene bottle. So the absorption contrast of two solutions was 3:1. The high µ
a solution was continuously pumped from the tank into an incoming latex tubing with an internal diameter (ID) of ¼ inch. This tubing was connected to one end of a T-shape connector. The second end of the T-shape connector was connected to an outgoing clear vinyl tube with the same ID and directed the high absorbing solution back to the tank. A white balloon was tied up to the third end of the connector and sealed with water-proof glue. The connector was fixed on a three-jaw chuck with the balloon submerged into the low µ
a solution in the Nalgene bottle. The incoming tubing was periodically compressed by the rotational part of a pump (503U/RL, Watson-Marlow L.L.C.) at a preset frequency. The pump was driven by a frequency-adjustable 0~10 Volt square wave voltage generator. Therefore the high µ
a solution was pumped through the T-shape connector in pulsatile mode at the same frequency, and the balloon expanded and contracted at the driving frequency.
Fig. 5. System setup of the pulsatile phantom experiment. High µa solution was continuously pumped through a balloon at 0.5 Hz. The balloon was submerged in a 72 mm thick slab container filled with low µa solution. The absorption contrast of the solution inside the balloon against outside was 3 to 1. Eight LDs were launched as sources, and 6 spectrometer systems were set as detectors to acquire data at 10 frames per second. Slit widths on all spectrographs were set to 100 µm.
Fig. 6. Signals of different S-D pairs in frequency domain.
Figure 6(a) is detected signals of 3 different optical paths in frequency domain after being normalized to their individual mean intensity. Only 0 to 1 Hz of these spectra were displayed.
Figure 6(b) illustrates the 3 different optical paths through the phantom with straight lines. D1 to D6 were 6 detection spots connected to spectrographs through optical fibers. S1 to S8 were 8 laser beam input spots connected to LDs through optical fibers.
This pumping system was designed to mimic the cyclical change in absorption seen in vaso-activity within tissue. The video-rate tomography system was used to image and recover this pulsatile change of the absorbing content, and quantify the magnitude of the effect upon the transmitted signal. The experiment was designed such that the elasticity of the balloon would only allow response below pump frequencies of 1 Hz. So a 0.5 Hz drive was chosen as the experimental pumping frequency for this study. Most pumped solution went through the T-shape connector directly, and only a small amount of the solution was to be pushed into the balloon corresponding to the change of pressure difference of the incoming and outgoing tubes. The volume of high µa solution in the balloon and thus the size of the balloon would vary according to the pumping frequency, but only on a very subtle scale which was barely apparent and not significantly measurable. So the variation of µa in the balloon was expected to be very small, and was quantified by the effect it had on the transmission of the light.
The detected signals for 3 different optical paths in frequency domain (FD) are displayed in
Fig. 6. After subtracting the mean value of each signal to remove DC component, a Fourier transform was performed on each signal and then the signal was normalized to the intensity of DC component separately, such that the change is relative to 1.0. The relative intensity of every frequency component on the FD spectrum reflected the portion of signal caused by the varying portion of absorbing property versus that by the homogeneous background. It is obvious from
Fig. 6(a) that path 1, which was the shortest and went through the balloon, received strongest 0.5 Hz pulsatile information, increasing from 0.005 up to 0.015, or 0.5% to 1.5%. Meanwhile, path 3 received weakest 0.5 Hz pulsatile signal since it was the longest, although it also went through the balloon. The largest effect of the pulsatile flow was set at approximately 1% change in the DC signal, as can be seen in the peak of
Figs. 6(a) and
6(b). This low magnitude in pulsatile flow is approximately what is seen in pulse oximetry, and so was chosen as a good example of typical physiological signal level.
After imaging the pulsatile phantom at 10 frames per second for 20 seconds, the raw data set was calibrated and put into the reconstruction program. Mean µ
a of the ROI using direct reconstruction algorithm versus time is shown in
Fig. 7(a). There was a blurry periodical trend in curve µ
a, but the noise made it not easy to distinguish. A Fourier transform was applied to the µ
a curve, and the result was shown in
Fig. 7(b). Since the experimental imaging speed was 10 frames per second, the highest distinguishable frequency after Fourier transform was 5 Hz. There was a peak right at 0.5 Hz, but several other peaks with similar intensity at higher frequencies made it not convincing to conclude the system had picked out the desired pulse. The result of difference reconstruction was not listed here, because the reconstructed µ
a difference was very small and could not be clearly discerned from reconstruction artifacts around.
Fig. 7. Reconstructed µ
a of the ROI in the time-domain and frequency-domain. In
Fig. 7(a) the mean recovered µ
a value in the ROI is plotted versus time with a direct diffuse reconstruction method. In
Fig. 7(b) the Fourier transform of the mean µ
a in
Fig. 7(a) is shown. In
Fig. 7(c) the mean µ
a recovered in the ROI versus time is shown using region-guided direct reconstruction methods.
Figure 7(d) is the Fourier transformed amplitude data from
Fig. 7(c), showing the dominance of the 0.5 Hz signal in this frequency spectrum.
A region-based direct reconstruction method was also applied to this data set. Classically the imaging plane can be segmented into several homogeneous regions based on the structural information obtained in advance with other imaging modalities, e.g. MRI and/or x-ray CT. The iterative update is then simplified from the situation where individual nodes are reconstructed to the problem where only the parameters for each homogenous region are estimated in the inversion [
10
B. Brooksby, H. Dehghani, B. W. Pogue, and K. D. Paulsen, “Near infrared (NIR) tomography breast image reconstruction with a priori structural information from MRI: algorithm development for reconstructing heterogeneities,” IEEE J. STQE
9, 199–209 (2003).
]. In this case, instead of inverting a matrix that is thousands of rows, it has just 2 rows, one for each material. The advantage of this method is that the contrast is enhanced between the inclusion and the background. The two regions defined in this experiment were from direct measurements of the phantom size and shape and location of the balloon. Since the balloon size change was too small to measure between the pulsatile motion, this was assumed to be constant. In
Fig. 7(c) the recovered µ
a of the ROI versus time showed a clear periodical trend, and after Fourier transformation a sharp peak at 0.5 Hz was recovered in
Fig. 7(d), which is much higher than the background. This peak verifies that this region-based method does allow recovery of the pulsatile varying absorption content in the balloon, likely due to noise suppression through the combined effect of using all the measurements together. As such the inversion matrix problem is highly over-conditioned, and so is a well posed problem.
4. Discussion
This work has shown the feasibility of using this spectrally-encoded tomography system to do fast speed imaging through thick tissues, of the size required in clinical breast imaging [
7
S. Jiang, B. W. Pogue, S. Srinivasan, S. Soho, S. P. Poplack, T. D. Tosteson, and K. D. Paulsen, “Assessment of the menstrual cycle upon total hemoglobin, water concentration and oxygen saturation in the female breast,” Proc. SPIE
4955, 342–348 (2003).
[CrossRef]
–
10
B. Brooksby, H. Dehghani, B. W. Pogue, and K. D. Paulsen, “Near infrared (NIR) tomography breast image reconstruction with a priori structural information from MRI: algorithm development for reconstructing heterogeneities,” IEEE J. STQE
9, 199–209 (2003).
]. The goal of the system was to achieve high speed imaging, of blood pulsation dynamics, and show the validation that the system could be used to recover subtle signals in realistic tissue phantoms. The system was built around CCD detection which allows only CW detection, and does not have the capability to image frequency domain data. This choice was intentionally made because of the high cost that PMT or APD detectors would incur in developing high-frame rate systems with a similar dynamic range. The tradeoff is the lack of the ability to accurately recover exact absorption and scattering coefficients independently from one another, because there is no clear ability to estimate photon traveled pathlength in tissue. However the application of the system has focused on vascular pulsatile flow imaging, and in this case the scatter coefficient remains largely static. Thus the errors in assumed scattering coefficient will mostly contribute to an overall offset of recovered absorption coefficient. Over the years, changes in absorption have been repeatedly been shown to be quantifiable with CW systems, such that this design was chosen for this application. The errors in absolute absorption from the unknown value of the scattering coefficient are a problem, but the error incurred is low if we have a reasonable value of scattering coefficient to use in the model fitting. If the estimated scattering coefficient is within 10% of the true value, then our error in absorption is less than this level. Still, changes in absorption are always recoverable at some level, even if the scatter signal is inaccurate by 20%, so the system does have highly robust response for dynamic imaging applications.
The use of spectrally-encoded sources in the tomography system has the implicit assumption that they are treated as interacting with the same tissue optical properties, even though they are at slightly different wavelengths. Strictly speaking, their wavelengths are different, but if treated as a multi-wavelength tomography system, and if we reconstructed the wavelengths separately, we would be able to potentially recover more information about the spectral features of the tissue. But this approach would have obviously fewer measurements at each individual wavelength, and would make a multispectral reconstruction possible, but highly ill-posed at any given wavelength. When assuming most of the absorption is due to whole blood, it can be seen that its absorption coefficient varies only by about 5% across the used wavelength range. Thus, by only incurring 5% error, we assume the attenuation is the same at all sources, and thereby solve the image reconstruction problem with a better posed inversion problem, and avoid the need to get into spectral reconstruction. However the idea of doing multispectral reconstruction at video rate is something to consider in future research in this area.
The dynamic blood diffusion phantom experiments illustrated the ability of the system to image heterogeneities through 64 mm of tissue-like media at high speed. The direct reconstruction method worked well, and so did difference reconstruction methods, since the absorption content in the inclusion changed much during the process. In the pulsatile phantom experiments, the system was used to illustrate the problem of accurately detecting very subtle periodic variations, on the order of 1% intensity signal change. The volume change in the phantom was approximately 0.5%, which resulted in the observed transmission signal change by 1%. The absorption properties of the inclusion, could not be usefully recovered with diffuse tomography imaging, as might be expected for a system where the signal change is only 1%, and the noise level of the system is near 1%. Thus, in a realistic tissue imaging situation, where SNR is close to unity, the application of diffuse tomography alone will be limited without improvements in signal amplification somehow. However, image-guided recovery of the absorption, using a region-fitting based method, did recover useful signal changes, and the amplitude of the periodic signal was found with a factor of 4 between the amplitude peak and the background noise at other frequencies (
Fig. 7(d)). Image-based recovery of broad regions turns the recovery problem from an ill-posed under-conditioned situation, to a well-posed over-conditioned inversion problem. Thus in the future, this type of imaging will likely be utilized in the context of MRI-guided hemodynamic imaging.
The repeatability and linearity of this system were verified in static blood phantom experiments, to demonstrate that the recovered values are linear with concentration of absorber, and that the accuracy can be consistent with our previous systems. There can be bias offsets of assumed scattering coefficients which prohibit absolute recovery of absorption coefficients, however that system is used with continuous wave data where the scattering coefficient is assumed anyway. So the goal of the system is not to do absolute tissue spectroscopy, instead the level of real-time absorption property change is of interest. The use of an assumed scattering coefficient is unfortunate, but earlier studies have found that reconstructing differential changes in tissue is feasible with assumed homogeneous scattering distributions. If the scattering coefficient is significantly erroneous, it can lead to significant errors in the reconstruction process, yet if it is within 10–20% of the true value, recovery of images is not significantly impeded.
Results of both direct and difference reconstruction methods were listed and compared. The approach of using difference-data based reconstruction generates images with fewer artifacts, especially at boundaries, but the ROI profile biases a little towards higher sensitivity region, and this is more of a problem towards the closest boundary. Besides, difference reconstruction doesn’t perform well when the difference of nearby images is subtle. So for applications like hemodynamics monitoring, direct reconstruction, maybe image-guided by other modalities, is more suitable.
The slab geometry was adopted here to decrease the differences of optical pathlengths of source-detector pairs, since dynamic range is such a problem in diffuse tomography. To have a very high imaging rate, all sources and detectors have to be kept on during the whole imaging process, and so the system needs to have a reasonably small dynamic range between individual source signals to make the signals useful. The slab geometry has a much smaller dynamic range than the circular geometry, and solves some of this problem. For tissue that is compressed into the slab shape of say 6–7 cm thickness, the 8 sources are placed at one side with 8 mm interval, and 8 detectors are oppositely placed in the similar way. Then the maximum path length difference is near 2 cm, which is much smaller than the circular geometry design, where differences in pathlengths can be on the same order as the diameter of the tissue. Thus, the slab approach solves this seemingly trivial but yet quite important issue of maintaining the dynamic ranges within a few orders of magnitude from each other between signals.
Real tissue imaging applications of this system are ongoing, using MRI imaging to visualize organs, and then rapid imaging with the spectrally encoded system to recover hemoglobin changes as well as exogenous dye uptake in the organs being tracked. The results of this study show that temporal changes in signal intensity at the level of 1% are still reliable enough to provide useful data for this image-guided tomographic approach. These results then allow estimation of how concentration changes of hemoglobin or dye per unit volume would translate into detectable levels of transmitted signal. The system testing in rodent tissue imaging is ongoing.