Optical tomography is a noninvasive functional medical imaging modality, the purpose of which is to obtain three-dimensional (3D) images of the optical properties of tissue at near-infrared (NIR) or visible red wavelengths. The propagation of light in tissue is stochastic in nature, and the optical properties of tissue are defined accordingly. The absorption and scattering coefficients ( and ) describe the frequency of absorption and scatter events per distance traversed by the photon. The absorption at different wavelengths is of particular interest, as it depends strongly on the concentrations of oxygenated and deoxygenated hemoglobin ( and Hb). Changes in and Hb with respect to an arbitrary baseline can be determined using a time series of the intensity or amplitude data. If the absolute optical properties of the tissue can be determined, it is possible to estimate the oxygen saturation and blood volume of the tissue.
In optical tomography, several source and detector fibers are placed on the surface of the volume of tissue to be imaged.1. In most implementations, the source fibers are sequentially activated, and light is detected in parallel by all measurement channels. The optical properties of the tissue are then derived using an appropriate image reconstruction algorithm.2 Three basic instrument types are widely used in research. Intensity measurements can be made with a relatively simple instrument, based on large-area detectors with high quantum efficiency (QE) and therefore allowing the rapid acquisition of data. However, it has been proven that in order to obtain simultaneous reconstructions of scatter and absorption, intensity measurements alone are insufficient unless appropriate prior (e.g., spectral) information is available.3 Radiofrequency intensity-modulated systems4, 5 allow the measurement of the phase shift and amplitude attenuation of the photon density wave between the source and detector positions. Usually such systems employ only one modulation frequency to simplify the technical implementation. Intensity-modulated diffuse optical imaging was first demonstrated by Gratton 6 and O’Leary 7 The review of phase and amplitude measurement systems by Chance 4 includes detailed descriptions of the sources of error in such systems. Time-resolved systems8, 9 illuminate the tissue with picosecond pulses of light and measure the temporal distribution of the detected photons. By deconvolving the temporal response of the system, it is possible to recover the so-called temporal point spread function (TPSF) of the tissue.10 Using a subsequent Fourier transform, the TPSF can be converted into the frequency domain. This conversion is useful because the frequency-domain diffusion equation and the corresponding data types, amplitude and phase, allow a simpler implementation of image reconstruction algorithms than using temporal data types directly.2 They also provide relatively good sensitivity to the inner optical properties of the tissue. By using only data corresponding to a single modulation frequency, the full information content of the time-domain system is not used. The additional information available in higher-order moments of the TPSF may improve image quality over what is presented here, but the improvement is not dramatic in our experience. The time-domain method has the additional advantage that the early part of the TPSF can be filtered out, which reduces the effects of optical cross talk due to light leaks in the fiber-tissue interface and those due to the internal cross talk of the switch. In this paper, we reduced optical cross talk by lining the fiber-holding rings with black velvet. In nontomographic applications, such as cortical activation mapping, the late photons in a time-domain system can be used to improve the sensitivity of the measurements to deeper tissues.11
The technical implementations of time- and frequency-domain systems are quite different.1, 12 Time-domain systems require picosecond-pulsed lasers and appropriate ultrafast detectors, which are coupled to electronics able to record the arrival times of individual photons. Frequency-domain systems use continuous-wave light sources and modulate the intensity of light with a radiofrequency signal (typically 50 to ). The amplitude and phase of the detected signal are usually measured with a lock-in amplifier.
In deciding which system type is most appropriate for a given application, it is necessary to consider the system cost, the ease of use (particularly in a hospital environment), and the quality and variety of data that can be obtained. In this paper we compare two systems, a time-domain system developed at University College of London (UCL),13 and a frequency-domain system developed at Heklsinki University of Technology (HUT).14 The data are compared in the frequency domain. We study the quality of the data and performance of the systems using the difference imaging approach. Two solid cylindrical phantoms were constructed for this purpose: one homogeneous and a second that contains two small objects representing perturbations in either absorption or scatter. In this paper, “difference data” is used to mean the differences in the measured phase and logarithmic amplitude data between measurements made on the inhomogeneous and homogeneous phantoms. The term “absolute data” is used to mean the attenuation of the amplitude and the phase shift that takes place between the source and detector in the phantom. While it is possible to do absolute imaging using calibrated data without a homogeneous reference phantom, difference imaging is often used as it has better sensitivity to small changes in the optical properties and it produces images with fewer artifacts provided that a geometrically accurate reference phantom exists. Because imaging time, which is the time required to collect a complete data set using all sources and detectors, is often a vital criterion in functional medical imaging applications, the data were acquired using a series of measurement times.
The measured difference data are compared with simulated data generated by a 3D finite-element method (FEM) solver to the diffusion approximation (DA) of the radiative transfer equation (RTE). To validate the FEM model, a Monte Carlo (MC) simulation program was developed, and simulated data from the FEM and MC models were compared both in terms of absolute data using the homogeneous model parameters as well as the difference data generated by subtracting the homogeneous data from the inhomogeneous. The FEM simulated difference data and the measured data were compared by calculating the rms difference between them for each measurement time.
The measured difference data sets were used to generate images of the absorption and scattering coefficients in the inhomogeneous phantom. The images were reconstructed using a regularized Gauss-Newton inverse algorithm.15 The effects of measurement time on the image quality are discussed for both systems.
UCL Time-domain optical tomography system
The UCL optical tomography system (Fig. 1 ) consists of 32 parallel detectors that measure the times of flight of transmitted photons at two NIR wavelengths simultaneously. The sources are two synchronized fiber lasers (built by IMRA, Inc.) operating at 780 and , which illuminate the tissue via a 32-way optical fiber switch. Light transmitted across the tissue is collected by 32 detector fiber bundles, which are coupled to 4 eight-anode microchannel-plate photomultiplier tubes (MCP-PMTs). The output end of each source fiber (graded index, diameter ) is permanently integrated along the central axis of a corresponding detector bundle (approximately 3000 step index fibers with a core diameter of and a cladding diameter of ; 3.2-mm bundle diameter). This coaxial arrangement has two major benefits. First, it decreases the number of connectors required to couple 32 sources and 32 detectors to the tissue. Second, it enables back-reflected light at the surface to be used to calibrate the system in situ, as described in Ref. 16. The detection electronics consists of 32 parallel time-correlated single photon counting (TCSPC) units. A detailed description of the imaging system, known as MONSTIR (Multi-channel Opto-electronic Near-infrared System for Time-resolved Image Reconstruction), is provided in Ref. 13. Each channel can simultaneously detect up to around photons per second. Light is coupled to the MCP-PMTs via 32 computer-controlled variable optical attenuators (VOAs), which provide a wavelength-independent attenuation up to about 3 optical densities. This increases the dynamic range of the system and ensures that the intensity of detected light does not saturate or damage the MCP-PMTs. Arrival times of detected photons are measured with respect to a laser-generated reference signal, and histograms of photon flight times (TPSFs) are accumulated. The full set of TPSFs is subsequently transferred to a dedicated workstation for processing.
The principal clinical application for which the UCL system was designed is 3D imaging of the newborn infant brain. Due to the presence of hair and uncertainty in the skin reflectance, the efficiency of the coupling of light into and out of the infant head is generally almost impossible to quantify. As a consequence, absolute measurements of integrated transmitted intensity are generally unusable for image reconstruction of the infant head; although, in principle, differences in integrated intensity acquired on the same subject can be used in the absence of motion-induced variability in coupling. In acknowledgment of this problem, the UCL system was not originally designed to provide quantitative (or indeed highly repeatable) measurements of intensity. However, recent careful calibration for the variable losses in each component has enabled useful difference intensity data to be obtained for infant brain imaging.17
Frequency-domain system developed at HUT
The source fibers are time multiplexed using a fiber-optic switch (DiCon VX500-16) and one of the two wavelengths is selected using a DiCon prism switch. The signals from the PMTs are mixed to an intermediate frequency of , and the phase and amplitude are calculated in software from signals digitized using two synchronized PCI-6704E cards (National Instruments).
When the active source fiber is changed during the measurement, the gain of the detectors is set according to a table determined before the actual measurement. This extends the available dynamic range in a manner analogous to the VOAs used in the UCL system.13
Two laser diodes are used to provide a selection between 760 and . The laser diodes are temperature stabilized using a software-controlled thermoelectric cooler system. For this study, a laser diode operating at a wavelength of was installed to acquire data that could be directly compared with the 780-nm data from the UCL system (the phantom optical properties are practically constant within this wavelength range).
Specifications for the two systems are given in Table 1 , along with information on selected published optical tomography and phase measurement systems.
Characteristics of the two instruments with performance data from other published systems. TD=time-domain, FD=frequency-domain, CW=continuous-wave.
|System Feature||HUT Nissilä1||UCL Schmidt2||Schmitz3||McBride4||Chance5||Tromberg, Pham5, 6|
|Number of sources||16||32||25||16|
|Number of detectors||16||32||32||16|
|Wavelengths (nm)||760(785), 830||780,815||660,761 785,808 826,||744 to 859||674,811 849,956|
|Source power (mW)||8||40||0.02 to 3||25|
|Modulation frequency (MHz)||100||0 to 2000||100||1 to 500||0.3 to 1000|
|Amplitude drift (%/h)||1||0.2 to 1||0.6 to 20|
|Phase drift (deg/h)||0.05||0.18||0.01 to 0.2||0.007 to 7.5|
|Noise in (%)||0.5||1.7|
|Noise in (%),||0.1||0.03 to 1.9||0.6|
|Noise in (deg),||0.5||0.2||3|
|Noise in (deg),||0.08||0.05 to 0.34||0.3|
|Source cross talk||7|
|Detector cross talk|
|Detection limit at|
|Repeatability of (%)||0.4||4||0.5|
|Repeatability of (deg)||0.08||0.08||0.4|
|cross-talk (deg/dB)||negligible||negligible||0.04 to 0.9||0.05|
|Imaging time||to minutes||minutes|
|Hardware cost (euros)|
The UCL system has 32 source and 32 detector channels while the HUT system only has 16 of each. The VOAs of the UCL system reduce the optical power of the light reaching each detector so that the count rate does not exceed the maximum rate for photon counting detection (around photons per second per channel). This has the negative effect of reducing the available signal-to-noise ratio (SNR) at short source-detector separations. The VOAs also have a zero-transmission option that enables the detectors nearest to each activated source to be deactivated when the maximum attenuation of the VOA is insufficient to prevent saturation.
The multichannel MCP-PMTs used in the UCL system have more interchannel cross talk than the independent detectors used in the HUT system. Source cross talk in the UCL 32-way fiber switch is reduced to a negligible level by a piezoelectric shutter in series with each source fiber.
The powerful fiber laser of the UCL system allows measurement across thicker tissues than the laser diodes used in the HUT system. However, this is partially compensated by the higher QE of the detectors and the simpler and more efficient light collection used in the HUT system.
The timing of the sending and receiving electronics is synchronized by an optical reference pulse in the UCL system and by a phase-locked-loop in the HUT system. The coaxial fiber bundles allow frequent calibration of the detection channels in the UCL system by using the reflected pulse from the tissue as a reference. The HUT system has no equivalent of this but it has a smaller intrinsic drift.
The lower hardware cost of the HUT system makes it attractive from a commercial point of view, although increasing the channel count to 32 would double the cost of the system and the UCL system could be implemented at a reduced cost using newer components.
Modeling and Image Reconstruction
A standard approach for the modeling of light propagation in random media is the RTE.21 For a source modulated with an angular frequency it is written asat position and in direction . Above, is the speed of light, is the normalized scattering phase function representing the probability density function for scattering , and is the spatial and angular distribution of the source.
Numerical solutions to the RTE are computationally expensive, and to obtain a practical forward model for optical tomography, approximations need to be made. In most tissues , and in such cases, the DA is commonly used. In the frequency domain, the DA is written asis the scalar diffusion coefficient, in which is the mean of the cosine of the scattering angle, is the photon density, and presents an isotropic source term.
In this study, we compare measured data to simulated data from a FEM solution to the DA. For the boundary condition for the DA, we use the Robin condition, which can be derived from the assumption that the inward-directed photon current at each point on the boundary (except the source positions) is zero. The sources are represented as inward-directed diffuse photon currents over the illuminated area (covered by the source fibers) . The sources are then included directly in the boundary conditionsis the outward unit normal vector to the surface and the source current. In order to model the experimental situation accurately, we used a source current with a Gaussian spatial profile .
If the measured data is written as a complex amplitude, it is proportional to the exitance, which is defined on the boundary of the medium asof the position of the detecting bundle. Within the DA, using the boundary conditions described in Eq. 3, the exitance can be written as
The finite element model for the DA used in this study was implemented in 3D using second-order elements and basis functions. The mesh used to generate the simulated data closely follows the boundaries between the perturbations and the background, and it has nodes and second-order elements.
The image reconstructions were performed using a regularized Gauss-Newton algorithm and a mesh with nodes and linear elements to reduce reconstruction time. The homogeneous model baseline data were calculated using the DA model. The baseline and the measured difference data were added together to generate the perturbed absolute data that was used to generate the reconstructions. The objective function included a regularization term that is proportional to the sum of the -norms of the gradients of the and images. In addition to the use of explicit regularization, we selected the best image from the series of iterations using visual criteria. Images at the first iterations have a low spatial resolution while late iterations are usually quite noisy, and a compromise had to be made for each data set. The objective function for iteration and image vector (including both the scattering and absorption images) can be written asand are the source and detector indices; is the measured data; is the simulated data; is a parameter that controls the amount of regularization applied, , is the initial estimate, and is defined by is the number of neighbors of basis component .
The MC simulation program was implemented in C using the photon packet principle outlined in Ref. 22. The light source was modeled as a collimated beam with a radius of , and the detectors were assumed to collect all photons that hit the surface of the phantom at a distance of up to from the marked detector position. Internal reflections were modeled. photon packets were launched at each source position. The phase and amplitude at a modulation frequency of were calculated using the Fourier transform.
The two imaging systems were tested using a pair of solid cylindrical phantoms, which have the same external geometry and identical background optical properties. The phantoms were made from a mixture of particles and NIR-absorbing dye (ICI Projet 900NP) within epoxy resin.23 Each phantom has a diameter of and a height of .
Both phantoms were designed to have optical properties of and at a nominal wavelength of , and the expected optical properties of the background at the experimental wavelength are and . The index of refraction is 1.56, and the anisotropy factor is approximately 0.5.
One of the phantoms is homogeneous and the other contains two small cylindrical targets with optical properties of ( , ) and ( , ) relative to the background. Each target has a diameter of and a height of . Both targets are positioned within the same cross-sectional plane halfway between the top and bottom of the cylindrical phantom (Fig. 3 ). Due to uncertainties involved in the process of manufacture, the nominal optical properties of the phantom and the targets are expected to be accurate within about 10%. In order to insert the targets, cylindrical holes were drilled in the phantom with a flat-ended milling tool and the solid targets, which were manufactured earlier, were dropped into the holes. Liquid resin with the same optical properties as the background of the phantom was then poured into the holes. Because the cured targets are slightly more dense than the liquid resin, the targets are expected not to move much during the curing process. In some cases, x-ray imaging can be used to detect the position of scattering inhomogeneities; however, the contrast in our case was too small. Contrast agents that enhance x-ray detection but do not significantly affect NIR optical properties are under investigation. A small amount of metallic powder could be used for this effect.
In order to verify the correctness of the optical properties, the homogeneous phantom was measured with the HUT system and calibration was applied to the data.14, 24 The two parameters of the homogeneous model ( and ) and the coupling coefficient were varied to find an optimal fit with the experimental data. It is necessary to add the coupling coefficient to the measured data because the efficiency of the coupling of the light between the optical fibers and the phantom is unknown. The predicted optical properties were found to be correct within 3%.
The phantoms were also x-rayed to verify that they are free of air bubbles.
The optical fiber bundles were placed in two rings, each containing 16 bundles. The rings were apart in the axial direction. The optodes were placed at a short distance from the surface of the phantom to reduce the sensitivity of the measurements to small surface inhomogeneities and to make the FE model less sensitive to the properties of the mesh. The actual gap was in the HUT system and in the UCL system due to different connectors used. This slight discrepancy does not significantly affect the difference data. Otherwise, the fiber positions were the same in measurements carried out by both systems. In the presentation of the difference data in Sec. 3, we use the concept of measurement number frequently. The measurement number is simply an index to the data vector that is used to contain measured data. In the graphs of this paper, the measurement number can be mapped to the phantom geometry in the following way: The sources (1 to 16) measured by both systems have been numbered and marked with an “x” in Fig. 3, and the detectors have been marked with an “o.” The UCL system uses all positions as both source and detector but these additional data are only used in the reconstructions, not the direct data comparison. In the HUT data vector, the first 16 values correspond to source 1, detectors 1 to 16 in this order. The next 16 values (17 to 32) correspond to source 2 and detectors 1 to 16, and so on.
The phantoms were measured using seven different measurement times by both systems (1, 2, 5, 10, 15, 20, and per source position). The imaging time is the time that is required by the system to complete one full data acquisition sequence using all sources and detectors. The imaging time may be expressed asis the number of source positions, is the measurement time for each source position, and is the switching time. For the HUT system, consists of the switch delay and the time required to adjust the detector gains, and for the UCL system, it includes the time required to download the data from each of the 32 TCSPC units to a PC. is approximately for the HUT system and 3 to for the UCL system. The phantoms were repeatedly exchanged so that the measurements on the two phantoms were made consecutively for each measurement time.
The UCL system TPSFs were converted to frequency-domain data using the Fourier transform. Amplitude and phase were extracted at a frequency of . The differences in and phase between the phantoms with and without the targets were calculated and compared with the FEM data. The rms errors were calculated for each measurement time. Of the 32 detection channels of the UCL system, 7 exhibited significant temporal jitter and were excluded from the analysis. The model comparison was made only for those source-detector pairs that were measured with both instruments.
To test for stochastic errors independently of measurement accuracy, data from a homogeneous cylindrical phantom were measured repeatedly six times using 16 channels with .
The quality of the reconstructed images was evaluated using four different quantitative parameters. Although the images were reconstructed using a 3D algorithm, the image quality parameters were calculated based on the central [two-dimensional (2D)] cross section of the phantom. Both targets of interest were present in this plane, and due to the placement of the optical fibers, the reconstruction algorithm had limited control over the optical properties in regions far from the central cross section. A 3D analysis would be more appropriate if the optical fibers had been placed with even spacing throughout the surface of the phantom.
The cross section of the phantom was divided into two regions: the target region (different for and ) and the background. The target region was defined to be the area surrounding the center point of the target (the perturbation in either or ) with a radius of (the actual radius of the physical target is ). The background region was defined as the remaining part of the cross section, not part of the target region. The idea was to have a large enough search region to get a positive identification of the main objects to be reconstructed in the resulting image cross sections, while avoiding a misidentification of image artifacts as the targets. The effects of the targets on the reconstructed optical properties in the background region are minimal.
We define the relative contrast of the reconstructed target to be the peak value of the image in the target region divided by the mean of the image in the background region. The contrast-to-noise ratio (CNR) is defined to be (peak value of the image in the target region minus the mean of the background region) divided by the standard deviation of the background region.
The radius of the reconstructed target is calculated by thresholding the image in the target region at a value that is halfway between the peak value and the background mean and requiring eight-connectivity within the thresholded region. Two image pixels are eight-connected if they are either eight-neighbors of each other, or they have neighbors that are eight-connected with each other. Two pixels are eight-neighbors of each other if they are adjacent in one of eight directions (left, right, top, bottom, or diagonal). The effective radius of the image of the target is calculated fromis the area of the region that is higher than the threshold and is eight-connected to the peak value. The localization error in the image of the target is calculated as the Euclidian distance between the center of gravity of the thresholded target area and the true position of the target.
Results and Discussion
FEM Versus MC for the Reduced-Size Case
The FEM DA model was validated by calculating the simulated absolute and difference data for two phantoms that were otherwise similar to the phantoms used in the experiments, but they were reduced in size by 30% to enable the use of the MC method with sufficient statistical power. The simulated phase and amplitude data using the optical properties of the homogeneous phantom using the FEM DA and MC methods are shown in Fig. 4 , and simulated difference data using both models are shown in Fig. 5 for sources 1 to 8 [(a) and (b)]. The mismatch between the predictions between the two models is shown as a function of source-detector separation in Figs. 4c and 4d for absolute data and in Figs. 5c and 5d for difference data. The FEM data were calculated using the quadratic mesh in order to minimize numerical errors in the simulations. The agreement between the models is good, although there are subtle differences in the magnitudes of the effects of the perturbations in the data predicted by the two models. The differences between the two models increase as a function of increased source-detector separation for two reasons: (1) the contrast of the targets is greater at larger separations, which leads to larger discrepancies, and (2) the MC phase difference data is affected by stochastic noise, which is primarily visible in the data corresponding to the largest separations.
FEM versus Measured Data
The measured and FEM-predicted difference data at are shown in Figs. 6a and 6b . The experimental data corresponds to for both instruments. The mismatches between the measured and simulated difference data are shown in Figs. 6c and 6d. The rms differences between the simulated and measured data as a function of and are shown in Fig. 7 . The rms errors between the model and measurement averaged over all the imaging times are given for both instruments in Table 2 . Histograms of the errors between the measured and model-predicted difference data are given in Fig. 8 . The histograms include the data measured at all measurement times.
Summary of the results of the phantom measurements (averaged over all imaging times).
|HUT System||UCL System|
|rms error in phase|
|rms error in amplitude||3.4%||6.0%|
|Reproducibility of phase ( )|
|Reproducibility of amplitude ( )||0.4%||4%|
|Relative contrast in||1.19||1.15|
|Relative contrast in||1.25||1.15|
|Localization error in||1.5||1.6|
|Localization error in||2.3||2.0|
|Effective radius of the reconstructed target||8.7||9.3|
|Effective radius of the reconstructed target||7.0||6.9|
Differences between the measured and the model-predicted data may be partly explained by our lack of precise knowledge of the optical properties (the accuracy is estimated to be ) and the vertical positions of the targets. There is also a slight difference in the diameters of the two phantoms. The model-predicted phase difference data is more positive than the measured data, which is evident in the histograms of Fig. 8. Any mismatch between the model and physical reality prevents the error estimates from going to zero even in the case of a perfect measurement.
The phase error of the UCL system was quite independent of the measurement time at , while the phase error of the HUT system was between 0.42 and . At short measurement times , the accuracy of the phase data measured by the HUT system was limited by shot noise, which is approximately proportional to . The shot noise is due to the limited number of photons detected, as well as the dark current, which is significant since the PMTs used in the HUT system are not cooled. Phase drift, which causes an error approximately proportional to , reduced the accuracy significantly at longer measurement times . The best accuracy was found to be obtained at . Photon shot noise and thermal noise at short measurement times are less significant in the UCL system because of the higher power of the laser source and the thermoelectric cooling of the MCP-PMTs.
The amplitude difference data measured using the UCL system matched the model difference data to within 4.6 to 7.8%, while the HUT data matched the model to within 3.2 to 3.5%. The VOAs, the piezoelectric shutters, and the fiber laser all contribute to amplitude noise in the UCL system. The implementation of the VOAs using holes in a rotating disk makes the amplitude sensitive to the positional repeatability of the disks. This noise is independent of the measurement time. The piezoelectric shutters cause additional switching noise in the amplitude, which may reduce the quality of amplitude data especially at short measurement times. The laser intensity oscillates with a period of the order of a few minutes, which makes the amplitude accuracy sensitive to how the imaging time and the period of oscillation relate. In the HUT system, the most accurate amplitude data was obtained using the same measurement time, which optimized the quality of the phase data. Hysteresis may cause switching noise in the HUT system, but this is a subject for further research.
In summary, the measurement accuracy is a complex function of instrument design and implementation and our estimates of measurement accuracy can be understood by considering the different components of the system in detail and analyzing the corresponding noise sources theoretically.
A homogeneous cylindrical phantom was measured six times using both systems, and the standard deviations of amplitude and phase for each source-detector combination are shown in Fig. 9 . Note that due to the attenuation of light by the VOAs, the measured optical powers for the two instruments are of different orders of magnitude (Fig. 9).
The repeatability of phase data measured by the two systems was found to be similar at . The mean repeatability of the amplitude data from the HUT system was 0.4% and the corresponding value for the UCL system was 4% (Table 2).
In presenting the reconstructed images, we consider only the cross section through the plane , which crosses through the targets. In Fig. 10 , absorption and scatter images reconstructed from simulated data with no added noise are shown. In Figs. 11 and 12 , reconstructions calculated from the data recorded at three different imaging times for the HUT and UCL systems are shown. The true positions of the targets are indicated with black dashed circles in the images.
The image quality parameters were calculated for both systems at all measurement times and Fig. 13 shows the parameters as a function of imaging time. The relative contrast and CNR of the images reconstructed from the HUT data were somewhat higher than those from the UCL data at equivalent imaging times. We also reconstructed the optical properties based on noise-free simulated data with 16-source and 16-detector positions using a very low regularization parameter ( instead of the value used to generate all the images shown in this paper) and obtained a relative contrast of 1.3. The true contrasts of the targets (2:1) cannot be recovered due to the small volumes of the targets relative to the spatial resolution of the method.
In our measurement geometry, the differences in spatial resolution and contrast due to the different number of source and detector fibers in the two setups were not significant. The explanation for the lower contrast in the images reconstructed from the UCL data is that the iterative reconstruction algorithm produced severe artifacts on the boundary of the image (corresponding to outlier points in the data) if many iterations were calculated. In the case of the HUT data, a few additional iterations could be calculated without significant artifacts, which increased the contrast and spatial resolution of the images and also reduced the cross talk between the two optical properties. The number of iterations needed to get the best result from the UCL data ranged from 1 to 3, while 2 to 5 iterations were needed for the HUT data.
The localization accuracy of both and targets and the size of the target were found to be similar in the reconstructions calculated from data measured with both systems. The quality of the scattering images from the HUT data improved in localization accuracy and CNR as the imaging time was increased, while the quality of the absorption image was best at an intermediate measurement time of per source . No clear dependency between the image quality parameters and the imaging time was found in the images reconstructed from the UCL data. This is partly explained by the fact that a relatively small fraction of the noise in the UCL data is shot noise. The image quality parameters averaged over the imaging times are given in Table 2.
Several artifacts could be identified in the images. In particular, the area of the image concident with the position of the scattering target shows cross talk from the image. The interparameter cross talk is also present to some extent in the first iterations calculated from noise-free simulated data. The match between the FEM DA model and the physical reality is not perfect and there may be measurement errors such as phase-amplitude cross talk and noise that contribute to the interparameter cross talk. Artifacts on the boundary of the reconstructed image may be due to imperfections in the surfaces of the phantoms or changes in the phase or amplitude caused by the bending of the optical fibers.
Instrument Performance Relative to Other Systems
Table 1 includes performance figures of instruments described in the review paper by Chance, 4 a frequency-domain tomography system described by McBride, 5 and a fast continuous-wave optical tomography system by Schmitz 25 The system by the Tromberg group4, 26 is mentioned in a separate column because of its wide range of modulation frequencies.
In Table 1, all noise figures and the detection limits are given for a bandwidth of . Phase noise and drift are given for a modulation frequency of , with the exception of the systems in Ref. 4, for which the values are given at the respective modulation frequencies of each system. Phase and amplitude noise are given for two detected optical power levels. The first is and the second corresponds to a measurement with a source-detector separation of on a caucasian adult forehead. The latter corresponds to a detected power of when the HUT system is used. The source power in Table 1 is the optical power incident on the tissue. The repeatability figures for the McBride system were measured using a different phantom and a faster measurement time than what was used for the repeatability measurements in this paper .
The UCL and HUT systems compared in this paper have relatively good detection limits for instruments of their type, allowing the study of thick tissues. Another advantage the two systems share is accurate phase measurement due to the lack of phase-amplitude cross talk. Phase-amplitude cross talk in the HUT system is negligible in practical measurements when the anode current is below .14 Improvements in the imaging times in the two systems is a subject of ongoing research.
Advantages and Disadvantages of the Two Systems
Advantages and disadvantages of the technologies used in the two systems are summarized in Table 3 .
Summary of the advantages (+) and disadvantages (−) of the technologies used. OD=units of optical density.
|HUT System||UCL System|
|Light detection||Room-temperature PMTs Moderate detection limit No cross talk due to detectors||Cooled MCP-PMTs Excellent detection limit Cross talk within detector (maximum 0.3%) Risk of condensation|
|Dynamic range extension technique||Gain switching No light loss Hysteresis ( of change) Adjustable range: 4 OD||VOAs (implementation using holes) Light loss up to 3 OD Transmission variability Adjustable range: 2 to 3 OD|
|Light source||Laser diode Inexpensive ( euros/laser) Variety of wavelengths available||Fibrelaser Expensive ( euros) System bandwidth Optical power Amplitude may be unstable|
|Calibration||Amplitude and phase Offline only||Temporal data types Online or offline|
|Imaging time||Moderate ( to minutes) Long with very thick tissues||Long (minutes)|
The UCL system has a lower detection limit (noise equivalent power= versus in the HUT system, Table 1), the ability to record individual photons leading to a higher SNR at very low intensities, and a high laser power ( versus in the HUT system, Table 1), which enable it to be used to perform measurements across relatively large thicknesses of tissue. This is important in, for example, optical tomography of premature infants’ brains. In order to extend the dynamic range of the HUT system, the laser power should be increased and the detectors replaced by cooled PMTs. The UCL system has a greater variety of data types available, which may improve image quality in the future as reconstruction techniques advance. Thirty two channels also help to increase contrast in the reconstructions and provide a more even sampling of the tissue. Both the phase and amplitude data types measured using the HUT system can be calibrated using the procedure described in Refs. 14, 24. Temporal data types of the UCL system can be calibrated as described in Ref. 10 and temporal drifts in the detection channels can be cancelled on-line using the procedure described in Ref. 16.
One major disadvantage of the UCL system compared to the HUT system is the long switching time, which means that imaging times are significantly longer for equivalent measurement times. This becomes a severe limitation of the UCL system when attempting to acquire data quickly, such as to image dynamic processes in tissue, particularly in cortical activation imaging.27. However, an effort has recently begun to replace the system electronics with new TCSPC modules built by Becker and Hickl GmbH (Berlin, Germany).28 These modules use double buffering so that data can be recorded continuously while data from the previous source is being stored on a personal computer. When fully implemented, the UCL system will have a switching time comparable to the HUT system. The HUT system has the ability to record high-quality amplitude data even at short measurement times, and both amplitude and phase data recorded at short and intermediate distances have a better SNR. If an attempt is made to image a hemodynamic event that is faster than the imaging time of the system, it is possible that the spatial location of the event is inaccurately reconstructed as the different source fibers are active at different times and therefore they sample different optical property distributions in the tissue. An ideal system would have many channels and a fast imaging time so that both spatial and temporal undersampling are avoided. If the image acquisition time is long, modeling the physiological changes in the tissue is likely to improve the results and help separate the effects of systemic oscillations from the functional processes studied.29, 30
Recommendations for Instrument Designers
In general, if tissues significantly more attenuating than our test phantom are to be measured, cooling of the PMTs is needed to maintain good performance. In Fig. 9c, the sharp increase in phase noise at the lowest measured optical powers is indicative of the effects of the dark current. In addition to cooling of the detector, increasing the amplification of the intermediate frequency amplifiers after the PMTs in the HUT system would improve the SNR slightly; however, a programmable amplifier would be needed to implement this in a practical tomographic measurement system. To our knowledge, no group has implemented a frequency-domain optical tomography system using cooled PMTs, and some technical challenges are expected there. A factor of 7 reduction in the dark current is expected if the R7400U-02 PMTs are cooled to . A time-domain system with cooled detectors is a good choice for highly attenuating subjects such as premature infants’ brains, while frequency-domain systems should be sufficient for mammography. Future improvements in frequency-domain systems will take them closer to the low-light performance of time-domain systems, yet their overall system cost will still remain lower than that of equivalent time-domain systems. In practice, it is easier to achieve good low-light performance using the time-domain principle because the pulsed light source is not emitting photons simultaneously with the detection of the photons from the previous pulse and thus electrical isolation is not as critical. In contrast, the source and detection electronics of frequency-domain systems are simultaneously active, and careful shielding is necessary to achieve a good detection limit. Increasing the number of wavelengths is very expensive for time-domain systems, while it can be done with a minimal increase in the cost of a frequency-domain system if time multiplexing of the wavelengths is used.
Photon counting systems have the disadvantage that the detected light may need to be attenuated and the implementation of VOAs can both reduce the amount of light entering the detector and cause problems with the reproducibility of the measurements. The implementation of VOAs with neutral density filters instead of holes is preferable to reduce random variations in the transmission. In transmission-only time-domain mammography with breast compression, it is not necessary to use VOAs.8 Gain switching with modern low-hysteresis PMTs works well and can be extended by changing the gain in further amplifier stages, since the gain range of the PMT itself is limited.
The stability of the light source and detection system are of critical importance in optical tomography systems. In frequency-domain systems, temperature stabilization of the laser diodes helps, and careful design of the electronics to minimize temperature dependency reduces drift.
Either time- or frequency-domain systems can be used for dynamic imaging. The use of dynamic range extension techniques such as gain switching and VOAs limit the rate at which images can be obtained. A smaller dynamic range may be acceptable in certain applications, such as cortical activation imaging. In this case, a compromise may need to be made between spatial and temporal resolution. In dynamic imaging applications, SNR becomes more important than stability. The sensitivity of measurements to deeper parts of the tissue in a back-reflection geometry can be enhanced by selecting “late photons” in a time-resolved system,11, 31, 32 or by using higher-order moments of the TPSF.33 Optical cross-talk within the fiber-tissue interface is easier to detect and eliminate using time-domain measurements.
In typical tomographic applications, implementation details have a greater role in determining the quality of data and reconstructed images than the choice between the time- and frequency-domain techniques. Future reconstruction techniques may provide improvements in image quality based on time-domain data types, but the lower cost as well as better SNR at short separations are likely to remain favorable to frequency-domain systems in applications where the detection limit and depth sensitivity of these systems are sufficient.
The authors would like to acknowledge the financial support by the Jenny and Antti Wihuri Foundation, the Finnish graduate school Functional Studies in Medicine, the Academy of Finland, the Foundation of Technology (TES), the National Technology Agency of Finland, the Emil Aaltonen Foundation, the Instrumentarium Science Foundation, and the Wellcome Trust. I. N. would like to thank Professor Simon Arridge, Professor Erkki Somersalo, and Petri Hiltunen for their advice. Resources of the Finnish IT Centre for Science were used in this study.