2.1.

Instrumentation

The MF-DOT system is based on a network analyzer (Agilent Technologies, Palo Alto, California). The network analyzer not only measures the amplitude and the phase of the detected signals, but also provides the rf signal to modulate the amplitude of the laser diode sources. Figure 1 shows the schematic of the DOT system. The system employs four different wavelengths: $665\mathrm{nm}$ ( $50\mathrm{mW}$ , Intelite, Genoa, Nevada), $785\mathrm{nm}$ ( $75\mathrm{mW}$ , Thorlabs, Newton, New Jersey), $800\mathrm{nm}$ ( $150\mathrm{mW}$ , Intelite, Genoa, Nevada), and $830\mathrm{nm}$ ( $100\mathrm{mW}$ , Intelite, Genoa, Nevada). The laser diodes are driven by a dc current source (ILX Lightwave, Bozeman, Montana). A bias-T allows the modulation of the lasers with the modulation signal generated by the network analyzer at the selected frequencies. A four-way power splitter (Mini-Circuits, Brooklyn, New York) is used to divide the rf signal into four, and a separate amplifier (Mini-Circuits, Brooklyn, New York) is used to amplify each signal prior feeding to the bias-Ts. A temperature-controlled laser diode base (Thorlabs, Newton, New Jersey) is used for each laser diode. A fiber optic switch (DICon Fiberoptic, Richmond, California) is used to multiplex the laser diode outputs through the source fibers. It has $1.2\text{-}\mathrm{dB}$ loss, less than $\pm 0.03\text{-}\mathrm{dB}$ repeatability, maximum $-80\text{-}\mathrm{dB}$ cross talk, and approximately $300\text{-}\mathrm{ms}$ switching time. The source fibers connected to the output of the fiber optic switch lead to the fiber optic interface and transfer the lasers’ outputs to the sample. The four laser mounts that hold the lasers are placed in a metal box.

Fig. 1

The schematic diagram of the multifrequency DOT setup.

Photomultiplier tubes (PMT) are chosen as detectors able to measure low-level signals in fan-beam geometry. Eight PMTs (R7400U-20 Hamamatsu, Japan) are used, one for each detection site. The R7400U series is a subminiature photomultiplier tube $16\mathrm{mm}$ in diameter and $12\mathrm{mm}$ in length. This PMT has an eight-stage electron multiplier composed of metal channel dynodes and features excellent response time. It has a rise time of $0.78\mathrm{ns}$ and ${\mathrm{f}}_{3\mathrm{DB}}$ around $200\mathrm{MHz}$ . Although its performance is degraded with increasing frequency, it has a better frequency response compared to other PMTs due to its faster rise and fall times. Typical spectral response of R7400U-20 covers the 300- to 900-nm range with 50-mA/W cathode radiant sensitivity at 800 nm. A compact on-board-type high-voltage power supply (Hamamatsu, Japan) is used for each PMT. The gain of PMTs can be adjusted separately by applying voltages between $1\text{to}4\mathrm{V}$ to the control input of these high-voltage power supplies. These eight control voltages are supplied by a computer via a DAC card (National Instruments, Austin, Texas). The gain of each PMT can be adjusted individually and is determined at the beginning of an experiment. With the help of a photodiode (Hamamatsu, Japan) mounted on the rotation stage in the detection unit, the approximate intensity collected by each fiber can be monitored. For each source position, the photodiode is positioned right under each detector fiber sequentially, and the photodiode output is measured by a transimpedance amplifier (Terahertz Technologies, Oriskany, New York). The output of the transimpedance amplifier is monitored and recorded by the computer with a general purpose DAQ card (National Instruments, Austin, Texas). Based on these approximate values, the initial gain setting of each individual PMT is determined. During an imaging session, the gains of the PMTs are kept constant. The outputs of the PMTs are multiplexed by an $8\times 1$ rf switch (Hitite, Chelmsford, Massachusetts) and amplified by a $65\text{-}\mathrm{dB}$ amplifier (Sonoma Instruments, Sonoma, California). The rf switch uses a GaAs switch and has 0.8-dB insertion loss and more than $55\text{-}\mathrm{dB}$ isolation at $100\mathrm{MHz}$ .

The detectors nearest to the source receive orders of magnitude more incident light than those that are the farthest due to the detection geometry. Thus, either the light level or the gain of the PMT must be adjusted because of the limited dynamic range of the PMT. A filter wheel is inserted between the PMTs and fibers in the coupling system to adjust the light level. The PMTs and the detector fibers are arranged in a circular geometry and placed in a detection unit. A rotation stage with high repeatability is used to change the positions of the filters and keep them fixed relative to the source, as shown in Fig. 2 . Therefore, the extensive dynamic range of the intensity measured from the sample is reduced to the dynamic range of the PMTs. The positions of the PMTs and the detector fibers are kept fixed and do not involve any mechanical motion. The gain of each PMT can even be set to the same value with the adjustment of the filters. This made the calibration procedure easier and reduced any error potentially inherent in other types of calibration techniques. Keeping the PMT’s bias voltage constant during the experiment decreases the voltage hysteresis, while the neutral density filters limit the dynamic range of the light detected by the PMT, thus reducing the light hysteresis.

Fig. 2

The detection unit. The PMTs and the detector fibers were arranged in a circular geometry and fixed. A rotation stage was used to change the positions of the filters.

The PMTs, their high-voltage power supplies, the rf switch, the filter wheel, and the rotation stage are all placed in a shielded aluminum box. This detection unit has three levels, as shown in Fig. 2. The upper level holds the detector fibers and the high-voltage power supplies. The middle level is the aluminum filter wheel mounted on a high-precision rotation stage (Phytron, Waltham, Massachusetts). The filter wheel is $20\mathrm{in.}$ in diameter and can accommodate up to 32 filter holders. The rotation stage is driven by a step motor power drive (National Instruments, Austin, Texas) and a motion controller card (National Instruments, Austin, Texas). The bottom level accommodates the PMT holders, the rf switch, and the amplifier unit. All parts of this detection unit were precision machined. The power supply lines for the circuits and the rotation stage, the control voltages for the high-voltage power supplies, and the rf switch are filtered by feed-through filters located on one side of the aluminum box.

2.2.

Fiber Optic Interface

Optical fibers are used both to conduct light from the sources to the tissue and to transfer the collected light from tissue to detectors. $62.5\text{-}\mu \mathrm{m}$ core diameter gradient-index fibers are used as the source fibers, while $1.1\text{-}\mathrm{mm}$ core diameter step-index fibers are used as the detector fibers. An adaptive interface was constructed to hold the sample and position the tips of source and detector fibers around the full circumference of the sample (Fig. 3 ). It consists of eight-source and eight-detector fiber probes with radically adjustable holders. Fiber probes were machined in a cylindrical form with a center hole that contains the fiber. Custom-made glass ferrules were used to fix and polish the tips of the fibers with the help of a polishing machine (Ultratec, Santa Ana, California). After polishing, each ferrule was positioned at the end of each probe and glued. A pair of molybdenum disulfide (MDS) filled nylon bearings (McMaster-Carr, Atlanta, Georgia) hold and guide each probe accurately in radial position. A compression spring is used for each probe to provide the required pressure to hold the fiber on the sample surface. The design allows having the same pressure level at each probe position (Fig. 4 ). The sample is placed at the center of the fiber optic interface, and the fiber probes are radially adjusted until they are in contact with the sample.

Fig. 3

The fiber optic adaptive interface, fiber optic probes, and the solid phantom with a hole.

Fig. 4

The spring-loaded fiber optic probe.

2.3.

Data Acquisition

A PC controls and synchronizes all of the equipment. Lab Windows CVI (National Instruments, Austin, Texas) is utilized for programming. The network analyzer measures the amplitude and the phase information from each detector. It provides heterodyne detection, and therefore gives a higher accuracy in phase detection by down-converting a received signal to an intermediate frequency (IF) that is much lower than the signal itself. In this study, the IF bandwidth of the network analyzer is set to $15\mathrm{Hz}$ . The dynamic amplitude error contribution is $0.1\mathrm{dBm}$ between $+10$ and $-50\mathrm{dBm}$ , while the dynamic phase error contribution is around $0.5\mathrm{deg}$ in the same range. The phase and amplitude of each PMT signal is acquired for approximately $2.5\mathrm{s}$ , averaged, and transferred to the computer via a GPIB card (National Instruments, Austin, Texas). For each of the eight source locations, measurements are done for eight detector locations, giving 64 measurements for each wavelength. The rotation stage accomplishes $45\text{-}\mathrm{deg}$ movements required for the filter positioning in less than $300\mathrm{ms}$ . Hence, this motion can be accomplished in parallel with the fiber optic switching and no extra time is required. The measurement time for one wavelength, including monitoring time for each PMT signal and $300\text{-}\mathrm{ms}$ switching time for source multiplexing, is around $2.5\min$ .

2.4.

Phantom Design

Two $63\text{-}\mathrm{mm}$ -diam solid phantoms simulating tissue optical properties were constructed using the method described by Firbank, Oda, and Delpy.¹⁷ The optical properties of the phantoms were $\mu \mathrm{a}=0.0132{\mathrm{mm}}^{-1}$ and $\mu {\mathrm{s}}^{\prime }=0.86{\mathrm{mm}}^{-1}$ at $785\mathrm{nm}$ . One of them was used for calibration purposes (homogeneous case), and a $15\text{-}\mathrm{mm}$ hole was drilled on the other one to simulate different embedded objects (heterogeneous case). The hole was filled with varying concentrations of Intralipid (10%, F.K. Clayton, Clayton, North Carolina) and Indian ink (Winsor and Newton, England) in water to test the performance of the system. A separate solid shell was prepared to characterize the optical properties of the Intralipid and the Indian ink. Different concentrations of Intralipid and Indian ink in water were placed in the solid shell phantom and measurements were acquired. The absorption coefficient of the solid phantom was higher than the breast tissue to evaluate the sensitivity of the system. Due to high absorption, the measured signals were around subpicowatt range when the separation between the source-detector probes was maximum, which was nearly the minimum detectable signal range, especially at high frequencies.

2.5.

Calibration

Calibration is the pivotal part of the data acquisition due to the variation in characteristics of each PMT, optical fiber, and neutral density filters. An accurate calibration is achieved in three steps: filter calibration, detector calibration, and homogeneous phantom calibration.

2.6.

Data Analysis

The diffusion approximation to the transport equation is given as

The finite element method (FEM) is used for the solution of the forward problem defined by the diffusion equation. A FEM mesh of 1761 nodes and 3264 first-order triangular elements is constructed with denser element distribution underneath the boundary. The source is place at a location $1∕{\mu }_s^{\prime }$ below the surface. A coarser mesh with 289 nodes and 512 elements is used for the solution of the inverse problem. For the solution of the inverse problem, the relation between the measurements and the unknown absorption and scattering coefficients are linearized around an initial as:

First, the absorption and scattering images were reconstructed at discrete frequencies separately. Later, a multifrequency approach was used, in which a single absorption and scattering image pair was obtained using the whole dataset obtained at different modulation frequencies. In this approach, the sensitivity information coming from each frequency coupled together to form a single Jacobian. One of the advantages of the new form of the Jacobian is that it makes the reconstruction more immune to noise than the single frequency counterpart.

3.1.

System Performance

The most critical elements of the system that affect its performance are the PMTs. The network analyzer, the amplifier, and the rf switch all have a bandwidth of $500\mathrm{MHz}$ , while the PMTs have a bandwidth of $200\mathrm{MHz}$ . Another important factor is the rf shielding of the detection unit. The noise level can be reduced by proper shielding and grounding to yield greater signal-to-noise ratio (SNR). The noise level of the system was measured while the rf modulation signal for the lasers and the bias voltage for the PMTs $\left(900\mathrm{V}\right)$ were kept applied, but no light is incident on the sample. The noise level was around $-90\mathrm{dBm}$ up to $200\mathrm{MHz}$ , $-80\mathrm{dBm}$ up to $250\mathrm{MHz}$ , $-70\mathrm{dBm}$ up to $300\mathrm{MHz}$ , and $-55\mathrm{dBm}$ up to $350\mathrm{MHz}$ . While the noise level increases with the modulation frequency, the signal decreases due to limited bandwidth of the PMT. To determine the useful bandwidth of the system, a homogeneous phantom was placed into the fiber optic interface, and the source-detector combination that provides the lowest signal was determined. Later, the modulation frequency was swept from $100\text{to}350\mathrm{MHz}$ for this source-detector combination. The bandwidth of the system is selected as the frequency at which the SNR drops below 30, and is determined as $300\mathrm{MHz}$ (Fig. 5 ). The minimum detectable signal at this frequency is $-67\mathrm{dBm}$ , which corresponds to light power of $0.4\mathrm{pW}$ incident on the detector. The maximum modulation frequency in the phantom experiments was $280\mathrm{MHz}$ . Figure 6 shows the measured 64 amplitude values for the modulation frequencies of 110 and $280\mathrm{MHz}$ .

Fig. 5

The amplitude of the measured signal and noise with respect to modulation frequency.

Fig. 6

The measured 64 amplitude values for the modulation frequencies of 110 and $280\mathrm{MHz}$ .

Measurement repeatability is also an important factor that affects performance of the imaging system. The DOT instrument and the lasers were turned on $30\min$ before the experiments to enable the system to warm up. To separate the errors arising from the fiber positioning and the instrument itself, repeatability tests were conducted with and without the repositioning of the homogeneous phantom.²³ First, the phantom was placed in the fiber optic interface and measurements were repeated ten times. The average rms errors in the amplitude and phase measurements were found to be 0.5% and $0.2\mathrm{deg}$ . The repeated measurements were also acquired for the case in which the phantom is removed and repositioned between each measurement. To eliminate the errors due to the inhomogeneity of the phantom, it was placed in the same orientation each time it was repositioned. In this case, the average rms error in the amplitude and phase measurements were found to be 4% and $0.4\mathrm{deg}$ . This error was caused only by the fiber probe repositioning. Later, the phantom was also rotated $45\mathrm{deg}$ between each measurement and the same error values were obtained. Therefore, it was concluded that the error due to the spatial variations of the phantom optical properties is less than the fiber repositioning errors. At the end of each experiment, a homogeneous phantom is placed in the fiber optic interface to acquire calibration data. The errors due to fiber repositioning would result in errors in the calibration and therefore affect the reconstructed image quality. These errors would cause fluctuations in the recovered optical properties in the imaging volume close to the source and detector probes, in addition to the inaccuracy in the recovered optical properties of the inclusion.

3.2.

Phantom Imaging

For calibration purposes, a homogeneous phantom was used. A second phantom with a hole $15\mathrm{mm}$ in diameter was used to evaluate the performance of the system. The hole was positioned half the distance between the center and the edge. First, the concentration of the Intralipid and Indian ink was set to match the background optical properties of the solid phantom ( ${\mu }_a=0.0132{\mathrm{mm}}^{-1}$ and ${\mu }_{s^{\prime }}=0.86{\mathrm{mm}}^{-1}$ ). Using this concentration as a reference, solutions with contrasts of 1.5:1, 2:1, and 2.5:1 in absorption coefficients are prepared. The data are acquired at $785\mathrm{nm}$ with six modulation frequencies, 110, 150, 175, 200, 240, and $280\mathrm{MHz}$ . After the calibration, the absorption and scattering maps were reconstructed using the inverse solver. First, six rows of Fig. 7 show the reconstructed absorption and scattering maps for the case with the contrast ratio of 2.5:1 at six different modulation frequencies. The seventh row shows the mean images obtained by averaging over six frequencies, while the last row shows the results obtained using multifrequency reconstruction. Figure 8 shows the mean of the reconstructed absorption values for all contrast cases in a region of interest (ROI) of $1\mathrm{cm}$ in diameter and the expected values. It is seen that the reconstructed absorption values increase with the modulation frequency and approach the expected value for all cases. On the other hand, the variations in the background region, which is homogeneous, increase as well. To take into account both of these effects, another metric for measuring the image quality, the contrast-to-noise ratio (CNR), was defined:²⁴

Fig. 7

First six rows show the reconstructed absorption and scattering maps for the case with the contrast ratio of 2.5:1 at six different modulation frequencies. The seventh row shows the mean image obtained by averaging over six frequencies, and the last row shows the results obtained using multifrequency reconstruction. The colormap was adjusted so that black and white colors correspond to 2.5 and 0.5 times the optical properties of the background, respectively.

Fig. 8

The mean of the reconstructed absorption values in a ROI of $1\mathrm{cm}$ in diameter for the contrast ratios of 1.5:1 (a), 2:1 (b), and 2.5:1 (c). The dotted lines show the expected values.

where ${\mu }_{a\mathrm{ROI}}$ is the mean absorption coefficient for the ROI, ${\mu }_{a\text{back}}$ is the mean absorption coefficient for the background, ${\sigma }_{\mathrm{ROI}}$ is the standard deviation of the absorption distribution in the ROI, and ${\sigma }_{\text{back}}$ is the standard deviation of the absorption distribution in the background. The coefficients ${\omega }_{\mathrm{ROI}}$ and ${\omega }_{\text{back}}$ are introduced due to the area difference of the ROI and the background. The background region was defined as the whole area except for a $2.5\text{-}\mathrm{cm}$ -diam circular region, as shown in Fig. 9 .

Fig. 9

The object and the background regions selected for calculating the image quality related parameters.

Table 2

The values of the parameters used for evaluation of the image quality for the contrast ratio of 1.5:1, 2:1, and 2.5:1. The first six columns show the values of parameters for the reconstructed images at six different modulation frequencies. The seventh column shows the same parameters for the mean image obtained by averaging over six frequencies. The last column, on the other hand, shows the same parameters for multifrequency reconstruction.

A number of parameters were defined and used along with CNR for the evaluation of the image quality:

As seen in Table 2 , the error in the recovered absorption coefficient reduces as the modulation frequency increases for all cases. For example, for the contrast case of 2.5:1, the error for the recovered absorption coefficient in the ROI reduces from 34% at $110\mathrm{MHz}$ down to 5% at $280\mathrm{MHz}$ . On the other hand, the fluctuations in the recovered optical properties in the background increase with the modulation frequency, as seen in Fig. 7. As a result, the related parameters Babstd (standard deviation of absorption coefficient distribution in the background) and Bscamean (the mean scattering coefficient in the background) increase as well.

These results show that evaluating the image quality only with error in the recovered optical parameters of the object would be misleading as expected, therefore the contrast ratio parameter (CNR) would help in assessing the image quality. As seen from Table 2, the CNR increases up to a certain modulation frequency and then starts decreasing. For cases with the contrast ratios of 2:1 and 2.5:1, the highest CNR was obtained at $200\mathrm{MHz}$ , while for the case with the smallest contrast ratio (1.5:1), the CNR obtained at $200\mathrm{MHz}$ was comparable to the one obtained at $150\mathrm{MHz}$ . The values in the seventh column of Table 2 show the values for the parameters obtained after averaging all of the images obtained at six different modulation frequencies (the mean image). Although the CNR values for the mean image are higher than the values obtained at discrete modulation frequencies, the CNR for the multifrequency reconstruction result is the highest for all cases (last column in Table 2).

The full width at half maximums (FWHMs) obtained from the profiles in the $x$ and $y$ directions along the object were averaged to get a single parameter, FWHM. The FWHMs obtained from the profiles at discrete frequencies were similar, except smaller FMHMs along the $x$ direction were observed at $280\mathrm{MHz}$ . It is believed to be due to the increased variations in the background.