Image-guided diffuse optical fluorescence tomography implemented with Laplacian-type regularization

A promising method to incorporate tissue structural information into the reconstruction of diffusion-based fluorescence imaging is introduced. The method regularizes the inversion problem with a Laplacian-type matrix, which inherently smoothes pre-defined tissue, but allows discontinuities between adjacent regions. The technique is most appropriately used when fluorescence tomography is combined with structural imaging systems. Phantom and simulation studies were used to illustrate significant improvements in quantitative imaging and linearity of response with the new algorithm. Images of an inclusion containing the fluorophore Lutetium Texaphyrin (Lutex) embedded in a cylindrical phantom are more accurate than in situations where no structural information is available, and edge artifacts which are normally prevalent were almost entirely suppressed. Most importantly, spatial priors provided a higher degree of sensitivity and accuracy to fluorophore concentration, though both techniques suffer from image bias caused by excitation signal leakage. The use of spatial priors becomes essential for accurate recovery of fluorophore distributions in complex tissue volumes. Simulation studies revealed an inability of the “no-priors” imaging algorithm to recover Lutex fluorescence yield in domains derived from T1 weighted images of a human breast. The same domains were reconstructed accurately to within 75% of the true values using prior knowledge of the internal tissue structure. This algorithmic approach will be implemented in an MR-coupled fluorescence spectroscopic tomography system, using the MR images for the structural template and the fluorescence data for region quantification. © 2007 Optical Society of America OCIS codes: (170.0110) imaging systems; (170.6280) spectroscopy, fluorescence and luminescence; (170.6960) tomography; 170.3010 image reconstruction techniques; (170.1610) clinical applications.


Introduction
Imaging the spatial distribution of fluorescence activity at depth in tissue is a challenging problem that has yet to impact clinical practice.The most popular approach is diffuse optical fluorescence tomography (DOFT), a model-based method which approximates photon propagation as a diffuse field and computationally matches model parameters to measured boundary data.The theoretical framework is well established [1][2][3][4][5] and feasibility studies have demonstrated fluorescence yield imaging in various phantom geometries [6][7][8][9].Much of the most recent research has focused on high-throughput, full body imaging of small animals [10,11] but no human images have been published to date.DOFT produces low resolution images compared to standard clinical modalities such as x-ray CT and MRI due to the highly scattered photon fields and relatively sparse measurement sampling of the tissue volume.The reconstruction problem is ill-posed and underdetermined and large heterogeneity in tissue optical properties caused by complex tissue morphology further challenges the imaging algorithm.Sensitivity drops with increasing depth resulting in a non-linear responsivity across the imaging field, making it more difficult to image larger tissue volumes.Experiments in larger volumes that more closely resemble a human breast typically consider unrealistically high fluorophore contrasts and simple geometries [9,12].Though these studies are important for advancing the understanding of the modality, improving image resolution and contrast sensitivity is critical for identifying a clinical role for DOFT.
DOFT is an extension of the more widely studied diffuse optical tomography (DOT) which suffers similarly from low resolution and depth-dependent contrast sensitivity in the absence of spatial or spectral prior information to guide the solution.However, methods to incorporate highly-resolved anatomical data obtained from standard clinical modalities have improved the quantification of the optically derived images [13][14][15].These hybrid approaches lead to a conceptually new application of optical tomography, one in which the highly resolved imaging system provides a structural template upon which volumetric optical spectroscopic images are constructed.This framework may be applied to either absorption and scatter spectroscopy, or fluorescence spectroscopy.Additional challenges specific to the fluorescence case include lower signal intensity, excitation source contamination of the fluorescence emission measurements due to filter leakage, and tissue optical property effects on both the excitation and emission photon propagation.Applying spatial guidance techniques to the more complicated DOFT problem may yield even larger gains in imaging capability.
This paper introduces a method to dramatically improve fluorescence imaging at depth in tissue by coupling MRI and fluorescence tomography.Tissue structural information determined from standard T1 and T2 MR images is encoded as a spatial filter in the DOFT reconstruction algorithm and used to guide the recovery of fluorescence activity.In this implementation, reconstruction parameters are loosely grouped into regions based on tissuetype determined from the MR images but are permitted to update independently, giving rise to the term "soft" spatial priors [16].The algorithm is tested with phantom data of Lutetium Texaphyrin, a near-infrared photosensitizer shown to accumulate preferential in a variety of malignancies [17][18][19][20], recorded from a newly developed MR-coupled spectroscopy-based fluorescence tomography system.A simulation study based on realistic geometries generated from MR images of a normal human breast serves as an initial illustration of expected improvements provided by the spatial priors approach in complex domains.

Theory
Assuming highly diffuse media, the transport of light at a given modulation frequency ω in the presence of fluorescence generated by an external source at the excitation wavelength (λ x ) of the fluorescing agent, is modeled by two diffusion equations, where the solution of the first equation provides the driving source term of the second [5]: where subscripts x and m represent the excitation and emission fluence at wavelengths λ x and λ m , respectively.The intrinsic optical parameters The most appropriate description of the air-tissue boundary is derived with an indexmismatched type III condition (also known as Robin or mixed), in which some fraction of the fluence at the external boundary of the tissue exits and does not return.The flux leaving the external boundary is equal to the fluence rate at the boundary weighted by a factor that accounts for the internal reflection of light back into the tissue.This relationship is described in the following equation: where ξ is a point on the external boundary, and A depends upon the relative refractive index (RI) mismatch between tissue Ω and air.A can be derived from Fresnel's law: ( ) , the angle at which total internal reflection occurs for photons moving from region Ω with RI n 1 to air with RI n AIR , and

Finite element implementation:
When the refractive index is homogenous (assumed to be 1.33 [21]), the finite element discretization of a volume Ω can be obtained by subdividing the domain into D elements joined at V vertex nodes.In the finite element method (FEM),

Φ
becomes one of sparse matrix inversion, and in this work, a bi-conjugate gradient stabilized solver is used.As developed previously [22,23], the coupled diffusion Eqs. ( 1) and (2) in the FEM framework can be expressed as a system of linear algebraic equations: where the matrices K x,m (κ , and F x,m have entries given by: and the source vector Q 0 has terms The source term is defined as a Gaussian distribution, matching the intensity profile at the tip of the optical fiber.Because the source is assumed spherically isotropic, modeling is more accurate when it is centered one scattering distance within the outer boundary.The source vector Q m for fluorescence re-emission is expressed as and is distributed throughout the domain.

No spatial priors
In the inverse problem, the goal is the recovery of optical properties at each FEM node using surface measurements of light fluence at both the excitation and emission wavelength sequentially (assuming the use of two externally applied sources, one at each wavelength), followed by fluorescence yield reconstruction at the emission wavelength.The computational approach is to minimize the difference between measured fluence,  ( ) where NM is the total number of measurements given by the imaging system, NN is the number of parameters representing the optical property distribution which corresponds to the number of nodes in the reconstruction mesh, and I is an NN × NN identity matrix.In general, [ ] ( ) where J is the Jacobian matrix, here calculated using the Adjoint method [24].Equation ( 13) is known as the Moore-Penrose generalized inverse and is found to be suitable for problems where the number of unknowns to be recovered is much larger than the amount of information (measurements) available.In standard practice, I is an identity matrix, and in this work λ is some fixed fraction multiplied by the maximum value on the diagonal of the Hessian matrix J J T , and is therefore updated at each iteration.To recover the optical properties at the emission wavelength,

(
) , the externally applied illuminating source is changed to one at the emission wavelength and the formulation presented in Eq. ( 13) is used.
The recovered optical properties are used in the fluorescence yield reconstruction which also adheres to the minimization formulation presented in Eq. (12).The unknown parameter in Eq. ( 13) becomes ) (r f a ημ , and the Jacobian can be calculated by similar Adjoint properties described above and in Ref. [5].

Spatial priors
Spatial prior information is incorporated by assuming a 'generalized Tikhonov' penalty term which is similar in structure to Eq. ( 12) except that the identity matrix is replaced with a Laplacian-type matrix, presented here for the excitation field: where NN is the number of unknowns in the model [16,25,26].The constant, β , balances the effect of the parameters with the model-data mismatch in the same manner as λ in Eq. ( 13).The dimensionless 'filter' matrix, L, is generated using MRI-derived priors and its construction is flexible.In this application, each node in the FEM mesh is labeled according to the region, or tissue type, with which it is associated (in the MR image).The L-matrix represents a Laplacian-type structure, the diagonal of which is L i,i =1 where i is the nodal index.When nodes i and j are in the same region containing n nodes, L i,j =-1/n, otherwise L i,j =0.This effectively relaxes the smoothness constraints at the interface between different tissues, in directions normal to their common boundary.The effect on image quality is similar to that achieved through total variation minimization schemes [27] but easily encodes internal boundary information from MR images.Following the same procedure as in the no-priors case, the parameter update is given by [ ] ( ) where much like λ in Eq. ( 12), β is a fixed fraction multiplied by the maximum value on the diagonal of J J T .This update formulation is also used to recover the optical properties at the emission wavelength as well as the fluorescence yield, as described in Section 2.2.1.

Reconstruction basis
A number of different strategies for defining the reconstruction basis are possible, including second mesh and pixel basis [28,29].The choice of reconstruction basis allows for computational efficiency which serves to reduce the number of unknowns in the algorithm.The problem at hand is twofold: The forward problem requires that the volume of interest is subdivided into adequate number of sub-domains which allow for an accurate description of the calculated fields, whereas a reduction in the number of unknowns improves the illposedness of the problem in the reconstruction algorithm.This is addressed by defining a separate reconstruction basis (different from the meshes used in the FEM implementation), upon which the unknowns are updated in Eqs. ( 13) and ( 15) [23].In cases where no prior structural information is available, a pixel basis which defines a set of regularly spaced pixels for the update of the quantities of interest is used.However, when spatial priors are involved, it is crucial to ensure that enough pixels are used to adequately define small structural regions.In this case, a set of regular pixel bases are introduced in each region of interest, which account for the region's individual shape and size.A semi-adaptive method which allows the number of pixels in each region to be selected was used in the studies presented here.The same region-based pixel basis is used throughout the iterative reconstruction process.

Simulation studies
Test domains for simulation studies were derived from a T1 weighted MR image of a human breast which measured approximately 10.5 cm in diameter.A 2-D slice of the breast volume was discretized into a mesh of approximately 2000 nodes and MR image intensity thresholds were used to assign adipose and fibro-glandular tissue volumes into mesh regions.Figure 1 shows the original MR breast image and its associated discretized, region-labeled mesh.The MR image originated from a clinical exam with our MR coupled frequency domain NIR system and therefore reveals an irregularly shaped breast boundary caused by the fiber optic array slightly compressing the breast at the contact positions.A cancerous tumor region was added as a target anomaly or region of interest for these studies and is depicted in the figure .A second test case with a tumor region located near the center of the domain was used to explore the nonlinear sensitivity of diffuse optical tomographic techniques.The sourcedetector positions were determined directly from the MR image and represent 16 optical fibers circumscribing the breast domain in a single plane.Light is detected at all non-source fiber positions providing a total of 240 measurements of intensity and phase per wavelength.This simulated configuration matches our experimental fluorescence tomography system.In this set of studies, regions were assigned values in terms of biologically relevant optical absorbing (HbO, dHb, Water, Lutetium Texaphyrin) as well as scattering parameters (Mie scattering amplitude and power).These values, provided in Table 1, represent typical known in vivo levels of endogenous chromophores and are consistent with previous clinical work [30,31].Though it is unknown exactly how Lutetium Texaphyrin (LuTex) will distribute in a human breast, the concentrations used are similar to those reported in ex vivo studies [17] and at least represent a reasonably complex distribution for demonstrating recovery of fluorescence yield.Tissue chromophore concentrations were used to calculate total optical absorption and scattering values at the excitation and emission wavelengths based on experimentally determined values of molar extinction coefficients.Simulated noisy data was generated based on these optical properties in the following manner: 1) frequency domain data for a light source at the excitation wavelength, 2) frequency domain data for a light source at the emission wavelength, and, 3) Continuous Wave (CW) fluorescence emission data at the emission wavelength.This represents a total of 3 data sets for a given imaging session: two frequency domain data sets for determining background optical properties and one CW fluorescence emission data set.CW fluorescence emission data was used to match the capabilities of our experimental system and results in a simplification of Eq. ( 2) which can be handled by setting ω = 0. Given the modest Stoke's shift of LuTex, depicted in Fig. 2, the choice of excitation wavelength has practical implications for experimental work.The difficulty in filtering the excitation signal precludes the use of Lutex's NIR absorption peak (about 735 nm) to excite the fluorophore.In this study, a 690 nm excitation wavelength was used for both simulated and experimental data acquisition.In addition to normally distributed random noise added to the frequency domain data (5% amplitude and 1 degree phase) and CW fluorescence emission (10% intensity), excitation signal leakage through the filter was added to the CW fluorescence emission intensity to simulate a typical 7 OD rejection of the excitation intensity.This number comes directly from experimentally measured rejection estimates for the filters in the tomography system used in the experiments performed here.The general image reconstruction protocol was as follows: 1) Reconstruct for optical properties at the excitation wavelength, μ ax and μ sx ', with frequency domain data, 2) Reconstruct for optical properties at the emission wavelength, μ am and μ sm ', with frequency domain data collected using a laser source at the emission wavelength, 3) Use the reconstructed optical properties and fluorescence intensity data to recover fluorescence yield.The same reconstruction algorithm was used to determine background optical properties in steps ( 1) and ( 2) and is based on previously reported work [32,33].Initial estimates for all parameters were generated using homogenous fitting algorithms which enforce a single value for all nodes.The Jacobian matrix was calculated on a fine mesh of approximately 2000 nodes and interpolated onto a course reconstruction pixel basis for inversion.A 30 by 30 pixel reconstruction basis was used for the no-priors case.The spatial priors reconstruction used a newly developed semi-adaptive pixel basis that redistributes the pixels based on the region information, as described in section 2.2.3.This method ensures that each region contains an adequate number of nodes to approximate the internal structure of the domain.Convergence was defined as less than a 2% change in projection error between successive iterations for the frequency domain optical properties algorithm and less than a 1% change in projection error between iterations for the fluorescence yield reconstruction.Similar algorithmic parameters were used for the phantom studies and are described in further detail below.

Phantom studies
A spectrometer based tomographic imaging system which couples directly into a Philips 3T MRI magnet was used to acquire fluorescence emission spectra.The system, depicted in Fig.  Lutetium Texaphyrin provided by Pharamcyclics was diluted in water and used as the imaging fluorophore.The test domain was a solid 5.5 cm diameter hardened epoxy phantom with scatterer and absorbers created by titanium dioxide powder and India ink, respectively [34,35].The phantom had a 14 mm hole located approximately 12 mm from the phantom center.The background optical properties were μ ax = 0.005 and μ sx ' = 1.0 mm -1 , measured with a frequency domain system near the excitation wavelength.Unlike the simulation experiments, the optical properties were assumed constant throughout the domain in this experiment.These values were also used as the optical properties at the emission wavelength.The hole was filled with a solution of 1% Intralipid to match the scattering value of the background, and varying concentrations of Lutex (0.3125 μM to 5 μM) were added.This represents a simple test case for investigating the imaging response to varying concentrations of fluorophore.The excitation source was a 690 nm laser diode which matches the wavelength used in the simulation studies.Total acquisition time for the fluorescence emission was less than 4 minutes (total of 240 data points).
Even after processing the collected light with a 720 nm long pass interference filter (Omega Optics) which provides 7 OD rejection of the excitation light as well as the filtering offered by the spectrograph grating, emission spectra recorded by the detector are composed of a sum of the pure fluorescence signal and excitation bleed-through.In order to decouple these signals, previously recorded "basis spectra" of the bleed-through signal and the pure fluorescence signal are fit to the data.The process is illustrated in Fig. 4 for one measurement at a single detector.A linear least squares algorithm operating on the basis spectra quantifies the amount of excitation bleed-through and true fluorescence signal that exists in each spectral recording, further reducing the effect of excitation bleed-through.This is accomplished by minimizing the summation To investigate the improved imaging capability of reconstructing data with spatial priors, each data set was reconstructed with and without spatial "soft" priors.Since the domain was easily characterized geometrically, spatial prior information was determined by direct manual measurement of the phantom.This information was encoded in the fine mesh used for reconstruction.Since background optical properties were previously determined with our frequency domain system, they were assumed known for image reconstruction.Convergence was defined as less than a 1 % change in projection error between iterations.All images were reconstructed using a 2GHz Centrino Duo laptop with 2GB RAM running Windows XP.

Simulation results
Figures 5 and 6 display the target values of the two test cases along with images reconstructed using no priors and spatial soft prior information.Qualitatively, image recovery using spatial priors produces significantly more accurate images.Spatial priors preserve the general internal structure of the Lutex distribution, detail that is lost almost entirely in the no-priors case.Images of fluorescence yield show the most dramatic difference between the spatial priors and no priors reconstructions.Without spatial priors, the algorithm appears to have no ability to recover "cancer" regions of elevated fluorescence yield for these complicated cases.However, incorporating spatial priors results in images that qualitatively appear accurate and quantitatively are reasonably close to the true values.Figure 7 provides 1-D cross-sections near the y-axis of each domain.These plots confirm an inability of the no-priors imaging algorithm to recover the simulated fluorescent tumor in either test field.Alternatively, the spatial prior-based imaging algorithm not only picks out the objects of interest, but provides fairly accurate reproductions of the complicated structure of simulated fibro-glandular layers.Mean values for the simulated cancer region near the domain center are 79% and 51% of the true values for the spatial priors and no-priors reconstructions, respectively.These numbers change to 75% and 45% for the case with the cancer region closer to the edge of the domain.They represent a significant improvement in imaging performance; however, they alone do not illustrate the full impact of incorporating spatial priors.The cross-sectional plots indicate that the no-priors images contain virtually no spatial discrimination of the cancer regions.Regional contrasts are depicted in Table 2 and further illustrate a dramatic overall improvement in cancer region quantification with spatial priors.s,x , μ a,m , μ / s,m , and fluorescence yield, ημ af , for reconstruction implementations using no prior information and with spatial prior information.In this case, the simulated cancer region is near the edge, which is known to be easier to recover without spatial priors.The image scales are at right.In addition to improved qualitative and quantitative accuracy, the spatial prior algorithm reduced the reconstruction time significantly.The full reconstruction time for the no-priors case, including background optical property estimation, was just under 9 minutes for both the central and superficial tumor cases.These times were reduced to less than 3 minutes and 90 seconds using spatial priors.In both cases, structural information guided the algorithm to a convergent solution in far fewer iterations than the no-priors cases.

Phantom results
Figure 8 shows fluorescence yield images recovered from phantom data using both no-prior and spatial prior image reconstruction approaches.A qualitative assessment of the images reveals a dramatic benefit of the spatial prior on image formation.The fluorescent object's borders are more clearly defined for all concentrations of Lutex when using spatial priors.Furthermore, the values of fluorescence yield throughout the region of interest are more homogeneous, and therefore, more similar to the actual distribution for spatially guided reconstructions.Incorporating spatial priors also suppresses edge artifacts significantly.This is most apparent in images of phantoms with low Lutex concentration.Figure 9 shows a full scale image of the 0.3125 μM phantom.Artifacts near the boundary virtually dominate the no-priors image and represent the largest fluorescent yield values in that imaging domain.None of these artifacts exist in the spatially guided image generated from the same fluorescence emission data.The spatially guided image also includes an easily discernable fluorescent object at the correct location which is difficult to define in the no-priors case.These results indicate that for low fluorophore concentrations in this imaging domain, the correct fluorescent object would be identified only if prior structural information were incorporated in the reconstruction.Both techniques suffer from a bias in the recovered fluorescence yield which results in a positive value in the region of interest for the case with no fluorophore.This is most likely caused by bleed-through of the excitation light in the emission measurements coupled with an imperfect spectral fitting technique which results in positive values for fluorescence emission intensity even in phantoms containing no Lutex fluorophore.Incorporating additional filtering, such as band-pass source filtering, and optimizing the spectral fitting routine should address the bias signal and decrease the potential for incorrect quantification.
In addition to improving image quality and fluorescence yield quantification, spatial priorbased reconstructions converged in substantially less time than the no-priors algorithms.Reconstruction times ranged from 50 -90 seconds when not using spatial priors and approximately 22 seconds for spatially guided reconstructions, marking an improvement of just over 75% in some cases.These numbers represent only the fluorescence yield reconstructions and not the recovery of background optical properties since μ a and μ s ' were assumed known from prior measurements.In cases requiring the recovery of background optical properties, the improvement in reconstruction time will be similar to that previously quoted in the simulation results.

Discussion
This study introduced an effective method to incorporate MR-derived tissue morphology for imaging fluorescence yield at depth in tissue.The conceptual assertion is that diffuse tomography will be more successful as an imaging modality when combined with pre-existing imaging systems which have higher spatial resolution, such as MRI.Simulation and phantom studies were used here to validate the method prior to implementation as a full scale system.Using Lutex phantom data, the soft spatial prior implementation improved qualitative and quantitative imaging response of fluorescence yield.Prior information also suppressed image artifacts and more accurately represented the internal distribution of fluorophore.In the nopriors case, edge artifacts dominate the image for lower concentrations of fluorophore, and provide a misleading interpretation of the internal distribution of the fluorescent agent.The phantom studies represented simple distributions of optical properties and fluorophore concentration and improvements in image quality were still substantial.It is expected that soft prior implementations will benefit imaging performance to an even greater extent in complex tissue domains, an assertion born out in the simulation studies.
A simulated breast mesh derived directly from an MR image of a human breast served as the test bed for a complex tissue domain.The simulations demonstrated no ability to recover the internal distribution of the complicated domain without the use of spatially guided reconstructions.Qualitatively, fluorescence yield images generated without spatial priors had little resemblance to the target domain and completely disregarded the 18mm simulated cancer region, as evidenced by 1-D cross-sectional plots of fluorescence yield.Breakdown of the images in the no-priors case is likely due to poorly recovered background optical properties as well as the complexity of the fluorescence yield distribution itself.Certainly, In these studies, it was assumed that the fluorophore distribution correlates directly to the fatty, fibro-glandular, and tumor tissue layers, which themselves exhibit no significant intraregion heterogeneity, in a given imaging domain.While this assumption appears reasonable based upon the biology of the tissue for endogenous chromophores [30,31], further studies of the uptake of Lutex in-vivo are required to determine a correlation between tissue type and fluorophore localization.Should this assumption prove to be unrealistic, the "soft" spatial prior approach offers some latitude in terms of correctly identifying structural prior information.The algorithmic implementation groups tissue regions together in a regionspecific regularization and allows individual nodes in those regions to update independently.As opposed to "hard" prior approaches where nodal values in a given region are assumed homogenous, a soft prior technique may recover positive fluorescent objects not directly encoded in the spatial prior information.This is a subject of further investigation.
The experimental system introduced here couples directly into a Philips 3T MRI to provide simultaneous MR and NIR fluorescence imaging.Simultaneous imaging simplifies co-registration of the MR image with the NIR reconstruction domains and reduces overall acquisition time.Optimization of this system for in vivo imaging is underway for both small animal and human breast imaging.

Conclusion
A spatially guided image reconstruction implementation based on prior knowledge of tissue morphology was shown to provide significant improvements in fluorescence yield recovery in complicated tissue volumes and to be highly beneficial for simple domains.Specifically, both phantom and simulation results demonstrated dramatic improvements in recovery and quantification of features in the fluorescence distribution.Structural guidance also reduced image reconstruction time substantially.A newly developed experimental system couples full-volume deep-tissue fluorescence spectroscopy capabilities directly into the MR bore for simultaneous MR-NIR fluorescence image acquisition.The results presented here show promise for this approach in all cases and tissue volumes considered.Extensions of this study are underway to determine the effect of incorrectly identifying the structural prior, especially in cases where the MR images produce false negative or false positive readings.A full characterization of the imaging limits of the experimental system will further complement the results presented here.It is clear that incorporating anatomical features derived from MR images in DOFT image reconstruction will improve sensitivity to lower concentrations of fluorophore, qualitative accuracy, and fluorescence yield quantification in-vivo.

Φ
is the photon fluence rate at position r.The diffusion coefficient is given by ) c(r) is the speed of light in the medium at any point, defined by c o /n(r), where n(r) is the index of refraction at the same location and c o is the speed of light in vacuum.The fluorescence parameters are the lifetime ) a product of the fluorophore's quantum efficiency η and its absorption coefficient

Φ
, from the model Eqs.(1) and (2) by adjusting the spatial distribution of the unknown parameters through minimization of the 'objective' function.The objective function for recovering the optical properties at the excitation wavelength,

2 χ
will not equal zero, but the values of x μ for which x zero can be determined, based on an initial estimate of x μ .Following the Taylor series method for deriving Newton's method, the second and higher order terms are ignored, leading to the iterative update equation:

Fig. 1 .
Fig. 1.An axial T1 weighted MR slice of a human breast is shown in (a), from which the test domain for simulation studies was derived.Darker regions indicate fibro-glandular tissue imbedded in adipose tissue indicated by lighter values.The image was acquired during a clinical exam with one of our MR-coupled NIR tomography systems and shows the indentations caused by the fiber optic probes.The test domain (b) was a discretization of the MR image thresholded into regions.The yellow anomaly was added to simulate a targeted cancerous tumor.

Fig. 2 .
Fig. 2. The normalized absorbance and fluorescence emission spectra of Lutetium Texaphyrin are shown.

#
78716 -$15.00USD Received 4 January 2007; revised 22 February 2007; accepted 13 March 2007 (C) 2007 OSA3, is composed of 16 spectrometers with low noise CCD cameras cooled to -70 C. Sixteen 13 meter long silica fiber bundles circumscribe the imaging domain in contact mode and couple into the spectrometers via custom designed input optics with a six position automated filter wheel.Fluorescence emission signals are processed with long pass interference filters before entering the spectrometer.Each fiber bundle contains eight 400 μm fibers, seven of which directly couple the tissue surface to the spectrometer input while the eighth couples the tissue surface to the light source.This detector configuration minimizes coupling losses and provides parallel detection of full spectra for each source position.Inter-fiber tissue coupling calibration factors were determined prior to imaging using a cylindrical homogenous phantom with a centrally located source/detector.The custom manufactured input optics were designed to optimize filtering and maximize the light detection efficiency of the spectrometers.Camera exposure times can be adjusted to the detected light level to maximize the signal to noise ratio.With 15 detectors per 16 source positions, a total of 240 measurements are recorded for a given acquisition.

Fig. 3 .
Fig. 3.The experimental spectrometer-based system depicted at left couples directly into the MR via 13 meter fiber optic bundles.Sixteen spectrometers are computer controlled for rapid image acquisition (left photograph).An animal interface (right photograph) is composed of a rodent MR coil custom built by Philips Research Hamburg to accommodate the optical fiber array for simultaneous MR and NIR fluorescence imaging.

Fig. 4 .
Fig. 4.An example of a pair of basis spectra for the excitation and fluorescence light (in counts/s as a function of CCD pixel number) is shown in (b).These spectra are recorded for each detector prior to imaging.In practice, the basis spectra are used to perform a least squares fit (a) to the spectrum measured for each source-detector pair to determine the relative contribution of the fluorescence and excitation light to the measured response.

Fig. 5 .
Fig. 5. Target and recovered values of μ a,x , μ /s,x , μ a,m , μ / s,m , and fluorescence yield, ημ af , for reconstruction implementations using no prior information and with spatial prior information.In this case, the simulated cancer region is near the edge, which is known to be easier to recover without spatial priors.The image scales are at right.

5 #Fig. 6 . 5 #
Fig.6.Target and recovered values of μ a,x , μ / s,x , μ a,m , μ / s,m , and fluorescence yield, ημ af , for reconstruction implementations using no priors and spatial priors.In this case, the simulated tumor region is near the center of the imaging domain, which is known to be more difficult to recover accurately.Image scales are at right.

Fig. 7 .
Fig. 7. Cross sectional plots of fluorescence yield are shown for the simulated imaging domains in (a) the case with an object near the edge and (b) the case with an object near the center.In both cases, the cross section is in the y-direction just off center from x = 0.The solid line represents the target value, the small dotted line the recovered value using a no-priors based algorithm, and the dashed line the recovered value using spatially guided reconstruction.

Fig. 8 .Fig. 9 .
Fig. 8. Recovered images of fluorescence yield are shown for varying concentrations of Lutex.The 14 mm diameter fluorescent inclusion was embedded in a 55 mm diameter solid epoxy tissue simulating phantom.Images were generated from the same data using algorithms based on no-priors and spatial soft prior implementations.

Table 1 .
Chromophore concentrations and scattering parameter values assigned to the mesh regions in the simulation studies

Table 2 .
Target and recovered fluorescence yield regional contrasts for the images in Figs.4 and 5.Reported ratios were calculated using the mean values in each region.