1. Introduction

Diffuse fluorescence tomography has been developed for over a decade with basic measurements and algorithms showing proof of concept [1–6], followed by groundbreaking studies in murine models showing the clinical value of these methods [7–10]. In recent years, several systems have been produced [11–16], and now commercial systems which achieve this are commonly available. Still optimization and performance analysis of these systems is not well developed, partly because of the proprietary nature of commercialized development, and partly because of the unregulated nature of pre-clinical imaging. Two of the most problematic areas in system development for diffuse fluorescence tomography are in the dynamic range of the instrument and how it performs, in terms of the effective noise level. In this study, we examine these features for a high performance academic system, developed to couple with MicroCT imaging.

The system examined here has been reported previously [16,17] and was designed to utilize a rotating gantry, allowing for the use of a single source with a fan beam of detectors that revolves around the surface of the animal in a non-contact manner. The signal generation and detection is based upon pulsed diode laser and time-correlated single photon counting approach, which has been used for well over two decades in a range of systems [18–20]. Currently this technology has been widely commercialized by Becker & Hickl, and adopted in a few fluorescence instruments. One of the single most problematic issues with this type of technology is that it is highly sensitive to measurement noise, and so appropriate attenuation and temporal constraints must be applied in order to get the signal levels in an appropriate range, where the generated amplitude is linear with incoming signal level. This issue is particularly problematic in diffuse tomography where signal changes by several orders of magnitude can easily be induced by variations in tissue thickness through which the measurements are taken. In this paper, we focus on an automated exposure control sub-system, which allows the tomographic signal measurement to proceed in an automated manner, and analysis of such effects of the performance of the imaging system and algorithms.

The second issue examined in this paper is the one of effective system noise level. Model-based reconstruction approaches are inherently difficult to assess imaging performance, because of the data calibration which is required to make the raw data match the model as well as is possible. In previous studies [21–23], it was demonstrated how diffuse tomography is almost always dominated by bias errors, rather than random errors, and the model-data mismatch which occurs in the calibration process determines how low the bias error can be. Unfortunately estimating bias error is complicated by the calibration process and the model and system performance, and so bias error is rarely used as a metric in assessing diffuse tomography instrumentation. In this paper, the ‘effective’ random noise level of the system is examined, which generates a single performance metric which integrates all the errors, but comparing the reconstruction accuracy to simulations with only random error. The approach is approximate, but yields an ‘effective noise’ level which is based upon the lowest concentration which can be accurately reconstructed for this system.

2. Materials and methods

2.1 Instrumentation

The system used in this study is described in detail in a previous paper [17]. Briefly, the system consists of a custom designed computer rack integrated with time-correlated single photon counting (TCSPC) instrumentation cards (Becker and Hickl, Berlin, Germany). Custom LabVIEW (National Instruments, Austin, TX) software programs have been developed, allowing automated control of the motion and data acquisition. The laser diode driver module (PicoQuant, Berlin, Germany) controls a diode at 635nm, pulsing at 80 MHz repetition frequency. The laser is coupled to a 50 μm fiber optic beamsplitter delivering 5% of the signal to a reference channel and the other 95% toward a collimated source fiber directed at the specimen. A motorized attenuator is inserted inline to the source delivery fiber, allowing dynamic adjustment to ensure that the detected transmission and fluorescence signals fall within the linear range of the photomultiplier tubes (PMTs). This is a critical part of the initial set up, and described in detail previously.

Diffuse fluorescence and excitation signals measurements are then collected, in a transmission configuration, by 5 optical focusers with an angular separation of 22.5°, having a 97mm working distance from the tissue. This fan beam geometry allows parallel detection at the 5 locations for each source location. The optical focusers are coupled to 400μm fibers (NA: 0.37), and the fibers split into 5% and 95%, with the lower intensity being directed to the transmission (Tr) channel and higher to the fluorescence (Fl) channel. In detection, each signal is collimated coming out of the fiber and filtered, with the fluorescence channels using a 650nm long pass interference filters in front of each PMT. Hamamatsu H7422P-50 PMTs are used for single photon counting and these generate a single analog pulse for each detected photon, on each of the 11 channels (including one reference channel). The pulses are amplified by approximately 26 dB, to ensure they fall within the suitable analog range of operation for the single photon counting instrumentation (SPC-134, B&H, Berlin Germany), and directed to routers (HRT-41, B&H, Berlin Germany) which are used to sequence more PMTs into a single counting module.

The underlying hardware for time-correlated single photon counting (TCSPC) is highlighted in Fig. 1(a) . This instrumentation is optimal for low-light applications such as diffuse optical tomography as it allows for more sensitive detection than the conventional continuous wave or frequency modulated approaches. The use of multiple TCSPC PMTs here allows for the parallel acquisition of fluorescence and transmission measurements. In this approach, individual photon arrivals are counted and binned into a temporal histogram. The 10% quantum yield of these PMTs limits the detection rate such that not all photons are detected, but the total cumulative signal is proportional to the detected intensity. This data is then used to generate a temporal point spread function (TPSF), or more simply a histogram of photon transport time through the medium under investigation. The fundamental ability to count photons using TCSPC dictates that a very fast repetition rate source be used in a low light level application to ensure that only a single photon be detected per period. As a rule of thumb it is suggested that light range be adjusted such that the probability of detecting a single photon be 0.1-0.01 per period. This histogram of data then requires that all detected photons are counted for each of N periods and then summed together for the final data set, shown in Fig. 1(b). Many papers have been published on this system, and the foundations can be read there [24].

 

Fig. 1 (a) The underlying principal of TCSPC is illustrated, adapted from Ref. [24], where individual photon pulses analog amplified, and time referenced from a separate PMT signal. The time to amplitude conversion (TAC) is used for analog to digital conversion (ADC) and the temporal data points are summed together and used to construct a histogram of the photon arrival time in TCSPC data memory. In (b) a typical TCSPC data set is shown where the secondary peaks are artifacts from the pulsed laser shape. Note that the y-axis is logarithmic, so that the first peak is actually substantially larger than the 2nd pear. After 500 on the x-axis, the signal is dark noise.

Download Full Size | PDF

Custom control software for the fluorescence tomography system was created using the National Instruments LabVIEW development environment. Several modules were developed to control the individual automation of each sub-system. Upon starting the system a series of semi-automated initialization and calibration steps occurs, as described in the following paragraphs.

Attenuator position calibration is done each time the system is powered up, where the motorized attenuator must be homed to the zero position. This routine communicates with the motorized attenuator’s optical position encoder in a closed loop fashion. Homing ensures that the device will be applying the appropriate amount of attenuation to the source during a scan and proper homing is very important for data calibration. Filter slide calibration occurs with fluorescence and transmission filters that are controlled by two individual motors, set to “home” at a closed or “blocking” position between the other two filter sets.

The Single Photon Counting (SPC) instrumentation starts up by initializing the cards and assigning control numbers and then determines if they are functioning properly, using custom software. As part of this the DCC Cards are initialized. These are 6 PCI-based DCC-100 cards used in the SPC system function as the power supplies controlling the gain of each of the 11 PMTs. These provide software protection against PMT overload and the subsequent damage in the event saturation occurs. The final step in the calibration routine entails verifying the proper setup of the other components. When setup is completed a single reference TPSF is acquired and displayed as a final verification, with a typical data set as shown in Fig. 1(b).

2.2 Automatic exposure control

In recent development, the system’s data acquisition routine was improved with an automatic exposure control (AEC) software routine, which enables high SNR, for arbitrary objects that are imaged. This is important because the non-contact imaging of small animals and correspondingly the imaging of irregular geometries, results in signal intensity changes over a many orders of magnitude in both fluorescence and transmission intensities. To ensure an adequate SNR will be maintained for each source position, the AEC routine was developed to target the intensity within the linear operating range of the PMTs.

The specifics of the AEC routine are highlighted in the flow chart depicted in Fig. 2 . Here, the motorized attenuator that controls the source intensity is initially set to attenuate the source by 25 dB for each of the 32 positions of laser source. Using a fixed integration time of 1 second, TPSFs are then acquired for the 10 channels. This continues sequentially, with an incremental reduction in the source attenuation, until the total counts exceeds 300 on the most sensitive channel. This value of 300 counts was chosen as at this level the standard error is 5%, and this was chosen as a lower limit on a reliable signal for use in tomography. Lower values could be chosen, but lead to less reliable images. At this point the attenuation is then reduced by another 2 dB and the integrated intensity for the same channel is again recorded. Using these two data points, a linear fit is performed and used to calculate the attenuation necessary to target 50% of the PMTs maximum linear counts.

 

Fig. 2 Flow chart illustrating the logical flow of steps in the automatic exposure control (AEC) process, to ensure the TPSF has optimal signal. In the lower move attenuator box, the target counts are subtracted by the measure counts, b, and divided by a known slope value, m, which relates counts to attenuator position.

Download Full Size | PDF

The source is finally adjusted to the optimum attenuation level and the final TPSF datasets are then acquired for the current source position. The procedure is repeated for each angle of the gantry. This approach was utilized to avoid issues related to PMT saturation as this is one of the greatest challenges in using the TCPSC instrumentation to perform non-contact imaging of an irregular geometry. As a benchmark, results illustrating the performance gains that can be achieved using the AEC algorithms are illustrated in the figures presented.

2.3 Data calibration: compensation for experimental limitations & model-fit

Following data acquisition, intensity calibration for each TPSF is achieved automatically through LabVIEW and by first integrating the detected counts for a total intensity value [16,17], using an approach illustrated in Fig. 3 . Each signal is corrected to subtract out the number of dark counts as well as normalize for changes in integration time and attenuation levels at each channel, as is listed below:

 

Fig. 3 Flowchart illustrating the steps required to generate images following data acquisition. The choice of calibration method, either difference-based or calibration simulating the in vivo situation where there is no data set which can be used as an exact reference. The choice of calibration approach affects the accuracy in recovery.

Download Full Size | PDF

Here, ${\phi }_{raw}^i$ represents the calibrated intensity measurement for a single channel and ${\phi }_{TPSF}^i$ is the fluence in counts, as detected by the TCSPC instrumentation. The dark counts for each channel are computed and subtracted from the integrated counts, ${\phi }_{TPSF}^i$. This calculation is performed for each detector and then scaled by the integration time and attenuation level, as applied to the source (OD_Source) and detection channels (OD_Detection).

Unfortunately the fluorescence signal of the system is not perfect, because interference optical filters do not reject all out of band light, and there is always some leakage signal from the transmission wavelengths which makes it through the filter [28]. This is typically less than 2 orders of magnitude attenuated, but if the transmission is large, and fluorescence is weak, than this can still be the most limiting factor in the system. This leakage is measured and can be effectively subtracted out of the signal by pre-measurement of a rejection factor for each filter, under the assumption that this is repeatable and stable. Specifically, the fluorescence crosstalk signal for a given fluence can be estimated based on the transmission signal intensity and the previously determined rejection for the fluorescence filters. This corrects for intensity errors in the signal, as follows:

Here, i denotes the measurement number and the subscripts fl and tr refer to the fluorescence and transmission signals respectively. OD_{rejection_fl} is the previously determined rejection of the filters and is considered on a channel specific basis. Much like the difference calibration the subscript hetero denotes heterogeneous (post-contrast) fluorescence data. It is significant to note that bias error is introduced into this adjusted fluorescence data, ${\phi }_{fl}^i$, calibrated using this method as the intrinsic tissue autofluorescence will not be removed.

To complete the data calibration, the measurements must be normalized and adjusted to better match the modeled data from the forward diffusion model, which is then used for image reconstruction. To do this, the adjusted fluorescence data, ${\phi }_{fl}^i$, is divided by the heterogeneous transmission dataset at the excitation wavelength,${\phi }_{meas(hetero\_tr)}^i$ and finally scaled to the transmission fluence generated by the forward model, ${\phi }_{calc(hetero\_tr)}^i$, to yield a calibrated fluorescence dataset, ${\phi }_{cal\_fl}^i$.

The values of ${\phi }_{calc(hetero\_tr)}^i$ are previously estimated on a set of cylindrical homogeneous phantoms with known optical properties, such that the calibrated data matches the forward model simulations to a good level of accuracy. Then it is this set of values, ${\phi }_{fl}^i$, which is used for reconstruction. This process has ample opportunity to produce bias errors, because of all the multiplying factors which go into the process. Despite this, phantom and animal reconstruction accuracy has been quite good in previously reported work [16]. Assessment of the effective noise introduced by this calibration proceedure is completed in this paper.

2.4 Model-based image reconstruction

The process of image reconstruction with diffuse tomography data has been published in many previous papers, using the NIRFAST algorithms [15,25–27]. The MicroCT imaging allows creation of customized finite element meshes to be generated for use in the data calibration and subsequent image formation for any objects, including irregular shapes from mice. The intensity calibration is performed automatically in acquiring data with the FT system and is used to perform either difference or model-based calibration. Finite-element model based image formation is achieved using NIRFAST, the research stage image reconstruction tool developed specifically for frequency domain tomographic image reconstruction, which is an iterative, non-linear solver, using diffusion light transport in tissue. The fluorescence yield is reconstructed based upon a previously published approach with NIRFAST [16,17,25,26].

In order for the sources and detectors to be placed in the model with very high accuracy it was necessary to ensure that a relationship between the center of rotation for the FT system and the microCT was established. After scanning the phantom with the MicroCT, an image of the entire 81.8mm x 81.8mm FOV is reconstructed. This produces DICOM images which are read into Mimics software (Materialize Inc). Binary masks are created using Mimics outlining the shape of the exterior, and when interior information is to be used, this is also output in the mask. Meshes are created with a custom program, and the sources and detectors are placed in a circle within the mesh origin at the previously determined center of rotation. Using a modified version of the NIRFAST each source and detector position is adjusted radially from this initial location to the corresponding boundary location. This approach allows the sources and detectors to be placed in an automated manner, with high accuracy.

Reconstruction is done with the calibrated data, in Eq. (3), using a fixed regularization procedure where the algorithm converges until a pre-determined accuracy is met in the objective function error value. The regularization approach used was to use a scalar value to multiply into the maximum value taken from the set given by the Hessian matrix diagonal values. As such that this regularization parameter was a single number in each iteration, and could decrease in value with subsequent iterations, as the maximum diagonal value decreased. Previous studies have found this process to be optimal in producing accurate recovery of known phantoms [15,26].

2.5 Phantom studies

The system was tested with cylindrical 25 mm diameter phantoms made from ceramic resin, embedded with titanium dioxide scatterer and ink as absorber, as described in previous work [17]. In one phantom used, a 8 mm diameter hole was added at a location 4 mm in from one surface, and this was filled with an intralipid solution which matched the scattering of the background which was absorption coefficient μ_a = 0.01 mm⁻¹ and reduced scattering coefficient, μ_s ^/ = 1.0 mm⁻¹. The protoporphyrin IX (PpIX) was titrated into this at varying concentrations. The phantom was imaged repeatedly with a range of concentrations. The reconstructed images are of the fluorophore absorption coefficient, μ_af, multiplied by the fluorescent quantum yield, η. This combined parameter is directly proportional to the fluorescence yield.

The effective system noise was estimated with and without the automatic-exposure control (AEC) feature. The primary difference between the data acquired with the two imaging sequences is that the same source strength was utilized in performing difference imaging of the data acquired without the AEC. To this end, no additional error resulting from the calibration of the motorized attenuator was introduced in calibrating these datasets. To perform a direct comparison between the two techniques, a serial dilution experiment was first performed over several concentrations, and the resultant image intensity in the ROI was used to generate a calibration curve for the recovered concentration.

3. Results

3.1 Automatic exposure control validation

The AEC was designed to ensure that an adequate signal level could be obtained for each source position by maintaining the fluorescence and transmission intensity within the linear operating range of the PMTs. To assess the improvement gains resulting from the implementation of the automatic control technique, a comparison was made between the “AEC” and “No AEC” data acquisition methods for a fixed PpIX concentration of 62.5 ng/ml. This data is shown in Fig. 4(a) , calibrated in 4(b), and with reconstructed images in 4(c). This concentration was selected as it was previously determined to be the point at which the standard control scheme became unreliable. Below this concentration, the data was sufficiently corrupt that image reconstruction would not converge. Additionally, without the AEC control scheme in place, 1ug/ml was the approximate upper bound for imaging in the setup being considered due to data exceeding the linear operating range of the PMTs, or even causing PMT saturation to occur. Interestingly the image in 4(c) with the AEC is larger and more physically accurate to the true size of the region in the phantom. The recovered values for the No AEC and with AEC are shown in the plot in Fig. 4(d). The no AEC data could not be recovered outside the range shown with the circles on the plot.

 

Fig. 4 In (a) the calibrated measured data are shown for both techniques (with and without the AEC) and clearly illustrates a breakdown in the SNR, without AEC. In (b) the AEC calibrated data and the reconstructed fit data are shown, illustrating a good fit in the reconstruction. In (c), the performance comparison between the “No AEC” and “AEC” data sets, reconstructing a single fluorescent region in the top center of the phantom with PpIX. This was at 62.5 ng/ml concentration, the previously determined floor with the “No AEC” technique. In the latter case, the recovered region at top is larger, and more representative of the true region size. In (d) the reconvened values in the region for different concentrations of PpIX is shown, with both AEC and no AEC in the data acquisition. For the no AEC data, useful images could not be obtained at concentrations above and below the range shown in the plot (circles).

Download Full Size | PDF

The specific detectability limits are of course a function of the optical properties of the medium, target size, and location within the domain. This analysis was merely designed to provide a direct comparison between the two techniques. Comparing the recovered values over the wider range of PpIX concentrations (10 µg/ml – 0.2µg/ml), the results indicate that the AEC allows a linear response to be achieved over an additional order of magnitude in terms of concentration. The lower linearity bound has been extended beyond what it was previously, though it clearly breaks down somewhere between the 0.039 µg/ml and 0.2µg/ml. It is significant to note that the upper bound of this plot could be extended substantially, but was capped at 10 µg/ml as this is an upper limit regarding the biologically feasible contrast that can be obtained in-vivo with PpIX.

3.2 Effective system noise estimation

Twenty datasets (N = 20) were acquired at a concentration of 1μg/ml for each of the two techniques and used to recover 40 fluorescence yield images. A representative set of images can be found in Fig. 5(a) and 5(b) for data collected without and with the AEC respectively. The recovered fluorescence yield in the ROI was recorded for each image and used to calculate the effective Pp-IX concentration, based on the substitution of this value into the experiment specific linear fit to the calibration curve relating recovered fluorescence yield and Pp-IX concentration. The reconstructed values of standard deviation were in the range of 12-14%.

 

Fig. 5 Experimental sets of 20 images were reconstructed from phantom data obtained without AEC and with AEC, and representative images from each are shown in (a) and (b), respectively. Estimating the error in recovering the fluorescence yield could then be compared to simulations done with a range of forward data error values. In simulations, repeated reconstructions were used to estimate the standard deviation in the recovery, which is plotted in (c). This graph shows the concentration error vs. forward modeled data noise that was used in the reconstruction, and can be used to estimate the effective system error based upon the reconstructed concentrations (shown by arrows).

Download Full Size | PDF

In transforming recovered fluorescence yield into recovered Pp-IX concentration, the effective system noise can be estimated following simulated data analysis of data with known levels of generated random noise. Thirty simulated fluorescence forward data was generated with 1%, 2%, 5%, and 10% Gaussian noise added and used to reconstruct images of fluorescence yield. The standard deviation in recovered concentration at each noise level was then used to generate another calibration curve representing the simulated system noise as a function of the recovered concentration error in the ROI as, as shown in Fig. 5(c). By performing a linear fit to this data it is then possible to know how the reconstructed parameter error is related to the forward data error in a general way. As shown by the arrows on Fig. 5(c), by extrapolating backwards from the experimentally recovered standard deviations value in the Pp-IX concentration, it is then possible to estimate the noise level of the system. Based upon the reconstructed phantoms, this was estimated to be 6.9% and 7.6% for reconstructions performed using data without and with the AEC routine respectively.

3.3 Model-based calibration analysis

The approach for data calibration does not incorporate difference imaging, but rather uses the transmission data to subtract out the influence of filter imperfection in the signal. While this is not an ideal approach, with filter based signal detection, this is an unfortunate necessary step in the process. Though this approach does not require a pre-contrast image, it does substantially limit the detectability, as random bias errors are inherent in the data set through uncompensated for variations in the crosstalk signal. In order to examine the magnitude of effect that this data calibration approach induces, the data was reconstructed in two approaches, as shown in Fig. 6 .

 

Fig. 6 A systematic series of reconstructions are shown with decreasing concentration of PpIX in the top region (a). Then the recovered values for images are shown in (b), generated from (i) simulated data with appropriate amount of system noise (ii) the difference calibration technique and (iii) the in-vivo calibration routine used with animals (IVC).

Download Full Size | PDF

Image results are presented in Fig. 6, and clearly illustrates a breakdown in performance beginning at low levels, approximately 0.63 µg/ml and transforming into ring artifacts with decreasing PpIX concentration. This breakdown translates into a deviation from linearity that worsens with decreases in PpIX concentration as shown. Here, the “simulated” and “difference” quantification results are in close accordance, again reconfirming the previously determined effective system noise. Inspection of this plot also confirms the ability to quantify the recovered ηµ_af is degraded using calibrated data. This effect would be compounded by reductions in target size, irregular geometries, and intrinsic optical property changes unaccounted for – such as what might be expected using a CCD based approach and the sequential acquisition of fluorescence and transmission data.

4. Discussion

The system analyzed here is highly sensitive in that it is based upon single photon counting in the near-infrared. From this standpoint, it has the most sensitive detection approach which could be achieved. This is critical for a non-contact based imaging system, as the loss in light from the lens-based pickup is substantial. Direct contact fiber based systems could achieve better light detection, but suffer from the contact affecting the tissue being imaged. There are also problems in accurately modeling contact-based systems which make the attraction of non-contact imaging quite important. The TCSPC approach provides high sensitivity but must be carefully designed to avoid saturation and also to allow detection within the useful part of its dynamic range. In this paper, an automated exposure control is developed and demonstrated as to how it improves recovery and widens the useful concentraition range which could be detected. Analysis of the images and data in Fig. 4 indicates that the most important benefits come at the higher concentration ranges, in that the useful detection limit is increased by an order of magnitude. Similarly, the data would support the conclusion that improvements in the lower limit detection are also present.

The effective system noise was estimated by comparison with simulations, where the input random noise in the forward data could be systematically increased, to create a graph where the relationship between noise in the forward data and recovered error in the fluorescent yield was known [see Fig. 5(c)]. Using this approach, the effective noise in the system appears to be near 7% if taken to be all randomly distributed. The application of the AEC reduces this noise level slightly, although the major benefit of the AEC appears to be in widening the effective recovery range of concentration values.

The data acquisition portion of the system allows data collection with 32 source positions corresponding to 160 measurements (32 x 5). In the current configuration, it takes 10 minutes to collect a single slice dataset with the AEC being utilized at each of the 32 source positions. It is projected that the acquisition time for a single slice could possibly be reduced to 4 minutes and 30 seconds through a series of optimizations and improvements including.

5. Conclusions

This revised system presented here has improvements to the basic function which extend the dynamic range and improve the accuracy of image reconstruction. Imaging 1 cm sized regions which have PpIX in them with the system allows recovery of concentrations down to near 50 ng/ml, which is useful enough for many tumor-imaging applications. Applications of this system can be determined from these bounding limitations, where imaging must be for applications where 4 minute acquisition is suitable, and where the property recovery to near 7% accuracy is suitable, and well above the lower detection limit of 50 ng/ml concentration. While the AEC and calibration approaches introduced here are not perfect, they do provide quantitative improvements to the system performance which are imperative for the system to be useful for mouse imaging where the range of fluorophore concentrations can be quite wide and the accuracy in quantification can make a difference in successful monitoring of tumor changes over time.