University of Southern California Dissertations and Theses

Hyperspectral and multispectral optical bioluminescence and fluorescence tomography in small animal imaging

(USC Thesis Other)

Hyperspectral and multispectral optical bioluminescence and fluorescence tomography in small animal imaging

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content HYPERSPECTRAL AND MULTISPECTRAL OPTICAL BIOLUMINESCENCE AND FLUORESCENCE TOMOGRAPHY IN SMALL ANIMAL IMAGING by Abhijit J. Chaudhari A Dissertation Presented to the FACULTY OF THE GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulfillment of the Requirements for the Degree DOCTOR OF PHILOSOPHY (ELECTRICAL ENGINEERING) May 2007 Copyright 2007 Abhijit J. Chaudhari Dedication To those who were and are To those who were and are not For him And also, for her... ii Acknowledgments ‘Interdependence is of higher value than independence. ’(Stephen Covey). In 2000, I got my first opportunity to work on the human genome project under the guidance of Dr. Deborah K. Van Alphen of California State University, Northridge. My task, then, was to analyze gene expression imaging data and develop tools that help in drawing statistical inferences about gene pathways. Dr. Van Alphen made me write my first grant proposal, and with the support of Dr. Mack Johnson, had me present my research work at conferences and state-wide research competitions. While working with imaging data, I got interested in the field of molecular imaging and decided to pursue higher studies in the field of medical imaging. I consider myself very lucky to be initiated into the field of imaging by such passionate people. Dr. Richard Leahy’s research group at USC was a perfect fit for my aspirations of pursuing a doctoral degree. While my knowledge of theory was sound, I lacked practi- cal experience. Dr. Leahy gave me the opportunity to get some hands-on experience at Dr. Desmond Smith’s lab at UCLA. Working at a pharmacology lab was quite unchar- acteristic for me as an engineer, however, this experience was extremely useful in the later part of my research work when I had to work with small animals. I wish to thank Dr. Desmond Smith for the opportunity and technician Arshad Khan for the extensive amount of practical knowledge he shared with me. iii Apart from being excellent at mentoring and as an adviser, Dr. Richard Leahy gave me a lot of freedom to try new things. He supported innovation and kept me focused. I have benefited greatly from his command on imaging concepts and his urge for perfec- tion, and I wish to thank him for the values that he instilled in us members of his research group. Outside research, Dr. Leahy’s humanitarian qualities and his advice have helped me a lot in my personal life. I would also like to thank Dr. Wlodek Proskurowski, my Masters Degree adviser from the Department of Mathematics at USC for his help and support. His courses at USC and discussions with him helped me tremendously in developing a good understanding of complex mathematical concepts. Anand Joshi, my colleague and officemate, was always very inspirational. His abil- ities in grasping difficult concepts easily and in problem solving were unparalleled. Dr. Felix Darvas’ clear explanation of concepts and his practical approach to things was very helpful and is highly appreciated. Discussions with Quanzheng Li, Belma Dogdas, Sangeetha Somayajula, Sangtae Ahn, all from the Department of Electrical Engineering at USC, and with Ryan Park, Michel Tohme and Dr. James Bading of the Department of Radiology at USC were very fruitful. Maitreyee Tripathi provided critical inputs and encouragement. I wish to thank all these people. I would also like to express my grati- tude to Dr. Alexander Sawchuk and Dr. Antonio Ortega whose courses at USC helped me a great deal in coming up with, and pursuing novel ideas presented in this thesis. Lastly, I would like to thank my parents Dr. Saudamini Chaudhari and Dr. Jayawant Chaudhari who have been a constant source of support, encouragement and inspiration. iv Table of Contents Dedication ii Acknowledgments iii List of Tables viii List of Figures ix Abstract xiv 1 Introduction 1 1.1 Review and scope . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Contributions of this work . . . . . . . . . . . . . . . . . . . . . . . . 4 1.3 Organization of this thesis . . . . . . . . . . . . . . . . . . . . . . . . 7 2 Forward models in multispectral and hyperspectral OBT and OFT 9 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.2 The diffusion equation: Approximations and boundary conditions . . . 10 2.3 Mouse-atlas based volumetric registration with surface constraints . . . 13 2.4 Monochromatic forward model formulations . . . . . . . . . . . . . . . 16 2.5 Comparison of the analytic model and FEM . . . . . . . . . . . . . . . 19 2.6 Formulation of the Multispectral and Hyperspectral Forward Model . . 20 2.7 Illumination in HOFT . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 2.8 Reconstruction for multiple probes . . . . . . . . . . . . . . . . . . . . 26 2.9 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 3 Investigation of fast forward model solvers 28 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 3.2 Sparsity of our FEM matrix and reordering . . . . . . . . . . . . . . . 30 3.3 Cholesky decomposition and the fill-in phenomenon . . . . . . . . . . 32 3.4 Solution schemes for forward models . . . . . . . . . . . . . . . . . . 34 3.4.1 Solution by sequential back-substitution . . . . . . . . . . . . . 34 3.4.2 The incomplete Cholesky factor as a preconditioner for PCG . . 35 v 3.4.3 Solution using lumped RHS . . . . . . . . . . . . . . . . . . . 37 3.5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 4 Reconstruction methods for multispectral and hyperspectral OBT and OFT 41 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 4.2 Mathematical formulation of the objective function . . . . . . . . . . . 42 4.2.1 Measure of the fit to the data . . . . . . . . . . . . . . . . . . . 42 4.2.2 Regularization penalty . . . . . . . . . . . . . . . . . . . . . . 43 4.2.3 Positivity Penalty . . . . . . . . . . . . . . . . . . . . . . . . . 43 4.2.4 The objective function . . . . . . . . . . . . . . . . . . . . . . 44 4.2.5 The optimization algorithms and preconditioner . . . . . . . . . 45 4.3 Using spatial basis functions for dimensionality reduction . . . . . . . 46 4.3.1 Choice of basis functions . . . . . . . . . . . . . . . . . . . . . 47 4.3.2 Simulation setup and methods . . . . . . . . . . . . . . . . . . 49 4.3.3 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . 51 4.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 5 Experimental setups and the imaging process 54 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 5.2 Comparison of single-view data acquisition versus using four orthogo- nal views . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 5.3 The mirror setup for tomographic data acquisition . . . . . . . . . . . . 56 5.4 Phantom preparation and materials for phantom studies . . . . . . . . . 58 5.5 An unified OBT and OFT process chain . . . . . . . . . . . . . . . . . 60 5.6 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 6 Simulation and experimental studies 62 6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 6.2 Simulation studies in a mouse atlas using hyperspectral data for investi- gating resolution limits and localization error . . . . . . . . . . . . . . 63 6.3 An investigation of the impact of skull thickness and skull optical prop- erties on reconstructions of brain tumors using OBT . . . . . . . . . . . 69 6.4 A multi-modality PET-OBT study for comparing the localization of tumors in an inhomogeneous mouse phantom . . . . . . . . . . . . . . 73 6.5 Simultaneous reconstruction of sources with different emission spectra and applications using quantum dots . . . . . . . . . . . . . . . . . . . 74 6.6 A feasibility study exploring the use of OFT in human breast cancer imaging . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 6.7 Phantom studies to analyze localization error and resolution in OBT . . 79 6.8 An in-vivo study in a mouse . . . . . . . . . . . . . . . . . . . . . . . 82 6.9 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 vi 7 Summary and future work 86 7.1 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 7.2 Future work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 7.2.1 Innovative methods for spatial encoding of probes in OFT . . . 88 7.2.2 Fast forward and inverse solvers . . . . . . . . . . . . . . . . . 89 Bibliography 90 vii List of Tables 3.1 Table showing the sample matricesF used for this thesis . . . . . . . . 29 3.2 Table shows the reordering and decomposition costs in seconds for three sample matrices in HOBT . . . . . . . . . . . . . . . . . . . . . . . . 32 3.3 Table shows the backsolving cost in seconds for a single sample RHS . 32 3.4 Table shows the backsolving cost for a single sample RHS for the matrix withv = 59332. MATLAB’s cholin operator was used for the incom- plete Cholesky decomposition. For the complete Cholesky decomposi- tion, the same function was used with droptol = 0. . . . . . . . . . . . . 33 3.5 Table shows the cost of repeated reordering, Cholesky decomposition, and solving versus pre-computation of the Cholesky factor and sequen- tial backsolving for the matrix withv = 59332. Also, see section 3.4.3. 35 3.6 Table shows the cost for PCG to converge as a function of droptol for the incomplete Cholesky decomposition for the matrix of sizev = 59322 36 3.7 Table shows the cpu for varying number of RHS using sequential oper- ations, PCG and lumped RHS for the matrix of sizev = 211372. . . . . 38 3.8 Table shows the CPU time for all three matrices under consideration and 50 RHS when using sequential operations, PCG and lumped RHS . . . 39 4.1 The table shows values that were obtained with noiseless data and using the DCT basis functions. . . . . . . . . . . . . . . . . . . . . . . . . . 51 4.2 The table shows values that were obtained with noiseless data and using the DaubechiesD 4 basis functions. . . . . . . . . . . . . . . . . . . . . 51 4.3 The table shows values that were obtained with noiseless data and using the DaubechiesD 4 basis functions. . . . . . . . . . . . . . . . . . . . . 53 6.1 Table showing the FWHM at various locations for point-like sources placed in the rectangular phantom . . . . . . . . . . . . . . . . . . . . 82 viii List of Figures 2.1 Labeled sagittal images of (a) the atlas, (b) the subject and (c) the warped atlas. The warping was carried out using only a surface matching con- straints. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.2 Optical coefficients estimated from the warped atlas on the tetrahedral mesh of the subject. Figures show (a) diffusion coefficients and (b) absorption coefficients for the warped atlas. . . . . . . . . . . . . . . . 16 2.3 Schematic showing the process of computation of the analytic model with the extrapolated boundary . . . . . . . . . . . . . . . . . . . . . . 17 2.4 Correlation analysis between the photon flux computed on the mouse surface by the FEM-based solver using inhomogeneous optical proper- ties defined based on the mouse atlas and the analytic model assuming homogeneous optical properties throughout the animal volume. (a,b): plots of the correlation coefficient between the FEM and analytic for- ward field sensitivities computed for each surface element. (c-f): plots of forward fields for a point source computed using analytic (c,e) and FEM (d,f) for worst case (c,d) and best case (e,f) correlation of surface fluence. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 2.5 Variation of spectra with depth and SVD; (a) Variation in the spectrum of the Vybrant DiD dye with change in depth in tissue from 1 mm to 8 mm; (b) The singular value decomposition of the spectrum vs depth matrix for 8 different depths. Clearly, only 3 singular values carry more than 97% of the energy in the SVD spectrum. . . . . . . . . . . . . . . 21 ix 2.6 Simulated data obtained from a shallow (depth 2 mm) bioluminescent point source and a deep (depth 7 mm) point source after projection onto the first three principal spectral basis functions. (a,b,c) The three orthog- onal spectral basis functions; (d,e,f) Forward fields on the surface gen- erated by a source 2 mm deep; (g,h,i) Forward fields from a source at 7 mm depth. For a shallow source, the contributions from the second and third basis functions are very small. For deeper sources, the contri- butions from the second and third basis functions are more significant and help to reduce the ill-posedness of the inverse problem. It is note- worthy that the forward fields resulting from projection onto the first principal basis vector (a) are not equivalent to a simple integration over the spectrum of the bioluminescent source. . . . . . . . . . . . . . . . . 23 3.1 Pattern of non-zeros in the FEM matrices F of (3.2) before and after reordering. Each row in the figure grid corresponds to matrices of the same size. The first row corresponds to the matrix of size 59332, the second to that of size 146872 and the third to that of size 211372. The columns correspond to the original FEM matrix, the FEM matrix reordered by AMD, and the FEM matrix reordered by RCM respectively. 31 3.2 Plot of the cost for computing the incomplete Cholesky factorization and for PCG, as a function of droptol for the matrix of sizev = 59322 . 37 3.3 Cost of solving multiple RHS using sequential operations, PCG and using lumped RHS for the matrix of size v = 211372 as a function of the number of RHS. . . . . . . . . . . . . . . . . . . . . . . . . . . 38 4.1 Coefficients of the 3D DCT and the 3D Daubechies wavelet; (a) A sam- ple fluorescence image plotted on coronal sections of a cube-shaped tis- sue phantom; (b) The 3D DCT coefficients; (c)The 3D DaubechiesD 4 coefficients. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 4.2 Noiseless data reconstructed using variable number of DCT basis func- tions; (a) Using 512 basis functions; (b) Using 1728 basis functions; (c) Using the complete basis . . . . . . . . . . . . . . . . . . . . . . . . . 50 4.3 Noiseless data reconstructed using variable number of Daubechies D 4 basis functions; (a) Using 512 basis functions; (b) Using 1728 basis functions; (c) Using the model as is . . . . . . . . . . . . . . . . . . . 51 x 5.1 The conditioning of the forward model and localization error. (a) The singular value spectra of the forward models. (b) The localization error if data from multiple views is used for reconstruction. (c) The localiza- tion error when using data from a single-view. The sources were placed along the line of intersection of the x = 20 and y = 20 planes. The detectors for case (b) are on the z = 40, z = 0, x = 40 and x = 0 planes, whereas those for case (c) are on the z = 40 . . . . . . . . . . . 55 5.2 A photograph of the mirror setup designed for acquisition of tomo- graphic data. The setup has no moving parts and can be accommodated by any of the 2D Xenogen IVIS bioluminescence imaging systems. . . . 57 5.3 (a) Schematic of the mirror setup for acquisition of tomographic data; (b) A sample image obtained from the CCD camera using the mirror setup. The bottom and side-views are scaled relative to the top view because of differing optical path lengths; this is accounted for in our forward model. The focusing of the bottom view can be corrected by moving the subject platform vertically to the same focal plane as the top-view. Since the setup has no moving parts, the relative distances remain constant enabling a fixed correction for perspective. . . . . . . . 58 5.4 Homogeneous mouse and slab phantoms and inhomogeneous slab phan- toms designed and fabricated for experiments described below. The phantoms were prepared following the recipe from Cubeddu [CPT + 97]. 60 6.1 SVD analysis of the achromatic versus hyperspectral forward models. The color scale here is such that black indicates no spatial contribu- tion and white indicates maximum contribution: (a) SVD spectra of the achromatic and hyperspectral bioluminescent forward models with matched detectors; (b) Overlay of the absolute value of the right singular vector corresponding to the fourth largest singular value of the achro- matic forward model on the horizontal sections through the 3D mouse volume; (c) Overlay of the absolute value of the right singular vector corresponding to the fourth largest singular value of the hyperspectral forward model on the same sections. (b) and (c) were gamma-corrected for better visibility. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 xi 6.2 Reconstruction of a deep point-like source in an atlas-mouse geometry using hyperspectral data. The color scale used here goes from light (no intensity) to dark (maximum intensity). (a) Horizontal sections travers- ing the torso of the mouse and indicating the location of a point-like bioluminescent source used for simulation. (b) The reconstructed source distribution using monochromatic data collected at 2645 detector loca- tions and average optical properties over the 600-800 nm window.(c) The reconstructed source distribution using achromatic data collected at 2645 detector locations. (d) The reconstructed source distribution using hyperspectral data collected at 2645 detector locations and 100 spectral bins (corresponding to FWHM of 1.5 mm). The reconstruction from hyperspectral data is able to accurately localize the source distribution, but reconstruction from monochromatic and achromatic data indicate a broadly distributed source. . . . . . . . . . . . . . . . . . . . . . . . . 67 6.3 A plot of the FWMH resolution obtained from hyperspectral data for a reconstructed point source located deep in the abdomen as a function of the spatial sampling on the surface of the 3D animal volume. . . . . . . 68 6.4 A plot of the FWMH resolution obtained from monochromatic data using multiple illumination patterns for a reconstructed point source located deep in the abdomen as a function of the spatial sampling on the surface of the 3D animal volume. . . . . . . . . . . . . . . . . . . . 69 6.5 Wavelength dependent optical properties of the skull and the skin. Fig- ures (a) and (b) shows the absorption coefficient and the reduced scat- tering coefficients respectively. . . . . . . . . . . . . . . . . . . . . . . 70 6.6 Localization error as a function of wavelength λ and of perturbation in optical properties. (a) and (c) show the error due to perturbation of optical properties alone, where as (b) and (d) show the combined effect of perturbation in the optical properties and change in skull thickness. Smoothing was performed by applying a simple spatial Gaussian fil- ter of width 0.4 mm to the solution, to eradicate effects of outliers and numerical errors. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 6.7 The optical properties for organs in the inhomogeneous mouse atlas; (a)Reduced scattering coefficientsμ ′ s and (b) absorption coefficientsμ a 73 6.8 Overlay of the PET and optical reconstruction in a inhomogeneous mouse with tumors in the head, spleen and the lungs . . . . . . . . . . . . . . 75 6.9 The animal geometry and the emission spectra for the probes, (a) the three tumor locations in 3D and (b) the corresponding emission spectra . 76 xii 6.10 Horizontal sections of the mouse torso showing simultaneous recon- struction of the Texas red, Red beads and Cy5.5 biological probes in- vivo. The first row represents the original probe location, while the second row shows the reconstruction result. . . . . . . . . . . . . . . . 77 6.11 The tessellated human breast phantom . . . . . . . . . . . . . . . . . . 78 6.12 HOFT reconstruction of a point-like source at a depth of 6 cm in the human breast: (a) The original source location and (b) The reconstructed source distribution estimated using noisy data. . . . . . . . . . . . . . 79 6.13 Diagram of the locations where the optical fiber tip was placed along two orthogonal planes for the localization and resolution studies. The points are ordered as 1: (12, 11, 6) mm, 2: (12, 18, 6) mm, 3: (12, 6, 6) mm, 4: (4, 11, 6) mm, 5: (20, 11, 6) mm, 6: (12, 11, 10) mm and 7: (12, 11, 4) mm. The two points in the trans-axial X directions are indicated by gray circles, those along the trans-axial Y direction are denoted by white circles, those along the axial direction are denoted by black circles and the location close to the center of the phantom is shown as a white square. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 6.14 Achromatic and multispectral reconstructions using the rectangular-block shaped phantoms for a single point source. The color scale used here goes from light (no intensity) to dark (maximum intensity); (a) displays the horizontal sections through the phantom showing the true location (12, 11, 6) mm of the tip of the orange LED-driven optical fiber; (b) horizontal sections showing reconstructions from achromatic data; (c) horizontal sections showing reconstruction from multispectral data. . . . 81 6.15 Bioluminescence data obtained with the IVIS 200 imaging system and mirror setup overlaid on white-light photographs of a mouse with an implanted brain tumor. (1,8) are achromatic (i.e. unfiltered) acquisi- tions obtained for calibration purposes.(2)-(7) were acquired using the six filters (560, 580, 600, 620, 640 and 660 nm) mounted on the system in order of increasing wavelength. Since no signal was observed in the bottom view, only the top and side-views were acquired. . . . . . . . . 83 6.16 Reconstruction of six horizontal planes of the bioluminescence data overlaid on co-registered MR slices. The red contour indicates the bound- ary of the tumor on the corresponding slices. The threshold was set to display bioluminescence data greater than or equal to10% of the maxi- mum reconstructed voxel value. . . . . . . . . . . . . . . . . . . . . . 85 xiii Abstract For bioluminescence and fluorescence imaging studies in small animals, it is important to be able to accurately estimate the 3D distribution of a light-emitting source within the animal volume in order to draw biological inferences about underlying processes. The spectrum of light produced by an emission source that escapes the subject varies with the depth of the source because of the wavelength dependence of the optical prop- erties of tissue. Consequently, multispectral or hyperspectral data acquisition should help in the 3D localization of deep sources. In this thesis, we develop a framework for fully 3D bioluminescence and fluorescent tomographic image acquisition and recon- struction that exploits spectral information. Singular value analysis is used both for data dimensionality reduction and to illustrate the advantage of using hyperspectral rather than monochromatic or achromatic data. Regularized tomographic reconstruction tech- niques that use analytic or numerical solutions to the diffusion approximation of photon transport through turbid media were implemented and their properties were studied in phantom and animal studies. A fixed arrangement of mirrors was used in conjunc- tion with a CCD camera of a bioluminescence imaging system (Xenogen Corporation’s IVIS200) for simultaneous acquisition of multispectral imaging data over most of the surface of the animal. Phantom studies conducted using this system demonstrate our ability to accurately localize deep point-like sources and show that a resolution of 1.5 xiv to 2.2 mm for depths up to 6 mm can be achieved. Phantom and simulation studies car- ried out using a fluorescence imaging system (Cambridge Research Instrument’s MAE- STRO) demonstrate that sources upto 6 mm depth can be reconstructed using only a single view. For a mouse atlas geometry, we present a method for reconstructing multi- ple fluorescent probes simultaneously. We also show an in-vivo study of a mouse with a brain tumor expressing firefly luciferase. Our results indicate good anatomical localiza- tion of the tumor when the 3D bioluminescent image was co-registered with magnetic resonance images. While Optical Bioluminescence Tomography (OBT) and Optical Fluorescence Tomography (OFT) are limited to small animal imaging at this stage, a feasibility study using OFT for breast tumor detection in humans shows the potential use of these methods for clinical in-vivo studies. A. C. Los Angeles, California December 2006. xv Chapter 1 Introduction 1.1 Review and scope The ability to perform non-invasive in-vivo molecular and cellular-level imaging stud- ies of gene and protein expression and other cellular-level events, is revolutionizing our understanding of the biology of normal and diseased tissue and is playing an increas- ingly important role in areas of gene therapy, immunology, drug discovery, cancer research and treatment [WN03, CR02, MTBW99]. For most in-vivo studies in humans, optical imaging is largely limited to superficial sites owing to the absorbing and scat- tering properties of tissue, and MRI or PET are preferred modalities [BBM + 01]. How- ever, in small animals, due to shorter path lengths, a large fraction of photons reaches the surface of the animal and tomographic 3D reconstruction of bioluminescent and fluorescent signals is possible at significantly lower cost compared to MRI and PET [CDB + 05, RCN01, CR02]. For the imaging of bioluminescent sources in-vivo, no external excitation source is needed, and in turn, background noise is low and sensitiv- ity is high [WN03]. For in-vivo fluorescence imaging, the measured signal additionally depends on the quantum efficiency of the fluorescent probe, and on the emission inten- sity and the strength of the excitation source [CDB + 05, RCN01]. We define Optical Bioluminescence Tomography (OBT) and Optical Fluorescence Tomography (OFT) as techniques for gathering data over the animal surface with the objective of reconstruct- ing the 3D distribution of the bioluminescent or fluorescent emission sources when the optical properties of tissue are known. 1 The problem of determining the photon density on the animal surface from a bio- luminescent or fluorescent source distribution within the animal volume (called the forward problem) requires accurate representation of the photon transport in biolog- ical tissue. In addition, fluorescence imaging would require a model of light prop- agation from the external source to the emission site. A deterministic description of these processes under the assumption of isotropic scattering is given by the diffu- sion equation, which allows modeling of light propagation in turbid media with inho- mogeneous optical properties and realistic geometries [Arr99]. Several approaches have been proposed for solving the diffusion equation, including methods based on the discretization of the volume into finite elements that allow modeling of inho- mogeneities [ASHD93, SAHD95, RSM01, WLJ04] and methods based on analyt- ical approximations using simplified geometries and assuming tissue homogeneity [HST + 94, Arr95, RCN01, GRWN03]. For OBT/OFT to be applicable to animal studies, fast solvers for the diffusion equation are essential. Three-dimensional source reconstruction from bioluminescent or fluorescent data and the light diffusion-based forward model is a challenging inverse problem due to it’s ill-posedness. For fluorescence imaging, the external illumination source can be used to elicit different mappings from the fluorescent source to the surface of the animal improving the conditioning of the inverse problem. Bioluminescence imaging, on the other hand, is passive and one observes only a single mapping from the source to the surface, which depends on the spatial distribution of the source and the optical prop- erties of the surrounding tissue. In both cases, the spectral dependence of the optical properties of biological tissue can be exploited to potentially improve the conditioning of the inverse problem in the following manner. Bioluminescent or fluorescent photons emitted at near-infrared wavelengths undergo less attenuation compared to those in the blue and green parts of the spectrum due to the dependence of optical properties on 2 wavelength [CPW90]. As a consequence, the spectrum of the source measured on the surface of the animal will be modified depending on the depth of the emission site in tissue. A CCD camera that integrates photons over the spectrum of the emission source and measures intensity alone, would not be able to distinguish between a broad emis- sion pattern produced on the skin by a deep focal source and a superficial distributed source. In contrast, if hyperspectral (i.e. resolving the spectrum of the emitted light at a high sampling rate, e.g. ∼ 100 spectral bins) or multispectral (i.e. resolving the spectrum at a few wavelengths only, e.g. ∼ 8 spectral bins) detection techniques are employed, the two sources may be distinguished due to the difference in spectral mea- surements on the surface. Thus, the changes in spectra encode depth and the use of hyperspectral or multispectral data is expected to reduce the ill-posedness of the inverse problem and yield improved depth reconstruction compared to achromatic (i.e. collec- tion of integrated data without filters) or monochromatic (i.e. data collection at a single wavelength) measurements. It is important to note that in both OBT and OFT, the reconstruction problem is linear with respect to the source distribution. This differs from the traditional, nonlinear Diffuse Optical Tomography (DOT) problem, where spatial reconstruction of the optical properties of tissue is the primary objective. Inverse methods have been developed for solving the problems of DOT [YWP + 97, Arr99, UOO + 01], OFT [RCN01, RSM01, GRWN03] and combined DOT and OFT [KH03, MOW + 03]. A model-based approach for bioluminescent tomography has been proposed [GZLJ04], however, this method does not take into account the wavelength dependence of the opti- cal properties of tissue. Average values for the optical properties are used and as a con- sequence, a monochromatic model is assembled. These methods did not use spectral data, and hence produced smooth solutions and were not able to localize deep sources. While wavelength information potentially adds a new dimension to the OFT inverse 3 problem, so can the choice of the excitation source. Methods that use modulated exci- tation light [RSM01, CBDT05] have been proposed that exploit the lifetimes of the fluorescent probes and can reduce the ill-posedness of the inverse problem. Further, these approaches may be helpful in distinguishing multiple probes. Hyperspectral imaging has been used extensively in the fields of remote sensing and geology for identification of natural and man-made materials that are indistinguishable using standard color imagery [Lan02]. Biomedical applications of hyperspectral imag- ing are mostly limited to 2D spatial data [LH00, SNZ + 01, FVM + 01] or restricted to reconstructions on a 2D focal plane [FVM + 01]. Commercial instruments are available for multispectral bioluminescence and fluorescence imaging (e.g. Xenogen Corpora- tion’s IVIS 200 R , CRI Inc.’s Maestro R ), but again are restricted to 2D planar imag- ing. A recent software addition to the IVIS 200 system allows for 3D reconstruction from spectral data but is restricted to data from a single view and is based on a sim- plified homogeneous slab model of light propagation. Bioluminescence or fluorescence imaging using planar detectors and a single view is generally insufficient for accurate 3D reconstruction of the distribution of bioluminescent sources [WLJ04]. Data from multiple views may be gathered by placing the subject on a rotating stage and imaging with a CCD camera [GZLJ04]. Following our published results [CDB + 05] and those shown in this thesis, there have been efforts to incorporate spectral imaging into these systems [WSD + 06]. 1.2 Contributions of this work The main contributions of this work are summarized as follows: 4 Fast and accurate optical forward problem solvers Time-resolved Optical Absorption and Scattering Tomography (TOAST) [Arr99] is a pre-compiled software for assembling the monochromatic optical forward model. It has several limitations because it was written to be used for the Diffuse Optical Tomography (DOT) problem and requires substantial computing resources. For example, the assem- bly of a optical forward model with approximately 60,000 nodes took about 18 hours. We successfully have written our own Finite Element code that solves for the forward model efficiently. The new code assigns optical properties element-wise in the volu- metric tessellation of the animal. This is especially important where rendering a thin layer of tissue or bone e.g. the skull. Finite Element-based models are large and solving them can be computationally time consuming. We implemented methods that use sparse matrix reordering and factorization techniques to reduce the computational cost substan- tially. Software has also been developed for the analytic solver using the extrapolated boundary condition. Further, methods allowing fast assembly of the hyperspectral and multispectral optical forward models have been derived and implemented. Atlas-based registration scheme: Development and evaluation For OBT/OFT, the animal’s surface topography can be reconstructed from structured light measurements, but internal anatomy is unavailable unless additional CT or MR images are acquired. Following the work of Anand A. Joshi [JSTL06] in the area of surface-constrained human brain registration, we developed an analogous method for estimating the internal organ structure of a mouse by warping a labeled 3D volumet- ric mouse atlas with the constraint that the surfaces of the two should match perfectly [CJDL07]. We demonstrate and evaluate the application of this warping scheme in OBT/OFT, where scattering and absorption coefficients of tissue are functions of the 5 internal anatomy. Hence, better estimates of the organ structures potentially lead to a more accurate forward model resulting in improved source localization. Design and fabrication of an optical setup and development of the OBT/OFT imaging process chain For enabling 3D multispectral and hyperspectral optical data collection from most of the mouse surface, we have successfully designed, fabricated and used a mirror setup for several phantom and in-vivo experiments along with the Xenogen IVIS200 2D optical imaging system. Software was written to read the image files directly in MATLAB and correct them for perspective, blur and spherical aberrations. The collected data were warped onto the mouse skin assuming the mouse surface to be lambertian and by stitch- ing together the top, side and bottom views. Additionally, Agar phantoms that exhibit optical properties were fabricated in the shape of cubes and of a mouse. Molds for curing the phantom material were prepared from plaster of Paris. From an animal lab- oratory point of view, it is important to have a process chain for OBT and OFT, that requires minimum human interaction and maximum automation. This involves individ- ual modules that perform data collection, image registration and data processing. We have integrated all these modules into a single process as a part of this work. 3D tomographic reconstruction of optical signal and experimental evaluation Iterative methods (based on conjugate gradients) for image reconstruction of monochro- matic, achromatic, multispectral and hyperspectral data were successfully derived and 6 implemented in phantom experiments and in-vivo studies. The process involved map- ping data from the camera plane to the surface of the animal, registration of the mul- tiple views of the subject obtained from the mirror setup, solving a forward and then an inverse problem. The complete reconstruction cycle was successfully completed for several studies. Currently, sophisticated inversion methods are being investigated based on their noise, resolution and convergence properties. 1.3 Organization of this thesis In this thesis, we describe the three essential components for 3D multispectral and hyper- spectral OBT and OFT, namely, forward models in OBT and OFT, the inverse methods, and the instrumentation necessary for optical tomography. This thesis is organized as follows. Chapter 2 presents methods of solving the forward problem in OBT and OFT i.e. for modeling light propagation in tissue. Wavelength dependence of optical prop- erties is used to assemble the multispectral and hyperspectral models. Dimensionality reduction techniques using singular value decomposition are presented for both, the OBT and the OFT problems. For accurate forward modeling, a precise definition of the underlying optical properties is essential. An atlas-based registration method is pre- sented and evaluated for use when the internal geometry of the animal is unknown. In Chapter 3, fast solvers for the OBT/OFT forward model are discussed. A detailed description of reconstruction methods is given in Chapter 4. A novel method using spa- tial basis functions (3D Discrete Cosine Transforms and 3D Daubechies D-4 wavelets) is proposed for speeding up the inverse computations. Chapter 5 describes the mirror setup designed and fabricated by us that allows for simultaneous acquisition of data on most of the animal skin. This setup can be used with any 2D optical imaging system (e.g. Xenogen Corporation’s IVIS200 R ). In chapter 6, we describe five simulation studies 7 and three experimental studies that validate methods presented in Chapters 2, 3 and 4 using the instrumentation setup described in Chapter 5. Chapter 7 includes concluding remarks and future work. 8 Chapter 2 Forward models in multispectral and hyperspectral OBT and OFT 2.1 Introduction In order to solve the inverse problem in OBT or OFT for the bioluminescent source or the fluorescent probe distribution (q), one has to solve the forward problem first. The forward problem is to predict the photon fluence (b) on the boundary dΩ of the ani- mal volume Ω, given the distribution of the bioluminescent source or the fluorescent probe (q) and the tissue optical properties. For OFT, the forward problem, further, must include parameters dependent on the external illumination (typically it’s intensity, wave- length spectrum and nature i.e. point-like, distributed etc.). Tissue optical properties viz. the reduced scattering coefficientμ ′ s (r,λ), the absorption coefficientμ a (r,λ), and the refractive indexη(r,λ) are all functions of the 3D locationr in tissue, and wavelengthλ. In the near-infrared (NIR) region of the electromagnetic spectrum, the reduced scatter- ing coefficient is much larger than the absorption coefficient, and as a consequence, the diffusion approximation to the radiative transport equation (RTE) can be used [Arr99]. The diffusion equation is a linear parabolic partial differential equation that maps the source strengthq(r,λ) to the measurementsb(r,λ). In this chapter, we first describe the diffusion approximation to the radiative trans- port equation and the assumptions made therein. Further, we describe two methods to solve the diffusion equation, the first that uses the Finite Element Method (FEM) with 9 a Robin-type boundary condition, and the second ‘analytic’ method that uses a semi- infinite slab approximation with an extrapolated boundary. The FEM-based approach can incorporate a inhomogeneous tissue model, whereas, the analytic model works only with the homogeneous assumption. As a consequence, FEM-based methods provide fairly accurate models when the animal internal geometry is known. However, for an all-optical study, there is no way of estimating the internal organ geometry of the animal. To address this issue, we discuss a registration method of warping a volumetric mouse atlas to the animal to be imaged [CJDL07]. In the subsequent sections, we present comparisons of the FEM-based model with the analytic model and describe the assembly of the multispectral and hyperspectral optical forward model that uses the spectral subspace projection idea. The FEM-based solution to the diffusion equation is very slow. A discussion of fast FEM-solvers for OBT/OFT that use sparse matrix reordering and decomposition to accelerate the forward model computation is given in Chapter 3. Comparisons between the multispectral model and the achromatic model (obtained by integrating over the spectrum of the emission source) based on the singular value decomposition will be shown in Chapter 6. 2.2 The diffusion equation: Approximations and boundary conditions In 3D, the steady-state RTE for monochromatic light as a function of the position r assumes the form [AH03] ˆ s∇φ(r,ˆ s)+(μ a +μ s )φ(r,ˆ s) =μ s Z 4π f(ˆ s·ˆ s ′ )φ(r,ˆ s)dˆ s ′ +q(r,ˆ s) (2.1) 10 where the scalar fieldφ(r,ˆ s) denotes the energy radiance in W mm −2 sr −1 ,μ a andμ s are the absorption and scattering coefficients of the medium in mm −2 , the unit vectorˆ s is the direction of scatter for the photon, whose initial direction was the unit vectorˆ s ′ , f(ˆ s·ˆ s ′ ) denotes the scattering phase function which is the probability of scattering of a photon with initial directionˆ s ′ in the directionˆ s, andq(r,ˆ s) is the directional source density. The dot product in the argument of the scattering phase functionf(ˆ s·ˆ s ′ ) implies its dependence only on the cosine of the angle betweenˆ s andˆ s ′ . With photon densityΦ defined as Φ(r) = Z 4π φ(r,ˆ s)dˆ s (2.2) the diffusion equation in steady state at a single wavelengthλ is given as −{c▽·κ(r)▽−μ a (r)c} −1 q(r) = Φ(r), (2.3) where κ(r) = 1 3[μ a (r)+μ ′ s (r)] (2.4) κ(r) is the diffusion coefficient, c is the velocity of light at wavelength λ, and μ ′ s = μ s (1−g) is the reduced scattering coefficient withg as the average cosine of the scatter angle. The detailed derivation of the diffusion equation from the RTE is given in [CZ67]. While the diffusion equation is simply the first angular moment (P-1) approximation to the RTE, it has been shown to be valid under the following circumstances; 1. The medium is dominated by scattering i.e. μ ′ s >> μ a , whereμ ′ s is the reduced scattering coefficient andμ a is the absorption coefficient 2. The point of interest is at least a few scatter lengths (1/μ ′ s ) away from the source 11 3. The propagating signal cannot be approximated correctly if measured at very early times While condition 1 is satisfied by most biological tissue, there are non-scattering regions e.g. the lungs, where the diffusion approximation is not valid. While there has been considerable debate about what diffusivity (the ratio μ ′ s /μ a ) is required for the validity of the diffusion approximation [AH03, DLA + 99], a value> 10 is widely accepted [DLA + 99]. Methods have been proposed to use the RTE directly for solv- ing the forward problem [AH03], however, these methods are very slow for realistic geometries. The usual boundary condition for the RTE is that no photons travel in the inward direction at the boundarydΩ, i.e. withν as the outward normal, φ(r,ˆ s) = 0, whenr∈dΩ,ˆ s·ν < 0 (2.5) While (2.5) cannot be satisfied exactly by the diffusion equation, the condition is replaced by requiring that the total inward photon current be 0 at the boundary [Arr99, CZ67], i.e. Z ˆ s·ν<0 ˆ sφ(r,ˆ s)dˆ s = 0,r∈dΩ (2.6) To account for the refractive index mismatch of the air-tissue interface (refractive index of 1.4 for tissue versus 1 for air), and total internal reflection at the boundary, a parameter R is introduced [SAHD95, Arr99], and the boundary condition is modified to Z ˆ s·ν<0 ˆ sφ(r,ˆ s)dˆ s = Z ˆ s·ν>0 ˆ sRφ(r,ˆ s)dˆ s,r∈dΩ. (2.7) 12 After simplifications, equation (2.7) transforms to [Arr99] Φ(r)+2κG ∂Φ ∂ν = 0,r∈dΩ. (2.8) where κ = 1/[3(μ a +μ ′ s )]. The parameter G may be found experimentally as G = (1+R)/(1−R) [GFB83] where R≈−1.440n −2 int +0.710n −1 int +0.668+0.0636n int (2.9) or using Fresnel’s laws and Keijzer’s method [KSS88] as G = 2/(1+R 0 )−1+|cosθ c | 3 1−|cosθ c | 2 (2.10) whereθ c = sin −1 (1/n) is the critical angle, andR 0 = (n−1) 2 /(n+1) 2 is the Fresnel’s reflection coefficient. The boundary condition in (2.8) is called the Robin-type boundary condition. 2.3 Mouse-atlas based volumetric registration with sur- face constraints Atlases are normalized representations of anatomy that provide a standard coordinate system for in vivo imaging studies. For OBT/OFT in small animals, the animal’s sur- face topography can be reconstructed from structured light measurements, but the inter- nal anatomy is unavailable unless additional CT or MR images are acquired. Here, we present a novel method of using a mouse atlas to generate an approximate internal anatomy for a mouse on which OBT/OFT studies are performed. This is especially useful because the scattering and absorption coefficients of tissue are functions of the 13 internal anatomy. As a result, better estimates of the organ structures can lead to a more accurate OBT/OFT forward model and hence, to improved source localization. Methods The objective of the method presented here is to estimate the internal organ structure of a mouse by warping a labeled 3D volumetric mouse atlas with the constraint that the surfaces of the two should match perfectly. Further, point landmarks corresponding to anatomical features may be defined for the two surfaces and are matched. The important steps in this warping procedure are as follows. The mice volumes are treated as 3D Riemannian manifolds with the outer skin surfaces as their boundaries. We first compute the initial mapping of these surfaces to unit disks by minimizing the covariant Cauchy- Navier elastic energy which simultaneously computes the flat maps and co-registers the point landmarks in the flat space. The surface coordinates satisfy the elastic equilibrium relation of the form Δφ +∇(∇·φ) = 0. We then compute a volumetric harmonic map between two mouse volumes, while constraining the two surfaces to map to each other such that the point landmarks remain registered. This is done by minimizing the harmonic energyE(u) of the mapu [JSTL06]. Intermediate spherical coordinates are used to facilitate the surface and volume alignment. This novel approach allows us to match mouse surfaces as well as vol- umes in such a way that the surfaces of the two mice are exactly matched. Segmented and labeled sagittal sections through the atlas, the ‘subject’ and the atlas warped to the subject using the surface constraint are shown in Figure 2.1. For evaluation of the warping scheme in OBT, the 3D animal volume was tessellated. Tissue labels from the warped atlas were transferred to the corresponding tetrahedrons in the space of the ‘subject’ mouse and optical properties were assigned element-wise to 11 tissue types (brain, liver, kidneys, lungs, heart, muscle, bladder, skeleton, stomach, 14 Figure 2.1: Labeled sagittal images of (a) the atlas, (b) the subject and (c) the warped atlas. The warping was carried out using only a surface matching constraints. spleen, and pancreas) based on published results as shown in Figure 2.2. The biolu- minescence sources were assumed to have the emission spectrum of firefly luciferase, and a multispectral optical forward model (A) was assembled using an efficient finite- element solver as derived in the next section. The forward model (H) corresponding to the homogeneous mouse volume was generated by assigning muscle optical properties to all internal organs. To evaluate the accuracy of the two forward models, a ‘true’ model (T) was computed. This was done using labeled MRI images of the subject. For a sam- ple source located in the stomach, the correlation coefficient for surface data between A and T was 0.90, as against 0.74 between H and T. Clearly, A presented a better approx- imation to the true data, and can be expected to yield improved reconstruction results. This proposed registration method is being evaluated by the computation of dice coeffi- cient between organs and by reconstruction of representative sources in various organs. 15 Figure 2.2: Optical coefficients estimated from the warped atlas on the tetrahedral mesh of the subject. Figures show (a) diffusion coefficients and (b) absorption coefficients for the warped atlas. 2.4 Monochromatic forward model formulations The 3D domain is discretized into T tetrahedral elements, connected at v vertex nodes. The solution Φ is discretized by finite element (FE) basis functions to Φ h . The prob- lem of solving forΦ h becomes a sparse matrix inversion problem in the Finite Element Method (FEM) framework [ASHD93]. The solution for photon flux at v discrete loca- tions is given by the solution to FΦ =w, F∈R v×v , Φ,w∈R v (2.11) where F =K(κ)+C(μ a )+ηB (2.12) 16 where w is the load vector for a specific source configuration. The FEM matrix F is assembled as an addition of matrices K and C that depend on the diffusion and absorp- tion coefficients respectively, and a matrix B that imposes the boundary conditions. Assume that the data are acquired for m surface nodes with n internal source locations. Then, the monochromatic forward modelA(λ) is given by A(λ) =DF −1 W, A(λ)∈R m×n (2.13) where D is a detector selection matrix andW = [w 1 w 2 ... w n ] is a matrix of the load vectors corresponding to the sources under consideration. Our implementation of the FEM followed the steps outlined in [ADSO00]. We also used Time-resolved Optical Absorption and Scatter Tomography (TOAST) [Arr99], a finite-element solver specifically tailored to solve the diffusion equation for preliminary studies. Figure 2.3: Schematic showing the process of computation of the analytic model with the extrapolated boundary For simplified geometries and assuming a homogeneous animal volume, it is pos- sible to analytically approximate the diffusion process [HST + 94]. For a homogeneous and infinite tissue model, and where the diffusion approximation is satisfied, the photon fluence Φ at distance d from a point source with source powerS decays exponentially away from the point source according to Φ(d,λ) = S 4πκ(λ)d exp −μ eff (λ)d (2.14) 17 whereμ eff (λ) = p 3μ a (λ)[μ a (λ)+(μ ′ s (λ)] andκ(λ) are as defined in equation (2.4). The animal is assumed to be piecewise homogeneous, thus, we have dropped the depen- dence of the optical parameters on d in equation (2.14). To account for the finite dimen- sion of the animal and consequently, the refractive index mismatch of the tissue-air inter- face, we use the extrapolated boundary condition. The extrapolated boundary condition is the same boundary condition as for the FEM and requires that the inward photon cur- rent equals zero at an extrapolated boundary at heightz b outside the tissue surface. To achieve this, the fluence at z b due to a point source inside the tissue can be forced to zero by introducing a negative image source of photons as is illustrated in Figure (2.3) [HST + 94]. Consequently, the fluence at any point on the tissue surface is the sum of the photon fluence from the point source and it’s image. Ifd 1 andd 2 represent the distance from a surface location to the point source inside dΩ and it’s image respectively, the photon fluence at that surface location is given by Φ(d 1 ,d 2 ,λ) = S 4πκ(λ) " exp −μ eff (λ)d 1 d 1 − exp −μ eff (λ)d 2 d 2 # (2.15) wherez b is computed using Fresnel’s reflection coefficientR 0 [Hec01] as z b = 1+R 0 1−R 0 · 2 3[μ a (λ)+μ ′ s (λ)] This analytic solution of the diffusion equation has been shown to be in reasonable agreement with Monte-Carlo simulations assuming a homogeneous mouse [RCN01]. 18 Figure 2.4: Correlation analysis between the photon flux computed on the mouse sur- face by the FEM-based solver using inhomogeneous optical properties defined based on the mouse atlas and the analytic model assuming homogeneous optical properties throughout the animal volume. (a,b): plots of the correlation coefficient between the FEM and analytic forward field sensitivities computed for each surface element. (c-f): plots of forward fields for a point source computed using analytic (c,e) and FEM (d,f) for worst case (c,d) and best case (e,f) correlation of surface fluence. 2.5 Comparison of the analytic model and FEM The assembly of the hyperspectral, multispectral or achromatic forward models required the computation of the monochromatic models first. We have evaluated two solutions of the monochromatic forward model. The first was based on the numerical solution of the diffusion equation using the FEM and the second ‘analytic’ approach used the extrapolated boundary condition. Using published optical properties for biological tis- sue [CPW90] and a mouse atlas [SCS + 02], we constructed an inhomogeneous model for light propagation at wavelength 620 nm using the TOAST FEM code. The ‘analytic’ model used for comparison with the FEM solution was computed assuming a homoge- neous geometry with average optical properties at 620 nm. Figures 2.4(a) and 2.4(b) 19 show the correlation coefficients between the FEM and analytic forward models, com- puted across the set of source locations, for each of the 24324 surface nodes. The figure shows generally high correlations other than in the limbs. We also compared the mod- els by computing correlations across the set of surface nodes for each source location. The average correlation coefficient between the inhomogeneous mouse and the analytic approximation was R = 0.84. The lowest correlations were found in the head where the light source is completely enclosed by the skull (e.g. figures 2.4(c) and (d), R = 0.67) and the highest correlations were found in the dorsal region of the mouse, which is mostly soft tissue and muscle (e.g. figures 2.4(e) and (f), R = 0.96). This comparison shows that the analytic model and the FEM agree well in regions where optical properties are mostly homogeneous, e.g. in the trunk of the mouse, but that significant differences can be expected for sources shielded by materials which have optical properties that deviate strongly from the bulk optical properties of the animal. 2.6 Formulation of the Multispectral and Hyperspectral Forward Model Either the inhomogeneous FEM model or the homogeneous analytic model can be used to compute the monochromatic forward model. If we assume n possible source locations inside the animal and m detector locations on the surface then the measured datab(λ) can be modeled as: b(λ) =A(λ)q (2.16) whereA(λ)∈R m×n is the monochromatic forward model at wavelengthλ andq∈R n is the source distribution. Ifs 0 = [ s 0 (λ 1 ) s 0 (λ 2 ) ... s 0 (λ h )] T ,s 0 ∈R h represents 20 (a) (b) Figure 2.5: Variation of spectra with depth and SVD; (a) Variation in the spectrum of the Vybrant DiD dye with change in depth in tissue from 1 mm to 8 mm; (b) The singular value decomposition of the spectrum vs depth matrix for 8 different depths. Clearly, only 3 singular values carry more than 97% of the energy in the SVD spectrum. the emission spectrum of the bioluminescent source forh wavelength bins, the achro- matic model (A achr ∈ R m×n ) was obtained by integrating the monochromatic models over the source spectrum, i.e. A achr = h X i=1 s 0 (λ i )A(λ i ) (2.17) It is noteworthy thatA achr is not equivalent to constructing a monochromatic model by using average values for the optical properties. Assume that optical data are collected in h spectral bins over the wavelength range of interest. The multispectral forward model A mult was formed by concatenating the individual weighted monochromatic forward models for each wavelength. The measurementsA mult were concatenated similarly to obtain a system of the form b mult =A mult q (2.18) where b mult = h b T (λ 1 ) b T (λ 2 ) ··· b T (λ h ) i T ,b mult ∈R hm 21 and A mult = h s 0 (λ 1 )A T (λ 1 ) s 0 (λ 2 )A T (λ 2 ) ··· s 0 (λ h )A T (λ h ) i T ,A mult ∈R hm×n . Providedh is small, the full data set can be used for reconstruction of the 3D image. This is not the case for hyperspectral data, where measurements from 100-300 spectral bands may make the inverse problem intractable. For example, assume a reconstruction grid spacing of 1 mm inside the mouse volume with 20,000 source nodes and a minimum spatial sampling of 1000 pixels per image for each of 4 orthogonal views of the animal. Each monochromatic bioluminescent forward model will then have 20000×4000 ele- ments. Hyperspectral instruments are routinely configured to collect ∼ 100 spectral bands, resulting in a concatenated forward model with on the order of 10 10 elements (about 32 GB). Inversion of these models would be extremely time consuming. How- ever, dimensionality reduction can be achieved by projecting each monochromatic for- ward model onto a subspace defined by a set of spectral basis functions selected to best exploit the depth-dependent spectral characteristics of the data. Applying Principal Components Analysis, we observed that a few principal basis functions, defining a sub- space much smaller thanh, carry most of the depth information, as is shown in Figure 2.5. To investigate the use of a reduced dimensionality search, we placed point sources at p different depths inside a computer simulated homogeneous mouse volume. Let h represent the number of wavelength bins in which hyperspectral data were col- lected. Let s 0 = [ s 0 (λ 1 ) s 0 (λ 2 ) ... s 0 (λ h )] T ,s 0 ∈R h represent the emission spectrum of the bioluminescent source as a function of wavelength λ and s i ∈ R h , i = 1,2,...,p, denote the spectral measurement at the surface for the source at depth d i . Each s i for a fixed source depth d i is a function of λ, i.e. with fixed 22 Figure 2.6: Simulated data obtained from a shallow (depth 2 mm) bioluminescent point source and a deep (depth 7 mm) point source after projection onto the first three principal spectral basis functions. (a,b,c) The three orthogonal spectral basis functions; (d,e,f) Forward fields on the surface generated by a source 2 mm deep; (g,h,i) Forward fields from a source at 7 mm depth. For a shallow source, the contributions from the second and third basis functions are very small. For deeper sources, the contributions from the second and third basis functions are more significant and help to reduce the ill-posedness of the inverse problem. It is noteworthy that the forward fields resulting from projection onto the first principal basis vector (a) are not equivalent to a simple integration over the spectrum of the bioluminescent source. d i ,s i = [ s i (λ 1 ) s i (λ 2 ) ... s i (λ h )] T ,s i ∈R h . Surface measurements were made collinear with the p source locations. We constructed a source depth versus normalized surface spectral measurement matrix S (normalization is necessary to avoid dependence on the absolute power of the spectra): S = h s 1 ks 1 k s 2 ks 2 k ... sp kspk i ,S∈R h×p . 23 A SVD of matrix S yielded the principal spectral basis vectors u i and the principal depth contributing vectorsv i i.e. S = P p i=1 σ i u i v T i . We then constructed a reduced- dimension subspace from the t left singular vectors corresponding to the set of singular values that constitute more than 95% of the energy. The t left singular vectorsU = [ u 1 u 2 ... u t ] T ,U∈R h×t form a basis for this subspace. Each principal basis vector can individually be written as a function of wavelengthλ as u j = [ u j (λ 1 ) u j (λ 2 ) ... u j (λ h )] T ,u j ∈R h ,j = 1,2,...,t. (2.19) The first three principal spectral basis functions are shown in Figure 2.6. Also shown on the mouse surface are the projections of the data onto three spectral basis functions for shallow (2mm) and deep (7mm) sources. For a shallow source, the contributions from the second and third basis functions are small compared to the first, while for deeper sources the contributions from the second and third basis functions are larger and will help to resolve source depth. Assume n voxels inside the animal volume and m detector locations on the surface of the animal. LetA(λ i ) andb(λ i ) represent the monochromatic forward model and the measurement at wavelengthλ i respectively. The monochromatic forward model from (2.16) at wavelength λ i , with b(λ i ) ∈ R m ,i = 1,2,...,h can be written as b(λ i ) =s 0 (λ i )A(λ i )q. (2.20) Projecting the measurements and the forward model onto the jth spectral basis vector (j = 1,2,...,t), with ˜ b j ∈R m , ˜ A j ∈R m×n , andq∈R n , we obtain ˜ b j = ˜ A j q (2.21) 24 where ˜ b j = h X i=1 b(λ i )u j (λ i ) and ˜ A j = h X i=1 s 0 (λ i )u j (λ i )A(λ i ). We concatenated the measurements and forward models corresponding to projections on the t spectral basis vectors to form the linear system b hyp =A hyp q (2.22) where b hyp = h ˜ b T 1 ˜ b T 2 ··· ˜ b T t i T ,b hyp ∈R tm and A hyp = h ˜ A T 1 ˜ A T 2 ··· ˜ A T t i T ,A hyp ∈R tm×n . In the context of the formulation of hyperspectral models presented here, the achromatic model is equivalent to generating a hyperspectral model using only one rectangular spec- tral basis function that is equal to 1 over the bandwidth of the source. The model up to this point describes the bioluminescence phenomenon. For OFT, the forward model must incorporate light propagation from the external illumination source to the fluo- rophore location. The diffusion equation can be used to model this. We understand that the spectral basis functions also depend on the tissue in which the source is located. For example, the spectral basis functions might be very different for sources in the brain or when placed in the lungs. The validity of this method could further be investigated for more complex geometries. 25 2.7 Illumination in HOFT Each spatial illumination pattern corresponds to a single hyperspectral forward model. If we denote b k hyp and A k hyp (k = 1, 2 ,..., w) as the hyperspectral measurement and the hyperspectral forward model for the kth spatial illumination pattern respectively, the fluorescent data can be modeled as b ill hyp =A ill hyp q (2.23) where b ill hyp = h b 1 hyp b 2 hyp ··· b w hyp i T ,b ill hyp ∈R wtm and A ill hyp = h A 1 hyp T A 2 hyp T ··· A w hyp T i T ,A ill hyp ∈R wtm×n By an optimum choice of the spatial illumination pattern, the conditioning of the for- ward model can be improved substantially. The hyperspectral optical system we are developing [ZVM + 06] employs an external laser as an illumination source. The sys- tem provides flexibility to step the laser over a 2D grid, each time, varying the spatial illumination pattern on the surface. 2.8 Reconstruction for multiple probes We assume that we have k probes that correspond to k different emission spectra. Let the hyperspectral model (2.22) (with illumination in the case of OFT (2.23)) for each source spectrumi be represented byA i . For anith source spectrum, the underlying equation is b =A i x i (2.24) 26 whereA i ∈R wtm×n ,x i ∈R n andb∈R wtm . Therefore, for ‘k’ sources with different spectra, we have the system b = h A 1 A 2 ··· A k ih x 1 x 2 ··· x k i T (2.25) or b =A mp x mp (2.26) where A mp ∈ R wtm×nk , x mp ∈ R nk andb ∈ R wtm . Using this technique, multiple source distributions can be estimated at once. An application of this method is demon- strated in Chapter 6. 2.9 Summary Methods for the assembly of the forward models in OBT and OFT have been presented. The hyperspectral model is on the order of 100 times larger than the monochromatic model, however this increase in size can be reduced to a factor of only three by using spectral dimensionality reduction. Thus, the dimensionality reduction idea provides considerable cost saving and better conditioning of the forward model. The forward model assembly can been extended further to incorporate multiple illuminations and multiple source distributions. 27 Chapter 3 Investigation of fast forward model solvers 3.1 Introduction The assembly of the monochromatic optical forward model by the FEM-based solution to the diffusion equation requires the solution of a system of linear equations represented in the matrix-vector notation as Fφ =w, (3.1) where the matrixF ∈ R v×v , and vectorsφ,w ∈ R v . Solving this general problem by the LU decomposition [GV96] would cost a total of O(v 3 ) flops (i.e. the cost of the decompositionF =LU plus the cost of back-substitution). It is important to note that for HOBT, the problem (3.1) needs to be solved repeatedly with the same matrix F on the left hand side, but several right hand sides (RHS)w i ,i = 1,...,p: FΦ =W, (3.2) where the matrixF∈R v×v ,Φ,W∈R v×p withΦ = [φ 1 φ 2 ... φ n ] representing a matrix whose columns are the solutions and W = [w 1 w 2 ... w n ] . Instead of solving the system directly, the decomposition of F can be performed only once, followed by repeated back-solving. The cost of solving (3.2) by this method will be 28 O(v 3 +pv 2 ) instead of O(pv 3 ) for the direct solver without preprocessing. A straight- forward implementation of this scheme, however, is not possible, because the matrix F is large, withv often of the order of several hundred thousands, v = O(10 5 ), and the number of RHS’sw i of orderp = O(10 3 ). It is important to note the following: 1. The matrix F is a symmetric, positive definite (SPD) and sparse matrix with the number of non-zero elements (denoted as nz) per row of the order of tens, i.e. nz/v = O(10). 2. We are interested in the solution φ m only at the surface nodes and not at all v nodes. Thus, we do not desire to have the entireΦ in the computer memory. 3. Since w is a finite grid representation of the source location and strength, W is sparse as well (about 4 nz elements per column). Matrix size Non-zeros Density as Density as Condition v nz nz/v nz/v 2 number 59332 782853 13.19 2.2×10 −4 6.8×10 4 146872 2004522 13.64 9.2×10 −5 7.8×10 5 211372 3156406 14.22 7×10 −5 9.7×10 7 Table 3.1: Table showing the sample matricesF used for this thesis In this chapter, we discuss reordering schemes that use the SPD property of F and evaluate them for use with realistic HOBT matrices as are indicated in Table 3.1. Fur- ther, we analyze complete and incomplete Cholesky decompositions of reordered matri- ces. These decompositions are employed for the following three solution schemes which are discussed in detail in this chapter; 1. Reordering, complete Cholesky factorization of the reordered matrix, and repeated back-substitution. 29 2. Reordering, incomplete Cholesky factorization of the reordered matrix, and using the Cholesky factor as a pre-conditioner for the Preconditioned Conjugate Gradi- ent (PCG) method. 3. Reordering, complete Cholesky factorization of the reordered matrix, and using lumped RHS’s for back-substitution. 3.2 Sparsity of our FEM matrix and reordering The cost associated with solving sparse matrix problems depends not only on their spar- sity but also on the distribution of the non-zero elements. An accurate estimate of com- putational cost, namely O(vk 2 ) for decomposition and O(vk) for backsolving, where k is the half bandwidth,k<<v, may be given only for banded SPD matrices. The matrix F for HOBT is not banded and has an irregular distribution of the non-zero elements, thus, the cost cannot be estimated accurately. To minimize the fill-in phenomenon due to the decomposition of F, the following two sparse matrix reordering schemes may be used: 1. Symmetric reverse Cuthill-McKee ordering (RCM) [GL81]: This scheme finds a permutation that produces a matrix that has its non-zero elements close to the main diagonal compared to the original matrix, thus, reducing the bandwidth. 2. Symmetric approximate minimum degree permutation (AMD) [ADD04] that pro- duces a matrix with a sparser Cholesky factor compared to the original matrix. The reordered matrices often have the ’inverted fir tree’ shape along their diago- nal. 30 (a) 0 1 2 3 4 5 x 10 4 0 1 2 3 4 5 x 10 4 nz = 782853 (b) (c) (d) (e) (f) (g) (h) (i) Figure 3.1: Pattern of non-zeros in the FEM matricesF of (3.2) before and after reorder- ing. Each row in the figure grid corresponds to matrices of the same size. The first row corresponds to the matrix of size 59332, the second to that of size 146872 and the third to that of size 211372. The columns correspond to the original FEM matrix, the FEM matrix reordered by AMD, and the FEM matrix reordered by RCM respectively. Typical distributions of the non-zero elements in our matrix F before and after reordering are shown in Figure 3.1. We have used MATLAB’s implementation of func- tions for sparse matrix operations in simulations, experiments and presentation. 31 Size (v) Decomposition RCM AMD without reordering Reordering Decomp. Reordering Decomp. 59332 99.4 0.34 10.6 0.6 2.5 146872 490 0.9 24.2 1.7 7.5 211372 2588 1.4 130 2.8 18 Table 3.2: Table shows the reordering and decomposition costs in seconds for three sample matrices in HOBT In Table 3.2, we show the cpu times for decomposition and reordering for the three sample matrices. In Table 3.3, the cost of backsolving is presented. We should note here that in our context, while the decomposition and the reordering will be performed just once, the efficiency of the back-solver is important, for this step is repeated several thousand times. Size (v) Without reordering RCM reordering AMD reordering 59332 0.142 0.064 0.065 146872 0.4 0.16 0.15 211372 0.7 0.39 0.28 Table 3.3: Table shows the backsolving cost in seconds for a single sample RHS Notice from tables 3.2 and 3.3 that the AMD scheme of reordering is superior in our case. Henceforth, all studies are carried out using the AMD reordering. 3.3 Cholesky decomposition and the fill-in phenomenon Consider the cost of solving the system ˜ F ˜ Φ = ˜ W (3.3) using the Cholesky decomposition ˜ F = ˜ R T ˜ R for ˜ F, where ˜ F is the AMD-reordered v×v matrix, ˜ Φ and ˜ W are v×p reordered matrices, and ˜ R is the v×v Cholesky 32 factor of ˜ F. Due to the fill-in phenomenon that occurs even in well-structured matrices, the number of non-zeros in the Cholesky factor may increase compared to that in the original matrix, thus, contributing significantly to the computation cost. In this section, we investigate the effects of incomplete Cholesky decomposition on the back-solving cost. MATLAB’s implementation of the incomplete Cholesky decomposition (cholinc) with a drop tolerance (droptol) option [Saa96] was used for all the studies. This factor- ization is computed by performing the incomplete LU decomposition (here ˜ R T ˜ R) with the pivot threshold option set to 0 (which forces diagonal pivoting) and then scaling the rows of the incomplete upper triangular factor, ˜ R, by the square root of the diagonal entries in that column. Setting the droptol parameter to 0 results in a complete Cholesky decomposition. In Table 3.4, we show the dependence of the cpu time for forward and back-solving on the non-zero entries in the Cholesky factor ˜ R as a function of the drop- tol parameter for a sample matrix. Droptol nz for Cost for Cost of ˜ R decomposition back-substitution 0 6874570 12.02 0.57 10 −5 2974210 10.06 0.26 10 −4 1980942 5.73 0.18 10 −3 1132135 2.55 0.11 10 −2 549798 1.03 0.06 10 −1 173272 0.23 0.02 Table 3.4: Table shows the backsolving cost for a single sample RHS for the matrix with v = 59332. MATLAB’s cholin operator was used for the incomplete Cholesky decomposition. For the complete Cholesky decomposition, the same function was used with droptol = 0. We note the following from the table: 1. As the drop-tolerance value increases, the number of non-zero entries in the Cholesky factor decreases, and the backsolving cost decreases. 33 2. The cost for computing the incomplete Cholesky decomposition shows a simi- lar trend by decreasing with increase in droptol. However, this cost is of lesser importance than the cost of the backsolver (because p, the number of RHS, is very large). 3. The cost of back-substitution is proportional to the number of non-zerosnz in ˜ R. This study was important for enabling us to make a decision about using the incomplete Cholesky factor as a preconditioner for PCG. Notice, however, that with the increase of droptol value, the Cholesky factors are becoming increasingly approximate (incom- plete), see 3.4.2. 3.4 Solution schemes for forward models 3.4.1 Solution by sequential back-substitution A direct method to solve the problem in 3.1 is to use MATLAB’s ‘\’ operator for one RHS at a time. This operator reorders the sparse matrix, computes a complete Cholesky decomposition and solves the problem by forward and back-substitution. The draw- backs of using this scheme are obvious. The Cholesky decomposition will be carried out every single time for each RHS and the process will be very costly. An alternate method is to reorder the matrix, compute the Cholesky factorization only once, and use forward and back substitution for each individual RHS. It is absolutely essential to reorder the matrix before the factorization. In our experiments, the Cholesky factoriza- tion of the unordered matrix took almost 30 times more time than an AMD-reordered 59332×59332 matrix. The ratio increased to about 120 times for a 211372×211372 matrix. With AMD-reordering for the FEM-matrices, the profile of the non-zero ele- ments changes so that they are now closer to the main diagonal. This implies smaller 34 fill-in and hence, the Cholesky factor is computed much faster than using an unordered matrix. Table 3.5 shows the cost of using the ‘\’ operator sequentially versus using the AMD-reordered FEM matrix and it’s Cholesky decomposition for solving the system with just 50 RHS. The saving by using AMD-reordering is considerable and increases with the number of RHS, see Section 3.4.3. It should also be noted that MATLAB’s ‘\’ operator includes AMD reordering and computation of the complete Cholesky decom- position by default. Solving sequentially using ‘\’ amounts to reordering the matrix and computing the Cholesky factorization for each RHS. Cost of using\ for solving Cost for reordering and Cost for forward with 50 RHS repeatedly Cholesky decomposition and back-substitution 128.5 2.6 28.5 Table 3.5: Table shows the cost of repeated reordering, Cholesky decomposition, and solving versus pre-computation of the Cholesky factor and sequential backsolving for the matrix withv = 59332. Also, see section 3.4.3. 3.4.2 The incomplete Cholesky factor as a preconditioner for PCG The rate of convergence of the conjugate gradient method can be strongly affected by the choice of the preconditioner. The preconditioner is a positive definite matrix chosen to approximate the inverse of the F matrix. If the eigenvalues of the product of the preconditioner andF are tightly clustered compared to those ofF itself, convergence is achieved in relatively few iterations. The incomplete Cholesky factors, ˜ R can be used as a preconditioner and the system in (3.1) can be solved by the preconditioned conjugate gradient (PCG) iterations. The product of cost per iteration and the total number of iterations required for convergence constitutes the total cost of PCG. There is a trade-off while selecting the droptol parameter. On one hand, for each iteration of PCG, a system of equations with the matrices R T R needs to be solved, i.e. two backsolving operations are performed at each step (Table 3.6 shows the cpu 35 droptol CPU for back # of PCG CPU for PCG Decomposition Total substitution iterations time CPU 10 −5 0.26 3 1.06 10.06 11.12 10 −4 0.18 4 1.08 5.73 6.81 10 −3 0.11 7 1.1 2.55 3.65 10 −2 0.06 16 1.9 1.03 2.93 10 −1 0.02 51 5.14 0.23 5.37 Table 3.6: Table shows the cost for PCG to converge as a function of droptol for the incomplete Cholesky decomposition for the matrix of sizev = 59322 cost for a single RHS as a function of droptol). This cost is the largest contribution to the total cpu time. On the other hand, the incomplete Cholesky decomposition ˜ R T ˜ R leads to an approximate matrix ˜ F = ˜ R T ˜ R, and the disparity between the original matrix F and it’s approximation ˜ F (in terms of the number of non-zero entries and of spectral properties) increases as droptol increases. Consequently, the incomplete Cholesky factor ˜ R progressively become worse as a preconditioner for PCG as an effect of increase in droptol, and thus, the number of PCG iterations required for convergence increases. In Table 3.6, the cost of computing the incomplete Cholesky decomposition, the # of iterations for PCG (it), and cost for PCG iterations are presented as a function of droptol for the matrix of sizev = 59322. A comparison of Tables 3.4 and 3.6 shows that increase in droptol implies decrease in cpu for backsolving, but increase in the number of PCG iterations. As a result, we should expect to find an optimal value for droptol that minimizes the total cost. Figure 3.2 shows a plot of the time required for computing an incomplete Cholesky factoriza- tion and finding the solution by PCG. As indicated by the plot and the table, droptol = 10 −3 is close to optimal. We also note here that if the Cholesky factor is computed instead of the incomplete one with droptol = 0 (i.e. using MATLAB’s chol instead of cholinc) for the system of Table 3.6, the cost was 2 seconds as against 12 seconds for cholinc. 36 10 −5 10 −4 10 −3 10 −2 10 −1 10 −1 10 0 10 1 10 2 droptol CPU in sec Time required by PCG to converge Time taken for the Cholesky decomposition Figure 3.2: Plot of the cost for computing the incomplete Cholesky factorization and for PCG, as a function of droptol for the matrix of sizev = 59322 Also, note that the PCG-based method is better than sequential backsolving, see 3.4.3. 3.4.3 Solution using lumped RHS Since F is sparse and symmetric positive definite (SPD), the Cholesky factorization requires less than half the time of a general factorization. The Cholesky factorization isF =R T R, where R is upper triangular. The solutionΦ is computed by solving two triangular systems,RΦ =K andR T K =W). We lump together the individual RHS in sub-matrices, each of 500 columns. Subse- quently, matrix reordering and Cholesky factorization is used to compute 500 solutions, all at once. This procedure can be repeated until all the required solutions are computed. We also note that matrix W transfers the source strength from a 3D location onto the closest tetrahedron nodes. As a result, W is sparse with an average of 4 non-zero ele- ments per column. Thus, W can be precomputed and does not require a lot of memory. Figure 3.3 and table 3.7 show the cost of using lumped RHS versus PCG and of solving each RHS individually for the matrix withv = 211372. These results clearly 37 50 100 150 200 250 300 350 400 450 500 10 1 10 2 10 3 10 4 10 5 10 6 Sequential operation PCG Lumped RHS Figure 3.3: Cost of solving multiple RHS using sequential operations, PCG and using lumped RHS for the matrix of sizev = 211372 as a function of the number of RHS. # of Sequential PCG Lumped RHS Ratio RHS (a) (b) (a/b) 50 1420 308 33.5 9.2 100 2840 616 38.6 16.0 200 5680 1232 50.5 24.4 300 8520 1848 64.3 28.7 400 11360 2464 82.6 29.3 500 14200 3080 102.6 30.0 Table 3.7: Table shows the cpu for varying number of RHS using sequential operations, PCG and lumped RHS for the matrix of sizev = 211372. indicate that there is a big advantage in using lumped RHS. Remark: there is a trade-off between memory and the number of solutions computed at once. The number 500 for the RHS was experimentally determined to be suitable for the largest matrix (v = 211372). This means that processing more thanp = 500 RHS simultaneously would require using auxiliary memory, thus increasing the CPU cost significantly, while using RHS less than 500 would not use the computing power of the machine efficiently. Table 3.7 examines the effect of the number of RHS ‘p’ on the cost of computing the solution. As expected, the cost is proportional to the number of RHS. There is a reduction in the cost of more than 100 times when using the ‘\’ operator with lumped 38 RHS for p = 500. Consequently, there is a big advantage when using the method of lumped RHS with a high enoughp that is supported by the computer memory. For our machine (memory of 2 GB), this wasp = 500 forv = 211372. Table 3.8 examines the dependence of the size of the matricesv on the cost of com- puting the solution. For the two largest matrices (v = 146872 and v = 211372), we observe that the CPU cost grows as the square of the non-zeros elements (nz). As indicated by the last column in Table 3.8, the increase in saving is proportional to the increase in the size of the matrices (v). Size Non-zeros Sequential PCG Lumped RHS Ratio nz (a) (b) (a/b) 59332 782853 125 58.5 3.57 16.4 146872 2004522 570 124.5 12.32 10.1 211372 3156406 1420 308 33.5 9.19 Table 3.8: Table shows the CPU time for all three matrices under consideration and 50 RHS when using sequential operations, PCG and lumped RHS 3.5 Summary In this chapter, we presented fast solvers for the assembly of forward models in HOBT. We examined three methods for solving (3.2), (a) using MATLAB’s \ operator in a loop for sequential solving, (b) using PCG with the incomplete Cholesky decomposi- tion of the FEM matrix as a pre-conditioner, and (c) using matrix reordering, complete Cholesky decomposition and forward and back-substitutions with lumped RHS. Fol- lowing is a summary of conclusions from the investigations undertaken. We refer to the largest FEM matrix (v = 211372) asH: 1. The combined cost for matrix reordering and incomplete Cholesky decomposition (withdroptol = 10 −3 ) forH is about a 100 times less than the cost of computing 39 the incomplete Cholesky decomposition of unorderedH with the samedroptol, see table 3.2. 2. Increase in droptol corresponds to decrease in the backsolving cost. The cost of computing the incomplete Cholesky factor decreases with increase indroptol, see table 3.4. It should be noted here that a choice of droptol> 0 is equivalent to approximating the system and will introduce errors in the solution. 3. The cost of backsolving is proportional to the non-zeros (nz) in ˜ R, the incomplete Cholesky factor, see table 3.6. 4. If PCG is used for solving (3.2),droptol = 10 −3 is close to optimal. 5. The cost of solving (3.2) by the method of lumped RHS is about 30 times less than the cost for PCG forH. The cost for the existing method is about 4 times higher than that of PCG for matrixH. 6. The cost for solving usingH for 500 RHS is at least 30 times less than the cost for the existing method of using PCG. This work clearly shows that the cost for solving the HOBT forward problem is reduced by more than 30 times by using the proposed method of lumped RHS (see Section 3.4.3) than the existing method of using PCG 40 Chapter 4 Reconstruction methods for multispectral and hyperspectral OBT and OFT 4.1 Introduction Large-scale linear systems of the typeAq =b whereA∈R m×n ,q∈R n andb∈R m , m > n occur naturally in many disciplines of science and mathematics, e.g. structural engineering, electrical engineering, biomedical engineering, computer vision, atmo- sphere studies, economics etc. For OBT and OFT, the forward problem was modeled similarly in (2.22) and (2.23), where the forward models A hyp or A ill hyp denoted the mapping from any possible internal 3D location of the bioluminescent source or a fluo- rescent probeq bound to a region of interest (e.g. a tumor), onto every point on the ani- mal skin, where surface measurementsb hyp orb ill hyp can be made [HST + 94, CDB + 05]. Our inverse problem is to estimate the 3D distribution of the bioluminescent source (q) inside the boundarydΩ, when the photon flux on the boundary and the optical scatter and absorption properties are known. For achieving this, we are required to invert the systemAq =b to findq (we have dropped the suffixes for convenience). For OBT and OFT, the forward operatorA is large and ill-posed and consequently, iterative techniques are used for inversion. Iterative methods use an initial guessq 0 and find a sequence of approximationsq k , that converge to a solutionq. This solutionq is not unique for a 41 rectangular system as is the case for OBT/OFT, thus, some a priori assumptions about the unknownq need to be incorporated while solving the system. In this chapter, we describe our formulation of the inversion as an optimization problem with objective function f(q). Depending on the type of the desired recon- structed distribution (smooth or sparse), the regularization penalty function is chosen. The update equations for minimizing f(q) using the Polak-Ribie` re form of the non- linear preconditioned conjugate gradient (PCG) method are shown. These equations can be directly extended to OS-PCG type methods where PCG is used on subsets of the data at each subiteration. For reducing the computational cost, dimensionality reduction ideas using spatial basis functions (3D Discrete Cosine Transform, and 3D Daubechies D-4 wavelets) are presented. 4.2 Mathematical formulation of the objective function The image reconstruction is formulated as an unconstrained minimization of the objec- tive functionf(q), where penalty functions are used to favor reconstruction of solutions close to the feasible set. The functionf(q) combines a measure of the fit to the data, a regularizing function, and a non-negativity penalty function. 4.2.1 Measure of the fit to the data As a measure of fit to the data, we used a simple squared error ψ(q) =kb−Aqk 2 2 . (4.1) 42 A more accurate model would take into account the spatially variant noise that arises in optical systems due to variations in the incident photon flux and properties of the CCD focal plane [MOW + 03]. 4.2.2 Regularization penalty If the forward modelA is rank-deficient (as is the case in OBT and OFT), there are an infinite number of reconstructions that minimizeψ(q). This is because adding any vector in the nullspace of A to the reconstruction will not affect the ψ(q). Instead, we may regard the data as defining a feasible and consistent set of reconstructions for whichψ(q) is less than some tolerance value depending on the noise level. One way of representing our bias to a particular class of solutions is by introducing additional information about the imageq through an image-dependent norm. Here we use aL 2 - norm or anL 1 -norm regularizing function to control the instabilities in the image: Ω(q) =kqk 2 2 = X i |q i | 2 . (4.2) or Ω(q) =kqk 1 = X i |q i |. (4.3) 4.2.3 Positivity Penalty Since bioluminescent or fluorescent source distributions are non-negative, the quality of the reconstructed images should be enhanced by constraining the estimated distributions to be nonnegative. Rather than explicitly constraining the search to the non-negative orthant, we instead used a penalty function φ(q) = X j q 2 j u(−q j ) (4.4) 43 where u(.) is the unit step function andj is the pixel index in image q [Lue84]. This type of penalty function has been successfully used by Mumcuoglu et al [MLCZ94] in the reconstruction of tomographic images using preconditioned conjugate gradient methods. 4.2.4 The objective function The combined objective function is: f(q) =kb−Aqk 2 2 +αkqk 2 2 +μ X i q 2 i u(−q i ) (4.5) whereα is the regularization parameter that controls the trade-off between agreement to the data and the smoothness of the solution. We selectedα using the L-curve technique [Han92]. An L-shaped curve is obtained from the plot of the logarithm ofΩ(q) against the logarithm ofψ(q) by varyingα. On the vertical part of the L-curve, there are still components in the image arising from small singular values and inadequate regulariza- tion. On the horizontal part of the curve, theΩ(q) does not change much but the residual norm increases sharply. The chosen value ofα is close to the knee of the curve. The parameterμ represents the weight attached to the non-negativity penalty. In theory, the solution to the constrained problem is given by the minimizer of f(q) in the limit as μ → ∞ [Lue84]. In practice, increasingμ can lead to ill-conditioning and poor con- vergence of the conjugate gradient algorithm. We found that a valueμ = 0.01×q (i) max , whereq (i) max is the maximum intensity in the image at iteration i, is sufficient to suppress negative values without adversely affecting convergence rates. 44 4.2.5 The optimization algorithms and preconditioner We used a preconditioned conjugate gradient (PCG) algorithm to minimizef(q). Since the penalty function is non-quadratic, we used a Newton-Raphson line search to deter- mine the optimal step size at each iteration [Lue84, She94]. The outline of the precon- ditioned PCG algorithm, withC as the preconditioner is as follows: d (0) =−g (0) =−f ′ (q (0) ) Findθ (i) that minimizesf(q (i) +θ (i) d (i) ) by line search q (i+1) =q (i) +θ (i) d (i) r (i) =C(g (i) ) β (i−1) = max ( r T i (g (i) −g (i−1) ) r T (i) g (i−1) ,0 ) , d (i) =r (i) +β (i−1) d (i) . The conditionβ = max(β PR ,0) whereβ PR = r T i+1 (g (i+1) −g (i) ) /(r T (i) g (i) ) is equiva- lent to restarting PCG ifβ PR < 0, and is required for convergence of the PCG algorithm [She94]. The rate of convergence of the conjugate gradient method can be strongly affected by the choice of the preconditioner. The preconditioner is a positive definite matrix chosen to approximate the inverse of the Hessian matrix. If the eigenvalues of the product of the preconditioner and the Hessian are tightly clustered compared to the Hessian itself, convergence is achieved in relatively few iterations. Kaufman [Kau87] noted that the EM maximum likelihood algorithm for image reconstruction of emission tomography data, in which the data follow a Poisson distribution, can be recast as a gradient descent method that uses the heuristic diagonal preconditioner of the form: C(q) = diag ( q i P n A hyp ni ) (4.6) 45 and proposed speeding up computation of a maximum likelihood solution by replacing the EM algorithm with a preconditioned conjugate gradient method using (4.6) as a preconditioner. Mumcuoglu et al. [MLCZ94] successfully applied this idea to penalized maximum likelihood with the same non-negativity penalty function as described above. We have adopted the same preconditioner here and found that it significantly speeds up convergence rates compared to a standard conjugate gradient search. Investigation of fast inverse methods is currently underway. Methods that are being evaluated are ordered subset methods, coordinate descent methods, incremental gradient methods and fast simultaneous forward and inverse solvers. 4.3 Using spatial basis functions for dimensionality reduction Consider the objective function after dropping the positivity penalty (an analogous method may be developed with the positivity penalty such thatq>0): f(q) =kb−Aqk 2 2 +αkqk 2 2 (4.7) Let B denote an orthonormal basis matrix [VK95] i.e. B = h e 0 e 1 ... e k i , B ∈ R n×k and BB T =I n . Using the same convention as in (4.7), we can rewrite the minimization as f(q) = b−A(BB T )q 2 2 +α (BB T )q 2 2 (4.8) LetP =AB,P∈R m×k andy =B T q,y∈R k . Assume a sparse representation such that most of the energy is contained only in a few basis coefficients, i.e. k<< n, the 46 minimization of the function in (4.8) is equivalent to solving a much smaller problem given by the following equations ˆ y = argmin y kb−Pyk 2 2 +αkByk 2 2 (4.9) ˆ q = Bˆ y (4.10) One way to solve the above system is by taking the gradient of the objective func- tion with respect to y and setting it to zero. Consequently, the above system becomes equivalent to solving the normal equations (P T P+αI k )y =P T b (4.11) Note that the bracketed term on the LHS of (4.11) belongs toR k×k , and that the matrix P is much smaller than A. The PCG algorithm, as proposed earlier, can now achieve convergence in less number of steps. 4.3.1 Choice of basis functions Any orthonormal basis should work in the above framework. However, to realize a substantial computational saving, it is important that only a few bases are required to accurately model the data up to a desirable tolerance value. We propose two families of basis functions and demonstrate their usage. The Discrete Cosine Transform (DCT) Emission images due to their smooth nature may be modeled as a highly correlated ran- dom fields such that they can be approximated using a separable model (with respect to the cardinal axes) in which each component is a 1st order Markov chain. It is known that 47 (a) (b) (c) Figure 4.1: Coefficients of the 3D DCT and the 3D Daubechies wavelet; (a) A sample fluorescence image plotted on coronal sections of a cube-shaped tissue phantom; (b) The 3D DCT coefficients; (c)The 3D DaubechiesD 4 coefficients. the DCT is equivalent to the Karhuenen-Lo` eve transform of a first-order Markov process [Jai88]. Therefore, most of the energy is compacted in the first few coefficients of the DCT. This property makes the DCT a potentially desirable choice for this application. The forward 3D DCT is given by S(w,v,u)=β 3D (w,v,u) N F −1 X z=0 N R −1 X y=0 N C −1 X x=0 s(z,y,x)cos(t 1 )cos(t 2 )cos(t 3 ) (4.12) where t 1 = (2x+1)uπ 2N C ,t 2 = (2y+1)vπ 2N R ,t 1 = (2z+1)wπ 2N F (4.13) β 3D (w,v,u)= r 2 N F 2 N R 2 N C C(w)C(v)C(u) (4.14) and C(k) =        1 √ 2 k=0 1 otherwise (4.15) Here,s(x,y,z) is the pixel value at location (x,y,z) andS(w,v,u) is a DCT coeffi- cient in the correspondingN F ×N R ×N C DCT cube. Analytical equations are available 48 for 3D DCTs and the basis matrix can be generated rapidly. Note that we do not require the full basis, since only a few coefficients will carry most of the energy. The Daubechies waveletD 4 basis An alternative choice for basis functions is to use a wavelet representation. The Daubecheis waveletD 4 is one choice. The coefficients of the Daubechies waveletD 4 [Dau92] are given by h 0 = 1+ √ 3 4 √ 2 (4.16) h 1 = 3+ √ 3 4 √ 2 (4.17) h 2 = 1− √ 3 4 √ 2 (4.18) h 3 = 3− √ 3 4 √ 2 (4.19) The scaling vectors = (h 0 ,h 1 ,h 2 ,h 3 ) and the wavelet vectorw = (h 3 ,−h 2 ,h 1 ,−h 0 ) are orthogonal. The 1DD 4 basis can be easily assembled using the scaling and wavelet vectors, with shifts of two. Since the basis is separable, the 3D orthonormal basis can be obtained by kronecker products of the 1D bases. 4.3.2 Simulation setup and methods The method is demonstrated here for OFT alone, but can be extended to OBT if no external laser illumination is used. We simulated a cube-shaped tissue geometry, with detectors placed on one side of the cube and a laser for illumination on the opposite side from the detectors (sometimes called an epifluorescence setup). The tissue slab of dimension 16 mm× 16 mm× 16 mm was divided into 4096 voxels, each 1 mm× 49 (a) (b) (c) Figure 4.2: Noiseless data reconstructed using variable number of DCT basis functions; (a) Using 512 basis functions; (b) Using 1728 basis functions; (c) Using the complete basis 1 mm× 1 mm. The optical fluence at 256 nodes was computed using our fast imple- mentation of the analytic monochromatic forward model solver A [HST + 94, CDS + 03]. The laser was stepped 4 times. Each laser location corresponded to one forward model. Multispectral data was assumed to be collected in four 20 nm bins between 600 nm and 680 nm. Consequently, one could view the system dimension as 256 detectors× 4096 sources× 4 illuminations× 4 spectra. The model A was projected onto a variable number of basis functions of the two families of bases. A fluorescent point source was simulated at a location 8 mm below the surface. The 3D location of the source was reconstructed using the PCG technique and the objective function derived in (4.11). The method presented here can be extended further by using a dictionary of basis functions and orthogonal matching pursuit [MZ93] methods to choose the optimal basis. The machine used to perform simulation studies was an AMD Opteron 250 64 bit, dual CPU machine. We evaluated the performance of the two families of basis functions by two metrics; CPU time (sec), and Full Width at Half Maximum (FWHM) (mm) resolution. 50 (a) (b) (c) Figure 4.3: Noiseless data reconstructed using variable number of DaubechiesD 4 basis functions; (a) Using 512 basis functions; (b) Using 1728 basis functions; (c) Using the model as is 4.3.3 Results and Discussion Reduction in computational cost We have tabulated below, the results obtained by reconstruction of a point source using the DCT basis functions and the Daubechies wavelet D 4 basis. Substantial reduction in cost is achieved by using 512 or 1728 basis functions at the expense of decrease in resolution (FWHM). No. of basis functions 512 1728 Full(4096) CPU Time (sec) 6.53 250.82 1668 FWHM (mm) 3.2 2.41 1.9 Table 4.1: The table shows values that were obtained with noiseless data and using the DCT basis functions. No. of basis functions 512 1728 Full(4096) CPU Time (sec) 9.615 320.51 1668 FWHM (mm) 2.42 2.24 1.9 Table 4.2: The table shows values that were obtained with noiseless data and using the DaubechiesD 4 basis functions. 51 Denoising The optical emission signal has spatial structure, whereas noise is random. Thus, the signal may be modeled very well using the basis functions, while the noise will tend to be suppressed in the reduced solution space. As a result, there is substantial gain in reconstructing noisy data using basis functions. The method proposed above has a natural property of denoising at the expense of obtaining smoother solutions. Improved conditioning While A is severely ill-conditioned, the projected model P has better conditioning since the bases act as regularizing functions. With 512 Daubechies basis functions, the con- dition number of P was in the order of 10 6 against that of A which was of the order of 10 11 . Thus, the solution using the above method is expected to be numerically much more stable compared to the direct solution. Because of this property and using a small number of basis functions, we may be able to directly invert the system. Computational advantage using the multiresolution/multigrid idea In a separate simulation, we projected the data on a small basis set first, and recon- structed the image on that coarse subspace. Subsequently, we projected the recon- structed image onto a larger set of basis functions and used this as a starting vector for the higher dimensional or finer subspace. Tabulated are the times required for direct computation and via using a two-step multi-resolution approach. Clearly, using multi- resolution reduces the cost by almost half for this example. 52 No. of basis functions 1728 Full(4096) Time required for individual computation (sec) 320.51 1668 Total time required using the multiresolution method (sec) 320.51 796 Table 4.3: The table shows values that were obtained with noiseless data and using the DaubechiesD 4 basis functions. 4.4 Summary The optimization scheme proposed in this chapter has been used extensively by us to generate the results for several tissue geometries as will be demonstrated in Chapter 6. The positively constrained PCG method has desirable properties for solving ill-posed problems. Work is currently in progress for comparing the convergence of PCG against the incremental gradient methods. The method of projection onto spatial basis functions offers dramatic reduction in computational cost, inherent denoising properties and better conditioning at the expense of a small loss in resolution. 53 Chapter 5 Experimental setups and the imaging process 5.1 Introduction Commercially available optical systems for bioluminescence and/or fluorescence imag- ing are largely limited to systems that collect data only from a single view. While these systems give satisfactory reconstruction results for source locations at about 1 cm deep, better depth resolution is required, particularly in the central region of the animal. The inversion problem using only a single view is highly ill-posed compared to when multi- ple views are used. We have designed a prototype mirror setup that allows for collection of data from four orthogonal views (top, two sides and bottom), and consequently fully 3D OBT or OFT. The system is designed to work with any 2D planar imaging system. In this chapter, we compare reconstructions from using simulated data from a single view versus using tomographic data acquisition, and clearly demonstrate that the latter is necessary for accurate localization of deep sources. Next, we describe the design of the optical setup that allows for tomographic viewing, and describe the fabrication of tissue-simulating phantoms to be used for experiments described in Chapter 6. Further, we describe an OBT/OFT imaging process chain that incorporates methods presented in Chapters 2, 3 and 4. 54 5.2 Comparison of single-view data acquisition versus using four orthogonal views 0 2000 4000 6000 8000 10000 12000 10 −14 10 −12 10 −10 10 −8 10 −6 10 −4 10 −2 10 0 # of sources Singular values Detectors on one side Detectors on two sides (a) Figure 5.1: The conditioning of the forward model and localization error. (a) The singu- lar value spectra of the forward models. (b) The localization error if data from multiple views is used for reconstruction. (c) The localization error when using data from a single-view. The sources were placed along the line of intersection of the x = 20 and y = 20 planes. The detectors for case (b) are on the z = 40, z = 0, x = 40 and x = 0 planes, whereas those for case (c) are on the z = 40 plane. A 40 mm× 40 mm× 40 mm homogeneous tissue cube was used for the simulation study. The source grid-spacing was 1.75 mm. For the single-view experiment, detectors were assumed to be on the z = 40 plane with grid spacing of 0.75 mm (corresponding to 55 the medium binning setting of Xenogen’s IVIS 200 system). For data acquisition using four orthogonal views, detectors were placed on the z = 40, z = 0, x = 40 and x = 0 planes with a grid-spacing of 1.5 mm (corresponding to the low resolution binning of the IVIS 200 system). The total number of detectors in each case were 2916. The tissue cube was assumed to have skin optical properties and multispectral data were simulated. Images were reconstructed using the truncated pseudo-inverse technique owing to the large number of inverses used. The singular value decomposition plots in Figure 5.1a clearly indicate that the model corresponding to tomographic viewing is far better conditioned than using single-view detection, and consequently, the inversion must be more stable than for a model using a single view. Figures 5.1b and c demonstrate that sources throughout the tissue volume can be reconstructed with sufficient accuracy using tomographic data, while data using only a single view is incapable of reconstructing sources deeper than 10 mm from the detector plane. 5.3 The mirror setup for tomographic data acquisition In order to perform fully 3D reconstruction of bioluminescent sources, we needed to gather data over the whole surface of the animal. A photograph and schematic of the mirror setup that we have designed for this task are shown in Figures 5.2 and 5.3. The setup consists of a 13 cm×13 cm×1 mm non-fluorescent transparent plexiglass animal bed that may be placed on the subject platform of a 2D imaging system. The bed is elevated by 8 cm on two sides by 4 mm thick plexiglass planks, leaving the other two sides open. A pair of L-shaped benches is designed to hold 7 cm×4.5 cm mirrors at a 45 0 angle with respect to the bed. This pair is mounted on the animal bed and is used to acquire two opposing side views of the subject. The animal is placed along the centerline 56 of the bed. A second pair of L-shaped benches is fabricated to hold mirrors of size 7 cm×8.5 cm, again at a 45 0 angle. The mirrors are placed on the subject platform such that one is directly below the plexiglass bed, and the other is displaced to the side. This mirror pair acquires the bottom view of the animal. All results presented here used the IVIS 200 2D imaging system (Xenongen, Alameda, CA). The subject platform of the IVIS 200 system is capable of vertical translation thus enabling us to focus individually on the top, side and bottom views. Figure 5.2: A photograph of the mirror setup designed for acquisition of tomographic data. The setup has no moving parts and can be accommodated by any of the 2D Xenogen IVIS bioluminescence imaging systems. This setup requires that the imaging system have a field of view of at least 17 cm×17 cm in order to capture the four orthogonal views. One of the standard platform posi- tions of the IVIS 200 (level D, demagnification x7.5) produces a field of view of 19.5 cm×19.5 cm and is well suited to the mirror setup. The Xenogen Living Image software corrects for cosmic noise and lens distortion. The lens system of the IVIS 200 allows f- stops of f/1, f/2, f/4 and f/8. Four orthogonal multispectral views may be acquired using the mirror setup and the IVIS 200’s back-thinned, back-illuminated Grade 1 26 mm×26 57 mm CCD camera without the need to move the phantom or mouse under study. Multi- spectral data for phantoms and experiments were acquired using a filter wheel mounted with six filters (response peaks at 560 nm, 580 nm, 600 nm, 620 nm, 640 nm and 660 nm respectively; FWHM = 20 nm) that are matched to the spectrum of firefly luciferase, which we used in the in vivo study reported in Chapter 5. Figure 5.3: (a) Schematic of the mirror setup for acquisition of tomographic data; (b) A sample image obtained from the CCD camera using the mirror setup. The bottom and side-views are scaled relative to the top view because of differing optical path lengths; this is accounted for in our forward model. The focusing of the bottom view can be cor- rected by moving the subject platform vertically to the same focal plane as the top-view. Since the setup has no moving parts, the relative distances remain constant enabling a fixed correction for perspective. 5.4 Phantom preparation and materials for phantom studies Our objectives in conducting the phantom studies were: 1. to have a simple homogeneous model with known optical properties and geome- try; 58 2. to validate the accuracy of localization of the point source in a real experimental setup; and 3. to perform resolution analysis at various depths in the volume. Rectangular block-shaped phantoms (24 mm×24 mm×14 mm) were prepared from material that mimics the average optical properties of biological tissue [CPT + 97]. Every 100 ml of the material solution contained a mixture of 95 ml of distilled water as a sol- vent, 5 ml of 20% Intralipid (I-141, Sigma-Aldrich, USA) as a scattering medium, 1μl of India ink (Rotring, Germany) as an absorbing medium, and 1.05 g of highly purified Agar powder (A7049, Sigma-Aldrich, USA) as a solidifying agent. The material cured within an hour and was easily cast into the desired size and shape. As a bioluminescent point source, we used the cleaved tip of a multimodal optical fiber (Siecor, USA) with a spot diameter of 100 microns coupled to an orange LED (#276-273, RadioShack, USA). The orange LED was chosen because its wavelength spectrum is similar to that of firefly luciferase. A glass capillary tube was used to guide and localize the fiber within the phantom. Inhomogeneous phantoms were prepared by varying the concentrations of Intralipid and India ink. As a preliminary experiment, we prepared rectangular inhomogeneous phantoms (dimension: 4 cm x 4 cm x 4 cm) having two-layered (skull enclosed in skin) and three-layered (one skull layer embedded in the skin) geometries. The skull absorption and reduced scattering coefficients deviated from those of skin by +30% and +20% respectively. Figure 5.4 shows the molds and phantoms that were prepared using the method described here. For fluorescence studies, we inserted fluorescence probes into the phantoms by means of a needle. After the experiment, the phantoms were sliced open and the exact location of the probe was measured by vernier calipers. 59 Figure 5.4: Homogeneous mouse and slab phantoms and inhomogeneous slab phantoms designed and fabricated for experiments described below. The phantoms were prepared following the recipe from Cubeddu [CPT + 97]. 5.5 An unified OBT and OFT process chain From an animal laboratory point of view, it is important to have an imaging process chain for OBT and OFT, that requires minimum human interaction and maximum automa- tion. The process chain we implemented is as follows. The animal was anesthetized and placed on a plexiglas bed so as to keep the same coordinate system throughout the experiment. The option for structured light was turned on, images were captured and the surface of the animal was constructed. A mouse atlas was warped to the animal sur- face. Optical properties were assigned to the internal organs and the computation of the forward model began. Concurrently, multispectral or hyperspectral data was acquired using the mirror setup described here. Every set of spectral images were interlaced with an image without using any filter. This set of unfiltered images was used to estimate the time activity curve of the underlying source. The collected multispectral or hyper- spectral data was registered to the mouse surface, and calibration factors and correction 60 were applied. After completion of the computation of the forward model, the inverse was computed using the methods presented in this thesis. The results are shown in the warped atlas coordinate system so that anatomical regions can be identified. We automated and integrated individual methods presented in Chapters 2, 3 and 4 into an imaging process chain as a part of this thesis work. 5.6 Summary Tomographic data acquisition is necessary for accurate localization of source distribu- tions in small animals. The mirror setup facilitates this, and may be used with any standard 2D imaging system. The method of phantom preparation described allows flexibility, since the phantoms can be cast in desired shapes and sizes, and with desired optical properties. A unified OBT/OFT process chain was developed and implemented. 61 Chapter 6 Simulation and experimental studies 6.1 Introduction In order to validate approaches proposed in this thesis for the three essential components for 3D multispectral and hyperspectral OBT or OFT, viz, the assembly of the multispec- tral and hyperspectral forward models, the 3D tomographic reconstruction technique, and the optical setup that enables tomographic data acquisition, we conducted several simulation studies, phantom experiments and animal studies. In this chapter, we present results from the following studies: Simulation studies 1. Simulation study in a mouse atlas using hyperspectral data where we compare bioluminescent achromatic and hyperspectral models based on SVD, and present resolution studies in OBT and OFT [CDCL05]. 2. Investigation of the impact of skull thickness and skull optical properties on recon- struction of brain tumors [DCL05]. 3. A multi-modality PET-OBT tumor study for evaluating localization accuracy in an inhomogeneous mouse phantom derived from a mouse atlas [SCS + 02]. 4. Simultaneous reconstruction of sources with different emission spectra and appli- cations using quantum dots. 62 5. A feasibility study exploring applications of OFT to human breast cancer imaging. Experimental studies 1. Phantom studies for analysis of localization accuracy and resolution for biolumi- nescent sources using the IVIS system [CDCL05]. 2. In-vivo study in a mouse with an implanted brain tumor [CDB + 05]. 6.2 Simulation studies in a mouse atlas using hyper- spectral data for investigating resolution limits and localization error An anatomical mouse atlas based on x-ray CT and cryosection images [SCS + 02] was used to define the surface of a 3D mouse volume which was then tessellated using a tetrahedral mesh. 16177 source locations were defined on a uniform grid in the thorax and abdomen with a 0.8 mm grid spacing. For the purposes of this study the mouse was assumed to have homogeneous optical scattering and absorption properties. While in practice, data are acquired in the focal plane of the CCD camera, for convenience we assumed here that measurements are collected directly on the surface of the ani- mal. The photon flux at each surface node was computed for each of 100 frequency bins uniformly spaced in wavelength from 600 nm to 800 nm. The bioluminescent source was assumed to have the frequency spectrum of Vybrant DiD R dye. The ana- lytic approximation was used to generate each monochromatic forward model. The monochromatic forward models were projected onto three depth-dependent principal spectra for dimensionality reduction. The achromatic model used for comparison was 63 generated as described in (2.17). Noiseless achromatic and hyperspectral data were sim- ulated. To facilitate the large number of inverses computed in this study, images were reconstructed using a truncated pseudo-inverse of the system matrixA in (2.22), where the truncation parameter was set to retain the singular vectors corresponding to singular values with magnitude greater than 10 −6 of the maximum singular value. In a limited number of comparisons to the PCG-based minimization of the regularized function (4.5) without the non-negativity constraint (again, so as to not over-estimate resolution), we found very similar performance for these noiseless data. An SVD-based analysis of the bioluminescent achromatic versus hyperspectral models The forward modelA∈R m×n can be decomposed using the SVD as A = min(n,m) X i=1 σ i u i v T i =UΣV T (6.1) whereU andV are orthogonal matrices andΣ is a diagonal matrix of singular values. The right-singular vectors (columns ofV) and left-singular vectors (columns ofU) may be viewed as a natural basis forA and represent respectively an orthogonal decomposi- tion of the source and data spaces. For comparison of the achromatic and hyperspectral models, we used three times as many surface nodes or measurements in the achromatic model as in the hyperspectral model to account for the additional hyperspectral mea- surements arising from the projection of the spectral data onto the first three spectral basis vectors. Figure 6.1a shows the singular value spectra of these achromatic and hyperspectral models. Typically, one would expect a superior imaging system to show a slower rate of drop off in the singular values. This is not the case here: if anything the spectra would lead to the erroneous conclusion that the achromatic system is preferable. 64 0 1000 2000 3000 4000 5000 6000 7000 8000 10 −12 10 −10 10 −8 10 −6 10 −4 10 −2 10 0 # of detectors Singular values Achromatic Hyperspectral (a) (b) (c) Figure 6.1: SVD analysis of the achromatic versus hyperspectral forward models. The color scale here is such that black indicates no spatial contribution and white indicates maximum contribution: (a) SVD spectra of the achromatic and hyperspectral biolumi- nescent forward models with matched detectors; (b) Overlay of the absolute value of the right singular vector corresponding to the fourth largest singular value of the achromatic forward model on the horizontal sections through the 3D mouse volume; (c) Overlay of the absolute value of the right singular vector corresponding to the fourth largest sin- gular value of the hyperspectral forward model on the same sections. (b) and (c) were gamma-corrected for better visibility. However, if one compares the 3D spatial distribution of the right singular vectors which span the source space of the larger singular values for both achromatic and hyperspec- tral data, it is evident that the hyperspectral data do indeed lead to a superior ability to detect deep sources. We have illustrated this observation in Figure 6.1b and c, where the absolute value of the fourth source basis vectors are shown overlaid onto horizontal sections through the mouse. These basis functions show the energy for the achromatic model to be concentrated in the superficial regions of the animal, which explains why the reconstructions in figures 6.2 and 6.14 have the point source reconstructed near the 65 boundary of the phantom. In contrast, the hyperspectral basis functions show substantial power in these vectors running through the deepest parts of the animal. Resolution limits with hyperspectral bioluminescent data The purpose of this study was to examine the effect on resolution of increasing the spatial sampling density at which photon flux is measured. We also used this study to determine an upper bound on achievable resolution for a realistic mouse geometry in the case where the optical properties are homogeneous and known. While this limit cannot be achieved practically, it is useful to compare these ideal resolutions to those measured in the phantom study presented in Section 6.7. The monochromatic, achromatic and the hyperspectral models were computed for detector sampling rates corresponding to 30, 146, 301, 628, 995, 1624, 2645, 3319, 4026, 5146 and 6129 surface nodes, and data were simulated for a point source with an emission spectrum of Vybrant DiD R dye located deep in the abdomen, as shown in Figure 6.2a. Average values of optical properties (μ a = 0.05 mm −1 ,μ ′ s = 2.3 mm −1 ,η = 1.4) were used to generate the monochromatic models. Figures 6.2b, c and d show that the monochromatic or the achromatic model, using same number of detectors as the hyperspectral model, reconstructs source distri- butions that are less focal and extend to the surface. The maxima of the reconstructed distribution for monochromatic data and achromatic data were at distances of 8.8 mm and 2.7 mm from the true 3D location of the simulated point source respectively. The reconstruction results using hyperspectral data show accurate localization of the 3D dis- tribution of the point-source. Figure 6.3 shows the effect on resolution of increasing the spatial sampling rate for hyperspectral data. Results indicated that sub-millimeter resolution may be achieved if more than 4600 surface nodes are used. We also observed that increasing the sampling rate beyond 5500 does not impact the resolution, indicating 66 Figure 6.2: Reconstruction of a deep point-like source in an atlas-mouse geometry using hyperspectral data. The color scale used here goes from light (no intensity) to dark (max- imum intensity). (a) Horizontal sections traversing the torso of the mouse and indicating the location of a point-like bioluminescent source used for simulation. (b) The recon- structed source distribution using monochromatic data collected at 2645 detector loca- tions and average optical properties over the 600-800 nm window.(c) The reconstructed source distribution using achromatic data collected at 2645 detector locations. (d) The reconstructed source distribution using hyperspectral data collected at 2645 detector locations and 100 spectral bins (corresponding to FWHM of 1.5 mm). The reconstruc- tion from hyperspectral data is able to accurately localize the source distribution, but reconstruction from monochromatic and achromatic data indicate a broadly distributed source. that at this point sampling requirements corresponding to the resolution limit have effec- tively been met. A similar study for the monochromatic and achromatic models using the same number of surface nodes as were used for the hyperspectral model revealed that the source distribution is not localized accurately. 67 0 1000 2000 3000 4000 5000 6000 7000 0.5 1 1.5 2 2.5 3 3.5 # of actual detectors FWHM in mm Figure 6.3: A plot of the FWMH resolution obtained from hyperspectral data for a reconstructed point source located deep in the abdomen as a function of the spatial sampling on the surface of the 3D animal volume. Using external illumination with monochromatic fluorescent data For the mouse atlas geometry, we computed the monochromatic forward models cor- responding to a wavelength of 700 nm using the optical properties at this wavelength. A point-like external illumination source (a laser with emission wavelength of 640 nm) was assumed to be placed at 4, 5, 6 7, and 8 laser positions on a concentric ring around the source location of interest. The source spectrum was assumed to be that of the Vybrant DiD dye. Figure 6.4 shows the change in the FWHM of a point source placed in the bulk of the mouse as a response to changing the spatial sampling of the detec- tors. Clearly, using external illumination reduces the illposedness of the HOFT inverse problem and adds independent information. 68 Figure 6.4: A plot of the FWMH resolution obtained from monochromatic data using multiple illumination patterns for a reconstructed point source located deep in the abdomen as a function of the spatial sampling on the surface of the 3D animal volume. 6.3 An investigation of the impact of skull thickness and skull optical properties on reconstructions of brain tumors using OBT In order to analyze the localization error due to mis-estimation of the skull optical prop- erties and the skull thickness, we performed several simulations using a mouse head model and a sphere. Both geometries contained a closed skull layer inside an otherwise homogeneous medium, for which we assumed skin optical properties. For each of these geometries, the skull optical properties and thickness were varied and multispectral data for wavelengths ranging from 600 nm to 700 nm (steps of 20 nm) were assumed to be collected for a point source inside the skull layer. The multispectral bioluminescent forward model was computed using the true optical properties just once. The simu- lated data reflected the perturbations in optical properties of the skull and the changes in thickness in individual simulations. Images were reconstructed using the regularized 69 600 620 640 660 680 700 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 Wavelength (nm) Absorption coefficient (1/mm) Skin Skull (no perturbation) Skull (+25%) Skull (−25%) 600 620 640 660 680 700 1.5 2 2.5 3 3.5 4 Wavelength (nm) Reduced scattering coefficient (1/mm) Skin Skull (no perturbation) Skull (+25%) Skull (−25%) (a) (b) Figure 6.5: Wavelength dependent optical properties of the skull and the skin. Fig- ures (a) and (b) shows the absorption coefficient and the reduced scattering coefficients respectively. pseudo-inverse owing to large number of inverses required for this study. The wave- length dependent optical properties for skin were taken from [CPW90] and for skull were from [FHED93] as shown in Figure 6.5. The skull properties were varied from the original published value in -25%, -10%, 10% and 25% steps. Sphere models Two spheres (sphere 0 and sphere 1), each with radius 11.8 mm, were used for the simulation. The skull layer thickness in sphere 0 was 0.76 mm and in sphere 1 was 1.16 mm. Sphere 0 was tessellated into 625989 tetrahedral elements and 116415 nodes. Sphere 1 was tessellated into 666158 tetrahedrons and 123283 nodes. A regular source grid of 7360 source points with a grid distance if 0.84 mm was placed inside the volume enclosed by the skull layer. 2989 surface nodes that were equally spread across the surface acted as detectors. A point source was placed at coordinates (0.85 0.85 5.41) mm with the center of the sphere at (0 0 0) mm. The source was placed off-center, because a source at the center would produce the same pattern at the surface detectors independent of the material properties due to symmetry. For the reconstruction of the 70 point source, we used different values for the regularization parameter, ranging from 10 −2 to10 −9 for studying the effect of regularization on the solution. Mouse head model Two mouse head geometries (head 0 and head 1) were used, which had an average skull thickness of 0.48± 0.074 mm and 0.67± 0.1 mm. The geometries were generated from a segmented CT image of a 26g athymic nude (nu/nu) mouse. The dimensions of the enclosing cuboid were 27.58 x 18.76 x 31.08 mm. The skull compartment was generated by dilating the outer hull of the segmented brain. Head 0 was tessellated into 274211 tetrahedrons and 59335 nodes. Head 1 was tessellated into 272454 tetrahedrons and 59075 nodes. A regular source space of 9296 nodes with a grid distance if 0.8 mm was used in both models with 2516 surface detectors. A single point source located to the frontal part of the brain was simulated inside the skull cavity. The regularization parameters were varied as described in the sphere model. Results In the absence of noise, the regularized pseudo-inverse gives accurate reconstruction results (see Figure 6.6a), if the true forward model is used for data generation and source reconstruction and the regularization parameterα is small (10 −9 ). If the model is per- turbed and a small regularization is used, the solution becomes distorted. If the regu- larization parameter is increased, the resolution drops. For Figure 6.6, the localization error was defined as the spatial distance between the maximum of the smoothed recon- structed solution and the true source location. Simulations in both geometries show significant changes in reconstruction results, when the optical properties of the skull layer are changed. In the spherical models, this leads to a mislocalization of the source of up to 10% of the size of the total geometry. In general, if the optical properties of 71 1e−2 1e−3 1e−4 1e−5 1e−6 1e−7 1e−8 1e−9 0 1 2 3 4 5 6 7 Regularization parameter Localization error in mm Minus 25% Minus 10% As Is Plus 10% Plus 25% 1e−2 1e−3 1e−4 1e−5 1e−6 1e−7 1e−8 1e−9 0 1 2 3 4 5 6 7 Regularization parameter Localization error in mm Minus 25% Minus 10% As Is Plus 10% Plus 25% (a) (b) 1e−2 1e−3 1e−4 1e−5 1e−6 1e−7 1e−8 1e−9 0 1 2 3 4 5 6 7 Regularization parameter Localization error in mm Minus 25% Minus 10% As Is Plus 10% Plus 25% 1e−2 1e−3 1e−4 1e−5 1e−6 1e−7 1e−8 1e−9 0 1 2 3 4 5 6 7 Regularization parameter Localization error in mm Minus 25% Minus 10% As Is Plus 10% Plus 25% (c) (d) Figure 6.6: Localization error as a function of wavelengthλ and of perturbation in opti- cal properties. (a) and (c) show the error due to perturbation of optical properties alone, where as (b) and (d) show the combined effect of perturbation in the optical properties and change in skull thickness. Smoothing was performed by applying a simple spa- tial Gaussian filter of width 0.4 mm to the solution, to eradicate effects of outliers and numerical errors. the skull are overestimated, the source tends to move deeper into the structure, whereas an underestimation of the optical properties leads to more superficial localizations. The choice of the regularization parameter has a significant impact on reconstruction results. While the increased regularization can compensate for the misestimation of the optical properties, it cannot reduce the error completely, if the skull thickness is also overesti- mated. Correct estimation of the optical properties and the geometry are crucial for accurate 3D localization of bioluminescent probes. However, if the deviations are not too large (i.e. in the range of -25% to 25%), the true model and the estimated model share a subspace of sufficiently high dimension to allow for a correct localization of the source at the cost of resolution. 72 6.4 A multi-modality PET-OBT study for comparing the localization of tumors in an inhomogeneous mouse phantom 600 620 640 660 680 700 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 Wavelength (nm) Reduced Scattering coefficient (mm −1 ) Skin Other organs Heart Liver Lungs Skeleton Kidneys Stomach 600 620 640 660 680 700 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 Wavelength (nm) Absorption coefficient (mm −1 ) Skin Other organs Heart Liver Lungs Skeleton Kidneys Stomach (a) (b) Figure 6.7: The optical properties for organs in the inhomogeneous mouse atlas; (a)Reduced scattering coefficientsμ ′ s and (b) absorption coefficientsμ a An anatomical mouse atlas based on x-ray CT and cryosection images [SCS + 02] was used to define the surface of a 3D mouse and of 17 internal organs. The vol- ume was tessellated using the constrained delaunay tetrahedralization method, such that organ boundaries were preserved. The resulting mesh had 109994 elements and 20989 nodes (7338 surface nodes and 13651 internal nodes). 10345 source locations were defined on an uniform grid inside the mouse volume with a 1.2 mm grid spacing. Opti- cal properties were assigned element-wise for all 17 organs , and for 6 wavelength bins uniformly spaced at every 20 nm, between 600 nm to 700 nm, as shown in Figures 6.7a and 6.7b [CPW90, ARC05]. While in practice data are acquired in the focal plane of the CCD camera, for convenience we assumed here that measurements are collected directly on the surface of the animal. The photon flux at 2448 detector locations (a subset 73 of the total surface nodes) was computed for each of the 6 wavelength bins using soft- ware developed in-house. Bioluminescence activity was simulated in 3 tumors located in the brain, spleen and lungs respectively. The tumors were assumed to express fire- fly luciferase (FFL). The multispectral forward model was assembled as is described in [CDB + 05]. Noise was added to the data for realistic simulations. Reconstructions were computed using the positively-constrained preconditioned Conjugate Gradient (PCG) method [CDB + 05]. For the same geometry, a PET simulation was carried out. A system model of the microPET Focus 220 scanner was generate using the method in [QLC + 98]. The sinogram was with max ring difference (MRD) 47 and span 3, the image size was 128× 128× 97, and the voxel size was 0.4 mm in transverse plane and 0.8 mm in axial direction. For the PET reconstruction, the atlas mouse was used as the background and activity in the three tumors was simulated. Poisson noise was added to the sino- gram so that the average count rate is 20 count/LOR. The reconstruction was done by 2 iterations of the OSEM3D method and 18 iterations of MAP with uniform resolu- tion [QL00, LAQL04]. Overlay of the reconstructed images from PET and MOBT are shown in Figure 6.8. The locations and distributions of the tumors were reconstructed satisfactorily by MOBT. 6.5 Simultaneous reconstruction of sources with differ- ent emission spectra and applications using quan- tum dots The theory corresponding to the reconstruction of multiple probes was developed in section 2.8. In this section, we present a study where we simulated three tumors at different locations inside the mouse volume. The tumors were assumed to be tagged 74 Figure 6.8: Overlay of the PET and optical reconstruction in a inhomogeneous mouse with tumors in the head, spleen and the lungs by molecules that corresponded to the emission spectra of Texas Red, Red beads and Cy5.5 respectively. Results are shown for the simultaneous reconstruction of the tumors from noisy hyperspectral bioluminescent data. Further, the same idea can be applied for HOFT by using quantum dots as tags. Quantum dots are advantageous for HOFT because the same illumination source can excite multiple quantum dots whose emission spectra may be different [MPB + 05]. The method described here can be used for the localization of all dots simultaneously. Simulation setup The mouse geometry and tumor locations used for simulation, and the probe spectra are shown in Figures 6.9a and 6.9b respectively. Data corresponding to 5272 sources placed on a 3D regular grid (grid spacing 1.2 cm) were simulated and were assumed to 75 Figure 6.9: The animal geometry and the emission spectra for the probes, (a) the three tumor locations in 3D and (b) the corresponding emission spectra be measured at 2878 detectors placed on the surface of the mouse in the central region. Gaussian noise was added to the data. The hyperspectral forward model was assem- bled using the extrapolated boundary method and images were reconstructed using the truncated pseudo-inverse of the forward model. Results and discussion Results of tumor reconstruction from the simulation study are shown in Figure 6.10. These reconstructions clearly indicate our ability to resolve multiple probes that could reflect different underlying biological processes or cells. Here we simulated an HOBT- like set up in which a single view is acquired corresponding to HOFT with a single uniform illumination within the animal. Quantum dots may be used for tagging genes or proteins in HOBT and the method described here can be used for estimating the position of quantum dots in 3D. 76 Figure 6.10: Horizontal sections of the mouse torso showing simultaneous reconstruc- tion of the Texas red, Red beads and Cy5.5 biological probes in-vivo. The first row represents the original probe location, while the second row shows the reconstruction result. 6.6 A feasibility study exploring the use of OFT in human breast cancer imaging HOFT has an additional advantage compared to HOBT that the illumination profile can be varied to elicit a different mapping from the sources to the detectors (i.e. a different HOFT forward model is computed for each illumination). Consequently, simultane- ous use of data from multiple illuminations increases the conditioning of the forward problem. This is especially important for the reconstruction of deep sources, e.g. in human breasts where tumors can be as deep as 6 cm. The fluorescent dye Indocyanine Green (ICG) has been approved by US Food and Drug Administration (FDA) for use in humans [NYSC00, GHA05]. The emission spectrum of ICG peaks at 805 nm, where the absorption coefficient of tissue is very small. Consequently, light may travel to the sur- face even from a depth of 6-8 cm. Metal-enhanced ICG has been developed to increase 77 the quantum yield by a factor of 25 and with additional resistance to photobleaching [MGGL03]. These advancements motivate this feasibility study of the application of HOFT to breast tumor imaging in humans. Figure 6.11: The tessellated human breast phantom The simulation setup Human breasts were assumed to be hemispherical with a diameter of 12 cm. The tes- sellated breast volume used for this simulation is shown in Figure 6.11. The tessellation had 501,324 tetrahedrons connected at 89,026 nodes. Hyperspectral data corresponding to 8003 sources placed on a regular grid inside the volume were simulated and assumed to be measured at 1413 detectors placed on the surface of the breast as shown in Figure 6.11b. Five point-like illumination profiles were used. The forward model was cal- culated by our fast FEM-based technique using published optical properties of human breasts. Realistic noise was added to the data and images were reconstructed using the truncated pseudo-inverse of the forward model. 78 Results and discussion Reconstruction results from simulation studies are shown in Figure 6.12. A point-like source at about 6 cm depth can accurately be localized by using HOFT. This feasibility study demonstrates that HOFT can potentially be a powerful modality for early breast cancer detection because of it’s high sensitivity, non-ionizing and non-radioactive prop- erties, specificity, and accuracy of 3D tumor localization. Figure 6.12: HOFT reconstruction of a point-like source at a depth of 6 cm in the human breast: (a) The original source location and (b) The reconstructed source distribution estimated using noisy data. 6.7 Phantom studies to analyze localization error and resolution in OBT Seven identical rectangular-block shaped phantoms (size 24 mm×24 mm×14 mm) were made, and a single phantom was placed on the mirror setup for each experiment. The optical fiber driven by the orange LED was inserted into the phantom at a particular location by means of a glass capillary tube marked to indicate depth of insertion. The experiment was repeated 7 times, with the optical fiber tip placed at a different location each time. Of the 7 locations, two were along the trans-axial X direction with the (x, 79 Figure 6.13: Diagram of the locations where the optical fiber tip was placed along two orthogonal planes for the localization and resolution studies. The points are ordered as 1: (12, 11, 6) mm, 2: (12, 18, 6) mm, 3: (12, 6, 6) mm, 4: (4, 11, 6) mm, 5: (20, 11, 6) mm, 6: (12, 11, 10) mm and 7: (12, 11, 4) mm. The two points in the trans-axial X directions are indicated by gray circles, those along the trans-axial Y direction are denoted by white circles, those along the axial direction are denoted by black circles and the location close to the center of the phantom is shown as a white square. y, z) coordinates (4, 11, 6) mm and (20, 11, 6) mm, two were along the trans-axial Y direction at (12, 6, 6) mm and (12, 18, 6) mm, one source was close to the center of the phantom at (12, 11, 6) mm and two sources were placed along the axial direction at locations (12, 11, 4) mm and (12, 11, 10) mm. Figure 6.13 shows these 7 locations along thez = 6 mm andy = 11 mm planes. The locations were verified after the experiments by cutting open the phantoms along the horizontal and vertical planes containing the fiber. Multispectral data were acquired in six spectral bins for all four orthogonal views with an acquisition time of 15 sec for each spectral bin. The achromatic data were collected without using any filter. The forward model was computed using the analytic model by defining a uniform mesh with 11745 nodes (sample spacing 0.85 mm) inside the rectangular block and 2668 surface nodes (sample spacing 0.85 mm). 80 Figure 6.14: Achromatic and multispectral reconstructions using the rectangular-block shaped phantoms for a single point source. The color scale used here goes from light (no intensity) to dark (maximum intensity); (a) displays the horizontal sections through the phantom showing the true location (12, 11, 6) mm of the tip of the orange LED-driven optical fiber; (b) horizontal sections showing reconstructions from achromatic data; (c) horizontal sections showing reconstruction from multispectral data. Reconstructions were computed using the PCG method. Since a positivity constraint can have the effect of artificially enhancing resolution for point sources in a zero back- ground, we dropped the non-negativity constraint from the cost function (4.5) for this study. Figure 6.14 shows a reconstruction of a point source located at (12, 11, 6) mm using multispectral data from 6 spectral bins versus reconstruction using achromatic data. The reconstruction using the latter produces superficial sources, whereas the for- mer reconstructs the source at the correct depth. Localization error for all seven sources was observed to be≤ 0.6 mm. Table 6.1 shows the measured FWHM resolution for the points. In theX andY directions, we observe a worst case resolution of 2.53 mm and 2.43 mm respectively at a depth of 6 mm near the center of the phantom. In theZ direc- tion we achieve a worst case resolution of 2.24 mm at a depth of 6 mm. For achromatic data, resolution cannot be defined, since the sources are not localized correctly. 81 True 3D location Resolution X (mm) Resolution Y (mm) Resolution Z (mm) (12, 11, 6) 2.53 2.43 2.24 (12, 18, 6) 2.04 1.91 2.14 (12, 6, 6) 2.10 2.21 2.11 (12, 11, 4) 1.46 1.47 1.58 (12, 11, 10) 1.58 1.61 1.57 (20, 11, 6) 2.01 2.38 2.08 (4, 11, 6) 1.96 2.22 2.13 Table 6.1: Table showing the FWHM at various locations for point-like sources placed in the rectangular phantom 6.8 An in-vivo study in a mouse The study was performed on a 26g athymic nude (nu/nu) mouse with a U87MG brain tumor (human glioblastoma) stereotactically implanted in the right cerebral hemisphere. Computed tomography was performed on the animal using a MicroCAT R scanner (Imtek, Inc.) with a 3.5 cm transaxial field of view and spatial resolution of 0.072 mm. MRI scans of the animal were collected using a 7.05T Bruker Biospin scan- ner(no. B-C-70/16US). The imaging protocol provided a transaxial resolution of 150 to 230 microns and an axial resolution of 670 microns. For the bioluminescence exper- iment, the mouse was anesthetised by an intraperitoneal (i.p.) injection of a mixture of ketamine and xylazine (75 and 15 mg/kg, respectively). Additional half doses were administered as needed to maintain anesthesia over the course of the imaging procedure. Luciferin (Xenogen Biosciences; 5μg/g) was injected i.p.. The mouse was placed prone on the mirror setup inside the imaging chamber of the IVIS 200 so as to enable imaging from the dorsal side. The subject platform height was adjusted to bring the dorsal surface of the mouse into focus. A series of unfiltered biolu- minescence images with integration time of 2 min each were collected from 2 min to 30 min and 37 min to 50 min after the luciferin injection in order to characterize the inten- sity variation of the bioluminescent signal over time. Beginning at 50 min, a sequence 82 Figure 6.15: Bioluminescence data obtained with the IVIS 200 imaging system and mir- ror setup overlaid on white-light photographs of a mouse with an implanted brain tumor. (1,8) are achromatic (i.e. unfiltered) acquisitions obtained for calibration purposes.(2)- (7) were acquired using the six filters (560, 580, 600, 620, 640 and 660 nm) mounted on the system in order of increasing wavelength. Since no signal was observed in the bottom view, only the top and side-views were acquired. of eight 2 min images (1 unfiltered, 6 using the filters on the filter wheel, 1 unfiltered) were acquired. The unfiltered images were used to correct for variation in biolumines- cence intensity over the course of the spectral image sequences by linear interpolation between the peak intensity values from the two images. The height of the platform was then adjusted to place the side views in focus, and the eight-image sequence was repeated. Because the tumor was located close to the dorsal surface of the brain, no signal was obtained in the bottom view. Each unfiltered bioluminescence acquisition and spectral imaging sequence was preceded by an white-light photograph to record the position of the mouse. The region of interest (the animal head) from the photograph of the top and side-views was aligned with the microCT images by a five parameter affine fit (3D translation plus rotations about two orthogonal axes) to minimize differences between the projected co-registered microCT and the photograph. For the microCT, the volume was integrated over the corresponding axis to form the silhouette image, which 83 was then matched to the respective photograph. The 3D CT data set was summed along the z and x axes, respectively to obtain the dorsal and lateral views. Only silhouettes were matched, i.e. the images used for the match were binary. The optical data were mapped onto the surface of the anatomic volume using the same affine transformation as for the photographs. The skull and soft tissue were segmented from the microCT using BrainSuite 2.0 [SL02], and custom software was used to generate a mesh as input for TOAST. Optical properties for skull and soft tissue at each node of the tessellated volume were assigned as per published data [CPW90](Skull: absorption coefficients in mm −1 of 0.0625, 0.0593, 0.0561, 0.0528 and reduced scattering coefficients in mm −1 of 2.1947, 2.1589, 2.1230, 2.0871 corresponding to wavelengths 580 nm, 600 nm, 620 nm and 640 nm, Soft tissue: absorption coefficients in mm −1 of 0.0283, 0.0260, 0.0237, 0.0222 and reduced scattering coefficients in mm −1 of 2.7085, 2.5000, 2.3128, 2.1597, again, corresponding to wavelengths 580 nm, 600 nm, 620 nm and 640 nm). The for- ward model was computed using TOAST for four wavelength bins (corresponding to data acquired by filters with peaks at 580 nm, 600 nm, 620 nm and 640 nm). Each for- ward model consisted of 10,814 internal and 58,000 surface nodes and was computed in approximately 30 hours on a AMD OPTERON 250 dual CPU computer. Reconstruc- tions were computed using the PCG method to minimize the cost function in (4.5), in this case including the non-negativity penalty. Figure 6.15 shows the sequence of bioluminescent data collected by the IVIS 200 system overlaid on the photograph of the mouse with the camera focused on the top- view. The results of the 3D reconstruction are shown in Figure 6.16, where we have overlaid the bioluminescence reconstruction on the horizontal MRI slices. The MRI images were aligned with the microCT images using RVIEW [SHH99]. The overlaid images demonstrate good anatomic localization of the bioluminescence image of the tumor, which is centered slightly dorsal to the MRI image of the tumor 84 Figure 6.16: Reconstruction of six horizontal planes of the bioluminescence data over- laid on co-registered MR slices. The red contour indicates the boundary of the tumor on the corresponding slices. The threshold was set to display bioluminescence data greater than or equal to10% of the maximum reconstructed voxel value. 6.9 Summary Results from simulation studies in the mouse atlas clearly indicate advantages of using hyperspectral/multispectral data over achromatic data. This is further confirmed by phantom studies. A resolution of about 2.4 mm is achieved in practical settings, though simulation studies indicate that better resolution will be achieved if hyperspectral data is used. The in-vivo experiment clearly indicates tremendous potential for 3D HOBT and HOFT for molecular imaging in small animals. 85 Chapter 7 Summary and future work 7.1 Conclusions We have demonstrated a method for 3D multispectral bioluminescence optical tomog- raphy (MOBT) and hyperspectral bioluminescence optical tomography (HOBT) and presented the merits of using spectral information for optical reconstruction. While the FEM-based method allowed for incorporation of tissue optical properties and yielded accurate solutions to the diffusion equation, it was very time consuming. The develop- ment of fast forward model solvers presented in this thesis was essential for the practi- cality of MOBT/HOBT for in-vivo studies. Work is currently in progress for the evalu- ation of fast inversion schemes for optical tomography. A combination of fast forward and inverse solvers for MOBT/HOBT potentially leads to more accurate results with reasonable computation time. The spectrum observed on the animal’s skin due to an underlying source distribution encodes depth information. Thus, hyperspectral and multispectral models have more information content, especially about deep sources, as compared with achromatic data. This is very clear from the 3D distributions of the significant right singular vectors of the achromatic and hyperspectral models as shown in Figures 6.1b and 6.1c respectively. While the singular-value spectrum is an effective measure of information content and ill- conditioning of the problem, it is the spatial distribution of the principal right singular vectors that governs the ability to localize deep or superficial sources. Each spectral measurement may be viewed as adding some amount of independent information to the 86 problem, thus making the inversion less ill-posed. As a direct consequence, hyperspec- tral and multispectral inversions are more robust to measurement noise. This was shown in the phantom studies, where the multispectral data produced superior results as com- pared to the achromatic measurements. The achromatic model is a projection of the monochromatic models onto a single depth-independent spectral basis function (a rect- angular function) and thus, has less information content compared to the hyperspectral or multispectral models. From the reconstruction results in Figures 6.2b and 6.2c, it is clear that the achromatic model is already superior to the monochromatic models com- puted using average values for optical properties. However, it should also be noted that the current published results for bioluminescent tomography [GZLJ04] use the latter. Hyperspectral data acquisition requires more sophisticated hardware [ZVM + 06] than multispectral data collection (which can be accomplished with bandpass filters). However, the hyperspectral forward model is expected to be better conditioned than the multispectral forward model because it incorporates more information about wavelength and should perform better for lower signal-to-noise ratios. Simulation and phantom studies conducted and presented in this thesis are useful for the description of resolution limits and of localization accuracy for setups and optical instruments currently used in research laboratories. While our resolution studies have been limited so far to phantoms and simulated mice, the reconstruction of the brain tumor in a live mouse clearly indicates great promise for applications of our method in functional imaging of small animals. Results from the in-vivo mouse study show the reconstructed bioluminescent image of the tumor to be slightly more dorsal than the anatomical tumor. This small mis-localization of the tumor might have been caused by mis-estimation of the skull optical properties. Studies using quantum dots for multiple 87 probe detection and those exploring the feasibility of OFT in human breast tumor imag- ing show potential of the approaches presented in this thesis to complex pre-clinical and clinical studies. 7.2 Future work 7.2.1 Innovative methods for spatial encoding of probes in OFT Bioluminescence imaging is limited in the sense that light is produced only from cellular-level interactions unlike fluorescence imaging, where external illumination can be used to produce different mappings of photon fluence on the animal surface. The propagation of the photons in tissue is a function of the wavelength. Thus, if the illumi- nation source emitted at wavelengths below 550 nm, only superficial sources will receive excitation radiation and will emit. On the other hand, external illumination at wave- lengths greater than 750 nm will correspond to larger penetration depth, and effectively, deeper sources will receive radiation. Thus, control over the excitation wavelength may be useful in order to spatially encode the location of the fluorescent probes. We are in the process of doing some phantom and simulation studies investigating this. Of great interest is further investigation into optimal illumination schemes for OFT. Methods of using lifetime information about the fluorescent probe have been proposed by using a modulated external source [RSM01, CBDT05]. This class of methods will especially be useful for distinguishing multiple probes and is expected to reduce the illposedness of the inverse problem. These methods in conjunction with spectral information could be further investigated. 88 7.2.2 Fast forward and inverse solvers Ordered subset methods e.g. the globally convergent relaxed OSEM (Ordered subset expectation maximization) methods have been proposed in PET [AF03] that use a subset of the measurement data for each update instead of the total data. The convergence rate of these methods may further be largely affected by the choice of the preconditioner. Coordinate descent methods have been proposed for use in PET [MOW + 03]. These and several other emission computed tomography inversion methods [QL06] could be investigated for OBT and OFT in order to achieve fast convergence to the solution. 89 Bibliography [ADD04] P.R. Amestoy, T. A. Davis, and I. S. Duff. Algorithm 837: AMD, an approximate minimum degree ordering algorithm. ACM Transactions on Mathematical Software, 30(3):381–388, Sep 2004. [ADSO00] S. R. Arridge, H. Dehghani, M. Schweiger, and E. Okada. The finite ele- ment model for propagation of light in scattering media: A direct method for domains with nonscattering regions. Medical Physics, 27(1):252–264, 2000. [AF03] S. Ahn and J. A. Fessler. Globally convergent image reconstruction for emission tomography using relaxed ordered subset methods. IEEE Trans. Medical Imaging, 22:613–626, 2003. [AH03] G. S. Abdoulaev and A. H. Hielscher. Three-dimensional optical tomogra- phy with the equation of radiative transfer. Journal of Electronic Imaging, 12(4):594–601, October 2003. [ARC05] G. Alexandrakis, F. R. Rannou, and A. F. Chatziioannou. Tomographic bioluminescence imaging by use of a combined optical-pet (opet) system: a computer simulation feasibility study. Physics in Medicine and Biology, 50(17):4225–4241, 2005. [Arr95] S. R. Arridge. Photon-measurement density functions. part i: Analytical forms. Applied Optics, 34(31):7395–7409, 1995. [Arr99] S.R. Arridge. Optical tomography in medical imaging. Inverse Problems, 15(2):41–93, 1999. [ASHD93] S. R. Arridge, M. Schweiger, M. Hiroaka, and D. T. Delpy. A finite ele- ment approach for modeling photon transport in tissue. Phys. Med. Biol., 20(2):299–309, 1993. [BBM + 01] D.A. Boas, D. H. Brooks, E. L. Miller, C. A. DiMarzio, M. Kilmer, R. J. Gaudette, and Quan Zhang. Imaging the body with diffuse optical tomog- raphy. IEEE Signal Processing Magazine, 18(6):57–75, November 2001. 90 [CBDT05] D. J. Cuccia, F. Bevilacqua, A. J. Durkin, and B. J. Tromberg. Modu- lated imaging: quantitative analysis and tomography of turbid media in the spatial-frequency domain. Optics Letters, 30(11):1354–1356, June 2005. [CDB + 05] A. J. Chaudhari, F. Darvas, J. R. Bading, R. A. Moats, P. S. Conti, D. J. Smith, S. R. Cherry, and R. M. Leahy. Hyperspectral and multispectral bioluminescence optical tomography for small animal imaging. Physics in Medicine and Biology, 50(23):5421–5441, Nov 2005. [CDCL05] A. J. Chaudhari, F. Darvas, S. R. Cherry, and R. M. Leahy. Resolution stud- ies in bioluminescence and fluorescence optical tomography using hyper- spectral data acquisition. Molecular Imaging and Biology, 7(6):125–126, 2005. [CDS + 03] A. J. Chaudhari, F. Darvas, D. J. Smith, S. R. Cherry, and R. M. Leahy. Hyperspectral fluorescence optical tomography for small animal imaging. Molecular Imaging, 2:225–226, August 2003. [CJDL07] A. J. Chaudhari, A. A. Joshi, F. Darvas, and R. M. Leahy. A method for atlas-based volumetric registration with surface constraints for optical bioluminescence tomography in small animal imaging. In SPIE Medical Imaging. SPIE - International Society for Optical Engineering, 2007. [CPT + 97] R. Cubeddu, A. Pifferi, P. Taroni, A. Torricelli, and G Valentini. A solid tissue phantom for photon migration studies. Phys. Med. Biol., 42:1971– 1979, 1997. [CPW90] W. F. Cheong, S. A. Prahl, and A. J. Welch. A review of the optical proper- ties of biological tissues. IEEE Journal of Quantum Electronics, 26:2166– 2185, 1990. [CR02] C. H. Contag and B. D. Ross. It’s not just about anatomy: In vivo bio- luminescence imaging as an eyepiece into biology. Journal of Magnetic Resonance Imaging, 16(4):378–387, September 2002. [CZ67] M. C. Case and P. F. Zweifel. Linear Transport Theory. Addison-Wesley Inc., 1967. [Dau92] I. Daubechies. Ten Lectures on Wavelets. Society of Industrial and Applied Mathematics, Philadelphia, PA, 1992. [DCL05] F. Darvas, A. J. Chaudhari, and R. M. Leahy. Effect of skull optical prop- erties on hyperspectral optical tomography in small animals. Molecular Imaging, 5, 2005. 91 [DLA + 99] R. M. Doornbos, R. Lang, M. C. Aalders, F. W. Cross, and H. J. Sternborg. The determination of in vivo human tissue optical properties and absolute chromophore concentrations using spatially resolved steady-state diffuse reflectance spectroscopy. Physics in Medicine and Biology, 44:967–981, 1999. [FHED93] M Firbank, M Hiraoka, M Essenpreis, and D T Delpy. Measurement of the optical properties of the skull in the wavelength range 650-950 nm. Physics in Medicine and Biology, 38(4):503–510, 1993. [FVM + 01] B. K. Ford, C. E. V olin, S. M. Murphy, R. M. Lynch, and M. R. Descour. Computer tomography-based spectral imaging for fluorescence microscopy. Biophysical Journal, 80(2):986–993, 2001. [GFB83] R. A. J. Groenhius, H. A. Ferwerda, and J. J. Ten Bosch. Scattering and absorption of turbid materials determined from reflection measurements 1: Theory. Applied Optics, 22(16):2456–2462, 1983. [GHA05] A P Gibson, J C Hebden, and S R Arridge. Recent advances in diffuse optical imaging. Physics in Medicine and Biology, 50(4):R1–R43, 2005. [GL81] A. George and J Liu. Computer Solutions of Large Sparse Positive Definite Systems. Prentice-Hall, 1981. [GRWN03] E. E. Graves, J. Ripoll, R. Weissleder, and V . Ntziachristos. A submil- limeter resolution fluorescence molecular imaging system for small animal imaging. Medical Physics, 30(5):901–911, 2003. [GV96] G. H. Golub and C. F. Van Loan. Matrix Computations. The John Hopkins University Press, MD, 3 edition, 1996. [GZLJ04] X. Gu, Q. Zhang, L. Larcom, and H. Jiang. Three-dimensional biolumi- nescence tomography with model based reconstruction. Optics Express, 12(17), August 2004. [Han92] P. C. Hansen. Analysis of discrete ill-posed problems by means of the l-curve. SIAM Review, 32:561–580, 1992. [Hec01] E. Hecht. Optics. Pearson Addison Wesley, 4th edition edition, August 2001. [HST + 94] R. C. Haskell, L. O. Svaasand, T. T. Tsay, T. C. Feng, M. S. McAdams, and B. J. Tromberg. Boundary conditions for the diffusion equation in radiative transfer. J. Opt. Soc. Am. A, 11(10):2727–2740, 1994. 92 [Jai88] A.K. Jain. Fundamentals of Digital Image Processing. Prentice Hall, 1 edition, September 1988. [JSTL06] A. A. Joshi, D. Shattuck, P. M. Thompson, and R. M. Leahy. Simultaneous surface and volume registration using 2d and 3d mappings. In SIAM Imag- ing Sciences Abstracts, page 55. SIAM Conference on Imaging Sciences, Society of Industrial and Applied Mathematics, 2006. [Kau87] L. Kaufman. Implementing and accelerating the em algorithm for positron emission tomography. IEEE Transactions on Medical Imaging, 6(1):37– 51, March 1987. [KH03] A. D. Klose and A. H. Hielscher. Fluorescence tomography with sim- ulated data based on the equation of radiative transfer. Optics Letters, 28(12):1019–1021, 2003. [KSS88] M. Keijzer, W. M. Star, and P. R. M. Storchi. Optical diffusion in layered media. Applied Optics, 27(9):1820–1824, 1988. [Lan02] D. Landgrebe. Hyperspectral image data analysis. IEEE Signal Processing Magazine, 19(1):17–28, 2002. [LAQL04] Q. Li, E. Asma, J. Qi, and R. M. Leahy. Accurate estimation of the fisher information matrix for the pet image reconstruction problem. IEEE Trans. Medical Imaging, 23:1057–1064, 2004. [LH00] R. M. Levenson and C. C. Hoyt. Spectral imaging and microscopy. Amer- ican Laboratory, 32(22):26+, 2000. [Lue84] D. Luenberger. Linear and Nonlinear Programming. Addison-Wesley Inc., Massachusetts, 2 edition, 1984. [MGGL03] J. Malicka, I. Gryczynski, C. D. Geddes, and J. R. Lakowicz. Metal- enhanced emission from indocyanine green: a new approach to in vivo imaging. Journal of Biomedical Optics, 8(3):472–478, 2003. [MLCZ94] E. Mumcuoglu, R. M. Leahy, S. R. Cherry, and Z. Zhou. Fast gradient based methods for bayesian reconstruction of transmission and emission pet images. IEEE Transactions on Medical Imaging, 13:687–701, 1994. [MOW + 03] A. B. Milstein, S. Oh, K. J. Webb, C. A. Bouman, Q. Zhang, D. A. Boas, and R. P. Millane. Fluorescence optical diffusion tomography. Applied Optics, 42(16):3081–3094, 2003. 93 [MPB + 05] X. Michalet, F. F. Pinaud, L. A. Bentolila, J. M. Tsay, S. Doose, J. J. Li, G. Sundaresan, A. M. Wu, S. S. Gambhir, , and S. Weiss. Quantum dots for live cells, in vivo imaging, and diagnostics. Science, 307(5709):538–544, 2005. [MTBW99] U. Mahmood, C. H. Tung, A. Bogdanov, and R. Weissleder. Near- infrared optical imaging of protease activity for tumor detection. Radi- ology, 213(3):866–870, 1999. [MZ93] S. Mallat and Z. Zhang. Matching pursuit with time-frequency dictionar- ies. IEEE Transactions on Signal Processing, 41:3397–3415, 1993. [NYSC00] V . Ntziachristos, A. G. Yodh, M. Schnall, and B. Chance. Concurrent MRI and diffuse optical tomography of breast after indocyanine green enhance- ment. Proc Natl Acad Sci USA, 97(6):2767–2772, March 2000. [QL00] J. Qi and R. M. Leahy. Resolution and noise properties of map reconstruc- tions in fully 3d pet. IEEE Trans. Medical Imaging, 19(5):493–506, May 2000. [QL06] J. Qi and R. M. Leahy. Iterative reconstruction techniques in emission computed tomography. Physics in Medicine and Biology, 51(15):R541– R578, 2006. [QLC + 98] J. Qi, R. M. Leahy, S. R. Cherry, A. Chatziioannou, and T. H. Farquhar. High resolution 3d bayesian image reconstruction using the micropet small-animal scanner. Physics in Medicine and Biology, 43(4):1001–1013, 1998. [RCN01] B. W. Rice, M. D. Cable, and M. B. Nelson. In vivo imaging of light- emitting probes. Journal of Biomedical Optics, 6(4):432–440, 2001. [RSM01] R. Roy and E. M. Sevick-Muraca. Three-dimensional unconstrained and constrained image-reconstruction techniques applied to fluorescence, frequency-domain photon migration. Applied Optics, 40(13):2206–2215, 2001. [Saa96] Y . Saad. Iterative Methods for Sparse Linear Systems. PWS Publishing Company, 1996. [SAHD95] M. Schweiger, S. R. Arridge, M Hiroaka, and DT Delpy. The finite element method for the propagation of light in scattering media: Boundary and source conditions. Medical Physics, 22(11):1779–1792, 1995. 94 [SCS + 02] D. Stout, P. Chow, R. Silverman, R. Leahy, X. Lewis, S. Gambhir, and A. Chatziioannou. Creating a whole body digital mouse atlas with pet, ct and cryosection images. Molecular Imaging and Biology, 4(4), 2002. [She94] J. R. Shewchuk. An introduction to the conjugate gradient method without the agonizing pain. http://www-2.cs.cmu.edu/ jrs/jrspapers.html, August 1994. [SHH99] C. Studholme, D. L. G. Hill, and D. J. Hawkes. An overlap invariant entropy measure of 3d medical image alignment. Pattern Recognition, 32(1):71–86, 1999. [SL02] D. W. Shattuck and R. M. Leahy. Brainsuite: An automated cortical surface identification tool. Medical Image Analysis, 6(2):129–142, 2002. [SNZ + 01] R. A. Schultz, T. Nielsen, J. R. Zavaleta, R. Ruch, R. Wyatt, and H. R. Garber. Hyperspectral imaging: a novel approach for microscopic analysis. Cytometry, 43(4):239–247, 2001. [UOO + 01] Y . Ueda, K. Ohta, M. Oda, M. Miwa, Y . Tsuchiya, and Y . Yamashita. Three-dimensional imaging of a tissuelike phantom by diffusion optical tomography. Applied Optics, 40(34):6349–6355, 2001. [VK95] M. Vetterli and J. Kovecevic. Wavelets and subband coding. Prentice Hall Signal Processing Series. Prentice-Hall PTR., Englewood Cliff, NJ, 1st edition, 1995. [WLJ04] G. Wang, Y . Li, and M. Jiang. Uniqueness theorems in bioluminescence tomography. Medical Physics, 31(8):2289–2299, August 2004. [WN03] R. Weissleder and V . Ntziachristos. Shedding light onto live molecular targets. Nature Medicine, 9(1):123–128, 2003. [WSD + 06] G. Wang, H. Shen, K. Durairaj, X. Qian, and W. Cong. The first bio- luminescence tomography system for simulteneous acquisition of multi- view and multi-spectral data. International Journal of Biomedical Imag- ing, 2006. accepted. 95 [YWP + 97] Y . Yao, Y . Wang, W. Pei, W. Zhu, and R. L. Barbour. Frequency-domain optical imaging of absorption and scattering distributions by a born iter- ative method. J. Opt. Soc. Am. A-Optics Image Science and Vision, 14(1):352–342, 1997. [ZVM + 06] G. Zavattini, S. Vecchi, G. Mitchell, U.Weisser, R. M. Leahy, B. J. Pichler, D. J. Smith, and S. R. Cherry. A hyperspectral fluorescence system for 3d in vivo optical imaging. Physics in Medicine and Biology, 51(8):2029– 2043, 2006. 96

Abstract (if available)

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Fast solvers in hyperspectral optical bioluminescence tomography for small animal imaging

PDF

Image registration with applications to multimodal small animal imaging

PDF

Fast iterative image reconstruction for 3D PET and its extension to time-of-flight PET

PDF

Contributions to structural and functional retinal imaging via Fourier domain optical coherence tomography

PDF

Development of a multi-mode optical imaging system for preclinical applications in vivo

PDF

3D vessel mapping techniques for retina and brain as an early imaging biomarker for small vessel diseases

Asset Metadata

Creator Chaudhari, Abhijit J. (author)

Core Title Hyperspectral and multispectral optical bioluminescence and fluorescence tomography in small animal imaging

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Electrical Engineering

Publication Date 01/26/2007

Defense Date 12/13/2007

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag bioluminescence,fluorescence,hyperspectral,mutispectral,OAI-PMH Harvest,small animal imaging,tomography

Language English

Advisor Darvas, Felix (committee member), Leahy, Richard M. (committee member), Proskurowski, Wlodek (committee member)

Creator Email achaudhari@gmail.com

Permanent Link (DOI) https://doi.org/10.25549/usctheses-m232

Unique identifier UC186720

Identifier etd-Chaudhari-20070125a (filename),usctheses-m40 (legacy collection record id),usctheses-c127-158799 (legacy record id),usctheses-m232 (legacy record id)

Legacy Identifier etd-Chaudhari-20070125a.pdf

Dmrecord 158799

Document Type Dissertation

Rights Chaudhari, Abhijit J.

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Repository Name Libraries, University of Southern California

Repository Location Los Angeles, California

Repository Email cisadmin@lib.usc.edu