<%BANNER%>

Three-Dimensional Photoacoustic Tomography and Its Application to Detection of Joint Diseases in the Hand

Permanent Link: http://ufdc.ufl.edu/UFE0041907/00001

Material Information

Title: Three-Dimensional Photoacoustic Tomography and Its Application to Detection of Joint Diseases in the Hand
Physical Description: 1 online resource (128 p.)
Language: english
Creator: Sun, Yao
Publisher: University of Florida
Place of Publication: Gainesville, Fla.
Publication Date: 2010

Subjects

Subjects / Keywords: finger, multispectral, osteoarthritis, photoacoustic, three
Biomedical Engineering -- Dissertations, Academic -- UF
Genre: Biomedical Engineering thesis, Ph.D.
bibliography   ( marcgt )
theses   ( marcgt )
government publication (state, provincial, terriorial, dependent)   ( marcgt )
born-digital   ( sobekcm )
Electronic Thesis or Dissertation

Notes

Abstract: This thesis research presents the study of three-dimensional (3-D) photoacoustic tomography (PAT) and its application for the first time to detection of osteoarthritis (OA) in the hand. PAT is an emerging non-ionizing, non-invasive imaging modality that can visualize high optical contrast in biological tissues with high ultrasound resolution. Compared with other imaging modalities that have been conventionally used or recently investigated to visualize the structural abnormalities in the finger joints with OA, the 3-D PAT approach studied in this dissertation not only provides high-resolution anatomical structures, but also offers quantitative tissue optical property as well as physiological / functional information including concentrations of oxy-hemoglobin (HbO2), deoxy-hemoglobin (Hb) and water (H2O) that could be used to detect OA in an early stage. In this study, a 3-D high performance PAT reconstruction algorithm is developed based on parallel computing technique and finite element method. The optimal detector performance and scanning geometry for 3-D photoacoustic imaging of the finger joints are investigated with considerable phantom experiments. The results of phantom experiments show that a spherical scanning geometry appears to provide improved spatial solution over a cylindrical scanning geometry, and that 1mm thick ?cartilage? can be accurately differentiated from the ?bones? with a 1 MHz transducer in a spherical scanning geometry. In addition, the absorption coefficient of the ?cartilage? can be effectively recovered when this optical property varied from 0.015mm-1 to 0.04mm-1. A 3-D PAT system in a spherical scanning geometry has been constructed and optimized for in-vivo examination of human finger joints. A distal interphalangeal (DIP) joint from a female healthy subject was photoacoustically examined by our 3-D PAT system, and major anatomical structures of the examined DIP joint along with the side arteries were clearly reconstructed in high quality, where joint space was well differentiated from surrounding finger phalanges. The performance of the 3-D PAT system was further improved with an ultrasound detection array composed of eight 1 MHz transducers and a 16-channel pulse/receive board. The 8-channel 3-D PAT system was carefully calibrated with controlled tissue phantom experiments, and is capable of completing a finger joint examination within 5 minutes (at single optical wavelength). The 3-D PAT reconstruction algorithm and scanning system developed has also been applied to a pilot clinical study aiming to test the possibility of detecting OA in the hand joints using 3-D PAT. In this pilot clinical study, seven subjects (two OA patients and five healthy controls) were enrolled and photoacoustically examined. The image quality of the reconstructed finger joints was greatly improved with the 8-channel 3-D PAT system, and apparent differences, in both the reconstructed size of the joint space and the absorption coefficient of the joint cavity, has been observed between the OA and normal joints. The successful results obtained suggest the possibility of 3-D PAT as a potential clinical tool for early detection of OA in the finger joints. Major chromophore concentrations (HbO2 and Hb) of in-vivo finger joints have also been quantitatively imaged using multispectral 3-D PAT approach with six optical wavelengths from 730nm to 880nm. The multispectral results obtained further confirmed that the 3D PAT approach implemented in this thesis research is able to differentiate OA from normal joints. While we target the detection of OA as a testing base for validating the single- and multi-spectral 3-D PAT approaches developed in this thesis research, many aspects of our work are fundamental to imaging in general. For example, the 3-D PAT approaches implemented are applicable to other biomedical problems such as rheumatoid arthritis (RA) detection and functional brain imaging.
General Note: In the series University of Florida Digital Collections.
General Note: Includes vita.
Bibliography: Includes bibliographical references.
Source of Description: Description based on online resource; title from PDF title page.
Source of Description: This bibliographic record is available under the Creative Commons CC0 public domain dedication. The University of Florida Libraries, as creator of this bibliographic record, has waived all rights to it worldwide under copyright law, including all related and neighboring rights, to the extent allowed by law.
Statement of Responsibility: by Yao Sun.
Thesis: Thesis (Ph.D.)--University of Florida, 2010.
Local: Adviser: Jiang, Huabei.

Record Information

Source Institution: UFRGP
Rights Management: Applicable rights reserved.
Classification: lcc - LD1780 2010
System ID: UFE0041907:00001

Permanent Link: http://ufdc.ufl.edu/UFE0041907/00001

Material Information

Title: Three-Dimensional Photoacoustic Tomography and Its Application to Detection of Joint Diseases in the Hand
Physical Description: 1 online resource (128 p.)
Language: english
Creator: Sun, Yao
Publisher: University of Florida
Place of Publication: Gainesville, Fla.
Publication Date: 2010

Subjects

Subjects / Keywords: finger, multispectral, osteoarthritis, photoacoustic, three
Biomedical Engineering -- Dissertations, Academic -- UF
Genre: Biomedical Engineering thesis, Ph.D.
bibliography   ( marcgt )
theses   ( marcgt )
government publication (state, provincial, terriorial, dependent)   ( marcgt )
born-digital   ( sobekcm )
Electronic Thesis or Dissertation

Notes

Abstract: This thesis research presents the study of three-dimensional (3-D) photoacoustic tomography (PAT) and its application for the first time to detection of osteoarthritis (OA) in the hand. PAT is an emerging non-ionizing, non-invasive imaging modality that can visualize high optical contrast in biological tissues with high ultrasound resolution. Compared with other imaging modalities that have been conventionally used or recently investigated to visualize the structural abnormalities in the finger joints with OA, the 3-D PAT approach studied in this dissertation not only provides high-resolution anatomical structures, but also offers quantitative tissue optical property as well as physiological / functional information including concentrations of oxy-hemoglobin (HbO2), deoxy-hemoglobin (Hb) and water (H2O) that could be used to detect OA in an early stage. In this study, a 3-D high performance PAT reconstruction algorithm is developed based on parallel computing technique and finite element method. The optimal detector performance and scanning geometry for 3-D photoacoustic imaging of the finger joints are investigated with considerable phantom experiments. The results of phantom experiments show that a spherical scanning geometry appears to provide improved spatial solution over a cylindrical scanning geometry, and that 1mm thick ?cartilage? can be accurately differentiated from the ?bones? with a 1 MHz transducer in a spherical scanning geometry. In addition, the absorption coefficient of the ?cartilage? can be effectively recovered when this optical property varied from 0.015mm-1 to 0.04mm-1. A 3-D PAT system in a spherical scanning geometry has been constructed and optimized for in-vivo examination of human finger joints. A distal interphalangeal (DIP) joint from a female healthy subject was photoacoustically examined by our 3-D PAT system, and major anatomical structures of the examined DIP joint along with the side arteries were clearly reconstructed in high quality, where joint space was well differentiated from surrounding finger phalanges. The performance of the 3-D PAT system was further improved with an ultrasound detection array composed of eight 1 MHz transducers and a 16-channel pulse/receive board. The 8-channel 3-D PAT system was carefully calibrated with controlled tissue phantom experiments, and is capable of completing a finger joint examination within 5 minutes (at single optical wavelength). The 3-D PAT reconstruction algorithm and scanning system developed has also been applied to a pilot clinical study aiming to test the possibility of detecting OA in the hand joints using 3-D PAT. In this pilot clinical study, seven subjects (two OA patients and five healthy controls) were enrolled and photoacoustically examined. The image quality of the reconstructed finger joints was greatly improved with the 8-channel 3-D PAT system, and apparent differences, in both the reconstructed size of the joint space and the absorption coefficient of the joint cavity, has been observed between the OA and normal joints. The successful results obtained suggest the possibility of 3-D PAT as a potential clinical tool for early detection of OA in the finger joints. Major chromophore concentrations (HbO2 and Hb) of in-vivo finger joints have also been quantitatively imaged using multispectral 3-D PAT approach with six optical wavelengths from 730nm to 880nm. The multispectral results obtained further confirmed that the 3D PAT approach implemented in this thesis research is able to differentiate OA from normal joints. While we target the detection of OA as a testing base for validating the single- and multi-spectral 3-D PAT approaches developed in this thesis research, many aspects of our work are fundamental to imaging in general. For example, the 3-D PAT approaches implemented are applicable to other biomedical problems such as rheumatoid arthritis (RA) detection and functional brain imaging.
General Note: In the series University of Florida Digital Collections.
General Note: Includes vita.
Bibliography: Includes bibliographical references.
Source of Description: Description based on online resource; title from PDF title page.
Source of Description: This bibliographic record is available under the Creative Commons CC0 public domain dedication. The University of Florida Libraries, as creator of this bibliographic record, has waived all rights to it worldwide under copyright law, including all related and neighboring rights, to the extent allowed by law.
Statement of Responsibility: by Yao Sun.
Thesis: Thesis (Ph.D.)--University of Florida, 2010.
Local: Adviser: Jiang, Huabei.

Record Information

Source Institution: UFRGP
Rights Management: Applicable rights reserved.
Classification: lcc - LD1780 2010
System ID: UFE0041907:00001


This item has the following downloads:


Full Text





THREE-DIMENSIONAL PHOTOACOUSTIC TOMOGRAPHY AND ITS APPLICATION
TO DETECTION OF JOINT DISEASES IN THE HAND




















By

YAO SUN


A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL
OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT
OF THE REQUIREMENTS FOR THE DEGREE OF
DOCTOR OF PHILOSOPHY

UNIVERSITY OF FLORIDA

2010

































2010 Yao Sun





























And all who


To my wife, Xiaorong Li;
My mom and my dad;
have been supportive to me in my life









ACKNOWLEDGMENTS

First, I would like to thank Professor Huabei Jiang, who provided this great

education opportunity to me. I appreciate the full financial support and precious advice I

have received from Professor Jiang during my study for the PhD degree. His broad

views and deep understandings on the research subjects guided me all through this

research. His awareness to timetable and duty priorities, and his persistence to

research objectives deeply impressed me. What I have learned and experienced in

Professor Jiang's lab will definitely benefit my future career and life.

Secondly, my thanks go to my academic committee members: Dr. David Gilland,

Dr. Rosalind Sadleir, and Dr. Sihong Song. Their comments and suggestions were very

helpful to my study and research. Their valuable time spent on reading my dissertation

is highly appreciated.

Thirdly, I would like to thank Dr. Eric S. Sobel, an associate professor in the

College of Medicine at University of Florida. As our clinical research partner, Dr. Sobel

helps us to recruit sufficient patients with osteoarthritis. And I would like to thank the

research assistant professors in our lab, Dr. Qizhi Zhang and Dr. Zhen Yuan, for

valuable discussions on experiments and computations, respectively.

Last but not least, I would like to take this opportunity to express my thanks to Lei

Xi, Qiang Wang, Lei Yao, Yiyong Tan, Ruixin Jiang, Alexandria Grubbs, Xiaoping Liang,

Lu Yin, Lijun Ji, Bo Wang, and all members in our biomedical optics lab for their

assistance and the great time we spent together.









TABLE OF CONTENTS

page

ACKNOW LEDGMENTS .......... ....... .. .......... ... ......... .................... 4

LIST O F TA BLES .......... ..... ..... ............................................................. ........ 7

LIS T O F F IG U R E S .................................................................. 8

LIST OF ABBREVIATIONS ......................... ..................... 12

A BST RA C T ............... ... ..... ......................................................... ...... 13

CHAPTER

1 BACKGROUND AND SIGNIFICANCE .......... .... ...... .. ............... .............. 16

1.1 Introduction to Photoacoustic Tomography............... ...................... 16
1.2 Overview of Osteoarthritis...................... ....... ............................. 18

2 RECONSTRCUTION ALGORITHM IN PAT .......... ............................................ 24

2.1 Finite Element Based Reconstruction Algorithm in PAT.............................. 24
2.1.1 Forward Modeling in Acoustically Homogenous Media......................... 26
2.1.2 Inverse Modeling and Nonlinear Optimization in PAT Reconstruction .. 27
2.1.3 Regularization in PAT Reconstruction ............ .............. 30
2.1.4 Simulation in PAT Reconstruction ............ .................................. 33
2.2 Computing Strategies in PAT Reconstruction................................ ............. 33
2.2.1 A djoin Sensitivity M ethod .................................. .......... ...... ............ 34
2.2.2 D ual M eshing Schem e ............... ................... ................................... 34
2.2.3 Partial Reconstruction in Region of Interest ................ .................. 37
2.2.4 Parallel Com putting Technique .......................................... ............... 39

3 PRE-INVESTIGATION OF 3-D PAT SYSTEM WITH TISSUE PHANTOM STUDY47

3.1 Tissue M im picking Phantom ......... ............................................. ............... 48
3.2 System Description and Experiments .................................. ....................... 49
3 .3 M etho ds ........................ ..................... 50
3.4 Results and Discussion............................ ......... 52

4 IMAGING FINGER JOINTS IN-VIVO WITH 3-D PAT IN A SPHERICAL SCAN..... 63

4.1 System Development and Description ............... ........ .. ...... ............ 63
4.1 .1 Fiber Guided Near Infrared Lighting ............ .................................. 64
4.1.2 Arrangement for Hand and Finger Placement.................. ........... 65
4.1.3 Signal Detecting in a Spherical Scanning Geometry........................... 65
4.1.4 Units Control and Data Acquisition System............. ................... 67









4 .2 In-v ivo E x pe rim e nts .................................................................. 6 7
4 .3 M e th o d s .................................................. 6 8
4.4 Results and Discussion............................ ............... 69

5 CLINICAL EXAMINATION OF OA PATIENTS WITH 3-D PAT................................ 79

5.1 Enhancement in 3-D PAT Approach ............................ ...... ................ 80
5.1.1 U ltrasound Detection A rray .......... ...................................... ............ 81
5.1.2 M multiple Channel Data Acquisition................... ............. ............... 83
5.1.3 Quantitative 3-D PAT with 3-D Photon Diffusion Model ..................... 83
5.1.4 Tissue-like Phantom Test......... ........... .. ...... .................... 84
5.2 Patients and C clinical Data .......... ........................ ................ .............. 85
5.3 Results ............... ........ ................................................................... 86
5.4 Discussion ....................................................................................... ................... 89

6 IN-VIVO IMAGING of CHROMOPHORE CONCENTRATIONS in FINGER
JOINTS W ITH MULTISPECTRAL 3-D PAT...................................................... 100

6.1 M methods ............... ..... ...................................................... ........ 102
6.2 In-vivo Experiments .......................... ........................ 104
6.3 Results and D discussion ......................................... ................. 105

7 CONCLUSION AND FUTHER WORK............................... 115

7.1 Conclusion .... .................................... ................ 115
7.2 Future Studies ............. ................ ............... ......... 116
7.2.1 Evaluation of 3-D PAT in OA Detection with Small Clinical Samples.. 116
7.2.2 In-vivo Detection of Rheumatoid Arthritis in the Hand...................... 116
7.2.3 Recover Acoustic Property Together with Optical Properties by PAT
A approach ............. .. ..... ................. .... ....... ............... 117

LIST OF REFERENCES .................. ....... ......... .. ..... ............... 121

BIOGRAPHICAL SKETCH ................ ................... ........ 128









LIST OF TABLES

Table page

3-1 Phantoms used in the pre-investigation of 3-D PAT system............................. 56

5-1 Averaged absorption coefficient and size of joint tissues for 2 OA and 4
norm al joints ............. ...... ..................... ..... ..........99









LIST OF FIGURES


Figure page

1-1 Schematic of photoacoustic effect ............... ......... ....... ..... ............ 23

2-1 The reconstructed images for simulated finger joint (a) coronal section slice
along X=0mm (b) cross section slice along Z=15mm. ............................. 42

2-2 Geometry of the absorbed energy interpolation at fine node i in 2-D case (a)
and in 3-D case (b), and the calculation of the Jacobian matrix element in the
effective integration zone (c). ............................. ............................ ....... 42

2-3 Slice along the axis of the cylindrical target from reconstructed 3-D images
with full volume reconstruction (a), and partial reconstruction in region of
interest (b). .............. .................................... ............ ........... ..... 43

2-4 Reconstructed 2-D images of a circular target in a phantom experiment with
full volume reconstruction (a), and partial reconstruction in region of interest
(b ). ......................................................................... ................. ..... ...... 4 3

2-5 Schematic of parallel computing with a combination of both shared memory
and distributed m em ory. ............... ................................. ............... 44

2-6 Flow chart of the high performance photoacoustic tomography based on
parallel computing technique and finite element method................ ......... 45

2-7 The reconstructed images for simulated blood vessels (a) and crossed hairs
in phantom experiment (b)............. .. ........ .... ........... .... .......... 46

3-1 Schematic of two potential cylindrical scanning geometries of finger joints
imaging in 3-D PAT. ............ .......................................... 57

3-2 A slice through the axis of the reconstructed finger joint scanned by a
cylindrical geometry shown in Figure 3-1(b) ...... .... ....................................... 57

3-3 Schematic of the three-dimensional photoacoustic imaging system................. 58

3-4 Schematic of the scanning arc/path in the XY plane .................... ............... 58

3-5 Reconstructed absorbed energy density images in the coronal (z=0mm),
sagittal (x=0mm) and cross (y=4mm) sections for a joint-like phantom with
1mm thick cartilage in cylindrical (a) and spherical (b) scanning geometries,
and with 2mm thick cartilage in spherical scanning geometry (c) .................. ... 59

3-6 Reconstructed absorption coefficient images (coronal section) for the
phantom with 1mm (a) and 2mm (b) thick cartilage....... ........ .... .......... 60









3-7 The recovered absorption coefficient distributions along a transect (x=0mm)
for the images shown in Figure 3-6. ... ................................................... 60

3-8 Reconstructed absorbed energy density images in the coronal (z=0mm),
sagittal (x=0mm) and cross (y=4mm) sections for a joint-like phantom where
the absorption coefficient of cartilage was 0.015 (a), 0.025 (b), 0.03 (c) and
0.04 (d) m m -1 ................. .. ................................. ............. ... 61

3-9 Reconstructed absorption coefficient images (coronal section) for the
phantom with an absorption coefficient of cartilage of 0.015 (a), 0.025 (b),
0.03 (c) and 0.04 m m -1 (d) ........................................................ ... ........... 62

3-10 The recovered absorption coefficient distributions along a transect (x=0mm)
for the images shown in Figure 3-9. ... ................................................... 62

4-1 Schematic of the 3-D spherical scanning PAT system for finger joint imaging... 72

4-2 Test of lighting beam with a cylindrical phantom (a), a line phantom (b), and a
p e n c il le a d (c ). .............. .......... ..................................................... 7 2

4-3 Photograph of the arrangement for hand and finger Placement......................... 73

4-4 Schematic of the spherical scanning geometry (a), circular locus L1 (b) and
L 2 (c ). ............ ......... ... ......... .................. ................................ 7 3

4-5 Reconstruction mesh and possible model mismatch of 2-D circular scanning
(a) and 3-D spherical scanning (b). ........... ................................................ 74

4-6 Bubble level gauge (a) and homemade marker cap (b) for calibration of
mechanical errors in the photoacoustic spherical scanner ............................. 74

4-7 A selected coronal section from the reconstructed 3-D image (absorbed
energy density) at Z=-5mm (a), and MRI (coronal section) from a similar joint
(b ). ........ .... ........ ................................................................. 7 5

4-8 A selected coronal section from the recovered absorption coefficient image at
Z=-5mm (a), and its profile along the cut line X=-2mm (b). ............................ 75

4-9 A selected cross section from the reconstructed 3-D image (absorbed energy
density) at Y=7mm (a), and MRI (cross section) from a similar joint (b)............. 76

4-10 A selected cross section from the recovered absorption coefficient image at
Y=7mm (a), and its profile along cut lines 1 (b) and 2 (c), respectively. ............. 76

4-11 Cross section images (absorbed energy density) at Y=-2.5mm (a), Y=0mm
(b), Y=3mm (c), Y=5mm (d), Y=7mm (e) and Y=9mm (f)................................ 77









4-12 Coronal section images (absorbed energy density) at Z=-5mm (a), Z=-3.5mm
(b), Z=1 mm (c) and Z=4.5m m (d). ............................. ................... ............... 78

5-1 Photograph of the transducer array/finger interface of the three-dimensional
photoacoustic tomography system used for the current study......................... 92

5-2 Geometry of the detector array....... ............ ........... ... ............... 92

5-3 An example of possible mismatches between the center of the detector array
(marked as 0) and the rotary center of the rotary stage R2 (marked as 0')....... 93

5-4 Alignment variations (a) and spectrum differences in impulse responses of
four selected transducer elements in the ultrasound detecting array.................. 93

5-5 Block diagram of the multiple-channel 3-D PAT system ............. ............... 94

5-6 The central values of the 16 channels in PCIAD1650 card (a) and Labview
control interface (b) in multiple-channel 3-D PAT system............................. 95

5-7 Reconstructed optical absorption coefficient images in the coronal
(z=0.0mm) and sagittal (x=0mm) sections as well as in 3-D view for a joint-
like phantom where the absorption coefficient of "cartilage" was 0.01 (a) and
0.03mm-1 (b). (c) shows the recovered absorption coefficient distributions
along a transect crossing the "bones" and "cartilage" (x=0.5mm) for the
images in coronal section ............................... ........ ........................ 96

5-8 Photographic schematic of the coronal (a) and sagittal (b) sections for the
finger joint im aging. ................. ........... ..... ............. ........ ........ 97

5-9 Recovered images at selected coronal (a) / sagittal (b) planes for case 1
(H healthy) ......................................................... ............ .... ...... 97

5-10 Recovered images at selected coronal (a) / sagittal (b) planes for case 5
(O A ) ............... .......................... ................. .................................. 9 8

6-1 Absorption spectra of oxy-hemoglobin (HbO2), deoxy-hemoglobin (Hb) and
w ater (H20 ) in biological tissues ................... .......... ....... ...... ............... 110

6-2 Homemade 3-D positioning ruler. .................................. ......... ............... 110

6-3 Recovered optical absorption maps in coronal section (a) and sagittal section
(b) of an in-vivo finger joint for light source in six wavelengths...................... 111

6-4 Concentrations of oxy-hemoglobin (a) and deoxy-hemoglobin (b) in /M-, and
the water content (c) at coronal section (z=-8mm) and sagittal section (x=-
1 mm) of examined finger joint from an OA patient. ......................................... 112









6-5 Concentrations of oxy-hemoglobin (a) and deoxy-hemoglobin (b) in /M, and
the water content (c) at coronal section (z=2mm) and sagittal section (x=-
0.5mm) of examined finger joint from a healthy subject ................................. 113

6-6 Average Concentrations of oxy-hemoglobin(HbO2), total hemoglobin (Total-
Hb) and the oxygen saturation (SO2) in finger joint bones (a) and in joint
c a v ity (b ) ..................................................................................... 1 1 4

7-1 Reconstructed images at selected sagittal (a) / coronal (b) planes for an in-
vivo P IP joint............ .. ................. .................................................. 119

7-2 Simulation geometry to recover ultrasound speed with PAT method. (a) is the
simulation geometry of optical contrasts; (b) is the simulation geometry of
ultrasound contrasts; ............. .......... ................ 119

7-3 Reconstructed optical contrast images (A) and ultrasound speed images (B)
at 1st iteration (a), 3rd iteration (b) and 5th iteration (c)................. ................. 120









LIST OF ABBREVIATIONS

2-D Two-dimensional

3-D Three-dimensional

CT Computed tomography

DOT Diffuse Optical Tomography

DIP Distal Interphalangeal

FE Finite Element

MRI Magnetic Resonance Imaging

NIR Near-infrared

OA Osteoarthritis

PAT Photoacoustic Tomography

PIP Proximal Interphalangeal

qPAT Quantitative Photoacoustic Tomography

RA Rheumatoid Arthritis

US Ultrasonography









Abstract of Dissertation Presented to the Graduate School
of the University of Florida in Partial Fulfillment of the
Requirements for the Degree of Doctor of Philosophy

THREE-DIMENSIONAL PHOTOACOUSTIC TOMOGRAPHY AND ITS APPLICATION
TO DETECTION OF JOINT DISEASES IN THE HAND

By

Yao Sun

August 2010

Chair: Huabei Jiang
Major: Biomedical Engineering

This thesis research presents the study of three-dimensional (3-D) photoacoustic

tomography (PAT) and its application for the first time to detection of osteoarthritis (OA)

in the hand. PAT is an emerging non-ionizing, non-invasive imaging modality that can

visualize high optical contrast in biological tissues with high ultrasound resolution.

Compared with other imaging modalities that have been conventionally used or recently

investigated to visualize the structural abnormalities in the finger joints with OA, the 3-D

PAT approach studied in this dissertation not only provides high-resolution anatomical

structures, but also offers quantitative tissue optical property as well as physiological /

functional information including concentrations of oxy-hemoglobin (HbO2), deoxy-

hemoglobin (Hb) and water (H20) that could be used to detect OA in an early stage.

In this study, a 3-D high performance PAT reconstruction algorithm is developed

based on parallel computing technique and finite element method. The optimal detector

performance and scanning geometry for 3-D photoacoustic imaging of the finger joints

are investigated with considerable phantom experiments. The results of phantom

experiments show that a spherical scanning geometry appears to provide improved

spatial solution over a cylindrical scanning geometry, and that 1mm thick "cartilage" can









be accurately differentiated from the "bones" with a 1 MHz transducer in a spherical

scanning geometry. In addition, the absorption coefficient of the "cartilage" can be

effectively recovered when this optical property varied from 0.015mm-1 to 0.04mm-1. A

3-D PAT system in a spherical scanning geometry has been constructed and optimized

for in-vivo examination of human finger joints. A distal interphalangeal (DIP) joint from a

female healthy subject was photoacoustically examined by our 3-D PAT system, and

major anatomical structures of the examined DIP joint along with the side arteries were

clearly reconstructed in high quality, where joint space was well differentiated from

surrounding finger phalanges. The performance of the 3-D PAT system was further

improved with an ultrasound detection array composed of eight 1 MHz transducers and

a 16-channel pulse/receive board. The 8-channel 3-D PAT system was carefully

calibrated with controlled tissue phantom experiments, and is capable of completing a

finger joint examination within 5 minutes (at single optical wavelength).

The 3-D PAT reconstruction algorithm and scanning system developed has also

been applied to a pilot clinical study aiming to test the possibility of detecting OA in the

hand joints using 3-D PAT. In this pilot clinical study, seven subjects (two OA patients

and five healthy controls) were enrolled and photoacoustically examined. The image

quality of the reconstructed finger joints was greatly improved with the 8-channel 3-D

PAT system, and apparent differences, in both the reconstructed size of the joint space

and the absorption coefficient of the joint cavity, has been observed between the OA

and normal joints. The successful results obtained suggest the possibility of 3-D PAT as

a potential clinical tool for early detection of OA in the finger joints. Major chromophore

concentrations (Hb02 and Hb) of in-vivo finger joints have also been quantitatively









imaged using multispectral 3-D PAT approach with six optical wavelengths from 730nm

to 880nm. The multispectral results obtained further confirmed that the 3D PAT

approach implemented in this thesis research is able to differentiate OA from normal

joints.

While we target the detection of OA as a testing base for validating the single- and

multi-spectral 3-D PAT approaches developed in this thesis research, many aspects of

our work are fundamental to imaging in general. For example, the 3-D PAT approaches

implemented are applicable to other biomedical problems such as rheumatoid arthritis

(RA) detection and functional brain imaging.









CHAPTER 1
BACKGROUND AND SIGNIFICANCE

1.1 Introduction to Photoacoustic Tomography

Photoacoustic tomography (PAT), also referred as optoacoustic tomography

(OAT), is an emerging non-ionizing and non-invasive imaging modality in the field of

biomedical imaging based on the physical principal of photoacoustic effect, where light

energy is converted into acoustic energy due to optical absorption and localized thermal

expansion in biological tissues as shown in the schematic of photoacoustic effect in

Figure 1-1. In photoacoustic tomography in biomedicine, pulsed laser beam with

duration in nanoseconds is delivered into biological tissues, and some of the delivered

energy will be absorbed and converted into heat, leading to transient thermoelastic

expansion and thus wideband ultrasonic emission. The generated ultrasonic waves

(~MHz) are then detected by ultrasonic transducers to form images.

As a hybridized imaging technique between optical imaging modality and

ultrasound imaging modality, PAT combines both high optical contrast and high

ultrasound resolution in a single modality. Studies have shown that optical absorption

contrast between tumor and normal tissues in the breast can be as high as 3:1 in near-

infrared region due to the significantly increased vascularity in the tumors.1-3 Significant

absorption contrast between diseased and normal joints has also been observed in the

hand with osteoarthritis (OA) and rheumatoid arthritis (RA).4-7 For example, for an

osteoarthritic joint, the ratio of its cartilage absorption coefficient to that of the

associated bone is increased by 40% relative to the healthy joints. Unfortunately, the

high optical contrast in biological tissues was poorly imaged in the pure optical imaging

methodologies, since optical scattering in soft tissues degrades spatial resolution









significantly with depth. For example, for imaging depth greater than ~1mm, the spatial

resolution provided by optical diffusion tomography is no better than 3mm. With PAT

technique, the high optical absorption contrast in biological tissues is recovered by

measuring the megahertz acoustic pressure wave with ultrasound detectors. Since

ultrasound scattering is two to three orders of magnitude weaker than optical scattering

in biological tissues, photoacoustic tomography is capable of providing significantly

improved yet depth-independent spatial resolution over all-optical imaging techniques

such as diffusion optical tomography (DOT).

Thus far, PAT has shown its potential to detect breast cancer, to assess vascular

and skin diseases, to monitor epilepsy in small animals, to image finger joints, to sound

out fluorescent proteins, and to evaluate exogenous contrast agents in molecular

imaging.824 Recent studies further demonstrated the possibility for PAT to image hard

tissues such as bones and associated soft tissues.25-27 In these studies, rat tail joint or

cadaver human finger joints were imaged and gold nanoparticles were used as contrast

agent to help diagnose / monitor RA disease.

While simplified 2-D PAT detection geometries and reconstruction algorithms have

been used in most of the previous studies, light scattering and pressure wave

propagation in tissue is inherently 3-D, and 2-D approximation to a real 3-D problem (2-

D reconstruction model with 2-D detecting geometry surrounding 3-D targets) will

inevitably bring errors, blurring and missing structures in the reconstructed 2-D images.

Besides that, the structures of the examined biological tissues remain unidentified in

sections other than the detecting plane, which may be important in some biomedicine

cases. For example, a finger joint from a human subject is highly irregular in its









volumetric structure. Although it seemed that tissues around proximal interphalangeal

(PIP) and distal interphalangeal (DIP) joints were well differentiated in reported 2-D

cross-sectional images of a photoacoustically examined cadaver human finger 25, some

tissue structures shown in the obtained 2-D cross sectional images were not in

agreement with those from histological photographs due to the collapse of the signals

from the entire 3-D joint structures into the 2-D cross-sectional imaging planes. While

RA disease can be diagnosed from the cross section images alone, it is known that the

coronal and sagittal section images of the finger joints are necessary for OA diagnosis

which is impossible for 2-D PAT. Thus, full volumetric image reconstruction is essential

to capture the joint tissues completely and accurately for diagnosis of both OA and RA.

Although studies on conventional PAT have focused on obtaining high resolution

anatomical structure of biological tissues by imaging absorbed optical energy density

(i.e., the product of the intrinsic absorption coefficient and the external optical fluence or

photon density), PAT is also capable of providing the intrinsic absorption coefficient

distribution when a light transport model is combine with conventional PAT method.2834

Moreover, physiological / functional information such as oxygen saturation or the

concentration of hemoglobin, to which optical absorption is very sensitive, can also be

effectively recovered in high resolution when multispectral light is used.35-42 These

intrinsic tissue properties may be important for accurate diagnostic decision-making in

early stage of diseases.

1.2 Overview of Osteoarthritis

Osteoarthritis (OA) is the most common degenerative joint disease, affecting tens

of millions of Americans and involves over 500,000 individuals for joint replacement

annually in the United States.4344









Primary OA commonly affects small joints in the hand, including the distal

interphalangeal (DIP) joints and the proximal interphalangeal (PIP), and large weight-

bearing joints in the hips and knees. Although knee and hip OA bear most responsibility

for the burden of OA, hand OA may be a marker of a systemic predisposition toward

OA; incidence study of OA further revealed that the second DIP joint is most frequently

affected (57% in women, 37% in men) among all the joints in the hand.45

Pain, stiffness, tenderness, joint enlargement and limitation of motion are the

typical symptoms and signs of OA, and the pathological features of OA include erosion

of articular cartilage and associated bony changes.46 Synovial effusion and associated

enhanced blood vessel growth may be occasionally observed as well in OA joints,

although it may not be as severe as that in joints affected by rheumatoid arthritis (RA).

Thus far, clinical examination remains the principal approach to OA diagnosis,

relying on symptoms and signs OA patients suffered as well as experience of arthritis

physicians. Imaging techniques, including standard x-ray radiography, computed

tomography (CT), ultrasonography and magnetic resonance imaging (MRI), have been

conventionally used or investigated to visualize the physiologic/anatomic abnormalities

in the joint cavity when OA has been established.47-54 However, all these available

imaging techniques are used only as supplemental methods in OA diagnosis, especially

in hand OA, either because they are insensitive to the abnormal changes in the OA joint

cavity or because they are too costly to be used as a routine examination method. For

example, plain radiographs are able to visualize structural abnormality in the bone

compartment (joint space narrowing and osteophyte formation) in joint cavity with high

spatial resolution when OA has been well established, however they are insensitive to









changes in soft tissues (cartilage, synovial fluid, etc) and therefore incapable of

capturing the primary features when OA is establishing in earlier stage.

As such, sensitive and affordable imaging methods are urgently needed for the

detection of OA, especially in early stages. Moreover, the progress of effective imaging

methods in OA diagnosis may at the same time accelerate the advancement of medical

therapies for OA, which are currently effective to relieve OA symptoms or prevent the

worsening of OA to a certain degree and yet of limited effectiveness in modifying OA.

Up to now, no disease-modifying OA drugs (DMOADs) has been approved. The

development of the DMOADs may benefit greatly from the progress of imaging methods

in OA diagnosis, where sensitive imaging techniques can serve as surrogate markers

for clinically meaningful outcomes to economically and efficiently validate the candidate

DMOADs.43

New imaging techniques based on near-infrared (NIR) light, including pure-optical

imaging techniques and photoacoustic tomography (PAT), have been recently studied

to image finger joints and to effectively detect joint diseases. 4-7,25, 55 The high "color"

contrast (absorption coefficient) provided by diseased biological tissues (as high as 3:1

between tumor and normal tissues; > 40% increase in absorption for osteoarthritic

cartilage compared to normal cartilage), has made NIR based imaging techniques a

promising modality for early disease detection.2-3 While optical imaging techniques are

able to detect the highly sensitive optical absorption / scattering abnormalities

associated with soft tissues (cartilage, synovial fluid, etc) in OA and RA joints, its spatial

resolution is relatively low (about 3~5mm). Compared to all-optical imaging techniques,

PAT is able to visualize the same optical absorption contrast with significantly improved









spatial resolution (0.5mm or better, adjustable with ultrasound frequency) for deep-

tissue imaging. For example, tendon, aponeurosis, volar plate, subcutaneous tissue,

phalanx and other tissues around DIP joints were well differentiated in high resolution

when a cadaver human finger was photoacoustically scanned in cross section with a

10MHz transducer.

In this study, we develop a full 3-D PAT approach (high performance 3-D PAT

reconstruction algorithm and a 3-D PAT system in spherical scanning geometry), aiming

to study finger joints imaging and detect osteoarthritis disease in the hand. The high

performance 3-D PAT reconstruction algorithm is based on parallel computing

technique and the finite element method, which is validated phantom experiments. The

feasibility, the optimal detector performance, and the optimal scanning geometry for 3-D

photoacoustic finger joints imaging were pre-investigated with a series of finger-joint

phantom experiments.56 A 3-D PAT system in a spherical scanning geometry has then

been constructed and optimized for in-vivo examination of human finger joints, and the

distal interphalangeal (DIP) joint from a human volunteer was imaged.57 The

performance of the 3-D PAT system has been further improved with an ultrasound

detection array composed of eight 1 MHz PZT transducers and a 16-channel

pulse/receive PCI board. The 3-D PAT system with eight detecting channels was

carefully calibrated with controlled tissue phantom experiments, and is capable of

completing a finger joint examination within 5 minutes (at single optical wavelength).

Seven subjects (two OA patients and five healthy controls) enrolled in our study from

2008 to 2009 were photoacoustically examined with our 8-channel 3-D PAT system.58

Major chromophore concentrations (HbO2, Hb) of in-vivo finger joints were further









quantitatively imaged by using multispectral 3-D PAT approach with six optical

wavelengths from 730nm to 880nm.59












ncdelt Fulsed{NIF Liht


f/ '^^. \ Generated
.(- Ultrasonic Pulses

Absorber "A


Figure 1-1. Schematic of photoacoustic effect.









CHAPTER 2
RECONSTRCUTION ALGORITHM IN PAT

Thus far, several algorithms have been implemented to effectively reconstruct

photoacoustic images from measured acoustic waves, such as backprojection, Fourier

transform, P-transform, k-wave method and statistical approaches.60-64 While others rely

on analytical solutions to the photoacoustic wave equation in a regularly shaped

imaging domain and assume acoustical homogeneity in biological media, finite element

based algorithm seems provide unrivaled advantage to accommodate tissue

heterogeneity and geometric irregularity as well as allow complex boundary conditions

and source representations.65-67 Moreover, acoustic property in the biological tissues is

able to be recovered simultaneously with optical property by finite element based

algorithm, when acoustic heterogeneity is taken in consideration in photoacoustic wave

equation. In the following section, we will introduce the finite element based algorithm in

photoacoustic tomography in detail.

2.1 Finite Element Based Reconstruction Algorithm in PAT

For pulse mode PA generation propagating in soft tissues at room temperature,

the thermal diffusion and electrostrictive effects can be ignored. Therefore, only

considering the thermal expansion mechanism, the acoustic field by the PA generation

in tissue is described by the following wave equation:

V 1 21 02
V 2 p ,t)= H(F, t) (2.1)
c p 2 at

Where c is the speed of acoustic wave in the tissue; p(r, t) is the acoustic pressure; ,

is the thermal expansion coefficient; Cp is the specific heat at constant pressure; and









H(F,t) is the heating function defined as thermal energy per time and volume absorbed

in the tissue.

Assuming the laser light is a delta function in time and uniformly irradiated on the

surface of the tissue, and the heating function H(F,t) can be written as the product of

spatial absorption function D(7) and the temporal illumination as the form of delta

function S(t), we can write equation (2.1) as


V 2 2 p(,t)= = [(F)3(t)] (2.2)


Taking the Fourier transform of the above equation (2.2) in the temporal domain,

we have the following Helmholtz equation:

V2P(F, c)+k2P(F, ) =ik C W)) (2.3)
Cp

Where k is the acoustic wave number defined ask= r/c, = 2zf.

Although the acoustic speed c can be considered as a constant if we assume the

biological tissues are acoustically homogenous, it is indeed a variable quality in realistic

biological tissues. For variable acoustic speed in biological tissues, we can re-write

equation (2.3) as the following equation, in a reference of a constant acoustic speed in a

coupling medium

V2P(F, 0) + k2(1 + A)P(, )) = ik ( (2.4)
C
p

Where A= c/c2()-1

Choose yj as the set of basis functions, and then the weighted weak form of

equation (2.4) can be stated as:









[V2 +k,2(1+ A)P] d Yfik c/p I() dV= 0 (2.5)
P

Expanding P(rF,r) and (7F) as the sum of complex coefficients multiplied by the

basis functions: P = p,/,, = y A k= N A the finite element discretization

of equation (2.5) can be written as the following linear equations:

[A]{p} = {b (2.6)


Where A =J (VV -Vyi/)dV- vk02vIVdV- vkO kkAydV- (u, +y- )i)ds


b = -ik J dV
Cp 1-i

p = { P2,**, PNT

And the following absorbing boundary conditions have been applied:

82P
VP.ih == P + y2 (2.7)


Whrl -iko -3/2p+i3/8kop2 -i/2kop2
Where = -; /= and ( is the angular coordinate at a
1-i/kop 1-i/kop

radial position .

2.1.1 Forward Modeling in Acoustically Homogenous Media

Although the acoustic speed c is indeed a variable quality in realistic biological

tissues, we may assume that the biological tissues are acoustically homogenous for

simplicity in some cases and consider the acoustic speed as a constant in interested

biological tissues during the process of image reconstruction. In these cases, equations

(2.6) can be simplified as the following equations









[A]{p}={b} (2.8)


Where A = (V V, ) k2 k0( ) (17 +Y ,ds



CoP 1
pp i=1

P =PI, P2,**, PN T

In photoacoustic tomography with the assumption of acoustical homogeneity in

biological tissues, equation (2.8) serves as the forward model to compute the spatial

variation of acoustic pressure wave (external observable) based on a given optical

property distribution, or to compute the updated acoustic pressure wave in space based

on a initial guess or iterative value of an optical property distribution in image

reconstruction.

2.1.2 Inverse Modeling and Nonlinear Optimization in PAT Reconstruction

In the inverse problem of PAT reconstruction, the optical absorption distribution

(ultrasound source term in ultrasound propagation equation) was estimated given the

governing photoacoustic equations, boundary conditions, and measurements of the

pressure wave on the detecting locations (external observable). Gauss-Newton iterative

scheme combined with Marquardt and Tikhonov regularizations has been pursued in

our research to obtain stable inverse solution to the PA equation. This approach uses

the hybrid regularizations-based method to update an initial (guess) optical absorption

distribution iteratively to minimize an objective function composed of a weighted sum of

the squared difference between computed and measured data for all frequencies.

Unlike forward problem in PAT where the stiffness matrix A is well-conditioned and

the solution is well-defined, inverse problem in PAT has to estimate the spatial tissue








properties with limited observables on the boundary and most commonly results in ill-

posed situation, yielding meaningless solutions in image reconstruction. In ill-posed

situation, solutions in inverse problem may be non-unique, non-existent or unstable (a

tiny perturbation in the measurement results in large differences in the solution).68 To

improve the conditioning of the problem and enable a numerical solution in inverse

modeling, Tikhonov regularization was used to find an optimized solution in inverse

problem by minimizing the following modified objective function

F(z)= pc(,)- -+, 2 (2.9)

Where the second term is an added penalizing factor with a weighting A in a Euclidean

distance referenced to *, and f* is set equal to the parameter distribution at the

previous iteration ,' 'for simplicity in Levenberg-Marquardt algorithm.

At the ith iteration, we assume that the Taylor's expansion of the objective function

F(2) around j(') can be expressed as


F() = F(2'))+ VF(('))( (')+ HF (2) ())2 +... (2.10)
2

Where VF(f"') and HF (()) are the first and second order derivative

At the minimum of the objective functionF(f), we have stationary point where the

following identity is satisfied

dF(j)
=d ) 0 (2.11)
d2

Substituting (2.10) into (2.11) and ignoring the derivatives higher than second

order, we have the following Gauss-Newton iterative scheme

P+1) = P -[ HF(2())]1 VF(')) (2.12)









Consider the regularized objective function given in equation (2.9), where f* = ),

we have the first and second order derivative of the objective function F(f)


VF(2)= 2 P(-) dp ) +2A(f -2 ))
di


HF(2)= 2dp( + .( 22 (2.13)
d- di

Combining the above first and second order derivative of the objective

function F(2) with the Gauss-Newton iterative scheme (2.12), we have the iterative

equation

(JTJ + 2AI)A = JT(pm _pC) (2.14)

Where the Jacobian matrix J is defined over total M detecting locations as

ap, ap, aOp
08, 082 ON
ap2 ap2 2ap
J = a, O2 aN


aPM aPM aPM
a8, a82 aN

PC {pcP2c. ,PMc T

Pm =PIm,p2m, PMm T

Aj= () = {Ao,, A 2,...,AoN }T (2.15)

Since p is a linear function of Ofrom the linear equation set (2.8), we can calculate

Jjk by applying partial differentiation on equation (2.8) with respect to k,









a(Ap) (B( B -k COP VdS (2.16)
S(Pk (P Cp '

A, and B are independent on
variable if j # k, therefore we have

AuJ k =Bk = B=k (2.17)

In photoacoustic tomography, the final equation (2.14) with Jacobian matrix

determined in equation (2.17) is the core of the inverse computation to update the

optical property distribution from an initial estimate. In each iterative step, it compares

the measured with the solved observables, which is the acoustic pressure wave for

different frequencies obtained via the forward model (2.8) and current estimate of the

optical property distribution, and solves the inverse problem by nonlinear optimization

equation (2.14) to obtain a new estimate of the optical property distribution. This

procedure is iteratively repeated until optimal least-squares fit of computed and

measured acoustic data is reached.

2.1.3 Regularization in PAT Reconstruction

The purpose of introducing the Tikhonov regularization on the Hessian matrix H is

to ensure the diagonal dominance and enable a stable numerical solution in inverse

procedure. However, the optimal weighting factor A is problem-specific and sometime

quite difficult to determine. An effective approach for choosing optimal weighting

factorA is to iteratively modify its value based on the trace of the Hessian matrix and the

relative least-square error erel at each iterative step: 69

2 = a x ere x trace([J'J]) (2.18)

Where is an empirical coefficient. We normally usea = 0.5 in our PAT problem.









Since the weighting factor A relies on er, which iteratively decreases, the value of

A and the regularization effect will be reduced at each iterative step. The variation of

A in each step allows more accurate image to be revealed as the iterations progress,

while a blurred image is initially reconstructed with large A at the start of the iterations.

The empirical coefficienta is carefully chosen so that the image can be accurately

reconstructed with iterations while stability to a numerical solution is ensured at the

same time.

Although weighting factorA defined as identify (2.18) in Tikhonov regularization

enable the stable numerical solution in inverse problem, the Gauss-Newton iterative

scheme may experience problem in iterative convergence since the approximation

(2.10) of Gauss-Newton iterative scheme are iteratively valid only near the minimum

(actual solution). If the initial estimates are far from the actual solution, Gauss-Newton

iterative scheme may have difficulty in convergence.

Other than Gauss-Newton iterative scheme, gradient descent, also referred as

steepest descent, relies barely on the initial estimates although it may exhibit oscillation

around the minimum and the convergent process is slow,68 since the first order

derivative of objective function is used in gradient descent instead of second order

derivative used in Gauss-Newton iterative scheme. The real-valued objective function

F(2) decreases fastest from a given estimate along the direction of the negative

gradient of F(f) at the given estimate, and hence the iterative algorithm for gradient

descent is given by

(+1) = (+1) p'VF( (')) (2.19)

Where p(' is a positive value used as iterative step size.









Considering VF(i"') given by equation (2.13) from Tikhonov regularized objective

function (2.9), we have the following iterative algorithm for gradient descent

( + A)Aj = J' (pm p) (2.20)
2p

Comparing the gradient descent method given in (2.20) and the Gauss-Newton

iterative scheme given in (2.14), we may see some similarity between them. With a

hybrid technique proposed in Levenberg-Marquardt method, we can use a steering

factor A' to switch between the Gauss-Newton iterative scheme and gradient descent

method. The final update equation can be written as

(JTJ + I + I)A = J (pm -pC) (2.21)

Where A is a weighting factor used in Tikhonov regularization to relieve the ill-posedness

of the inverse problem; A' is a steering factor introduced in Levenberg-Marquardt

method to balance the convergent ability and the convergent speed.

When the steering factor A' is very big, the term jTj can be ignored and equation

(2.21) shows the iterative behavior of gradient descent method shown in equation

(2.20). When is small, then equation (2.21) reduces to equation (2.14) and exhibits

the behavior of Gauss-Newton iterative scheme. So practically, the value of A' is

chosen iteratively by the following rules:

1) The steering factorA' is initially set with a big value so that the initial estimates

can be chosen with less caution with the gradient descent behavior shown in (2.21).

2) In each iterative step, the steering factor A' reduces by a given ratio (such as

10) to accelerate the convergence process if the objective function decreases.

Otherwise, the steering factor A' is increased by a given ratio.









2.1.4 Simulation in PAT Reconstruction

The reconstruction algorithm developed for 3-D PAT is tested with simulated finger

joints, where a cylindrically simulated finger with a 2.5mm thickness "cartilage"

sandwiched between two bones is buried in background. The absorbed energy contrast

between the simulated cartilage and the bones is assumed as 3:1, and the two bones

are 10mm in diameter.

The mesh for algorithm testing is a cylindrical mesh with total 14784 tetrahedral

elements and 2907 nodes, which is generated from a 2-D triangle mesh slice layer by

layer. Each layer of the total 17 layered cylindrical mesh has 308 elements and 171

nodes, with layer interval 1.25mm in the axial direction (Z axis). The coronal section and

the cross section of the reconstructed finger joint from simulated data are shown in

Figure 2-1, where the actual shape of the simulated finger joint is outlined with black

dash lines. From the reconstructed finger joint, we can see the bones and the cartilage

are recovered with great quality. The structure and the size of the bones and the

cartilage are in great agreement with the actual shape and size of the simulated bones

and the cartilage. And the absorbed energy contrast between the simulated cartilage

and the bones is also recovered with high accuracy. The results shown in this

simulation validate the reconstruction algorithm of 3-D PAT for further finger joint

imaging.

2.2 Computing Strategies in PAT Reconstruction

In finite element based PAT reconstruction, the major task is solving the forward

equations (2.8), the Jacobian matrix (2.17) and the inverse equations (2.14), which are

normally time / memory costly depending on the scales of the linear system of equation

(the total number N of finite element nodes). As the scale increase, the cost in








computational memory and time grows dramatically, which may go beyond the

computation limitation of a single PC. Herein, some efficient computing strategies have

been adopted to bring the high resolution PAT in full potential.

2.2.1 Adjoin Sensitivity Method

Normally only the pressure on the boundary nodes (total M nodes) is observable,

therefore, in the calculation of Jk, only j = det(s), s = 1, 2,...,M is useful for updating q in

the inverse equation (2.14). Considering the Jacobian matrix defined in equation (2.17),

we have

[J=sk [g ] [J, j ] [ s, A[ [ ]=[S, ] [Bk] (2.22)

Where [,,] is defined as 8,, =1 forj = det(s), otherwise 3,, = 0; and S, can be

calculated from [8 ] = [S,, ][A]

It is worth to note that Bk = 0 if node iand k belong to different element as

shown in equation (2.16). As a result, we only need sum up the contribution of the

nodes in elements contain nodek, even though the sum index i in S,, varies from 1

to N.

2.2.2 Dual Meshing Scheme

In finite element based reconstruction algorithm, spatial sampling rate (mesh

density) that is used to discretize the target region is critical and greatly determines the

computational costs and the resolution of the reconstructed images. For example, the

rapid spatial variation in the field of acoustic pressure wave (megahertz frequencies)

governed by the photoacoustic wave equation demands high mesh resolution

(typically A/10, A is the ultrasound wavelength) with finely discretized mesh in order to









maintain forward calculation accuracy of the acoustic pressure field in the media.

However a finely discretized mesh with increased spatial sampling rate, and hence the

total numbers of mesh nodes, will inevitably results in dramatically growing

computational cost, which almost grows squarely and cubically with the total mesh

number in the forward calculation of the acoustic pressure and the inverse image

reconstruction of unknown absorbed energy profiles (or absorption coefficient profiles),

respectively. At the same time, the unknown absorbed energy profiles to be

reconstructed in PAT are relatively uniform and hence representable with fewer mesh

nodes (degrees-of-freedom). To balance the accuracy required for the calculation of the

acoustic pressure wave and the cost of computation, dual-mesh scheme seems an

optimal solution in the finite element based PAT algorithm, where a dense mesh is used

in the forward calculation of the acoustic pressure and a relatively coarse mesh is used

in the inverse image reconstruction.71

Implementation of the dual mesh scheme affects two components of the

reconstruction algorithm:7173 (1) the forward calculation of the acoustic pressure field in

the media at each iteration, where the updated absorbed energy profile to be used in

the forward calculation is defined on the coarse mesh while the forward calculation is

based on the fine mesh, and (2) calculation of the Jacobian matrix that is used to

update the estimates of the absorbed energy profile during the inverse image

reconstruction.

To calculate the acoustic pressure field in the media with the updated absorbed

energy defined on the coarse mesh, interpolation process is required to get the

absorbed energy profile over the fine mesh from the coarse mesh. Assume that the ith









fine mesh node is inside the coarse meshL, as shown in Figure 2-2(a) for a simplified

2-D case and Figure 2-2(b) for a 3-D case, the absorption energy value at the ith fine

mesh node can be interpolated by the following formula

NP
-= (X ) (2.23)
n=1

Where (o is the absorbed energy value at each node of the coarse meshL, N, is the

number of the element nodes. N = 3 for 2-D triangle mesh, and N = 4 for 3-D

tetrahedral mesh.

For a certain pair of coarse/fine mesh, the interpolation coefficient Y'I () is

always a constant, and therefore can be pre-calculated and stored in a mesh related file

before the reconstruction process, which saves the time of re-computing during each

iteration step.

With the dual mesh scheme, the elements of the Jacobian matrix will be calculated

by


A, = B'' = -ik cpy k'V dS (2.24)
a (p' CP

Where k is the node on the coarse mesh, ,k' is the basis function centered on node k in

this mesh, and the integration is still performed over the elements in the fine mesh.

Since yk' and y, are the basis functions defined over the coarse mesh and fine mesh,

respectively, the integration kernel k y is non zero only in the overlapping zone of the

coarse mesh the coarse nodek belongs to, and the fine mesh the fine node i belongs

to. Complications arise from the integration when the fine mesh spans more than one

coarse element. To simplifying this problem, we can generate the fine mesh from the









coarse mesh by splitting the coarse elements into fine elements so that each fine mesh

resides entirely within one coarse mesh element as shown in Figure 2-2(c).

2.2.3 Partial Reconstruction in Region of Interest

Unlike optical tomography which is an inverse field problem to quantitatively

reverse the optical scattering / absorption properties in the media by measuring the

transported light arriving optical detectors from light sources, photoacoustic tomography

is indeed an inverse source problem. In photoacoustic tomography, the goal is to trace

the ultrasound sources (optical absorbers excited by near-infrared pulse) by measuring

the ultrasound field on a boundary with pre-known ultrasound properties in the media.

The particular feature of this imaging technique gives us an inspiration that we may

confine our inverse reconstruction in the region of interest that may contain the

ultrasound sources, which could greatly reduce the computation cost in the inverse

procedure.

To confine our inverse reconstruction in the region of interest, a slight adjustment

on the inverse equation is required where inverse Equation (2.14) is rewritten as

(JT + 2A)I)A = JT(pm pc) (2.25)

Where the new Jacobian matrix J is defined as

ap, ap, ... ap,
a8, a 2 aK
Oap2 ap2 2a
J= ax, a82 aK


aPM aPM aPM
aX1, aX2 axK









The new Jacobian matrix J is therefore a first order derivative of the ultrasound

pressure at each detecting location to the unknown optical absorbed energy (the

intensities of ultrasound sources) at each finite element node within the region of

interest, which includes total K finite element nodes.

The partial reconstruction is tested by simulation, where a cylindrical target (10mm

in diameter with 20mm height) is embedded in a cylindrical background. The contrast

between the target and the background is 7:1. The section image along the axis of the

cylindrical target from the reconstructed 3-D images with full volume reconstruction is

shown in Figure 2-3(a), where total 2907 nodes in the whole volume are involved in the

inverse reconstruction. At the same time, the cylindrical target is reconstructed with

partial reconstruction in region of interest, which is a small cylindrical volume containing

the cylindrical target. The region of interest include total 801 nodes and is outlined with

dash black lines in the section image along the axis of the cylindrical target as shown in

Figure 2-3(b), which is a 2-D section image from the reconstructed 3-D images with

partial reconstruction in the region of interest. As observed from Figure 2-3(b), the

partial reconstruction in region of interest in our study obtained almost the same image

shown in Figure 2-3(a) with full volume reconstruction in term of structure and

quantitative values, while the computational cost in the partial reconstruction in region of

interest is only around 1/47 of the computational cost in full volume reconstruction.

The partial reconstruction is further validated with phantom experiment as shown

in Figure 2-4, where the full volume contains 1301 nodes and the region of interest

outlined with black dash line includes only 229 nodes. Compared with the reconstructed

images with full volume reconstruction (Figure 2-4(a)), images reconstructed with partial









reconstruction in region of interest maintain the same quality and quantitative results

(Figure 2-4(b)) while the computational cost in the partial reconstruction in region of

interest drops to around 1/32.

2.2.4 Parallel Computing Technique

In some cases where a fine mesh both in the forward calculation and inverse

reconstruction is required for high resolution PAT images, the computational cost may

go beyond the limitation of a single PC. For example, computational memory in tens of

gigabytes is required to store the matrix A in forward equations, the Jacobian matrix J

and the Hessian matrix JTJto image a blood vessel in subl00om spatial resolution. To

overcome the computational cost in high resolution PAT (as well as 3-D PAT, where the

problem scale grows cubically with the spatial resolution instead of squarely in 2-D

case), parallel computing technique has to be utilized. In a parallel computing scheme,

the memory request and the computation tasks for a computation assignment can be

spread on multiple processors as shown in Figure 2-5. In the following, we will give a

detailed description on the parallel computing based PAT algorithm.

In finite element based PAT algorithm, the forward calculation of acoustic field by

Equation (2.8) and the determination of Jacobian matrix J by Equation (2.17) and (2.22)

are based on the forward modeled matrix A at each frequency o ordered from o,

to oK with totalK frequency elements. The matrix A in forward equations is usually a

symmetric sparse matrix and stored in memory by banded storage or compressed

storage strategy, requiring less memory than the Jacobian matrix J and the Hessian

matrixJTJ, which is usually a full matrix. Herein, we spread the forward calculation of

acoustic field and the determination of Jacobian matrix on distributed processors (total









number is Q +1) by the frequency element o,, where the whole matrix A under certain

frequency elements is stored in the memory of the specifically assigned processor. As

shown in flow chart of the high performance photoacoustic tomography (Figure 2-6), the

matrix A at frequency is stored in the memory of a processor determined

byMod(i-1,L), whereL is the average number of frequency elements on each

processor (total frequency elements K divided by total processorsQ+1). In each

processor, the forward modeled acoustic field as well as the elements of Jacobian

matrix at the frequencies associated with the specified processor is calculated

independently, after which the Jacobian matrix J over the whole frequency range is

assembled and stored dispersedly in the memories of the distributed processors. In this

way of parallel computing and storage, the computation load and storage is evenly

assigned among the processors, and the computation assignment is maximally

parallelized since each processor can run the assigned task independently without

communication with other processors.

Jacobian matrix J and the Hessian matrix jTj (denoted asH) is full matrix,

requiring distributed storage over the processors. A wrapping storage over the columns

of the Jacobian matrix J is used to evenly store the Jacobian sub-matrix over the

processors and further efficiently calculate the Hessian matrixHwith minimally mutual

communication between the processors. The wrapping storage is described below.

} E CPU Myid, if Myid = Mod( -1,Q +1) (2.26)

Where { J, J ,... ,J ,J ,...,Ml) J J= 1,...,N.








With the wrapping storage of the Jacobian matrix J, the Hessian matrixH can be

calculated and stored by

{H} ={, i ... r {j CPU Myid, if Myid = Mod( -1, Q +1) (2.27)

Where {j H}= {HI,,...,H H j=1,...,N.

Because the Hessian matrix H is symmetric, we store only the upper triangle of the

Hessian matrix spread over the processors. In the calculation of elementH!y,

communication and data exchange may be involved among the distributed processors,

except that the vector {.j and {j are both stored in the memory of the same

processor. Since the Hessian matrix is stored dispersedly, the inverse reconstruction

by Equation (2.14) requires a parallel solving method, where parallelized Cholesky

decomposition is used in our high performance PAT.

With the high performance PAT, simulated blood vessels with different diameters

(0.6mm, 0.28mm and 0.14 mm) can be accurately reconstructed as shown Figure 2-

7(a), using a finite element mesh with 15200 triangle elements and 7721 nodes. The

0.14mm diameter vessel (red dash line in Figure 2-7(a)) is quantitatively recovered as

0.133mm in diameter based on FWHM analysis. Crossed hairs (~60pm in diameter) in a

phantom experiment can also be imaged in a 0.1mm resolution when a finite element

mesh with 14377 nodes is used, as shown in Figure 2-7(b).




















1C


I:

N
E'
E'i


fi-


25

2 5


n


-ra


1 (b)

(b)


i 1is


Figure 2-1. The reconstructed images for simulated finger joint (a) coronal section slice
along X=0mm (b) cross section slice along Z=15mm.


(a) (b)


Figure 2-2. Geometry of the absorbed energy interpolation at fine node i in 2-D case (a)
and in 3-D case (b), and the calculation of the Jacobian matrix element in the
effective integration zone (c).













rI


-I t


-ID a ID 13


Figure 2-3. Slice along the axis of the cylindrical target from reconstructed 3-D images
with full volume reconstruction (a), and partial reconstruction in region of
interest (b).


/
20 0 I


-10 0 10
(a)


-10 0 10
(b)


Figure 2-4. Reconstructed 2-D images of a circular target in a phantom experiment with
full volume reconstruction (a), and partial reconstruction in region of interest
(b).


6
4


r


I


- I f0












A Computation Assignment


UP P
CP 0P


UPUCP


U


UP P


Figure 2-5. Schematic of parallel computing with a combination of both shared memory
and distributed memory.


* 0 0


* 0 0













CPU Q


* 0


Solving [ S, ] by Solving [S, ] by

[1s, ]=[.1Mi] [ ]= 1 ]Mi 1
s=1 -M,i =1 N s =1 -M,i =1 N


Solving [S,] by

s 1 [ Mi ]= Ns.]
s=I~M,i=I-N


Figure 2-6. Flow chart of the high performance photoacoustic tomography based on
parallel computing technique and finite element method.

















2.5 2000
5 -
to
0 2 0 1000ooo

-2 1.5 -5 0
-10
1 -1000
-4 -2 0 2 4 -15 -10 0 10
X(mm) 10 0 10
(a)



Figure 2-7. The reconstructed images for simulated blood vessels (a) and crossed hairs
in phantom experiment (b).









CHAPTER 3
PRE-INVESTIGATION OF 3-D PAT SYSTEM WITH TISSUE PHANTOM STUDY

In this section, we aim to systematically evaluate the possibility of 3-D

photoacoustic imaging of the finger joints, based on a series of tissue-like phantom

experiments. The 3-D structures and absorption coefficient images of the joint-

mimicking phantoms are obtained using finite element based PAT reconstruction

algorithm and its enhancement.

Among different 3-D scanning geometries, cylindrical scanning geometries (shown

in Figure 3-1) adopted from DOT seem straightforward and are easiest to implement,

however they give poor resolution along the cylindrical axis due to the limited aperture

size of a transducer. For example, a joint-like phantom with cartilage thickness in 1mm

was scanned by a 1MHz transducer in a cylindrical scanning geometry shown in Figure

3-1(a), and the joint space is missing in the reconstructed images as shown in Figure 3-

2. A transducer with small aperture size can theoretically offer enhanced resolution;

however, such a transducer often gives poor signal-to-noise ratio (SNR) with added

cost. Compared to a cylindrical scanning geometry, a spherical scanning geometry with

a modest size of transducer aperture may provide an optimal solution for 3-D imaging

as the spatial resolution with such a scanning depends little on the aperture size of the

transducer, especially around the center of the imaging zone. In this study, phantom

experiments based on both cylindrical scanning geometry and spherical scanning

geometry are conducted for comparison, the results of which show that spherical

scanning geometry provides obvious spatial resolution improvement over cylindrical

scanning geometry and might be suitable for finger joint imaging. The spherical

scanning geometry is adopted in further phantom experiments to systematically









evaluate the differentiability in controlled phantoms with varied thickness and varied

optical absorption values in "cartilage".

3.1 Tissue Mimicking Phantom

To evaluate the reconstruction algorithm of 3-D PAT and the possibility of 3-D

photoacoustic imaging of the finger joints with phantom experiments, solid tissue

phantoms that well approximate optical characters of in-vivo finger joints are fabricated.

The solid tissue phantoms used here are made of Intralipid, India ink, distilled water and

agar powder, among which Intralipid acts as the optical scatter, India ink supplies color

contrast (optical absorber) and the agar (1%-2%) serves as the coagulator. The joint-

like tissue phantoms used in our experiments are optically well-characterized, in term of

the optical absorption/scattering coefficients. The highly purified agar powder (A-7049,

Sigma, USA) is dissolved in distilled water, and is liquefied after being boiled to the

melting temperature of 950C with a regular microwave oven. The solution has to be

stirred continuously in order to obtain good uniformity in fabricated phantom while

cooling down to 400C, when the solution must be poured into a mould and left there for

some time (normally about 10 minutes for 400mL solution) to reach a proper hardening

and stable optical properties. The optimal temperature for adding Intralipid and India ink

to the agar solution is around 600C, although it is non-critical between 800C and 400C.

The detailed fabrication procedure and the performance in controlled optical properties

of the tissue phantom have been fully described before.74

The two phantom "bones" fabricated in our study to mimick real finger bones had:

/u, = 0.07mm-1, /, =2.5mm-1 and a diameter of 6mm, and the "cartilage" between the

bones had: p = 0.015~0.04mm-, /, =1.5mm-1 and a thickness of 1 or 2mm, as shown









in Table 3-1. All the optical properties of the fabricated phantom are based on real finger

joints.

3.2 System Description and Experiments

The joint-like tissue phantom experiments in this study were conducted with a 3-D

PAT imaging system as shown in Figure 3-3, which is capable to scan the phantoms

along a cylindrical surface as well as along an equivalent spherical surface.

The Ti: Sapphire laser (LOTIS, Minsk, Belarus) generates a pulsed beam with a

pulse repetition rate of 10 Hz and a pulse width <10ns. The wavelength of the laser was

tunable from 600 to 950nm and was locked at 820nm in our experiments for deep

medium penetration. After the reflection mirror, the laser is then expanded by a lens and

a ground glass so that the phantom can be uniformly illuminated by an area source of

light. The laser energy used at the phantom surface was about 10mJ/cm2, which is far

below the safety standard of 22mJ/cm2.

To detect the acoustic field generated by the pulsed laser, 1 MHz transducer

(Valpey Fisher, Hopkinton, MA) was used. The transducer has a bandwidth of 80% at -

6dB and an aperture size of 6mm. To reduce the attenuation of the acoustic field, the

transducer and the phantom were all immersed in a water tank, in which the sound

speed was specified as 1495m/s.

During an experiment, the phantom was placed at the center of the rotary stage,

and the transducer was attached to the rotary stage via an arm having a length of

~20mm. To acquire the acoustic field along a cylindrically scanning surface, the

transducer was stepped along the Z axis as well as rotated along an arc facing the

center of the joint-like phantom in XY plane. In our experiments, the arc scanning in the

XY plane covered 3000 around the phantom as shown in Figure 3-4, allowing the









collection of signals at 50 locations with a step interval of 60; the linear scanning along

the Z axis permitted data collection at 650 positions along the cylindrical surface using

13 steps with an interval of 1mm.

For the spherical scanning mode, the arc scanning in the XY plane covered the

same 3000 around the phantom as the cylindrical scanning, and the phantom was also

rotated along the Y axis so that the transducer was equivalently rotated around the

phantom. The phantom was rotated 150 per step along the Y axis until a full spherical

surface was covered, which allowed the data collection at 600 scanning positions.

The detected acoustic signals were amplified and filtered by a Pulser/Receiver

5058PR (Panametrics, Waltham, MA). The bandwidth of 5058PR ranges from 10kHz to

10MHz, and five scales from 0.01 to 1.0MHz are selectable in high pass filtering. The

preamplifier integrated in 5058PR is in 30dB gain and the gain of the 5058PR itself is

40/60dB with attenuation tunable from OdB to 80dB in 1dB steps. The Labview

programming was used to control the entire data acquisition procedure.

3.3 Methods

Our finite element based 3-D PAT reconstruction algorithm and its enhancement

for extracting quantitative absorption includes two steps. The first is to obtain the 3-D

images of absorbed optical energy density through our 3-D PAT reconstruction

algorithm that is based on finite element solution to the photoacoustic wave equation in

frequency domain subject to the radiation or absorbing boundary conditions (BCs),

which is described in detail in Chapter 2. The second step is to recover the optical

absorption coefficient distribution using the photon diffusion model coupled with the

absorbed optical energy density images obtained from the first step.









In the first step, the 3-D PAT reconstruction algorithm uses the regularized Newton

iterative strategy to update an initial absorbed optical energy density distribution to

minimize an object function, which is composed of a weighted sum of the squared

difference between computed and measured acoustic data.

To recover the optical absorption coefficient from absorbed energy density, ,, the

photon diffusion equation as well as the Robin boundary conditions can be written in

consideration of D =,u ,

V.D(F)V(E(F)(F))- O() = -S(F) (3.1)
-D(i)V (E(F)(F)) n = E()ac()) (3.2)

WhereE(F)=1/ ,(F) D(F) is the diffusion coefficient, D(F) =/[3(U, (F)+ U; ())] and pu;(F) is

the reduced scattering coefficients, a is a boundary condition coefficient related to the

internal reflection at the boundary, and S(F) is the incident point or distributed source

term. For the inverse computation, the so-called Tikhonov-regularization sets up a

weighted term as well as a penalty term in order to minimize the squared differences

between computed and measured absorbed energy density values,

Min (F)- 2 + (F IL(F)[E(F)- Eo (F)]12 (3.3)

Where L is the regularization matrix or filter matrix, / is the regularization parameter,

so =N (oo....)T and o, =(v, _.... ,)T where l is the absorbed energy density

obtained from PAT, and 0C is the absorbed energy density computed from Equations.

(3.1) and (3.2), for i=1, 2..., N locations within the entire PAT reconstruction domain.

The initial estimate of absorption coefficient can be updated based on iterative Newton

method as follows with 3=1,

A(E) = (JTJ + + 1+LTL)- [J ((io ()] (3.4)









In addition to the usual Tikhonov regularization, the PAT image (absorbed optical

energy density) is used both as input data and as a priori structural information to

regularize the solution so that the ill-posedness associated with such inversion can be

reduced. The whole imaging zone is segmented by Amira (Indeed-Visual Concepts

GmbH, Berlin, Germany) into several sub-regions, and the segmented a priori spatial

information based on the PAT image can be incorporated into the iterative process

using Laplacian-type filter matrix, L ,

i if 1=j
L, = -1INN if i, cone region (3.5)
0 if i, c region

Where NN is the total node number within one region or tissue; the optical absorption

coefficient distribution is then reconstructed through the iterative procedures described

by equation (3.1) and (3.3), where a priori spatial information from the PAT image is

incorporated in through the matrix L.

3.4 Results and Discussion

The 3-D images were reconstructed using a finite element mesh of 41323

tetrahedral elements and 7519 nodes which required around 120 minutes per iteration

for a total of 10 iterations on a computation cluster with 10 processors. Each of the

processors worked at 2.2GHz with 2GB memory.

The comparison between the cylindrical and spherical scanning geometries can be

made from the images shown in Figure 3-5a and Figure 3-5b for a typical joint-like

phantom with 1mm-thick cartilage. Figure 3-5a display the coronal, sagittal and cross

sections of the 3-D images from the reconstructed absorbed energy density distributions

with the cylindrical scanning, while Figure 3-5b show the corresponding sectional

images with the spherical scanning geometry. To compare the reconstructed images









with the phantom, the actual shape of the phantom is outlined in Figure 3-5 (Dashed

rectangular or circle). As can be observed from Figure 3-5a, the phantom structures in

the coronal section (XY plane) are well reconstructed, while both the structures in the

sagittal and cross sections (YZ and XZ planes) are clearly distorted especially for the

cross section. Here we see that the structures along the circular plane (XY plane) can

be recovered accurately and that the structures are erroneously stretched along the

cylindrical axis (Z axis). The images presented in Figure 3-5b with the spherical

scanning show clear improvement over that using the cylindrical scanning. We see that

there is almost no distortion along the Z axis especially for the cross section. These

images with the spherical scanning are quantitatively accurate in terms of the recovered

size, shape and location of the objects.

Since the spherical scanning gives us the best image quality for 3-D joint imaging,

the remaining results were all obtained based on this scanning geometry. Figure 3-5c

presents the coronal, sagittal and cross sections of the reconstructed 3-D absorbed

energy density for a joint-like phantom with 2mm-thick cartilage. It is observed from

Figure 3-5b and 3c that the structures along each section are well detected and that the

thickness of cartilage for both 1mm and 2mm cases are correctly recovered. The

recovered absorption coefficient images and their profiles along a transect in the sagittal

section are displayed in Figure 3-6 and 3-7, respectively. We see that the reconstructed

absorption coefficients of the bones and cartilage are quantitatively accurate relative to

their actual values. By estimating the full width at half maximum (FWHM) of the

absorption coefficient profiles, the recovered cartilage thickness was found to be









around 1.3 and 2.3mm, which are in good agreement with the actual object size of 1

and 2mm used in the experiments.

We also conducted experiments to evaluate the ability of resolving different optical

contrasts when the absorption coefficient of cartilage varied from 0.015 to 0.04 mm-1.

The reconstructed results from these experiments are shown in Figure 3-8, where

Figure 3-8a to 3-6d display the coronal, sagittal and cross sections of the reconstructed

3-D absorbed energy density when the cartilage had an absorption coefficient of 0.015

mm-1, 0.025 mm-1, 0.03mm-1 and 0.04 mm-1, respectively. We can observe that the

bones and cartilage are clearly distinguishable when the absorption contrast between

the bones and cartilage is large (Figure 3-8a and 3-8b) and that it is increasingly difficult

for these two types of structures to be separated (Figure 3-8c and 3-8d). The recovered

absorption coefficient images and their profiles along a transect in the sagittal section

are displayed in Figure 3-9 and Figure 3-10. We note that the absorption coefficient of

the cartilage is quantitatively recovered-it was found to be around 0.016, 0.026, 0.031

and 0.043 mm-1, which are in very good agreement with the actual values of 0.015,

0.025, 0.03 and 0.04 mm1.

It is worth to note that homogeneous acoustic speed was assumed in this study

which is certainly an approximation to the possible heterogeneous acoustic medium

used in the experiments. Previous studies have shown that this assumption may

generate blurring in the reconstructed images. However, Compared to the blurring

effects from an inappropriate scanning geometry and detection configuration, the

blurring due to the assumption of homogeneous acoustic speed is significantly smaller,

especially for larger dimension targets. In addition, the blurring due to the assumption of









homogeneous acoustic speed can be eliminated by the use of advanced reconstruction

algorithms without such an assumption.

In summary, in this study we have shown that a spherical scanning geometry is

probably an optimal way for us to perform 3-D PAT imaging of the finger joints. Our

results also indicate that the joint-like structures and their absorption coefficients can be

quantitatively reconstructed using our 3-D PAT approach.











Table 3-1. Phantoms used in the pre-investigation of 3-D PAT system
No Ca y Optical absorption Optical scattering Thickness
egory coefficient (mm-1) coefficient (mm-) (mm)

1 Bone 0.07 2.5 N/A

2 Cartilage 0.015 1.5 1.0

3 Cartilage 0.015 1.5 2.0

4 Cartilage 0.025 1.5 1.0

5 Cartilage 0.03 1.5 1.0

6 Cartilage 0.04 1.5 1.0












Light source


Detectors


Light source







4-


Figure 3-1. Schematic of two potential cylindrical scanning geometries of finger joints
imaging in 3-D PAT.









1000
-25
-15-

-15 -10 -5 0 5 10 15
Y (mm)


Figure 3-2. A slice through the axis of the reconstructed finger joint scanned by a
cylindrical geometry shown in Figure 3-1(b)


CD C70













Ti:Sapphire


Figure 3-3. Schematic of the three-dimensional photoacoustic imaging system.



4---..


I Ie


X

I--0


\L ''S ~


Figure 3-4. Schematic of the scanning arc/path in the XY plane.














58












Sagittal Section


IF






Fl
I iGl.


,, a 0 5 1o


Y4on




42
-6j


Y *'
-rt -B 0
Y(Hmt)


Coronal Section


Figure 3-5. Reconstructed absorbed energy density images in the coronal (z=0mm),
sagittal (x=0mm) and cross (y=4mm) sections for a joint-like phantom with
1mm thick cartilage in cylindrical (a) and spherical (b) scanning geometries,
and with 2mm thick cartilage in spherical scanning geometry (c).


Cross Section






20



410 4 0 6 10



4



-4

4 I
-10 4 D 6 10


S10N0





8Om)
44


0-.10 JS 0F 10:


a ig
















0.06
0.06
0.06 5
0.06
0.00.05
0.04 0 0.04

0.03 0.03

-5 0.02 -5 0.02

0.01 0.01
-10 -10
-10 -5 0 5 10 -10 -5 0 5 10
(a) (b)


Figure 3-6. Reconstructed absorption coefficient images (coronal section) for the
phantom with 1mm (a) and 2mm (b) thick cartilage.


4 B 8 10


-2mm thckness
-*-lmm thickness


0.07


E
E



: 9
04

05
0.o3

o 0.02

001-

0i --------------------
-8 -8 4 .-2 0 2
Location along X-0 (mni)


Figure 3-7. The recovered absorption coefficient distributions along a transect (x=Omm)
for the images shown in Figure 3-6.



















60










Sagittal Section


4 Im)


I:




















1o
10!
ii







I-






tn


Figure 3-8. Reconstructed absorbed energy density images in the coronal (z=0mm),
sagittal (x=0mm) and cross (y=4mm) sections for a joint-like phantom where
the absorption coefficient of cartilage was 0.015 (a), 0.025 (b), 0.03 (c) and
0.04 (d) mm1.










61


x Mmn4


Y kih


,i, 1[
I e


I


ii





I


Ay4


(c)


*w 4 1 B 10

-(I

(d)
I
.100


'10 64
Yhimr


Coronal Section


Cross Section


Y mm,


-1P ,@ ..I


Xlmn1














0.07 10
0.06
0.05


003
0.02 -

0 5 10 10
0 5 10 -10 .5


0.07


0.05


.03
.02
.01
0 5 10


10 1U
0.07 0.06

5 0.06 5 0.05
0.05
So 0.04
0.03
-5 0.02 -5 0.02

0.01 0.01
-10 -10
-10 -5 0 5 10 -10 0 10
(c) (d)


Figure 3-9. Reconstructed absorption coefficient images (coronal section) for the
phantom with an absorption coefficient of cartilage of 0.015 (a), 0.025 (b),
0.03 (c) and 0.04 mm-1 (d).


Figure 3-10. The recovered absorption coefficient distributions along a transect
(x=0mm) for the images shown in Figure 3-9.










CHAPTER 4
IMAGING FINGER JOINTS IN-VIVO WITH 3-D PAT IN A SPHERICAL SCAN

In this section, we develop a 3-D PAT system in a spherical scanning geometry

indicated in finger-like phantom study, and attempt to image a distal interphalangeal

(DIP) joint in vivo, which is most frequently affected (57% in women, 37% in men)

among all the joints in the osteoarthritic hand indicated in an incidence study of OA.45 In

vivo joint imaging represents a challenge because unlike the joint from a phantom or a

cadaver finger, an in vivo joint has abundant blood vessels located along the skin which

strongly absorb light. The strong absorption from the blood vessels and bones give rise

to significantly reduced signal-to-noise ratio in an in vivo setting. The investigation in this

section indicate that major features/components in the DIP finger joint as well as the

absorption coefficient profiles of these joint structures can be successfully achieved in-

vivo, by using Ti: Sapphire laser based photoacoustic system in a spherical scanning

geometry.

4.1 System Development and Description

The joint-like phantom experiments conducted in Chapter 3 indicates that a 3-D

PAT approach in an equivalently spherical scanning geometry based on a 1 MHz

transducer with an aperture size of 6mm seems potential to be applied for imaging

human finger joints. Here, we develop a 3-D PAT system in a real spherical scanning

geometry so that the ultrasound detector can be jointly controlled by two rotary stepper

motors to move on a spherical detecting surface surrounding examined in-vivo finger

joints, instead of rotating the examined targets in tissue phantom experiments for

equivalently spherical scanning. The developed 3-D PAT system for in-vivo finger joints

imaging in a spherical scanning geometry consists of a fiber guided near-infrared









lighting subsystem, spherical scanning stepper motors, a ultrasound detector and an

arrangement for hand and finger placement, as shown in Figure 4-1.

4.1.1 Fiber Guided Near Infrared Lighting

Instead of using a reflection mirror to guild the pulsed light beam illuminating the

examined subject from the top of the water tank, the pulsed light beam at 820nm

generated from the Ti: Sapphire laser is guided by an optical fiber to the bottom of the

water tank in in-vivo examination setting. With this illumination approach, there is

enough room to optimize the light beam without interfering with the spherical scanning

mechanical arms. The optical fiber used in our 3-D PAT allows more flexibility in the

lighting approach of the examined finger joints. Although a fixed area source of NIR light

area source is used in our study to visualize the joint cavity, a rotary area source can

also be delivered by our lighting system, which light approach may differentiate the skin

and the inside structures in cross-section of finger joints when a high frequency

transducer is used.25 The fiber guided light beam is expanded by a lens and a ground

glass before reaching the finger so that the DIP joint can be uniformly illuminated by an

area source of light. The laser energy is measured by a power meter on the surface a

finger joint is placed so that the energy we use for the in-vivo examination is well below

the safety standard of 22mJ/cm2 and will not do any damage to human subjects. The

laser energy is also measured after it comes out of the Ti: Sapphire laser to ensure the

energy safety and good working condition of the laser device itself. Before being used

in-vivo, the lighting system is carefully adjusted so that the light beam is centered and

uniform for further finger joints examination. The adjusted light system is verified with

phantoms as shown in Figure 4-2, where the cylindrical phantom is illuminate from the

bottom surface and the line phantom as well as the pencil lead is illuminated from side.









From the testing results shown in Figure 4-2, we see that the light beam is almost

perfectly uniform and centered, and is ready for in-vivo examination.

4.1.2 Arrangement for Hand and Finger Placement

A mechanical arrangement for hand and finger placement was designed and

shown in Figure 4-3, where the hand and the proximal end of the examined finger rest

on the plastic holders cushioned by blue and pink sponges in Figure 4-3(a). The distal

tip of the examined finger slips into a thin plastic ring as shown in Figure 4-3(b). Multiple

ring sizes are available for examined finger joints from different subjects so that the

plastic ring well fits the individual finger joint. Although the plastic ring is close to the

examined DIP joint, it has a little impact on the light illumination and on the propagation

of the ultrasound wave from the light heated DIP joint. The arrangement for hand and

finger placement is screwed to the bottom surface of the water tank, and can be

adjusted for the placement of little finger and index finger of both hands. If the finger

joint is to be illuminated from the dorsal side, the hand should face up resting on the

plastic holder; otherwise, the hand should face down on the plastic holder, which might

be a more comfortable position for human subjects.

4.1.3 Signal Detecting in a Spherical Scanning Geometry

The acoustic field generated by the pulsed laser beam is detected by the same 1

MHz transducer used in phantom experimental studies in Chapter 3, which has a

bandwidth of 80% at -6dB and an aperture size of 6mm. To reduce the attenuation of

the acoustic field, the transducer and the phantom were all immersed in a water tank.

As can be seen from Figure 4-1, the transducer was attached to the rotary stage

R2 whose position was controlled by the rotary stage R1. The collection of acoustic

signals along a spherical scanning surface was realized through the combined rotations









of rotary stages R1 and R2 (Note that R1 and R2 were attached to the same L-shape

arm, see Figure 4-4). As shown in Figure 4-4, the examined finger joint was aligned in

the direction of Y axis, and located at the rotation center, 0. In the scanning, the

transducer was first rotated by R1 to a position along the circular locus L1 and then

rotated by R2 along the circular locus L2 (Figure 4-4a). This process was repeated until

the spherical scanning was completed. Since the finger/hand blocked the scanning of

the transducer in certain angles, the scanning along L1 and L2 covered only 2400 and

2520, respectively (Figure 4-4b and 4-4c).

To run the transducer along the above designed spherical locus, not only the joint

motion of the two rotary stages should be well controlled, but also the position of the

transducer and the rotary centers of the two stepper motors should be well calibrated.

Unlike the 2-D circular scanning shown in Figure 4-5(a) where the error of absolute

starting point only results to slight rotation or slight translation of the reconstructed

targets, mechanical errors in the spherical scanning as shown in Figure 4-5(b) could

lead to distorted scanning locus and big errors in the reconstruction because of model

mismatch. To calibrate the mechanical parts for the spherical scanning geometry, a

bubble level gauge purchased from Sears and a marked paper cap for the transducer

as shown in Figure 4-6 are used. By using the bubble level gauge, we can make sure

that the rotary axes of the stage R1 and R2 are perpendicular to each other. With the

center of the transducer marked on the paper cap, we will know that the circular trace of

the transducer controlled by the rotary stage R2 is on the same plane with the rotary

axis of the stage R1, since the center of the transducer should be a stationary point at

point S shown in Figure 4-4(c) when the rotary stage R1 rotates. After the calibration,









the transducer is able to run the designed spherical locus when it is jointly controlled by

the two rotary stages.

4.1.4 Units Control and Data Acquisition System

The received acoustic signal on each detecting point along the spherical scanning

surface is pre-amplified, amplified and filtered by the same Pulser/Receiver 5058PR

used in phantom experimental studies in Chapter 3, and further feeds to PCI A/D card

(CS121000, Gage Applied Technologies, Lockport, IL, USA) where the received analog

signal is converted into digital signal and stored for further signal processing and image

reconstruction. The CS121000 is in a 12-bit resolution and its sampling rate is 100MS/s.

A Labview program was used to control the entire data acquisition procedure.

4.2 In-vivo Experiments

A female healthy volunteer entered our in-vivo study, and was photoacoustically

examined by our developed 3-D PAT system. The volunteer sat on a chair behind the

water tank, the height of which is adjustable for comfortable positioning of the examined

subject. The low end of the arm and the examined hand were all immersed in a water

tank, which is filled with warm water so that the human subjects will not feel discomfort

during the whole examination procedure. During the examination, the finger joint was

placed at the center of the two rotary stages, and the palmar side of the examined DIP

finger joint faced up allowing the finger joint to be illuminated from the dorsal side of the

finger. Our experience indicates that this way of light illumination can give us maximized

tissue penetration, resulting in clearly differentiable joint cavity in the reconstructed

images. The proximal end and the distal end of the examined finger are secured on the

plastic holder to reduce the motion error during the in-vivo examination. A complete









scanning allowed the collection of signals at 387 positions along the spherical surface

which took about 40 minutes.

4.3 Methods

The collected photoacoustic signals in temporal domain were converted into

signals in frequency domain by Fourier Transforming for our finite element based 3-D

PAT reconstruction algorithm, which is based on finite element solution to the

photoacoustic wave equation in frequency domain subject to the radiation or absorbing

boundary conditions (BCs). The 3-D images of absorbed optical energy density were

reconstructed through our 3-D PAT reconstruction algorithm. The second step is to

recover the optical absorption coefficient distribution by using a forward fitting procedure

based on the photon diffusion model and the reconstructed absorbed optical energy

density images in the first step.29

The following photon diffusion equation as along with the Robin boundary

conditions can be used in the forward fitting procedure to recover the optical absorption

coefficient using the absorbed optical energy density ( reconstructed,

V -D(F)V (Y(F)) u (F)Y(F) = -S(F) (4.1)
-DVY(F). = aY(F) (4.2)

Where Y(F) is the photon density, Y(F) = 0(F)/l ,(f), #u,() is the optical absorption

coefficient, D(F) is the diffusion coefficient, D= 1/(3(#u +,j')) and j' is the reduced

scattering coefficients, a is a boundary condition coefficient related to the internal

reflection at the boundary, and S(F) is the incident point or distributed source term.

To recover the optical absorption coefficient from the reconstructed absorbed

optical energy density (D from an initial distribution of u(F), the optical fluence Y(F)









and the absorbed energy density Dc) are iteratively calculated through the photon

diffusion equation and (c) = /u(F)Y, respectively. If the error between D and V(c) is

not small, then ,(rF) is updated by u,(F) = D(F)/ (F) and the above procedure is

repeated until a small error between D and V(c) is reached, resulting in a stable

quantitative distribution of p (r).

4.4 Results and Discussion

In this section, we present the in vivo 3-D reconstruction results and make

observations based on these results.

A coronal section of the reconstructed 3-D image (absorbed energy density) is

shown in Figure 4-7a, while a similar section of the DIP joint from a typical human finger

by MRI is given in Figure 4-7b for comparison. As can be observed from Figure 4-7a,

the cartilage/joint space is well differentiated from the adjacent distal phalanx (DP) and

intermediate phalanx (IP). The recovered absorption coefficient image and the

absorption coefficient profile along a transect X=-2mm corresponding to the coronal

slice shown in Figure 4-7a are displayed in Figure 4-8a and 4-8b, respectively. We see

that the absorption coefficient values of the phalanxes and cartilage/joint space are

quantitatively consistent with that reported in the literature. By estimating the full width

at half maximum (FWHM) of the absorption coefficient profiles, the recovered thickness

of cartilage/joint space was found to be 1.6mm, which is in agreement with the actual

size.

Figure 4-9a presents a cross section of the reconstructed 3-D image (absorbed

energy density). Again, a similar cross section of the DIP joint from a typical human

finger by MRI is shown in Figure 4-9b for comparison. Comparing both Figure 4-9a and









4-9b, it appears that several joint tissue types are visible including phalanx (PX), lateral

artery (LA) and median artery (MA), and tendon (TE). The recovered absorption

coefficient image and the absorption coefficient profiles along two transects are

displayed in Figure 4-10a to 4-10c. The first transect (line 1 in Figure 4-10a) goes

through the phalanx and tendon, while the second one (line 2 in Figure 4-10a) crosses

the tendon and two arteries. The recovered absorption coefficients of the phalanx,

tendon, median and lateral arteries were found to be 0.07mm-1, 0.074mm-1, 0.064mm-

and 0.058mm-1, respectively. Again, these values are generally consistent with that

reported in the literature. We note that blurs around the two blood vessels are

noticeable (Figure 4-10a), largely due to the relatively low-frequency (1MHz) transducer

used in this study. We believe the resolution can be enhanced with a wideband

transducer having higher frequency (e.g., 5-10MHz) which will make the millimeter-

submillimeter sized arteries accurately imaged.

To view the volumetric structures of the tissues around the joint, a series of cross

and coronal section images are displayed in Figure 4-11 and 4-12, respectively. As

shown in Figure 4-11 a to 1-11f, the phalanx is clearly visible in all the cross section

images from the proximal end (Y=-2.5mm) to the distal end (Y=9mm). The tendon is

differentiated from the phalanx (Figure 4-11 a, 4-11e and 4-11f). The lateral artery is

visible almost in all the cross section images except for that shown in Figure 4-11f at the

distal end of the finger. The median artery is seen in Figure 4-11e and Figure 4-11f, and

weakly visible in Figure 4-11d. Other joint structures including the phalanges (DP and IP)

and the cartilage (CL) are clearly imaged from the coronal section images close to the

dorsal side of the finger (Figure 4-12a and 4-12b). Close to the palmar side of the finger,









the tendon and two arteries begin to show up (Figure 4-12c), and the two arteries

eventually become clearly visible (Figure 4-12d). Negative values are also observed in

the reconstructed absorbed energy density images. It is likely due to the homogenous

acoustic approximation (constant acoustic speed) to the actual heterogeneous acoustic

media of the finger joints and limited signal-to-noise ratio because of the strong light

scattering of joint tissues.

In summary, we have presented a 3-D PAT technique that is able to image finger

joints in vivo. Although it seems impossible for PAT to provide image quality (in terms of

spatial resolution) comparable to MRI for joint or breast imaging, PAT is capable of

obtaining absorption coefficient or functional information. In addition, PAT is portable

and low in cost. While our current experimental setup is not optimized, it allows us to

demonstrate the possibility of 3-D in vivo joint imaging for the first time. The results

obtained indicate that major joint structures and their absorption coefficients can be

quantitatively reconstructed using our 3-D PAT approach.











........ Ti:Sapphir


Figure 4-1. Schematic of the 3-D spherical scanning PAT system for finger joint
imaging. L: lens; T: detector/transducer; R1/R2: rotary stages.


Figure 4-2. Test of lighting beam with a cylindrical phantom (a), a line phantom (b), and
a pencil lead (c).





















(a) (b)


Figure 4-3. Photograph of the arrangement for hand and finger Placement.


z
S \ I
(72
I '


Y

S" S -.L2
25201
T
""'
__-__


\ L1
%


Figure 4-4. Schematic of the spherical scanning geometry (a), circular locus L1 (b) and
L2 (c).

















.-. -,.


.~ li..
... .. .. -._

r--t ---,~





(a) (b)

Figure 4-5. Reconstruction mesh and possible model mismatch of 2-D circular scanning
(a) and 3-D spherical scanning (b).


Figure 4-6. Bubble level gauge (a) and homemade marker cap (b) for calibration of
mechanical errors in the photoacoustic spherical scanner.











DP
10


DIP


-10


-16 -10 -6


0
X (mm)
(a)


5 10 15


Figure 4-7. A selected coronal section from the reconstructed 3-D image (absorbed
energy density) at Z=-5mm (a), and MRI (coronal section) from a similar joint
(b). DP: distal phalanx; IP: intermediate phalanx; DIP: distal interphalangeal
joint.


-15 -10 -5 0 5 10 15
X(mm)
(a)


4
Y(mm)
(b)


Figure 4-8. A selected coronal section from the recovered absorption coefficient image
at Z=-5mm (a), and its profile along the cut line X=-2mm (b).












in-


d 6 -
rG I r


(am)
(a)


Figure 4-9. A selected cross section from the reconstructed 3-D image (absorbed
energy density) at Y=7mm (a), and MRI (cross section) from a similar joint (b).
PX: Phalanx; MA: median artery; TE: tendon; LA: lateral artery.


0.07
E
. 0.065

o 0.06
U
o
| 0.055

S0.05


-15 -10 -5 0 5 10 15
Z(mm)

(a)
0.065

I \
0.06/


0.055 / /


0.05


0.045
2 4 6 8 1(
Distance along the transect(mm)

(c)


4 5 6 7
Distance along the transect(mm)
(b)


Figure 4-10. A selected cross section from the recovered absorption coefficient image at
Y=7mm (a), and its profile along cut lines 1 (b) and 2 (c), respectively.


-- /
I' I

I
I


8










-15 800 -15 800

-10 P 600 -10 P 600



1 200 0 200

10
5 0 5 -200

10 -200 10 M -400

15 15
-10 0 10 -10 0 10
Z (mm) Z mm)
(a) (d)
-15 -15


J-1 -1
S500



50 5
10 I10 20D
15 1-500 15 I4oo
-10 0 10 -10 0 10
Z (mm) Z (mm)

-15 -15
1000 600



-1 i P S 1 .0



10 10
1000 -200

0 10 -400
15 ".,oo 15
-10 0 10 -10 0 10
Z (mml Z (mm)
(c) (f)



Figure 4-11. Cross section images (absorbed energy density) at Y=-2.5mm (a), Y=0mm
(b), Y=3mm (c), Y=5mm (d), Y=7mm (e) and Y=9mm (f). PX: Phalanx; MA:
median artery; TE: tendon; LA: lateral artery.























IP /


15
-16 -10 -6 0 8
Xlmm)
(a)







0 t






-15 -10 4- 0 5
X fmrm)
(b)


10 I6
10o If


10

6




4

-10


-1B
-1


I 1o
1,I

00

200





1-0
10 1 -15


-10i -5 0 6
X (m)I
(c)


E -1 0
-300

10 1i


S1 200
100

S1300



.10 -6 0 6 10 15
X (mmO
(d)


Figure 4-12. Coronal section images (absorbed energy density) at Z=-5mm (a), Z=-
3.5mm (b), Z=lmm (c) and Z=4.5mm (d). DP: distal phalanx; IP: intermediate
phalanx; CL: cartilage; PX: Phalanx; MA: median artery; TE: tendon; LA:
lateral artery.


40I
3 0



!- :i,









CHAPTER 5
CLINICAL EXAMINATION OF OA PATIENTS WITH 3-D PAT

Osteoarthritis (OA) is the most common degenerative joint disease, which affects

tens of millions of Americans. Although knee and hip OA bear most responsibility for the

burden of OA, hand OA may be a marker of a systemic predisposition toward OA;

Incidence study of OA further revealed that the second DIP joint is most frequently

affected (57% in women, 37% in men) among all the joints in the hand.45

Thus far, clinical examination remains the principal approach of OA diagnosis,

relying on symptoms and signs OA patients suffered as well as experience of arthritis

physicians. Imaging techniques, including standard x-ray radiography, computed

tomography (CT), ultrasonography and magnetic resonance imaging (MRI), are used

only as supplemental methods in OA diagnosis, especially in hand OA, either because

they are insensitive to the abnormal changes in the OA joint cavity or because they are

too costly to be used as a routine examination method. As such, sensitive and

affordable imaging methods are urgently needed for the detection of OA, especially in

early stage. Moreover, the progress of effective imaging methods in OA diagnosis may

at the same time accelerate the advancement of medical therapies for OA, which are

currently effective to relief OA symptoms or prevent the worsening of OA in certain

degree and yet of limited effectiveness in modifying OA. Up to now, no approved

disease-modifying OA drugs (DMOADs) are available, development of which may

benefit greatly from the progress of imaging methods in OA diagnosis, where sensitive

imaging techniques can serve as surrogate markers for clinically meaningful outcomes

to economically and efficiently validate the candidate DMOADs.43









New imaging techniques based on near-infrared (NIR) light, including pure-optical

imaging techniques and photoacoustic tomography (PAT), has been recently studied to

image finger joints and to effectively detect joint diseases.4-7'25 55-57 While optical

imaging techniques are able to detect the highly sensitive optical absorption / scattering

abnormalities associated with soft tissues (cartilage, synovial fluid, etc) in OA and RA

joints, its spatial resolution is relatively low (about 3~5mm). Compared to all-optical

imaging techniques, PAT is able to visualize the same optical absorption contrast with

significantly improved spatial resolution (0.5mm or better, adjustable with ultrasound

frequency) for deep-tissue imaging.

In this section, we take the advantage of the high spatial resolution provided by

PAT technique in imaging the optical absorption property, continue previous study in

finger joints imaging using 3-D PAT in a spherical scanning geometry, and present a

pilot clinical study of detecting OA using 3-D PAT on both osteoarthritis patients and

healthy controls.

5.1 Enhancement in 3-D PAT Approach

In the in-vivo finger joints imaging conducted in Chapter 5, a single 1MHz

transducer is used to receive the generated ultrasound signals, and then is feed into a

single channel signal processing unit (preamplifier, amplifier, filter, and A/D conversion).

The whole procedure of examination lasts about 40 minutes with single channel

detecting system, which is too long for OA patients. To speed up the examination

process, an ultrasound detection array together with multiple channel signal processing

unit is designed, calibrated and tested with tissue-like phantom experiments. With

multiple channel detection system, the examination of a finger joint can be completed in

5 minutes, allowing the collection of photoacoustic signals at 360 positions surrounding









the examined DIP joint. DIP finger joints from OA patients and healthy volunteers are

further imaged with our enhanced 3-D PAT system in a spherical scanning configuration,

which are comprised of a pulsed NIR light source as same as the one in in-vivo finger

joint imaging conducted in Chapter 4, a spherical scanning subsystem, and an

ultrasound detection array and associated data acquisition system as shown in Figure

5-1.

5.1.1 Ultrasound Detection Array

The ultrasound detection array (photograph shown in Figure 5-1 and the design

shown in Figure 5-2) is composed of eight 1 MHz transducers (Valpey Fisher,

Hopkinton, MA) arranged uniformly along a 2100 circular arc. The transducer has a

bandwidth of 80% at -6dB and an aperture size of 6mm. The position and performance

of each transducer in the ultrasound array has to be calibrated carefully before

experimental study. Ideally, the rotation of ultrasound detecting array in certain degrees

in the plane where all the transducer elements are arranged should result in the overlap

of the positions and received acoustic signals between different transducer elements.

For example, positions and received signals from transducer T1, T2 until T7 after being

counterclockwise rotated in 300 (together with rotary stage R2) should coincide with

those from transducer T2, T3 until T8 shown in Figure 5-2 before the ultrasound

detecting array was rotated.

In calibrating the system errors in mechanical parts, the center of the ultrasound

detector array must coincide with the rotation axis of the rotary stage R2, and the plane

where all the transducer elements are securely arranged should be perpendicular with

the rotation axis of the rotary stage R2, so that the rotary motion of R2 will not cause

evident position error between transducer elements shown in Figure 5-3. After the









calibration between the center of the ultrasound detecting array and the rotary center of

the stage R2, the position of each transducer element has to be precisely aligned and

secured in the ultrasound detecting array so that each transducer element is in the

same distance from the center of the ultrasound detecting array. The position calibration

of each transducer element is conducted by measuring the pulse-echo response of the

each transducer element from a fixed reflection surface, as shown in Figure 5-4(a)

where pulse-echo responses of fours selected transducer elements are measured in

50MHz sampling rate. From the responses shown in Figure 5-4 (a), location variations

between the four transducer elements are observed and require fine alignment so that

the four pulse-echo responses are of the same flight of time (in the same dashed

vertical line as shown in Figure 5-4(a)). The alignment errors of the transducer elements

can be limited within 0.2mm (about 7 sampling points' variation) by observing the above

pulse-echo responses. The sensitivities and impulse responses of the transducer

elements may also exhibit some differences, which result in the amplitude fluctuation

and slight spectrum shift shown in Figure 5-4(b). We measured the received acoustic

signals of each element on the same location with a well-controlled phantom experiment

and calculated the performance calibration coefficients of all the transducer elements,

which are used in the signal calibration of each experiment conducted with this system.

The same method as employed in Chapter 4 is used to calibrate the altitude of the

plane where all the transducer elements are securely arranged. With the marked cap

shown in Figure 4-6(b), the transducer elements can be finely elevated on the same

plane with the rotary axis of the stage R1 since the center of the transducer element at

point F shown in Figure 5-2 should be a stationary point when the stage R1 rotates.









5.1.2 Multiple Channel Data Acquisition

The detected acoustic signals were feed into preamplifiers and converted to digital

signals by an A/D board (PREAMP2-D and PCIAD1650, US Ultratek, Concord, CA), as

shown in the block diagram of the multiple-channel 3-D PAT system (Figure 5-5). Two

preamplifier sets (PREAMP2-D) are used in our 3-D PAT system, each including four

pre-amplifying channels. The preamplifier PREAMP2-D has a bandwidth from 20 KHz to

30MHz at -3dB, with a gain range from 10 to 60dB. The signal to Noise Ratio of the

preamplifier is 67dB. The PCI card PCIAD1650 consists of 16 pulse/receiver channels,

each of which is integrated with an amplifier (0 dB to 80 dB in 0.01dB increments) and a

signal filter (low pass from 5.9MHz to 48MHz, and high pass from 16 KHz to 4.8MHz).

The sampling rate of the PCIAD1650 is tunable from 390.625 KHz up to 50MHz (8

scales). The central values of the pulse/receiver channels have some variation as

shown in Figure 5-6(a), and require relative calibration between the detecting channels.

The multiple channels A/D board and the entire examination procedure are controlled

by a Labview interface shown in Figure 5-6(b), which include the initiation of step motor

controllers, procession of the spherical scanner(ultrasound detecting array),

communication between the interface with the PCIAD1650, and signal display/storage.

5.1.3 Quantitative 3-D PAT with 3-D Photon Diffusion Model

A forward fitting procedure based on the 3-D photon diffusion model and the

reconstructed 3-D absorbed optical energy density images in 3-D PAT is used to

recover the three-dimensional images of optical absorption coefficient distribution, which

is important for further quantitative analysis of the full joint. This forward fitting method

requires no segmentation and is relatively easy to implement. 29To recover the optical

absorption coefficient from the reconstructed absorbed optical energy density ((rF) from









an initial distribution of/ u(F), the same updating strategy is used in the 3-D forward

fitting procedure as the one used in 2-D forward fitting procedure. 29 Briefly, the optical

fluenceY(F), which relates with /u(F) and D(F) by the identity Y(F)= =(F)/ (F), and

the absorbed energy density D(c)(F) are iteratively calculated through the 3-D photon

diffusion equation and D(c)(F) ,u0(F)Y(F), respectively. If the error between ((rF) and

D(c)(F) is not small, then /u,(F)is updated by u,(F)= '(F)/(F) and the above

procedure is repeated until a small error between (D(F) and OD(c(F) is reached, resulting

in a stable quantitative distribution of pu(F).

5.1.4 Tissue-like Phantom Test

The 3-D quantitative PAT (qPAT) system is validated on tissue mimicking

phantoms before the clinical study was performed, to ensure quantitatively accurate in

vivo results. The joint-like tissue phantoms used in our experiments were optically well-

characterized. The two "bones" had: absorption coefficient = 0.07mm-1, reduced

scattering coefficient u; =2.5mm-1 and a diameter of 6mm, and the "cartilage" between

the bones had: u, = 0.01 or 0.03mm-1, p; =1.5mm-1 and a thickness of 1.5mm. The

optical absorption / scattering coefficients of the joint-like phantoms well approximate

the corresponding coefficients of in-vivo finger joints. Figure 5-7(a) and 5-7(b) show the

reconstructed phantom joints in three-dimensional view, and in sagittal and coronal

sections with our 3-D qPAT. The distributions of the recovered absorption coefficient

along a transect crossing the "bones" and "cartilage" for the images shown in Figure 5-

7(a) and 5-7(b) are shown in Figure 5-7(c). The average absorption coefficient of the

joint cavity when the "cartilage" had absorption coefficient of 0.01 or 0.03mm-1 was









found to be 0.016mm-1 and 0.030mm-1, respectively, while the absorption coefficient of

the "bones" was recovered as 0.067mm-1 and 0.070 mm-1, respectively, for the two

phantom cases. The recovered "cartilage" thicknesses were found to be 1.5mm and

1.7mm, which is in agreement with the exact size.

5.2 Patients and Clinical Data

Seven subjects (two OA patients and five healthy subjects) were enrolled in the

current study between 2008 and 2009. All participants (white female; mean age 61

years, range 45-71 years) provided informed consent as part of the protocol approved

by the Institutional Review Board of University of Florida. Participants were recruited

from the Rheumatology Clinic at the Shands Health Center of the University of Florida

(UF), Gainesville, FL USA. One distal interphalangeal (DIP) finger joint, joint II at the left

hand, from each subject was examined clinically and photoacoustically in this study,

because it was mostly vulnerable to OA disease and at the same time easily accessible

with our current 3-D photoacoustic scanning system.

Clinical examination of each patient was performed independently by a single

rheumatologist (E.S. Sobel). Patients with OA were identified by clinical history and

main clinical features, including symptoms (predominantly pain and stiffness), functional

impairment and signs (joint enlargement and redness). The healthy controls had no

known OA or other joint diseases.

In photoacoustic examination, the examined DIP finger joints were placed at the

center of the two rotary stages R1 and R2, and the palmar side of the examined DIP

finger joint faced up, allowing the finger joint to be illuminated from the dorsal side of the

finger. Both the examined finger and the ultrasound array were all immersed in a water

tank filled with warm water for acoustic coupling. The proximal end and the distal end of









the examined finger were secured on the plastic holder. A complete scanning allowed

the collection of signals at 360 positions along the spherical surface which took less

than 5 minutes.

Finite element based reconstruction algorithm of 3-D quantitative photoacoustic

tomography (qPAT) (i.e., 3-D conventional reconstruction algorithm coupled with 3-D

photon diffusion equation) was used to quantitatively recover the optical absorption

coefficient images of the DIP joints. Full-width at 30% maximum method (FW30%M)

was used in pure-optical imaging modalities to quantitatively analyze the recovered

finger joints, which is able to compensate the inadequate spatial resolution of pure-

optical imaging method in finger joints imaging. While PAT is able to provide high

resolution in imaging the joint cavity (1~2mm) and hence the interfaces between bones

and cartilages, it seems that PAT can only image the synovial fluid (~0.5mm) in poor

resolution. Herein, we combine the FW30%M criteria and the traditional full-width at half

maximum (FWHM) method to quantitatively analyze the recovered joints in our study,

which enables us to get the best quantitative results on different tissues (cartilage and

synovial fluid) in the joint cavity.

Data from six of the seven recruited subjects (2 OA and 4 healthy) resulted in

successful photoacoustic image reconstruction, while data from one healthy control was

not useable due to the severe motion error of finger during the photoacoustic

examination.

5.3 Results

The 3-D photoacoustic images of the DIP joints from both OA patients and healthy

controls were reconstructed using a finite element mesh of 7519 nodes and 41323

tetrahedral elements. Based on the recovered 3-D images, the optical absorption









coefficient of different joint tissues (bone, cartilage and synovial fluid) and the structural

size of joint cavity (cartilage and synovial fluid) are extracted. These parameters are

used as indicators/markers for differentiation of OA from normal joints in this study.

Consecutive coronal and sagittal slices from the recovered 3-D images from the

joint cavity as indicated in Figure 5-8 are presented here. The exact locations of these

selected slices may be slightly different from subject to subject. The coronal and sagittal

slices for a normal joint are shown in Figures 5-9(a) and 5-9(b). As can be observed

from Figures 5-9(a) and 5-9(b), the bones are clearly delineated from the adjacent

tissues (cartilage and fluid) in the joint cavity. While there is no sharp boundary between

the bones and cartilage/fluid, the joint tissue/space is clearly identified in both the

coronal and sagittal sections. Further observation finds that the synovial fluid seems to

be differentiable from the surrounding cartilages in some slices (e.g., Z=7.5mm in the

coronal section and X=-3.0mm in the sagittal section). It is worth noting that the

recovered size of the bones appears to be smaller than the actual anatomic size. We

believe it is mostly due to the limited bandwidth of the ultrasound detectors we used in

the photoacoustic examination where the loss of low frequency signals resulted in the

reduced size of the recovered DIP joint.

Figures 5-10(a) and 5-10(b) display selected coronal and sagittal slices from the

recovered 3-D images for an OA joint. Again, the bones are clearly differentiable from

the adjacent joint tissues (cartilage and fluid), and the joint cavity is easily identified in

both the coronal and sagittal sections. It is interesting to note that the spatial distribution

of absorption coefficient appears to be heterogeneous for the OA joint (Figures 5-10),

while it is relatively homogeneous for the healthy joint (Figures 5-9). This observation is









consistent with the findings from the previous optical imaging study. Compared with the

healthy joint (Figures 5-9), the joint space narrowing seems apparent for the OA joint

(Figures 5-10).

The average absorption coefficient and structural size of joint tissues for all the

subjects are calculated and presented in Tables 5-1. We note from Table 1 that the

recovered absorption coefficient of the diseased synovial fluid ranges from 0.022 to

0.023mm-1, while it is below 0.017mm-1 (0.0089 to 0.017mm-1) for the normal fluid.

Difference between OA joints and healthy controls is also apparent, in terms of the

absorption coefficient of the synovial fluid in the joint cavity. The difference is more

striking when the ratio of the absorption coefficient of the synovial fluid to that of the

bone is considered (0.58~0.63 for OA joints verse 0.25~0.46 for normal joints). For the

diseased cartilage, the absorption coefficient (0.026~0.028mm-1) is generally larger than

that of normal joints (0.015~0.024mm-1). The ratio of the absorption coefficient of the

cartilage to that of the bone further confirms the difference between OA and normal

joints (0.72-0.76 for OA joints verse 0.43-0.61 for normal joints).

We also note from Table 1 that the measured thickness of the synovial fluid for OA

joints ranges from0.3 to 0.5mm, with a mean value of 0.4mm, whereas the range of this

size for healthy joints is from 0.4 to 0.9mm, with a mean a mean value of 0.67mm. For

the diseased cartilage, the measured thickness ranges from 0.4 to 0.5 mm, with a mean

value of 0.45mm, whereas the range of this size for healthy joints is from 0.5 to 0.8 mm,

with a mean value of 0.675mm. Compared with the normal joints, the dimensional

narrowing of both synovial fluid and cartilage is obvious for OA joints. A more obvious

difference is observed in the thickness of the joint cavity (synovial fluid and cartilage)









between OA and normal joints (1.1~1.5mm for OA joints verse 1.7~2.5mm for healthy

joints), where the observable joint space narrowing for OA joints is in agreement with

the findings from x-ray radiography.

5.4 Discussion

In this study, for the first time we in vivo imaged distal interphalangeal (DIP) joints

and associated osteoarthritis in sub-mm resolution using our three-dimensional

quantitative photoacoustic tomography (qPAT) approach. Apparent differences, in both

the reconstructed size and optical absorption coefficient of the joint cavity, are observed

between osteoarthritic and normal joints. The photoacoustic imaging of the joints offers

a spatial resolution comparable to x-ray radiography, and has the potential to be further

enhanced when larger amount of transducer elements and higher frequency

transducers are used. The joint space narrowing observed photoacoustically is in

agreement with the findings from x-ray radiography, and the observed optical changes

associated with diseased joints are in consistent with documented biological roots. For

example, it is known that with the onset of OA, the synovial membrane/fluid in articular

cavity becomes increasingly turbid and has an enhanced vascularization.75 The

increased turbidity and vascularization would accompany increased optical absorption

coefficient (as well as optical scattering coefficient) in the diseased synovial

membrane/fluid. In fact, this optical increase can be as large as 100% at certain

wavelengths in the NIR region. 76-77 In addition, there is increasing evidence that OA is

a disease involving a metabolic dysfunction of bone.78-79 It is likely that this metabolic

dysfunction of bone, often associated with high metabolism of subchondral bone, will

cause changes in its tissue optical absorption due to the changes in oxygenation, Hb,

Hb02 and water content.









An ultrasound array composed of eight 1MHz transducers was used in the current

study, which was sufficient to image the abnormalities in the joint cavity and the

adjacent bones where cartilage erosion has been generally believed to be the major

cause of OA disease. By using 1 MHz transducers, the joint cavity can be clearly

identified and be measured. As shown previously in photoacoustic imaging of cadaver

joints,25 PAT was capable of imaging other interesting joint structures as well, including

volar plate, tendon and aponeurosis when a high frequency transducer was used. In a

recent study of whole-body mouse imaging using 3-D PAT,80 both relatively large

organs such as liver, kidney and spleen, and small structures such as spinal cords,

ovarian vessel, abdominal aorta and femoral veins were clearly imaged when multi-

scale wavelet filter was used along with a 64-element 3.1 MHz ultrasound array. With

these advancements, we believe that our 3-D PAT approach can be further improved to

capture other interesting joint tissues such as ligaments, tendon and synovium in

addition to the articular cavity and the bones. This adds the potential for our 3-D PAT

technique to detect OA that is initiated by these other joint tissues.81'82

It is also worth to note that PAT imaging is governed by ultrasound propagation

behavior in tissue while ultrasound signals are excited by optical pulses. For simplicity,

in this study we have assumed an acoustically homogenous medium in image

reconstruction. The actual inhomogeneous nature of joint tissues will definitely affect the

ultrasound propagation behavior in computation. In this case an advanced

reconstruction algorithm considering acoustical heterogeneity is needed for improved

image reconstruction. Such an advanced algorithm can also reconstruct tissue acoustic

properties such as ultrasound speed in addition to tissue optical properties.38 65 While









the advanced reconstruction algorithm requires increased computation cost, the

recovered tissue acoustic properties together with tissue optical properties and

anatomical structure of joints tissues will offer more accurate detection of OA disease.

For example, decrease of ultrasound speed and increase of ultrasound attenuation has

been observed in osteoarthritic cartilage.83' 84

In summary, this study represents the first attempt to in vivo detect osteoarthritis in

the finger joints using three-dimensional quantitative photoacoustic tomography. We

show that apparent differences, in both the reconstructed size of the joint space and the

absorption coefficient of the joint cavity, exist between OA and normal joints. The results

reported here suggest the possibility of 3-D qPAT as a potential clinical tool for early

detection OA in the finger joints. As an emerging imaging technique, the advantages of

3-D qPAT for early detection of OA lie in its capability of quantitatively imaging optical

absorption, acoustical properties, and physiological/functional parameters such as

oxygenation and water content as well as revealing joint abnormality in high resolution.

In addition, PAT is portable and low in cost, and hence may play an important role in

long term monitoring of OA.






















Figure 5-1. Photograph of the transducer array/finger interface of the three-dimensional
photoacoustic tomography system used for the current study.
T4


ST2


T7


T8


I -K.




-- 0


Figure 5-2. Geometry of the detector array. T1~T8 are eight detectors, O is the center of
the detector array.


ii-I
~ o~-o
i;W,
















ST2


/o
__0 O


L


Figure 5-3. An example of possible mismatches between the center of the detector
array (marked as 0) and the rotary center of the rotary stage R2 (marked as
0').


AII
---4-D-------


(a) (b)

Figure 5-4. Alignment variations (a) and spectrum differences in impulse responses of
four selected transducer elements in the ultrasound detecting array.











93


14. .00


/,








Rotator RI Rotator R2



I I
Controller 1 Controller 2


PC Controlling
Interface


16U


Figure 5-5. Block diagram of the multiple-channel 3-D PAT system.


F^i: Saph ire H










































VISA of Nemark


VISA of Oriel
r 11 4


Sample Rate

Waveform Length

Post Delay
.|f^~


0 2 4


Set Low Pass Fllte

Set High Pass Fllt,


Trigger Source
, ,,


Ol ennahC pened


Gain (dB) DC Offset rmV) Channel Status Calibration Facts











FY s










r i


Card Number


6 8 10 12 14 16 18
Channel No.


(a)


Memory Size Serial Number
I,, I


Original Phall

Phal No. Gamnma No.


.i-i Sbep- -I







.Di Faa File



function (open:O)
rreate or
append to file? (new File:F)
I append to Fle


Figure 5-6. The central values of the 16 channels in PCIAD1650 card (a) and Labview

control interface (b) in multiple-channel 3-D PAT system.


4 ,


I,,.


r


':' '"
':' "










Coronal Section

0.08
0.06
E
0.04


-5 0 5
X (mm)


0.0:


-5 0 5
X(mm)


E
0.04

10.02
0


Sagittal Section

0,08
0.06
0.04


-5 0 5
Z (mm)


i0 0608
ral 10.06


0
Z (mm)


0.04
0.02
0


Location along the transect (mm)


Figure 5-7. Reconstructed optical absorption coefficient images in the coronal
(z=0.0mm) and sagittal (x=0mm) sections as well as in 3-D view for a joint-
like phantom where the absorption coefficient of "cartilage" was 0.01 (a) and
0.03mm-1 (b). (c) shows the recovered absorption coefficient distributions
along a transect crossing the "bones" and "cartilage" (x=0.5mm) for the
images in coronal section .


3-D view

























Figure 5-8. Photographic schematic of the coronal (a) and sagittal (b) sections for the
finger joint imaging.


Z=5.0mm Z=5.5mm Z=6.0mm Z=6.5mm Z=7.0mm Z=7.5mm



(a) Coronal section (dorsal palm)
(an)0 x oonnS X on EafI, '[m' X (II:
(a) Coronal section (dorsal 4 palm)


X=-5.0mm X=-4.5mm X=-4.0mm


X=-3.5mm X=-3.0mm X=-2.5mm

RW: 3

IUri


(b) Sagittal section (Medial Lateral)
Figure 5-9. Recovered images at selected coronal (a) / sagittal (b) planes for case 1
(Healthy).











Z=-9.0mm


X=-3.0mm


* 1

*i


Z=-8.5mm Z=-8.0mm Z=-7.5mm Z=-7.0mm Z=-6.5mm






(a) Coronal section (dorsal palm)
X=-2.5mm X=-2.0mm X=-1.5mm X=-1.0mm X=-0.5mm




it:Z Z ihZa) ZlrIt- Z zn


(b) Sagittal section (Medial Lateral)
Figure 5-10. Recovered images at selected coronal (a) / sagittal (b) planes for case 5
(OA).









Table 5-1. Averaged absorption coefficient and size of joint tissues for 2 OA and 4


normal joints
Average
Joint
N Tissues ,
(mm-1)
Bone 0.035
1 Cartilage 0.015
Fluid 0.008
Bone 0.039
2 Cartilage 0.024
Fluid 0.017
Bone 0.034
3 Cartilage 0.020
Fluid 0.015
Bone 0.037
4 Cartilage 0.023
Fluid 0.017
Bone 0.036
5 Cartilage 0.028
Fluid 0.023
Bone 0.037
6 Cartilage 0.026
Fluid 0.021


Mean
Thickness
(mm)
N/A
0.5
0.7
N/A
0.7
0.7
N/A
0.8
0.9
N/A
0.7
0.4
N/A
0.5
0.5
N/A
0.4
0.3


Joint
Thickness
(mm)

1.7



2.1



2.5



1.8



1.5



1.1


U,(b)


Ula(b)
'Ua(b)


0.43 0.25



0.61 0.45



0.59 0.45



0.61 0.46



0.76 0.63



0.72 0.58









CHAPTER 6
IN-VIVO IMAGING OF CHROMOPHORE CONCENTRATIONS IN FINGER JOINTS
WITH MULTISPECTRAL 3-D PAT

Since last decade, near-infrared spectroscopy and NIR based multispectral

imaging techniques (such as diffuse optical tomography) has emerged as promising

methods to noninvasively image or measure tissue chromophores and associated

hemodynamics in living tissues. In those techniques, concentration distributions of

physiologically / functionally significant tissue chromophores (such as oxy-hemoglobin ,

deoxy-hemoglobin, water, lipid) and their variations due to disease development,

metabolic changes, cerebral hemorrhage, and brain activity can be measured or imaged

by taking the advantage of the absorption spectra of different tissue chromophores.85-94

Thus far, studies has shown that the measurement of these tissue chromophores and

their physiological / functional changes are of great significance both in basic science

(such as neuroscience, molecular imaging of small animals) and clinical applications

(such as breast cancer detection and classification). For example, the concentrations of

oxy-hemoglobin and deoxy-hemoglobin in breast tumors are found to be higher than

normal tissues.96

In recently years, multispectral photoacoustic tomography (PAT) has been

investigated aiming to quantify these same tissue chromophores and the oxygenation

for molecular imaging and functional imaging with high spatial resolution, which can be

as high as subl00/um for several mm penetration depths in highly resolved PAT

images.3542 Although the challenges remain since the conventional PAT images are the

distributions of absorbed optical energy density and not actual images of the absorption

coefficients that can be immediately ready for quantitative spectroscopic analysis,35


100









several models have been proposed and developed to quantitatively resolve the tissue

chromophores with multispectral PAT.36-40 For example, maps of optical absorption

coefficients can be recovered when the traditional PAT is couple with a photon diffusion

model under certain assumptions, such as pre-known scattering property or wavelength

dependency of the scattering property in the medium. The recovered maps of optical

absorption coefficients in multi-wavelengths further go through spectroscopic analysis

for quantitatively resolved concentrations of the tissue chromophores. Some groups

made further attempts to in-vivo measure the concentrations of oxy-hemoglobin and

deoxy-hemoglobin within blood vessels in a human finger,41 and to image the

hemoglobin and oxygenation in a rat brain,42 although more accurate results require

more accurate models of multispectral PAT.

In this section, our goal is to advance the application of multispectral PAT to in-

vivo imaging of chromophore concentrations in finger joints. Although previous studies

focused on 2-D modeled multispectral PAT, a multispectral 3-D PAT approach is used

in our study, which is more accurate and is essential to capture the full volumetric

structures in high resolved PAT images with less artifacts and blurring than an approach

with 2-D approximation, particularly for tissues with high structural heterogeneity (such

as finger joints). Optical wavelengths from 730nm to 880nm are used in our

multispectral 3-D PAT approach since studies have indicated that multispectral light with

wavelengths located on the two sides of 805nm is required to distinguish the

concentrations of oxy-hemoglobin and deoxy-hemoglobin from the combined optical

absorption values.97


101









6.1 Methods

The algorithm used in our multispectral 3-D PAT is a fitting method, which includes

three steps. The first step is to obtain the 3-D images of absorbed optical energy density

D(A,F) in each optical wavelength through our 3-D finite element based PAT

reconstruction algorithm. The second step is to recover the 3-D optical absorption

coefficient distribution for each specific optical wavelength by using a forward fitting

procedure based on the reconstructed 3-D absorbed optical energy density images

D(A,F)in the first step and the following 3-D photon diffusion equation at each specific

optical wavelength A,

V D(A, F)V((A, F)) ,ua (, F))(A, F) = -S(A, F) (6.1)

-D(A, F)VP(A, F). f- = aP(A, F) (6.2)

Where Y(A,F)is the photon density, /u(A,F)is the optical absorption coefficient, D(A, F)

is the diffusion coefficient, D(A, F)= 1/(3(u (A, F)+ ,u(A, F))) and u/(A, F) is the reduced

scattering coefficients, a is a boundary condition coefficient related to the internal

reflection at the boundary, and S(A,F)is the incident point or distributed source term.

To recover the optical absorption coefficient from the reconstructed absorbed

optical energy density, (A,F), the forward fitting procedure starts from an estimated

distribution of u,(A,F)which is an optimized initial based on the results of an searching

scheme for optimal initial parameters. The reduced scattering coefficient /,j (A, ) is

regarded as pre-known constant for simplicity. The optical fluence Y(A,F) and the

absorbed energy density'(c)(A,F) are iteratively calculated through the 3-D photon

diffusion equation andD(c)(A,F)= ,,(A,F)Y(,F), respectively. If the error between


102









(D(2,F) and ODc>(A,F) is not small for a specific optical wavelength, then the

current (A, F) is updated by u, (A, F) = ((A, F) /(A, F) and the above procedure is

repeated until a small error between (D(2,)and O')(A,F) is reached, resulting in a

stable quantitative distribution of /u(A,F).

A third step is to recover the concentrations of major tissue chromophores (oxy-

hemoglobin, deoxy-hemoglobin, and water) from the resolved optical absorption

coefficients maps at each wavelength, based on the absorption spectra of these tissue

chromophores and the following Beer's law in consideration of the optical absorption

contributions from L chromophores.

L
/(A ( F) = ,-(F)c(A) (6.3)
i=i

Where c (F) is the concentration of ithchromophore in the unit of molarM(mole/L);

s,(A) is the absorption extinction coefficient of the ih chromophore at wavelength A. The

chromophores numberL = 2 for the recovery of two tissue chromophores (oxy-

hemoglobin, deoxy-hemoglobin) in our study, which comprise the major contributions to

the optical absorption in the wavelength region we used.

To recover the concentrations of these tissue chromophores, an inverse equation

is used, which is written as the following

[jT + ]{ ()} [ ] T { (F)} (6.4)

Where J, = s (,); the elements {c(r)}and { t,(f)} are cj(F) and i (, ,F) respectively.

p is a regularization factor to stabilize the inverse solution.


103









6.2 In-vivo Experiments

The in-vivo experiments are conducted by our 3-D PAT system described in detail

in Chapter 5, which are comprised of a pulsed NIR light source with wavelength tunable

from 600nm to 950nm, a spherical scanning subsystem, an ultrasound detection array

and associated data acquisition system. Six optical wavelengths from 730nm to 880nm

(730, 760, 805, 825, 850 and 880nm) are chosen for the multispectral 3-D PAT scanner,

which may give the improved distinguishability between HbO2 and Hb. Both the

ultrasound array and the examined finger are immersed in a water tank, which is filled

with warm water, for minimized ultrasound attenuation. During the examination, the

palmar side of the examined DIP finger joint faced up at the center of the two rotary

stages R1 and R2, allowing the finger joint to be illuminated from the dorsal side of the

finger. The collection of the acoustic signals at 360 locations along the spherical

scanning surface surrounding the examined DIP joint requires about 5 minutes in each

wavelength, and the entire procedure of multispectral 3-D photoacoustic examination

lasts for about 30 minutes. Again, the proximal end and the distal end of the examined

finger are secured on the plastic holder to reduce the motion error during the entire in-

vivo examination procedure. To reduce the impact of possible motion errors on our

multispectral PAT, our multispectral 3-D PAT scanner operates in the following

wavelength order: 730nm, 805nm, 850nm, 880nm, 825nm, 760nm. A homemade 3-D

positioning ruler (Figure 6-2) is used to check the location of the examined finger joint

before and after the whole examination procedure, to eliminate failed examinations with

obvious finger motion during the examination procedure.

Seven subjects (two OA patients and five healthy subjects) entered the current

study between 2008 and 2009. All participants (white female; mean age 61 years, range


104









45-71 years) provided informed consent as part of the protocol approved by the

Institutional Review Board of University of Florida. One distal interphalangeal (DIP)

finger joint, joint II at the left hand, from each subject was examined clinically and

photoacoustically in this study, because it was mostly vulnerable to OA disease and at

the same time easily accessible with our current multispectral 3-D PAT scanner. Clinical

examination of each patient was performed independently by a single rheumatologist

(E.S. Sobel). The healthy controls had no known OA or other joint diseases.

Four of the seven recruited subjects (two OA patients and two healthy controls)

completed the whole procedure of 3-D multispectral photoacoustic scanning resulted in

successful photoacoustic image reconstruction and quantitative measurement of oxy-

hemoglobin and deoxy-hemoglobin in our study.

6.3 Results and Discussion

The 3-D photoacoustic images of the DIP joints from all four human subjects were

reconstructed by using a spherical mesh of 7519 finite element nodes and 41323

tetrahedral elements. The optical absorption maps of reconstructed DIP joints are

further recovered under each optical wavelength by coupling reconstructed 3-D

photoacoustic images (deposited optical absorption energy) of the examined DIP joints

with 3-D photon diffusion model. Among the six optical wavelengths where the DIP

joints are reconstructed and their optical absorption maps are quantitatively recovered,

three optical wavelengths are selected, based on the relatively structural consistencies

of the recovered optical absorption maps, for further spectroscopic analysis to

quantitatively resolve the concentrations of HbO2 and Hb.

Figure 6-3(a) and 6-3(b) are the recovered optical absorption maps in coronal

section and sagittal section, respectively, of an in-vivo finger joint at wavelengths from


105









730nm to 880nm. As can be observed in Figure 6-3(a) and 6-3(b), the joint tissue/space

is differentiated from the adjacent bones in the coronal and/or sagittal sections at

several optical wavelengths (760nm, 805nm, 850nm and 880nm). Relatively structural

consistencies of the recovered optical absorption maps are observed between the

wavelengths in 760nm, 805nm and 850nm, which are selected to further quantitatively

resolve the tissue chromophore concentrations. The reason that the joint tissue/space is

missing or fused with the adjacent bones at some optical wavelengths (such as 730nm

and 825nm) are most likely due to motion error during the 30 minutes examination

procedure, although an optimized wavelength order is followed in our 3-D multispectral

PAT scanning. We believe that the motion error and the structural inconsistencies of the

recovered optical absorption maps at some wavelengths can be significantly reduced by

using an ultrasound array integrated with much more ultrasound detectors (for example,

64 detector elements as being used in high resolution whole-body mouse imaging).80

The quantitatively resolved concentrations of HbO2 and Hb from a typical OA

finger joint, which is the same finger joint shown in Figure 6-3, are displayed in Figure 6-

4. The concentrations of HbO2 (Figure 6-4(a)) in different joint tissues are clearly

differentiated in high resolution both in coronal section and sagittal section, and it is

relatively higher in the bones than in the joint cavity. Figure 6-4(b) shows the

concentrations of Hb both in coronal section and sagittal section, where the joint cavity

is differentiated from the adjacent phalanges in the sagittal section and is barely

identifiable in the coronal section. We believe that the relatively poor quality of the

resolved Hb concentration images in coronal section is most likely due to the relatively

low quality image at the wavelength 760nm, where the absorption spectrum of Hb has a


106









high peak and motion error may have happened during the 3-D multispectral PAT

scanning.

Figure 6-5 shows the quantitatively resolved concentrations of HbO2 and Hb from

a typical healthy finger joint. Again, the structure of joint space is clearly identified in the

images of each resolved chromophore, and the concentrations of HbO2 and Hb in the

joint cavity are found to be relatively lower than that in the adjacent phalanges. In

comparison with the OA finger joint shown in Figure 6-4, the joint cavity in healthy

subjects seems more distinctive.

For all the four examined in-vivo finger joints, the average concentrations of oxy-

hemoglobin (HbO2), total hemoglobin (Total-Hb) and the oxygen saturation (SO2) and

H20) in the joint cavity and the adjacent bones are calculated and presented in Figure

6-6. As can be observed from Figure 6-6(a), the average concentrations of HbO2 and

Total-Hb in the finger phalanges are in the range from 44 to 58 uM, and from 68 to

96 ,uM, respectively. The quantitatively resolved concentration levels of HbO2 and

Total-Hb are in agreement with the measurements by diffuse near-Infrared

spectroscopy. 98 Further observations indicate a difference in the levels of total-Hb

between finger phalanges with OA disease and healthy phalanges, which are in range

of 86 ~ 96/,M and 68 ~ 69uM respectively. The finding of enhanced hemoglobin levels

in OA joints has the pathologic roots revealed by previous studies, where angiogenesis

has been observed in diseased tissues. The oxygen saturations of the four in-vivo

finger joints are observed from 58% to 66% in the finger phalanges, while finger

phalanges with OA disease seems to have lower oxygen saturations levels (55%~60%)

than healthy phalanges (65%~66%). Hypoxia in osteoarthritic phalanges has been


107









revealed and seemed to appear in pre-arthritic stage by previous studies of finger joints

diseases. 99 Our results are in great agreement with previous findings of hypoxia in

osteoarthritic phalanges.

As shown in Figure 6-6(b), the average concentrations of HbO2 and Total-Hb in

the finger joint cavity are in the range from 16 to 29 u/M, and from 27 to 53 /uM,

respectively. An obvious difference in the levels of total-Hb between joint cavities with

OA disease and healthy joint cavities can be observed in Figure 6-6(b), which are in

range of 46 ~ 53/,M and 27 ~ 42,uM, respectively. The oxygen saturations in the joint

cavities are observed from 55% to 69%, and osteoarthritic fingers seem to have lower

oxygen saturations levels (55%~59%) than healthy phalanges (59%~69%) in the joint

cavities. Again, the enhanced hemoglobin levels and dropped oxygen saturations in OA

joint cavities are in agreement with the pathologic angiogenesis and hypoxia related

with finger joint diseases.

In summary, major chromophore concentrations (HbO2 and Hb) of finger joints

have also been imaged in-vivo using multispectral 3-D PAT approach with six optical

wavelengths from 730nm to 880nm. The total hemoglobin levels of finger joints are

found to range in 68 96/IM and 27 53/IM for joint phalanges and joint cavities,

respectively. The oxygen saturations vary from 58% to 66% in joint phalanges, and from

55% to 69% in joint cavities. The Enhanced hemoglobin levels and dropped oxygen

saturations in osteoarthritic phalanges and joint cavities have been observed, and the

findings are in agreement with the pathologic angiogenesis and hypoxia related with

finger joint diseases. The multispectral results obtained in this study indicated that the

multispectral 3- PAT approach might be capable of detecting the angiogenesis and


108









hypoxia in pre-arthritic stage of OA disease, and is promising as a potential clinical tool

for early detection OA in the finger joints.


109



















E0.31



I0.25

0 0.
a
J 0.15


Wavelength (nm)


Figure 6-1. Absorption spectra of oxy-hemoglobin (HbO2), deoxy-hemoglobin (Hb) and
water (H20) in biological tissues.


Figure 6-2. Homemade 3-D positioning ruler.


110











730nm

0.07
0.06
0.05S
004
0.03
0.02


-5 0 5
X (mm)
825nm



10.03

0 02



I o
-5 0 5
X (mm)


730nm


-10 -5
Z (mm)
825nm


0.06
0.05
0.04
'0.03
10.02
0.01
0


7Rnnm


0.06
0.05
0.04


0.03
0.02
0.01

-5 0 5
X (mm)
850nm

0.05
'0.04
0.03 E


S.01

-5 0 5
X (mm)


760nm


-10 -5
Z(mm)
RFfnnm


I0.04
0.03 _
E
0.02 -

0.01

0


805nm
0.07
.06
.05
0.04
0.03



-5 0 5
X(mm)
880nm





003
I 0.02

0 5 001
-5 0 5
X (mm)


805nm


II::
0.05

0.04

0.03





-10 -5 0 0
Z (mm)
880nm


10 0.03 10

S.0 I00 10.03
E 0.02 0
> 0 0 015 0 0.02

r 10.01
.5 0 .0! o,00 .5

-10 -5 0 -10 -5 0
Z (mm) Z (mm)


-10 -5
Z(mm)


Figure 6-3. Recovered optical absorption maps in coronal section (a) and sagittal
section (b) of an in-vivo finger joint for light source in six wavelengths.


111


I0.03
0.03
0.02:
0.02
001!
001
0.001













10 10 70
80 60

560 5 50
E 6 E
Sm 40 ((a)
> 0 40 "0 o 30
20
20
-5 -5 10

-5 0 5 0 -10 -5 0
X (mm) Z (mm)

10 .50 10 50


5 -40 5 40
E E
E 30 E 30
0 o (b)
20 20

-5 10 -5 10

-5 0 5 0 -10 -5 0 0
X (mm) Z (mm)

Figure 6-4. Concentrations of oxy-hemoglobin (a) and deoxy-hemoglobin (b) in uM-, and
the water content (c) at coronal section (z=-8mm) and sagittal section (x=-
1mm) of examined finger joint from an OA patient.


112


Z=-Rmm


X=-lmm













X=-O.5mm


-5 0
X (mm)


-5 0 5
X (mm)


30

20

10

0

35

30

25

20
E
15

10

5

0


-5 0
Z (mm)


3 Z (
Z (mm)


Figure 6-5. Concentrations of oxy-hemoglobin (a) and deoxy-hemoglobin (b) in pM, and
the water content (c) at coronal section (z=2mm) and sagittal section (x=-
0.5mm) of examined finger joint from a healthy subject


113


Z=2mm


60

50

40

30 (a)

20

10

0











IHbO ~
90 ZTolal-Hb 90


70 70
"1
r 0*
~0 60 m
0
40 40



1 1 10


Subject No.


(a)


100 100
mHbO2
90 Total-Hb 90
M SO2

S70 70 S
So s So
0 l50




01: 1 I I I I -I I -I 1
: 40 40



0


1 2 3 4
Subject No.


(b)

Figure 6-6. Average Concentrations of oxy-hemoglobin(HbO2), total hemoglobin (Total-
Hb) and the oxygen saturation (802) in finger joint bones (a) and in joint
cavity (b).


114









CHAPTER 7
CONCLUSION AND FURTHER WORK

7.1 Conclusion

In this dissertation, we have investigated the feasibilities of 3-D PAT approach to

image in-vivo finger joints and to detect the OA disease in the hand. Also multispectral

3-D PAT is studied to measure the major tissue chromophore concentrations (HbO2, Hb

and H20) of in-vivo finger joints.

Results based on finger-like phantom experiments, in-vivo finger joint imaging and

clinical cases study show that 3-D PAT in a spherical scanning geometry is able to

image in-vivo finger joints in high resolution and is promising to detect the finger OA in

the hand. Major anatomical structures of an examined distal interphalangeal (DIP) joint

along with the side arteries are clearly imaged, where joint space is well differentiated

from surrounding finger phalanges. Seven subjects enrolled in our study are

photoacoustically examined with our 8-channel 3-D PAT system, and apparent

differences, in both the reconstructed size of the joint space and the absorption

coefficient of the joint cavity, are observed between OA and normal joints. Major

chromophore concentrations (HbO2, Hb) of in-vivo finger joints are quantitatively imaged

by using 3-D multispectral PAT approach with six optical wavelengths from 730nm to

880nm, and the quantitative results are consistent with previous studies.

As our studies indicate, the parallel computing technique, 3-D PAT approach in a

spherical scanning geometry, efficient optical lighting method and calibrated multiple

ultrasound detecting system are essential for effective finger joints reconstruction, OA

detection and chromophores recovery. An ultrasound detecting system with much more

detectors is expected to further improve the images quality and quantitative results.


115









Although the image quality (in terms of spatial resolution) provided by our 3-D PAT

seems not comparable to MRI, our 3-D PAT is able to obtain optical absorption

coefficient and physiological / functional information of the joint tissues (cartilage, fluid,

phalanges, etc) in higher spatial resolution over pure optical techniques. In addition,

PAT is portable and low in cost. While our current 3-D PAT approach is not optimized, it

allows us to demonstrate the possibility of 3-D in in-vivo finger joint imaging and clinical

cases study of osteoarthritis for the first time, which suggests the possibility of 3-D PAT

as a potential clinical tool for early detection of OA in the finger joints.

7.2 Future Studies

7.2.1 Evaluation of 3-D PAT in OA Detection with Small Clinical Samples

The ultimate goal of this dissertation and our study is to detect OA joints from

normal joints at an early stage. The advancement of the 3-D PAT both in the hardware

and algorithms introduced in this dissertation serve the goal, and studies on limited

clinical OA cases have already indicated that it is promising to detect the finger

osteoarthritis in the hand by imaging abnormal rise of the optical absorption coefficients

in joint cavity (cartilage and synovial fluid) and the narrowing of the joint space with3-D

PAT. However further evaluation with small clinical samples plays a key role to validate

our methods, where the specificity and sensitivity of the 3-D PAT in the detection of OA

patients will be investigated.

7.2.2 In-vivo Detection of Rheumatoid Arthritis in the Hand

Rheumatoid arthritis is another common joint disease, which is characterized as

inflamed joint. Although the 3-D PAT is initially developed to image DIP joint for

detection of OA disease in the hand, it is adjustable for PIP joint imaging and RA

disease detection in the hand. Compared to the DIP joint imaging, PIP joint imaging


116









may involve more ultrasound attenuation in the bones, resulting in further drop of signal

to noise ratio (SNR). A PIP joint from a healthy human subject is photoacoustically

examined by the adjusted 3-D PAT system, and the reconstructed joint cavity seems

differentiated from the surrounding phalanges consecutively in both coronal section and

sagittal section as shown in Figure 7-1. A photoacoustic study on another healthy

human subject indicates that our 3-D PAT seems capable of identifying the synovial

fluid in the joint, which may be interesting and important for rheumatoid arthritis

detection since inflammation in the synovial fluid / synovial membrane is severe in

rheumatoid joints. More in-vivo studies of PIP finger joints imaging is planned for further

clinical cases study, and human subjects (RA patients and healthy controls) has being

recruited.

7.2.3 Recover Acoustic Property Together with Optical Properties by PAT
Approach

Homogeneous acoustic speed was assumed in our study which is certainly an

approximation to the heterogeneous acoustic medium (finger joints) in reality.

Considering the acoustic heterogeneity in photoacoustic equation as shown in Equation

(2.6), the variation of sound speed in the acoustically heterogeneous medium may be

recovered. Figure 7-2 is the simulation geometry of a testing case for advanced 2-D

PAT with consideration of acoustic heterogeneity in photoacoustic equation, where the

optical contrast between the targets and the background is 2:1 and the ultrasound

contrast (ultrasound speed) is 1550m/s Vs.1485.5m/s. The lower target has both optical

contrast and ultrasound contrast, while the upper left target and upper right target have

only optical contrast and ultrasound contrast, respectively. As shown in Figure 7-3,

strong crosstalk between the optical contrast and ultrasound contrast is observed in first


117









several iterations (Figure 7-3(a) and 7-3(b)), however the ultrasound speed and its

contrast can be effectively recovered with more iteration in inverse reconstruction. We

are planning to use advanced 3-D PAT to recover the ultrasound speed and its variation

in in-vivo finger joints, which may further help quantify the OA disease in finger joints

since the rise in ultrasound speed has been observed in OA finger joints from previous

ultrasound studies of OA disease.


118











X=-2.6mm Pltne

'iI

SmmI1

a mm


t.- mm. n


I20
10
i


1:
is 'i

'0


r..lmmr4n.


Z [mm) Z tmm e Z {mim
(a) Sagittal section (medial -4 lateral)


X(mm) X ln > {.-
(b) Coronal section (dorsal palm)


5 0
X{(mm)


Figure 7-1. Reconstructed images at selected sagittal (a) / coronal (b) planes for an in-
vivo PIP joint.









V Y I r o=14W 87.Wis
Sbackgrot d =l 1.0 Y \ 6rgt=155..l .


D=4nm

O O



(a) (b)

Figure 7-2. Simulation geometry to recover ultrasound speed with PAT method. (a) is
the simulation geometry of optical contrasts; (b) is the simulation geometry of
ultrasound contrasts;


119


X-.6mm PlaRn


Z imn,


- K.-m 01l-
































-10 0 10


15 2 15 ..
15

10 1.8 1


6 1.6 5

1.4


-10 1
1 1
-15
-15 -10 -6 0 5 10 15 -15


1


15

10o

5

0

-5

-10

-15


-10 0 10

















.10 0 10


-15 -10 -5 0 5 10


Figure 7-3. Reconstructed optical contrast images (A) and ultrasound speed images (B)
at 1st iteration (a), 3rd iteration (b) and 5th iteration (c).


120


-10 -5 0 5 10


xlO

1.65


1.6


1.55 (a)


1.5


1.45


x 10

1.58

1.56

1.54

1.52

1.5

1.48


x 10
1.56


1.54


1.52 (C)


1.5


11.48


15

10

5

0



-10


2

1.8

* 1.6

1.4

1.2


'"









LIST OF REFERENCES


1. H. Jiang, N. Iftimia, Y. Xu, J. Eggert, L. Fajardo, and K. Klove, "Near-infrared optical
imaging of the breast with model-based reconstruction," Acad. Radiol. 9, 186-194
(2002).
2. X. Gu, Q. Zhang, M. Bartlett, L. Schutz, L. Fajardo, and H. Jiang, "Differentiation of
cysts from solid tumors in the breast with diffuse optical tomography," Acad. Radiol.
11, 53-60 (2004).
3. A. E. Cerussi, A. Berger, F. Bevilacqua, N. Shah, D. Jakubowski, J. Butler, R.
Holcombe, and B. Tromberg, "Sources of absorption and scattering contrast for
nearinfrared optical mammography," Acad. Radiol. 8, 211-218 (2001).
4. A. H. Hielscher, A. D. Klose, A. K. Scheel, B. Moa-Anderson, M. Backhaus, U. J.
Netz, et al, "Sagittal laser optical tomography for imaging of rheumatoid finger
joints," Phys. Med. Biol. 49, 1147-1163 (2004).
5. A. Pifferi, A. Torricelli, P. Taroni, A. Bassi, E. Chikoidze, E. Giambattistelli, "Optical
biopsy of bone tissue: a step toward the diagnosis of bone pathologies," J. Biomed.
Opt. 9, 474-480 (2004).
6. Z. Yuan, Q. Zhang, H. Jiang, "3D diffuse optical tomography imaging of
osteoarthritis: Initial results in finger joints," J. Biomed. Opt. 12, 034001 (2007).
7. A. K. Scheel, M. Backhaus, A. D. Klose, B. Moa-Anderson, U. J. Netz, K. G.
Hermann, et al, "First clinical evaluation of sagittal laser optical tomography for
detection of synovitis in arthritic finger joints," Ann. Rheum. Dis. 64, 239-245 (2005).
8. A. A. Oraevsky, A. A. Karabutov and S. V. Solomatin, "Laser optoacoustic imaging
of breast cancer in vivo," Proc. SPIE 4526, 6-11 (2001).
9. S. Manohar, A. Kharine, J. C. G. van Hespen, W. Steenbergen and T. G. van
Leeuwen, "Photoacoustic mammography laboratory prototype: imaging of breast
tissue phantoms," Phys. Med. Biol. 50, 2543-2557(2005).
10.R. G. M. Kollkman, J. H. G. M. Klaessens, E. Hondebrink, J. C. W. Hopman, F. F. M.
de Mul, W. Steenbergen, J. M. Thijssen and T. G. V. Leeuwen, "Photoacoustic
determination of blood vessel diameter," Phys. Med. Biol. 49, 4745-4756(2004).
11.R. I. Siphanto, K. K. Thumma, R. G. M. Kolkman, T. G. van Leeuwen, F. F. M. de
mul, J.W. van Neck, L.N.A. van Adrichem, and W. Steenbergen, "Serial noninvasive
photoacoustic imaging of neovascularization in tumor angiogenesis," Opt. Express
13, 8995 (2005).
12.S. Yang, D. Xing, Q. Zhou, L. Xiang and Y. Lao, "Functional imaging of
cerebrovascular activities in small animals using high-resolution photoacoustic
tomography," Med. Phys. 34, 3294-3301 (2007).
13.X. Wang, Y. Pang, G. Ku, X. Xie, G. Stoica, and L. V. Wang, "Noninvasive laser-
induced photoacoustic tomography for structural and functional in vivo imaging of
the brain," Nat. Biotechnol. 21, 803-806 (2003).


121









14. E. Z. Zhang, J. Laufer and P. Beard, "Three dimensional photoacoustic imaging of
vascular anatomy in small animals using an optical detection system," Proc. SPIE
6437, 643710S (2007).
15.A. A. Karabutov, E. Savateeva, A. Oraevsky, "Imaging of layered structures in
biological tissues with opto-acoustic front surface transducer," Proc. SPIE 3601,
284-295 (1999).
16.C. G. A. Hoelen, F. F. M. de Mul, R. Pongers, and A. Dekker, "Three-dimensional
photoacoustic imaging of blood vessels in tissue," Opt. Lett. 23, 648-650 (1998).
17.J. A. Viator, G. Au, G. Paltauf, S. Jacques, S. Prahl, H. Ren, Z. Chen, J. Nelson,
"Clinical testing of a photoacoustic probe for port wine stain depth determination,"
Lasers Surg. Med. 30, 141-148 (2002).
18. S. Y. Emelianov, P. Li and M. O'Donnell, "Photoacoustics for molecular imaging and
therapy," Phys. Today. 05, 34-39 (2009).
19.D. Razansky, M. Distel, C. Vinegoni, R. Ma, N. Perrimon, R. W. Ko"ster and V.
Ntziachristos, "Multispectral opto-acoustic tomography of deep-seated fluorescent
proteins in vivo," Nat. Photonics 3, 379-379 (2009).
20. H. F. Zhang, K. Maslov, G. Stoica and L. V. Wang, "Functional photoacoustic
microscopy for high-resolution and noninvasive in vivo imaging," Nat. Biotechnol. 24,
848-851 (2006).
21. H. F. Zhang, K. Maslov, G. Stoica and L. V. Wang, "In vivo imaging of subcutaneous
structures using functional photoacoustic microscopy," Nat. Protoc. 2, 797-804
(2007).
22.A. D. L. Zerda, C. Zavaleta, S. Keren, et al., "Carbon nanotubes as photoacoustic
molecular imaging agents in living mice," Nat. Nanotechnol. 3, 557-562 (2008).
23. L. S. Boucharda, M. S. Anwarb, G. L. Liu, et al., "Picomolar sensitivity MRI and
photoacoustic imaging of cobalt nanoparticles," PNAS 106, 4085-4089 (2009).
24. Q. Zhang, Z. Liu, P. R. Carney, Z. Yuan, H. Chen, S. N Roper and H. Jiang, "Non-
invasive imaging of epileptic seizures in vivo using photoacoustic tomography,"
Phys. Med. Biol. 53, 1921-1931 (2008).
25.X. Wang, D. L. Chamberland, D. A. Jamadar, "Noninvasive photoacoustic
tomography of human peripheral joints toward diagnosis of inflammatory arthritis,"
Opt. Lett. 32, 3002-3004 (2007).
26. D. L. Chamberland, X. Wang, B. J. Roessler, "Photoacoustic tomography of
carrageenan-induced arthritis in a rat model," J. Biomed. Opt. 13, 011005 (2008).
27. D. L. Chamberland, A. Agarwal, N. Kotov, J. B. Fowlkes, P. L Carson and X. Wang,
"Photoacoustic tomography of joints aided by an Etanercept-conjugated gold
nanoparticle contrast agent-an ex vivo preliminary rat study," Nanotechnology 19,
095101 (2008).
28. J. Ripoll, V. Ntziachristos, "Quantitative point source photoacoustic inversion
formulas for scattering and absorbing medium," Phys. Rev. E 71, 031912 (2005).


122









29.Z. Yuan, H. Jiang, "Quantitative photoacoustic tomography: recovery of optical
absorption coefficient map of heterogeneous medium," Appl. Phys. Lett. 88, 231101
(2006).
30. B. Cox, S. Arridge, K. Kostli and P. Beard, "2D quantitative photoacoustic image
reconstruction of absorption distributions in scattering medium using a simple
iterative method," Appl. Opt. 45, 1866-1875 (2006).
31. L. Yao, Y. Sun, and H. Jiang, "Quantitative photoacoustic tomography based on the
radiative transfer equation," Opt. Lett. 34, 1765-1767 (2009).
32. L. Yin, Q. Wang, Q. Zhang, H. Jiang, Tomographic imaging of absolute optical
absorption coefficient in turbid medium using combing photoacoustic and diffusing
light measurements," Opt. Lett. 32, 2556-2558 (2007).
33. L. Yao, Y. Sun, and H. Jiang, "Transport-based quantitative photoacoustic
tomography: Simulations and experiments," Phys. Med. Biol. 55, 1917-1934 (2010).
34.Z. Yuan, Q. Wang, H. Jiang, "Reconstruction of optical absorption coefficient maps
of heterogeneous media by photoacoustic tomography coupled with diffusion
equation based regularized Newton Method," Opt. Express 15, 18076-18081
(2007).
35. B. T. Cox, J. G. Laufer and P. C. Beard, "The challenges for quantitative
photoacoustic imaging," Proc. SPIE 717713, 717713-1-9 (2009).
36. J. Laufer, D. Delpy, C. Elwell and P. Beard, "Quantitative spatially resolved
measurement of tissue chromophore concentrations using photoacoustic
spectroscopy: application to the measurement of blood oxygenation and
haemoglobin concentration," Phys. Med. Biol. 52, 141-168 (2007).
37. B. T. Cox, S. R. Arridge, and P. C. Beard, "Estimating chromophore distributions
from multiwavelength photoacoustic images," J. Opt. Soc. Am. A 26, 443-455
(2009).
38.Z. Yuan and H. Jiang, "Simultaneous recovery of tissue physiological and acoustic
properties and the criteria for heterogeneous media by quantitative photoacoustic
tomography," Opt. Lett. 34, 1714-1716 (2009).
39. Z. Yuan and H. Jiang, "Quantitative photoacoustic tomography," Phil. Trans. R. Soc.
A 367, 3043-3054 (2009).
40. J. Laufer, B. Cox, E. Zhang, and P. Beard, "Quantitative determination of
chromophore concentrations from 2D photoacoustic images using a nonlinear
model-based inversion scheme," Appl. Opt. 10, 1219-1233 (2010).
41. J. Laufer, E. Zhang, P. Beard, "Quantitative in vivo measurements of blood oxygen
saturation using multiwavelength photoacoustic imaging," Proc. SPIE 6437,
64371Z-1-9 (2007).
42.X. Wang, X. Xie, G. Ku, and L. V. Wang, "Noninvasive imaging of hemoglobin
concentration and oxygenation in the rat brain using high-resolution photoacoustic
tomography," J. Biomed. Opt. 11, 024015-1-9 (2006).


123









43. C. V. Oddis, "New Perspective on Osteoarthritis," Am. J. Med. 100 2A, 10S-15S
(1996).
44. R. W. Moskowitz, R. D. Altman, M. C. Hochberg, J. A. Buchwalter, V. M. Goldberg,
"Osteoarthritis: diagnose and medical/surgical management," Lippincott Williams &
Wilkins, Philadelphia, USA (2007).
45. P. Sarzi-Puttini, M. A. Cimmino, R. Scarpa, R. Caporali, F. Parazzini, A. Zaninelli, et
al, "Osteoarthritis: an overview of the disease and its treatment strategies," Semin.
Arthritis. Rheum. 35, 1-10 (2005).
46. S. B. Abramson, M. Attur, Y. Yazici, "Prospect for disease modification in
osteoarthritis," Nat. Clin. Prat. Rheumatol. 2, 304-312 (2006).
47. P. Peloschek, G. Langs, M. Weber, J. Sailer, M. Reisegger, H. Imhof, H. Bischol, F.
Kainberger, "An automatic model-based system for joint space measurements on
hand radiographs: Initial experience," Radiology 243, 855-862 (2007).
48. J. T. Sharp, "Assessment of radiographic abnormalities in rheumatoid arthritis: what
have we accomplished and where should we go from here," J. Rheumatol. 22,
1787-1791 (1995).
49.T. Nishii, H. Tanaka, N. Sugano, H. Miki, M. Takao, H. Yoshikawa, "Disorders of
acetabular labrum and articular cartilage in hip hysplasia: evaluation using isotropic
high-resolution CT arthrography with sequential radial reformation," Osteoarthritis.
Cart. 15, 251-57 (2007).
50. H. I. Keen, R. J. Wakefield, A. J. Grainger, E. M. A. Hensor, P. Emery, P. G.
Conaghan, "An Ultrasonographic Study of Osteoarthritis of the Hand: Synovitis and
Its Relationship to Structural Pathology and Symptoms," Arthritis. and Rheumatism.
59, 1756-1763 (2008).
51.J. Chao and K. Kalunian, "Ultrasonography in osteoarthritis: recent advances and
prospects for the future," Current Opinion in Rheumatology 20, 560-564 (2008).
52.J. A. Tyler, P. J. Watson, H. L. Herrod, M. Robson, L. D. Hall, "Detection and
monitoring of progressive cartilage degenerationof osteoarthritic cartilage by MRI,"
Acta. Orthoped. Scand. Suppl 266, 130-138 (1995).
53. H. Sugimoto, A. Takeda, K. Hyodoh, "Early-stage rheumatoid arthritis: Prospective
study of the effectiveness of MR imaging for diagnosis," Radiology 216, 569-575
(2000).
54. H. Sugimoto, A. Takeda, J. Masuyama, M. Furuse, "Early-stage rheumatoid arthritis:
diagnostic accuracy of MR imaging," Radiology 198, 185-192 (1996).
55.Z. Yuan, Q. Zhang, E. Sobel, H. Jiang, "Tomographic x-ray guided three-
dimensional diffuse optical tomography of osteoarthritis in the finger joints," J.
Biomed. Opt. 13, 044006 (2008).
56.Y. Sun, and H. Jiang, "Quantitative three-dimensional photoacoustic tomography of
the finger joints: Phantom studies in a spherical scanning geometry," Phys. Med.
Biol. 54, 5457-5467 (2009)


124









57.Y. Sun, E. Sobel and H. Jiang, "Quantitative three-dimensional photoacoustic
tomography of the finger joints: an in-vivo study," J. Biomed. Opt. 14, 064002-1-5
(2009).
58. Y. Sun, E. Sobel and H. Jiang, "First assessment of three-dimensional quantitative
photoacoustic tomography for in vivo detection of osteoarthritis in the finger joints,"
Radiology (under review).
59. Y. Sun, E. Sobel and H. Jiang, "In vivo imaging of chromophore concentrations in
finger joints by three-dimensional multispectral photoacoustic tomography," Opt.
Exp. (under submission).
60. K. P. Kostli, M. Frenz, H. Bebie, and H. P. Weber, "Temporal backward projection of
optoacoustic pressure transients using Fourier transform methods," Phys. Med. Biol.
46, 1863-1872 (2001).
61. C. G. A. Hoelen, F. F. de Mul, R. Pongers, and A. Dekker, "Three dimensional
photoacoustic imaging of blood vessels in tissue," Opt. Lett. 23, 648-650 (1998).
62. P. Y. Liu, "The P-transform and photoacoustic image reconstruction," Phys. Med.
Biol. 43, 667-674 (1998).
63.Y. V. Zhulina, "Optical statistical approach to optoacoustic image reconstruction,"
Appl. Opt. 39, 5971-5977 (2000).
64. K. P. Kostli, D. Frauchiger, J. Niederhauser, G. Paltauf, H. Weber, and M. Frenz,
"Optoacoustic imaging using a three-dimensional reconstruction algorithm," IEEE J.
Sel. Top. Quantum Electron. 7, 918-923 (2001).
65. H. Jiang, Z. Yuan, and X. Gu, "Spatially varying optical and acoustic property
reconstruction using finite-element-based photoacoustic tomography," J. Opt. Soc.
Am. A. 23, 878-888 (2006).
66.Z. Yuan, Q. Zhang, H. Jiang, "Simultaneous reconstruction of acoustic and optical
properties of heterogeneous media by quantitative photoacoustic tomography," Opt.
Exp. 14, 6749-6754 (2006).
67.Z. Yuan and H. Jiang, "Three-dimensional finite-element-based photoacoustic
tomography: Reconstruction algorithm and simulations," Med. Phys. 34, 538-546
(2007).
68. Q. Fang, "Computational Methods for Microwave Medical Imaging," PhD
dissertation, Dartmouth College (2004).
69. K. D. Paulsen, P. M. Meaney, and L. Gilman, "Alternative Breast Imaging: Four
Model-Based Approaches," Springer Science+Business Media (2005).
70. K. D. Paulsen and H. Jiang, "Spatially varying optical property reconstruction using a
finite element diffusion equation approximation," Med. Phys. 22, 691-701 (1995).
71. K. D. Paulsen, "A Dual Mesh Scheme for Finite Element Based Reconstruction
Algorithms," IEEE Trans. Med. Imag. 14, 504-514 (1995).


125









72. H. Jiang, K. D. Paulseny, U. L. Osterbergy and M. S. Pattersonz, "Improved
continuous light diffusion imaging in single- and multi-target tissue-like phantoms,"
Phys. Med. Biol. 43, 675-693 (1998).
73.X. Gu, Y. Xu, and H. Jiang, "Mesh-based enhancement schemes in diffuse optical
tomography," Med. Phys. 30, 861-869 (2003).
74. 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).
75.A. D. Meisel, P.G. Bullough, "Atlas of Osteoarthritis," Gower Medical Publishing,
New York (1984).
76. J. Beuthan, V. Prapavat, R. Naber, O. Minet, G. Muller, "Diagnostic of inflammatory
rheumatic diseases with photon density waves," Proc. SPIE 2676, 43-53 (1996).
77.V. Prapavat, W. Runge, J. Mans, A. Krause, J. Beuthan, G. Muller, "The
development of a finger joint phantom for the optical simulation of early inflammatory
rheumatic changes," Biomedizinische. Technik. 42, 319-326 (1997).
78. J. Mansell, et al, "Bone, not cartilage, should be the major focus in osteoarthritis,"
Nat. Clin. Prat. Rheumatol. 3, 306-307 (2007).
79. L. Knott, A. Bailey, "Collagen cross-links in mineralizing tissues: a review of their
chemistry, function and clinical relevance," Bone 22, 181-187 (1998).
80. H. Brecht, R. Su, M. Fronheiser, S. A. Ermilov, A. Conjusteau, A. A. Oraevsky,
"Whole-body three-dimensional optoacoustic tomography system for small animals,"
J. Biomed. Opt. 064007, 064007-1-8 (2009).
81.A. L. Tan, H. Toumi, M. Benjamin, A. J. Grainger, S. F. Tanner, P. D. Emery, D.
McGonagle, "Combined high-resolution magnetic resonance imaging and
histological examination to explore the role of ligaments and tendons in the
phenotypic expression of early hand osteoarthritis," Ann. Rheum. Dis. 65, 1267-
1272 (2006).
82. K. D. Brandt, E. L. Radin, P. A. Dieppe, L. van de Putte, "Yet more evidence that
osteoarthritis is not a cartilage disease," Ann. Rheum. Dis. 65, 1261-1264 (2006).
83. S. L. Myers, K. Dines, D. A. Brandt, K. D. Brandt, M. E. Albrecht, "Experimental
assessment by high frequency ultrasound of articular cartilage thickness and
osteoarthritic changes," J. Rheumatol. 22, 109-116 (1995).
84. G. A. Joiner, E. R. Bogoch, K. P. Pritzker, M. D. Buschmann, A. Chevrier, F. S.
Foster, "High frequency acoustic parameters of human and bovine articular cartilage
following experimentally-induced matrix degradation," Ultrason. Imaging 23, 106-
116(2001).
85. F. F. Jobsis, "Noninvasive, infrared monitoring of cerebral and myocardial oxygen
sufficiency and circulatory parameters," Science 198, 1264-1267 (1977).
86. M. S. Patterson, B. Chance, and B. C. Wilson, "Time resolved reflectance and
transmittance for the noninvasive measurement of tissue optical properties," Appl.
Opt. 28, 2331-2336 (1989).


126









87.M. Tamura, O. Hazeki, S. Nioka, B. Chance, and D. Smith,"Simultaneous
Measurements of Tissue Oxygen Concentration and Energy State by Near-Infrared
and Nuclear Magnetic Resonance Spectroscopy," Adv. Exp. Med. Biol. 222, 359-63
(1988).
88. M. Cope and D. T. Delpy, "System for Long-Term Measurement of Cerebral Blood
and Tissue Oxygenation on Newborn Infants by Near Infrared Transillumination,"
Med. Biol. Eng. Comput. 26, 289-294 (1988).
89. P. W. McCormick, M. Stewart, M. G. Goetting, M. Dujovny, G. Lewis, and J. I.
Ausman, "Noninvasive cerebral optical spectroscopy for monitoring cerebral oxygen
delivery and hemodynamics," Crit. Care. Med. 19, 89-97 (1991).
90. S. P. Gopinath, C. S. Robertson, R. G. Grossman, and B. Chance, "Nearinfrared
spectroscopic localization of intracranial hematomas," J. Neurosurg. 79, 43-47
(1993).
91.T. Kato, A. Kamei, S. Takashima, and T. Ozaki, "Human visual cortical function
during photic stimulation monitoring by means of near-infrared spectroscopy," J.
Cereb. Blood. Flow. Metab. 13, 516-520 (1993).
92. B. Chance, Z. Zhuang, C. UnAh, C. Alter, and L. Lipton, "Cognitionactivated low-
frequency modulation of light absorption in human brain," Proc. Natl. Acad. Sci.
U.S.A. 90, 3770-3774 (1993).
93.A. Villringer, J. Planck, C. Hock, L. Schleinkofer, and U. Dirnagl, "Near infrared
spectroscopy (NIRS): a new tool to study hemodynamic changes during activation of
brain function in human adults," Neurosci. Lett. 154, 101-104 (1993).
94.Y. Hoshi and M. Tamura, "Dynamic multichannel near-infrared opticalimaging of
human brain activity," J. Appl. Physiol. 75, 1842-1846 (1993).
95.A. Yodh, B. Chance, "Spectroscopy and imaging with diffusing light," Phys. Today
48, 34-40 (1995).
96. S. Srinivasan, B. W. Pogue, S. D. Jiang, H. Dehghani, C. Kogel, S. Soho, J. J.
Gibson, T. D. Tosteson, S. P. Poplack and K. D. Paulsen, "Interpreting hemoglobin
and water concentration, oxygen saturation, and scattering measured in vivo by
near-infrared breast tomography," Proc. Natl. Acad. Sci. USA 100, 12349-12354
(2003).
97.Y. Yamashita, A. Maki, and H. Koizumi, "Wavelength dependence of the precision of
noninvasive optical measurement of oxy-, deoxy-, and total-hemoglobin
concentration," Med. Phys. 28, 1108-1114 (2001).
98. B. W. Pogue, S. P. Poplack, T. O. McBride, W. A. Wells, K. S. Osterman, U. L.
Osterberg, K. D. Paulsen, "Quantitative Hemoglobin Tomography with Diffuse Near-
Infrared Spectroscopy: Pilot Results in the Breast," Radiology 218, 261-266 (2001).
99. C. H. Jeon, J. K. Ahn, J. Y. Chai, H. J. Kim, E. K. Bae, S. H. Park, E. Y. Cho, H. S.
Cha, K. S. Ahn, E. M. Koh, "Hypoxia appears at pre-arthritic stage and shows co-
localization with early synovial inflammation in collagen induced arthritis," Clin. Exp.
Rheumatol. 26, 646-648 (2008).


127









BIOGRAPHICAL SKETCH

Yao Sun received his B.S. in electrical engineering from Northeast Institute of

Electric Power in July 2000 and his M.S. in physics, with specialty in acoustics, from

Tsinghua University in 2003. After two years' study in the PhD program of physics at

Washington State University, he decided to devote himself to the society of biomedical

engineering and entered The University of Florida in July 2005 for his PhD education in

biomedical engineering. His research focused on high performance photoacoustic

tomography and its application to the imaging of finger joints and detection of

osteoarthritis in the hand.


128





PAGE 1

1 THREE-DIMENSIONAL PHOTOACOUSTIC TOMOGRAPHY AND ITS APPLICATION TO DETECTION OF JOINT DISEASES IN THE HAND By YAO SUN A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORID A IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 2010

PAGE 2

2 2010 Yao Sun

PAGE 3

3 To my wife, Xiaorong Li; My mom and my dad; And all who have been supportive to me in my life

PAGE 4

4 ACKNOWLEDGMENTS First, I would like to thank Professor Huabei J iang, who provided this great education opportunity to me. I appreciate the full financial support and precious advice I have received from Professor Jiang during my study for the PhD degree. His broad views and deep understandings on the research subjects guided me all through this research. His awareness to timetable and duty priorities, and his persistence to research objectives deeply impressed me. What I have learned and experienced in Professor Jiangs lab wil l definitely benefit my future career and life. Secondly, my thanks go to my academic committee members: Dr. David Gilland, Dr. Rosalind Sadleir, and Dr. Sihong Song. T heir comments and suggestions were very helpful to my study and research. Their va luable time spent on reading my dissertation is highly appreciated. Thirdly, I would like to thank Dr. Eric S. Sobel, an associate professor in the College of Medicine at University of Florida. As our clinical research partner, Dr. Sobel helps us to recruit sufficient patients with os teoarthritis. And I woul d like to thank the research assistant professors in our l ab, Dr. Qizhi Zhang and Dr. Zhen Yuan, for valuable discussions on experiments and computations, respectively. Last but not least, I would like to take this opportunity to expre ss my thanks to Lei Xi, Qiang Wang, Lei Yao, Yiyong Tan, Ruixin Jiang, Alexandria Gr ubbs, Xiaoping Liang, Lu Yin, Lijun Ji, Bo Wang, and all members in our biomedical optics lab for their assistance and the great ti me we spent together.

PAGE 5

5 TABLE OF CONTENTS page ACKNOWLEDG MENTS .................................................................................................. 4 LIST OF TABLES............................................................................................................ 7 LIST OF FIGURES .......................................................................................................... 8 LIST OF ABBR EVIATIONS ........................................................................................... 12 ABSTRACT................................................................................................................... 13 CHA PTER 1 BACKGROUND AND S IGNIFIC ANCE ................................................................... 16 1.1 Introduction to P hotoacoustic Tomogr aphy ....................................................... 16 1.2 Overview of Os teoarthritis................................................................................. 18 2 RECONSTRCUTION AL GORITHM IN PAT ........................................................... 24 2.1 Finite Element Based Rec onstruction Algo rithm in PAT ................................... 24 2.1.1 Forward Modeling in Ac oustically Ho mogenous M edia ......................... 26 2.1.2 Inverse Modeling and Nonlinear Optimization in PAT Recons truction.. 27 2.1.3 Regularization in PAT Recons truction ................................................... 30 2.1.4 Simulation in PAT Recons truction ......................................................... 33 2.2 Computing Strategies in PAT Rec onstruc tion ................................................... 33 2.2.1 Adjoin Sensitivity Method ...................................................................... 34 2.2.2 Dual Meshing Sc heme .......................................................................... 34 2.2.3 Partial Reconstructi on in Region of Interest .......................................... 37 2.2.4 Parallel Computing Tec hnique .............................................................. 39 3 PRE-INVESTIGATION OF 3-D PAT SY ST EM WITH TISSUE PHANTOM STUDY 47 3.1 Tissue Mimi cking Phantom ............................................................................... 48 3.2 System Descripti on and Expe riments ............................................................... 49 3.3 Me thods ............................................................................................................ 50 3.4 Results and Discu ssion ..................................................................................... 52 4 IMAGING FINGER JOINTS INVIVO WITH 3-D PAT IN A SPHERICAL SCAN..... 63 4.1 System Developm ent and Descr iption .............................................................. 63 4.1.1 Fiber Guided Near Infrared Lighting ...................................................... 64 4.1.2 Arrangement for Hand and Finger Placement ....................................... 65 4.1.3 Signal Detecting in a Spherical Scann ing Geometry ............................. 65 4.1.4 Units Control and Da ta Acquisiti on System ........................................... 67

PAGE 6

6 4.2 In-vivo Experiments .......................................................................................... 67 4.3 Me thods ............................................................................................................ 68 4.4 Results and Discu ssion ..................................................................................... 69 5 CLINCAL EXAMINATION OF OA PATIENTS WITH 3-D PAT ................................ 79 5.1 Enhancement in 3-D PAT Ap proach ................................................................. 80 5.1.1 Ultrasound De tection Array ................................................................... 81 5.1.2 Multiple Channe l Data Acqu isition ......................................................... 83 5.1.3 Quantitative 3-D PAT with 3-D P hoton Diffusion Model ........................ 83 5.1.4 Tissue-like Phantom Test ...................................................................... 84 5.2 Patients and Clinical Data ................................................................................. 85 5.3 Results .............................................................................................................. 86 5.4 Discussion ........................................................................................................ 89 6 IN-VIVO IMAGIING of CHROMOPHOR E CONCENTRATIONS in FINGER JOINTS WITH MULT ISPECTRAL 3-D PAT .......................................................... 100 6.1 Me thods .......................................................................................................... 102 6.2 In-vivo Experim ents ........................................................................................ 104 6.3 Results and Discu ssion ................................................................................... 105 7 CONCLUSION AND FUTHER WORK .................................................................. 115 7.1 Conc lusion ...................................................................................................... 115 7.2 Future Studi es ................................................................................................ 116 7.2.1 Evaluation of 3-D PAT in OA De tectio n with Small Clinical Samples.. 116 7.2.2 In-vivo Detection of Rheum atoid Arthritis in the Hand ......................... 116 7.2.3 Recover Acoustic Property Toget her with Optical Properties by PAT Approach ........................................................................................... 117 LIST OF RE FERENCES ............................................................................................. 121 BIOGRAPHICAL SKETCH .......................................................................................... 128

PAGE 7

7 LIST OF TABLES Table page 3-1 Phantoms used in the pre-in vestigation of 3-D PAT system ............................... 56 5-1 Averaged absorption coefficient and size of joint tissues for 2 OA and 4 normal jo ints ....................................................................................................... 99

PAGE 8

8 LIST OF FIGURES Figure page 1-1 Schematic of photoacoustic effect. ..................................................................... 23 2-1 The reconstructed images for simulated finger joint (a) coronal sec tion slice along X=0mm (b) cross sect ion slice al ong Z=15mm......................................... 42 2-2 Geometry of t he absorbed energy interpolation at fine node i in 2-D case (a) and in 3-D case (b), and the calculation of the Jacobian matrix element in the effective integr ation zo ne (c).............................................................................. 42 2-3 Slice along the axis of the cylindric al target from rec onstructed 3-D images with full volume reconstruction (a), and parti al reconstruction in region of interest (b).......................................................................................................... 43 2-4 Reconstructed 2-D images of a circul ar target in a phantom experiment with full volume reconstruction (a), and partial reconstruction in region of interest (b). ...................................................................................................................... 43 2-5 Schematic of parallel computing with a combination of both shared memory and distribut ed memo ry...................................................................................... 44 2-6 Flow chart of the high perform ance photoacoustic tomography based on parallel computing te chnique and finite element meth od.................................... 45 2-7 The reconstructed images for simula ted blood vessels (a) and crossed hairs in phantom experiment (b).................................................................................. 46 3-1 Schematic of two potential cylindric al scanning geometries of finger joints imaging in 3-D PAT. ........................................................................................... 57 3-2 A slice through the axis of the reconstructed finger joint scanned by a cylindrical geometry show n in Figure 3-1(b)....................................................... 57 3-3 Schematic of the three-dim ensional photoacoustic imaging system. .................. 58 3-4 Schematic of the scanni ng arc/path in the XY plane. ......................................... 58 3-5 Reconstructed absorbed energy densit y images in the coronal (z=0mm), sagittal (x =0mm) and cross (y=4mm) sect ions for a joint-like phantom with 1mm thick cartilage in cylindrical (a) and spherical (b) scanning geometries, and with 2mm thick cartilage in s pherical scanning geometry (c)....................... 59 3-6 Reconstructed absorption coeffici ent images (coronal section) for the phantom with 1mm (a) and 2 mm (b) thick cartil age. ........................................... 60

PAGE 9

9 3-7 The recovered absorption coefficient distributions along a transect (x =0mm) for the images show n in Figur e 3-6.................................................................... 60 3-8 Reconstructed absorbed energy densit y images in the coronal (z=0mm), sagittal (x =0mm) and cross (y=4mm) secti ons for a joint-like phantom where the absorption coefficient of cartilage was 0.015 (a), 0.025 (b), 0.03 (c) and 0.04 (d) mm-1...................................................................................................... 61 3-9 Reconstructed absorption coeffici ent images (coronal section) for the phantom with an absorption coefficient of cartilage of 0.015 (a), 0.025 (b), 0.03 (c) and 0.04 mm-1 (d).................................................................................. 62 3-10 The recovered absorption coefficient distributions along a transect (x =0mm) for the images show n in Figur e 3-9.................................................................... 62 4-1 Schematic of the 3-D spherical s canning PAT system for finger joint imaging... 72 4-2 Test of lighting beam with a cylindrical phan tom (a ), a line phantom (b), and a pencil l ead (c)..................................................................................................... 72 4-3 Photograph of the arrangement for hand and finger Placement. ........................ 73 4-4 Schematic of the spherical scanning geometry (a), circular locus L1 (b) and L2 (c). ................................................................................................................. 73 4-5 Reconstruction mesh and possible model mismatch of 2-D circular scanning (a) and 3-D spherical scanning (b)..................................................................... 74 4-6 Bubble level gauge (a) and homemade marker cap (b) for calibration of mechanical errors in the phot oacoustic spher ical sc anner................................. 74 4-7 A selected coronal section from the reconstructed 3-D image (absorbed energy density) at Z=-5mm (a), and MRI (c oronal section) from a si milar joint (b)....................................................................................................................... 75 4-8 A selected coronal section from the re covered absorption coefficient image at Z=-5mm (a ), and its profile al ong the cut li ne X=-2mm (b)................................. 75 4-9 A selected cross secti on from the reconstructed 3-D image (absorbed energy density) at Y=7mm (a), and MRI (cross se ction) from a sim ilar joint (b)............. 76 4-10 A selected cross section from the recovered absorpt ion coefficient image at Y=7mm (a), and its profile along cut li nes 1 (b) and 2 (c), respectively. ............. 76 4-11 Cross section images (absorbed energy density) at Y=-2 .5mm (a), Y=0mm (b), Y=3m m (c), Y=5mm (d), Y=7mm (e) and Y=9mm (f)................................... 77

PAGE 10

10 4-12 Coronal section images (absorbed energy density) at Z=-5mm (a), Z=-3.5mm (b), Z=1 mm (c) and Z=4.5 mm (d)....................................................................... 78 5-1 Photograph of the transducer array/finger interf ace of the three-dimensional photoacoustic tomography system used for the cu rrent study. ........................... 92 5-2 Geometry of t he dete ctor array........................................................................... 92 5-3 An example of possible mismatches bet ween the center of the detector array (marked as O) and the rotary center of the rotary stage R2 (marked as O)....... 93 5-4 Alignment variations (a) and spectrum differences in impulse resp onses of four selected transducer elements in the ultrasound de tecting array.................. 93 5-5 Block diagram of the mult iple-channel 3D PAT system. .................................... 94 5-6 The central values of the 16 channels in PCI AD1650 card (a) and Labview control interface (b) in mult iple-channel 3D PAT system................................... 95 5-7 Reconstructed optical absorption coefficient images in the coronal (z=0.0mm) and sagittal (x =0mm) sections as well as in 3-D view for a jointlike phantom where the absor ption coefficient of cartilage was 0.01 (a) and 0.03mm-1 (b). (c) shows the recovered abs orption coefficient distributions along a transect crossing the bones and cartilage (x=0.5mm) for the images in coronal section .................................................................................. 96 5-8 Photographic schematic of the coronal (a) and sagittal (b) sections for the finger joint imagi ng. ............................................................................................ 97 5-9 Recovered images at selected coronal (a) / sagittal (b) planes for case 1 (Health y). ............................................................................................................ 97 5-10 Recovered images at selected coronal (a) / sagittal (b) planes for case 5 (OA). ................................................................................................................... 98 6-1 Absorption spectra of oxy-hemoglobin (HbO2), deoxy-hemoglobin (Hb) and water (H2O) in biologic al tiss ues....................................................................... 110 6-2 Homemade 3-D positioning ruler. ..................................................................... 110 6-3 Recovered optical absorpt ion maps in coronal section (a) and sagittal section (b) of an in-vivo finger joint for li ght source in si x wavelengths......................... 111 6-4 Concentrations of oxy-hemoglob in (a) and deoxy-hemo globin (b) in M and the water content (c) at coronal se ction (z=-8mm) and sagittal section (x=1mm) of examined finger jo int from an OA patient........................................... 112

PAGE 11

11 6-5 Concentrations of oxy-hemoglob in (a) and deoxy-hemo globin (b) in M and the water content (c) at coronal se ction (z=2mm) and sagittal section (x=0.5mm) of examined finger joint from a heal thy subject................................... 113 6-6 Average Concentrations of oxy-hemoglobin(HbO2), total hemoglobin (TotalHb) and the oxygen saturation (SO2) in finger joint bones (a) and in joint cavity (b)........................................................................................................... 114 7-1 Reconstructed images at selected sagittal (a) / coronal (b) planes for an invivo PIP joint. .................................................................................................... 119 7-2 Simulation geometry to recover ultr asound speed with PAT method. (a) is the simulation geometry of optical contrast s; (b) is the simulation geom etry of ultrasound c ontrasts;........................................................................................ 119 7-3 Reconstructed optical contrast images (A) and ultrasound speed images (B) at 1st iteration (a), 3rd iteration (b) and 5th iteration (c)....................................... 120

PAGE 12

12 LIST OF ABBREVIATIONS 2-D Two-dimensional 3-D Three-dimensional CT Computed tomography DOT Diffuse Optical Tomography DIP Distal Interphalangeal FE Finite Element MRI Magnetic Resonance Imaging NIR Near-infrared OA Osteoarthritis PAT Photoacoustic Tomography PIP Proximal Interphalangeal qPAT Quantitative Phot oacoustic Tomography RA Rheumatoid Arthritis US Ultrasonography

PAGE 13

13 Abstract of Dissertation Pr esented to the Graduate School of the University of Florida in Partial Fulf illment of the Requirements for t he Degree of Doctor of Philosophy THREE-DIMENSIONAL PHOTOACOUSTIC TOMOGRAPHY AND ITS APPLICATION TO DETECTION OF JOINT DISEASES IN THE HAND By Yao Sun August 2010 Chair: Huabei Jiang Major: Biomedical Engineering This thesis research presents the study of three-dimensiona l (3-D) photoacoustic tomography (PAT) and its application for the firs t time to detection of osteoarthritis (OA) in the hand. PAT is an emerging non-ionizi ng, non-invasive imaging modality that can visualize high optical contrast in biologic al tissues with high ultrasound resolution. Compared with other imaging moda lities that have been conventionally used or recently investigated to visualize the structural abno rmalities in the finger joints with OA, the 3-D PAT approach studied in this dissertation not only provides high-resolution anatomical structures, but also offers quantitative tissue optical property as well as physiological / functional information including conc entrations of oxy-hemoglobin (HbO2), deoxyhemoglobin (Hb) and water (H2O) that could be used to detect OA in an early stage. In this study, a 3-D high performance PAT reconstruction algorithm is developed based on parallel computing technique and finite element method. The optimal detector performance and scanning geometry for 3-D photoac oustic imaging of the finger joints are investigated with consi derable phantom experiments. The results of phantom experiments show that a spherical scanni ng geometry appears to provide improved spatial solution over a cylindrical scanning geometry, and that 1mm thick cartilage can

PAGE 14

14 be accurately differentiated from the bones with a 1 MHz transducer in a spherical scanning geometry. In addition, the absorption coefficient of the cartilage can be effectively recovered when this opt ical property varied from 0.015mm-1 to 0.04mm-1. A 3-D PAT system in a spherical scanning geometry has been cons tructed and optimized for in-vivo examination of hum an finger joints. A distal interphalangeal (DIP) joint from a female healthy subject was photoacoustically examined by our 3-D PAT system, and major anatomical structures of the examined DIP joint along with the side arteries were clearly reconstructed in high quality, wher e joint space was well differentiated from surrounding finger phalanges. The performanc e of the 3-D PAT system was further improved with an ultrasound detection array co mposed of eight 1 MHz transducers and a 16-channel pulse/receive board. The 8-channel 3-D PAT system was carefully calibrated with controlled tiss ue phantom experiment s, and is capable of completing a finger joint examination within 5 minutes (at single optical wavelength). The 3-D PAT reconstruction algorithm and scanning system developed has also been applied to a pilot clinical study aiming to test the possibility of detecting OA in the hand joints using 3-D PAT. In this pilot clinic al study, seven subjects (two OA patients and five healthy controls) were enrolled a nd photoacoustically examined. The image quality of the reconstructed finger joints was greatly im proved with the 8-channel 3-D PAT system, and apparent differences, in both the reconstructed size of the joint space and the absorption coefficient of the joint cavity, has been observed between the OA and normal joints. The successful results obtained suggest the possibility of 3-D PAT as a potential clinical tool for ear ly detection of OA in the fi nger joints. Major chromophore concentrations (HbO2 and Hb) of in-vivo finger joints have also been quantitatively

PAGE 15

15 imaged using multispectral 3-D PAT approach with six optical wavelengths from 730nm to 880nm. The multispectral results obtaine d further confirmed that the 3D PAT approach implemented in this thesis research is able to differentiate OA from normal joints. While we target the detection of OA as a testing base for validating the singleand multi-spectral 3-D PAT approaches developed in this thesis research, many aspects of our work are fundamental to imaging in gene ral. For example, t he 3-D PAT approaches implemented are applicable to ot her biomedical problems such as rheumatoid arthritis (RA) detection and functional brain imaging.

PAGE 16

16 CHAPTER 1 BACKGROUND AND SIGNIFICANCE 1.1 Introduction to P hotoacoustic Tomography Photoacoustic tomography (PAT), also referred as optoacoustic tomography (OAT), is an emerging non-ionizing and non-in vasive imaging modality in the field of biomedical imaging based on the physical principal of photoac oustic effect, where light energy is converted into acoustic energy due to optical absorption and localized thermal expansion in biological tissues as shown in the schematic of photoacoustic effect in Figure 1-1. In photoacoustic tomography in biomedicine, pulsed laser beam with duration in nanoseconds is delivered into bi ological tissues, and some of the delivered energy will be absorbed and converted into hea t, leading to transient thermoelastic expansion and thus wideband ultrasonic em ission. The generated ultrasonic waves (~MHz) are then detected by ultrasonic transducers to form images. As a hybridized imaging technique between optical imaging modality and ultrasound imaging modality, PAT combines both high optical contrast and high ultrasound resolution in a single modality. St udies have shown that optical absorption contrast between tumor and normal tissues in the breast can be as high as 3:1 in nearinfrared region due to the significantly increased vascularity in the tumors.13 Significant absorption contrast between diseased and norma l joints has also been observed in the hand with osteoarthritis (OA) and rheumatoid arthritis (RA).47 For example, for an osteoarthritic joint, the ratio of its cartil age absorption coefficient to that of the associated bone is increased by 40% relative to the healthy joints. Unfortunately, the high optical contrast in biological tissues was poorly imaged in the pure optical imaging methodologies, since optical scattering in soft tissues degrades spatial resolution

PAGE 17

17 significantly with depth. For example, for imaging depth greater t han ~1mm, the spatial resolution provided by optical diffusion to mography is no better than 3mm. With PAT technique, the high optical absorption contrast in biological tissues is recovered by measuring the megahertz acoustic pressure wave with ultrasound detectors. Since ultrasound scattering is two to three orders of magnitude weaker than optical scattering in biological tissues, photoacoustic tomography is capable of providing significantly improved yet depth-independent spatial resolution over alloptical imaging techniques such as diffusion optical tomography (DOT). Thus far, PAT has shown its potential to detect breast cancer, to assess vascular and skin diseases, to monitor epilepsy in sma ll animals, to image finger joints, to sound out fluorescent proteins, and to evaluat e exogenous contrast agents in molecular imaging.8-24 Recent studies further demonstrated t he possibility for PAT to image hard tissues such as bones and associated soft tissues.25-27 In these studies, rat tail joint or cadaver human finger joints were imaged and gold nanoparticles were used as contrast agent to help diagnose / m onitor RA disease. While simplified 2-D PAT detection geomet ries and reconstruction algorithms have been used in most of the previous studi es, light scattering and pressure wave propagation in tissue is inherently 3-D, and 2-D approximation to a real 3-D problem (2D reconstruction model with 2-D detecting geometry surrounding 3-D targets) will inevitably bring errors, blurring and missing st ructures in the reconstructed 2-D images. Besides that, the structures of the examined biological tissues remain unidentified in sections other than the detecting plane, whic h may be important in some biomedicine cases. For example, a finger joint from a human subject is highly irregular in its

PAGE 18

18 volumetric structure. Although it seemed t hat tissues around proximal interphalangeal (PIP) and distal interphalangeal (DIP) joints were well differentiated in reported 2-D cross-sectional images of a photoacousti cally examined cadaver human finger 25, some tissue structures shown in the obtained 2D cross sectional images were not in agreement with those from hist ological photographs due to t he collapse of the signals from the entire 3-D joint stru ctures into the 2-D cross-se ctional imaging planes. While RA disease can be diagnosed from the cross section images alone, it is known that the coronal and sagittal section im ages of the finger joints ar e necessary for OA diagnosis which is impossible for 2-D PAT. Thus, full vo lumetric image reconstruction is essential to capture the joint tissues completely and accurately for diagnosis of both OA and RA. Although studies on convent ional PAT have focused on obtaining high resolution anatomical structure of biol ogical tissues by imaging absorbed optical energy density (i.e., the product of the intrinsi c absorption coefficient and the external optical fluence or photon density), PAT is also capable of prov iding the intrinsic absorption coefficient distribution when a light tr ansport model is combine with conventional PAT method.28-34 Moreover, physiological / functional informa tion such as oxygen saturation or the concentration of hemoglobin, to which optical absorption is very sensitive, can also be effectively recovered in high resoluti on when multispectral light is used.35-42 These intrinsic tissue properties may be important for accurate diagnostic decision-making in early stage of diseases. 1.2 Overview of Osteoarthritis Osteoarthri tis (OA) is the most common degenerative joint disease, affecting tens of millions of Americans and involves over 500,000 individuals for joint replacement annually in the United States.43-44

PAGE 19

19 Primary OA commonly affects small join ts in the hand, including the distal interphalangeal (DIP) joints and the proximal interphalangeal (PIP), and large weightbearing joints in the hips and knees. Although knee and hip OA bear most responsibility for the burden of OA, hand OA may be a marker of a systemic predisposition toward OA; incidence study of OA fu rther revealed that the second DIP joint is most frequently affected (57% in women, 37% in men) among all the joints in the hand.45 Pain, stiffness, tenderness, joint enlar gement and limitation of motion are the typical symptoms and signs of OA, and the pathological featur es of OA include erosion of articular cartilage and associated bony changes.46 Synovial effusion and associated enhanced blood vessel growth may be occasi onally observed as well in OA joints, although it may not be as severe as that in jo ints affected by rheumatoid arthritis (RA). Thus far, clinical examination remain s the principal approach to OA diagnosis, relying on symptoms and signs OA patients su ffered as well as experience of arthritis physicians. Imaging techniques, includi ng standard x-ray radiography, computed tomography (CT), ultrasonography and magnet ic resonance imaging (MRI), have been conventionally used or investigated to vi sualize the physiologic/anatomic abnormalities in the joint cavity when OA has been established.47-54 However, all these available imaging techniques are used only as supplemental methods in OA diagnosis, especially in hand OA, either because they are insensitiv e to the abnormal changes in the OA joint cavity or because they are too costly to be used as a routine examination method. For example, plain radiographs are able to visu alize structural abnormality in the bone compartment (joint space narrowing and osteophyte formation) in joint cavity with high spatial resolution when OA has been well est ablished, however they are insensitive to

PAGE 20

20 changes in soft tissues (cart ilage, synovial fluid, etc) and therefore incapable of capturing the primar y features when OA is establishing in earlier stage. As such, sensitive and affordable im aging methods are urgently needed for the detection of OA, especially in early stages. Mo reover, the progress of effective imaging methods in OA diagnosis may at the same ti me accelerate the advancement of medical therapies for OA, which are currently effect ive to relieve OA symptoms or prevent the worsening of OA to a certain degree and ye t of limited effectiveness in modifying OA. Up to now, no disease-modifying OA drugs (DMOADs) has been approved. The development of the DMOADs may benefit greatly from the progress of imaging methods in OA diagnosis, where sensitive imaging te chniques can serve as surrogate markers for clinically meaningful outcomes to economic ally and efficiently validate the candidate DMOADs.43 New imaging techniques based on near-infra red (NIR) light, including pure-optical imaging techniques and photoacoustic tom ography (PAT), have been recently studied to image finger joints and to effe ctively detect joint diseases. 4-7,25, 55 The high color contrast (absorption coefficient) provided by diseased biological tissues (as high as 3:1 between tumor and normal tissues; > 40% increase in absorption for osteoarthritic cartilage compared to normal cartilage), has made NIR based imaging techniques a promising modality for early disease detection.23 While optical imaging techniques are able to detect the highly sensitive optic al absorption / scattering abnormalities associated with soft tissues (cartilage, synovial fl uid, etc) in OA and RA joints, its spatial resolution is relatively low (about 3~5mm). Compared to all-optical imaging techniques, PAT is able to visualize the same optical abs orption contrast with significantly improved

PAGE 21

21 spatial resolution (0.5mm or better, adjustable with ultrasound frequency) for deeptissue imaging. For example, tendon, aponeurosis, volar plate, subcutaneous tissue, phalanx and other tissues around DIP joints were well differ entiated in high resolution when a cadaver human finger was photoacoustically scanned in cross section with a 10MHz transducer. In this study, we develop a full 3-D PAT approach (high performance 3-D PAT reconstruction algorithm and a 3-D PAT system in spherical scanning geometry), aiming to study finger joints imaging and detect ost eoarthritis disease in the hand. The high performance 3-D PAT reconstruction algorithm is based on parallel computing technique and the finite element method, wh ich is validated phantom experiments. The feasibility, the optimal detector performanc e, and the optimal scanni ng geometry for 3-D photoacoustic finger joints imagi ng were pre-investigated with a series of finger-joint phantom experiments.56 A 3-D PAT system in a spheric al scanning geometry has then been constructed and optimized for in-vivo ex amination of human finger joints, and the distal interphalangeal (DIP) joint from a human volunteer was imaged.57 The performance of the 3-D PAT system has been further improved with an ultrasound detection array composed of eight 1 MHz PZT transducers and a 16-channel pulse/receive PCI board. The 3-D PAT syst em with eight detecting channels was carefully calibrated with controlled tissue phantom experiments, and is capable of completing a finger joint examination within 5 minutes (at single optical wavelength). Seven subjects (two OA patients and five heal thy controls) enrolled in our study from 2008 to 2009 were photoacoustically exam ined with our 8-channel 3-D PAT system.58 Major chromophore concentrations (HbO2, Hb) of in-vivo finger joints were further

PAGE 22

22 quantitatively imaged by using multispectr al 3-D PAT approach with six optical wavelengths from 730nm to 880nm.59

PAGE 23

23 Figure 1-1. Schematic of photoacoustic effect. Incident Pulsed NIR Light Absorber Generated Ultrasonic Pulses

PAGE 24

24 CHAPTER 2 RECONSTRCUTION ALGORITHM IN PAT Thus far, several algorithms have been implemented to effectively reconstruct photoacoustic images from measured acous tic wa ves, such as backprojection, Fourier transform, P-transform, k-wave method and statistical approaches.60-64 While others rely on analytical solutions to the photoacousti c wave equation in a regularly shaped imaging domain and assume acoustical homogeneit y in biological media, finite element based algorithm seems provide unriv aled advantage to accommodate tissue heterogeneity and geometric irregularity as we ll as allow complex boundary conditions and source representations.65-67 Moreover, acoustic property in the biological tissues is able to be recovered simultaneously with optic al property by finite element based algorithm, when acoustic heterogeneity is tak en in consideration in photoacoustic wave equation. In the following section, we will intr oduce the finite element based algorithm in photoacoustic tomography in detail. 2.1 Finite Element Based Reco nstruction Algorithm in PAT For pulse mode PA generation propagating in soft tissues at room temperature, the thermal diffusion and electrostrictive effects can be ignored. Therefore, only considering the thermal expansion mechanism the acoustic field by the PA generation in tissue is described by the following wave equation: 2 2 221 (,)(,)P p rtHrt ctCt (2.1) Where c is the speed of acoustic wave in the tissue;(,) prt is the acoustic pressure; is the thermal expansion coefficient; Cp is the specific heat at constant pressure; and

PAGE 25

25 (,) H rt is the heating function defined as thermal energy per time and volume absorbed in the tissue. Assuming the laser light is a delta func tion in time and uniformly irradiated on the surface of the tissue, and the heating function (,) H rt can be written as the product of spatial absorption function () r and the temporal illuminatio n as the form of delta function() t we can write equation (2.1) as )()( ),( 12 2 2 2tr tC trp tcP (2.2) Taking the Fourier transform of the above equation (2.2) in the temporal domain, we have the following Helmholtz equation: 22() (,)(,)pcr PrkPrik C (2.3) Where k is the acoustic wave number defined as/ kc 2 f Although the acoustic speed c can be consider ed as a constant if we assume the biological tissues are acoustically homogenous, it is indeed a variable quality in realistic biological tissues. For variable acoustic spee d in biological tissues, we can re-write equation (2.3) as the following equation, in a reference of a constant acoustic speed in a coupling medium 22 0() (,)(1)(,)pcr PrkPrik C (2.4) Where 22 0/()1 ccr Choose j as the set of basis functions, and then the weighted weak form of equation (2.4) can be stated as:

PAGE 26

26 22 0 () (1)0jj VV pcr PkPdVikdV C (2.5) Expanding (,) Pr and () r as the sum of complex coefficients multiplied by the basis functions:1 N ii iPp ,1 N ii i ,1 N kk i the finite element discretization of equation (2.5) can be written as the following linear equations: A pb (2.6) Where 2 22 00 2() ()j ijjijikkjiji VVVA dVkdVkdV ds 0 1 N iiji V i Pc bikdV C T Npppp ,,,21 And the following absorbing boundar y conditions have been applied: 2 2 P PnP (2.7) Where 0 2 0 0/1 8/32/3 ki ki ik ; 0 2 0/1 2/ ki ki and is the angular coordinate at a radial position 2.1.1 Forward Modeling in Acoustically Homogenous Media Although the acoustic speed c is indeed a vari able qualit y in realistic biological tissues, we may assume that the biological tissues are acoustically homogenous for simplicity in some cases and consider the acoustic speed as a constant in interested biological tissues during the process of image reconstruction. In these cases, equations (2.6) can be simplified as the following equations

PAGE 27

27 A pb (2.8) Where 2 2 0 2()j ijjijiji VVAk d s 0 1 N iiji V i Pc bik C T Npppp ,,,21 In photoacoustic tomography with the assu mption of acoustical homogeneity in biological tissues, equation (2.8) serves as the forward model to compute the spatial variation of acoustic pressure wave (ext ernal observable) based on a given optical property distribution, or to compute the updat ed acoustic pressure wave in space based on a initial guess or iterative value of an optical pr operty distribution in image reconstruction. 2.1.2 Inverse Modeling and Nonlinear Optimiza tion in PAT Reconstruction In the inverse problem of PAT reconstruc tion, the optical absorption distribution (ultrasound source term in ultrasound propa gation equation) was estimated given the governing photoacoustic equatio ns, boundary conditions, and measurements of the pressure wave on the detecting locations (ext ernal observable). Gauss-Newton iterative scheme combined with Marquardt and Tikhonov regularizations has been pursued in our research to obtain stable inverse solution to the PA equation. This approach uses the hybrid regularizationsbased method to update an initial (guess) optical absorption distribution iteratively to minimize an objecti ve function composed of a weighted sum of the squared difference between computed and measured data for all frequencies. Unlike forward problem in PAT where the st iffness matrix A is well-conditioned and the solution is well-defined, inverse problem in PAT has to estimate the spatial tissue

PAGE 28

28 properties with limited observables on t he boundary and most commonly results in illposed situation, yielding meaningless soluti ons in image reconstruction. In ill-posed situation, solutions in inverse problem may be non-unique, non-existent or unstable (a tiny perturbation in the measurement result s in large differences in the solution).68 To improve the conditioning of the problem and enable a numer ical solution in inverse modeling, Tikhonov regularizat ion was used to find an opt imized solution in inverse problem by minimizing the followi ng modified objective function 22()()cmFpp (2.9) Where the second term is an added pena lizing factor with a weighting in a Euclidean distance referenced to and is set equal to the parameter distribution at the previous iteration () i for simplicity in LevenbergMarquardt algorithm. At the thiiteration, we assume that the Taylors expansion of the objective function ()F around () i can be expressed as () ()()()()2() ()()()()() 2i iiii FH FFF (2.10) Where ()()iF and ()()i FH are the first and second order derivative At the minimum of the objective function()F we have stationary point where the following identity is satisfied () 0dF d (2.11) Substituting (2.10) into (2.11) and ignoring the derivatives higher than second order, we have the following Gauss-Newton iterative scheme 1 (1)()()()()()iiii FHF (2.12)

PAGE 29

29 Consider the regularized objecti ve function given in equation (2.9), where() i we have the first and second order der ivative of the objective function()F ()() ()2()2()c cmidp Fpp d ()() ()22cc Fdpdp H dd (2.13) Combining the above first and second or der derivative of the objective function()F with the Gauss-Newt on iterative scheme (2.12), we have the iterative equation (2)()TmcJJJTI pp (2.14) Where the Jacobian matrix J is defined over total M detecting locations as 111 12 222 12 12 N N MMM N p pp p pp J p pp 12,,,T cccc Mpppp 12,,,T mmmm Mpppp (1)() 12,,,T ii N (2.15) Since p is a linear function of from the linear equation set (2.8), we can calculate jkJby applying partial di fferentiation on equation (2.8) with respect tok

PAGE 30

30 ()()ijj ijj kkApB 0 ijj S i Pc B ikdS C (2.16) ijAandij B are independent onk and jkij sincej and k are independent variable if j k, therefore we have ijjkijjkik A JBB (2.17) In photoacoustic tomography, the final equation (2.14) with Jacobian matrix determined in equation (2.17) is the core of the in verse computation to update the optical property distri bution from an initial es timate. In each iterativ e step, it compares the measured with the solved observables, whic h is the acoustic pressure wave for different frequencies obtained via the forward model (2.8) and current es timate of the optical property distri bution, and solves the inverse problem by nonlinear optimization equation (2.14) to obtain a new estimate of th e optical property distribution. This procedure is iteratively r epeated until optimal least-squares fit of computed and measured acoustic data is reached. 2.1.3 Regularization in PAT Reconstruction The purpos e of introducing the Tikhonov regularization on the Hessian matrix H is to ensure the diagonal dominance and enable a stable numerical solution in inverse procedure. However, the opt imal weighting factor is problem-specific and sometime quite difficult to determine. An effect ive approach for choosing optimal weighting factor is to iteratively modify its value based on the trace of the Hessian matrix and the relative least-square error rele at each iterative step: 69 ([])T reletraceJJ (2.18) Where is an empirical coefficient. We normally use0.5 in our PAT problem.

PAGE 31

31 Since the weighting factor relies on relewhich iteratively decreases, the value of and the regularization effect will be reduced at each iterative step. The variation of in each step allows more accurate image to be revealed as the iterations progress, while a blurred image is initially reconstructed with large at the start of the iterations. The empirical coefficient is carefully chosen so that the image can be accurately reconstructed with iterations while stability to a numerical solution is ensured at the same time. Although weighting factor defined as identify (2.18) in Tikhonov regularization enable the stable numerical solution in inve rs e problem, the Gau ss-Newton iterative scheme may experience problem in iterat ive convergence sinc e the approximation (2.10) of Gauss-Newton iterative scheme are iteratively valid only near the mi nimum (actual solution). If t he initial estimates are far from the actual solution, Gauss-Newton iterative scheme may have di fficulty in convergence. Other than Gauss-Newt on iterative scheme, gradient descent, also referred as steepest descent, relies barely on the initial es timates although it may exhibit oscillation around the minimum and the convergent process is slow,68 since the first order derivative of objective function is used in gradient descent instead of second order derivative used in Gauss-Newt on iterative scheme. The real-valued objective function ()F decreases fastest from a given esti mate along the direction of the negative gradient of ()F at the given estimate, and hence t he iterative algorithm for gradient descent is given by (1)(1)()()()iiiiF (2.19) Where() i is a positive value used as iterative step size.

PAGE 32

32 Considering ()()iF given by equation (2.13) from Tikhonov regularized objective function (2.9), we have the following iterat ive algorithm for gradient descent 1 ()() 2TmcJ pp (2.20) Comparing the gradient de scent method given in (2.20) and the Gauss-Newton iterative scheme given in (2.14), we may see some sim ilarity between them. With a hybrid technique proposed in Lev enberg-Mar quardt method, we can use a steering factor to switch between the Gauss-Newton iterative scheme and gradient descent method. The final update equat ion can be written as ( )()TmcJJ J TII pp (2.21) Where is a weighting factor used in Tikhonov regularization to relieve the ill-posedness of the inverse problem; is a steering factor introduced in Levenberg-Marquardt method to balance the convergent ability and the convergent speed. When the steering factor is very big, the termJJT can be ignored and equation (2.21) shows the iterative behavior of gr adient descent method shown in equation (2.20). When is small, then equation (2.21) reduces to equation (2.14) and exhibits the behavior of Gauss-Newton iterative scheme. So practically, the value of is chosen iteratively by the following rules: 1) The steering factor is initially set with a big value so that the initial estimates can be chosen with less caution with t he gradient descent behavior shown in (2.21). 2) In each i terative step, the steering factor reduces by a given ratio (such as 10) to accelerate the convergence proce ss if the objective function decreases. Otherwise, the steering factor is increased by a given ratio.

PAGE 33

33 2.1.4 Simulation in PAT Reconstruction The reconstruction algorithm developed for 3-D PAT is tested with simulated finger joints, where a cylindrically simulated finger with a 2.5mm thickness cartilage sandwiched between two bones is buried in background. The absorbed energy contrast between the simulated cartilage and the bones is assumed as 3:1, and the two bones are 10mm i n diameter. The mesh for algorithm testing is a cy lindrical mesh with total 14784 tetrahedral elements and 2907 nodes, which is generated from a 2-D triangle mesh slice layer by layer. Each layer of the total 17 layer ed cylindrical mesh has 308 elements and 171 nodes, with layer interval 1.25mm in the axial di rection (Z axis). The coronal section and the cross section of the reconstructed finger joint from simulated data are shown in Figure 2-1, where the actual shape of the simulated finger joint is outlined with black dash lines. From the reconstructed finger joint, we can see the bones and the cartilage are recovered with great quality. The structure and the size of the bones and the cartilage are in great agreement with the actual shape and size of the simulated bones and the cartilage. And the abs orbed energy contrast betwe en the simulated cartilage and the bones is also recovered with high accuracy. The results shown in this simulation validate the recons truction algorithm of 3-D PAT for further finger joint imaging. 2.2 Computing Strategies in PAT Reconstruction In finite element based PAT reconstruction, the major task is solving the forw ard equations (2.8), the Jacobian matrix (2.17) and the inverse equations (2.14), which are normally time / memo ry costly depending on the scales of the linear system of equation (the total number N of finite element nodes). As t he scale increase, the cost in

PAGE 34

34 computational memory and time grows dr amatically, which may go beyond the computation limitation of a single PC. Herein some efficient computing strategies have been adopted to bring the high resolution PAT in full potential. 2.2.1 Adjoin Sensitivity Method Normally only the pressure on the boundary nodes (total M nodes) is observable, therefore, in the calculation ofjkJ, only det(),1,2,, j ssM is useful for updating in the inverse equation (2.14). Considering the Jac obian matrix defined in equation (2.17), we have 1 ,, s ksjjksjijiksiikJJABSB (2.22) Where, s j is defined as,1sj fordet() j s otherwise,0sj ; and s iScan be calculated from s jsiijSA It is worth to note that 0ikB if node iand kbelong to different element as shown in equation (2.16). As a result, we only need sum up the contribution of the nodes in elements contain nodek, even though the sum index iin s iS varies from 1 toN. 2.2.2 Dual Meshing Scheme In finite element based reconstruction algorit hm, spatial sampling rate (mesh density) that is used to discretize the target region is critical and greatly determines the computational costs and the resolution of the reconstructed images. For example, the rapid spatial variation in the field of ac oustic pressure wave (megahertz frequencies) governed by the photoacoustic wave equation demands high mesh resolution (typically/10 is the ultrasound wavelength) with fi nely discretized mesh in order to

PAGE 35

35 maintain forward calculation accuracy of the acoustic pressure field in the media. However a finely discretized mesh with in creased spatial sampling rate, and hence the total numbers of mesh nodes, will inevitabl y results in dramatically growing computational cost, which almost grows squarely and cubically with the total mesh number in the forward calculation of t he acoustic pressure and the inverse image reconstruction of unknown absorbed energy prof iles (or absorption coefficient profiles), respectively. At the same time, th e unknown absorbed energy profiles to be reconstructed in PAT are relatively uni form and hence representable with fewer mesh nodes (degrees-of-freedom). To balance the accuracy required for the calculation of the acoustic pressure wave and the cost of computation, dual-mesh scheme seems an optimal solution in the finite element based PAT algorithm, where a dense mesh is used in the forward calculation of the acoustic pressure and a relatively coarse mesh is used in the inverse image reconstruction. 71 Implementation of the dual mesh sc heme affects two components of the reconstruction algorithm:71-73 (1) the forward calculation of the acoustic pressure field in the media at each iteration, where the updated absorbed ener gy profile to be used in the forward calculation is defined on the coar se mesh while the forward calculation is based on the fine mesh, and (2) calculation of the Jacobian matrix that is used to update the estimates of the absorbed energy profile during the inverse image reconstruction. To calculate the acoustic pressure field in the media with the updated absorbed energy defined on the coarse mesh, interpolation process is required to get the absorbed energy profile over the fine mesh from the coarse mesh. Assume that the thi

PAGE 36

36 fine mesh node is inside the coarse meshL, as shown in Figure 2-2(a) for a simplified 2-D case and Figure 2-2(b) for a 3D case, the absorption energy value i at the thifine mesh node can be interpolated by the following formula 1()p nnN iLLi nr (2.23) Where nL is the absorbed energy value at each node of the coarse meshL,pNis the number of the element nodes. 3pN for 2-D triangle mesh, and 4pN for 3-D tetrahedral mesh. For a certain pair of coarse/fine me sh, the interpolation coefficient ()nLir is always a constant, and therefor e can be pre-calculated and stor ed in a mesh related file before the reconstruction process, which sa ves the time of re-computing during each iteration step. With the dual mesh scheme, the elements of the Jacobian matrix will be calculated by 0, j ijikikk S i kPp c ABBikdS C (2.24) Where kis the node on the coarse mesh,k is the basis function centered on node kin this mesh, and the integration is still perfo rmed over the elements in the fine mesh. Since k and i are the basis functions defined over the coarse mesh and fine mesh, respectively, the integration kernelki is non zero only in the overlapping zone of the coarse mesh the coarse nodek belongs to, and the fine mesh the fine node i belongs to. Complications arise from the integrat ion when the fine mesh spans more than one coarse element. To simplifying this problem we can generate the fine mesh from the

PAGE 37

37 coarse mesh by splitting the coarse element s into fine elements so that each fine mesh resides entirely within one coarse mesh element as shown in Figure 2-2(c). 2.2.3 Partial Reconstructi on in Region of Interest Unlike optical tomography whic h is an inverse field problem to quantitatively reverse the optical scattering / absorption pr operties in the media by measuring the transported light arriving optical detectors from light sour ces, photoacoustic tomography is indeed an inverse source problem. In photoacoustic tomography, the goal is to trace the ultrasound sources (optical absorbers excited by near-infrared pulse) by measuring the ultrasound field on a boundary with pre-k nown ultrasound properties in the media. The particular feature of this imaging techni que gives us an inspiration that we may confine our inverse reconstruction in the region of interest that may contain the ultrasound sources, which coul d greatly reduce the computat ion cost in the inverse procedure. To confine our inverse reconstruction in t he region of interest, a slight adjustment on the inverse equation is requi red where inverse Equation (2.14) is rewritten as (2)()TmcJJJTI pp (2.25) Where the new Jacobian matrix J is defined as 111 12 222 12 12 K K MMM K p pp p pp J p pp

PAGE 38

38 The new Jacobian matrix J is therefore a first order derivative of the ultrasound pressure at each detecting location to t he unknown optical absorbed energy (the intensities of ultrasound sour ces) at each finite elem ent node within the region of interest, which includes total K finite element nodes. The partial reconstruction is tested by simu lation, where a cylindrical target (10mm in diameter with 20mm height) is embedded in a cylindrical background. The contrast between the target and the background is 7:1. The section image along the axis of the cylindrical target from the reconstructed 3D images with full volume reconstruction is shown in Figure 2-3(a), where total 2907 nodes in the whole volume are involved in the inverse reconstruction. At the same time, th e cylindrical target is reconstructed with partial reconstruction in region of interest, wh ich is a small cylindrical volume containing the cylindrical target. The region of inte rest include total 801 nodes and is outlined with dash black lines in the section image along the ax is of the cylindrical target as shown in Figure 2-3(b), which is a 2-D section im age from the reconstructed 3-D images with partial reconstruction in the region of inte rest. As observed from Figure 2-3(b), the partial reconstruction in region of interest in our study obtained almost the same image shown in Figure 2-3(a) with full volume reconstruction in term of structure and quantitative values, while the computational cost in the partial reconstruction in region of interest is only around 1/47 of the computational cost in fu ll volume reconstruction. The partial reconstruction is further validated with phantom experiment as shown in Figure 2-4, where the full volume contains 1301 nodes and the region of interest outlined with black dash line includes only 229 nodes. Compared with the reconstructed images with full volume reconstruction (Figure 2-4(a)), images recons tructed with partial

PAGE 39

39 reconstruction in region of interest main tain the same quality and quantitative results (Figure 2-4(b)) while the computational cost in the partial reconstruction in region of interest drops to around 1/32. 2.2.4 Parallel Computing Technique In some cases where a fine mesh both in the forward calculation and inverse reconstruction is required for high resolution PAT images, the computational c ost may go beyond the limitation of a single PC. For exam ple, computational memory in tens of gigabytes is required to store the matrix A in forward equations, the Jacobian matrix J and the Hessian matrix JJTto image a blood vessel in sub100m spatial resolution. To overcome the computational cost in high re solution PAT (as well as 3-D PAT, where the problem scale grows cubically with the spatial resolution instead of squarely in 2-D case), parallel computing technique has to be utilized. In a parallel computing scheme, the memory request and the computation tasks for a computation assignment can be spread on multiple processors as shown in Fi gure 2-5. In the follo wing, we will give a detailed description on the parallel computin g based PAT algorithm. In finite element based PAT algorithm, the forward calculation of acoustic field by Equation (2.8) and the determination of Jacobian matrix J by Equation (2.17) and (2.22) are based on the forw ard modeled matrix A at each frequencyi ordered from1 to K with totalK frequency elements. The matrix A in forward equations is usually a symmetric sparse matrix and stored in memory by banded storage or compressed storage strategy, requiring less memory than the Jacobian matrix J and the Hessian matrixJJT, which is usually a full matrix. Herein we spread the forward calculation of acoustic field and the determination of Jacobi an matrix on distributed processors (total

PAGE 40

40 number is1Q ) by the frequency elementi where the whole matrix A under certain frequency elements is stored in the memory of the specifically assigned processor. As shown in flow chart of the high performanc e photoacoustic tomography (Figure 2-6), the matrix A at frequencyi is stored in the memory of a processor determined by(1,) M odiL whereLis the average number of frequency elements on each processor (total frequency elements K divided by total processors1Q ). In each processor, the forward modeled acoustic fi eld as well as the elements of Jacobian matrix at the frequencies associated with the specified processor is calculated independently, after which the Jacobian matrix J over the whole frequency range is assembled and stored dispersedly in the memori es of the distributed processors. In this way of parallel computing and storage, the computation load and storage is evenly assigned among the processors, and the co mputation assignment is maximally parallelized since each processor can r un the assigned task independently without communication with other processors. Jacobian matrix J and the Hessian matrix JJT(denoted asH) is full matrix, requiring distributed storage over the processors. A wrapping storage over the columns of the Jacobian matrix J is used to evenly store the Jacobian sub-matrix over the processors and further efficient ly calculate the Hessian matrixHwith minimally mutual communication between the processors. The wrapping storage is described below. (1,1)jJCPUMyidifMyidModjQ (2.26) Where 1( 1 )(),,,,,,1,,T jjMjMjMKjJJJJJjN

PAGE 41

41 With the wrapping storage of the Jacobian matrix J, the Hessian matrixHcan be calculated and stored by 12 (1,1)T jjjHJJJJCPUMyidifMyidModjQ (2.27) Where 1,,,1,,T jjjjHHHjN Because the Hessian matrix H is symmetric, we store only the upper triangle of the Hessian matrix spread over the proc essors. In the calculation of elementijH, communication and data exchange may be involved among the distributed processors, except that the vector jJ and jJ are both stored in the memory of the same processor. Since the Hessian matrix is st ored dispersedly, the inverse reconstruction by Equation (2.14) requires a parallel solv ing method, w here parallelized Cholesky decomposit ion is used in our high performance PAT. With the high performance PAT, simulated blood vessels with different diameters (0.6mm, 0.28mm and 0.14 mm) c an be accurately reconstructed as shown Figure 27(a), using a finite elem ent mesh with 15200 triangle elements and 7721 nodes. The 0.14mm diameter vessel (red dash line in Fi gure 2-7(a)) is quantitatively recovered as 0.133mm in diameter based on FWHM analysis. Crossed hairs (~60 m in diameter) in a phantom experiment can also be imaged in a 0.1mm resolution wh en a finite element mesh with 14377 nodes is used, as shown in Figure 2-7(b).

PAGE 42

42 (a) (b) Figure 2-1. The reconstructed images for simulated finger joint (a) coronal section slice along X=0mm (b) cross section slice along Z=15mm. (a) (b) (c) Figure 2-2. Geometry of the absorbed energy interpolation at fine node i in 2-D case (a) and in 3-D case (b), and the calculation of the Jacobian matrix element in the effective integration zone (c). i L1 L2 L3 L4

PAGE 43

43 (a) (b) Figure 2-3. Slice along the axis of the cylindrical target from reconstructed 3-D images with full volume reconstruction (a), and partial reconstruction in region of interest (b). (a) (b) Figure 2-4. Reconstructed 2-D images of a ci rcular target in a phantom experiment with full volume reconstruction (a), and partial reconstruction in region of interest (b).

PAGE 44

44 Figure 2-5. Schematic of para llel computing with a combination of both shared memory and distributed memory. CPU CPU CPU CPU Memor y CPU CPU CPU CPU Memor y CPUCPU CPU CPU Memor y CPUCPU CPU CPU Memor y A Computation Assignment

PAGE 45

45 Figure 2-6. Flow chart of the high perfo rmance photoacoustic tomography based on parallel computing technique and finite element method. (1~)iiQLK (1~2)iiLL () r initial values (1~)iiL Solvingpby Apb 1~,1~si sjsiijSolvingSby SA sMiN Solvingpby Apb Solvingpby Apb 1~,1~si sjsiijSolvingSby SA sMiN CPU0 CPU1 CPUQ Fill sub-matrix 0J Fill sub-matrix 1J Fill sub-matrix QJ CPU communication and organize sub-matrix (0~)iJiQ Calculate error function ()()()()cmeprpr Is error small enough? Stop Yes No CPU communication and Fill Hessian sub-matrix Solve the inverse equati ons in parallel and update () r 1~,1~si sjsiijSolvingSby SA sMiN

PAGE 46

46 (a) (b) Figure 2-7. The reconstructed images for si mulated blood vessels (a) and crossed hairs in phantom experiment (b).

PAGE 47

47 CHAPTER 3 PRE-INVESTIGATION OF 3-D PAT SY ST EM WITH TISSUE PHANTOM STUDY In this section, we aim to systematic ally evaluate the possibility of 3-D photoacoustic imaging of the finger joints based on a series of tissue-like phantom experiments. The 3-D stru ctures and absorption coeffi cient images of the jointmimicking phantoms are obtained using finite element based PAT reconstruction algorithm and its enhancement. Among different 3-D scanning geometries, cylindrical scanning geometries (shown in Figure 3-1) adopted from DOT seem straightforward and are easiest to implement, however they give poor resolution along the cy lindrical axis due to the limited aperture size of a transducer. For example, a join t-like phantom with cartilage thickness in 1mm was scanned by a 1MHz transducer in a cyli ndrical scanning geometry shown in Figure 3-1(a), and the joint space is missing in the reconstructed images as shown in Figure 32. A transducer with small aperture size can theoretically offer enhanced resolution; however, such a transducer often gives poor signal-to-noise ratio (SNR) with added cost. Compared to a cylindrical scanning ge ometry, a spherical scanning geometry with a modest size of transducer aperture may provide an optimal solution for 3-D imaging as the spatial resolution with such a scanni ng depends little on the aperture size of the transducer, especially around the center of the imaging zone. In this study, phantom experiments based on both cylindrical scanning geometry and spherical scanning geometry are conducted for comparison, the results of which show that spherical scanning geometry provides obvious spatial resolution improvement over cylindrical scanning geometry and might be suitable fo r finger joint imaging. The spherical scanning geometry is adopted in further pha ntom experiments to systematically

PAGE 48

48 evaluate the differentiability in controlled phantoms with varied thickness and varied optical absorption values in cartilage. 3.1 Tissue Mimicking Phantom To evaluat e the reconstruction algorithm of 3-D PAT and the possibility of 3-D photoacoustic imaging of the finger joints with phantom experiments, solid tissue phantoms that well approxim ate optical characters of in-viv o finger joints are fabricated. The solid tissue phantoms used here are made of Intralipid, India ink, distilled water and agar powder, among which Intralipid acts as t he optical scatter, India ink supplies color contrast (optical absorber) and the agar (1%2%) serves as the coagulator. The jointlike tissue phantoms used in our ex periments are optically well-ch aracterized, in term of the optical absorption/scattering coefficients. The highly purified agar powder (A-7049, Sigma, USA) is dissolved in distilled wate r, and is liquefied afte r being boiled to the melting temperature of 95oC with a regular microwave oven. The solution has to be stirred continuously in order to obtain good uniformity in fabricated phantom while cooling down to 40oC, when the solution must be poured into a mould and left there for some time (normally about 10 minutes for 4 00mL solution) to r each a proper hardening and stable optical properties. The optimal tem perature for adding Intralipid and India ink to the agar solution is around 60oC, although it is non-critical between 80oC and 40oC. The detailed fabrication proc edure and the performance in c ontrolled optical properties of the tissue phantom have been fully described before.74 The two phantom bones fabricated in our study to mimick real finger bones had: a = 0.07mm-1, s =2.5mm-1 and a diameter of 6mm, and the cartilage between the bones had: a = 0.015~0.04mm-1, s =1.5mm-1 and a thickness of 1 or 2mm, as shown

PAGE 49

49 in Table 3-1. All the optical properties of the fabricated phantom are based on real finger joints. 3.2 System Description and Experiments The joint-lik e tissue phantom experiments in this study were conducted with a 3-D PAT imaging system as shown in Figure 3-3, which is capable to scan the phantoms along a cylindrical surface as well as along an equivalent spherical surface. The Ti: Sapphire laser (LOT IS, Minsk, Belarus) generates a pulsed beam with a pulse repetition rate of 10 Hz and a pulse wi dth <10ns. The wavelength of the laser was tunable from 600 to 950nm and was locked at 820nm in our experiments for deep medium penetration. After the reflection mirror, t he laser is then expanded by a lens and a ground glass so that the phant om can be uniformly illuminated by an area source of light. The laser energy used at t he phantom surface was about 10mJ/cm2, which is far below the safety standard of 22mJ/cm2. To detect the acoustic field generated by the pulsed laser, 1 MHz transducer (Valpey Fisher, Hopkinton, MA) was used. T he transducer has a bandwidth of 80% at 6dB and an aperture size of 6mm. To reduce the attenuation of the acoustic field, the transducer and the phantom were all immers ed in a water tank, in which the sound speed was specified as 1495m/s. During an experiment, the phantom was placed at the center of the rotary stage, and the transducer was attached to the ro tary stage via an arm having a length of ~20mm. To acquire the acoustic field al ong a cylindrically scanning surface, the transducer was stepped along the Z axis as well as rotated along an arc facing the center of the joint-like phantom in XY plane. In our experim ents, the arc scanning in the XY plane covered 300 around the phantom as shown in Figure 3-4, allowing the

PAGE 50

50 collection of signals at 50 loca tions with a step interval of 6; the linear scanning along the Z axis permitted data collection at 650 pos itions along the cylindrical surface using 13 steps with an interval of 1mm. For the spherical scanning mode, the arc scanning in the XY plane covered the same 300 around the phantom as the cylindr ical scanning, and the phantom was also rotated along the Y axis so that the tr ansducer was equivalently rotated around the phantom. The phantom was rotated 15 per st ep along the Y axis until a full spherical surface was covered, whic h allowed the data collecti on at 600 scanning positions. The detected acoustic signals were amplif ied and filtered by a Pulser/Receiver 5058PR (Panametrics, Waltham, MA ). The bandwidth of 5058PR ranges from 10kHz to 10MHz, and five scales from 0.01 to 1.0MHz are selectable in high pass filtering. The preamplifier integrated in 5058PR is in 30dB gain and the gain of the 5058PR itself is 40/60dB with attenuation tunable from 0dB to 80dB in 1dB steps. The Labview programming was used to control the entire data acquisition procedure. 3.3 Methods Our finite element based 3-D PAT rec onstruction algorit hm and its enhancement for extracting quantitative absorption includes two steps. The first is to obtain the 3-D images of absorbed optical energy densit y through our 3-D PAT reconstruction algorithm that is based on finite element so lution to the photoacoustic wave equation in frequency domain subject to the radiation or absorbing boundary conditions (BCs), which is described in detail in Chapter 2. The second step is to recover the optical absorption coefficient distribution using the photon diffusion model coupled with the absorbed optical energy density images obtained from the first step.

PAGE 51

51 In the first step, the 3-D PAT reconstructi on algorithm uses the regularized Newton iterative strategy to update an initial absor bed optical energy density distribution to minimize an object function, which is co mposed of a weighted sum of the squared difference between computed and measured acoustic data. To recover the optical absorption c oefficient from absorbed energy density, the photon diffusion equation as well as the Robin boundary conditions can be written in consideration of a, D(r)E(r)(r)(r)S(r) (3.1) ()()()()() DrErrnErr (3.2) Where()1/()aErr ,() Dr is the diffusion coefficient, ()1/3(()())asDrrr and ()sr is the reduced scattering coefficients, is a boundary condition coefficient related to the internal reflection at the boundary, and S(r) is the incident point or distributed source term. For the inverse computation, the so -called Tikhonov-regularization sets up a weighted term as well as a penalty term in or der to minimize the squared differences between computed and measured abs orbed energy density values, 2 2 0()()()[()()]coMinrrrrr LEE (3.3) Where L is the regularization ma trix or filter matrix, is the regularization parameter, ooooT 12N(,,...,) and ccccT 12N(,,...,) where o iis the absorbed energy density obtained from PAT, and c i is the absorbed energy dens ity computed from Equations. (3.1) and (3.2), for i =1, 2, N locations within the entire PAT reconstruction domain. The initial estimate of absorption coeffi cient can be updated based on iterative Newton method as follows with =1, 1()( )[()]TTToc EJJILLJ (3.4)

PAGE 52

52 In addition to the usual Tikhonov regularization, the PAT image (absorbed optical energy density) is used both as input data and as a priori structural information to regularize the solution so t hat the ill-posedness associated with such inversion can be reduced. The whole imaging zone is segm ented by Amira (Indeed-Visual Concepts GmbH, Berlin, Germany) into seve ral sub-regions, and the segmented a priori spatial information based on the PAT image can be in corporated into the iterative process using Laplacian-type f ilter matrix, L 1 1/, 0,ijifij LNNifijoneregion ifijdiffrentregion (3.5) Where NN is the total node number within one regi on or tissue; the optical absorption coefficient distribution is then reconstruct ed through the iterative procedures described by equation (3.1) and (3.3), where a priori spatial infor mation from the PAT image is incorporated in through the matrix L 3.4 Results and Discussion The 3-D images were reconstructed us ing a finite element mesh of 41323 tetrahedral elements and 7519 nodes which required around 120 minutes per iteration for a total of 10 iterations on a computation cluster with 10 processors. Each of the processors worked at 2.2GHz with 2GB memory. The comparison between the cy lindrical and s pherical scanning geometries can be made from the images shown in Figure 35a and Figure 3-5b for a typical joint-like phantom with 1mm-thick cartilage. Figure 35a display the coronal, sagittal and cross sections of the 3-D images fr om the reconstructed absorbed energy density distributions with the cylindrical scanning, while Figure 3-5b show the corresponding sectional images with the spherical sc anning geometry. To compare the reconstructed images

PAGE 53

53 with the phantom, the actual s hape of the phantom is outlined in Figure 3-5 (Dashed rectangular or circle). As can be observed fr om Figure 3-5a, the p hantom structures in the coronal section (XY plane) are well reconstructed, while both the structures in the sagittal and cross sections (YZ and XZ planes) are clearly distorted especially for the cross section. Here we see that the struct ures along the circular plane (XY plane) can be recovered accurately and that the stru ctures are erroneously stretched along the cylindrical axis (Z axis). The images presented in Figure 3-5b with the spherical scanning show clear improvement over that using the cylindrical scanning. We see that there is almost no distortion along the Z ax is especially for the cross section. These images with the spherical scanning are quantitatively accurate in terms of the recovered size, shape and location of the objects. Since the spherical scanning gives us t he best image quality for 3-D joint imaging, the remaining results were all obtained based on this scanning geometry. Figure 3-5c presents the coronal, sagittal and cross se ctions of the reconstructed 3-D absorbed energy density for a joint-like phantom with 2mm-thick cartilage. It is observed from Figure 3-5b and 3c that the st ructures along each section are well detected and that the thickness of cartilage for both 1mm and 2m m cases are correctly recovered. The recovered absorption coefficient images and t heir profiles along a transect in the sagittal section are displayed in Figure 3-6 and 3-7, re spectively. We see th at the reconstructed absorption coefficients of the bones and cartilage are quantitatively accurate relative to their actual values. By estimating the full width at half maximum (FWHM) of the absorption coefficient profile s, the recovered cartilage thickness was found to be

PAGE 54

54 around 1.3 and 2.3mm, which ar e in good agreement with the actual object size of 1 and 2mm used in the experiments. We also conducted experiments to evaluate t he ability of resolving different optical contrasts when the absorption coefficient of cartilage varied fr om 0.015 to 0.04 mm-1. The reconstructed results from these expe riments are shown in Figure 3-8, where Figure 3-8a to 3-6d display the coronal, sagi ttal and cross sections of the reconstructed 3-D absorbed energy density when the cartilage had an absorption coefficient of 0.015 mm-1, 0.025 mm-1, 0.03mm-1 and 0.04 mm-1, respectively. We can observe that the bones and cartilage are clearly distinguishable when the absorption contrast between the bones and cartilage is large (Figure 3-8a and 38b) and that it is in creasingly difficult for these two types of structures to be s eparated (Figure 3-8c and 3-8d). The recovered absorption coefficient images and their profiles along a trans ect in the sagittal section are displayed in Figure 3-9 and Figure 3-10. We note that the absorption coefficient of the cartilage is quantitatively recoveredi t was found to be around 0.016, 0.026, 0.031 and 0.043 mm-1, which are in very good agreement wit h the actual values of 0.015, 0.025, 0.03 and 0.04 mm-1. It is worth to note that homogeneous acoustic speed was assumed in this study which is certainly an approximation to the possible heterogeneous acoustic medium used in the experiments. Previous studies have shown that this assumption may generate blurring in the reconstructed im ages. However, Compared to the blurring effects from an inappropriate scanning geometry and detection c onfiguration, the blurring due to the assumption of homogeneous acoustic speed is significantly smaller, especially for larger dimension targets. In addition, the blu rring due to the assumption of

PAGE 55

55 homogeneous acoustic speed can be eliminated by the use of advanced reconstruction algorithms without such an assumption. In summary, in this study we have shown that a spherical scanning geometry is probably an optimal way for us to perform 3-D PAT imagi ng of the finger joints. Our results also indicate that the joint-like stru ctures and their absorption coefficients can be quantitatively reconstructed using our 3-D PAT approach.

PAGE 56

56 Table 3-1. Phantoms used in the pr e-investigation of 3-D PAT system No Category Optical absorption coefficient (mm-1) Optical scattering coefficient (mm-1) Thickness (mm) 1 Bone 0.07 2.5 N/A 2 Cartilage 0.015 1.5 1.0 3 Cartilage 0.015 1.5 2.0 4 Cartilage 0.025 1.5 1.0 5 Cartilage 0.03 1.5 1.0 6 Cartilage 0.04 1.5 1.0

PAGE 57

57 (a) (b) Figure 3-1. Schematic of two potential cylindrical scanning geometries of finger joints imaging in 3-D PAT. Figure 3-2. A slice through t he axis of the reconstruct ed finger joint scanned by a cylindrical geometry show n in Figure 3-1(b)

PAGE 58

58 Figure 3-3. Schematic of the three-dimensional photoacoustic imaging system. Figure 3-4. Schematic of the sca nning arc/path in the XY plane. PC Rotator Sta g e Mirro r Rotator Controller Lens Transducer Pulser & Receive r Phantom Rotator Transparent Linear Stage Z Y X Ti:Sa pp hire Y X

PAGE 59

59 Coronal Section Sagittal Section Cross Section (a) (b) (c) Figure 3-5. Reconstructed absorbed energy density images in the coronal (z=0mm), sagittal (x=0mm) and cross (y=4mm) sect ions for a joint-like phantom with 1mm thick cartilage in cylindrical (a) and spherical (b) scanning geometries, and with 2mm thick cartilage in s pherical scanning geometry (c).

PAGE 60

60 (a) (b) Figure 3-6. Reconstructed absorption coe fficient images (coronal section) for the phantom with 1mm (a) and 2 mm (b) thick cartilage. Figure 3-7. The recovered absorption coeffi cient distributions along a transect (x=0mm) for the images shown in Figure 3-6.

PAGE 61

61 Coronal Section Sagittal Section Cross Section (a) (b) (c) (d) Figure 3-8. Reconstructed absorbed energy density images in the coronal (z=0mm), sagittal (x=0mm) and cross (y=4mm) secti ons for a joint-like phantom where the absorption coefficient of cartilage was 0.015 (a), 0.025 (b), 0.03 (c) and 0.04 (d) mm-1.

PAGE 62

62 (a) (b) (c) (d) Figure 3-9. Reconstructed absorption coe fficient images (coronal section) for the phantom with an absorption coefficient of cartilage of 0.015 (a), 0.025 (b), 0.03 (c) and 0.04 mm-1 (d). Figure 3-10. The recovered absorption coe fficient distributions along a transect (x=0mm) for the images shown in Figure 3-9.

PAGE 63

63 CHAPTER 4 IMAGING FINGER JOINTS INVIVO WITH 3-D PAT IN A SPHERICAL SCAN In this section, we develop a 3-D PAT system in a spherical scanning geometry indicated in finger-like phantom study, and a ttempt to image a distal interphalangeal (DIP) joint in vivo, which is most frequently affected (57% in wo men, 37% in men) among all the joints in the osteoarthritic hand indicated in an incidence study of OA.45 In vivo joint imaging represents a challenge bec ause unlike the joint from a phantom or a cadaver finger, an in vivo joint has abundant blood vessels located along the skin which strongly absorb light. The strong absorption from the blood ve ssels and bones give rise to significantly reduced signal-to-noise ratio in an in vivo setting. The investigation in this section indicate that major features/components in the DIP finger joint as well as the absorption coefficient profiles of these joint structures can be successfully achieved invivo, by using Ti: Sapphire laser based phot oacoustic system in a spherical scanning geometry. 4.1 System Development and Description The joint-lik e phantom experiments conduct ed in Chapter 3 indicates that a 3-D PAT approach in an equivalently spheric al scanning geometry based on a 1 MHz transducer with an aperture si ze of 6mm seems potential to be applied for imaging human finger joints. Here, we develop a 3-D PAT system in a real spherical scanning geometry so that the ultrasound detector can be jointly controlled by two rotary stepper motors to move on a spherical detecting surface surrounding examined in-vivo finger joints, instead of rotating the examined targets in tissue phantom experiments for equivalently spherical scanning. The developed 3-D PAT system for in-vivo finger joints imaging in a spherical scanning geometry co nsists of a fiber guided near-infrared

PAGE 64

64 lighting subsystem, spherical scanning stepper motors, a ultrasound detector and an arrangement for hand and finger placem ent, as shown in Figure 4-1. 4.1.1 Fiber Guided Near Infrared Lighting Instead of using a reflection mirror to gu ild the pulsed light beam illuminating the examined s ubject from the top of the wate r tank, the pulsed light beam at 820nm generated from the Ti: S apphire laser is guided by an optic al fiber to the bottom of the water tank in in-vivo examination setting. With this illumination approach, there is enough room to optimize the light beam without interfering with the spherical scanning mechanical arms. The optical fiber used in our 3-D PAT allows more flexibility in the lighting approach of the examined finger join ts. Although a fixed area source of NIR light area source is used in our study to visualize the joint cavi ty, a rotary area source can also be delivered by our lighting system, whic h light approach may di fferentiate the skin and the inside structures in cross-secti on of finger joints when a high frequency transducer is used.25 The fiber guided light beam is expanded by a lens and a ground glass before reaching the finger so that the DIP joint can be uniformly illuminated by an area source of light. The laser energy is measured by a power meter on the surface a finger joint is placed so that the energy we use for the in-v ivo examination is well below the safety standard of 22mJ/cm2 and will not do any damage to human subjects. The laser energy is also measured after it comes out of the Ti: Sapphire laser to ensure the energy safety and good working condition of t he laser device itself. Before being used in-vivo, the lighting system is carefully adjus ted so that the light beam is centered and uniform for further finger joints examination. The adjusted li ght system is verified with phantoms as shown in Figure 4-2, where the cylindrical phantom is illuminate from the bottom surface and the line phantom as well as the pencil lead is illuminated from side.

PAGE 65

65 From the testing results show n in Figure 4-2, we see that the light beam is almost perfectly uniform and centered, and is ready for in-vivo examination. 4.1.2 Arrangement for Hand and Finger Placement A mechanical arrangement for hand and finger placement was designed and shown in Figure 4-3, w here the hand and the pr oximal end of the ex amined finger rest on the plastic holders cushioned by blue and pi nk sponges in Figure 4-3(a). The distal tip of the examined finger slips into a thin pl astic ring as shown in Figure 4-3(b). Multiple ring sizes are available for examined finger joints from different subjects so that the plastic ring well fits the individual finger jo int. Although the plastic ring is close to the examined DIP joint, it has a little impact on the light illumination and on the propagation of the ultrasound wave from the light heated DIP joint. The arrangement for hand and finger placement is screwed to the botto m surface of the water tank, and can be adjusted for the placement of little finger and index finger of both hands. If the finger joint is to be illuminated from the dorsa l side, the hand should face up resting on the plastic holder; otherwise, the hand should face down on the plastic holder, which might be a more comfortable pos ition for human subjects. 4.1.3 Signal Detecting in a Spherical Scanning Geometr y The acoustic field generated by the pulsed laser beam is detected by the same 1 MHz transducer used in phantom experimental studies in Chapter 3, which has a bandwidth of 80% at -6dB and an aperture si ze of 6mm. To reduce the attenuation of the acoustic field, the transducer and the phantom were all immersed in a water tank. As can be seen from Figure 41, the transducer was attached to the rotary stage R2 whose position was controlled by the ro tary stage R1. The collection of acoustic signals along a spherical scanning surface was realized through the combined rotations

PAGE 66

66 of rotary stages R1 and R2 (Note that R1 and R2 were attached to the same L-shape arm, see Figure 4-4). As show n in Figure 4-4, the examined finger joint was aligned in the direction of Y axis, and located at t he rotation center, O. In the scanning, the transducer was first rotated by R1 to a position along the circular locus L1 and then rotated by R2 along the circ ular locus L2 (Figure 4-4a). Th is process was repeated until the spherical scanning was completed. Si nce the finger/hand blocked the scanning of the transducer in certain angles, the sc anning along L1 and L2 covered only 240 and 252, respectively (Figure 4-4b and 4-4c). To run the transducer along the above desig ned spherical locus, not only the joint motion of the two rotary stages should be well controlled, but also the position of the transducer and the rotary cent ers of the two stepper motors should be well calibrated. Unlike the 2-D circular scanning shown in Fi gure 4-5(a) where the error of absolute starting point only results to slight rotation or slight translation of the reconstructed targets, mechanical errors in the spherical scanning as s hown in Figure 4-5(b) could lead to distorted scanning locus and big errors in the reconstruction because of model mismatch. To calibrate the mechanical parts for the spher ical scanning geometry, a bubble level gauge purchased from Sears and a ma rked paper cap for the transducer as shown in Figure 4-6 are used. By usi ng the bubble level gauge, we can make sure that the rotary axes of t he stage R1 and R2 are perpendicular to each other. With the center of the transducer marked on the paper cap, we will know that the circular trace of the transducer controlled by the rotary stage R2 is on the same plane with the rotary axis of the stage R1, since t he center of the transducer should be a stationary point at point S shown in Figure 4-4(c) when the rota ry stage R1 rotates. After the calibration,

PAGE 67

67 the transducer is able to run the designed spherical locus when it is jointly controlled by the two rotary stages. 4.1.4 Units Control and Data Acquisition Sy stem The received acoustic signal on each detec ting point along the spherical scanning surface is pre-amplified, am plified and filtered by the same Pulser/Receiver 5058PR used in phantom experimental studies in Chapt er 3, and further feeds to PCI A/D card (CS121000, Gage Applied Technologies, Lock port, IL, USA) where the received analog signal is converted into digital signal and stored for further signal processing and image reconstruction. The CS121000 is in a 12-bit resolution and its sampling rate is 100MS/s. A Labview program was used to control the entire data acquisition procedure. 4.2 In-vivo Experiments A female healthy volunteer entered our in -vivo study, and was photoacoustically examined by our developed 3D PAT system. The volunteer sat on a chair behind the water tank, the height of which is adjustable for comfortable positio ning of the examined subject. The low end of the arm and the exam ined hand were all immersed i n a water tank, which is filled with warm water so that the human subjects will not feel discomfort during the whole examination procedure. During the examinat ion, the finger joint was placed at the center of the two rotary stages, and the palmar side of the examined DIP finger joint faced up allowing the finger joint to be illuminated from the dorsal side of the finger. Our experience indi cates that this way of light illumination can give us maximized tissue penetration, resulting in clearly differ entiable joint cavity in the reconstructed images. The proximal end and the distal end of the examined finger are secured on the plastic holder to reduce the motion error dur ing the in-vivo examination. A complete

PAGE 68

68 scanning allowed the collecti on of signals at 387 positions along the spherical surface which took about 40 minutes. 4.3 Methods The collect ed photoacoustic signals in te mporal domain were converted into signals in frequency domain by Fourier Trans forming for our finite element based 3-D PAT reconstruction algorithm, which is based on finite element solution to the photoacoustic wave equation in frequency domain subject to the radiation or absorbing boundary conditions (BCs). The 3-D images of absorbed optical energy density were reconstructed through our 3-D PAT reconstr uction algorithm. The second step is to recover the optical absorption coefficient dist ribution by using a forward fitting procedure based on the photon diffusion model and the re constructed absorbed optical energy density images in the first step.29 The following photon diffusion equati on as along with the Robin boundary conditions can be used in the forward fitting pr ocedure to recover the optical absorption coefficient using the absorbed optical energy density reconstructed, ()()()()()aDrrrrSr (4.1) ()()Drnr (4.2) Where () r is the photon density,()()/()arrr ()ar is the optical absorption coefficient, () Dr is the diffusion coefficient, 1/(3())asD and s is the reduced scattering coefficients, is a boundary condition coeffici ent related to the internal reflection at the boundary, and () Sr is the incident point or distributed source term. To recover the optical absorption coefficient from the re constructed absorbed optical energy density from an initial distribution of ()ar the optical fluence ()r

PAGE 69

69 and the absorbed energy density()c are iteratively calculated through the photon diffusion equation and ()()c ar respectively. If the error between and ()c is not small, then ()ar is updated by ()()/()arrr and the above procedure is repeated until a small error between and ()c is reached, resulting in a stable quantitative distribution of()ar 4.4 Results and Discussion In this secti on, we present the in vivo 3-D reconstruction results and make observations based on these results. A coronal section of the reconstruct ed 3-D image (absorbed energy density) is shown in Figure 4-7a, while a si milar section of the DIP joint from a typical human finger by MRI is given in Figure 4-7b for comparis on. As can be observed from Figure 4-7a, the cartilage/joint space is well differentiated from the adjacent di stal phalanx (DP) and intermediate phalanx (IP). The recover ed absorption coefficient image and the absorption coefficient profile along a transect X=-2mm co rresponding to the coronal slice shown in Figure 4-7a are displayed in Figure 4-8a and 4-8b, respectively. We see that the absorption coefficient values of the phalanxes and cartilage/joint space are quantitatively consistent with that reported in the lit erature. By estimating the full width at half maximum (FWHM) of the absorption coe fficient profiles, the recovered thickness of cartilage/joint space was found to be 1.6m m, which is in agreement with the actual size. Figure 4-9a presents a cro ss section of the reconstructed 3-D image (absorbed energy density). Again, a similar cross sect ion of the DIP joint from a typical human finger by MRI is shown in Figure 4-9b for comparison. Comparing both Figure 4-9a and

PAGE 70

70 4-9b, it appears that several joint tissue types are visible includi ng phalanx (PX), lateral artery (LA) and median artery (MA), and tendon (TE). The recovered absorption coefficient image and the absorption coeffi cient profiles along two transects are displayed in Figure 4-10a to 4-10c. The firs t transect (line 1 in Figure 4-10a) goes through the phalanx and tendon, while the second one (line 2 in Figure 4-10a) crosses the tendon and two arteries. The recovered absorption coefficients of the phalanx, tendon, median and lateral arteries were found to be 0.07mm-1, 0.074mm-1, 0.064mm-1 and 0.058mm-1, respectively. Again, these values are generally cons istent with that reported in the literature. We note that blurs around the two blood vessels are noticeable (Figure 4-10a), largely due to the relatively low-frequency (1MHz) transducer used in this study. We believe the resolution can be enhanced with a wideband transducer having higher frequency (e.g., 510MHz) which will make the millimetersubmillimeter sized arteri es accurately imaged. To view the volumetric stru ctures of the tissues around t he joint, a series of cross and coronal section images are displayed in Figure 4-11 and 4-12, respectively. As shown in Figure 4-11a to 1-11f, the phalanx is clearly visible in all the cross section images from the proximal end (Y=-2.5mm) to the dist al end (Y=9mm). The tendon is differentiated from the phalan x (Figure 4-11a, 4-11e and 4-11f). The lateral artery is visible almost in all the cross section images except for that shown in Figure 4-11f at the distal end of the finger. The median artery is seen in Figure 4-11e and Figure 4-11f, and weakly visible in Figure 4-11d. Other joint structures including the phalanges (DP and IP) and the cartilage (CL) are clearly imaged from the coronal section images close to the dorsal side of the finger (Figur e 4-12a and 4-12b). Close to the pa lmar side of the finger,

PAGE 71

71 the tendon and two arteries begin to show up (Figure 4-12c), and the two arteries eventually become clearly visible (Figure 4-12d). Negative values are also observed in the reconstructed absorbed energy density im ages. It is likely due to the homogenous acoustic approximation (constant acoustic speed) to the actual heterogeneous acoustic media of the finger joints and limited signal-t o-noise ratio because of the strong light scattering of joint tissues. In summary, we have presented a 3-D PAT technique that is able to image finger joints in vivo. Although it seems impossible fo r PAT to provide image quality (in terms of spatial resolution) comparable to MRI for jo int or breast imaging, PAT is capable of obtaining absorption coefficient or functional information. In addition, PAT is portable and low in cost. While our curre nt experimental setup is not optimized, it allows us to demonstrate the possibility of 3-D in vivo joint imaging for the fi rst time. The results obtained indicate that major joint structures and their abs orption coefficients can be quantitatively reconstructed using our 3-D PAT approach.

PAGE 72

72 Figure 4-1. Schematic of the 3-D spherical scanning PAT system for finger joint imaging. L: lens; T: detector/tr ansducer; R1/R2: rotary stages. (a) (b) (b) Figure 4-2. Test of lighting beam with a cylindrical phantom (a), a line phantom (b), and a pencil lead (c). Z Y X Pulser & Receiver Rotator Controller Ti:Sa pp hir

PAGE 73

73 (a) (b) Figure 4-3. Photograph of the arrangement for hand and finger Placement. (a) (b) (c) Figure 4-4. Schematic of the spherical scanning geometry (a), circular locus L1 (b) and L2 (c). S

PAGE 74

74 (a) (b) Figure 4-5. Reconstruction mesh and possible model mismatch of 2-D circular scanning (a) and 3-D spherical scanning (b). (a) (b) Figure 4-6. Bubble level gauge (a) and homemade marker cap (b) for calibration of mechanical errors in the phot oacoustic spherical scanner. 1 1

PAGE 75

75 (a) (b) Figure 4-7. A selected coronal section from the reconstructed 3-D image (absorbed energy density) at Z=-5mm (a), and MRI (c oronal section) from a similar joint (b). DP: distal phalanx; IP: intermediate phalanx; DIP: distal interphalangeal joint. (a) (b) Figure 4-8. A selected coronal section from the recovered absorption coefficient image at Z=-5mm (a), and its profile along the cut line X=-2mm (b). DP IP DIP DP IP DIP

PAGE 76

76 (a) (b) Figure 4-9. A selected cross section from the reconstructed 3-D image (absorbed energy density) at Y=7mm (a), and MRI (cross section) fr om a similar joint (b). PX: Phalanx; MA: median artery; TE: tendon; LA: lateral artery. (a) (b) (c) Figure 4-10. A selected cross section from the recovered absor ption coefficient image at Y=7mm (a), and its profile along cut lines 1 (b) and 2 (c), respectively. PX LA MA TE PX LA MA TE

PAGE 77

77 (a) (d) (b) (e) (c) (f) Figure 4-11. Cross section images (absorbed energy density ) at Y=-2.5mm (a), Y=0mm (b), Y=3mm (c), Y=5mm (d), Y=7mm (e) and Y=9mm (f). PX: Phalanx; MA: median artery; TE: tendon; LA: lateral artery.

PAGE 78

78 (a) (c) (b) (d) Figure 4-12. Coronal secti on images (absorbed energy dens ity) at Z=-5mm (a), Z=3.5mm (b), Z=1mm (c) and Z=4.5mm (d). DP: distal phalanx; IP: intermediate phalanx; CL: cartilage; PX: Phalanx; MA: median artery; TE: tendon; LA: lateral artery. LA MA DP IP CL DP IP CL LA MA TE

PAGE 79

79 CHAPTER 5 CLINCAL EXAMINATION OF OA PATIENTS WITH 3-D PAT Osteoarthri tis (OA) is the most common degenerative joint disease, which affects tens of millions of Americans. Although knee and hip OA bear mo st responsibility for the burden of OA, hand OA may be a marker of a systemic predisposition toward OA; Incidence study of OA further revealed that the second DIP joint is most frequently affected (57% in women, 37% in men) among all the joints in the hand.45 Thus far, clinical examination remain s the principal approach of OA diagnosis, relying on symptoms and signs OA patients su ffered as well as experience of arthritis physicians. Imaging techniques, includi ng standard x-ray radiography, computed tomography (CT), ultrasonography and magnet ic resonance imaging (MRI), are used only as supplemental methods in OA diagnosi s, especially in hand OA, either because they are insensitive to the abnormal changes in the OA joint cavity or because they are too costly to be used as a routine exam ination method. As such, sensitive and affordable imaging methods are urgently needed for the detection of OA, especially in early stage. Moreover, the progress of effe ctive imaging methods in OA diagnosis may at the same time accelerate the advancemen t of medical therapies for OA, which are currently effective to relief OA symptoms or prevent the worsening of OA in certain degree and yet of limited effectiveness in modifying OA. Up to now, no approved disease-modifying OA drugs (DMOADs) are available, development of which may benefit greatly from the progr ess of imaging methods in OA diagnosis, where sensitive imaging techniques can serve as surrogate ma rkers for clinically meaningful outcomes to economically and efficiently validate the candidate DMOADs.43

PAGE 80

80 New imaging techniques based on near-infra red (NIR) light, including pure-optical imaging techniques and photoacoustic tomography (PAT), has been recently studied to image finger joints and to effectively detect joint diseases.4-7,25, 5557 While optical imaging techniques are able to detect the hi ghly sensitive optical absorption / scattering abnormalities associated with soft tissues (cart ilage, synovial fluid, etc) in OA and RA joints, its spatial resolution is relatively low (about 3~5mm). Compared to all-optical imaging techniques, PAT is able to visualize the same optical absorption contrast with significantly improved spatial resolution (0 .5mm or better, adjustable with ultrasound frequency) for deep-tissue imaging. In this section, we take the advantage of the high spatial resolution provided by PAT technique in imaging the optical absorpt ion property, continue previous study in finger joints imaging using 3-D PAT in a spherical scanning geometry, and present a pilot clinical study of detecting OA usi ng 3-D PAT on both osteoarthritis patients and healthy controls. 5.1 Enhancement in 3-D PAT Approach In the in-vivo finger joints imaging c onducted in Chapter 5, a single 1MHz transducer is used to receive the generated ultrasound signals, and then is feed into a single channel signal processing unit (preamplifier, amplifier, filter, and A/D conversion). The whole procedure of examination la sts about 40 minutes with single channel detecting system, which is too long for OA patients. To speed up the examination process, an ultrasound detection array together with multiple channel signal processing unit is designed, calibrated and tested with tissue-lik e phantom experiments. With multiple channel detection system the examination of a finger joint can be completed in 5 minutes, allowing the collect ion of photoacoustic signals at 360 positions surrounding

PAGE 81

81 the examined DIP joint. DIP fi nger joints from OA patie nts and healthy volunteers are further imaged with our enhanced 3D PAT system in a spherical scanning configuration, which are comprised of a pulsed NIR light source as same as the one in in-vivo finger joint imaging conducted in Chapter 4, a spherical scanning subsystem, and an ultrasound detection array and associated data acquisition system as shown in Figure 5-1. 5.1.1 Ultrasound Detection Array The ultrasound detection array (photograph s hown in Figure 5-1 and the design shown in Figure 5-2) is composed of ei ght 1 MHz transducers (Valpey Fisher, Hopkinton, MA) arranged uniformly along a 210o circular arc. The transducer has a bandwidth of 80% at -6dB and an aperture size of 6mm. The position and performance of each transducer in the ultrasound arra y has to be calibrated carefully before experimental study. Ideally, t he rotation of ultrasound detecting array in certain degrees in the plane where all the tr ansducer elements are arranged s hould result in the overlap of the positions and received acoustic si gnals between different transducer elements. For example, positions and received signals from transducer T1, T2 until T7 after being counterclockwise rotated in 30o (together with rotary stage R2) should coincide with those from transducer T2, T3 until T8 shown in Figure 5-2 before the ultrasound detecting array was rotated. In calibrating the system errors in mec hanical parts, the center of the ultrasound detector array must coincide with the rotati on axis of the rotary stage R2, and the plane where all the transducer elements are securely arranged should be perpendicular with the rotation axis of the rotary stage R2, so that the rotary motion of R2 will not cause evident position error between transducer el ements shown in Figure 5-3. After the

PAGE 82

82 calibration between the center of the ultrasound detecting arra y and the rotary center of the stage R2, the position of each transducer element has to be precisely aligned and secured in the ultrasound detecting array so that each transducer element is in the same distance from the center of the ultrasound detecting arra y. The position calibration of each transducer element is conducted by measuring the pulse-echo response of the each transducer element from a fixed reflec tion surface, as shown in Figure 5-4(a) where pulse-echo responses of fours sele cted transducer element s are measured in 50MHz sampling rate. From the responses show n in Figure 5-4 (a), location variations between the four transducer el ements are observed and require fine alignment so that the four pulse-echo responses are of the sa me flight of time (in the same dashed vertical line as shown in Figure 5-4(a)). T he alignment errors of the transducer elements can be limited within 0.2mm (about 7 sampling poi nts variation) by observing the above pulse-echo responses. The sensitivities and impulse responses of the transducer elements may also exhibit some differences, which result in the amplitude fluctuation and slight spectrum shift shown in Figur e 5-4(b). We measured the received acoustic signals of each element on the same location with a well-controlled phantom experiment and calculated the performance calibration coe fficients of all the transducer elements, which are used in the signal calibration of each experiment conducted with this system. The same method as employed in Chapter 4 is used to calibrate the altitude of the plane where all the transducer elements are securely arranged. With the marked cap shown in Figure 4-6(b), the transducer elements can be fi nely elevated on the same plane with the rotary axis of the stage R1 si nce the center of the transducer element at point F shown in Figure 5-2 should be a stat ionary point when the stage R1 rotates.

PAGE 83

83 5.1.2 Multiple Channel Data Acquisition The detected acoustic signals were feed into preamplifiers and converted to digital signals by an A/D board (PREAMP2-D and PCIA D1650, US Ultratek, Concord, CA), as shown in the block diagram of the multip le-channel 3-D PAT system (Figure 5-5). Two preamplifier sets (PREAMP2-D ) are used in our 3-D PAT system, each including four pre-amplifying channels. The preamplifier PREAMP2-D has a bandwidth from 20 KHz to 30MHz at -3dB, with a gain range from 10 to 60dB. The signal to Noise Ratio of the preamplifier is 67dB. The PCI card PCIAD1650 consists of 16 pulse/receiver channels, each of which is integrated with an amplifier (0 dB to 80 dB in 0.01dB increments) and a signal filter (low pass from 5.9MHz to 48MHz, and high pa ss from 16 KHz to 4.8MHz). The sampling rate of the PCIAD1650 is tunable from 390.625 KHz up to 50MHz (8 scales). The central values of the pulse/r eceiver channels have some variation as shown in Figure 5-6(a), and require relative calibration between the detecting channels. The multiple channels A/D board and the entir e examination proce dure are controlled by a Labview interface shown in Figure 5-6(b), which include the initiation of step motor controllers, procession of the spherical scanner(ultrasound detecting array), communication between the interface with the PCIAD1650, and signal display/storage. 5.1.3 Quantitative 3-D PAT with 3-D Photon Diffusion Model A forward fi tting procedure based on the 3-D photon diffusion model and the reconstructed 3-D absorbed optical energy density images in 3-D PAT is used to recover the three-dimensional im ages of optical absorption c oefficient distribution, which is important for further quantit ative analysis of the full join t. This forward fitting method requires no segmentation and is relatively easy to implement. 29To recover the optical absorption coefficient from the recons tructed absorbed optical energy density () r from

PAGE 84

84 an initial distribution of()ar the same updating strategy is used in the 3-D forward fitting procedure as the one used in 2-D forward fitting procedure. 29 Briefly, the optical fluence() r which relates with ()ar and () r by the identity()()/()arrr and the absorbed energy density()()cr are iteratively calcul ated through the 3-D photon diffusion equation and()()()()c arrr respectively. If the error between () r and ()()cr is not small, then ()ar is updated by ()()/()arrr and the above procedure is repeated until a small error between () r and ()()cr is reached, resulting in a stable quantitative distribution of()ar 5.1.4 Tissue-like Phantom Test The 3-D quantitative PAT (qPAT) system is validated on tissue mimicking phantoms before the c linical study was performed, to ens ure quantitatively accurate in vivo results. The joint-l ike tissue phantoms used in our experiments were optically wellcharacterized. The two bones had: absorption coefficient a = 0.07mm-1, reduced scattering coefficient s =2.5mm-1 and a diameter of 6mm, and the cartilage between the bones had: a = 0.01 or 0.03mm-1, s =1.5mm-1 and a thickness of 1.5mm. The optical absorption / scattering coefficients of the joint-like phantoms well approximate the corresponding coefficients of in-vivo finge r joints. Figure 5-7(a) and 5-7(b) show the reconstructed phantom joints in three-dimensional view, and in sagittal and coronal sections with our 3-D qPAT. The distributio ns of the recovered absorption coefficient along a transect crossing the bones and cart ilage for the images shown in Figure 57(a) and 5-7(b) are shown in Figure 5-7(c). The average absorption coefficient of the joint cavity when the cartilage had abs orption coefficient of 0.01 or 0.03mm-1 was

PAGE 85

85 found to be 0.016mm-1 and 0.030mm-1, respectively, while the absorption coefficient of the bones was recovered as 0.067mm-1 and 0.070 mm-1, respectively, for the two phantom cases. The recovered cartilage thicknesses were found to be 1.5mm and 1.7mm, which is in agreement with the exact size. 5.2 Patients and Clinical Data Seven subj ects (two OA patients and five healthy subjects) were enrolled in the current study between 2008 and 2009. All part icipants (white female; mean age 61 years, range 45 years) prov ided informed consent as part of the protocol approved by the Institutional Review Board of Universi ty of Florida. Participants were recruited from the Rheumatology Clinic at the Shands Health Center of the University of Florida (UF), Gainesville, FL USA. O ne distal interphalangeal (DIP) fi nger joint, joint II at the left hand, from each subject was ex amined clinically and photoacoustically in this study, because it was mostly vulnerable to OA dis ease and at the same time easily accessible with our current 3-D phot oacoustic scanning system. Clinical examination of each patient was performed independently by a single rheumatologist (E.S. Sobel). Patients with OA were identified by clinical history and main clinical features, including symptoms (predominantly pain and stiffness), functional impairment and signs (joint enlargement and redness). The healthy controls had no known OA or other joint diseases. In photoacoustic examination, the examined DIP finger joints were placed at the center of the two rotary stages R1 and R2 and the palmar side of the examined DIP finger joint faced up, allowing the finger joint to be illuminated from the dorsal side of the finger. Both the examined finger and the ultr asound array were all immersed in a water tank filled with warm water for acoustic couplin g. The proximal end and the distal end of

PAGE 86

86 the examined finger were secured on the pl astic holder. A complete scanning allowed the collection of signals at 360 positions along the spherical surface which took less than 5 minutes. Finite element based reconstruction algor ithm of 3-D quantitative photoacoustic tomography (qPAT) (i.e., 3-D conventional reconstruction algorithm coupled with 3-D photon diffusion equation) was used to quantitat ively recover the optical absorption coefficient images of the DIP joints. Fu ll-width at 30% maximum method (FW30%M) was used in pure-optical imaging modalities to quantitatively analyze the recovered finger joints, which is able to compensate the inadequate spatial resolution of pureoptical imaging method in finger joints imaging. While PAT is able to provide high resolution in imaging the joint cavity (1~2 mm) and hence the interfaces between bones and cartilages, it seems that PAT can only im age the synovial fluid (~0.5mm) in poor resolution. Herein, we combine the FW30%M crit eria and the traditional full-width at half maximum (FWHM) method to quantitatively analyze the recovered joints in our study, which enables us to get the best quantitative results on different tissues (cartilage and synovial fluid) in the joint cavity. Data from six of the sev en recruited subjects (2 OA and 4 healthy) resulted in successful photoacoustic image reconstruction, while data from one healthy control was not useable due to the severe motion e rror of finger during the photoacoustic examination. 5.3 Results The 3-D photoacoustic images of the DIP joints from both OA patients and healthy controls were reconstructed using a finite element mesh of 7519 nodes and 41323 tetrahedral elements. Based on the recove red 3-D images, the optical absorption

PAGE 87

87 coefficient of different joint tissues (bone, cartilage and synovial flui d) and the structural size of joint cavity (cartilage and synovial fluid) are extracted. These parameters are used as indicators/markers for differentiation of OA from normal joints in this study. Consecutive coronal and sagittal slices fr om the recovered 3D images from the joint cavity as indicated in Figure 5-8 are presented here. The exac t locations of these selected slices may be slightly different from subject to subject. The coronal and sagittal slices for a normal joint are shown in Figur es 5-9(a) and 5-9(b). As can be observed from Figures 5-9(a) and 59(b), the bones are clearly de lineated from the adjacent tissues (cartilage and fluid) in the joint ca vity. While there is no sharp boundary between the bones and cartilage/fluid, the joint tissue/space is clearly identified in both the coronal and sagittal sections. Further observa tion finds that the sy novial fluid seems to be differentiable from the su rrounding cartilages in some slices (e.g., Z=7.5mm in the coronal section and X=-3.0mm in the sagitta l section). It is worth noting that the recovered size of the bones appears to be smaller than the actual anatomic size. We believe it is mostly due to the limited bandwidth of the ultrasound dete ctors we used in the photoacoustic examination where the loss of low frequency signals resulted in the reduced size of the recovered DIP joint. Figures 5-10(a) and 5-10(b) display selected coronal and sagittal slices from the recovered 3-D images for an OA joint. Again, the bones are clearly differentiable from the adjacent joint tissues (cart ilage and fluid), and the joint cavi ty is easily identified in both the coronal and sagittal sections. It is inte resting to note that th e spatial distribution of absorption coefficient appears to be heterogeneo us for the OA joint (Figures 5-10), while it is relatively hom ogeneous for the healthy joint (Figur es 5-9). This observation is

PAGE 88

88 consistent with the findings from the previous optical imaging study. Compared with the healthy joint (Figures 5-9), the joint spac e narrowing seems apparent for the OA joint (Figures 5-10). The average absorption coeffi cient and structural size of joint tissues for all the subjects are calculated and presented in Tabl es 5-1. We note from Table 1 that the recovered absorption coefficient of the di seased synovial fluid ranges from 0.022 to 0.023mm-1, while it is below 0.017mm-1 (0.0089 to 0.017mm-1) for the normal fluid. Difference between OA joints and healthy cont rols is also apparent, in terms of the absorption coefficient of the synovial fluid in the joint cavity. The difference is more striking when the ratio of the abs orption coefficient of the sy novial fluid to that of the bone is considered (0.58~0.63 for OA joints verse 0.25~0.46 for normal joints). For the diseased cartilage, the absorpt ion coefficient (0.026~0.028mm-1) is generally larger than that of normal joints (0.015~0.024mm-1). The ratio of the absor ption coefficient of the cartilage to that of the bone further confirms the difference between OA and normal joints (0.72~0.76 for OA joints verse 0.43~0.61 for normal joints). We also note from Table 1 that the measur ed thickness of the synovial fluid for OA joints ranges from0.3 to 0.5m m, with a mean value of 0.4mm, whereas the range of this size for healthy joints is from 0.4 to 0. 9mm, with a mean a mean value of 0.67mm. For the diseased cartilage, the measured thickne ss ranges from 0.4 to 0.5 mm, with a mean value of 0.45mm, whereas the range of this size for healthy joints is from 0.5 to 0.8 mm, with a mean value of 0.675mm. Compared with the normal joints, the dimensional narrowing of both synovial fluid and cartilage is obvious for OA joints. A more obvious difference is observed in the thickness of th e joint cavity (synovial fluid and cartilage)

PAGE 89

89 between OA and normal joints (1.1~1.5mm fo r OA joints verse 1.7~2.5mm for healthy joints), where the observable joint space na rrowing for OA joints is in agreement with the findings from x-ray radiography. 5.4 Discussion In this study, for the first time we i n vivo imaged distal interphalangeal (DIP) joints and associated osteoarthritis in sub-mm re solution using our three-dimensional quantitative photoacoustic tomography (qPAT) approach. Apparent differences, in both the reconstructed size and optic al absorption coefficient of the joint cavity, are observed between osteoarthritic and normal joints. The phot oacoustic imaging of the joints offers a spatial resolution comparable to x-ray r adiography, and has the potential to be further enhanced when larger amount of transduc er elements and higher frequency transducers are used. The joint space narrowing observed photoacoustically is in agreement with the findings fr om x-ray radiography, and t he observed optical changes associated with diseased joints are in consistent with documented bi ological roots. For example, it is known that wit h the onset of OA, the synovia l membrane/fluid in articular cavity becomes increasingly turbid and has an enhanced vascularization.75 The increased turbidity and vascularization would accompany increased optical absorption coefficient (as well as optical scatteri ng coefficient) in the diseased synovial membrane/fluid. In fact, this optical increase can be as large as 100% at certain wavelengths in the NIR region. 76-77 In addition, there is incr easing evidence that OA is a disease involving a metabolic dysfunction of bone.7879 It is likely that this metabolic dysfunction of bone, often associated with high metabolism of subchondral bone, will cause changes in its tissue optical absorpt ion due to the changes in oxygenation, Hb, HbO2 and water content.

PAGE 90

90 An ultrasound array composed of eight 1M Hz transducers was used in the current study, which was sufficient to image the abnormalities in the joint cavity and the adjacent bones where cartilage erosion has been generally believed to be the major cause of OA disease. By using 1MHz tr ansducers, the joint cavity can be clearly identified and be measured. As shown previously in photoacoustic imaging of cadaver joints,25 PAT was capable of imaging other interesti ng joint structures as well, including volar plate, tendon and aponeurosis when a hi gh frequency transducer was used. In a recent study of whole-body mouse imaging using 3-D PAT,80 both relatively large organs such as liver, kidney and spleen, and sm all structures such as spinal cords, ovarian vessel, abdominal aorta and femora l veins were clearly imaged when multiscale wavelet filter was used along with a 64-element 3.1MHz ultrasound array. With these advancements, we believe that our 3-D PAT approach can be further improved to capture other interesting joint tissues such as ligaments, tendon and synovium in addition to the articular cavity and the bones This adds the potential for our 3-D PAT technique to detect OA that is initiated by these other joint tissues.81, 82 It is also worth to note that PAT im aging is governed by ultrasound propagation behavior in tissue while ultrasound signals are ex cited by optical pulses. For simplicity, in this study we have assumed an ac oustically homogenous medium in image reconstruction. The actual inhomogeneous nature of joint tissues will definitely affect the ultrasound propagation behavior in computation. In this case an advanced reconstruction algorithm considering acoustical heterogeneity is needed for improved image reconstruction. Such an advanced algori thm can also reconstruct tissue acoustic properties such as ultrasound speed in addition to tissue optical properties.38, 65 While

PAGE 91

91 the advanced reconstruction algorithm requires increased computation cost, the recovered tissue acoustic properties together with tissue optical properties and anatomical structure of joints tissues will offer more accura te detection of OA disease. For example, decrease of ultrasound speed and increase of ultrasound attenuation has been observed in osteoarthritic cartilage.83, 84 In summary, this study represents the first a ttempt to in vivo detect osteoarthritis in the finger joints using three-dimensiona l quantitative photoacoustic tomography. We show that apparent differences, in both the re constructed size of the joint space and the absorption coefficient of the joint cavity, exist between OA and normal joints. The results reported here suggest the possibility of 3-D qPAT as a potential clinical tool for early detection OA in the finger joints. As an em erging imaging techniqu e, the advantages of 3-D qPAT for early detection of OA lie in its capability of quantitat ively imaging optical absorption, acoustical properties, and physi ological/functional par ameters such as oxygenation and water content as well as rev ealing joint abnormality in high resolution. In addition, PAT is portable and low in cos t, and hence may play an important role in long term monitoring of OA.

PAGE 92

92 Figure 5-1. Photograph of the transducer array/finger interf ace of the three-dimensional photoacoustic tomography system used for the current study. Figure 5-2. Geometry of the detector array. T1~T8 are eight detectors, O is the center of the detector array. O T1 T2 T3 T4 T5 T6 T7 T8 F

PAGE 93

93 Figure 5-3. An example of possible mismat ches between the center of the detector array (marked as O) and the rotary cent er of the rotary stage R2 (marked as O). (a) (b) Figure 5-4. Alignment variations (a) and spec trum differences in impulse responses of four selected transducer elements in the ultrasound detecting array. O T1 T2 T3 T4 T5 T6 T7 T8 O D1 D2 D3 D4

PAGE 94

94 Figure 5-5. Block diagram of t he multiple-channel 3-D PAT system. PC Controlling Interface Rotator R2 Rotator R1 16 Channel PCIAD Card T1 T2 T3 T4 T5 T6 T i: Sapphire Synchronize Controller 1 Controller 2

PAGE 95

95 (a) (b) Figure 5-6. The central values of the 16 channels in PCIAD1650 card (a) and Labview control interface (b) in mu ltiple-channel 3-D PAT system.

PAGE 96

96 Coronal Section Sagittal Section 3-D view (a) (b) Figure 5-7. Reconstructed optical absorpt ion coefficient images in the coronal (z=0.0mm) and sagittal (x=0mm) sections as well as in 3-D view for a jointlike phantom where the absor ption coefficient of cartilage was 0.01 (a) and 0.03mm-1 (b). (c) shows the recovered abs orption coefficient distributions along a transect crossing the bones and cartilage (x=0.5mm) for the images in coronal section (c)

PAGE 97

97 (a) (b) Figure 5-8. Photographic schemat ic of the coronal (a) and sagittal (b) sections for the finger joint imaging. Z=5.0mm Z=5.5mm Z=6.0mm Z=6.5mm Z=7.0mm Z=7.5mm (a) Coronal section (dorsal palm) X=-5.0mm X=-4.5mm X=-4.0mm X=-3.5mm X=-3.0mm X=-2.5mm (b) Sagittal section (Medial Lateral) Figure 5-9. Recovered images at selected coronal (a) / sagittal (b) planes for case 1 (Healthy). Y Z X Y X ZSection Line Section Line

PAGE 98

98 Z=-9.0mm Z=-8.5mm Z=-8 .0mm Z=-7.5mm Z=-7.0mm Z=-6.5mm (a) Coronal section (dorsal palm) X=-3.0mm X=-2.5mm X=-2.0mm X=-1.5mm X=-1.0mm X=-0.5mm (b) Sagittal section (Medial Lateral) Figure 5-10. Recovered images at selected co ronal (a) / sagittal (b ) planes for case 5 (OA).

PAGE 99

99 Table 5-1. Averaged absorption coefficient and size of joint tissues for 2 OA and 4 normal joints No Joint Tissues Average a (mm-1) Mean Thickness (mm) Joint Thickness (mm) () ()c a b a () () f a b a Bone 0.035 N/A Cartilage 0.015 0.5 1 Fluid 0.008 0.7 1.7 0.43 0.25 Bone 0.039 N/A Cartilage 0.024 0.7 2 Fluid 0.017 0.7 2.1 0.61 0.45 Bone 0.034 N/A Cartilage 0.020 0.8 3 Fluid 0.015 0.9 2.5 0.59 0.45 Bone 0.037 N/A Cartilage 0.023 0.7 Healthy 4 Fluid 0.017 0.4 1.8 0.61 0.46 Bone 0.036 N/A Cartilage 0.028 0.5 5 Fluid 0.023 0.5 1.5 0.76 0.63 Bone 0.037 N/A Cartilage 0.026 0.4 OA 6 Fluid 0.021 0.3 1.1 0.72 0.58

PAGE 100

100 CHAPTER 6 IN-VIVO IMAGIING OF CHROMOPHORE CO NCENTRATIONS IN FINGER JOINTS WITH MULTISPECT RAL 3-D PAT Since last decade, near-infrared spec troscopy and NIR based multispectral imaging techniques (such as diffuse optic al tomography) has emerged as promising methods to noninvasively image or meas ure tissue chromophores and associated hemodynamics in living tissues. In those techniques, concentration distributions of physiologically / functionally significant tissue chromophores (such as oxy-hemoglobin deoxy-hemoglobin, water, lip id) and their variations due to disease development, metabolic changes, cerebral hemorrhage, and brain activity can be measured or imaged by taking the advantage of the absorption spectra of different tissue chromophores.85-94 Thus far, studies has shown that the m easurement of these ti ssue chromophores and their physiological / functional changes are of great significance both in basic science (such as neuroscience, molecular imaging of small animals) and clinical applications (such as breast cancer detection and classifica tion). For example, the concentrations of oxy-hemoglobin and deoxy-hemoglobin in breast tumors are found to be higher than normal tissues.96 In recently years, multispectral photoacoustic tomography (PAT) has been investigated aiming to quantify these same tissue chromophores and the oxygenation for molecular imaging and functional imaging with high spatial resolution, which can be as high as sub100m for several mm penetration depths in highly resolved PAT images.3542 Although the challenges remain since the conventional PAT images are the distributions of absorbed optical energy densit y and not actual images of the absorption coefficients that can be immediately r eady for quantitative spectroscopic analysis,35

PAGE 101

101 several models have been proposed and developed to quantitatively resolve the tissue chromophores with multispectral PAT.3640 For example, maps of optical absorption coefficients can be recovered when the traditional PAT is couple with a photon diffusion model under certain assumptions, such as pre-known scattering property or wavelength dependency of the scattering proper ty in the medium. The recovered maps of optical absorption coefficients in mu lti-wavelengths further go through spectroscopic analysis for quantitatively resolved concentrations of the tissue chromophores. Some groups made further attempts to in -vivo measure the concentrations of oxy-hemoglobin and deoxy-hemoglobin within blood vessels in a human finger,41 and to image the hemoglobin and oxygenati on in a rat brain,42 although more accurate results require more accurate models of multispectral PAT. In this section, our goal is to advance the application of multispectral PAT to invivo imaging of chromophore concentrations in finger joints. Although previous studies focused on 2-D modeled multispectral PAT, a mu ltispectral 3-D PAT approach is used in our study, which is more accurate and is essential to capture the full volumetric structures in high resolved PAT images with less artifacts and blurring than an approach with 2-D approximation, particu larly for tissues with high structural heterogeneity (such as finger joints). Optical wavelengths from 730nm to 880nm are used in our multispectral 3-D PAT approach since studies have indicated that mult ispectral light with wavelengths located on the two sides of 805nm is required to distinguish the concentrations of oxy-hemoglobin and deoxyhemoglobin from the combined optical absorption values.97

PAGE 102

102 6.1 Methods The algorit hm used in our multispectral 3D PAT is a fitting method, which includes three steps. The first step is to obtain the 3-D images of absorbed optical energy density (,) r in each optical wavelength through our 3-D finite element based PAT reconstruction algorithm. The second step is to recover the 3-D optical absorption coefficient distribution for eac h specific optical wavelength by using a forward fitting procedure based on the reconstructed 3-D absorbed optical energy density images (,) r in the first step and the following 3-D photon diffusion equation at each specific optical wavelength (,)(,)(,)(,)(,)aDrrrrSr (6.1) (,)(,)(,)Drrnr (6.2) Where (,) r is the photon density, (,)ar is the optical absorption coefficient,(,) Dr is the diffusion coefficient, (,)1/(3((,)(,)))asDrrr and (,)sr is the reduced scattering coefficients, is a boundary condition coeffici ent related to the internal reflection at the boundary, and (,) Sr is the incident point or distributed source term. To recover the optical absorption coefficient from the re constructed absorbed optical energy density,(,) r the forward fitting procedure starts from an estimated distribution of (,)ar which is an optimized initial based on the results of an searching scheme for optimal initial parameter s. The reduced scattering coefficient(,)sr is regarded as pre-known constant for simplicity. The optical fluence (,)r and the absorbed energy density()(,)cr are iteratively calculat ed through the 3-D photon diffusion equation and()(,)(,)(,)c arrr respectively. If the error between

PAGE 103

103 (,) r and ()(,)cr is not small for a specific optical wavelength, then the current(,)ar is updated by (,)(,)/(,)arrr and the above procedure is repeated until a small error between (,) r and ()(,)cr is reached, resulting in a stable quantitativ e distribution of(,)ar A third step is to recover the concentrations of major tissue chromophores (oxyhemoglobin, deoxy-hemoglobin, and water) from the resolved optical absorption coefficients maps at each wavelength, bas ed on the absorption spec tra of these tissue chromophores and the following Beers law in consideration of the optical absorption contributions from L chromophores. 1(,)()()L aii ircr (6.3) Where ()icr is the concentration of thichromophore in the unit of molar(/) M moleL; ()i is the absorption extinction coefficient of the thichromophore at wavelength The chromophores number2 L for the recovery of two tissue chromophores (oxyhemoglobin, deoxy-hemoglobin) in our study, which comprise the major contributions to the optical absorption in the wavelength region we used. To recover the concentrations of thes e tissue chromophores, an inverse equation is used, which is written as the following ()()aJJcrJr TTI (6.4) Where ()ijjiJ ; the elements () cr and ()ar are ()jcr and(,)air respectively. is a regularization factor to stabilize the inverse solution.

PAGE 104

104 6.2 In-vivo Experiments The in-vivo experiments are cond ucted by our 3-D PAT system described in detail in Chapter 5, which are comprised of a pulsed NIR light source with wavelength tunable from 600nm to 950nm, a spherical scanning subsystem, an ultrasound detection array and associated data acquisition system. Six opt ical wavelengths from 730nm to 880nm (730, 760, 805, 825, 850 and 880nm) are chosen for the multispectral 3-D PAT scanner, which may give the improved distinguishability between HbO2 and Hb. Both the ultrasound array and the examined finger are imme rsed in a water tank, which is filled with warm water, for minimized ultrasound a ttenuation. During the examination, the palmar side of the examined DIP finger joint faced up at the center of the two rotary stages R1 and R2, allowing the finger joint to be illuminated from the dorsal side of the finger. The collection of the acoustic si gnals at 360 locations along the spherical scanning surface surrounding the examined DIP joint requires about 5 minutes in each wavelength, and the entire procedure of mu ltispectral 3-D photoacoustic examination lasts for about 30 minutes. Again, the proxim al end and the distal end of the examined finger are secured on the plastic holder to reduce the motion error during the entire invivo examination procedure. To reduce the impact of possible motion errors on our multispectral PAT, our multispectral 3D PAT scanner operates in the following wavelength order: 730nm, 805nm, 850nm 880nm, 825nm, 760nm. A homemade 3-D positioning ruler (Figure 6-2) is used to che ck the location of the examined finger joint before and after the whole exam ination procedure, to elimi nate failed examinations with obvious finger motion during the examinatio n procedure. Seven subjects (two OA patients and five healthy subjects) entered the current study between 2008 and 2009. All participants (whi te female; mean age 61 years, range

PAGE 105

105 45 years) provided informed consent as part of the protocol approved by the Institutional Review Board of University of Florida. O ne distal interphalangeal (DIP) finger joint, joint II at the left hand, from each subject was examined clinically and photoacoustically in this study, because it wa s mostly vulnerable to OA disease and at the same time easily accessible with our cu rrent multispectral 3-D PAT scanner. Clinical examination of each patient was performed independently by a single rheumatologist (E.S. Sobel). The healthy controls had no known OA or other joint diseases. Four of the seven recruited subjects (t wo OA patients and two healthy controls) completed the whole procedure of 3-D multis pectral photoacoustic scanning resulted in successful photoacoustic image reconstructi on and quantitative m easurement of oxyhemoglobin and deoxy-hemoglobin in our study. 6.3 Results and Discussion The 3-D photoacoustic images of the DIP joints from all f our human subjects were reconstructed by using a spherical mesh of 7519 finite element nodes and 41323 tetrahedral elements. The optical absorption maps of reconstructed DIP joints are further recovered under each optical wave length by coupling reconstructed 3-D photoacoustic images (deposited optical absorption energy) of the examined DIP joints with 3-D photon diffusion model. Among the si x optical wavelengths where the DIP joints are reconstructed and their optical absorption maps are quantit atively recovered, three optical wavelengths are selected, based on the relatively structural consistencies of the recovered optical absorption maps, for further spectroscopic analysis to quantitatively resolve the concentrations of HbO2 and Hb. Figure 6-3(a) and 6-3(b) ar e the recovered optical absorption maps in coronal section and sagittal section, respectively, of an in-vivo finger joint at wavelengths from

PAGE 106

106 730nm to 880nm. As can be observed in Figure 6-3(a) and 6-3(b), t he joint tissue/space is differentiated from the adjacent bones in the coronal and/or sagittal sections at several optical wavelengths (760nm, 805nm, 850nm and 880nm). Relatively structural consistencies of the recovered optical absorption maps are observed between the wavelengths in 760nm, 805nm and 850nm, which ar e selected to further quantitatively resolve the tissue chromophore concentrations. The reason that the joint tissue/space is missing or fused with the adjacent bones at some optical wavelengths (such as 730nm and 825nm) are most likely due to motion er ror during the 30 minutes examination procedure, although an optimized wavelength or der is followed in our 3-D multispectral PAT scanning. We believe that the motion erro r and the structural inconsistencies of the recovered optical absorption maps at some wavelengths can be significantly reduced by using an ultrasound array integrated with much more ultrasound detect ors (for example, 64 detector elements as being used in high resolution whole-body mouse imaging).80 The quantitatively resolved concentrations of HbO2 and Hb from a typical OA finger joint, which is the same finger joint sh own in Figure 6-3, are displayed in Figure 64. The concentrations of HbO2 (Figure 6-4(a)) in different joint tissues are clearly differentiated in high resolution both in co ronal section and sagittal section, and it is relatively higher in the bones than in t he joint cavity. Figure 6-4(b) shows the concentrations of Hb both in coronal section and sagittal section, where the joint cavity is differentiated from the adjacent phalanges in the sagittal section and is barely identifiable in the coronal sect ion. We believe that the re latively poor quality of the resolved Hb concentration images in coronal se ction is most likely due to the relatively low quality image at the wavele ngth 760nm, where the absorpti on spectrum of Hb has a

PAGE 107

107 high peak and motion error may have happ ened during the 3-D multispectral PAT scanning. Figure 6-5 shows the quantitatively resolved concentrations of HbO2 and Hb from a typical healthy finger joint. Again, the structur e of joint space is clearly identified in the images of each resolved chromophore, and the concentrations of HbO2 and Hb in the joint cavity are found to be relatively lo wer than that in the adjacent phalanges. In comparison with the OA finger joint shown in Figure 6-4, the joint cavity in healthy subjects seems more distinctive. For all the four examined in-vivo finger joints, the average concentrations of oxyhemoglobin (HbO2), totoal hemoglobin (Total-Hb) and the oxygen saturation (SO2) and H2O) in the joint cavity and the adjacent bones are calculated and presented in Figure 6-6. As can be observed from Figure 66(a), the average concentrations of HbO2 and Total-Hb in the finger phalanges are in the range from 44 to 58 M and from 68 to 96 M respectively. The quantitatively re solved concentration levels of HbO2 and Total-Hb are in agreement with the measurements by diffuse near-Infrared spectroscopy. 98 Further observations indicate a difference in the levels of total-Hb between finger phalanges with OA diseas e and healthy phalanges, which are in range of 86~96 M and68~69 M respectively. The finding of enhanced hemoglobin levels in OA joints has the pathologic roots reveal ed by previous studies, where angiogenesis has been observed in diseased tissues. The ox ygen saturations of the four in-vivo finger joints are observed from 58% to 66% in the finger phalanges, while finger phalanges with OA disease seems to have lo wer oxygen saturations levels (55%~60%) than healthy phalanges (65%~66%). Hypoxia in osteoarthritic phalanges has been

PAGE 108

108 revealed and seemed to appear in pre-arthritic stage by previous studies of finger joints diseases. 99 Our results are in great agreement wit h previous findings of hypoxia in osteoarthritic phalanges. As shown in Figure 6-6(b), t he average concentrations of HbO2 and Total-Hb in the finger joint cavity are in the range from 16 to 29 M and from 27 to 53 M respectively. An obvious difference in the levels of total-Hb between joint cavities with OA disease and healthy joint cavities can be observed in Figure 6-6(b), which are in range of 46~53 M and27~42 M respectively. The oxygen saturations in the joint cavities are observed from 55% to 69%, and osteoarthritic fingers seem to have lower oxygen saturations levels (55%~59%) than healthy phalanges (59%~69%) in the joint cavities. Again, the enhanced hemoglobin le vels and dropped oxygen saturations in OA joint cavities are in agreement with the pathologic angiogenesis and hypoxia related with finger joint diseases. In summary, major chrom ophore concentrations (HbO2 and Hb) of finger joints have also been imaged in-vivo using multis pectral 3-D PAT approa ch with six optical wavelengths from 730nm to 880nm. The total hemoglobin levels of finger joints are found to range in68~96 M and27~53 M for joint phalanges and joint cavities, respectively. The oxygen saturations vary fr om 58% to 66% in joint phalanges, and from 55% to 69% in joint cavities. The E nhanced hemoglobin leve ls and dropped oxygen saturations in osteoarthritic phalanges and joint cavities have been observed, and the findings are in agreement with the pathologic angiogenesis and hypoxia related with finger joint diseases. The multispectral result s obtained in this study indicated that the multispectral 3PAT approach might be capable of detecting the angiogenesis and

PAGE 109

109 hypoxia in pre-arthritic stage of OA disease, and is promising as a potential clinical tool for early detection OA in the finger joints.

PAGE 110

110 Figure 6-1. Absorption spec tra of oxy-hemoglobin (HbO2), deoxy-hemoglobin (Hb) and water (H2O) in biological tissues. Figure 6-2. Homemade 3-D positioning ruler.

PAGE 111

111 730nm 760nm 805nm 825nm 850nm 880nm (a) 730nm 760nm 805nm 825nm 850nm 880nm (b) Figure 6-3. Recovered optical absorption maps in coronal section (a) and sagittal section (b) of an in-vivo finger joint for light source in six wavelengths.

PAGE 112

112 Z=-8mm X=-1mm (a) (b) Figure 6-4. Concentrations of oxy-hem oglobin (a) and deoxy-hemoglobin (b) in M and the water content (c) at coronal se ction (z=-8mm) and sagittal section (x=1mm) of examined finger joint from an OA patient.

PAGE 113

113 Z=2mm X=-0.5mm (a) (b) Figure 6-5. Concentrations of oxy-hem oglobin (a) and deoxy-hemoglobin (b) in M and the water content (c) at coronal se ction (z=2mm) and sagittal section (x=0.5mm) of examined finger jo int from a healthy subject

PAGE 114

114 (a) (b) Figure 6-6. Average Concentra tions of oxy-hemoglobin(HbO2), total hemoglobin (TotalHb) and the oxygen saturation (SO2) in finger joint bones (a) and in joint cavity (b).

PAGE 115

115 CHAPTER 7 CONCLUSION AND FUTHER WORK 7.1 Conclusion In this dissertation, we have investigated the feasibilities of 3-D P AT approach to image in-vivo finger joints and to detect the OA disease in the hand. Also multispectral 3-D PAT is studied to measure the majo r tissue chromophore concentrations (HbO2, Hb and H2O) of in-vivo finger joints. Results based on finger-like phantom experiments, in-viv o finger joint imaging and clinical cases study show that 3-D PAT in a spherical scanning geometry is able to image in-vivo finger joints in high resolution and is promising to detect the finger OA in the hand. Major anatomical structures of an ex amined distal inter phalangeal (DIP) joint along with the side arteries are clearly imaged, where joint space is well differentiated from surrounding finger phal anges. Seven subjects enrolled in our study are photoacoustically examined with our 8channel 3-D PAT system and apparent differences, in both the reconstructed size of the joint space and the absorption coefficient of the joint cavity, are obs erved between OA and normal joints. Major chromophore concentrations (HbO2, Hb) of in-vivo finger jo ints are quantitatively imaged by using 3-D multispectral PAT approach with six optical wavelengths from 730nm to 880nm, and the quantitative results are c onsistent with previous studies. As our studies indicate, the parallel com puting technique, 3-D PAT approach in a spherical scanning geometry, efficient optic al lighting method and calibrated multiple ultrasound detecting system are e ssential for effective finger joints reconstruction, OA detection and chromophores recovery. An ul trasound detecting system with much more detectors is expected to further improve the images quality and quant itative results.

PAGE 116

116 Although the image quality (in terms of spatial resolution) provided by our 3-D PAT seems not comparable to MR I, our 3-D PAT is able to obtain optical absorption coefficient and physiological / functional informat ion of the joint tissues (cartilage, fluid, phalanges, etc) in higher spatial resolution over pure optical techniques. In addition, PAT is portable and low in cost. While our curr ent 3-D PAT approach is not optimized, it allows us to demonstrate the po ssibility of 3-D in in-vivo finger joint imaging and clinical cases study of osteoarthritis for the first time, which sugges ts the possibility of 3-D PAT as a potential clinical tool for early detection of OA in the finger joints. 7.2 Future Studies 7.2.1 Evaluation of 3-D PAT in OA Detection with Small Clinical Samples The ultimate goal of this dissertation and our study is to detect OA joints from normal joints at an early stage. The advancement of the 3-D PAT bot h in the hardware and algorithms introduced in this dissertation serve the goal, and studies on limited clinical OA cases have already indicated t hat it is promising to detect the finger osteoarthritis in the hand by imaging abnormal ri se of the optical absorption coefficients in joint cavity (cartilage and synovial fluid) and the narrowing of the joint space with3-D PAT. However further evaluation with small clin ical samples plays a key role to validate our methods, where the specificity and sensitivit y of the 3-D PAT in the detection of OA patients will be investigated. 7.2.2 In-vivo Detection of Rheumatoid Arthritis in the Hand Rheumatoid arthritis is another common joint disease, which is characterized as inflamed joint. Although the 3-D PAT is in itially developed to image DIP joint for detection of OA disease in the hand, it is adjustable for PIP joint i maging and RA disease detection in the hand. Compared to the DIP joint imaging, PIP joint imaging

PAGE 117

117 may involve more ultrasound attenuation in the bones, resulting in further drop of signal to noise ratio (SNR). A PIP joint from a healthy human subject is photoacoustically examined by the adjusted 3D PAT system, and the reconstructed joint cavity seems differentiated from the surr ounding phalanges consecutively in both coronal section and sagittal section as shown in Figure 7-1. A photoacoustic study on another healthy human subject indicates that our 3-D PAT seems capable of identifying the synovial fluid in the joint, which may be interesting and important for rheumatoid arthritis detection since inflammation in the synovial fluid / synovial membrane is severe in rheumatoid joints. More in-vivo studies of PIP finger joints imaging is planned for further clinical cases study, and human subjects (RA patients and healthy controls) has being recruited. 7.2.3 Recover Acoustic Property Togethe r with Op tical Properties by PAT Approach Homogeneous acoustic speed was assumed in our study which is certainly an approximation to the heterogeneous acoustic medium (finger joints) in reality. Considering the acoustic heterogeneity in photoacoustic equation as shown in Equation (2.6), the variation of sound speed in the acoustically heterogeneous medium may be recovered. Figure 7-2 is the simulation geom etry of a testing case for advanced 2-D PAT with c onsideration of acoustic heterogeneity in photoacoustic equation, where the optical contrast between the targets and the background is 2:1 and the ultrasound contrast (ultrasound speed) is 1550m/s Vs.1 485.5m/s. The lower target has both optical contrast and ultrasound contra st, while the upper left target and upper right target have only optical contrast and ultrasound contrast respectively. As shown in Figure 7-3, strong crosstalk between the optic al contrast and ultrasound contrast is observed in first

PAGE 118

118 several iterations (Figure 7-3(a) and 7-3(b)), however the ultrasound speed and its contrast can be effectively recovered with more iteration in inverse reconstruction. We are planning to use advanced 3-D PAT to reco ver the ultrasound speed and its variation in in-vivo finger joints, which may further help quantify the OA disease in finger joints since the rise in ultrasound speed has been obs erved in OA finger joints from previous ultrasound studies of OA disease.

PAGE 119

119 (a) Sagittal section (medial lateral) (b) Coronal section (dorsal palm) Figure 7-1. Reconstructed images at selected sagittal (a) / coronal (b) planes for an invivo PIP joint. (a) (b) Figure 7-2. Simulation geometry to recove r ultrasound speed with PAT method. (a) is the simulation geometry of optical contrast s; (b) is the simulation geometry of ultrasound contrasts;

PAGE 120

120 A B (a) (b) (c) Figure 7-3. Reconstructed optical contrast images (A) and ultrasound speed images (B) at 1st iteration (a), 3rd iteration (b) and 5th iteration (c).

PAGE 121

121 LIST OF REFERENCES 1. H. Jiang, N. Iftimia, Y. Xu J. Eggert, L. Fajardo, and K. Klove, Near-infrared optical imaging of the breast with m odel-based reconstruction, Acad. Radiol. 9 186 (2002). 2. X. Gu, Q. Zhang, M. Bartlett, L. Schutz, L. Fajardo, and H. Jiang, Differentiation of cysts from solid tumors in the breas t with diffuse optical tomography, Acad. Radiol. 11, 53 (2004). 3. A. E. Cerussi, A. Berger, F. Bevilac qua, N. Shah, D. Jakubowski, J. Butler, R. Holcombe, and B. Tromberg, Sources of absorption and scattering contrast for nearinfrared optical mammography, Acad. Radiol. 8 211 (2001). 4. A. H. Hielscher, A. D. Klose, A. K. Scheel, B. MoaAnderson, M. Backhaus, U. J. Netz, et al, Sagittal laser optical tomography for imaging of rheumatoid finger joints, Phys. Med. Biol. 49, 1147 (2004). 5. A. Pifferi, A. Torricelli, P. Taroni, A. Bassi E. Chikoidze, E. Giambattistelli, Optical biopsy of bone tissue: a step toward the diagnosis of bone pathologies, J. Biomed. Opt. 9 474 (2004). 6. Z. Yuan, Q. Zhang, H. Jiang, D diffuse optical tomography imaging of osteoarthritis: Initial results in finger joints, J. Biomed. Opt. 12, 034001 (2007). 7. A. K. Scheel, M. Backhaus, A. D. Klose, B. Moa-Anderson, U. J. Netz, K. G. Hermann, et al, First clinical evaluatio n of sagittal laser optical tomography for detection of synovitis in arthritic finger joints, Ann. Rheum. Dis. 64, 239 (2005). 8. A. A. Oraevsky, A. A. Ka rabutov and S. V. Solomatin, Laser optoacoustic imaging of breast cancer in vivo, Proc. SPIE 4526 6 (2001). 9. S. Manohar, A. Kharine, J. C. G. van Hespen, W. Steenbergen and T. G. van Leeuwen, Photoacoustic mammography laboratory prototype: imaging of breast tissue phantoms, Phys. Med. Biol. 50, 2543(2005). 10. R. G. M. Kollkman, J. H. G. M. Klaessens, E. Hondebrink, J. C. W. Hopman, F. F. M. de Mul, W. Steenbergen, J. M. Thijssen and T. G. V. Leeuwen, Photoacoustic determination of blood vessel diameter, Phys. Med. Biol. 49, 4745(2004). 11. R. I. Siphanto, K. K. T humma, R. G. M. Kolkman, T. G. van Leeuwen, F. F. M. de mul, J.W. van Neck, L.N.A. van Adrichem, and W. Steenbergen, Serial noninvasive photoacoustic imaging of neovascula rization in tumor angiogenesis, Opt. Express 13, 8995 (2005). 12. S. Yang, D. Xing, Q. Zhou, L. Xiang and Y. Lao, Functional imaging of cerebrovascular activities in small animals using high-resolution photoacoustic tomography, Med. Phys. 34, 3294 (2007). 13. X. Wang, Y. Pang, G. Ku, X. Xie, G. Stoica, and L. V. Wang, Noninvasive laserinduced photoacoustic tomography for struct ural and functional in vivo imaging of the brain, Nat. Biotechnol. 21, 803 (2003).

PAGE 122

122 14. E. Z. Zhang, J. Laufer and P. Beard, Three dimensional photoacoustic imaging of vascular anatomy in small animals us ing an optical detection system, Proc. SPIE 6437 643710S (2007). 15. A. A. Karabutov, E. Savateeva, A. Or aevsky, Imaging of layered structures in biological tissues with opto-acoustic front surface transducer, Proc. SPIE 3601 284 (1999). 16. C. G. A. Hoelen, F. F. M. de Mul, R. Pongers, and A. Dekke r, Three-dimensional photoacoustic imaging of blood vessels in tissue, Opt. Lett. 23, 648 (1998). 17. J. A. Viator, G. Au, G. Palt auf, S. Jacques, S. Prahl, H. Ren, Z. Chen, J. Nelson, Clinical testing of a photoacoustic probe for port wine stain depth determination, Lasers Surg. Med. 30, 141 (2002). 18. S. Y. Emelianov, P. Li and M. ODonne ll, Photoacoustics for molecular imaging and therapy, Phys. Today. 05, 34 (2009). 19. D. Razansky, M. Distel, C. Vinegoni, R. Ma, N. Perri mon, R. W. Koster and V. Ntziachristos, Multispectral opto-acousti c tomography of deep-seated fluorescent proteins in vivo, Nat. Photonics 3 379 (2009). 20. H. F. Zhang, K. Maslov, G. Stoica and L. V. Wang, Functional photoacoustic microscopy for high-resolution and noninvasive in vivo imaging, Nat. Biotechnol. 24, 848 (2006). 21. H. F. Zhang, K. Maslov, G. Stoica and L. V. Wang, In vivo imag ing of subcutaneous structures using functional photoacoustic microscopy, Nat. Protoc. 2 797 (2007). 22. A. D. L. Zerda, C. Zavaleta, S. Ker en, et al., Carbon nanot ubes as photoacoustic molecular imaging agents in living mice, Nat. Nanotechnol. 3 557 (2008). 23. L. S. Boucharda, M. S. Anwarb, G. L. Liu, et al ., Picomolar sensitivity MRI and photoacoustic imaging of cobalt nanoparticles, PNAS 106 4085 (2009). 24. Q. Zhang, Z. Liu, P. R. Carney, Z. Yuan, H. Chen, S. N Roper and H. Jiang, Noninvasive imaging of epileptic seizures in vivo using photoacoustic tomography, Phys. Med. Biol. 53, 1921 (2008). 25. X. Wang, D. L. Chamberland, D. A. Jamadar, Noninvasive photoacoustic tomography of human peripheral joints toward diagnosis of infla mmatory arthritis, Opt. Lett. 32 3002 (2007). 26. D. L. Chamberland, X. Wang, B. J. Roessler, Photoacoustic tomography of carrageenan-induced arthriti s in a rat model, J. Biomed. Opt. 13, 011005 (2008). 27. D. L. Chamberland, A. Agarwal, N. Kotov, J. B. Fowlkes, P. L Carson and X. Wang, Photoacoustic tomography of joints aided by an Etanerc ept-conjugated gold nanoparticle contrast agent-an ex vi v o preliminary rat study, Nanotechnology 19, 095101 (2008). 28. J. Ripoll, V. Nt ziachristos, Quantitative point source photoacoustic inversion formulas for scattering and absorbing medium, Phys. Rev. E 71, 031912 (2005).

PAGE 123

123 29. Z. Yuan, H. Jiang, Quant itative photoacoustic tomograph y: recovery of optical absorption coefficient m ap of heterogeneous medium, Appl. Phys. Lett. 88, 231101 (2006). 30. B. Cox, S. Arridge, K. Kostli and P. Beard, D quantitative photoacoustic image reconstruction of absorption distributions in scattering medium using a simple iterative method, Appl. Opt. 45, 1866 (2006). 31. L. Yao, Y. Sun, and H. Jiang, Quantitative photoacoust ic tomography based on the radiative transfer equation, Opt. Lett. 34, 1765 (2009). 32. L. Yin, Q. Wang, Q. Z hang, H. Jiang, Tomographic imaging of absolute optical absorption coefficient in turbid medium using combing photoacoustic and diffusing light measurements, Opt. Lett. 32 2556 (2007). 33. L. Yao, Y. Sun, and H. Jiang, Transport-based quantitative photoacoustic tomography: Simulations and experiments, Phys. Med. Biol. 55, 1917 (2010). 34. Z. Yuan, Q. Wang, H. Jiang, Reconstruc tion of optical absorption coefficient maps of heterogeneous media by photoacoustic tomography coupled with diffusion equation based regularized Newton Method, Opt. Express 15, 18076 (2007). 35. B. T. Cox, J. G. L aufer and P. C. Beard, The challenges for quantitative photoacoustic imaging, Proc. SPIE 717713, 717713 (2009). 36. J. Laufer, D. Delpy, C. Elwell and P. Beard, Quantitative spatially resolved measurement of tissue chromophore concentrations using photoacoustic spectroscopy: application to the m easurement of blood oxygenation and haemoglobin concentration, Phys. Med. Biol. 52, 141 (2007). 37. B. T. Cox, S. R. Arridge, and P. C. Beard, Estimating chromophore distributions from multiwavelength photoacoustic images, J. Opt. Soc. Am. A 26 443 (2009). 38. Z. Yuan and H. Jiang, Simultaneous reco very of tissue physiological and acoustic properties and the criteria for heterogeneous media by quantitative photoacoustic tomography, Opt. Lett. 34, 1714 (2009). 39. Z. Yuan and H. Jiang, Quant itative photoacoustic tomography, Phil. Trans. R. Soc. A 367 3043 (2009). 40. J. Laufer, B. Cox, E. Zhang, and P. Beard, Quantitativ e determination of chromophore concentrations from 2D photoacoustic images using a nonlinear model-based inversion scheme, Appl. Opt. 10, 1219 (2010). 41. J. Laufer, E. Zhang, P. B eard, Quantitative in vivo m easurements of blood oxygen saturation using multiwavelength photoacoustic imaging, Proc. SPIE 6437 64371Z (2007). 42. X. Wang, X. Xie, G. Ku, and L. V. Wang, Noninv asive imaging of hemoglobin concentration and oxygenation in the rat br ain using high-resolution photoacoustic tomography, J. Biomed. Opt. 11, 024015 (2006).

PAGE 124

124 43. C. V. Oddis, New Pe rspective on Osteoarthritis, Am. J. Med. 100 2A, 10SS (1996). 44. R. W. Moskowitz, R. D. Altman, M. C. Hochberg, J. A. Buchwalter, V. M. Goldberg, Osteoarthritis: diagnose and me dical/surgical management, Lippincott Williams & Wilkins Philadelphia, USA (2007). 45. P. Sarzi-Puttini, M. A. Ci mmino, R. Scarpa, R. Caporali, F. Parazzini, A. Zaninelli, et al, Osteoarthritis: an over view of the disease and it s treatment strategies, Semin. Arthritis. Rheum. 35, 1 (2005). 46. S. B. Abramson, M. Attur, Y. Yazi ci, Prospect for disease modification in osteoarthritis, Nat. Clin. Prat. Rheumatol. 2 304 (2006). 47. P. Peloschek, G. Langs, M. Weber, J. Sa iler, M. Reisegger, H. Imhof, H. Bischol, F. Kainberger, An automatic model-based system for join t space measurements on hand radiographs: Initial experience, Radiology 243 855 (2007). 48. J. T. Sharp, Assessment of radiographic abnormalities in rheumatoid arthritis: what have we accomplished and where should we go from here, J. Rheumatol. 22, 1787 (1995). 49. T. Nishii, H. Tanak a, N. Sugano, H. Miki, M. Tak ao, H. Yoshikawa, Disorders of acetabular labrum and articular cartilage in hip hyspla sia: evaluation using isotropic high-resolution CT arthrography with sequential radial reformation, Osteoarthritis. Cart. 15, 251 (2007). 50. H. I. Keen, R. J. Wakefield, A. J. Gr ainger, E. M. A. Hens or, P. Emery, P. G. Conaghan, An Ultrasonographic Study of Osteoarthritis of the Hand: Synovitis and Its Relationship to Structur al Pathology and Symptoms, Arthritis. and Rheumatism. 59, 1756 (2008). 51. J. Chao and K. Kalunian, Ultrasonography in osteoarthri tis: recent advances and prospects for the future, Current Opinion in Rheumatology 20, 560 (2008). 52. J. A. Tyler, P. J. Watson, H. L. He rrod, M. Robson, L. D. Hall, Detection and monitoring of progressive cartilage degenerationof osteoarth ritic cartilage by MRI, Acta. Orthoped. Scand. Suppl 266 130 (1995). 53. H. Sugimoto, A. Takeda, K. Hyodoh, Early-stage rheumatoid arthritis: Prospective study of the effectiveness of MR imaging for diagnosis, Radiology 216 569 (2000). 54. H. Sugimoto, A. Takeda, J. Masuyama, M. Furuse, Early-stage rheumatoid arthritis: diagnostic accuracy of MR imaging, Radiology 198 185 (1996). 55. Z. Yuan, Q. Zhang, E. Sobel, H. Jiang, Tomographic x-ray guided threedimensional diffuse optical tomography of osteoarthritis in the finger joints, J. Biomed Opt 13, 044006 (2008). 56. Y. Sun, and H. Jiang, Quantitative three-dimensional photoacoustic tomography of the finger joints: Phantom studies in a spherical scanning geometry, Phys. Med. Biol. 54, 5457-5467 (2009)

PAGE 125

125 57. Y. Sun, E. Sobel and H. Jiang, Quant itative three-dimens ional photoacoustic tomography of the finger joints: an in-vivo study, J Biomed Opt 14, 064002 (2009). 58. Y. Sun, E. Sobel and H. Jiang, First assessment of three-dimensional quantitative photoacoustic tomography for in vivo detection of osteoarthritis in the finger joints, Radiology (under review). 59. Y. Sun, E. Sobel and H. Jiang, In vi vo imaging of chromophor e concentrations in finger joints by three-dimensional mu ltispectral photoacoustic tomography, Opt. Exp. (under submission). 60. K. P. Kostli, M. Frenz, H. Bebie, and H. P. Weber, T emporal backward projection of optoacoustic pressure transients us ing Fourier transform methods, Phys. Med. Biol. 46, 1863 (2001). 61. C. G. A. Hoelen, F. F. de Mul, R. Pongers, and A. Dekker, Three dimensional photoacoustic imaging of blood vessels in tissue, Opt. Lett. 23, 648 (1998). 62. P. Y. Liu, The P-transform an d photoacoustic image reconstruction, Phys. Med. Biol. 43, 667 (1998). 63. Y. V. Zhulina, Optical statistical approach to optoacoustic image reconstruction, Appl. Opt. 39, 5971 (2000). 64. K. P. Kostli, D. Frauchi ger, J. Niederhauser, G. Paltau f, H. Weber, and M. Frenz, Optoacoustic imaging using a three-dime nsional reconstruction algorithm, IEEE J. Sel. Top. Quantum Electron. 7 918 (2001). 65. H. Jiang, Z. Yuan, and X. Gu, Spatially varying optical and acoustic property reconstruction using finite-eleme nt-based photoacoustic tomography, J. Opt. Soc. Am. A. 23, 878 (2006). 66. Z. Yuan, Q. Zhang, H. Jiang, Simult aneous reconstruction of acoustic and optical properties of heterogeneous media by qua ntitative photoacoustic tomography, Opt. Exp. 14, 6749 (2006). 67. Z. Yuan and H. Jiang, Three-dimensio nal finite-element-based photoacoustic tomography: Reconstruction algorithm and simulations, Med. Phys. 34, 538 (2007). 68. Q. Fang, Computational Methods for Microwave Medical Imaging, PhD dissertation, Dartmouth College (2004). 69. K. D. Paulsen, P. M. Meaney, and L. Gilman, Alternative Br east Imaging: Four Model-Based Approaches, Springer Science+Business Media (2005). 70. K. D. Paulsen and H. Jiang, Spatially varying optical property reconstruction using a finite element diffusion equation approximation, Med. Phys. 22, 691 (1995). 71. K. D. Paulsen, A Dual Mesh Scheme for Finite Element Based Reconstruction Algorithms, IEEE Trans. Med. Imag. 14, 504 (1995).

PAGE 126

126 72. H. Jiang, K. D. Paulseny, U. L. Os terbergy and M. S. Pattersonz, Improved continuous light diffusion imag ing in singleand multi-tar get tissue-like phantoms, Phys. Med. Biol. 43, 675 (1998). 73. X. Gu, Y. Xu, and H. Jiang, Mesh-based enhancement schemes in diffuse optical tomography, Med. Phys. 30, 861 (2003). 74. R. Cubeddu, A. Pifferi, P. Taroni, A. Torricelli and G. Valentini, A solid tissue phantom for photon migration studies, Phys. Med. Biol. 42, 1971 (1997). 75. A. D. Meisel, P.G. Bullough, Atlas of Osteoarthritis, Gower Medical Publishing, New York (1984). 76. J. Beuthan, V. Prapavat, R. Naber, O. Minet G. Muller, Diagnos tic of inflammatory rheumatic diseases with photon density waves, Proc. SPIE 2676 43 (1996). 77. V. Prapavat, W. Runge, J. Mans, A. Krause, J. Beuthan, G. Muller, The development of a finger joint phantom for the optical simula tion of early inflammatory rheumatic changes, Biomedizinische. Technik. 42, 319 (1997). 78. J. Mansell, et al, Bone, not cartilage, should be the major focus in osteoarthritis, Nat. Clin. Prat. Rheumatol. 3 306 (2007). 79. L. Knott, A. Bailey, Collagen cross-links in mineralizing tissues: a review of their chemistry, function and clinical relevance, Bone 22, 181 (1998). 80. H. Brecht, R. Su, M. Fr onheiser, S. A. Ermilov, A. Conjusteau, A. A. Oraevsky, Whole-body three-dimensi onal optoacoustic tomography syst em for small animals, J. Biomed. Opt. 064007, 064007 (2009). 81. A. L. Tan, H. Toumi, M. Benjamin, A. J. Grainger, S. F. Tanner, P. D. Emery, D. McGonagle, Combined high-resolution magnetic resonance imaging and histological examination to explore the role of ligaments and tendons in the phenotypic expression of early hand osteoarthritis, Ann. Rheum. Dis. 65, 1267 1272 (2006). 82. K. D. Brandt, E. L. Radi n, P. A. Dieppe, L. van de Pu tte, Yet more evidence that osteoarthritis is not a cartilage disease, Ann. Rheum. Dis. 65, 1261 (2006). 83. S. L. Myers, K. Dines, D. A. Brandt, K. D. Brandt, M. E. Albrecht, Experimental assessment by high frequency ultrasound of articular cartilage thickness and osteoarthritic changes, J. Rheumatol. 22, 109 (1995). 84. G. A. Joiner, E. R. Bogoch, K. P. Pritzker, M. D. Bu schmann, A. Chevrier, F. S. Foster, High frequency acoustic parameter s of human and bovine articular cartilage following experimentally-induc ed matrix degradation, Ultrason. Imaging 23 106 116 (2001). 85. F. F. Jobsis, Noninvasive, infrared m onitoring of cerebral and myocardial oxygen sufficiency and circulatory parameters, Science 198 1264 (1977). 86. M. S. Patterson, B. Chance, and B. C. Wilson, T ime resolved reflectance and transmittance for the noninvasive measur ement of tissue optical properties, Appl. Opt. 28, 2331 (1989).

PAGE 127

127 87. M. Tamura, O. Hazeki, S. Nioka, B. Chance, and D. Smith,"Simultaneous Measurements of Tissue Oxygen Concentration and Energy State by Near-Infrared and Nuclear Magnetic Resonance Spectroscopy," Adv. Exp. Med. Biol. 222 359 (1988). 88. M. Cope and D. T. Delpy, "System fo r Long-Term Measurement of Cerebral Blood and Tissue Oxygenation on Newborn Infants by Near Infrared Transillumination," Med. Biol. Eng. Comput. 26, 289 (1988). 89. P. W. McCormick, M. Stew art, M. G. Goetting, M. Dujovny, G. Lewis, and J. I. Ausman, Noninvasive cerebral optical spectroscopy for monitoring cerebral oxygen delivery and hemodynamics, Crit. Care. Med. 19, 89 (1991). 90. S. P. Gopinath, C. S. Robertson, R. G. Grossman, and B. Chance, Nearinfrared spectroscopic localization of intracranial hematomas, J. Neurosurg. 79, 43 (1993). 91. T. Kato, A. Kamei, S. Takashima, and T. Ozaki, Human visual cortical function during photic stimulation monitoring by m eans of near-infrared spectroscopy, J. Cereb. Blood. Flow. Metab. 13, 516 (1993). 92. B. Chance, Z. Zhuang, C. UnAh, C. Al ter, and L. Lipton, Cognitionactivated lowfrequency modulation of light absorption in human brain, Proc. Natl. Acad. Sci. U.S.A. 90, 3770 (1993). 93. A. Villringer, J. Planck, C. Hock, L. Schleinkofer, and U. Dirnagl, Near infrared spectroscopy (NIRS): a new tool to study hemodynamic changes during activation of brain function in human adults, Neurosci. Lett. 154 101 (1993). 94. Y. Hoshi and M. Tamura, Dynamic mult ichannel near-infrared opticalimaging of human brain activity, J. Appl. Physiol. 75 1842 (1993). 95. A. Yodh, B. Chance, Spectro scopy and imaging with diffusing light, Phys. Today 48, 34 (1995). 96. S. Srinivasan, B. W. P ogue, S. D. Jiang, H. Dehghani, C. Kogel, S. Soho, J. J. Gibson, T. D. Tosteson, S. P. Poplack and K. D. Paulsen, Interpreting hemoglobin and water concentration, oxygen saturation, and scattering measured in vivo by near-infrared breast tomography, Proc. Natl. Acad. Sci. USA 100 12349 (2003). 97. Y. Yamashita, A. Maki, and H. Koizumi, Wavelength d ependence of the precision of noninvasive optical meas urement of oxy-, deoxy-, and total-hemoglobin concentration, Med. Phys. 28, 1108 (2001). 98. B. W. Pogue, S. P. Poplack, T. O. McBride, W. A. We lls, K. S. Osterman, U. L. Osterberg, K. D. Paulsen, Quantitative Hemoglobin To mography with Diffuse NearInfrared Spectroscopy: Pilot Results in the Breast, Radiology 218 261 (2001). 99. C. H. Jeon, J. K. Ahn, J. Y. Chai, H. J. Kim, E. K. Bae, S. H. Park, E. Y. Cho, H. S. Cha, K. S. Ahn, E. M. K oh, Hypoxia appears at prearthritic stage and shows colocalization with early synovial infla mmation in collagen induced arthritis, Clin. Exp. Rheumatol. 26, 646 (2008).

PAGE 128

128 BIOGRAPHICAL SKETCH Yao Sun received his B.S. in elec trical engineering from No rtheast Institute of Electric Power in July 2000 and his M.S. in physics, with specialty in acoustics, from Tsinghua University in 2003. After two years study in the PhD pr ogram of physics at Washington State University, he decided to devote himself to the society of biomedical engineering and entered The Univer sity of Florida in July 2005 for his PhD education in biomedical engineering. His research focused on hi gh performance photoacoustic tomography and its application to the imagi ng of finger joints and detection of osteoarthritis in the hand.