1 November 2006 Comparison between a time-domain and a frequency-domain system for optical tomography
Author Affiliations +
The quality of phase and amplitude data from two medical optical tomography systems were compared. The two systems are a 32-channel time-domain system developed at University College London (UCL) and a 16-channel frequency-domain system developed at Helsinki University of Technology (HUT). Difference data measured from an inhomogeneous and a homogeneous phantom were compared with a finite-element method (diffusion equation) and images of scattering and absorption were reconstructed based on it. The measurements were performed at measurement times between 1 and 30 s per source. The mean rms errors in the data measured by the HUT system were 3.4% for amplitude and 0.51 deg for phase, while the corresponding values for the UCL data were 6.0% and 0.46 deg, respectively. The reproducibility of the data measured with the two systems was tested with a measurement time of 5 s per source. It was 0.4% in amplitude for the HUT system and 4% for the UCL system, and 0.08 deg in phase for both systems. The image quality of the reconstructions from the data measured with the two systems were compared with several quantitative criteria. In general a higher contrast was observed in the images calculated from the HUT data.



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 ( μa and μs ) 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 ( HbO2 and Hb). Changes in HbO2 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 500MHz ). 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 815nm , 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μm125μm ) is permanently integrated along the central axis of a corresponding detector bundle (approximately 3000 step index fibers with a core diameter of 40μm and a cladding diameter of 50μ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×105 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.17


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 wavelengths14 (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×2 prism switch. The signals from the PMTs are mixed to an intermediate frequency of 5kHz , 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 830nm . The laser diodes are temperature stabilized using a software-controlled thermoelectric cooler system. For this study, a laser diode operating at a wavelength of 785nm 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).


Technical comparison

Specifications for the two systems are given in Table 1 , along with information on selected published optical tomography and phase measurement systems.

Table 1

Characteristics of the two instruments with performance data from other published systems. TD=time-domain, FD=frequency-domain, CW=continuous-wave.

System FeatureHUT Nissilä1UCL Schmidt2Schmitz3McBride4Chance5Tromberg, Pham5, 6
Instrument typeFDTDCWFDFDFD
Number of sources16322516
Number of detectors16323216
Wavelengths (nm)760(785), 830780,815660,761 785,808 826,744 to 859674,811 849,956
Source power (mW)840 30 310 0.02 to 325
Modulation frequency (MHz)1000 to 2000 0.01 1001 to 5000.3 to 1000
Amplitude drift (%/h)1 <2 0.2 to 10.6 to 20
Phase drift (deg/h) to 0.20.007 to 7.5
Noise in A (%) 1pW 0.51.7
Noise in A (%), 3cm 0.10.03 to 1.90.6
Noise in φ (deg), 1pW 0.50.23
Noise in φ (deg), 3cm 0.080.05 to 0.340.3
Source cross talk 108 7 1014 104
Detector cross talk <103 <3×103 106
Detection limit at 1Hz 1fW 20aW 1pW 0.5pW
Repeatability of A (%)0.440.5
Repeatability of φ (deg)
φA cross-talk (deg/dB)negligiblenegligible0.04 to 0.90.05
Imaging time 15s to minutesminutes <1s <30s
Detector QE (800nm) 7%84%9
Hardware cost (euros) 80,000k 800,000k


Reference 12.


Reference 11.


Reference 23.


Reference 5.


Reference 4.


Reference 26.


Reference 18.


Reference 19.


Reference 20.

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 3×105 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 as


The RTE describes the change in radiance L(r,ŝ;ω) at position r and in direction ŝ . Above, c is the speed of light, f(ŝ,ŝ) is the normalized scattering phase function representing the probability density function for scattering ŝŝ , and q(r,ŝ;ω) 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 μaμs , and in such cases, the DA is commonly used. In the frequency domain, the DA is written as


where κ=13(μa+(1g)μs)1 is the scalar diffusion coefficient, in which g is the mean of the cosine of the scattering angle, Φ is the photon density, and Q 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) Ωs . The sources are then included directly in the boundary conditions


where n̂ is the outward unit normal vector to the surface and Γs the source current. In order to model the experimental situation accurately, we used a source current with a Gaussian spatial profile (σ=1.4mm) .

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 exitance is averaged on the boundary within 2.5mm of the position of the detecting bundle. Within the DA, using the boundary conditions described in Eq. 3, the exitance can be written as


For comparison with the measured data, the exitance was calculated as the average of the nodal values within the boundary area between the detecting fiber bundles and the phantom.

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 L2 -norms of the gradients of the lnμa and lnμs 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 xk (including both the scattering and absorption images) can be written as


where s and d are the source and detector indices; y is the measured data; F is the simulated data; τ is a parameter that controls the amount of regularization applied,


where Δxk=xkx0 , x0 is the initial estimate, and L is defined by


where nn is the number of neighbors of basis component i .

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.5mm , and the detectors were assumed to collect all photons that hit the surface of the phantom at a distance of up to 2.5mm from the marked detector position. Internal reflections were modeled. 1×109 photon packets were launched at each source position. The phase and amplitude at a modulation frequency of 100MHz 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 TiO2 particles and NIR-absorbing dye (ICI Projet 900NP) within epoxy resin.23 Each phantom has a diameter of 69.25mm and a height of 110mm .

Both phantoms were designed to have optical properties of μs=1±0.1mm1 and μa=0.01±0.001mm1 at a nominal wavelength of 800nm , and the expected optical properties of the background at the experimental wavelength 780nm are μa=0.0097mm1 and μs=1.04mm1 . 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 ( μa , 2μs ) and ( 2μa , μs ) relative to the background. Each target has a diameter of 9.5mm and a height of 9.5mm . 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 ( μa and μs ) and the lnA coupling coefficient were varied to find an optimal fit with the experimental data. It is necessary to add the lnA 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 12mm 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 6mm in the HUT system and 10mm 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 30s per source position). The imaging time tI 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


where NS is the number of source positions, tM is the measurement time for each source position, and tS is the switching time. For the HUT system, tS 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. tS is approximately 0.6s for the HUT system and 3 to 5s 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 100MHz . The differences in lnA 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 tM=5s .


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 μa and μs ) and the background. The target region was defined to be the area surrounding the center point of the target (the perturbation in either μa or μs ) with a radius of 13.5mm (the actual radius of the physical target is 4.25mm ). 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 reff is calculated from


where A is 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.

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


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



FEM versus Measured Data

The measured and FEM-predicted difference data at 780nm are shown in Figs. 6a and 6b . The experimental data corresponds to tM=10s 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 tM and tI 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) lnA , as a function of source-detector separation.


Fig. 7

The rms error between measured and model-predicted difference data: (a) Phase difference versus tM , (b) amplitude difference versus tM , (c) phase difference versus tI , (d) amplitude difference versus tI . 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) lnA (UCL system), (c) phase (HUT system), and (d) lnA (HUT system).


Table 2

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

HUT SystemUCL System
rms error in phase 0.51deg 0.46deg
rms error in amplitude3.4%6.0%
Reproducibility of phase ( tM=5s ) 0.08deg 0.08deg
Reproducibility of amplitude ( tM=5s )0.4%4%
Relative contrast in μa 1.191.15
Relative contrast in μs 1.251.15
CNR in μa 6.25.6
CNR in μs 6.25.2
Localization error in μa 1.51.6
Localization error in μs 2.32.0
Effective radius of the reconstructed μa target8.79.3
Effective radius of the reconstructed μs target7.06.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 ±10% ) and the vertical positions (±2mm) of the targets. There is also a slight difference (0.1mm) 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.45deg , while the phase error of the HUT system was between 0.42 and 0.6deg . At short measurement times (tM<5s) , the accuracy of the phase data measured by the HUT system was limited by shot noise, which is approximately proportional to 1tM . 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 tI , reduced the accuracy significantly at longer measurement times (tM10s) . The best accuracy was found to be obtained at tM=5s . 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).

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 tM=5s . 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.08deg . 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).


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 (--) in the images.

Fig. 10

Reconstructions from simulated noise-free data: (a) μa , (b) μs .


Fig. 11

Reconstructed images from the HUT data: (a) μa at tM=1s , (b) μs at tM=1s , (c) μa at tM=5s , (d) μs at tM=5s , (e) μa at tM=10s , and (f) μs at tM=10s .


Fig. 12

Reconstructed images from the UCL data: (a) μa at tM=1s , (b) μs at tM=1s , (c) μa at tM=5s , (d) μs at tM=5s , (e) μa at tM=10s , and (f) μs at tM=10s .


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 ( τ=107 instead of the value τ=105 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) μa , (b) μs ; CNR in (c) μa , (d) μs ; localization error in (e) μa ; (f) μs ; effective radius of the reconstructed (g) μa , and (h) μs 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 μs 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 μa and μs targets and the size of the μs 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 5s per source (tI90s) . 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 μa image concident with the position of the scattering target shows cross talk from the μs 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 1Hz . Phase noise and drift are given for a modulation frequency of 100MHz , 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 1pW and the second corresponds to a measurement with a source-detector separation of 3cm on a caucasian adult forehead. The latter corresponds to a detected power of 50pW 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 (tM<2s) than what was used for the repeatability measurements in this paper (tM=5s) .

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 70nA .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 .

Table 3

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

HUT SystemUCL System
Light detectionRoom-temperature PMTs Moderate detection limit (1fW) + No cross talk due to detectorsCooled MCP-PMTs + Excellent detection limit (20aW) Cross talk within detector (maximum 0.3%) Risk of condensation
Dynamic range extension techniqueGain switching + No light loss Hysteresis ( 0.2% of change) + Adjustable range: 4 ODVOAs (implementation using holes) Light loss up to 3 OD Transmission variability ± Adjustable range: 2 to 3 OD
Light sourceLaser diode + Inexpensive ( 1000 euros/laser) + Variety of wavelengths availableFibrelaser Expensive ( 200000 euros) + System bandwidth (2GHz) + Optical power (40mW) Amplitude may be unstable
Calibration + Amplitude and phase Offline only ± Temporal data types + Online or offline
Imaging time ± Moderate ( 15s to minutes) Long with very thick tissues Long (minutes)

The UCL system has a lower detection limit (noise equivalent power= 20aW versus 1fW 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 ( 40mW versus 8mW 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 0°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.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.



Hebden J. C., Arridge S. R., Delpy D. T., “Optical imaging in medicine I: Experimental techniques,” Phys. Med. Biol., 42 825 –840 (1997). 10.1088/0031-9155/42/5/007 0031-9155 Google Scholar


Arridge S. R., “Optical tomography in medical imaging,” Inverse Probl., 15 R41 –R93 (1999). 10.1088/0266-5611/15/2/022 0266-5611 Google Scholar


Arridge S. R., Lionheart W. R. B., “Non-uniqueness in diffusion-based optical tomography,” Opt. Lett., 23 882 –884 (1998). 0146-9592 Google Scholar


Chance B., Cope M., Gratton E., Ramanujam N., Tromberg B., “Phase measurement of light absorption and scatter in human tissue,” Rev. Sci. Instrum., 69 3457 –3481 (1998). 10.1063/1.1149123 0034-6748 Google Scholar


McBride T. O., Pogue B. W., Jiang S., Österberg U. L., Paulsen K. D., “A parallel-detection frequency-domain near-infrared tomography system for hemoglobin imaging of the breast in vivo,” Rev. Sci. Instrum., 72 1817 –1824 (2001). 10.1063/1.1344180 0034-6748 Google Scholar


Gratton E., Mantulin W. M., van de Ven M. J., Fishkin J. B., Maris M. B., Chance B., “A novel approach to laser tomography,” Bioimaging, 1 40 –46 (1993). 0966-9051 Google Scholar


O’Leary M. A., Boas D. A., Chance B., Yodh A. G., “Experimental images of heterogeneous turbid media by frequency domain diffusing photon tomography,” Opt. Lett., 20 426 –428 (1995). 0146-9592 Google Scholar


Ntziachristos V., Ma X., Chance B., “Time-correlated single photon counting imager for simultaneous magnetic resonance and near-infrared mammography,” Rev. Sci. Instrum., 69 4221 –4233 (1998). 10.1063/1.1149235 0034-6748 Google Scholar


Eda H., Oda I., Ito Y., Wada Y., Oikawa Y., Tsunazawa Y., Takada M., Tsuchiya Y., Yamashita Y., Oda M., Sassaroli A., Yamada Y., Tamura M., “Multichannel time-resolved optical tomographic imaging system,” Rev. Sci. Instrum., 70 3595 –3602 (1999). 10.1063/1.1149965 0034-6748 Google Scholar


Hillman E. M. C., Hebden J. C., Schmidt F. E. W., Arridge S. R., Schweiger M., Deghani H., Delpy D. T., “Calibration techniques and datatype extraction for time-resolved optical tomography,” Rev. Sci. Instrum., 71 3415 –3427 (2000). 10.1063/1.1287748 0034-6748 Google Scholar


Selb J., Stott J. J., Franceschini M. A., Sorenson A. G., Boas D. A., “Improved sensitivity to cerebral dynamics during brain activation with a time-gated optical system: Analytical model and experimental validation,” J. Biomed. Opt., 10 011013-1 –011013-12 (2005). 10.1117/1.1852553 1083-3668 Google Scholar


Nissilä I., Noponen T., Heino J., Kajava T., Katila T., Lin J. C., “Diffuse optical imaging,” Advances in Electromagnetic Fields in Living Systems, 4 77 –130 Springer-Verlag, New York (2005). Google Scholar


Schmidt F. E. W., Fry M. E., Hillman E. M. C., Hebden J. C., Delpy D. T., “A 32-channel time-resolved instrument for medical optical tomography,” Rev. Sci. Instrum., 71 256 –265 (2000). 10.1063/1.1150191 0034-6748 Google Scholar


Nissilä I., Noponen T., Kotilahti K., Katila T., Lipiäinen L., Tarvainen T., Schweiger M., Arridge S., “Instrumentation and calibration methods for the multichannel measurement of phase and amplitude in optical tomography,” Rev. Sci. Instrum., 76 044302-1 –044302-10 (2005). 10.1063/1.1884193 0034-6748 Google Scholar


Schweiger M., Arridge S. R., Nissilä I., “Gauss-Newton method for image reconstruction in diffuse optical tomography,” Phys. Med. Biol., 50 2365 –2386 (2005). 10.1088/0031-9155/50/10/013 0031-9155 Google Scholar


Hebden J. C., Gonzalez F. M., Gibson A., Hillman E. M. C., Yusof R., Everdell N., Delpy D. T., Zaccanti G., Martelli F., “Assessment of an in situ temporal calibration method for time-resolved optical tomography,” J. Biomed. Opt., 8 87 –92 (2003). 10.1117/1.1528206 1083-3668 Google Scholar


Hebden J. C., Gibson A., Austin T., Yusof R., Everdell N., Delpy D. T., Arridge S. R., Meek J. H., Wyatt J. S., “Imaging changes in blood volume and oxygenation in the newborn infant brain using three-dimensional optical tomography,” Phys. Med. Biol., 49 1117 –1130 (2004). 10.1088/0031-9155/49/7/003 0031-9155 Google Scholar


Dicon Fiberoptics VX500 data sheet, (2001) http://www.diconfiberoptics.com Google Scholar


Metal package PMT photosensor modules, H5773/H5783/H6779/H6780 series data sheet, Hamamatsu Photonics. Google Scholar


Schmidt F. E. W., “Development of a time-resolved optical tomography system for neonatal brain imaging,” University College London, (1999). Google Scholar


Chandrasekhar S., Radiative Transfer, Dover, New York (1960). Google Scholar


Prahl S. A., Keijzer M., Jacques S. L., Welch A. J., “A Monte Carlo model of light propagation in tissue,” 102 –111 (1989) Google Scholar


Firbank M., Oda M., Delpy D. T., “An improved design for a stable and reproducible phantom material for use in near-infrared spectroscopy and imaging,” Phys. Med. Biol., 40 955 –961 (1995). 10.1088/0031-9155/40/5/016 0031-9155 Google Scholar


Tarvainen T., Kolehmainen V., Vauhkonen M., Vanne A., Gibson A. P., Schweiger M., Arridge S. R., Kaipio J. P., “Computational calibration method for optical tomography,” Appl. Opt., 44 1879 –1888 (2005). 10.1364/AO.44.001879 0003-6935 Google Scholar


Schmitz C. H., Löcker M., Lasker J. M., Hielscher A. H., Barbour R. L., “Instrumentation for fast functional optical tomography,” Rev. Sci. Instrum., 73 429 –439 (2002). 10.1063/1.1427768 0034-6748 Google Scholar


Pham T. H., Coquoz O., Fishkin J. B., Anderson E., Tromberg B. J., “Broad bandwidth frequency domain instrument for quantitative tissue optical spectroscopy,” Rev. Sci. Instrum., 71 2500 –2513 (2000). 10.1063/1.1150665 0034-6748 Google Scholar


Kotilahti K., Nissilä I., Huotilainen M., Mäkelä R., Gavrielides N., Noponen T., Björkman P., Fellman V., Katila T., “Bilateral hemodynamic responses to auditory stimulation in newborn infants,” 1373 –1377 (2005). Google Scholar


Becker W., Bergmann A., Hink M. A., König K., Biskup C., “Fluorescence lifetime imaging by time-correlated single-photon counting,” Microsc. Res. Tech., 63 58 –66 (2004). 10.1002/jemt.10421 1059-910X Google Scholar


Prince S., Kolehmainen V., Kaipio J. P., Franceschini M. A., Boas D., Arridge S. R., “Time-series estimation of biological factors in optical diffusion tomography,” Phys. Med. Biol., 48 1491 –1504 (2003). 10.1088/0031-9155/48/11/301 0031-9155 Google Scholar


Kolehmainen V., Prince S., Arridge S. R., Kaipio J. P., “State estimation approach to non-stationary optical tomography problem,” J. Opt. Soc. Am. A, 20 876 –889 (2003). 0740-3232 Google Scholar


Steinbrink J., Wabnitz H., Obrig H., Villringer A., Rinneberg H., “Determining changes in NIR absorption using a layered model of the human head,” Phys. Med. Biol., 46 879 –896 (2001). 10.1088/0031-9155/46/3/320 0031-9155 Google Scholar


Montcel B., Chabrier R., Poulet P., “Detection of cortical activation with time-resolved diffuse optical methods,” Appl. Opt., 44 1942 –1947 (2005). 10.1364/AO.44.001942 0003-6935 Google Scholar


Liebert A., Wabnitz H., Steinbrink J., Obrig H., Möller M., Macdonald R., Villringer A., Rinneberg H., “Time-resolved multidistance near-infrared spectroscopy of the adult head: Intracerebral and extracerebral absorption changes from moments of distribution of times of flight of photons,” Appl. Opt., 43 3037 –3047 (2004). 10.1364/AO.43.003037 0003-6935 Google Scholar
© (2006) Society of Photo-Optical Instrumentation Engineers (SPIE)
Ilkka T. Nissilä, Ilkka T. Nissilä, Jeremy C. Hebden, Jeremy C. Hebden, David Jennions, David Jennions, Jenni Heino, Jenni Heino, Martin Schweiger, Martin Schweiger, Kalle M. Kotilahti, Kalle M. Kotilahti, Tommi E.J. Noponen, Tommi E.J. Noponen, Adam P. Gibson, Adam P. Gibson, S. Jaervenpaeae, S. Jaervenpaeae, Lauri Lipiäinen, Lauri Lipiäinen, Toivo E. Katila, Toivo E. Katila, "Comparison between a time-domain and a frequency-domain system for optical tomography," Journal of Biomedical Optics 11(6), 064015 (1 November 2006). https://doi.org/10.1117/1.2400700 . Submission:

Back to Top