xml version 1.0 encoding UTF-8 standalone no
record xmlns http:www.loc.govMARC21slim xmlns:xsi http:www.w3.org2001XMLSchema-instance xsi:schemaLocation http:www.loc.govstandardsmarcxmlschemaMARC21slim.xsd
leader nam 22 Ka 4500
controlfield tag 007 cr-bnu---uuuuu
008 s2010 flu s 000 0 eng d
datafield ind1 8 ind2 024
subfield code a E14-SFE0004595
Fluorescence molecular tomography :
b a new volume reconstruction method
h [electronic resource] /
by Stephen Shamp.
[Tampa, Fla] :
University of South Florida,
Title from PDF of title page.
Document formatted into pages; contains X pages.
Thesis (MSEE)--University of South Florida, 2010.
Includes bibliographical references.
Text (Electronic thesis) in PDF format.
Mode of access: World Wide Web.
System requirements: World Wide Web browser and PDF reader.
ABSTRACT: Medical imaging is critical for the detection and diagnosis of disease, guided biopsies, assessment of therapies, and administration of treatment. While computerized tomography (CT), magnetic resonance imaging (MRI), positron emission tomography (PET), and ultra-sound (US) are the more familiar modalities, interest in yet other modalities continues to grow. Among the motivations are reduction of cost, avoidance of ionizing radiation, and the search for new information, including biochemical and molecular processes. Fluorescence Molecular Tomography (FMT) is one such emerging technique and, like other techniques, has its advantages and limitations. FMT can reconstruct the distribution of fluorescent molecules in vivo using near-infrared radiation or visible band light to illuminate the subject. FMT is very safe since non-ionizing radiation is used, and inexpensive due to the comparatively low cost of the imaging system. This should make it particularly well suited for small animal studies for research. A broad range of cell activity can be identified by FMT, making it a potentially valuable tool for cancer screening, drug discovery and gene therapy. Since FMT imaging is scattering dominated, reconstruction of volume images is significantly more computationally intensive than for CT. For instance, to reconstruct a 32x32x32 image, a flattened matrix with approximately 1010, or 10 billion, elements must be dealt with in the inverse problem, while requiring more than 100 GB of memory. To reduce the error introduced by noisy measurements, significantly more measurements are needed, leading to a proportionally larger matrix. The computational complexity of reconstructing FMT images, along with inaccuracies in photon propagation models, has heretofore limited the resolution and accuracy of FMT. To surmount the problems stated above, we decompose the forward problem into a Khatri-Rao product. Inversion of this model is shown to lead to a novel reconstruction method that significantly reduces the computational complexity and memory requirements for overdetermined datasets. Compared to the well known SVD approach, this new reconstruction method decreases computation time by a factor of up to 25, while simultaneously reducing the memory requirement by up to three orders of magnitude. Using this method, we have reconstructed images up to 32x32x32. Also outlined is a two step approach which would enable imaging larger volumes. However, it remains a topic for future research. In achieving the above, the author studied the physics of FMT, developed an extensive set of original computer programs, performed COMSOL simulations on photon diffusion, and unavoidably, developed visual displays.
Advisor: Vijay K. Jain, Ph, D,
Diffuse optical tomography
Singular value decomposition
High spatial sampling
x Electrical Engineering
t USF Electronic Theses and Dissertations.
Fluorescence Molecular Tomography: A New Volume Rec onstruction Method by Stephen Joseph Shamp A thesis submitted in partial fulfillment of the requirements for the degree of Master of Science in Electrical Engineering Department of Electrical Engineering College of Engineering University of South Florida Major Professor: Vijay K. Jain, Ph.D. Andrew Hoff, Ph.D. Gokhan Mumcu Ph.D. Date of Approval: July 6, 2010 Keywords: fmt, khatri-rao, diffuse optical tomograp hy, singular value decomposition, high spatial sampling, overdetermined Copyright 2010, Stephen Joseph Shamp
Acknowledgements I wish to express my gratitude to my major professo r, Dr. Vijay K. Jain for his guidance and inspiration throughout this research and my gra duate studies. I also thank my committee members Dr. Andrew Hoff a nd Dr. Gokhan Mumcu for their helpful comments and suggestions. Finally, I would like to thank my wife, parents, an d family for their love and encouragement.
i Table of Contents List of Tables iiiList of Figures ivAbstract vi1. Introduction 12. Review of Literature 62.1. Experimental Techniques for Acquiring FMT Imag ing Data 62.2. Model for Photon Diffusion in Scattering Media 82.2.1. Dirac Delta Point Sources 82.2.2. Photon Diffusion Boundary Conditions 122.2.3. Sinusoidal Amplitude Modulated Point Sources 142.3. Normalized Born Field 182.4. Reconstruction Methods 203. Challenges in FMT Reconstruction 244. FMT Reconstruction Via Khatri-Rao Decomposition 264.1. Theory 264.1.1. Weight Matrix in Tensor and Array Forms 264.1.2. Decomposition of Weight Matrix to Extract No rmalizing Term 274.1.3. Decomposition of Weight Matrix to Khatri-Rao Product 294.1.4. Pseudoinverse of Khatri-Rao Product 314.2. Reconstruction Algorithms 324.2.1. Method 1: Row-Wise SVD-KR Reconstruction 334.2.2. Method 2: Column-Wise SVD-KR Reconstruction 344.2.3. Minimizing Reconstruction Errors from Noisy Measurements 36
ii 4.3. Results for the New FMT Reconstruction 405. Modeling of Photon Diffusion 475.1. Photon Propagation Model for Homogeneous Media 475.1.1. Finite Element Model 485.1.2. Method of Sources with Green's Function 495.1.3. Model Comparison 525.2. Finite Element Model for Heterogeneous Media 5 45.2.1. CT Image Segmentation 555.2.2. Finite Element Modeling of Photon Propagati on 565.3. Precession of Normalized Born Field in Heterog eneous Media 586. Administering and Imaging Multiple Fluorochromes Simultaneously 606.1. Results 637. Parallel Processing Implementations 667.1. Parallel Reconstruction by SVD-KR 667.2. Two-Stage Approach for Larger Imaging Volumes 677.2.1. Implementing SVD Region of Interest Enhancem ent 697.2.2. SVD Column Removal Update 697.2.3. SVD Column Addition Update 728. Conclusions 759. List of References 77 About the Author End Page
iii List of Tables Table 1. Optical Properties of Tissue 10Table 2. SVD-KR and SVD Reconstruction Time 42Table 3. SVD-KR and SVD Calculated Reconstruction M emory Usage 43Table 4. SVD-KR and SVD Reconstruction Errors 44Table 5. SVD-KR Reconstruction Time and Error 45Table 6. Optical Properties of Selected Tissue 58Table 7. SVD-KR Parallel Reconstruction Time and Re lative Speed 67
iv List of Figures Figure 1. Photon-Matter Interactions 2Figure 2. Simplified Jablonski Energy Diagram 2Figure 3. Simplified FMT Imaging System 3Figure 4. Reconstruction Using SVD-KR 5Figure 5. Simplified Non-contact FMT Imaging Syste m 7Figure 6. Solution to (1) Using COMSOL 10Figure 7. Photon Density vs. Time Given by (6) 11Figure 8. Model of Incident Collimated Light as Poi nt Sources 12Figure 9. Magnitude of Photon Density vs. Time for Example (31) 17Figure 10. Calculation of Weight Matrix for Simplif ied 2-D FMT System 21Figure 11. Method of Projections for Two Equations and Two Unknowns 22Figure 12. Flattening the Weight Matrix from a Tens or to an Array 27Figure 13. Decomposition of Weight Matrix to Remove Normalizing Term 28Figure 14. KR Product Decomposition of Wsd 30Figure 15. Row-Wise SVD-KR Reconstruction Algorith m 34Figure 16. Column-Wise SVD-KR Reconstruction Algor ithm 35Figure 17. Plot of Signal and Noise Contributions v s. Singular Value Index 39Figure 18. Reconstruction Error vs. Rank 40Figure 19. Diffuse Non-ellipsoidal Phantom Used in Reconstruction 41
v Figure 20. SVD-KR and SVD Reconstruction Time 42Figure 21. SVD-KR and SVD Calculated Reconstruction Memory Usage 43Figure 22. Reconstruction Time vs. RMSE 44Figure 23. Reconstruction Using SVD-KR 46Figure 24. Imaging Chamber Schematic for COMSOL Sim ulation 48Figure 25. Source Boundary Conditions 50Figure 26. 2-D Method of Sources to Enforce Boundar y Conditions 51Figure 27. Surface Fit of Finite Element Model by M ethod of Sources Model 53Figure 28. Residuals of Surface Fit 53Figure 29. Segmented Digimouse Skeleton Showing Reg ion of Interest (ROI) 56Figure 30. Segmented Mouse Leg 56Figure 31. Imaging Chamber Schematic for Heterogene ous COMSOL Simulation 57Figure 32. Synthetic Intensity vs. Wavelength Measu red at Detector 61Figure 33. Multi-fluorochrome Imaging with Contrast Enhancement 65Figure 34. SVD-KR Parallel Reconstruction Relative Speed 67Figure 35. Mesh Refinement in Region of Interest 68Figure 36. Removing Columns from the Weight Matrix 70Figure 37. Appending Columns to the Weight Matrix 73
vi Fluorescence Molecular Tomography: A New Volume Rec onstruction Method Stephen Joseph Shamp Abstract Medical imaging is critical for the detection and d iagnosis of disease, guided biopsies, assessment of therapies, and administration of trea tment. While computerized tomography (CT), magnetic resonance imaging (MRI), positron emission tomography (PET), and ultra-sound (US) are the more familiar m odalities, interest in yet other modalities continues to grow. Among the motivation s are reduction of cost, avoidance of ionizing radiation, and the search for new informat ion, including biochemical and molecular processes. Fluorescence Molecular Tomogra phy (FMT) is one such emerging technique and, like other techniques, has its advan tages and limitations. FMT can reconstruct the distribution of fluorescent molecul es in vivo using near-infrared radiation or visible band light to illuminate the subject. FM T is very safe since non-ionizing radiation is used, and inexpensive due to the compa ratively low cost of the imaging system. This should make it particularly well suite d for small animal studies for research. A broad range of cell activity can be identified by FMT, making it a potentially valuable tool for cancer screening, drug discovery and gene therapy. Since FMT imaging is scattering dominated, reconstr uction of volume images is significantly more computationally intensive than f or CT. For instance, to reconstruct a
vii 323232 image, a flattened matrix with approximate ly 1010, or 10 billion, elements must be dealt with in the inverse problem, while re quiring more than 100 GB of memory. To reduce the error introduced by noisy measurement s, significantly more measurements are needed, leading to a proportionally larger matr ix. The computational complexity of reconstructing FMT images, along with inaccuracies in photon propagation models, has heretofore limited the resolution and accuracy of F MT. To surmount the problems stated above, we decompose the forward problem into a Khatri-Rao product. Inversion of this model is sho wn to lead to a novel reconstruction method that significantly reduces the computational complexity and memory requirements for overdetermined datasets. Compared to the well known SVD approach, this new reconstruction method decreases computatio n time by a factor of up to 25, while simultaneously reducing the memory requirement by u p to three orders of magnitude. Using this method, we have reconstructed images up to 323232. Also outlined is a two step approach which would enable imaging larger volumes. However, it remains a topic for future research. In achieving the above, the author studied the phys ics of FMT, developed an extensive set of original computer programs, performed COMSOL simulations on photon diffusion, and unavoidably, developed visual displays.
1 1. Introduction Fluorescence Molecular Tomography (FMT) is a medica l imaging technology which uses near-infrared radiation or visible band light to il luminate and reconstruct the distribution of fluorescent molecules in deep tissue . FMT s ystems provide functional medical imaging using non-ionizing radiation and relatively inexpensive components. A broad range of cell activity can be identified, making FM T a potentially valuable tool for cancer detection, drug discovery and gene therapy . Ho wever, obstacles must be overcome before this potential can be realized; the resoluti on, imaging volume and accuracy of FMT imaging is currently limited by the computation al complexity of reconstruction and inaccuracies in photon propagation models. It is e xpected that improvements in these areas would facilitate high resolution imaging of d eep structures in vivo making FMT suitable for many applications. To acquire imaging data, the tissue of interest is illuminated with laser light with a wavelength of 650-900 nm; These 'excitation' photon s interact with the tissue as they travel though it, undergoing scattering and absorpt ion interactions ; see Figure 1 for a simple schematic of these photon-matter interaction s. For many biological tissues in this spectral window photon propagation is "scattering d ominated", typically with scattering interactions several orders of magnitude more commo n than absorption interactions . Since absorption is low in this spectral window, a significant proportion of excitation
photons travel several cent high scattering and low absorption, the transport o f photons may be modeled by the diffusion approximation to radioactive transport th eory As the excitation photons diffuse molecules in the tissue ; subject or be produced by excitation photon interact and a 'fluorescent photon of longer These fluorescent photons tissue . Figure 2 photons travel several cent imeters or more into the tissue . Under such conditions of high scattering and low absorption, the transport o f photons may be modeled by the diffusion approximation to radioactive transport th eory . Figure 1. Photon-Matter Interactions As the excitation photons diffuse through the media, they also interact with fluoresc ent ; These fluorescent molecules can either be injected into the or be produced by the subject due to genetic modification excitation photon interact s with a fluorescent molecule, the excitation photon is absorbed photon of longer wavelength is emitted this is shown in These fluorescent photons then undergo a similar diffusionlike transport Figure 2. Simplified Jablonski Energy Diagram Under such conditions of high scattering and low absorption, the transport o f photons may be modeled by the through the media, they also interact with fluoresc ent These fluorescent molecules can either be injected into the subject due to genetic modification . When an the excitation photon is absorbed this is shown in Figure 2. like transport through the
3 Some of these fluorescent photons diffuse to the s urface of the tissue, where they can be detected by direct contact detectors such as optica l fibers  or non-contact detectors such as CCD or CMOS cameras  . By using appr opriate filters, both the intensity of fluorescent photons and excitation photons at the t issue surface can be measured. The normalized born approximation can then be calculate d by dividing the intensity of fluorescent photons by that for excitation photons, leading to better experimental data by canceling out the effects of detector quantum effic iency and source strength, and reducing deviations from the model introduced by he terogeneous optical properties of the tissue . Figure 3. Simplified FMT Imaging System By illuminating the subject with a sufficient numbe r of sources, one at a time, and measuring the resulting fluorescent and excitation light distribution with a sufficient number of detectors, the distribution of fluorochro me inside the subject can be reconstructed. The first step in this reconstructi on is to discretize the tissue volume into a set of three-dimensional volumetric pixels, also k nown as voxels. A model is then used Imaging Plane Imaging Chamber with Animal or Plant Inside Excitation Light Source Fluorescent Tissue
4 to estimate the photon propagation from each source to each voxel to each detector, this is called the forward model . The data from the forward model and imaging measurements can then be written as a system of lin ear equations, where is the flattened weight matrix generated by the fo rward model, is a column vector of normalized experimental measurements, and is a column vector of the unknown fluorochrome concentrations in the voxels [ 2]. Solving this system of linear equations is called the inverse problem, or reconst ruction, and results in an estimate of the fluorochrome concentration, The inverse problem for FMT is significantly mor e computationally intensive than for CT; the photon p ropagation in this wavelength window is scattering dominated, therefore filtered back projection cannot be used. Instead, typical reconstruction methods include Moo re-Penrose pseudoinverse of the matrix by singular value decomposition , and iterative methods such as the algebraic reconstruction technique with randomized projection order (R-ART)  . To reduce the computational complexity of reconstru ction, we decompose the forward problem into a Khatri-Rao product. Inversion of th is model leads to a novel reconstruction method that significantly reduces th e computational complexity and memory requirements for overdetermined datasets. C ompared to the well known singular value decomposition based approach, this new recons truction method decreases computation time by a factor of up to 25, while sim ultaneously reducing the memory requirement by up to three orders of magnitude. An example of an image reconstructed using this novel reconstruction method is shown in Figure 4.
Figure 4. Reconstruction Using SVD 5 Reconstruction Using SVD -KR. Montage on left, 3-D rend er Tumor er ing on right.
6 2. Review of Literature 2.1. Experimental Techniques for Acquiring FMT Imag ing Data The first step in collecting in vivo experimental imaging data typically involves introducing a fluorescent agent into the subject. This is done in one of two ways. The first method is to inject a solution containing one or more different fluorochromes into the subject shortly before imaging . The fluoro chrome gets distributed throughout the subject by bulk transport and diffusion, and in the process the fluorochrome interacts with and binds to its biological target. A variety of f luorochromes are commercially available; they are divided into families with different biolo gical targets, including antibody conjugated molecular probes, nucleic acid probes, f luorescent proteins, reactive probes, and cell function probes. Alternatively, a gene th at expresses a fluorescent protein can be introduced into an organism, allowing for the measu rement of gene expression and regulation . The next step is to collect imaging data using an F MT imaging system. Current FMT imaging systems fall into two major categories, dir ect contact and non-contact. In direct contact imaging systems, optical fibers used both t o illuminate the tissue and measure fluorescence are in direct contact with the surface of the tissue, or paired to the surface by use of a matching fluid with optical properties sim ilar to the tissue . The advantage of direct contact imaging systems is that the shape of the tissue surface does not have to be
7 known or measured since the tissue is either compre ssed into a fixed geometry or placed in a closed container filled with liquid . There are a few disadvantages to this imaging geometry; compressing the tissue into a fixed geome try distorts the natural shape, using a liquid filled chamber precludes the use of live lab oratory animals, and the matching fluid induces additional photon scattering and attenuatio n . In non-contact imaging systems, the tissue is illum inated with a laser and fluorescence is measured with a CCD camera located around the tissu e. The advantages of this geometry are that large datasets from multiple angles can be acquired, and the tissue does not have to be compressed or immersed in fluid. The disadva ntage of this method is that the surface shape of the tissue must be known, and inac curacies in this measurement induce errors in the forward model . Figure 5. Simplified Non-contact FMT Imaging Syste m Before FMT can become a clinically useful tool, sev eral challenges must be surmounted. First, large datasets and extensive computations ar e necessary to obtain a reconstruction
8 with sufficient resolution to be clinically meaning ful. Second, these large datasets need to be collected quickly to minimize errors from mov ement of the subject and physiological changes . 2.2. Model for Photon Diffusion in Scattering Media In order to accurately and quickly reconstruct an F MT image, the forward model of photon propagation has to accurately predict actual photon density in the tissue while being computationally simple. Typically, the forwa rd model of choice is the diffusion approximation to the radiative transfer equation so lved for the type of light source used in the imaging setup. The light source used in FMT is laser light either directly incident on the tissue or coupled to the tissue with an optical fiber in direct contact. In the second case the light source may be modeled as a point sou rce. The light intensity can either be constant for steady-state measurements, frequency m odulated for frequency domain measurements, or a short pulse for time-resolved me asurements. The following sections review literature related to the modeling of photon diffusion in scattering media. In section 2.2.1 we discuss Dirac delta point sources, in section 2.2.2 we discuss boundary conditions, and in section 2.2.3 we discuss sinusoi dally modulated point sources. 2.2.1. Dirac Delta Point Sources Photons propagating through a scattering medium, su ch as tissue, scatter and attenuate as they travel. The radiative transport equation that describes this process is difficult to solve directly, so it is typically approximated by diffusion equations. However, traditional diffusion equations assume isotropic sc attering, while photon scattering is
9 actually anisotropic. Consequently, it is assumed that after numerous anisotropic scattering events the photon distribution will be a pproximately isotropic. To incorporate this assumption, the scattering coefficient, is modified by the average cosine of the angle of a scattering event, resulting in the reduced scattering coefficient, The photon density in the tissue n r induced by a source of photons n may be found by the diffusion approximation to the radiative transfer equation  n n n n n (1) where, for the 3-D case, n (2) and for the 2-D case, n (3) In (1), v is the speed of light in the tissue and has a typi cal value of r which corresponds to an index of refraction of n = 1.4, and D is the diffusion coefficient given by  (4) where is the absorption coefficient, is the scattering coefficient of the media, and is the average cosine of the angle of a scattering event . Values for these coefficients in selected components can be found in Table 1.
10 Table 1. Optical Properties of Tissue Tissue Source mm-1 mm-1 mm-1 mm Bone Pig Skull  0.04 2.625 35 0.925 0.125 Muscle Chicken  0.017 0.33 0.41 0.20 0.961 Skin Albino Murine Dermis  0.28 6.2 23.9 0.74 0. 051 Lung Human  0.81 8.1 32.4 0.75 0.037 Prostate Tumor Rat  0.049 0.81 27.0 0.97 0.388 Blood Human  0.13 0.611 124.6 0.995 0.450 The diffusion equation (1) is then solved for n at points far from sources and boundaries and for media where For wavelengths in the = 650 to 900 nm range in soft tissue, this second condition is g enerally true . The solution to (1) plotted for the 2-D case is shown in Figure 6. Figure 6. Solution to (1) Using COMSOL (mm) Collimated Incident Light 1.e 1 1.e 2 1.e 3 1.e 4 1.e 5 1.e 6 1.e 7 n
In an infinite medium, and for an isotropic point s ource solution to (1) becomes [ 4 n r where has units of from a source at As an example, an isotropic source defined to be at the location illuminates a sample of and rn would the source would b simplify to n This can be plotted as a function of time, resultin g in Figure 11 In an infinite medium, and for an isotropic point s ource 4 ] r n r r has units of and is the distance of point As an example, an isotropic source defined to be at the location illuminates a sample of muscle, with optical properties found in !". At a location in the muscle, n would the source would b e # $ n % and n n &n This can be plotted as a function of time, resultin g in Figure 7. Figure 7. Photon Density vs. Time Given by (6) the (5) is the distance of point As an example, an isotropic source defined to be at the location muscle, with optical properties found in Table 1, the distance and (5) would (6)
12 In the case of fluorescence molecular tomography, t he light source is anisotropic, so modifications must be made to (5). To model this s ituation, it is assumed that all source photons are initially scattered at a depth n. As shown in Figure 8, to specify a boundary condition of n r at the tissue surface, a negative point source is added at . Figure 8. Model of Incident Collimated Light as Poi nt Sources The photon density from an anisotropic source can t hen be written as  n n n (7) 2.2.2. Photon Diffusion Boundary Conditions For the diffusion approximation to the radiative tr ansfer equation, given by (1), absorbing boundaries can be modeled with a Dirichlet conditio n  Legend: Positive Point Source Zero U Boundary Negative Point Source Tissue Collimated Incident Light
13 (8) where In contact imaging systems, boundaries between t he imaging chamber and the sources or detectors are tissue-gla ss-air boundaries, and may be modeled with a modified Robin condition. In non-contact im aging systems, boundaries are tissueair interfaces, which also may be modeled with a mo dified Robin condition   n (9) where !" is the local normal vector at the boundary, and   # (10) # $ % & & ( & & (11) % ! (12) ! r r (13) n ) (14) where is the critical angle at the boundary, !r is the refractive index of the tissue, !r is the refractive index of air, is the relative refraction coefficient across the boundary, % is the power reflection coefficient, and # is a coefficient to describe photon propagation across boundary derived from FresnelÂ’s equations, and where for the 3-D case, n + , , , (15)
14 Alternatively, for the 2-D case, n + , , (16) 2.2.3. Sinusoidal Amplitude Modulated Point Sources In subsection 2.2.1, we studied the photon diffusio n approximation for Dirac delta point sources. That model is helpful for time-gated meas urements, which are typically used for reflection measurements to resolve shallow structur es. On the other hand, deep structures are typically resolved using transmissio n measurements. Transmission measurements are typically frequency domain or stea dy state measurements, acquired using a sinusoidally modulated or steady state ligh t source respectively. Both situations can be modeled in highly scattering media by the re sults of this chapter, (29) and (30). The photon diffusion from a isotropic sinusoidal am plitude modulated point source in a homogeneous isotropic media satisfies the Boltzmann transport equation , n n n r n n (17) n n r n r n (18) where n is the photon density with units r, r n is the photon current density with units and D is the diffusion coefficient given above by (4). The relationship between n and r n will be given later. n and r n can be accurately determined for points far from sources o r boundaries in media where . A sinusoidal isotropic point photon source i s given by 
15 n n / # # 01 2 (19) where # is the amplitude of frequency-dependant source mod ulation, #is the amplitude of dc source intensity where # #, and is the angular frequency of the source modulation. It is assumed that n and r n have the following forms : n n n 0 1 3 (20) r n 4 r n 5 4 r n 5 0 1 6 (21) where 3 and 6 are phase angles. Substituting (19) into (17), an d using the above forms for n and r n results in the steady-state and frequency-dependen t equations for n and r n , n n 4 r n 5 # n (22) 4 r n 5 n n (23) 01 n n 4 r n 5 # n (24) 4 r n 5 7 0 1 1, 8 n n (25) By making the assumption that 1, which is equivalent to making the assumption that the mean free path between scattering events i s much shorter than the wavelength of the sinusoidal modulation, (25) reduces to  4 r n 5 n n (26) where r n has units of and n has units of r By combining (22) and (23) to eliminate 4r n 5, results in 
16 n n n # n (27) Likewise, combining (25) and (26) to eliminate 4r n 5, results in  n n ) 01 n # n (28) It is worth noting that (27) and (28) are the stead y-state and frequency domain equivalents of (1) respectively. (27) and (28) can be solved for an infinite media, resulting in  n # n n 9 # n : n ; 1 < =>7 n) 1 *8? :0n ; 1 < @!7 n) 1 *8 0 1 (29) This equation is the frequency domain equivalent of (5). For the case where (29) simplifies to  n # n # n n 9 1 0n 9 1 0 1 (30) As an example, an isotropic steady-state source, wi th #, #, 1 =, and defined to be at the location n illuminates a sample of muscle, with optical properties found in Table 1 an d r At a
location in the muscle, # $ n % and '( Plotting ) ) as a function of time results i Figure 9. Magnitude of Photon Alternatively, when + n n r As an example, an isotropic steady location illuminates a sample of muscle, with optical prop erties found in Table 1 and r distance to the source would be simplify to n n 17 n, the dis tance to the source would be and (30) would simplify to $ '' n& r as a function of time results i n Figure 9. Magnitude of Photon Density vs. Time for Example (29) simplifies to the steady state solution, / 0 1 As an example, an isotropic steady -state source, with +n and defined to be at th illuminates a sample of muscle, with optical prop erties found in n !". At a location in the muscle, the source would be 2 # $ n % and n 3 (& % 4 % 5 (& 6 (7 tance to the source would be 8 (31) Example (31) (32) and defined to be at th e illuminates a sample of muscle, with optical prop erties found in n, the and (32) would 6 (7 (33)
18 2.3. Normalized Born Field The forward problem for FMT involves not only photo n propagation through a media, but also photon interaction with fluorochromes. To model this combined forward problem a Born-field approximation can be used. No ting from (30) that the photon density attenuates with the form 3n n, and can be represented in the frequency domain in the following simplified form : n n 1 A n 0.n n (34) where n n 1 has units of r and is the photon density at position n due to a source at n with modulation angular frequency 1, A n is the source gain factor, 01 n is the scalar propagation constant, n B n n B is the distance from the source, and n Building upon this equation, the detected intensity of excitation wave length light at a detector located at ndue to a source located at nis given by  r n n C A n A n D n n E (35) where C is the quantum efficiency of the detector for at w avelength A n is the detector gain factor, is the wave propagation constant for excitation wa velength light, and is the wavelength of excitation light and has a ty pical value of 672 nm . The intensity of fluorescence wavelength light detected at a detector located at n due to the fluorochrome distribution n illuminated by a source located at n, is given by 
19 n n F A n A C A n D n n E n 01G n D n n E n (36) where n is the distribution of fluorophore in the media wi th units n, G n is the fluorescence time constant with units of seconds, is the wave propagation vector for the fluorescent wavelength, A is the attenuation of the filter used to select fo r detection only fluorescent light, and C is the quantum efficiency of the detector for fluo rescent wavelength light. is the wavelength of fluorescent light and has a t ypical value of 710 nm . Since A n and A n, the source and detector gain, are position depend ent, finding a solution to (35) and (36) would require m easurement of these values for each pair of sources and detectors. To reduce the expe rimental measurements required, (35) and (36) can be combined to eliminate these terms r esulting in the normalized born field equation , n n A n n r n n C C (37) By assuming that C C due to the similarity of these wavelengths in FMT imaging systems, and combining (37) with (35) and (36), yie lds  n n n n FDn n .E n 01G n D n n E n (38)
20 This equation is an approximation of the normalized born field. In the section that follows, this equation will be used to generate th e weight matrix in the forward problem. 2.4. Reconstruction Methods The inverse problem for FMT requires solving the fo llowing equation : (39) where is the measurement vector given by (38), is the weight matrix relating each measurement to each voxel of fluorochrome, and is the fluorochrome concentration column vector. In (39), has units r has units r and has units n; while this might seem unbalanced at first glance, the matrix multiplication is performing a summation over three spatial dimension s, leading to the correct balance of units. For the steady-state case, where 1 is calculated using the right side of (38), as shown in the following equation : n n n D n n E D n n E n n (40) where W is typically a 2-D matrix, and for each element each i is a unique source detector pair, and each j is a unique voxel in the volume to be imaged.
Figure 10 Calculation of Weight Matrix for Simplified 2 An example of how one element system is shown in Figure the location of the annotated The system of equations concentration vector, when (40). Ther e are multiple methods for used an iterative method called the method of proje ctions reconstruction technique with randomized projection order (R 21 Calculation of Weight Matrix for Simplified 2 D FMT System one element of this matrix would be calculated for a simple Figure 10; it is important to note that the tip of each vector represents the location of the annotated shown in (39) can be solved for 9: the estimated when is determined experimentally and ; is calculated using e are multiple methods for solving this system of equations used an iterative method called the method of proje ctions also known as the algebraic reconstruction technique with randomized projection order (R -ART). In the method of D FMT System would be calculated for a simple 2-D the tip of each vector represents the estimated fluorochrome is calculated using solving this system of equations Graves (2003) also known as the algebraic In the method of
22 projections, each row vector of and its corresponding value of defines an affine hyperplane. The guess for the value of r is updated each iteration by projecting the previous guess, rn onto an affine hyperplane using the following for mula : r r n r n H r r H r H r H r (41) where r is the fluorochrome vector after the projection, rn is the fluorochrome vector before the projection, H r is the ith row of the matrix, and r is the ith element in the measurement vector. Figure 11 shows an example of the method of projections for solving a system of two equations a nd two unknowns. An initial guess is projected onto a hyperplane, and the resulting poin t is projected onto the other hyperplane; this is repeated until convergence to w ithin a specified tolerance is achieved. Figure 11. Method of Projections for Two Equations and Two Unknowns H H H H Initial Guess Projections Solution
23 Alternatively the Moore-Penrose pseudoinverse of th e matrix W can be found using singular value decomposition to obtain a direct sol ution for . This method requires finding the singular value decomposition of the mat rix W , I (42) where and I are orthogonal basis matrices and is a diagonal matrix that contains the singular values. Using the above singular value de composition, a solution for (39) can be obtained: I n (43) where is the Moore-Penrose pseudoinverse of To enhance the accuracy of the reconstruction from noisy measurements, singular va lues below a certain threshold are typically discarded . For a reconstruction of an N x N x N voxel volume, the minimum size of the matrix W is N3 by N3. As an example, for N = 32 the minimum size of the matrix W would be 32,768 by 32,768 for N = 64 this increases to 262,144 by 262,144 To obtain better reconstructions with noisy measurements more imagin g data than this minimum value is typically collected, making the system overdetermin ed. Solving systems of equations with such a large number of equations and unknowns requires a considerable number of computations and large amounts of memory.
24 3. Challenges in FMT Reconstruction Past reconstructions of FMT images have been limite d by memory and computational constraints to a relevantly small number of sourcedetector pairs, and consequently are restricted to relatively low spatial sampling. To o vercome these limitations, the formulation of the forward problem was decomposed i nto a diagonal matrix multiplied by a Khatri-Rao product. For overdetermined cases, th is decomposition significantly reduces the computational and memory complexity of reconstruction; this allows for larger imaging datasets with more source-detector p airs and higher spatial sampling. Higher spatial sampling has been shown to allow hig her reconstruction resolution, improve the signal to noise ratio of measurements, and improve reconstruction image quality . Additionally, for overdetermined dat asets those with more source-detector pairs than reconstruction voxels this method sign ificantly reduces the computational complexity of reconstruction compared to method of projections reconstructions and the memory complexity of reconstruction when compared t o SVD based Moore-Penrose pseudoinverse reconstructions. The performance of t his algorithm will be compared against these reconstruction methods. The reconstruction of FMT images is an ill-posed pr oblem, with a poorly conditioned weight matrix. Consequently, small errors in the f orward model can create significantly larger errors in the reconstructed image. Improvin g the prediction accuracy of photon
25 propagation in the forward model would therefore al low for significant increases in the accuracy and image quality of FMT reconstructions.
26 4. FMT Reconstruction Via Khatri-Rao Decomposition Large measurement datasets are necessary to obtain reconstruction with sufficient resolution to be clinically meaningful. However, t he computational complexity of reconstructing large datasets can be prohibitive. To reduce the computational complexity of reconstruction, the forward problem weight matri x, (40), was decomposed into a product of smaller matrices. This reformulation of the forward problem significantly reduces the computational complexity of reconstruct ing overdetermined datasets. 4.1. Theory 4.1.1. Weight Matrix in Tensor and Array Forms The weight matrix relates the experimental measurem ents of each source and detector pair in the imaging system to the concentration of fluorochrome in each voxel. The weight matrix is generated by a forward model. The forward model used here is the normalized Born equation, given by (40); however th e following method can be used for any forward model in imaging systems where every so urce contributes to the detected signal at every detector. Each element of the weig ht matrix specifies a unique sourcevoxel-detector combination, and is calculated as sh own in Figure 10 and (40). Although the matrix W is shown as a three dimensional tensor in (40), fo r reconstruction purposes the matrix is typically flattened to a two dimensional array of size m by n where
27 m is the number of source-detector pairs and n is the number of voxels being reconstructed, n n n J K K K K L n n (44) Flattening the matrix can be represented graphically: Figure 12. Flattening the Weight Matrix from a Tens or to an Array 4.1.2. Decomposition of Weight Matrix to Extract No rmalizing Term When the matrix W is flattened, the normalizing term Dn n .E becomes a constant for each row, n Consequently, the matrix W can be expressed as a constant times the matrix product of a diagonal matrix to the left of the non-normalized Detectors Voxels Sources W1 . WN W1 WN Voxels SourceDetector Pairs 3rd Order Tensor Array .
28 matrix W (45) where the diagonal entries of are n n n (46) The non-diagonal entries are equal to zero. This d ecomposition can be represented graphically as: Figure 13. Decomposition of Weight Matrix to Remove Normalizing Term In general, if this method is used for a different forward model which does not have a normalization term, would be equal to the identity matrix, M, and this step can be omitted entirely. . . . N . . . N N !" Full matrix (arrows represent rows n ) Full matrix (arrows represent rows n ) Diagonal matrix (arrow represents nonzero elements)
29 4.1.3. Decomposition of Weight Matrix to Khatri-Rao Product the non-normalized matrix W can be decomposed to be expressed as a Kroneckerlike product of two smaller matrices. Before discu ssing this decomposition in detail, some basic information on the Kronecker product and Khatri-Rao product is provided below. Given two matrices A and B each of size 2 by 2, the Kronecker product of A and B is equal to # O 7 P O P O POPO 8 Q P R P R P R P R P R P R P R P R PRPR P R P R PRPR P R P R S (47) where is the Kronecker operator. More generally, if A and B are of size ma by na and mb by nb respectively, the Kronecker product of A and B results in a matrix of size ma mb by na nb formed by multiplying each element in A by every element in B . The Khatri-Rao (KR) product, also called the column wise Kronecker product, is a matrix operation that is related to the Kronecker product. Given two matrices C and D of size mc by n and md by n the KR product is defined as  T = = # = (48) where is the KR operator. The resulting matrix is of si ze mc md by n and is a subset of the columns of T . For matrices C and D each of size 2 by 2, t heir KR product is equal to
30 T Q = = = = =,=, = = S (49) Note that the matrices C and D must have the same number of columns. can be decomposed to be expressed as a KR product of the matrices and (50) where n n D n n E (51) n n D n n E (52) (50) can be represented graphically as shown in Fig ure 14, for a reconstruction of an n voxel image from an imaging system with ms sources and md detectors. Figure 14. KR Product Decomposition of Wsd . N Full matrix Size: msmd by n (arrows represent columns n ) . N N . Full matrix Size: ms by n (arrows represent columns n ) Full matrix Size: md by n (arrows represent columns n )
31 By combining (45) and (50), W can be rewritten as (53) 4.1.4. Pseudoinverse of Khatri-Rao Product From (53), the Moore-Penrose pseudoinverse of can be written as n (54) where is given by (46), is given by (51), is given by (52). The inverse of since it is a diagonal matrix, is the element-wis e inverse of each non-zero element in the matrix. A simple example of the inverse of a diagonal matrix such as is U r V r# W K K K L n U r r V (55) The following identity is then applied to the pseud oinverse of a KR product shown in (54); given matrices A and B , # O # O # # $ O O (56) where $ is the Hadamard operator, also known as the elemen t-wise multiplication operator, defined as X $ Y $ ( $ 7 Z Z Z Z 8 7 Z Z Z Z 8 (57) (56) can be rewritten as  # O / # # $ O O 2 # O (58) (58) can be substituted into (54) for # and O resulting in
32 $ [ $ \ ( n (59) This is a new result to the best of my knowledge. As will be discussed shortly, this result can significantly reduce the reconstruction computa tional complexity and memory requirements for reconstruction when the matrix W is overdetermined. 4.2. Reconstruction Algorithms By combining (42) and (59), the following equation for the reconstructed fluorochrome concentration, is obtained: [ $ \ n (60) Recall that is given by (46), is given by (51), is given by (52), and $ is the Hadamard operator given by (59). To simplify the n otation in the following discussion, (59) will be rewritten as # R (61) where # [ $ \ (62) R n (63) For an imaging system with sources and detectors, for a total of source-detector pairs, and with voxels in the forward model, both the matrix and are of size by For overdetermined systems with % these matrices can become prohibitively large, requiring significant amounts of memory to store. Consequently, a reconstruction method that does not require these matrices to be
33 stored in memory could significantly reduce the mem ory complexity of the problem. Two different reconstruction algorithms that satisf y this criteria were investigated. Both methods involve reconstruction by taking the singul ar value decomposition of a matrix generated from a KR product, this method will be ca lled SVD-KR reconstruction for short. 4.2.1. Method 1: Row-Wise SVD-KR Reconstruction In this reconstruction method, one row at a time of the matrix is calculated and multiplied by the pre-computed vector $# n and stored as the ith element of the column vector R The result after iterating through all n rows of the matrix is an n by 1 column vector, R # is then calculated and multiplied by R yielding the reconstructed fluorochrome concentration, A flowchart of this algorithm is shown in Figure 15, with acquired data shown in dark gray, c alculated data in textured light gray, and methods in white.
34 Figure 15. Row-Wise SVD-KR Reconstruction Algorith m 4.2.2. Method 2: Column-Wise SVD-KR Reconstruction In this reconstruction method, one column at a time of the matrix is calculated and multiplied by the corresponding elem ent of the vector $# n Iterating through all m columns of the matrix, the resulting vectors from each calculation Calculate and store: $ # n Calculate row i of: Last row? Multiply and store as R r No Yes Imaging Measurements Normalized Born Field Generate Forward Model: and Calculate and store # matrix Multiply and store as Reconstructed Fluorochrome Concentration Limiting step memory i = i + 1
35 are summed resulting in the n by 1 column vector, R # is then calculated and multiplied by R yielding the reconstructed fluorochrome concentra tion, A flowchart of this algorithm is shown in Figure 16, with acqui red data shown in dark gray, calculated data in textured light gray, and methods in white. Figure 16. Column-Wise SVD-KR Reconstruction Algor ithm Calculate i th element of: $ # n Calculate column i of: Last column? Multiply then sum result with R No Yes Imaging Measurements Normalized Born Field Generate Forward Model: and Calculate and store # matrix Multiply and store as Reconstructed Fluorochrome Concentration Limiting step memory i = i + 1
36 4.2.3. Minimizing Reconstruction Errors from Noisy Measurements The reconstruction of FMT images is an ill-posed pr oblem with a poorly conditioned weight matrix. Consequently, small errors in meas urements or models can cause significant errors in the reconstruction. The lea st squares solution to the linear equation # R is given by (60); the reconstruction error from t his least squares solution is given by  ] ] ] ] & B # B B # n B B # B B # B (64) ] ] ] ] & B # B B # n B B R B B R B (65) Combining (64) and (65) results in ] ] ] ] & B # B B # n B ] R ] ] R ] B # B B # B (66) where is the reconstruction error given by R is the measurement error, # is the forward model error, and B # BB #nB is the condition number of #. To better understand this equation and its implications, we n eed to first state the relationship between the singular value decomposition of # and the condition number of #. As shown in (42), A can be expressed as # % % I % (67) where % and I% are orthonormal matrices and S is a diagonal matri x of singular values arranged in a descending order. The norm of #, B # B, can then be written as B # B B % B B % B ] I % ] (68) Since the norm of an orthonormal matrix is 1 and th e norm of a diagonal matrix is equal
37 to the largest entry on the diagonal, B # B is therefore equal to the largest singular value i n %, ^&'. Likewise, B #nB is equal to the inverse of the smallest singular v alue, ^&r n. Therefore the condition number of # is equal to B # B B # n B ^ &' ^ &r (69) Inspecting (66), we can see that by reducing the co ndition number of #, the relative error in the reconstruction, ] ]] ] _, can be reduced as well. This can be accomplished by replacing # with a reduced rank approximation; assuming that # has n singular values, for n a rankr approximation to #, # `, is equal to # ` I (70) where is the first r columns of %, I is the first r columns of I%, and is diagonal matrix with the first r singular values of %. By using only the r largest singular values ^&r is increased, which decreases the condition number of # `, leading to reduced reconstruction error. However, this approximation also introduces error by decreasing the accuracy of the forward model; the approximati on error increases B # BB # B in (66), leading to increased reconstruction error. Inspect ing (66), we can see that for values of r that are not significantly less than n, ]R ]]R ] % B # BB # B and so the effect of increasing B # BB # B will be less than the associated decrease of B # BB #nB. However, if r is decreased further, B # BB # B can become a significant term when ]R ]]R ] B # BB # B introducing a significant source of reconstructio n error. Consequently, in selecting a value for r a balance must be struck, since lowering the value of r reduces the noise in the reconstruction but if lowered too far can introduce significant errors.
38 Previous research has addressed the issue of select ing r by analyzing the singular value decomposition of the weight matrix  . As sh own in (42), the singular value decomposition of the weight matrix is I, which combined with (39) can be rewritten as  I (71) can be decomposed into two components, r( and )r by analyzing multiple measurement sets . The vectors r( and )r are calculated, filtered, and plotted; the value of r is selected as the index at which r( crosses and becomes smaller than )r . At this index, addition of additional singul ar values to the approximation will contribute more noise to the rec onstruction than signal, increasing the reconstruction error. This method of selecting r can be extended for use in the reconstruction meth ods presented here. (67) and (60) can be combined and r ewritten to obtain % R % I % (72) This equation can be combined with (63) to yield 4 % n 5 % I % (73) Like before, can be decomposed into two components, r( and )r resulting in two equations, 4 % n 5 r( % I % r( (74)
r < ; = The two resulting vectors from theses equations are filtered with a median filter and then plotted against each other. The highest index for which gives the value of r to u se for the rank example for a reconstruction of 512 voxels using synthetic meas urements from 256 sources and 1024 detectors with a 40 dB SNR. the value of r was determined to be 262. To confirm this result, were performed for > fluorochrome concentration. This analysis confirme d that root mean squared error ( in Figure 18 It is worth noting that the value of measurement noise and the number These plots were created using a synthetically gene rated FMT imaging data. Figure 17 Plot of Signal and Noise Contributions 39 ; ; ? @ 9 The two resulting vectors from theses equations are filtered with a median filter and then plotted against each other. The highest index for which @ 9 se for the rank r approximation of the matrix + a reconstruction of 512 voxels using synthetic meas urements from 256 sources and 1024 detectors with a 40 dB SNR. From the data used to generate was determined to be 262. To confirm this result, SVDKR > >% and the RMS error measured against the known fluorochrome concentration. This analysis confirme d that & has the minimum root mean squared error ( RMSE), with a value of 0.6370 A plot of the results It is worth noting that the value of r depends on many variables, such as measurement noise and the number and placement of voxels, sources and detectors. These plots were created using a synthetically gene rated FMT imaging data. Plot of Signal and Noise Contributions vs. Singular Value Index (75) The two resulting vectors from theses equations are filtered with a median filter and then A@ 9 Figure 17 is an a reconstruction of 512 voxels using synthetic meas urements from 256 From the data used to generate Figure 17, KR reconstructions and the RMS error measured against the known has the minimum A plot of the results is shown on many variables, such as and placement of voxels, sources and detectors. These plots were created using a synthetically gene rated FMT imaging data. Singular Value Index
40 Figure 18. Reconstruction Error vs. Rank 4.3. Results for the New FMT Reconstruction To solve (60) for the Moore-Penrose pseudoinverse of the matrix # must be calculated. The matrix # is of size n by n for any matrix W of size m by n where n is the number of voxels in the forward model and m is the number of source-detector pairs. Since mat rix A is of size n by n the computational and memory complexity of calcul ating its pseudoinverse is dependent solely on the number of voxels in the forward problem and not on the number of source detector pairs; since t his step is the limiting step in terms of memory usage, the reconstruction algorithm is able to solve for systems with a large number of sources and detectors. Previous imaging systems have typically used an und erdetermined system of equations for W because of the computational limitations of existi ng reconstruction techniques.
41 Compared to these existing techniques, the SVD-KR r econstruction method significantly reduces the computational complexity and memory req uirements of reconstruction when the matrix W is overdetermined. To test the performance of this method, synthetic v olume image data was generated using a diffuse non-ellipsoidal phantom designed to emula te structures that may be encountered in vivo A slice of this phantom is shown in Figure 19. Figure 19. Diffuse Non-ellipsoidal Phantom Used in Reconstruction It was then reconstructed using three different rec onstruction techniques: row-wise SVDKR, column-wise SVD-KR, and singular value decompos ition. For simplicity, these tests used an equal number of sources and detectors for e ach dataset. Results of this test are shown in Figure 20, Figure 21, Figure 22, Table 2, Table 3, Table 4, and Table 5.
42 Table 2. SVD-KR and SVD Reconstruction Time Voxels SourceDetector Pairs Row-Wise SVD-KR Reconstruction Time (sec) Column-Wise SVDKR Reconstruction Time SVD Reconstruction Time 1,000 10,000 4.7 4.7 11.4 1,000 50,625 6.4 7.2 31.4 1,000 160,000 14.2 13.8 86.1 1,000 390,625 29.7 28.3 709.6 1,000 810,000 62.2 52.9 Out of memory 1,000 2,560,000 211.0 165.2 Out of memory 1,000 6,250,000 454.4 378.3 Out of memory 1,000 12,960,000 1,022.6 803.2 Out of memory Figure 20. SVD-KR and SVD Reconstruction Time nnnnrnr r r r
43 Table 3. SVD-KR and SVD Calculated Reconstruction M emory Usage Voxels SourceDetector Pairs Row-Wise SVD-KR Reconstruction Memory Usage (MB) Column-Wise SVDKR Reconstruction Memory Usage (MB) SVD Reconstruction Memory Usage (MB) 1,000 10,000 37 37 305 1,000 50,625 44 44 1,545 1,000 160,000 55 55 4,883 1,000 390,625 69 69 11,921 1,000 810,000 85 85 24,719 1,000 2,560,000 128 128 78,125 1,000 6,250,000 183 183 190,735 1,000 12,960,000 250 250 395,508 Figure 21. SVD-KR and SVD Calculated Reconstruction Memory Usage nrnr r r
44 Table 4. SVD-KR and SVD Reconstruction Errors Voxels SourceDetector Pairs Row-Wise SVD-KR Reconstruction RMSE Column-Wise SVDKR Reconstruction RMSE SVD Reconstruction RMSE 1,000 10,000 0.472 0.472 0.423 1,000 50,625 0.402 0.402 0.339 1,000 160,000 0.332 0.332 0.293 1,000 390,625 0.301 0.301 0.283 1,000 810,000 0.272 0.272 Out of memory 1,000 2,560,000 0.243 0.243 Out of memory 1,000 6,250,000 0.234 0.234 Out of memory 1,000 12,960,000 0.220 0.220 Out of memory Figure 22. Reconstruction Time vs. RMSE From these results, we can see that the SVD-KR reco nstructions were significantly faster than SVD reconstruction, with speed increases of up to 25X. The column-wise SVD-KR reconstruction method performed slightly faster tha n the row-wise SVD-KR. SVD-KR was able to reconstruct four large datasets that co uld not be reconstructed by the SVD reconstruction algorithm; the largest dataset would have needed approximately 386 GB nrnrn r r r
45 of memory for reconstruction with SVD. Error was s lightly higher for the SVD-KR reconstructions compared to the SVD reconstructions when solving for the same measurement dataset. However, Figure 22 shows that SVD-KR can reconstruct to the same accuracy as SVD in less time. Alternatively, SVD-KR can reconstruct more accurately than SVD in the same time. Since the previous datasets used an equal number of sources and detectors, new datasets were generated to determine the effect of the sourc e to detector ratio on reconstruction. The results are shown in Table 5. Table 5. SVD-KR Reconstruction Time and Error Voxels SourceDetector Pairs Sources Detectors Row-Wise SVDKR Reconstruction Column-Wise SVDKR Reconstruction Time Error Time Error 1000 640000 6400 100 45.4 0.541 37.1 0.541 1000 640000 1600 400 47.7 0.441 41.8 0.441 1000 640000 400 1600 48.0 0.493 44.0 0.493 1000 640000 100 6400 53.0 0.523 38.2 0.523 These results show that reconstruction error is min imized when the ratio between sources and detectors is near one. Many current FMT setups use more detectors than sources; in order to make best use of SVD-KR reconstruction, th e number of sources in new systems needs to be increased. Additionally, since measurem ent sets from large detector arrays can be reconstructed using this method, many additi onal CCD cameras could be arrayed around the sample to collect additional measurement s. The improved detector spatial diversity and spatial sampling in such a system wou ld reduce error in the reconstruction.
As an example of the image reconstructions obtained using SVD voxel image was reconstructed using imaging data fr om 9216 sources and with a signal to noise ratio of 40 dB. The results are shown in Figure 23 Reconstruction Using SVD 46 As an example of the image reconstructions obtained using SVD KR, a 24 voxel image was reconstructed using imaging data fr om 9216 sources and with a signal to noise ratio of 40 dB. The results are shown in Figure 23 . Reconstruction Using SVD -KR. Montage on left, 3D rendering on right. Tumor KR, a 24 x 24 x 24 voxel image was reconstructed using imaging data fr om 9216 sources and 9216 detectors D rendering on right.
47 5. Modeling of Photon Diffusion The reconstruction of FMT images is an ill-posed pr oblem, with a poorly conditioned weight matrix; consequently, small errors in the fo rward model can create significantly larger errors in the reconstructed image. Improvin g the prediction accuracy of the photon propagation forward model would therefore allow for significant increases in the accuracy and image quality of FMT reconstructions. In this chapter, we will study photon diffusion in media with homogeneous and hete rogeneous optical properties. A finite element model is developed in COMSOL, and is validated using existing models. COMSOL is a finite element based multi-physics simu lation program. Finite element models discretize the volume into a set of nodes an d elements. This allows a continuous function to be approximated. Additionally, the fin ite element mesh is able to adapt the size of elements where necessary in order to better approximate the function near boundaries and small features such as point sources This finite element model is then used to study the effects of heterogeneous optical properties on the accuracy of the normalized Born field in models that assume homogen eous optical properties. 5.1. Photon Propagation Model for Homogeneous Media Photon diffusion in homogeneous media was modeled u sing two different approaches: finite element modeling in COMSOL, and method of so urces using Green's function.
48 5.1.1. Finite Element Model The diffusion approximation to the radiative transf er equation, given by (1), was modeled for a 3-D imaging system using COMSOL 3.4a. The ph oton density inside of the imaging chamber, n was solved for using stationary analysis of the p artial differential equation given by (1). A schematic of the imaging chamber setup is shown in Figure 24. Figure 24. Imaging Chamber Schematic for COMSOL Sim ulation The imaging chamber was modeled as a cube. Two different boundary conditions were applied to the chamber; Ab sorbing boundaries were modeled with the Dirichlet condition from (8), while tissue -glass-air interfaces located at the source plane and detector plane were modeled with t he modified Robin condition from (9) though (15). The glass was assumed to be treat ed with an anti-reflective coating, giving a value of a r for the modified Robin conditions. The incident c ollimated 50 mm 50 mm 50 mm Tissue-Glass-Air Boundary (Shaded) Absorbing Boundary Point Source at Depth 1/t
49 laser light was modeled as an isotropic point sourc e at depth n . Optical properties inside the imaging chamber were assumed homogeneous with r, r(( and r(( 5.1.2. Method of Sources with Green's Function To validate the precision of the above COMSOL finit e element model, a model was constructed based on the method of sources  and Green's function  given by (34). Like the COMSOL model above, the incident collimate d light is modeled as an isotropic point source at a depth of . Zero photon density boundary conditions were assumed at absorbing boundaries and at an extrapolated boun dary offset from glass-tissue-air interfaces , % % (76) where is the diffusion coefficient, and % is the Fresnel reflection at the boundary. For anti-reflection coating glass, it was assumed t hat %. These boundary conditions were met by adding additional positive a nd negative isotropic point sources . Figure 25 shows how these boundary condition s were enforced the at the source.
50 Figure 25. Source Boundary Conditions In general, to enforce a boundary condition a copy of the point source(s) with opposite sign needs to be mirrored across the boundary. Mu ltiple boundary conditions can be enforced by mirroring the point sources across each boundary. Figure 26 shows how this would be accomplished for the 4 boundary conditions from the 4 sides of a square imaging system. * Glass Boundary: Absorbing Boundary: Collimated Incident Light Legend: Positive Point Source Zero U Boundary Negative Point Source Tissue
51 Figure 26. 2-D Method of Sources to Enforce Boundar y Conditions Although hard to visualize, the method of sources p resented in Figure 26 can be extended to 3-D to enforce the 6 boundary conditions from th e 6 sides of the cubic imaging chamber. For the 3-D case, Figure 26 can be though t of as a perpendicular slice though the imaging chamber that contains the source. Whil e 4 of the 6 sides of the imaging Legend: Positive Point Source Zero U Boundary Negative Point Source Tissue
52 chamber intersect the page, 2 sides of the imaging chamber would lie parallel to the page, above and below; to enforce these 2 boundary condit ions, all of the sources in Figure 26 would need to be mirrored over each of these bounda ries into and out of the page respectively. 5.1.3. Model Comparison A 3-D method of sources model with 36 sources was c onstructed in MATLAB. The photon density in the imaging chamber induced by ea ch source is given by (34). To compensate for differences in A between the two models, a MATLAB symbolic equation was built with the form: Z b ) P 0. ) (77) where 01 n, is the distance from the source, and P is a symbolic variable to allow for the determinat ion of an optimum value of A to compensate for differences in source intensity betw een the two models. The resulting equation was fit to the COMSOL data using MATLAB su rface fitting tool, this is shown in Figure 27. Referring to Figure 25(a), for a 1-D system with 2 sources (77) becomes Z P D 0. E D 0. E (78) Recall that is the speed of light in the media, is given by (4), n, and 01
53 Figure 27. Surface Fit of Finite Element Model by M ethod of Sources Model Figure 28. Residuals of Surface Fit ln(Photon Density) X Y
54 A value of P r was found to minimize error between the two models resulting in a RMSE of 0.004991. Further than 1 mm from the sou rce, the models had a maximum error of 8%. Much of this error was likely the res ult of interpolating values from a finite element mesh used in COMSOL into a regular grid. 5.2. Finite Element Model for Heterogeneous Media Many previous FMT imaging systems have used forward models which do not explicitly take into account the heterogeneous optical propert ies of the various tissues being imaged    . Instead, it is assumed that the optical properties in the tissue are homogenous and equal to the average optical propert ies of the tissues; errors introduced by heterogeneous optical properties are partially c anceled out by dividing fluorescent measurements by intrinsic measurements to obtain th e normalized Born field. This provides an accurate approximation when small, simp le heterogeneities are present, such as those found in deliberately constructed phantoms However, to resolve deep structures with high resolution in vivo despite the complex anatomical structures and the diverse optical properties of tissue, a forward model that more accurately takes into account heterogeneous optical properties is necessary. One promising forward model involves using a hybrid CT-FMT system; In such a system, anatomical data would be collected by the C T system concurrently with fluorescent data by the FMT system. The CT dataset could be used to determine the shape of anatomical structures, classify their comp osition, and then lookup their known optical properties. The shape and composition of th e anatomical structures extracted from
55 the CT data could then be used to accurately model photon propagation using finite element methods. The photon density data generated by this method could then be used in a reconstruction method such as the SVD-KR recon struction method presented in the previous chapter. 5.2.1. CT Image Segmentation From the images produced by the CT scanner, anatomi cal structures can be extracted by one of many image segmentation techniques . Fo r the purposes of this paper, kmeans clustering based image segmentation was used, however other image segmentation techniques would be suitable as well. Image segmen tation was performed on the Digimouse CT dataset  . The right hind leg of the mouse was cropped from th e dataset and segmented using kmeans clustering based on the voxel intensity. Bas ed on the major structures present in the leg, the voxel intensities were grouped into 5 clusters: bone, skin/fat, fast twitch muscle, slow twitch muscle, and the surrounding air Initial values were set manually using typical values for these structures. The vox el intensity values for the 5 clusters were obtained using a k-means clustering algorithm; these values were used to segment the dataset into the tissue types represented by th e clusters based on the intensity of each voxel. Each tissue was then converted into a 3-D m esh. The resulting 3-D structures are shown in Figure 30.
56 Figure 29. Segmented Digimouse Skeleton Showing Reg ion of Interest (ROI) Bone Only With Muscle Added With Skin and Fat Added Figure 30. Segmented Mouse Leg 5.2.2. Finite Element Modeling of Photon Propagati on To analyze photon propagation through heterogeneous media, the 3-D meshes shown in Figure 30 were imported into COMSOL. The photon de nsity inside of the imaging ROI
57 chamber, n was solved using stationary analysis of the follo wing partial differential equation: n n n n n n n n (79) This is a more general form of (1) which is valid f or heterogeneous optical properties. A schematic of the imaging chamber setup is shown in Figure 31. Figure 31. Imaging Chamber Schematic for Heterogene ous COMSOL Simulation A small amount of fluorochrome, modeled as a point source, was placed between the fat and muscle on the lateral face of the right lower l eg. The optical properties for each of the tissues were assigned the following values: 10 mm 7.6 mm 12 mm Tissue-Glass-Air Boundary (Shaded) Absorbing Boundaries Point Source at Depth 1/t 101 by 101 Detector Array Mouse Leg Surrounded by Matching Fluid
58 Table 6. Optical Properties of Selected Tissue Tissue Source mm-1 mm-1 mm-1 mm Bone Pig Skull  0.04 2.625 35 0.925 0.125 Muscle Chicken  0.017 0.33 0.41 0.20 0.961 Skin Murine Dermis (Albino)  0.28 6.2 23.9 0.74 0.051 Matching Fluid Intralipid and Ink Solution  0.03 1.0 0.1 0.90 0.324 By building a similar model with the correct physic al layout of the sources, detectors and imaging chamber, this model can be used to calculat e the , and matrices from (53) for many FMT imaging systems. This method off ers the advantage of explicitly taking into account heterogeneous optical propertie s, leading to increased model accuracy and more accurate reconstructions. 5.3. Precession of Normalized Born Field in Heterog eneous Media The model from chapter 5.2.2 was used to generate t wo sets of synthetic measurement data, differing only in the optical properties of t he mouse leg; one with the heterogeneous optical properties listed in Table 6, the other wit h homogeneous optical properties equal to the matching fluid. Fluorescent and intrinsic m easurements were taken with each system for one source by a 101 by 101 array of dete ctors. The normalized Born field was calculated for each m easurement set, and the error calculated. The root mean squared error (RMSE) was found to be 0.2985. While the accuracy of these heterogeneous and homogeneous mea surements cannot be determined
59 without experimental measurements, the low precessi on in the normalized Born field measurements between two measurement sets that diff er only in the assumption of homogeneity shows that the normalized Born field ca n have errors of at least this magnitude introduced by this assumption. However, t he normalized Born field did significantly reduce the errors introduced by heter ogeneous optical properties compared to non-normalized measurements, which were found to have an RMSE over 50 times larger. To validate the accuracy of the above homogeneous m easurement set, a third measurement set was synthetically generated using h omogeneous optical properties and the method of sources from chapter 5.1.2. The erro r between the normalized Born fields of the two homogeneous measurement sets was found t o be 0.0180.
60 6. Administering and Imaging Multiple Fluorochromes Simultaneously There are a variety of commercially available fluor ochromes that vary in biological target, excitation wavelength and fluorescent wavel ength. Different classes of fluorochromes include antibody conjugate probes, nu cleic acid probes, cell function probes and fluorescent proteins. The ability to i mage multiple distinct probes simultaneously, and independently reconstruct the d istribution of each probe could increase the visibility of low contrast targets, an d allow researchers to better understand relationships between biological processes. Existing FMT reconstruction techniques can be exten ded to reconstruct the distribution of multiple fluorochromes. A mixture of two or more f luorochromes with affinities for different biological targets would be injected into the subject by a single syringe, and be distributed by blood flow and diffusion until inter acting with their respective targets. These fluorochromes would be chosen to emit light o f different wavelength when excited by a single excitation wavelength. Figure 32 shows a synthetic example of intensity vs. wavelength for an imaging system with a source of l ight with peak wavelength 1, and two fluorochromes that fluoresce with peak waveleng ths 2 and 3.
61 Figure 32. Synthetic Intensity vs. Wavelength Measu red at Detector Transmission measurements in a typical single fluor ochrome system are taken with one CCD camera fitted with a band-pass filter. Transmi ssion measurements in a multifluorochrome system would still be taken by one CCD camera, however this camera would have a changeable band-pass filter for each d istinct fluorochrome. Images would be acquired using the same procedures for acquiring single fluorochrome FMT images . For each source, three separate images would be acquired by a CCD camera with a different filter used for each image; One band-pass filter for each peak frequency, 1, 2, and 3. Alternatively, three CCD cameras could be used, each with a different filter, resulting in faster acquisition times but potential ly introducing errors due to their different physical positions. Either of the above methods would result in three sets of measurements: one for the intrinsic illumination at the excitation wavelength, 1, and one for each fluorochrome, 2 and 3. These measurements are collected into a matrix, one column for each wavelength, 1 2 3
62 n 4 n n n 5 (80) This matrix, is of size m by 3 where m is equal to the number of source-detector pairs. Since band-pass filters for each peak wavelength do not exclude all of the light from the source or other fluorochromes, this matrix is a mix ture of the columns of a more fundamental matrix, that describes the true photon intensity at the d etector from the source and each of the fluorochromes. The matrix can be calculated by multiplying by the Moore-Penrose pseudoinverse of the mixing ma trix, R where R is given by % U * * * * V (81) where *+, is the coefficient that corrects for the transmitt ance of r thru the filter for -. It could have appropriately been called however this simpler notation is being used for convenience. If no filter is used for 1 , the matrix R can be simplified to % U * * ** V (82) Multiplying by the Moore-Penrose pseudoinverse of the matrix R results in the unmixed Born field for each fluorescent wavelength, % (83) This expression can be expanded to 4 5 c % % % % % % % % % d (84) where %ris the element of R+ in the ith row and jth column, is the ith column of
63 and is the ith column of The matrix can be normalized to obtain an approximation to the normalized Born field by perfo rming an element-wise division of the first column of into the other columns of n + n n n n (85) The matrix is of size m by 2 where m is equal to the number of source-detector pairs. This matrix can then be used to solve for the conce ntration of fluorochrome using many standard methods of FMT reconstruction, including t he method of projections, singular value decomposition, and SVD-KR. More generally, t he solution can be written as n (86) The resulting matrix will be of size n by 2, where n is equal to the number of voxels in the forward problem. Each column of r, is a vector representing the concentration of the fluorochrome with fluorescent wavelength r. 6.1. Results A phantom was created from the Digimouse dataset, a nd two different types of fluorochromes were synthetically injected. To test the ability of imaging multiple fluorochromes to enhance the image quality for low contrast targets, the fluorochromes were assumed to have a concentration 4 times greate r in their target than the surrounding tissue. The low contrast makes identification of t he target difficult in the reconstructed image; however, when two or more fluorochromes are imaged simultaneously, the background can be canceled out to increase the cont rast of the target.
64 A measurement dataset with signal mixing was genera ted using the phantom. The dataset was reconstructed using row-wise SVD-KR reconstruct ion. The results for a slice of the reconstruction that contains both targets is shown in Figure 33. The results indicate that imaging using multiple fluorochromes can allow for enhanced detection of low contrast targets.
65 Original Reconstruction Fluorochrome 1 Fluorochrome 2 Fluorochrome 1 and 2 Combined Variations Colored Figure 33. Multi-fluorochrome Imaging with Contrast Enhancement
66 7. Parallel Processing Implementations Parallel computing allows for multiple processors t o simultaneously carry out calculations. This allows for increased reconstruc tion speed and for reconstructions of large numbers of voxels using large imaging dataset s. Parallel reconstruction algorithms have the potential to increase reconstruction resol ution, reconstruction image quality and imaging volume. 7.1. Parallel Reconstruction by SVD-KR SVD-KR reconstruction can be implemented in paralle l environments, allowing for decreased reconstruction times. To evaluate the p erformance of SVD-KR algorithms in a parallel environment, reconstructions were perfo rmed with different numbers of processor cores in a shared memory environment. Th e results of these reconstructions are shown in Table 7 and Figure 34 below. These re sults show that both row-wise SVDKR and column-wise SVD-KR see performance gains in a parallel environment. However, because of communication overhead and a se rial SVD algorithm, the gains are not proportional to the number of processor cores u sed. If implemented using a parallel SVD algorithm, larger speed increases would be exp ected.
67 Table 7. SVD-KR Parallel Reconstruction Time and Re lative Speed Voxels SourceDetector Pairs Sources Detectors Processor Cores Row-Wise SVD-KR Reconstruction Column-Wise SVD-KR Reconstruction Time Speed-up Time Speed-up 1728 331776 576 576 1 66.6 100% 52.9 100% 1728 331776 576 576 2 51.4 129% 39.2 135% 1728 331776 576 576 3 48.0 139% 35.2 150% Figure 34. SVD-KR Parallel Reconstruction Relative Speed 7.2. Two-Stage Approach for Larger Imaging Volumes A two-stage approach could be utilized to reconstru ct larger imaging volumes. In this approach, the imaging volume is initially reconstru cted with low resolution, and the resolution of a selected region is subsequently enh anced by a second reconstruction. A simple 2-D example of region of interest enhancemen t is shown in Figure 35. n !n r r
Figure This algorithm could be parallelized by having each computer independently refine a different region of the imaging chamber. Once all of the regions in where enhanced, they could be combined into a Attempts made to create a working algorithm ran int o problems resolution in the region of interest created a twof old problem: the condition n weight matrix significantly increased, while the regions increased the model erro error, reconstruction error quality. This algorithm is presented here it will allow someone to improve 68 Figure 35. Mesh Refinement in Region of Interest This algorithm could be parallelized by having each computer independently refine a different region of the imaging chamber. Once all of the regions in the imaging chamber where enhanced, they could be combined into a single high resolution image. Attempts made to create a working algorithm ran int o problems resolution in the region of interest created a twof old problem: the condition n weight matrix significantly increased, while the low resolution reconstruction of other increased the model erro r. With the increase in both condition number and reconstruction error s increased significantly leading to an overall decrease in image quality. This algorithm is presented here for the interest of the reader with the hope improve upon it in the future. This algorithm could be parallelized by having each computer independently refine a the imaging chamber high resolution image. Increasing the resolution in the region of interest created a twof old problem: the condition n umber of the reconstruction of other number and model overall decrease in image with the hope that
69 7.2.1. Implementing SVD Region of Interest Enhancem ent To enhance the resolution in the region of interest the matrix W will need to be updated; consequently, the SVD of W will need to be updated for reconstruction. Calcu lating the singular value decomposition of the matrix W as shown in (30) can be intensive in both computational and memory requirements for the large overdetermined matrices needed for a high resolution FMT reconstruction. Since th e SVD of the matrix W is already known from the initial low resolution reconstructio n, the SVD of W can be updated directly, as opposed to updating W and recalculating its SVD; this can significantly reduce the computational and memory complexity of r econstruction when W is overdetermined. 7.2.2. SVD Column Removal Update Refining the reconstruction resolution in a small r egion of interest requires the modification of the matrix W after initially calculating the singular value dec omposition. To remove the original voxels in this region from t he forward model, columns of data need to be removed from the matrix W and the singular value decomposition of this new matrix WÂ’ is calculated. Instead of recalculating the singu lar value decomposition of the matrix WÂ’ from the entire matrix W' the known singular value decomposition of W can be updated with the newly appended data, significan tly reducing the computational complexity. This method will work for removal of an y set of arbitrary columns, however to simplify the formalism without loss of generalit y a set of sequential columns will be removed in this example. Given a matrix W of size m by n, for which the rank-r singular value decomposition is known,
70 I (87) If an array of columns of size m by is removed from W (88) (89) Figure 36. Removing Columns from the Weight Matrix The known singular value decomposition of W can be updated to obtain the singular value decomposition of matrix WÂ’ I (90) First, remove the columns of VT that correspond to the columns of D these columns are at the same indexes in the matrix VT as the matrix W I I I $ I (91) m n n
71 I I I (92) This operation causes a loss of orthogonality in th e matrix I.. To restore orthogonality, the matrix I. can be decomposed into an orthogonal and an upper triangular matrix by the QR decomposition, e % /0 f g I (93) where Q is an orthogonal basis of I., and R is an upper triangular matrix. WÂ’ can then be rewritten as I e% % e # e (94) Where # %. U and QT are both orthogonal matrices however, A is not diagonal. To make A diagonal, its singular value decomposition is calc ulated, # % % I % (95) WÂ’ can then be written as % % I % e % % e I % (96) The final formulation for WÂ’ is therefore I (97) where % (98) % (99) I e I % (100) This method is to the best of my knowledge a new re sult. The advantages of this method is that it allows the singular value decomposition of WÂ’ to be determined to a high degree of accuracy by calculating the singular value decom position of A Since the size of A can
72 be significantly smaller than WÂ’ this method can significantly decrease the comput ational and memory complexity of updating the singular valu e decomposition of a matrix when columns of data are removed. 7.2.3. SVD Column Addition Update Refining the reconstruction resolution in a small r egion of interest requires the modification of the matrix W after initially calculating the singular value dec omposition. To incorporate the new voxels in this region into t he forward model, columns of data need to be appended to the matrix W and the singular value decomposition of this new matrix WÂ’ is calculated. Instead of recalculating the singu lar value decomposition of the matrix WÂ’ from the entire matrix WÂ’ the known singular value decomposition of W can be updated with the newly appended data, significa ntly reducing the computational and memory complexity. Given a matrix W of size m by n for which the rank-r singular value decomposition is known, I (101) If an array of new columns T of size m by is appended to W T (102)
73 Figure 37. Appending Columns to the Weight Matrix the known singular value decomposition of W can be updated to obtain the singular value decomposition of matrix WÂ’ I (103) First, let  h T (104) i M T T h (105) L is then the projection of C onto the orthonormal basis U and H is the component of C orthogonal to U Next, find the QR decomposition of H , r j /0 f g i (106) where J is an orthogonal matrix, and K is an upper triangular matrix. WÂ’ is equal to  m n n n'
74 T r $ h j ( $ I M ( (107) Upon inspection, the left and right matrices in the matrix product are orthogonal, however the middle matrix, denoted Q is not diagonal. To make Q diagonal requires finding its singular value decomposition , / / I / 12$ f k g e where e $ h j ( (108) The updated singular value decomposition is therefo re  I (109) where  r / (110) / (111) I $ I M ( I / (112) This procedure takes lD n E time, most of which is for the matrix multiplications that rotate the subspace shown in ( 101), (102) and (103) .
75 8. Conclusions Since SVD-KR reconstruction significantly reduces t he computational complexity and memory requirements for reconstruction of overdeter mined imaging datasets, imaging systems with larger detector arrays, higher spatial sampling, and a larger number of sources can be reconstructed; this improves the inf ormation content of the measurements and decreases the ill-posedness of the inverse prob lem, leading to increased resolution and accuracy for in vivo FMT imaging . Compared to reconstruction with SVD, the SVD-KR reconstruction method decreased reconstructi on time up to 25 times and decreased memory usage by up to three orders of mag nitude. Consequently, SVD-KR reconstruction allows for fast, high resolution rec onstructions with low reconstruction error. To make best use of this new reconstruction method, FMT imaging systems would be designed to have a large number of sources and d etectors. Future research will be necessary to design a FMT imaging system that takes advantage of SVD-KR for reconstructions of large detector arrays with a lar ge number of sources. The SVD-KR reconstruction method can be used for a variety of forward models. In one compatible model, CT imaging data from a combined C T/FMT system could be segmented into separate tissues; the known optical properties these tissues along with their 3-D shape could allow for more accurate appro ximations to the photon density in vivo Paired with improvements in the accuracy of the forward model, SVD-KR could
76 allow for FMT reconstructions of sufficient resolut ion and quality to be clinically meaningful, and significantly expand the applicatio ns of FMT.
77 9. List of References  R. B. Schulz, J. Ripoll, and V. Ntziachristos, "Exp erimental Fluorescence Tomography of Tissues With Noncontact Measurements, IEEE Transactions on Medical Imaging vol. 23, no. 4, pp. 492-500, 2004.  E.E Graves, J. Ripoll, R. Weissleder, and V. Ntziachr istos, "A submillimeter resolution fluorescence molecular imaging system fo r small animal imaging," Medical Physics vol. 30, no. 5, pp. 901-911, 2003.  W. F. Cheong, S. A. Prahl, and A. J. Welsh, "A Review of the Optical Properties of Biological Tissues," IEEE Journal of Quantum Electronics pp. 2166-2185, 1990.  M. S. Patterson, B. Chance, and B. C. Wilson, "Time resolved reflectance and transmittance for the noninvasive measurement of ti ssue optical properties," Applied Optics vol. 28, no. 12, pp. 2331-2336, 1989.  V. Ntziachristos and R. Weissleder, "Experimental t hreedimensional fluorescence reconstruction of diffuse media by use of a normali zed Born approximation," Optics Letters vol. 26, no. 12, pp. 893-895, 2001.  J. P. Culver et al., "Threedimensional diffuse optical tomography in the paral lel plane transmission geometry: Evaluation of a hybrid frequency domain continuous wave clinical system for breast imaging," Medical Physics vol. 30, no. 2, pp. 235247, 2003.  M. Firbank, M. Hiraoka, M. Essenpreis, and D. T. De lpy, "Measurement of the optical properties of the skull in the wavelength r ange 650-950 nm," Phys. Med. Biol. vol. 38, pp. 503-510, 1993.  T.C. Zhu and J. Lee, "Finiteelement modeling of light fluence distribution in prostate during photodynamic therapy," in Excerpt from the Proceedings of the COMSOL Multiphysics User's Conference Boston, 2005.  M. Schwieiger, S. R. Arridge, M. Hiraoka, and D. T. De lpy, "The finite element method for the propagation of light in scattering m edia: Boundary and source conditions," Med. Phys. vol. 22, no. 11, pp. 1779-1792, 1995.
78  J.B. Fishkin and E. Gratton, "Propagation of photon density waves in strongly scattering media containing an absorbing semiinfinite plane bounded by a straight edge," J. Opt. Soc. Am. A. vol. 10, no. 1, pp. 127-140, 1993.  A. C. Kak and M. Slaney, Principles of Computerized Tomographic Imaging New York: IEEE Press, 1988.  J. P. Culver, V. Ntziachristos, M. J. Holboke, and A. G. Yodh, "Optimizationzation of optode arrangements for diffuse optical tomograp hy: A singular value analysis," Opt. Lett. vol. 26, pp. 701-703, 2001.  T. Kolda, "Multilinear operators for higher-o rder decompositions," Sandia National Laboratories, Albuquerque, New Mexico, 2006.  A. K. Kaw, Introduction to Matrix Algebra Second Edition ed. Tampa, 2008.  J. P. Culver, V. Ntziachristos, M. J. Holboke, and A. G. Yodh, "Optimization of opto de arrangements for diffuse optical tomography: A s ingular value analysis," Opt. Lett. vol. 26, pp. 701-703, 2001.  R. C. Haskell et al., "Boundary conditions for the diffusion equation in radiative transfer," J. Opt. Soc. Am. A vol. 11, no. 10, pp. 2727-2741, Oct. 1994.  A. Sourbret, J. Ripoll, and V. Ntziachristos, "Accu racy of Fluorescent Tomography in the Presence of Heterogeneities: Study of the No rmalized Born Ratio," IEEE Transactions on Medical Imaging vol. 24, no. 10, pp. 1377-1389, 2005.  D. L. Pham, C. Xu, and J. L. Prince, "Current Metho ds in Medical Image Segmentation," Annual Review of Biomedical Engineering vol. 2, pp. 315337 August 200.  B. Dogdas, D. Stout, A. Chatziioannou, and R. M. Le ahy, "Digimouse: A 3D Wh ole Body Mouse Atlas from CT and Cryosection Data," Phys Med Biol. vol. 52, no. 3, pp. 577-87, Feb 2007.  D. Stout et al., "Creating a whole body digital mou se atlas with PET, CT and cryosection images.," Molecular Imaging and Biology vol. 4, no. 4, p. S27, 2002.  M. Brand, "Incremental Singular Value Decomposition Of Uncertain Data With Missing Values," in In ECCV 2002, pp. 707-720.
79  B. C. Wilson and S. L. Jacques, "Optical Reflectanc e and Transmittance of Tissues: Principles and Applications," IEEE Journal of Quantum Electronics vol. 26, no. 12, pp. 2186-2197, 1990.  W. Cong, K. Durairaj, L. V. Wang, and G. Wang, "A B orntype approximation method for bioluminescence tomography," Medical Physics vol. 33, no. 3, p. 679, 2006.
About the Author Stephen J. Shamp was born in Akron, OH and earned a B.S. degree in electrical engineering from the University of Notre Dame. He is currently a candidate for a M.S. degree in electrical engineering at the University of South Florida. After graduating, he is going to pursue a M.S. degree in Physiology at t he University of Cincinnati, with hopes of continuing thereafter to medical school.