2.1.

Instrumentation

2.1.1.

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 $815\mathrm{nm}$ , 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 $62.5\mu \mathrm{m}∕125\mu \mathrm{m}$ ) is permanently integrated along the central axis of a corresponding detector bundle (approximately 3000 step index fibers with a core diameter of $40\mu \mathrm{m}$ and a cladding diameter of $50\mu \mathrm{m}$ ; 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 $3\times {10}^5$ 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.

Fig. 1

The UCL time-resolved system.

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.¹⁷

2.1.2.

Frequency-domain system developed at HUT

The frequency-domain optical tomography system used in this paper has 16 source fibers, 16 parallel detection channels, and 2 wavelengths¹⁴ (Fig. 2 ).

Fig. 2

The HUT frequency-domain system.

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 $1\times 2$ prism switch. The signals from the PMTs are mixed to an intermediate frequency of $5\mathrm{kHz}$ , 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.¹³

Two laser diodes are used to provide a selection between 760 and $830\mathrm{nm}$ . The laser diodes are temperature stabilized using a software-controlled thermoelectric cooler system. For this study, a laser diode operating at a wavelength of $785\mathrm{nm}$ 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).

2.2.

Modeling and Image Reconstruction

A standard approach for the modeling of light propagation in random media is the RTE.²¹ For a source modulated with an angular frequency $\omega$ it is written as

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 ${\mu }_a\gg {\mu }_s^{\prime }$ , and in such cases, the DA is commonly used. In the frequency domain, the DA is written as

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 $\partial \Omega$ (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) $\partial {\Omega }_s$ . The sources are then included directly in the boundary conditions

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 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 $279403$ nodes and $194703$ second-order elements.

The image reconstructions were performed using a regularized Gauss-Newton algorithm and a mesh with $35228$ nodes and $181278$ 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 ${\mathrm{L}}_2$ -norms of the gradients of the $\ln {\mu }_a$ and $\ln {\mu }_s^{\prime }$ 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 $k$ and image vector ${\mathbf{x}}_k$ (including both the scattering and absorption images) can be written as

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 $2.5\mathrm{mm}$ , and the detectors were assumed to collect all photons that hit the surface of the phantom at a distance of up to $2.5\mathrm{mm}$ from the marked detector position. Internal reflections were modeled. $1\times {10}^9$ photon packets were launched at each source position. The phase and amplitude at a modulation frequency of $100\mathrm{MHz}$ were calculated using the Fourier transform.

2.3.

Phantoms

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 ${\mathrm{TiO}}_2$ particles and NIR-absorbing dye (ICI Projet 900NP) within epoxy resin.²³ Each phantom has a diameter of $69.25\mathrm{mm}$ and a height of $110\mathrm{mm}$ .

Both phantoms were designed to have optical properties of ${\mu }_s^{\prime }=1\pm 0.1{\mathrm{mm}}^{-1}$ and ${\mu }_a=0.01\pm 0.001{\mathrm{mm}}^{-1}$ at a nominal wavelength of $800\mathrm{nm}$ , and the expected optical properties of the background at the experimental wavelength $780\mathrm{nm}$ are ${\mu }_a=0.0097{\mathrm{mm}}^{-1}$ and ${\mu }_s^{\prime }=1.04{\mathrm{mm}}^{-1}$ . The index of refraction $n$ is 1.56, and the anisotropy factor $g$ is approximately 0.5.

One of the phantoms is homogeneous and the other contains two small cylindrical targets with optical properties of ( ${\mu }_a$ , $2{\mu }_s^{\prime }$ ) and ( $2{\mu }_a$ , ${\mu }_s^{\prime }$ ) relative to the background. Each target has a diameter of $9.5\mathrm{mm}$ and a height of $9.5\mathrm{mm}$ . 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.

Fig. 3

The cylindrical phantom with perturbations: (a) 3D projection, (b) 2D top view with source and detector positions.

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 ( ${\mu }_a$ and ${\mu }_s^{\prime }$ ) and the $\ln A$ coupling coefficient were varied to find an optimal fit with the experimental data. It is necessary to add the $\ln A$ 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.

2.4.

Measurements

The optical fiber bundles were placed in two rings, each containing 16 bundles. The rings were $12\mathrm{mm}$ 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 $6\mathrm{mm}$ in the HUT system and $10\mathrm{mm}$ 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 $30\mathrm{s}$ per source position). The imaging time $t_I$ 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 as

The UCL system TPSFs were converted to frequency-domain data using the Fourier transform. Amplitude and phase were extracted at a frequency of $100\mathrm{MHz}$ . The differences in $\ln A$ 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 $t_M=5\mathrm{s}$ .

2.5.

Image Analysis

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 $(z=0)$ 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 ${\mu }_a$ and ${\mu }_s^{\prime }$ ) and the background. The target region was defined to be the area surrounding the center point of the target (the perturbation in either ${\mu }_a$ or ${\mu }_s^{\prime }$ ) with a radius of $13.5\mathrm{mm}$ (the actual radius of the physical target is $4.25\mathrm{mm}$ ). 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 $r_{eff}$ is calculated from

3.1.

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.

Fig. 4

MC and FEM DA absolute data for source fibers 1 to 8: (a) phase, (b) ln $A$ ; 3D FEM with a quadratic mesh (blue, –) and MC (red, - -). Model mismatch between 3D FEM and MC as a function of source-detector separation for (c) phase and (d) $\ln A.$

Fig. 5

MC and FEM DA difference data for (a) phase and (b) ln $A$ for source fibers 1 to 8. 3D FEM with a quadratic mesh (blue, –) and MC (red, - -). Model mismatch between 3D FEM and MC as a function of source-detector separation for (c) phase and (d) $\ln A$ .

3.2.

FEM versus Measured Data

The measured and FEM-predicted difference data at $780\mathrm{nm}$ are shown in Figs. 6a and 6b . The experimental data corresponds to $t_M=10\mathrm{s}$ 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 $t_M$ and $t_I$ 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.

Fig. 6

Measured and model-predicted difference data for the two phantoms using source fibers 1 to 8: UCL system (black), HUT system (red), 3D FEM (magenta). The mismatch between each measurement (HUT: red dots, UCL: black circles) and the 3D FEM model for (c) phase and (d) $\ln A$ , as a function of source-detector separation.

Fig. 7

The rms error between measured and model-predicted difference data: (a) Phase difference versus $t_M$ , (b) amplitude difference versus $t_M$ , (c) phase difference versus $t_I$ , (d) amplitude difference versus $t_I$ . UCL system (black, –), HUT system (red, - -).

Fig. 8

Histograms of the difference between the measured and model-predicted difference data: (a) Phase (UCL system), (b) $\ln A$ (UCL system), (c) phase (HUT system), and (d) $\ln A$ (HUT system).

Table 2

Summary of the results of the phantom measurements (averaged over all imaging times).

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 $\pm 10\%$ ) and the vertical positions $(\pm 2\mathrm{mm})$ of the targets. There is also a slight difference $\left(0.1\mathrm{mm}\right)$ 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 $0.45\mathrm{deg}$ , while the phase error of the HUT system was between 0.42 and $0.6\mathrm{deg}$ . At short measurement times $(t_M<5\mathrm{s})$ , the accuracy of the phase data measured by the HUT system was limited by shot noise, which is approximately proportional to $1∕\sqrt{t_M}$ . 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 $t_I$ , reduced the accuracy significantly at longer measurement times $(t_M\le 10\mathrm{s})$ . The best accuracy was found to be obtained at $t_M=5\mathrm{s}$ . 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.

3.3.

Repeatability

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).

Fig. 9

The standard deviation of the measured data over six repeated trials as a function of the optical power of the detected light with $t_M=5\mathrm{s}$ . Each dot corresponds to the repeatability of a different source-detector combination: (a) phase (UCL system), (b) amplitude (UCL system), (c) phase (HUT system), and (d) amplitude (HUT system).

The repeatability of phase data measured by the two systems was found to be similar at $0.08\mathrm{deg}$ . 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).

3.4.

Reconstructed Images

In presenting the reconstructed images, we consider only the cross section through the plane $z=0$ , 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 $\left(\text{-}\text{-}\right)$ in the images.

Fig. 10

Reconstructions from simulated noise-free data: (a) ${\mu }_a$ , (b) ${\mu }_s^{\prime }$ .

Fig. 11

Reconstructed images from the HUT data: (a) ${\mu }_a$ at $t_M=1\mathrm{s}$ , (b) ${\mu }_s^{\prime }$ at $t_M=1\mathrm{s}$ , (c) ${\mu }_a$ at $t_M=5\mathrm{s}$ , (d) ${\mu }_s^{\prime }$ at $t_M=5\mathrm{s}$ , (e) ${\mu }_a$ at $t_M=10\mathrm{s}$ , and (f) ${\mu }_s^{\prime }$ at $t_M=10\mathrm{s}$ .

Fig. 12

Reconstructed images from the UCL data: (a) ${\mu }_a$ at $t_M=1\mathrm{s}$ , (b) ${\mu }_s^{\prime }$ at $t_M=1\mathrm{s}$ , (c) ${\mu }_a$ at $t_M=5\mathrm{s}$ , (d) ${\mu }_s^{\prime }$ at $t_M=5\mathrm{s}$ , (e) ${\mu }_a$ at $t_M=10\mathrm{s}$ , and (f) ${\mu }_s^{\prime }$ at $t_M=10\mathrm{s}$ .

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 ( $\tau ={10}^{-7}$ instead of the value $\tau ={10}^{-5}$ 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.

Fig. 13

Image quality parameters as a function of imaging time. Relative contrast in (a) ${\mu }_a$ , (b) ${\mu }_s^{\prime }$ ; CNR in (c) ${\mu }_a$ , (d) ${\mu }_s^{\prime }$ ; localization error in (e) ${\mu }_a$ ; (f) ${\mu }_s^{\prime }$ ; effective radius of the reconstructed (g) ${\mu }_a$ , and (h) ${\mu }_s^{\prime }$ targets; UCL system (black, –), HUT system (red, - -).

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 ${\mu }_s^{\prime }$ 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 ${\mu }_a$ and ${\mu }_s^{\prime }$ targets and the size of the ${\mu }_s^{\prime }$ 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 $5\mathrm{s}$ per source $(t_I\approx 90\mathrm{s})$ . 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 ${\mu }_a$ image concident with the position of the scattering target shows cross talk from the ${\mu }_s^{\prime }$ 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.

3.5.

Instrument Performance Relative to Other Systems

Table 1 includes performance figures of instruments described in the review paper by Chance, ⁴ a frequency-domain tomography system described by McBride, ⁵ and a fast continuous-wave optical tomography system by Schmitz ²⁵ The system by the Tromberg group^{4, 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 $1\mathrm{Hz}$ . Phase noise and drift are given for a modulation frequency of $100\mathrm{MHz}$ , with the exception of the systems in Ref. ⁴, 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 $1\mathrm{pW}$ and the second corresponds to a measurement with a source-detector separation of $3\mathrm{cm}$ on a caucasian adult forehead. The latter corresponds to a detected power of $\sim 50\mathrm{pW}$ 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 $(t_M<2\mathrm{s})$ than what was used for the repeatability measurements in this paper $(t_M=5\mathrm{s})$ .

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 $70\mathrm{nA}$ .¹⁴ Improvements in the imaging times in the two systems is a subject of ongoing research.

3.6.

Advantages and Disadvantages of the Two Systems

Advantages and disadvantages of the technologies used in the two systems are summarized in Table 3 .

Table 3

Summary of the advantages (+) and disadvantages (−) of the technologies used. OD=units of optical density.

The UCL system has a lower detection limit (noise equivalent power= $20\mathrm{aW}$ versus $1\mathrm{fW}$ 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 ( $40\mathrm{mW}$ versus $8\mathrm{mW}$ 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.²⁷. However, an effort has recently begun to replace the system electronics with new TCSPC modules built by Becker and Hickl GmbH (Berlin, Germany).²⁸ 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}

3.7.

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 $0°\mathrm{C}$ . 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.⁸ 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.³³ 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.