. nun
an!
A.

.

uvfidi.
: 3%..

 

XII
l. y

.651 .. .

,
.3 . .a
.1315"... .
94.5.19.

3».

gni'g
. I e e .31..
. rowan... . ,

‘ :33»...-
.a 3*.)"Nldh‘muutuuuu.
1.4;:

.v ..§..ufl.

0‘9 undilllbtrr l’i.t~Y.l
IAN! 1.5 I ‘ 37,- ‘A. "u\-
2 , .- lin i
6!

.l ‘IGAU... . .lil.
0"}... I!- a

\t‘vsio ‘
11311.25!“ R I! 31...!!!

AI

:.._.. s... .4 V. Fifi
3 it $1.1; A5,. .5 _
$§§a ..

x
If. 13:»; .31.... k...

 

Emggé,

E

'7') \_\x

l” LIBRARY 7
I Michigan §WI°
I... UniveISIlV

 

 

__.l

This is to certify that the
dissertation entitled

FORWARD AND INVERSE PROBLEMS IN NONINVASIVE

IMAGING TECHNIQUES

presented by

YIMING DENG

has been accepted towards fulfillment
of the requirements for the

Ph.D. degree in Electrical Engineering

 

saga/MM

Major Professor's Signatfire
6*? / <94 Am 7

Date

MSU is an Affinnative Action/Equal Opportunity Employer

 

PLACE IN RETURN BOX to remove this checkout from your record.
TO AVOID FINES return on or before date due.
MAY BE RECALLED with earlier due date if requested.

 

DATE DUE DATE DUE DATE DUE
0535121 592%111

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

5/08 K:IProj/Acc&Pres/ClRC/DateDue.indd

FORWARD AND INVERSE PROBLEMS IN NONINVASIVE
IMAGING TECHNIQUES

By

Yiming Deng

A DISSERTATION

Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of

DOCTOR OF PHILOSOPHY

Electrical Engineering
2009

ABSTRACT

FORWARD AND INVERSE PROBLEMS IN NONINVASIVE IMAGING
TECHNIQUES

By

Yiming Deng

Noninvasive imaging techniques are used in both engineering and clinical fields to
determine the state or internal conditions of structures and human body on the basis of
information contained in measured signals and images without the use of invasive
approaches. Development of the forward and inverse models for non-invasive
imaging not only helps to understand the physics and biological processes better, but
also enables development of efficient, low-cost and accurate algorithms for early
detection and diagnosis of anomalies. Recent developments in two noninvasive
imaging applications, Electromagnetic imaging including Magneto Optic (MO) and
Giant Magneto Resistive (GMR) methods and Positron Emission Tomography (PET)
imaging in biomedical field are covered in this dissertation. For MO and GMR
Imaging, the element free forward model as well as fast image reconstruction and
automated classification approaches are presented. In biomedical imaging application,
the PET modality is described and in particular the degradation of PET images due to
the patients’ thoracic (including cardiac and respiratory) motion is addressed. In
contrast to the other approaches for correcting motion artifacts in the image space, an
innovative framework for motion correction in the raw projection space (Sinogram) is

introduced based on the forward Monte Carlo model.

To my wife Dr. Xin Liu, daughter Joanna and parents

iii

 

 

ACKNOWLEDGEMENTS

I would like to express my sincere appreciation to my major advisor and mentor, Dr.
Lalita Udpa, for her guidance, support, and encouragement throughout my Ph.D.
study She also encouraged me in teaching and other academic activities. Dr. Lalita
Udpa inspires me so much through her own passion in research, excellence in
teaching and professionalism in service.

I would like to express my sincere gratitude and thanks to Dr. Satish S. Udpa for his
continuous support and guidance as my professor within the NDEL and as the Dean
in the College of Engineering. I would like to thank Dr. Kevin L. Berger for his
invaluable help during my biomedical research and continuous support during my
PhD study. I would like to thank Dr. David C. Zhu for his guidance in medical
imaging research and encouragement in my professional development. I would like to
thank Dr. Robert J. McGough for his support and help; especially his comprehensive
and interesting Medical Imaging and Acoustics classes which help me discover my
Passion in medical imaging and choose my future research areas afterwards. I would
also like to thank Dr. V. Mandrekar for being my committee member, and providing
valuable help in completing my study.

I would particularly like to thank my wife and colleague, Dr. Xin Liu for her
numerous support, encouragement and love. Finally, I would like to thank my
beloved parents, fiiends, Ms. Linda Clifford and the other colleagues of NDEL.
Thank you!

iv

TABLE OF CONTENTS

LIST OF TABLES ..................................................................................... viii
LIST OF FIGURES ..................................................................................... ix
CHAPTER 1. INTRODUCTION ................................................................. l
1.1 Motivation & Objectives .................................................................................... 1
1.2 Scope and Organization of the Dissertation ....................................................... 4
CHAPTER 2.NONINVASIVE IMAGING ................................................ 5
2.1 Introduction ........................................................................................................ 5
2.2 Electromagnetic imaging .................................................................................... 6
2.3 Biomedical imaging ........................................................................................... 7
2.3.1 Nuclear Medicine Imaging ............................................................................. 8
2.3.2 Positron Emission Tomography Imaging ..................................................... 10
2.4 Forward and inverse problems ......................................................................... 11
2.4.1 Forward problems in imaging ...................................................................... 12
2.4.2 Inverse problems in imaging ........................................................................ 13
2.4.3 Mathematical description of inverse problems ............................................ 13
CHAPTER 3.ELECTROMAGNETIC IMAGING ................................. 15
3.1 Introduction ...................................................................................................... 15
3.2 Fundamental of Electromagnetism ................................................................... 16
3.3 Principle of Magneto optic imaging ................................................................. 21
3.4 Principle of Giant Magneto Resistive (GMR) imaging ................................... 23
3.5 Issues in Electromagnetic Imaging .................................................................. 30
3.5.1 MO Imaging ................................................................................................. 30
3.5.2 GMR Imaging ............................................................................................... 31

CHAPTER 4. FORWARD MODELS FOR ELECTROMAGNETIC

MAGING....... ........................................................................................... 32
4.1 Introduction ...................................................................................................... 32
4.2 Finite element methods .................................................................................... 35
4.2.1 Sample Geometry and Mesh Generation ...................................................... 37
4.2.2 Interpolation/Shape functions and Global Matrix Assembly ....................... 38
4.2.3 Boundary and Interface Conditions .............................................................. 41
4.2.4 Matrix Solution ............................................................................................. 41
4.2.5 Post-Processing ............................................................................................. 42
4.2.6 Modeling M0 and GMR sensors with Induction Foil Excitation ................ 44
4-3 Forward model results ...................................................................................... 45

4.3.1 Validation Results ........................................................................................ 45

4.3.2 3D Results .................................................................................................... 55
4.4 Electromagnetic Imaging Parameters ............................................................... 64
4.4. 1 Frequency ..................................................................................................... 64
4.4.2 Sensor Liftoff ................................................................................................ 72
4.4.3 Conductivity of Top Layer ........................................................................... 78
4.4.4 Conductivity of Bottom Layer ..................................................................... 81
4.4.5 Conductivity of Fastener .............................................................................. 85
4.4.6 Fastener to Edge Distance ............................................................................ 87
4.4.7 Crack Dimension .......................................................................................... 94
4.4.8 Parameters Study Conclusions ..................................................................... 99

CHAPTER 5.1NVERSE PROBLEMS IN ELECTROMAGNETIC

 

IMAGING ........ - 101

5.1 Introduction .................................................................................................... 101

5.2 Data Analysis for M0 Imaging ...................................................................... 101
5.2.1 Automatic Rivet Detection ......................................................................... 102
5.2.2 Skewness Definition and automated classification .................................... 103

5.3 MO Image Classification ............................................................................... 110

5.4 Inverse Problems of MO imaging Conclusion ............................................... 119

5.5 Inversion and Data analysis for GMR imaging .............................................. 119
5.5.1 Optimization of Frequency Parameter in model predition - Skewness
Calculations ............................................................................................................ 120
5.5.2 Optimization of fi'equency parameter in model prediction-SNR Calculations

122

5.5.3 Magnitude Based Approach (MAG) data analysis for experimental data .126
5.5.4 Optimum Detection Angle (ODA) Based Approach for experimental data
133
5.6 Discussion and Conclusion for Inversion of GMR imaging .......................... 137

CHAPTER 6. POSITRON EMISSION TOMOGRAPHY IMAGING 140

6.1 Introduction .................................................................................................... 140
6.2 History of PET in Medicine ........................................................................... 141
6.3 Physics and instrumentation of PET .............................................................. 145
6.4 Challenges in PET imaging ............................................................................ 151

CHAPTER 7. FORWARD MODELS FOR POSITRON EMISSION

 

TOMOGRAPHY IMAGING 152
7.1 Introduction .................................................................................................... 152
7.2 Monte Carlo Methods in PET ........................................................................ 154
7.2.1 System Model ............................................................................................. 154
7.2.2 Random Number Generator ....................................................................... 155
7.2.3 Sampling Techniques ................................................................................. 157

vi

7.3 T wo dimensional (2D) model ........................................................................ 158

7.4 Three dimensional (3D) model ...................................................................... 163
7.5 System Data Acquisition ................................................................................ 164
7.5.1 Sinogram ..................................................................................................... 164
7.5.2 List-Mode Acquisition ................................................................................ 166
7.6 PET Image Performance and Quality ............................................................. 167
7.6.1 Scatter Noise ............................................................................................... 169
7.6.2 Random Noise ............................................................................................ 170
7.7 Results and Discussion ................................................................................... 170
7.7.1 2D results .................................................................................................... 170
7.7.2 3D Results .................................................................................................. 177
7.8 Motion Models ............................................................................................... 178
7.9 Motion Effects on Sinogram .......................................................................... 182
7.10 Motion Synthesis Results ............................................................................... 185

CHAPTER 8. INVERSE PROBLEMS IN POSITRON EMISSION

 

 

TOMOGRAPHY IMAGING - 189
8.1 Introduction .................................................................................................... 189
8.2 PET imaging with motion .............................................................................. 190
8.3 Current motion correction approaches ........................................................... 192
8.4 Motion Correction in Projection Space .......................................................... 194
8.5 Motion Tracking, Estimation and Extraction ................................................. 197
8.6 Motion Correction Implementation and Results ............................................ 202
8.7 Discussion and Conclusion ............................................................................ 206
8.8 Scaling-mapping method pseudo codes ......................................................... 208
CHAPTER 9. CONCLUSIONS AND FUTURE WORK- 210
9.1 Accomplishments and Conclusions ............................................................... 210
9.2 Future Work ................................................................................................... 210

vii

LIST OF TABLES

Table 4.1 Electromagnetic GMR imaging sensor parameters ..................................... 64

Table 4.2 Peak to peak values of Magnetic flux density (2 component) for various
sensor lifioffs ........................................................................................................ 77

Table 4.3 Peak values of magnetic flux density (2 component) Magnitude for various

top-layer conductivity .......................................................................................... 81
Table 4.4 Peak value of the signal magnitude verses bottom layer plate conductivity85
Table 4.5 Effect of fastemer conductivity on peak values of signal magnitude .......... 86
Table 4.6 Peak to peak signal magnitude vs. fastener to edge distance ...................... 93
Table 4.7 Peak to peak values of Magnetic flux density (2 component) for various

crack dimensions after substracting the crack-free signal .................................... 98
Table 5.1 Signal to Noise Ratio Comparison between conventional ECT and GMR

approaches .......................................................................................................... 139
Table 6.1 Comparison between PET and SPECT ......................................... 150

viii

LIST OF FIGURES

Figure 2.1 Noninvasive imaging techniques .................................................................. 5
Figure 2.2 A generic EM imaging system ..................................................................... 7
Figure 2.3 Schematic of SPECT/PET imaging system ............................................... 10
Figure 2.4 Relationship between the test object space and image space ..................... 12
Figure 3.1 Geometry of Crack Under Fastener (CUF) problem (Multilayer geometry
with subsurface defect) .......................................................................... 15
Figure 3.2 A boundary between two materials ............................................................ 18
Figure 3.3 Faraday rotation effect in MO sensor ......................................................... 22
Figure 3.4 A schematic of MO imaging instrument .................................................... 23
Figure 3.5 A simple GMR imaging sensor configuration ........................................... 24
Figure 3.6 Layout of the excitation current sheet ........................................................ 24
Figure 3.7 Response of GMR sensor to a field applied in the direction of “easy” axis
that is obtained experimentally in GMR system .................................................. 26
Figure 3.8 GMR sensor array with sheet current excitation foil ................................. 28

Figure 3.9 Typical GMR output images: defect-free (Boxl) and defective (Box2) and
central line plots: defect-free (dashed) and defective (solid) ............................... 29

Figure 4.1 Modeling geometry for the test sample in electromagnetic imaging for
NDE ............................................................................................. 37

Figure 4. 2 Photographs of (a) the sensor head at the zero-balancing calibration
position and (b) an overhead view of the calibration sample with slot down the
middle. The imaging sensors and excitation current are aligned parallel to the slot

.............................................................................................................................. 46
Figure 4.3 Sample Geometry used in the model .......................................................... 46
Figure 4.4 (a) Side view and (b) top view of the finite element mesh ........................ 47

Figure 4.5 Modeling results of the normal component of the magnetic flux density (a)
real part (b) imaginary part ................................................................................... 49

ix

 

Figure 4.6 Geometry of scan plan with excitation current parallel to the slot ............. 50
Figure 4.7 Comparison of demodulated and calibrated signal across the calibration

slot ........................................................................................................................ 51
Figure 4.8 Geometry of the test sample ....................................................................... 56
Figure 4.9 FE mesh for the test sample: (a) side view (b) top view ............................ 57
Figure 4.10 Magnetic flux density for the test sample geometry in Figure 4.9 at 400

Hz (a) Magnitude (b) Real component (c) Imaginary component ....................... 58
Figure 4.11 GMR data for crack free fastener in the test sample at 400 Hz (The line

scan is marked by the red dashed line) ................................................................. 59
Figure 4.12 Comparison between experimental and modeling signal with various

detection phases ........................................................................... 61
Figure 4.13 Frequency Parameter: simulation results at 100 Hz ................................. 66
Figure 4.14 Frequency Parameter: simulation results at 400 Hz ................................. 67
Figure 4.15 Frequency Parameter: simulation results at 500 Hz ................................. 68
Figure 4.16 Frequency Parameter: simulation results at 700 Hz ................................. 69
Figure 4.17 Frequency Parameter: simulation results at 2 kHz ................................... 70
Figure 4.18 Frequency Parameter: simulation results at 7 kHz ................................... 71

Figure 4.19 Liftoff Parameter: 0.0050 inch liftoff (a) real, (b) imaginary, (c)
magnitude ............................................................................................................. 72

Figure 4.20 Liftoff Parameter: 0.0095 inch liftoff (a) real, (b) imaginary, (c)
magnitude ............................................................................................................. 72

Figure 4.21 Liftoff Parameter: 0.0150 inch lifioff (a) real, (b) imaginary, (c)
magnitude ............................................................................................................. 73

Figure 4.22 Liftoff Parameter: 0.1 inch liftoff (a) real, (b) imaginary, (c) magnitude 73

Figure 4.23 Liftoff Parameter: 0.15 inch liftoff (a) real, (b) imaginary, (c) magnitude
.............................................................................................................................. 73

Figure 4.24 Liftoff Parameter: 0.2 inch liftoff (a) real, (b) imaginary, (c) magnitude 74

Figure 4.25 Line scans of real component (across the center of fastener) .................. 75
Figure 4.26 Line scans of imaginary component (across the center of fastener) ........ 75

Figure 4.27 Mixed line scans using optimum detection angle (across the center of
fastener) ................................................................................................................ 76

Figure 4.28 Peak values of Real, Imaginary, and Mixed MR signals vs. liftoff ......... 77

Figure 4.29 Magnetic flux density (2 component) for 28% IACS top layer
conductivity (a) real, (b) imaginary, (c) magnitude ............................................. 78

Figure 4.30 Magnetic flux density (2 component) for 29% IACS t0p layer
conductivity (a) real, (b) imaginary, (c) magnitude ............................................. 78

Figure 4.31 Magnetic flux density (2 component) for 30% IACS top layer
conductivity (a) real, (b) imaginary, (c) magnitude ............................................. 79

Figure 4.32 Magnetic flux density (2 component) for 31% IACS top layer
conductivity (a) real, (b) imaginary, (c) magnitude ............................................. 79

Figure 4.33 Magnetic flux density (2 component) for 32% IACS top layer
conductivity (a) real, (b) imaginary, (c) magnitude ............................................. 79

Figure 4.34 Magnetic flux density (2 component) for 33% IACS top layer
conductivity (a) real, (b) imaginary, (c) magnitude ............................................. 80

Figure 4.35 Peak value of the signal magnitude verses top layer plate conductivity.. 80

Figure 4.36 Magnetic flux density (2 component) for 30% IACS bottom layer

conductivity (a) real, (b) imaginary, (c) magnitude ............................................. 82
Figure 4.37 Magnetic flux density (2 component) for 31% IACS bottom layer
conductivity (a) real, (b) imaginary, (e) magnitude ............................................. 82
Figure 4.38 Magnetic flux density (2 component) for 32% IACS bottom layer
conductivity (a) real, (b) imaginary, (c) magnitude ............................................. 83
Figure 4.39 Magnetic flux density (2 component) for 33% IACS bottom layer
conductivity (a) real, (b) imaginary, (c) magnitude ............................................. 83
Figure 4.40 Magnetic flux density (2 component) for 34% IACS bottom layer
conductivity (a) real, (b) imaginary, (c) magnitude ............................................. 83
xi

 

Figure 4.41 Magnetic flux density (2 component) for 35% IACS bottom layer
conductivity (a) real, (b) imaginary, (c) magnitude ............................................. 84

Figure 4.42 Magnetic flux density (2 component) for 36% IACS bottom layer
conductivity (a) real, (b) imaginary, (c) magnitude ............................................. 84

Figure 4.43 Peak value of the signal magnitude vs. bottom layer plate conductivity .84

Figure 4.44 Real, Imaginary and Magnitude of magnetic flux density for conductivity
values of Ti fastener (1.0% IACS, 2.2% IACS, and 3.1% IACS) as indicated on
the left ................................................................................................................... 86

Figure 4.45 Real(a), Imaginary(b), Magnitude(c), Demodulated(d) magnetic flux
density (2 component) for fastener edge distance equals to 0.4 inch ................... 88

Figure 4.46 Real(a), Imaginary(b), Magnitude(c), Demodulated(d) magnetic flux
density (2 component) for fastener edge distance equals to 0.5 inch ................... 88

Figure 4.47 Real(a), Imaginary(b), Magnitude(c), Demodulated(d) magnetic flux
density (2 component) for fastener edge distance equals to 0.6 inch ................... 89

Figure 4.48 Real(a), Imaginary(b), Magnitude(c), Demodulated(d) magnetic flux
density (2 compenent) for fastener edge distance equals to 0.7 inch ................... 89

Figure 4.49 Real(a), Imaginary(b), Magnitude(c), Demodulated(d) magnetic flux
density (2 component) for fastener edge distance equals to 0.8 inch ................... 90

Figure 4.50 Line scans across the center of the fastener for real part of magnetic flux
density (2 component) .......................................................................................... 91

Figure 4.51 Line scans across the center of the fastener for imaginary part of magnetic

flux density (2 component) ................................................................................... 91
Figure 4.52 Line scans across the center of the fastener for mixed signal using ODA
.............................................................................................................................. 92
Figure 4.53 Edge effect on defect signal amplitude: Real, Imaginary, and Mixed
signals vs. fastener-to-edge distance .................................................................... 93
Figure 4.54 Magnetic flux density (2 component) for 0.2 inch crack (a) real, (b)
imaginary, and (c) magnitude ............................................................................... 94
Figure 4.55 Magnetic flux density (2 component) for 0.22 inch crack (a) real, (b)
imaginary, and (c) magnitude ............................................................................... 95

xii

 

Figure 4.56 Magnetic flux density (2 component) for 0.25 inch crack (a) real, (b)

imaginary, and (c) magnitude ............................................................................... 95
Figure 4.57 Magnetic flux density (2 component) for 0.3 inch crack (a) real, (b)

imaginary, and (c) magnitude ............................................................................... 95
Figure 4.58 Magnetic flux density (2 component) for no crack (a) real, (b) imaginary,

and (c) magnitude ................................................................................................. 96
Figure 4.59 Line scans across the center of the fastener for real image data .............. 96

Figure 4.60 Line scans across the center of the fastener for imaginary image data ....97
Figure 4.61 Mixed line scans using ODA of 70 degree .............................................. 97

Figure 4.62 Peak-to-peak values of Real, Imaginary, and Mixed MR signals vs. crack
area ....................................................................................................................... 99

Figure 5.1 The overall approach for automated rivet classification in magneto optic

Imaging ........................................................................................ 102
Figure 5.2 Feature Parameters definition in skewness function S1 ............................ 104
Figure 5.3 A set of boundary distances defined in skewness function Sz .................. 105

Figure 5.4 Parameters definition in skewness function S2 and corresponding results
............................................................................................................................ 107

Figure 5.5 Classification results in multi-dimensional skewness space: (a) l-D Space
(b) 2-D Space (c) 3-D Space .............................................................................. 110

Figure 5.6 Classification results in normalized central moment space: (a) M20-M40
Space (b) M20-M40-M22 Space ........................................................................ 113

Figure 5.7 POD curves from response data of automated MOI inspections of surface
layer cracks around fasteners ............................................................................. 115

Figure 5.8 ROC curves for Automated MOI on surface cracks with different flaw size
a .......................................................................................................................... 115

Figure 5.9 POD fits for Automated MOI inspections of surface cracks. (POD for SI

(dark, red), for S; (green, blue), both 2-parameter and 4-parameter fits were
performed) .......................................................................................................... 117

xiii

Figure 5.10 Three-step characterization for surface interior MO images (a) raw
images (b) enhanced images using motion-based filtering (0) detection and
classification of images ...................................................................................... 117

Figure 5.11 Characterization for structural surface edge MO images (a) bad edge rivet
image detection and classification (b) good edge rivet image detection and
classification ....................................................................................................... 1 18

Figure 5.12 Characterization for subsurface interior MO images obtained under 5 kHz
frequency and various intensities of magnetic field: (a)-(b) are raw and classified
images (high intensity), (c)-(d) are raw and classified images (low intensity) .. l 18

Figure 5.13 Surface plot of the image data showing the asymmetry in two lobes of

fastener image .................................................................................................... 121
Figure 5.14 Skewness vs. Frequency for results with 0.3” crack .............................. 121
Figure 5.15 Simulated line scans at different frequencies for 0.3” subsurface cracks

after ODA processing ......................................................................................... 123
Figure 5.16 2D Feature Space for simulated data at different frequencies ............... 125
Figure 5. l7 SNR vs. Frequency for simulated signal for test sample with 0.3 inch

subsurface crack ................................................................................................. 126
Figure 5.18 The test standard schematic .................................................................... 128

Figure 5.19 Raw GMR output images for the test standard at 400Hz: (fiom left to
right) Steel in-phase, Steel quadrature, Titanium in-phase and Titanium

quadrature components ...................................................................................... 129
Figure 5.20 Raw output images for the test standard at 400Hz: (a) defect-free, (b)
0.20”, (c) 0.22”, (d) 0.25” and (e) 0.30” subsurface cracks ............................... 129
Figure 5.21 MAG images at 400Hz: (a) defect-free, (b) 0.20”, (c) 0.22”, (d) 0.25” and
(e) 0.30” subsurface cracks ................................................................................ 129
Figure 5.22 Definitions of features to quantify the MAG signals: (a) F, and F 2, (b) F3,
(0) F4 and ((1) F5 ................................................................................................. 130
Figure 5.23 MAG signals for the test standard Ti rivets and 2-D classification scatter
plots .................................................................................................................... 131
Figure 5.24 GMR mixed images using ODA = 70 deg.: (a) defect-free, (b) 0.20”, (c)
0.22”, (d) 0.25” and (e) 0.30” subsurface cracks ............................................... 133
xiv

 

 

Figure 5.25 Central line plots of images in Figure 5.24 (a) raw image-inside rivet (b)
raw image-outside rivet, (c) mixed image-inside rivet (d) mixed image-outside

rivet ..................................................................................................................... 134
Figure 5.26 2D classification scatter plots for ODA approach .................................. 136
Figure 5.27 Illustration of SNRL definition ................................................................ 137
Figure 6.1 Schematic diagram of the PET imaging system with physics of
PET. . . .. ...................................................................................................................... 141
Figure 6.2 Three generations of PET scanner: the last is the PET/CT fusion

system ..................................................................................... 144
Figure 6.3 Illustration of radioactive decay in which the positron is generated ...... 146
Figure 6.4 Illustrations of the overall photon emission and detection process of PET

scanner .................................................................................... 148
Figure 6.5 Schematic of the PET detector structure ...................................... 149

Figure 7.1 Illustration of sampling the cumulative distribution function (cdf) for
simulating the annihilation event ....................................................................... 158

Figure 7.2 (a) 2D MC Model Geometry, (b) scatter-free (c) single scatter and ((1)
random noise ...................................................................................................... 161

Figure 7.3 Projection data acquisition and sorting in Sinogram: dashed lines and solid
lines represent the projections at an angle of -90 deg. and an arbitrary angle,
where the angle is between -90 deg. and 90 deg.. Each projection reprents one
complete view of the object at one particular angle by lines of response (LORs)

............................................................................................................................ 165
Figure 7.4 Illustration of list-mode acquisition ......................................................... 167
Figure 7.5 (a) True, (b) scatter and (0) random coincidences affecting the PET Image

Quality ................................................................................................................ 168
Figure 7.6 2D Monte Carlo Model for patient motion study with discretization grids

illustrated ............................................................................................................ 171

Figure 7.7 (a) 2D phantom object and detectors, (b) discretized solution domain with
detection tube(red dotted line) and (c) attenuation coefficients map ................. 172

XV

 

Figure 7.8 Simulation results with one million trials: (a) true image (b) simulated
sonogram and (c) reconstructed image ............................................................... 172

Figure 7.9 Reconstructed images at different scatter noise level from (a) to (d):
reconstructed images at alpha = 0, 0.2, 0.4 and 0.6 ........................................... 174

Figure 7.10 Reconstructed images at different random noise level from (a) to (d):
reconstructed images at alpha = 0, 0.2, 0.4 and 0.6 ........................................... 175

Figure 7.11 Motion-free PET images modeling and image reconstructions procedure
............................................................................................................................ 176

Figure 7. 12 3D view of the human torso model in the motion model, the chest wall,
lungs, heart and spine are simulated with 3D objects with varying shape at
various frequencies ............................................................................................. 178

Figure 7. 13 Motion signals for both pulmonary (circle & solid line) and cardiac (star

& solid line) motion ........................................................................................... 179
Figure 7. 14 Cross section of the 3D human torso modeled as cylindrical and elliptical

objects with different sizes ...... ’ ........................................................................... 180
Figure 7.15 Illustration of motion encoded sinograrn with motion detection and

estimation ........................................................................................................... 181
Figure 7.16 Mathematical parameters illustration for motion models ...................... 183

Figure 7. 17 Six respiratory phases with different lungs sizes for input phantom
images: exhalation minima (upper left) and inhalation maxima (lower right) .. 186

Figure 7.18 Sliced simulation results from 3D model: (a) Sinogram and (b)

corresponding reconstructed images using FBP ................................................ 186
Figure 8.1 Motion artifacts in PET in comparison to CT images .............................. 191
Figure 8.2 Overview of motion correction scheme ................................................... 195

Figure 8.3 Schematic of projection and image space motion correction in PET
imaging ............................................................................................................... 196

Figure 8.4 Fiducial points allocation for motion tracking ......................................... 197

Figure 8.5 Reconstructed image with fiducial points for motion tracking and
estimation ........................................................................................................... 198

xvi

Figure 8.6 Motion extraction in Sinogram based on prior information obtained using
fast Radon transform result in (a) and high intensity extraction in simulated

Sinogram (b) ....................................................................................................... 200

Figure 8.7 Motion effect estimation in Sinogram for right lung PET image in
noise-free case .................................................................................................... 201

Figure 8.8 Illustration of the line plots of motion-corrupted and motion-free Sinogram
intensity with boundaries estimation labeled ..................................................... 202

Figure 8.9 Sinogram before and after the motion correction. The cumulative motion
effects with widen ROI can be seen in (a) due to five cycles of respiration
simulated ............................................................................................................ 205

Figure 8.10 Reconstructed PET images: (a) target image at one particular respiratory
phase, (b) motion blurred image and (c) corrected image using the
scaling-mapping algorithm ................................................................................. 205

xvii

CHAPTER 1. INTRODUCTION

1.1 Motivation & Objectives

In the past decades, there has been a tremendous increase in interest in imaging, the
formation of images, and related physics, image processing and analysis for both
engineering and clinical applications. Due to the inter-disciplinary or
multi-disciplinary research fields that imaging covers, there are many unsolved
problems. Numerous imaging methodologies have been developed, viz, medical
imaging, electromagnetic imaging, molecular imaging, optical imaging, geophysical
imaging, etc. based on the technologies and underlying imaging physics. Imaging
techniques can also be simply divided into two categories: invasive and noninvasive
methods. Noninvasive imaging methods reveal the internal conditions of object of
interest without destroying or cutting/opening the test object. Such techniques are
widely used in both engineering and medicine where the objects are
structures/materials or in-vivo subjects, like human body. To understand the
noninvasive imaging better, both the forward models that simulate the image
formation and the inverse techniques that process and analyze the imaging data are
crucial. This dissertation comprises two parts, electromagnetic (EM) imaging for
metallic aircraft body inspection and positron emission tomography (PET) imaging

for in-vivo human subject imaging and diagnosis.

 

The first application is related to EM methods for nondestructive inspection of
airframe structures. Eddy current (EC) testing is the most widely used method for
aluminum structures. However it is limited by skin depth in detecting defects in
multi-layered complex structures. Electromagnetic noninvasive imaging techniques in
this dissertation, namely, magneto optic (MO) imaging and giant magneto resistive
(GMR) imaging offer much higher sensitivity that allow detection of defects located
deep in the test object. MO imaging is a technology, based in part on the Faraday
rotation effect that uses eddy current induction techniques along with a MO sensor,
which generates real-time analog images of the magnetic fields associated with
induced eddy currents interactions with structural anomalies. GMR imaging sensor is
unipolar in nature and is sensitive to the magnetic fields along the “easy” axis of the
sensor. The local magnetic fields measured by the GMR sensors are converted into
electric voltages, which are also associated with the induced current interactions with
structural anomalies. However, MO imaging measurement suffers severe background
noise and lack of quantification and automated analysis. GMR imaging data suffers
fi'om relative lower signal to noise ratio due to lower frequencies applied for
overcoming the skin depth effect; and lack of systems and sensors configuration
optimization due to little knowledge of actual electromagnetic field distribution and
sensor-field interaction. Positron Emission Tomography (PET) imaging is both a
clinical and a research imaging modality first introduced in 19708. However, patient
motion, including body movement, head movement and thoracic motion due to
respiratory and cardiac cycles, is a major source of artifacts in PET imaging. In the

motion-corrupted PET imaging data, various phases of motion information are

superposed and make the reconstructed images blurred. The motivation for motion
correction was not initiated, however, until recently when the spatial resolution of
PET scanners was greatly improved. The motion-corrupted blurring is relatively
larger and unbearable in the new generation of PET scanners. Also, the clinical
applications of PET have been emphasized more in thorax scans and tumor detection,
in contrast to the early research of the brain. The effects of the thoracic motion,
including cardiac and respiratory motion, are much more severe than the head
movement. Current compensation algorithms are performed after or during image
reconstruction process, e. g. multiple frame acquisition, gating methods, and image
registration. Also due to the trends in the combination of PET and Computed
Tomography (CT), the correction of anatomical motion in PET is essential for the
image fusion of PET and CT modalities and accurate attenuation corrections for PET
fiom CT acquisition.

The principle objective of this dissertation is to develop forward mathematical models
for noninvasive imaging and also develop inversion techniques and image processing
methods to interpret the measurement and reconstruct the true object accurately. More
specifically, in MO and GMR imaging application, a real-time reconstruction and
automated image classification approaches are presented. In the PET imaging
application, the specific aim is to develop a motion correction technique based on the
Sinogram data prior to the image reconstruction. In order to study the effects of
motion on PET image, we need to develop an accurate and reliable computational

model that allows controlled variation of motion due to respiratory and cardiac cycles.

1.2 Scope and Organization of the Dissertation

There are nine chapters in this dissertation and is composed two parts. Part I related to
noninvasive imaging of engineering structures comprises Chapter 3 to Chapter 5, and
is referred to as nondestructive evaluation (NDE). Chapter 1 and Chapter 2 are the
introduction and background of the dissertation, including the basic ideas of forward
and inverse problems in noninvasive imaging. Chapter 3 presents the fundamentals of
electromagnetic imaging including magneto optic and giant magneto resistive
imaging principles and challenges. Chapter 4 introduces the forward finite element
models for electromagnetic imaging techniques. System optimization and sensor
design are also investigated based on the numerical models. Chapter 5 discusses
image processing, characterization, classification and enhancement algorithms for
inversion techniques. Part 11 comprises Chapter 6 to Chapter 8 is devoted to
biomedical imaging, particularly PET. Chapters 6 to 8 discusses the biomedical
imaging application, namely, positron emission tomography imaging. Its history,
physical principles, instrumentation and challenges are discussed in Chapter 6.
Chapter 7 presents the forward 2D and 3D emission tomography models. The patient
motion is also considered and introduced in the models in this Chapter. Chapter 8
proposes the innovative motion correction methodology in the projection space, i.e.
Sinogram. The advantage and disadvantage of this method are discussed. Major
accomplishments and future work are summarized in the last chapter, Chapter 9,

followed by the Appendix and References.

CHAPTER 2. NONINVASIVE IMAGING

2.1 Introduction

Imaging, simply speaking, is the formation of images. Imaging science is concerned
with the formation, collection, duplication, analysis, modification, and visualization
of images.

A typical NDE/Noninvasive imaging system is illustrated in Figure 2.1, where the
detector/sensor makes measurement on the surface of the test objects, i.e. structures
and human body. Noninvasive imaging techniques are used in both engineering and
medicine to determine the state or internal conditions of the objects on the basis of

information contained in measured images or signals.

Detector/

- Sensor

Object
{Structures, Human body, ...}

 

 

 

 

 

Figure 2.1 Noninvasive imaging techniques

Noninvasive imaging is an interdisciplinary research field that combines
methodologies from Electrical Engineering, Biomedical Engineering, Medicine,

Radiology and Statistics. This dissertation covers two major research fields: 1)

electromagnetic imaging methods that deal with the materials and structures,
commonly referred to as NDE, and 2) biomedical imaging of human body or small
animals. Although the physical principles for these imaging techniques are different,
they can be unified by a common goal: both of them search for solutions of inverse
problems fi'om the observed images/signals. A brief description of the theory and
applications of various imaging modalities are in the following sections of this

Chapter.

2.2 Electromagnetic imaging

EM imaging is essential for detecting anomalies or defects in conducting materials
based on the electromagnetic principles. To understand EM imaging techniques, this
section is organized as following: the fundamental of electromagnetism is first
introduced; then the introduction of image processing is briefly discussed followed by
a general approach to processing and understanding the electromagnetic images. The
details of two electromagnetic imaging techniques, MO Imaging and GMR imaging,
are discussed in Chapter 3.

A generic electromagnetic (EM) imaging system is shown in Figure 2.2. The
excitation transducer couples the EM energy into the test objects. The receiving
sensors measure the response of energy/material interaction. This is the forward
imaging procedure. Different EM sensors are used for different applications, e. g. eddy

current, microwave imaging, Terahertz imaging, etc. After acquiring and storing the

EM images, those data are passed through the inversion techniques, which involve the

object reconstruction.

   

Inversion, Image
processing

\
I
\ I i /

\ EM Images J

Forward Imaging I
Using different sensors I

U
9.
m
>
S
C
E:
2:".
o

 

 

 

 

Figure-2.2 A generic EM imaging system

2.3 Biomedical imaging

Biomedical imaging is the imaging of human body (or parts of body) for clinical
purposes and is essential for medical diagnosis procedures to detect and treat diseases.
The different modes of generating biomedical images are referred to as modalities and
each modality is characterized by its own form of energy, data acquisition process,
physical principles of generating images, instruments and applications in medicine.

Commonly used biomedical imaging modalities include: conventional radiography

 

such as X-rays, computed tomography (CT), ultrasound, magnetic resonance imaging
(MRI) and nuclear medicine applications such as single photon emission computed
tomography (SPECT) and positron emission tomography (PET).

Thus biomedical imaging represents a series of noninvasive techniques that can
visualize and interpret the internal condition of human body in a clinical perspective.
Biomedical imaging is in some restricted mathematical sense seen as the solution of
inverse problems because the internal properties of human body is evaluated by from
the observed signals or images using different modalities. More introduction of

biomedical imaging can be found in Webb’s book [1].

2.3.1 Nuclear Medicine Imaging

Biomedical imaging application in this dissertation focuses on the PET imaging,
which is one of the nuclear medicine imaging modalities. In contrast to the
conventional planar X-ray imaging, CT, Ultrasound and MRI, nuclear medicine
imaging techniques focus on the distribution of radiopharrnaceuticals injected or
ingested into the body and physiologic or metabolic interactions between the
radioactive compounds and living tissues are imaged. These “radiopharmaceuticals”,
also termed radiotracers or radioisotopes, are compounds consisting of a chemical
substrate linked to a radioactive element. The way to introduce radioactive
compounds into human body include: inhalation into the lungs, direct injection into
the bloodstream, subcutaneous administration or oral administration [1]. The
abnormal rate of change at which the radiopharmaceutical accumulates and distributes

strongly indicates diseases in some particular tissues, which usually leads to a higher

contrast in the final images for diagnosis and distinguishes diseased from healthy
tissue. In nuclear medicine imaging, there are planar projection imaging techniques,
e.g. the anger scintillation camera and emission tomography imaging, like PET and
SPECT that depict the activity distribution in single cross sections of the subjects.

The nuclear medicine images are typically generated as follows: first, the radioactive
compounds undergo decay and the radiation, usually in the form of y rays is detected
using a gamma camera. Secondly, rather than using film as in X-ray imaging,
scintillation crystal is widely used to convert the energy of the y rays into light. Third,
the light photons are converted into electrical signals by photomultiplier tubes (PMT)
[1]. The signals are stored and analyzed by computers and the final nuclear medicine
images are produced and visualized. Positron emission tomography produces
emission images as opposed to transmission images, i.e. X-ray imaging, because the
energy of y rays emitted by the radioisotopes from inside the human body. A general

schematic of the nuclear medicine imaging system is shown in Figure 2.1 [7] and the

schematic for PET system is illustrated in Figure 2.3.

 

I Imaging Display

i

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

, Sinogram I
Computers/CPUS i
(Data binning, storage) List-mode Data I
I 7 Position Circuits I
I Electronic Circuits Energy Summing Circuits I
I PMTs
(Photomultiplier Tubes)
. Chystals
I Detector (Planar/Rings) (Na, BOO etc.) I
Collimator I

 

 
   

Patient/Object

Figure 2.3 Schematic of SPECT/PET imaging system

2.3.2 Positron Emission Tomography Imaging

Positron emission tomography imaging belongs to the class of nuclear medicine
imaging and provides functional information in patients by imaging the distribution of
positron emitting nuclides. In contrast to another nuclear medicine emission
tomography, SPECT, the PET system incorporates annihilation coincidence detection
(ACD) instead of collimation. The popular nuclides include 11C, 150, 18F, and 13N that

are incorporated into metabolically relevant compounds. Although there are many

10

radiopharmaceuticals suitable for PET imaging, for the applications discussed in this
dissertation such as cardiac imaging and lung tumor detection, 18F-FDG,
fluorine-18-fluorodeoxyglucose is more widely adopted. 18F-FDG is a glucose analog
that can differentiate not only malignant tumors from benign lesions, but also
differentiate severely hypo-perfused but viable myocardium from scar. More details

of PET imaging are discussed in Chapter 6 and Chapter 7.

2.4 Forward and inverse problems

The common goal of the noninvasive imaging applications is to interpret the internal
condition of the test object. This is referred to as the inverse problems in both
Engineering and Medicine. The measurement and image acquisition is referred to as
the forward problems. A schematic of the forward and inverse problems in imaging
techniques is depicted in Figure 2.4. There are two spaces involved in this loop, one is
the subject space and the other is the imaging space. For various imaging techniques,
the subject space includes all kinds of subjects with different properties to be imaged;
while the imaging space is a set of data, i.e. images, signals that are directly measured
or acquired by simulation. After understanding the concepts of subject space and
imaging space, it is' straightforward to define the forward problems and inverse
problems in imaging techniques. .

Forward problems involve the procedures and methods that the imaging space data
are obtained fiom test objects with various properties. The procedures can involve

direct measurement, forward system modeling and simulation based on underlying

ll

governing equations and so on. On the other hand, the inverse problems are the
procedures that the true properties of test objects in subject space are revealed from
the imaging space dataset. The concept of inversion techniques is actually very broad
in imaging, which includes image reconstruction, image classification, pattern
recognition, image visualization, etc. As we can see from Figure 2.4, the forward and
inverse procedures are crucial in imaging techniques because they are the bridges that

connect the test object space and the imaging space.

 

 

 

 

 

 

Forwud System Models
’ Underlylng Governlng Equations %
Subjects with 1 Measured or
various properties Sinulsted inages

 

 

 

 

 

 

 

 

lnverslon tectmlques, %
’ ImageISIgnal Reconstruction,

Vlsuallzatlon, Classification

 

 

 

Figure 2.4 Relationship between the test object space and image space

2.4.1 Forward problems in imaging

As we discussed earlier, forward problems in imaging are the procedures that acquire
images or signals from subjects in subject space, either by direct measurement or
modeling. In this dissertation, we focus on the forward system models that can be
represented as mathematical equations and explain the underlying physics. Simply

speaking, the imaging process is modeled by the forward models.

12

2.4.2 Inverse problems in imaging

The inverse problems and the forward problems inverse of one another, which means
the forward models input becomes the inverse output and vice versa. For historical
reason, the forward problems, not only in imaging but also in mathematics, physics or
mechanics, have been studied more extensively than the inverse problems. In this
chapter, we will spend more time on discussing the general ideas, issues and
techniques of inversion in imaging. More general introduction on inverse problems in
imaging, the readers can refer to the book by M. Bertero [4]. The inverse problems in
noninvasive imaging basically include reconstructing and/or interpreting an image of

conditions inside the test objects from noninvasive or nondestructive evaluation.

2.4.3 Mathematical description of inverse problems

The system of interest can be simply specified by input variables X and output
measurement Y. The measurement Y are functions of the input variables, which are
also called the system state variables, and their relationship is shown as equation 2.1.
Y = A(X)

(2.1)
The map A: X —) Y can be linear or nonlinear depending on the applications. In
this chapter, we will mainly discuss the linear problems and the inversion of equation
3.1 becomes simply the inversion of A. By assuming we are working on a finite
dimensional system and the map is linear and can be simplified into matrix form, we

have the following equation:

13

iii

l .‘r
.u

(2.2)
By assuming A is invertible, so the inverse problem is easily described as:
X=A4Y
(2.3)
Although equation 2.3 is straightforward mathematically, it is not always realistic due

to the ill-posed problem in inversion.

14

CHAPTER 3. ELECTROMAGNETIC IMAGING

3.1 Introduction

Detection and characterization of cracks under fastener holes (CUF) in aging aircraft
is a significant challenge in nondestructive evaluation (NDE) applications. Figure 3.1
shows the geometry of the multi-layered structure in the CUF problem. Eddy current
(EC) testing has been widely used for detecting such cracks for decades. However
ECT has its own limitations: difficulty in detecting deeper defects due to skin depth,
time consuming inspection process and cumbersome process of interpreting the
complex impedance plane signal. Advanced electromagnetic imaging methods, such
as eddy current Magneto Optic imaging and Giant Magneto Resistive imaging
technique have recently found widespread use in inspection of airframe structures.
Both these methods use eddy current excitation with the difference being in the

sensors used to detect the magnetic field associated with the induced eddy currents.

 

 

 

 

 

 

 

Figure 3.1 Geometry of Crack under Fastener (CUF) problem

In this chapter, we start fi'om the fundamental of electromagnetism, which helps us to
understand the physics of Electromagnetic (EM) imaging. Principles of the two EM

imaging, MO imaging and GMR imaging are then discussed. Challenges for these

'15

two EM imaging are then brought up with the proposed solutions using the forward

models and inversion techniques are presented in the following Chapters 4 and 5.

3.2 Fundamental of Electromagnetism

This section begins with a brief introduction of Electromagnetism, Electrical Field
and Magnetic Field. Then the Maxwell’s equations and related topics are discussed.
For more detailed introduction for fundamental of Electromagnetism, readers may
refer to [5] and [6].

The relationship between the magnetic field and electric field can be described by the

Maxwell’s equations, whose differential form are given by

VX5913
at
(3.1)
V xH = Js +J+§2
at
(3.2)
V-B=O
(3.3)
V-D=p
(3.4)

16

with the constitutive relations:

B = #H

(3.5)
D = (IE

(3.6)
J = of

(3.7)

In Equation (3.1) through (3. 7),

E = electric field intensity (volts/meter)

D = electric flux density (coulombs/square meter)
H = magnetic filed intensity (amperes/meter)

B = magnetic flux density (webers/square meter)

J 5.: source current density (amperes/square meter)

J = conduction current density (amperes/ square meter)
p = electric charge density (coulombs/cubic meter)

.11 = permeability (henrys/peter)

5 = permittivity (farads/meter)

‘7 = conductivity (Siemens/meter)

l7

III.)-
awe

ii] I»
«CI

28';

52,112 1 II

 

 

 

 

 

 

 

81,,ul I

Figure 3.2 A boundary between two materials

Maxwell’s equations should be satisfied everywhere in space, including the interface
joining different materials. Applying equations (3.1) — (3.4) to the material interface

shown in Fig 3.1 gives the boundary conditions: equations (3.8)-(3.11),

(E11 "E1)Xfi=0

(3.8)
(H11 ‘H1)Xfi =Js
(3.9)
(B,,-B,).a=0
(3.10)
(D11 “131)”? =Ps
(3.11)

where I”: is defined as the unit vector normal to the interface pointing from medium

I to medium 11. This implies that only the tangential component of the electric field

18

III:

‘I‘
bit

Ir.-
l 1

u...“

f )

II,

4.1

intensity E and the normal component of the magnetic flux density B are always
continuous on inter-material boundaries. The discontinuity of the tangential

component of H and the normal component of D are characterized by the surface

current J S and the surface charge p3.

Time-harmonic analysis plays an important role in various engineering applications.
A time-harmonic quantity refers to a variable which varies sinusoidally with time. A

time-harmonic quantity can be expressed in exponential form
_ 1mt
E (t) E pe

(3.12)

So, for time-harmonic fields, the Maxwell-Faraday’s law (3.1) and

Maxwell-Ampere’s law (3.2) become
V x E = —jcoB
(3.13)
V x H = J + jcoD
(3.14)

In this way, a time domain problem can be considered in the frequency domain.

Both MO and GMR imaging applications are quasi-static problems. As can be seen

from equations (3 .2) or (3.14), there are two kinds of currents: conduction current and

19

displacement current. The conduction current is proportional to the electric field

intensity, as stated by Ohm’s law,

(3.15)

The displacement current is defined as the time varying rate of electric flux density.

_aD_.

J———a)aE
datj

(3.16)

In many circumstances, the time varying rate is low enough such that 0' >> 608. In
such a case, the displacement current can be neglected and the quasi-static

approximation is said to apply to the Maxwell-Ampere’s law, i.e.
V x H = J
(3.17)
This equation implies that

V-J

ll
0

(3.18)

Consequently, in steady state, electric current density is divergence free, since there is
no accumulation of net charges. Two electromagnetic imaging applications: Magneto
optical imaging and Giant magneto resistive imaging are discussed in the following

sections.

20

3.3 Principle of Magneto optic imaging

Magneto optic (MO) imaging [27] [28] [29] is a relatively recent method used for
inspecting aging metallic airfi'ames for cracks and corrosion. MO imaging is a
technology, based in part on the Faraday rotation effect, uses eddy current induction
techniques along with a MO sensor, which generates real-time analog images of the
magnetic fields associated with induced eddy currents. The MO imaging is used in
aircraft inspection of rivet sites where data (images) from a large number of rivets
need to be analyzed.

MO imaging technology uses an induction foil for inducing eddy currents in the
conducting test sample and an M0 sensor for detecting the magnetic flux density
associated with the induced eddy currents. The MO sensor consists of a thin film of
bismuth-doped iron garnet grown on a substrate of gadolinium gallium garnet with its
easy axis perpendicular to the sensor surface.

The Faraday rotation effect observed by Faraday in 1845 states the plane of
polarization of a linearly polarized light transmitted through a MO material (MO

sensor) is rotated in the presence of a magnetic field as shown in Figure 3.3.

21

Sensor (garnet film)

Linear polarized light
before rotation

iill
IIII

Linear polarized light
after rotation

////
////

    
  
  

afl’r‘t‘i‘ifi" (Drum-4L 7 r a

: Applied magnetic filed
—.~_H_.____
Path Length

Figure 3.3 Faraday rotation effect in MO sensor

In a defect free specimen the induced current is uniform and the associated magnetic
field is in the plane of the sensor. Anomalies in the specimen, such as fasteners and
surface and subsurface cracks result in generating a normal magnetic flux density
along the easy axis of the sensor. If linearly polarized light is incident on the sensor,
the polarization plane of light is rotated by an angle 0 that is proportional to the

sensor thickness, given approximately by [28]:
9 z 9f (1? -M)1/(|IE||M|)
(3.19)
where k is the wave vector of the incident light, I is the sensor thickness, M is the

local state of magnetization of the sensor. Note that M is always directed parallel to

22

the ‘easy’ axis of the magnetization of the sensor which is perpendicular to the sensor
surface. When the reflected polarized light is viewed through the analyzer, the
presence of normal magnetic field is seen as a ‘dark’ area in the MO image. A
schematic of the MO imaging system using reflection mode geometry is illustrated in

Figure 3.4.

Light source
iii]

 

 

Polarizer
Bias c-oil /\ /
Sensor ® Induction foil

 

 

 

 

 

.---....................... . .........-...... .................... .-....
-!-.:.:.-.-.:.:.;.:.:.:.,.:.°.~ .-.- - .-.-.-.u .- - -. .- - - r. . . - - . .- 0 -. ~ - o -

g: ; .; . . .,.;.;.;.;.g.;.;.;.'.. ...‘.'._._..., ,

o'v'e'::n:o. .- l' I :-:I:o:.'-'o 0.. o 0...... I'q'l U s. . ' . .

 

Figure 3.4 A schematic of MO imaging instrument

An intrinsic property of MO materials is the domain structure [30]. The MO images
contain the images produced by these domain structures that must be accounted for in

automated image processing schemes.

3.4 Principle of Giant Magneto Resistive (GMR) imaging

Among sensors that are capable of measuring the field due to induced currents
directly, giant magneto-resistive (GMR) imaging sensors appear to be very promising

because GMR imaging sensors offer exceptional levels of sensitivity, small size and

23

low cost. A simple GMR imaging sensor configuration, in which a linear excitation
current sheet is applied, is shown in Figure 3.5, where the DC biased coil encircles
the GMR sensor. The layout of the current sheet for inducing uniform current in the
central region is shown in-Figure 3.6. The GMR imaging sensor with biasing coil is
located along the line of symmetry of current sheet, so that the output of GMR
imaging sensor without sample discontinuity is always zero. The lift-off of the GMR
sensor above the multi-layered structure with subsurface defects for the particular

problem is 0.0095 inches in this dissertation.

/

 

Figure 3.5 A simple GMR imaging sensor configuration

Line of symmetry '1 JVOltage feed points

 

 

 

 

 

 

Coil A

Figure 3.6 Layout of the excitation current sheet

24

The magnetic field generated by the excitation coil induces currents in the specimen.
The induced currents are distorted if a crack is encountered. The GMR imaging
sensor picks up perturbations in the fields associated with the induced eddy currents.
Signal processing algorithms may be required to demodulate the signal detected by

the GMR imaging sensor.

GNIR sensor is unipolar in nature and is sensitive to the magnetic fields along the
“easy” axis of the sensor. A typical response of GMR sensor to a field applied in the
direction of “easy” axis is obtained experimentally by the GMR system [54] and
shown in Figure 3.7. To quantitatively measure the magnetic flux density, the sensor
needs to be biased using the dc biasing coil (shown in Figure 3.5) to work around the
center of its linear range. The mathematical forms of the input and output signals of

the GMR system are given as:
Sr = 40 Sin(wt)
(3.20)
Sm = An sin(a)t + (on)
(3.21)
where S r is the excitation signal (also the reference signal) with the excitation
frequency a) and amplitude A0. For the sake of simplicity, the initial phase is

assumed to be zero.

25

GMR Calibration Curve

 

Output Voltage (V)
U) A M

N

 

 

 

 

-40 -30 -20 -10 0 10 20 30 40
Flux Density (G)

Figure 3.7 Response of GMR sensor to a field applied in the direction of “easy” axis that
is obtained experimentally in GMR system

Sm is the measurement GMR signal with the same frequency but different output
amplitude A,l and phase (on. The signal demodulation [59] first entails multiplying
the two signals together and results in a high frequency AC term and a DC component,
shown in Equation 3.22. Low pass filtering is applied to filter out the AC component

by integrating over a couple of cycles. More generally, an arbitrary initial phase 6

is introduced in this process as shown in Equation 3.23.

 

 

so) - 5.0) = 140141. sin<ax>sin<ax + a.) = - ’40:" cos<2ax + w.) + ”‘02": case.)

(3.22)

26

my my
0) a)

 

A041 j dtcos(2a1t+¢,, +5) j dt
- = — 0 n _ _9__
{5.0) Sm (01.; 2 ,% + 2 00800,. 5) ”a,

Idt I dt

0 0

= MOOS(¢n _ 6)
2
(3.23)

From Equation 3.23, the in-phase component S1 and quadrature component
S Q are easily obtained by setting the initial phase into 0° and 90°.

A0 All

 

 

nT
S 1 = {S,(t)« Sm (t)}o = AoAn I dt sin(a)t)sin(a)t + (pn) = cos(¢n)
0
(3.24)
A A .
SQ = {Sr(t)'Sm (t)}90 = 02 n Sln(¢n)
(3.25)

From the two components, both the amplitude An and phase q)" could be

extracted. The magnitude of the sensor signal can also be calculated as

(3.26)

27

 

Figure 3.8 GMR sensor array with sheet current excitation foil

Figure 3.8 shows an alternative approach to using the GMR imaging sensor where the
sinusoidal excitation current is applied to an induction foil. An array of GMR imaging
sensors can be used to pick up the normal component of magnetic flux density
generated by the induced currents. In the absence of any discontinuity, the magnetic
flux is tangential to the specimen surface resulting in a null signal. Anomalies in the
specimen result in a normal component of the magnetic flux density. The normal
component of the magnetic field is sensed by the GMR imaging sensor. Altemately, a
pulse excitation current can be applied to the induction foil and an array of GMR
imaging sensors can be used to measure the transient magnetic field at different
locations simultaneously. The output signals from the GMR sensor array contain
time and space information. The information can be used to estimate the defect depth
and location. Typical output magnitude images obtained using the GMR imaging
sensor are shown in Figure 3.9(a), where the Box 1 shows the defect-free image and
Box 2 shows the image of a fastener with a radial crack. The line plots across the
center of the images in Figure 3.9(b) clearly show the difference between the

defective fastener signal (solid line) and defect-free signal (dashed line) [54].

28

 

 

 

- - Without defect
— With defect

 

Field Magnitude a.u.

 

 

 

r r I l 1 l I l
SW 600 703 BED an "In 1100 1200

 

X-axis position index

Figure 3.9 Typical GMR output images: defect-free (Boxl) and defective (Box2) and
central line plots: defect-free (dashed) and defective (solid)

29

3.5 Issues in Electromagnetic Imaging

3.5.1 MO Imaging

The proposed research objective is electromagnetic imaging of this dissertation is the
development of a reliable, efficient, fully automated system for the inspection of
aging airplane skins. MO imaging is currently performed manually by scanning the
MO imager over the inspection surface. “Detection” using MOI depends on the
judgment of the inspector viewing the displayed MOI image. Even with training, this
is a subjective decision; therefore, this can produce varying results from inspector to
inspector depending on their ability to discern a defect. Another problem is the
inspector’s inevitable ennui as only a very small percentage of the images show any
anomalies. The MO images are 2D analog images of the 3D magnetic field
distribution. The image processing in MO imaging requires one to quantify the
asymmetry in the 2D images to determine the defect size and shape.

These problems can be mitigated with the application of modern image processing
and pattern recognition technologies to provide automated defect detection. An
automated real-time method for evaluating the MO images for structural defects can
reduce human error and increase accuracy and speed of the inspection. Such a method
will require the elimination of the background (image) noise due to magnetic domain

walls in the MO sensor and enhancing the image for subsequent interpretation.

30

3.5.2 GMR Imaging

GMR imaging technique offers much higher sensitivity enabling detection of defects
located deep in the test object. The interpretation and analysis of GMR data is crucial
to eliminate the operator variability and subjectivity and maximize the probability of
detection (POD). The GMR images are complex requiring higher computational
power for analyzing the images. This dissertation presents development of new
mathematical forward models, efficient image processing and crack detection
algorithms for GMR data. Through the forward and inverse techniques development,
the GMR imaging sensor can also be optimized to increase the GMR imaging system

detectability.

31

CHAPTER 4. FORWARD MODELS FOR

ELECTROMAGNETIC IMAGING

4.1 Introduction

Forward models using computational methods play an important role in
nondestructive evaluation (NDE). Numerical models serve many purposes such as
understanding the underlying physics, visualizing the fields, conducting a parametric

study and serving as a test bed to generate signals due to complex defect shapes.

Computational methods for modeling electromagnetic imaging problems fall into two
major classes, namely, analytical and numerical approaches. Analytical models are
capable of providing ideal, closed form solutions to problems. However, such
analytical models are, in general, limited in their application to simple specimen and
defect geometries, such as infinite half plane samples and regular, rectangular or
hemispherical crack geometry. Numerical models have gained popularity in recent
years due to their ability to simulate arbitrary shaped defects in complex sample
geometries. Numerical models can simulate problems with high dimensionality as

well as nonlinearity, anisotropy and inhomogeneous material properties.

With the development of faster and more powerful computers, the development and
implementation of the numerical models have become a major research focus in

nondestructive evaluation (NDE) [37] [3 8]. Several numerical modeling methods

32

have been developed for solving the governing equations describing various NDE
phenomena, such as finite difference method (F DM), finite element method (FEM),
boundary element methods (BEM), boundary integral methods (BIM) and volume
integral methods (VIM). Each method has its own advantages and limitations and is

appropriate for different kinds of numerical problems [39] [40] [41].

The finite difference method is a simple and intuitive numerical modeling method for
solving partial differential equations (PDE). The major advantage is its relative ease
for approximating partial derivatives by corresponding finite difference formulation.
It is widely applied for direct current, quasi-static, transient fields, and linear
problems. A limitation of the finite difference technique is its poor convergence
property, particularly in the case of irregular geometries. Also, finite difference
formulations do not lend themselves easily for modeling distributed parameters such

as the current densities, conductivities and permeability [42].

FEM evolved in the late 19503 as a numerical technique in structural analysis and was
quickly adapted as a major numerical modeling method in various engineering fields
to solve PDE. The finite element domain discretization scheme allows modeling of
complex shapes, making it efficient and relatively more accurate. Compared with
FDM, FEM has many advantages including its ease of imposing essential boundary
conditions and modeling complex geometries. FEM formulations can be easily
extended to handle higher order approximation and thereby lead to faster convergence
and better accuracy. The resulting matrix equations are in general sparse, banded and

diagonally dominant making the solutions process very stable [45]. FE models have

33

been developed and used for a variety of two and three-dimensional eddy current

NDE problems [43][46].

The drawbacks of FEM include large computer resources especially for nonlinear and
time-dependent problem. Also, the finite element method is not well suited for open
region problems [45] encountered in wave regimes. In the area of antenna and
electromagnetic wave propagation, integral equation based approaches are more
commonly used. BEM and VIM methods based on the integral equations [45] along
with absorbing boundary conditions are typically used in wave propagation problems

where the Green’s function is readily available.

Boundary integral methods solve Maxwell’s equations by integrating the appropriate
Green’s function over the boundary using the given boundary conditions. Hence there
is a need to discretize only the boundary surface in the integral equations, rather than
the entire solution domain. In the post-processing phase, the integral equation is used

again to calculate the required physical quantity at any point in the solution domain.

In the case of inhomogeneous problems, VIM is commonly applied to discretize the
volumetric source for solution. In VIM, the field is determined at a point by summing
the effects of the sources at all points. The volume integral computes a convolution of
the source points with the appropriate Green’s function. The advantage is it is only
necessary to construct a mesh over the test sample and solve for the currents in the
test sample. But this method also requires computation of the Green’s function which

in itself may not be trivial.

34

Integral equation based formulation is in general more efficient than differential
equation model with regard to computer resources but the resulting global stiffness
matrix equation is generally fully populated and ill conditioned leading to instabilities
[45]. More importantly VIM and BEM are inherently linear and can be used to solve
linear homogeneous problems. Lastly these methods are usefirl only when the Green’s

function for the problem is easily available.

In this dissertation, we use the finite element method for modeling the multi-layered
geometry largely due to the ease of modeling ferrite and non-ferrite fasteners and

complex defect shapes, as well as the inherent stability and accuracy of the solution.

4.2 Finite element methods

The finite element method (FEM) is one of the most widely used tools for solving
field problems. The method finds applications in diverse areas of engineering and
science largely due to its ability to model intricate geometries efficiently and
accurately. In the last few decades, FEM has been applied with great success in such
areas as the study of DC and low frequency electromagnetic fields in electrical
machines, design of structures, modeling thermal fields and nondestructive

evaluation.

A formulation of the problem from first principles derived from Maxwell’s equations
governing the differential equations to be solved in the model is provided by Jin [45].

Two- and three-dimensional electromagnetic finite element models for

35

 

;;1
\I

fly

.. 53‘

37‘
‘1

 

electromagnetic imaging applications have been developed [43][46][47][48][49].
More details of the MO imaging and GMR imaging modeling efforts can be found in
[55][62]. The major steps in finite element modeling are described briefly below.

Maxwell’s equations in differential form are given in Chapter 3, Equations 3.1 to 3.7.

The Vector potential formulation, A-V Formulation is used in this dissertation. From
Equation 3.3, the magnetic flux density is divergence fiee and can be expressed in

Equation 4.1 as the curl of a vector
B=VxA
(4.1)

where A is the magnetic vector potential. The magnetic field intensity can also be

expressed in terms of A, as

H = — V x A
,u
(4.2)
Substituting (4.1) into the Faraday’s Law (3.1), we get
V x E = —V x 6_A
at

(4.3)

V x [E + 241-) = 0

at

(4.4)

36

The electric field intensity can be express as
E = — j aIA — VV

(4.5)

The jw term is introduced since we are dealing with the time-harmonic problems.
Equation 4.5 tells us that the electric field intensity E can be expressed in terms of the

magnetic vector potential A and the scalar function V.
4.2.1 Sample Geometry and Mesh Generation

The detailed test geometry simulated in BM imaging is described in Section 3.3,
which has a multilayer geometry and consists of two layers of aluminum shown in
Figure 4.1. Only the general geometry is shown in this dissertation, without detailed

information

 

GIVIR imaging sensor I Scan direction it
——> I

I Finite multi—line foil excitation
!

 

 

 

   
 

 

 

Figure 4.1 Modeling geometry for test sample in BM imaging

37

The sample contains a top layer and second layer edges close to the fastener. Notches
are introduced at the fastener holes as shown in Figure 4.1. The three dimensional
domain is discretized using brick or hexahedral elements interconnected at 8 nodes.
The mesh size and complexity increases particularly when a small air gap of 0.005
inch is introduced between the layers and also between the fastener and plate. The
density of mesh elements were optimized and validated by comparing the model

prediction with experimental measurements.

4.2.2 Interpolation/Shape functions and Global Matrix Assembly

The magnetic vector A and electric scalar potential V in an element e is expressed in

terms of selected linear or quadratic interpolation fimctions as follows

8 24
A =Zleijx+Njijy+NjAzj =;N,Ak
J: =
(4.6)
8
e _ e e
V My.
j:
(4.7)

where

38

Njic k=3j—2
N;=<Nje.j} k=3j—1
LN-fé k=3j

 

(4.8)
and
A; k=3j—2
A,f= A; k=3j—1
A; k=3j
(4.9)

where ij, ij, AZ]. and VJ- are the three Cartesian components of the vector

potential and the scalar potential at node j; and N J- is the shape function associated

with node j that assumes the value of one at this node and the value of zero at any

other node. The vectors 55, )7, 2 are the Cartesian unit vectors.

Using the Galerkin formulation [45], and imposing the Coulomb gauge on Equation

4.1 and 4.6, we obtain

24 1 1
—VxN.e -VxNe. —V-N." V-N‘f (IV A"?
21%;: .)( ,)+(fl ,)( ,) } ,+

j:

39

24 8
Eff. W1" f WWI/M: +2149 0N: -VN§dV}Vf -
e j= e

j:

I Nf-(—1—Verxii)dS—-I Nf-(filV.Ae)dS=I Nf-deV
a0. .11 69. .U Q.

(4.10)

,for i=1,2,--°,24.

24 8
~2ti J’wVNi -N§dV}A;T + 21% WM: -VN;’dV}Vje +
1: e j= e

Ian 0N;(—jar4 —VV) ~1’idS = 0

(4.11)

for k = l,2,...,8 where it is the outward unit vector and 60,, is the boundary of
the sub-domain Qe . Combining these equations together, we obtain a

32 x 32 matrix equation,

[G16 [U ]. =IQIe
(4.12)

The global stiffness matrix is obtained by assembling each element matrix [Gr

together. The value at each entry is the sum of values contributed by all the connected

elements.

40

GU Q

(4.13)

where G is a complex, symmetric sparse matrix, U is the vector of unknowns
consisting of the electric scalar potential and the three components of the magnetic
vector potential at each node and Q is the load vector incorporating the current

source.
4.2.3 Boundary and Interface Conditions

In order to obtain a unique solution, appropriate boundary conditions need to be
imposed before solving the system of equations; Either Dirichlet boundary conditions
(values of A or V on the boundary) or Neumann boundary conditions (values of the
magnetic flux density B or the magnetic field intensity H on the boundary) need to be
specified. Since Neumann boundary conditions are usually included implicitly in the
finite element formulation, we only need to impose Dirichlet boundary conditions.
Current continuity conditions will be explicitly imposed at the interface boundaries to

avoid spurious solutions.
4.2.4 Matrix Solution

The finite element procedure results in a linear algebraic system of equations that
must be solved to determine the unknown coefficients of the shape fimctions. The
equations can be solved using either direct or iterative methods. The Gaussian

elimination method is a direct approach that can be used to solve either full or sparse

41

matrix equations. In the finite element method, the stiffness matrix is usually sparse
and banded. In such cases, it is more advantageous to employ iterative methods
since they are more efficient with respect to computational cost and data storage. A
number of iterative methods have been proposed over the years. However, only a
small number of them can be used for solving the complex-valued matrix equations

arising in quasi-static electromagnetic problems.

An important issue in iterative methods is its convergence properties. The
convergence of iterative process depends on the properties of the matrix, such as its
eigen values, the singular value or condition number. Usually a large condition
number will result in a large solution error or sometimes even failure to converge.
Thus the process of “preconditioning” is essential for the success of iterative methods.

The overall solution procedure requires 0(N) multiplications and additions [50].

4.2.5 Post-Processing

After we solve for the potentials, physical and measurable quantities of interest, such
as the magnetic flux density B, coil impedance and induced current density J, can be

calculated using:
B=VxA
(4.14)
E = — j arA — V V

(4.15)

42

 

’51?-

L51.“

J=0E

(4.16)

Where Equations 4.14 through 4.16 are the same as Equations 4.1, 4.5 and 3.7 and

presented here again for the reader’s convenience. GMR imaging sensors measure the

normal component of the magnetic flux density B. The different components of the

magnetic flux density can be obtained from the nodal values of magnetic vector

potential as:

 

6A 6A

 

 

 

 

(4.17)

(4.18)

(4.19)

The current density J inside element e is determined at the element center fi'om the

relation in Equation 4.20 and the components J; , J: can be calculated similarly.

43

 

 

6W" 8 8 EN?
J"=—'a)0'Ae—0' =—'a)0' NiAe.—0' ’Ve
x J x 6x J g]: l x' ,2]: 6x '

(4.20)

4.2.6 Modeling MO and GMR sensors with Induction Foil Excitation

The excitation source in this project is applied to an infinite planar induction foil
located parallel to the surface of the specimen. For a linear sinusoidal current

excitation with angular frequency of a) , the source current is
. '01
1,," = Re[JOe’ ]= JO coswt

(4.21)

where J 0 is the magnitude of the current. However a probe with a linear excitation

current flow pattern is sensitive only to cracks that are perpendicular to the current

direction. In this project, to best detect a crack oriented in the x direction the current

flow is assumed to be along the y direction, i.e., J in = 5410 where j} is the unit
vector along the y-axis. The normal component of magnetic flux density 32 is

measured by the GMR imaging sensor array.

4.3 Forward model results

Although difference imaging sensors are used in the MO and GMR imaging, the EM
fields generated by the excitation current are the same. Hence the same forward
model results can be used for simulating both MO and GMR signals. We start from a
two dimensional validation using a simple geometry and finally analysis the EM field

of the real testing sample.

4.3.1 Validation Results

The finite element model developed for electromagnetic imaging inspection of
multilayer test samples was first validated using experimental GMR measurements

from the calibration notch sample.

The validation sample is made of standard aluminum, which is a plate of dimensions
16 inches times 10 inches and of thickness 0.625 inch. A notch of dimensions 0.05
inch width and 0.25 inch depth is machined at the center of the sample. A photograph

of the sample is shown in Figure 4.2.

45

 

 

 

 

 

(a) (b)

Figure 4. 2 Photographs of (a) the sensor head at the zero-balancing calibration position
and (b) an overhead view of the calibration sample with slot down the middle. The
imaging sensors and excitation current are aligned parallel to the slot

The three dimensional domain of the finite element model and the corresponding

mesh are shown in Figures 4.3 and 4.4 respectively.

I U—O’ZIAluminumI B625"
11

 

 

0.05"

10"

Aluminum

 

 

 

 

 

1 6" VI

Figure 4.3 Sample Geometry used in the model

46

y-z plane at x = -127 mm
0.01

0.005

 

-0.005
-0.01
-0.015

 

-O. 02 -

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

-002 -001 0 0.01 0.02

(a)

x-y plane atz = -5.08 mm

 

-0.02 -0.01 0 0.01 0.02

(b)

Figure 4.4 (a) Side view and (b) top view of the finite element mesh

47

The FEM model solutions produce the electric scalar potential V and magnetic vector
potential A at each node in the solution region. The magnetic flux density B = V x A
is computed and its normal component is compared with the measured GMR signal.
Figure 4.5 shows the normal component of the magnetic flux density predicted by the
model. The 2D scan across the slotted sample is applied to obtain the sensor signal.

Figure 4.6 shows geometry of the scan.

48

I
-.O
0.00
Ls—xtn

 

Magnetic Flux Density (Gauss)-Real
O

 

 

 

U)

(U

E

1 0.4

3

8 0.4

g 0.2

3 0.2 I“ ..

8 0 0
-0.2 ,

’3

E -0.4 5 _0-2

.2 0.1

E

5:

m -o.4

E

 

(b)

Figure 4.5 Modeling results of the normal component of the magnetic flux density (a)
real part (b) imaginary part

49

 

 

 

Line Scam

o—b_o—o o—oi— —.— —.->

‘i’ .
Excrtatron

Current

 

 

 

 

 

 

 

 

 

 

I lurr inumJ

 

Figure 4.6 Geometry of scan plan with excitation current parallel to the slot

Experimental signals for the sample were available for various demodulation angles
as described in Section 3.4, principle of GMR imaging. The in-phase and quadrature
components were demodulated to generate a mixed signal using varying
demodulation angle 5. The mathematical form for the calibrated demodulated signal

for sensor 11 at arbitrary detection angle 6 is:

(SO-Sn)§=an((SO-SM) —--S<SO noo))cos(5)
+a ((SO-S ) fl —<SOS ,,°>/2 )sin(§)

(4.22)

where (S0 ~Sn)5 = —A£24’lcos(¢n — 6)

50

In finite element modeling, the zero balancing term <50 ~Sno>6 and the scaling

factor an are 0 and 1 respectively since these two terms are only useful in
experimental data analysis. From the validation results show in Figure 4.7, we can
conclude that the model result accurately predicts the experimental measurements.

Figure 4.7 Comparison of demodulated and calibrated signal across the calibration slot

m ..n ~-—~ and]. mm Inna—c

      

 

-l
I‘

r. .- .m.
m an...

‘ ' ' 4/ M...
u x as a a‘. I
..............

m ‘9“)

 

 

1...--. M... a... a...“

a
I
g

 

1. .. i 1'.
_- .. . no...“

 

51

Figure 4.7 continued
x10’5

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

-2 -I.5 -l.0 -0.5 0 0.5 1.0 1.5 2.0

Distance (inch)
(‘0)

 

 

x10‘5
6 r '
4 1.. ............ I .......................... 4 ................. I; ........ .—

 

 

8:15 ., ........... ..

 

 

 

 

 

 

 

 

 

 

oooooooooooooooooooooooo

 

 

 

 

 

 

 

 

-2 -1.5 -1.0 -0.5 0 0.5 1.0 1.52.0
Distance (inch)
(C)

52

Figure 4.7 continued
x 10-5
6 I

 

I I I I I I
E 5
I D O O : z 0
: : : : : : :
: : : : : : :
4 b IIIIIIIIIIII §0.-I.I....II‘3'...'OIOCCOOCOEDOOOOOCOOIOCBE ............. 3 IIIIIIIIIIIII iIOOOICOOOOOCI; .............

5 E E E i 3 i
: : : : 5 E i
s s s s . . I

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

.2 -1.5 -1.0 -0.5 0 0.5 1.0 1.5 2.0
Distance (inch)
(4)

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

-2 -1.5 .10 0.5 '0 0.5 1:0 1.5 2.0
Distance (inch)
(e)

53

Figure 4.7 continued
x l 0‘5

6

4
2

0

-2

 

-2 —l.5 -l.0 -0.5 0 0.5 1.0 1.5 2.0
Distance (inch)
(0
1110'5

 

-2 -1.5 -1.0 -0.5 0 0.5 1.0 1.5 2.0
Distance (inch)
(g)

54

Figure 4.7 continued
x10'5

6

4

 

-2 ~15 -1.0 -0.5 0 0.5 1.0 1.5 2.0
Distance (inch)
01)

4.3.2 3D Results

A schematic of the three dimensional (3D) geometry of the sample with subsurface
crack, imaging sensor configuration and scan direction are shown in Figure 4.8,

which is the same as Figure 4.1.

55

 

. Scan direction
GMR gin . I I
Finite multi—line foil excitation "“a g sens” .._.... .

_-.-._,_._,-.-..‘

    

 

—— Air boundary

 

Figure 4.8 Geometry of the test sample

The fastener material in this sample is Titanium. The side view shows top layer
truncated in accordance with the sample geometry. The model also simulates the gap
between the two top skins of the test sample. The current is parallel to the notch
length and the scan direction is perpendicular to it. The optimized 3D finite element

mesh for this sample with edge is shown in Figure 4.9.

56

y-z plane at x = -115 mm

 

-0.03 -0.02 -0.01 0 0.01 0.02 0.03 0.04

(a)

x-y plane atz = -O.79714 mm
0.04

O .03
0 .02

0.01

-0.01
-0.02
-0.03

-0.04

 

-0.05 -0.04 -0.03 -0.02 -0.01 O 0.01 0.02 0.03 0.04 0.05

(b)

Figure 4.9 FE mesh for the test sample: (a) side view (b) top view

57

The modeling results at 400Hz excitation frequency are shown in Figure 4.10. The
real and imaginary of the normal component of magnetic flux density are shown in
Figures 4.10 (b) and (c). The magnitude of the magnetic flux density is presented in

Figure 4.10

Figure 4.10 Magnetic flux density for the test sample geometry in Figure 4.9 at 400 Hz (a)
Magnitude (b) Real component (c) Imaginary component

 

(b)

58

Figure 4.10 continued

 

(C)

The corresponding experimental data from the rivet is shown in Figure 4.11.

 

Figure 4.11 GMR data for crack free fastener in the test sample at 400 Hz (The line scan
is marked by the red dashed line)

A comparison between simulation and experimental results is presented. Line scans in
the x-direction across the center of the fastener is extracted from the real and
imaginary image data (both the simulation and experiment) which are then mixed
using different demodulation or detection angles. The detection angle a and the mixed

GMR sensors signal is given by

Sa 2 SO cos(a) + S77 sin(a)
2

(4.23)

where S a , S0 and 53,2 are the mixed, in-phase and quadrature signals

respectively.

The mixed model and experiment signals are compared for a number of detection
angle values in Figures 4.12 with the detection angle varying from 80 to 190 degrees.
Both the magnitude and shape of the signals are remarkably similar. The results in
Figure 4.12 indicate the use of the model in conducting a parametric study as

described in the following sections.

60

Figure 4.12 Comparison between experimental and modeling signal with varying
detection angles

 

 

 

Detection I Experiment Model I
Angle

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

61

Figure 4.12 continued

 

 

Detection Experiment Model
Angle

 

”a“ on!!!)

 

 

 

 

 

 

130

 

 

140

 

 

 

x

“n...”

 

 

 

 

 

u u t; i m ms at a: o n ; -u u u I; I O. u u a: I
am

 

 

 

62

Figure 4.12 continued

 

 

 

 

 

 

 

   
 

 

       

 

 

 

 

 

 

 

  

 

 

     

 

 

 

 

O
Detection Experiment Model
O)-.. OM- “mi-Mutt! ‘1
I . . .7 ....... ,. T
3.... .m i - . .. .. . -
a an ., .. . - . . ' . .4. .
am. , I' . t: 4
a. ‘III- t v - b 9 v t u . ‘
m... .. -..-. . . ..‘. . ..
. . a... . . - . . .. . . . ..
.....i..i..i ...i
L A A A 1 x A J I 1 L A
u -u i: l a M n u I l.‘ J ‘5 M N ‘3 ' ‘1 H N ‘ 0
luv-an!
tram-”ma ”warm“:- w
y» l .4 -. v ' ’
on "I” . . . - .. i . <
or. ..... , . - . .. _ i. .
I, in e ' v a c . A .n 4
; .
170 ... _, . . . . ., .
w “n . a . .. n . . u
. .4
.n , . .
I. .l‘ '1 3 a. (I. ‘l {I D .7 0! 0
ID» ...L ..... y z. a...
is» .a v , - . ‘- . r
en ~ . I . a - .
nm- -» .. . .... v .. . t .
.l‘r"" -¢- ~ ~»-~--- .— ‘- - - ‘
‘flr - v v . - a _ . _ .
...- . A .. n i
.t u I. u r: ' u a. u n I
.~“
Mun-«emuw-m-o
nu ‘*"r=‘f'T“‘fi".—ir~"T“;—’?“"?"'f"

 

 

 

 

   

 

 

 

63

4.4 Electromagnetic Imaging Parameters

The validated FEM model was used to study the effect of various experimental
parameters on the GMR signal. Table 4.1 lists parameters along with the rage of

values.

Table 4.1 Electromagnetic GMR imaging sensor parameters

 

 

 

 

PARAMETER RANGE OF VALUES
Frequency (kHz) 0.1, 0.2, 0.4, 0.5, 0.7, 1.0, 2.0, 3.0, 5.0, 7.0
Sensor liftoff (inch) 0.0050, 0.0095, 0.015, 0.1, 0.15, 0.2
Conductivity of top layer (%IACS) Between 28 and 33

 

Conductivity of bottom layer (%IACS) Between 30 and 36

 

 

 

 

Conductivity of fastener (%IACS) Titanium (6A1, 4V): 1.0, 2.2, 3.1
Fastener to edge distance (inchO 0.40, 0.50, 0.60, 0.70, 0.80
Crack dimensions (inch) 0.20, 0.22, 0.25, 0.30

 

 

 

4.4.1 Frequency

Since most aircraft geometry contains multiple layers with cracks in second and third
layer the operating frequencies is largely at the low end. The range of frequencies
considered in this study is: 0.1, 0.2, 0.4, 0.5, 0.7, 1.0, 2.0, 3.0, 5.0 and 7.0 kHz. The
geometry of the test sample and sensor configuration and finite element mesh are

shown in Figures 4.8 and 4.9. A 0.3 inch second layer defect is introduced under the

64

fastener through the second layer on the side of the stringer edge. The normal

component of magnetic flux density Bz is calculated for each frequency, for both, the
crack free case and 0.3 inch crack case.

Figure 4.13 shows the simulation results at a frequency of 0.1 kHz. The images in
Figure 4. 13(a), 4.13 (c) and 4.13 (e) are the magnitude, real and imaginary parts of the
normal component of flux density for crack free case at a lift-off of 0.0095 inch.
Images in Figure 4.13 (b), 4.13 (d) and 4.13 (i) are the corresponding magnitude, real
and imaginary parts of the normal component of flux density for 0.3 inch crack
geometry. Figures 4.14 to 4.18 present similar model predicted results for frequencies
of 400, 500, 700, 2000 and 7000 Hz. The real, imaginary and absolute magnitude of
the normal component magnetic flux density for both defect free and 0.3 inch crack

are presented.

65

 

 

Figure 4.13 Frequency Parameter: simulation results at 100 Hz

66

 

 

 

(e) (0

Figure 4.14 Frequency Parameter: simulation results at 400 Hz

67

 

(e) (0

Figure 4.15 Frequency Parameter: simulation results at 500 Hz

68

 

(e) (0

Figure 4.16 Frequency Parameter: simulation results at 700 Hz

69

 

(a) (b)

   

 

(C) (d)

     

(e) ( 0

Figure 4.17 Frequency Parameter: simulation results at 2 kHz

70

 

(e) (t)

Figure 4.18 Frequency Parameter: simulation results at 7 kHz

71

4.4.2 Sensor Liftoff

The parametric effect of sensor lift-off on the signals was studied. Liftoff values
varying fi'om 0.0050, 0.0095, 0.0150, 0.1, 0.15, and 0.2 inches were considered. The
geometry and finite element mesh of the test sample and sensor configuration are

shown in Figures 4.8 and 4.9. A 0.3 inch second layer defect is introduced under the

fastener.

 

(a) (b) (C)
Figure 4.19 Liftoff Parameter: 0.0050 inch liftoff (a) real, (b) imaginary, (c) magnitude

 

(a) (b) (C)
Figure 4.20 Lifioff Parameter: 0.0095 inch lifioff (a) real, (b) imaginary, (c) magnitude

72

 

(a) (b) (C)
Figure 4.21 Lifioff Parameter: 0.0150 inch lifioff (a) real, (b) imaginary, (c) magnitude

 

(a) (b) (c)
Figure 4.22 Liftoff Parameter: 0.1 inch liftoff (a) real, (b) imaginary, (c) magnitude

 

(a) ’ (b) (c)

Figure 4.23 Liftofl‘ Parameter: 0.15 inch liftoff (a) real, (b) imaginary, (c) magnitude

73

 

(a) (b) (C)
Figure 4.24 Liftoff Parameter: 0.2 inch liftoff (a) real, (b) imaginary, (c) magnitude

The normal component of magnetic flux density 32 is calculated at excitation

frequency of 400 Hz. The real, imaginary and magnitude of the normal component of
magnetic flux density images are plotted for each individual liftoff value in Figures

4.19 through 4.24.

In order to see the variation in the signals with sensor liftoff, the peak to peak values
in line scans across the center of the fastener are considered. The real and imaginary
parts of the line scans are plotted in Figures 4.25 (Real) and 4.26 (Imaginary). The
real and imaginary data are demodulated using the optimum detection angle of 70

degrees at 400 Hz and the resulting mixed signals are shown in Figure 4.27. The peak

to peak values of 32 (real, imaginary, magnitude and mixed) at each lifioff are listed

in Table 4.2 and plotted in Figure 4.28.

74

x10'8 Sensor Liftoff Line Plots(Real)
2

 

 

 

 

 

-— '0. 005'
— 0.0095'
0.015"
-—-- 0.053'
— 0.1'
— 0.15'
0.2"

 

 

 

-2
-3

o.—

-2 -1

4
x104

Figure 4.25 Line scans of real component (across the center of fastener)

x10'9 Sensor Liftoff Line Plots(Imag)
8

 

 

A

  

 

—— 0.005'

— 0.0095'
0.015'

---- 0.053'

—"—0.1'

 

 

 

Figure 4.26 Line scans of imaginary component (across the center of fastener)

1:5 2
x10"4

 

x "3'9 Sensor Liftoff Mixed Line Plot using ODA

 

3 I l l l I l l

   
 

 

 

 

 

 

 

 

2 - A _ \
\
1 _
0 _
-1 -
— 0.005"
— 0.0095"
'2 ‘ 0.015"
—— 0.053"
-3 _ —- 0.1"
— 0.15" ’\
— 0.2"
_4 I l I 1 l l
-2 -1.5 -1 as o 0.5 1 1.5

Figure 4.27 Mixed line scans using optimum detection angle (across the center of
fastener)

76

Table 4.2 Peak to peak values of Magnetic flux density (2 component) for various sensor

 

 

 

 

 

 

 

 

 

 

 

 

 

 

lifioffs
SENSOR REAL IMAGINARY MIXED
LIFTOFF COMPONENT COMPONENT AFTER ODA
(INCH) (TESLA) (TESLA) (TESLA)
0.005 3.3484e-008 1.3586e-008 6.31 19e-009
0.0095 3.1232e-008 1.2620e-008 5.9523e-009
0.015 2.8929e-008 1.1657e-008 5.5312e-009
0.053 1.9075e-008 7.66676-009 3.2867e-009
0.1 1.2221e-008 5.3627e-009 1.5625e-009
0.15 7.7481e-009 4.5275e-009 7.3745e-010
0.2 4.7604e-009 4.2386e-009 1.7312e-009
x10'8 Peak to peak values vs. sensor liftoff
-0- Real Component
->- lmag Component
-I- Mixed at ODA=
‘U
m
m
r
6'
'0
m
m
2'
i
C
m
C.
m
9.
3’,

 

00 0.02 0.04 0.06 0.08
Sensor Liftoff (inch)

0.1 0.12 0.14 0.16 0.18 0.2
Figure 4.28 Peak values of Real, Imaginary, and Mixed MR signals vs. lifioff

77

4.4.3 Conductivity of Top Layer

The conductivity values for the top layer were selected according to the Parametric
Table presented earlier. The operation frequency is chosen as 400 Hz and liftoff is
chosen as 0.0050 inch. The nominal conductivity for top layer plate for test sample is
29.6% IACS, and the range of conductivity of the top layer plate is assigned values
from 28% to 33% IACS. Figures 4.29 through 4.34 present simulation results of real,
imaginary and magnitude values of the normal component of magnetic flux density
for different values of top layer conductivity. Table 4.3 shows the peak value of the

signal magnitude at different values of conductivity.

1
1'
(a) (b) (C)

Figure 4.29 Magnetic flux density (2 component) for 28% IACS top layer conductivity (a)
real, (b) imaginary, (c) magnitude

‘ ..
l
(b)

Figure 4.30 Magnetic flux density (2 component) for 29% IACS top layer conductivity (a)
real, (b) imaginary, (c) magnitude

     

(a)

78

 

(a) (b) (C)

Figure 4.31 Magnetic flux density (z component) for 30% IACS top layer conductivity (a)
real, (b) imaginary, (c) magnitude

   

(a) (b) (C)

Figure 4.32 Magnetic flux density (2 component) for 31% IACS top layer conductivity (a)
real, (b) imaginary, (c) magnitude

       

(a) (C)

Figure 4.33 Magnetic flux density (2 component) for 32% IACS top layer conductivity (a)
real, (b) imaginary, (c) magnitude

79

     

(a) (C)

Figure 4.34 Magnetic flux density (2 component) for 33% IACS top layer conductivity (3)
real, (b) imaginary, (c) magnitude

Peak Value of Magnitude vs. Top layer

-8
281 0 Conductivity

 

2.09 - -

2.08 - ~

2.07 ~ .

2.06 - .

2.05 - _

2.04

2.03

(e|se_L) epmgufiew need or need

2.02

 

2.01 - -

 

 

2 i I I I
28 29 30 31 32 33

Top Layer Plate Conductivity (%IACS)

 

Figure 4.35 Peak value of the signal magnitude verses top layer plate conductivity

80

Table 4.3 Peak values of magnetic flux density (2 component) Magnitude for various
top-layer conductivity

 

 

 

 

 

 

 

CONDUCTIVITY CONDUCTIVITY MAGNITUDE
(% IACS) (SIEMENS/M) (TESLAL
28 1.624067 2.014le-8
29 1.6820e7 2.0150e-8
30 1.739167 2.0206e-8
31 1.798067 2.0312e-8
32 1.8560e7 2.0358e—8
33 1.9140e7 2.0419e-8

 

 

 

 

 

Figure 4.35 shows a plot of peak value of magnitude signal magnitude verses top

layer plate conductivity.

From the results in Table 4.3 and the Figure 4.35, it is seen that the peak value of the
signal increases slightly when top layer plate conductivity increases. This is as
expected since increase in the top layer plate conductivity increases the induced eddy

currents of that layer thereby increasing the amplitude of the measured signal.

4.4.4 Conductivity of Bottom Layer

The conductivity values for the bottom layer were selected according to the
Parametric Table presented earlier. For bottom layer plate, the conductivity value
ranges from 30% IACS to 36% IACS. The operation frequency is chosen as 400 Hz

and liftoff is chosen as 0.0050 inch. The nominal conductivity for bottom layer is

81

taken as 33.5% IACS, and the range of the bottom layer conductivity is chosen to
vary from 30% IACS to 36% IACS. Figures 4.36 through 4.42 present results of real,
imaginary and magnitude value of the normal component of magnetic flux density for
different values of bottom layer conductivity. Table 4.4 shows the peak value of the
signal magnitude at different values of conductivity. Figure 4.43 shows the plot of

peak value of the signal magnitude verses bottom layer plate conductivity.

' (c)

Figure 4.36 Magnetic flux density (2 component) for 30% IACS bottom layer
conductivity (a) real, (b) imaginary, (c) magnitude

(c)

Figure 4.37 Magnetic flux density (2 component) for 31% IACS bottom layer
conductivity (a) real, (b) imaginary, (c) magnitude

 

(a) (b)

 

(a) (b)

82

     

(a)
Figure 4.38 Magnetic flux density (2 component) for 32% IACS bottom layer

(C)

conductivity (a) real, (b) imaginary, (c) magnitude

 
 

' (c)

Figure 4.39 Magnetic flux density (2 component) for 33% IACS bottom layer
conductivity (a) real, (b) imaginary, (c) magnitude

(a) (b) ' (c)

Figure 4.40 Magnetic flux density (2 component) for 34% IACS bottom layer
conductivity (a) real, (b) imaginary, (c) magnitude

(3) (b)

83

 

(a) (b) '0 ' (c)

Figure 4.41 Magnetic flux density (2 component) for 35% IACS bottom layer
conductivity (a) real, (b) imaginary, (c) magnitude

   

(a) (b) (c)

Figure 4.42 Magnetic flux density (2 component) for 36% IACS bottom layer
conductivity (3) real, (b) imaginary, (c) magnitude

_8 Peak Value of Magnitude vs. Bottom layer
05 X10 Conductivity

 

 

(else J.) epmgufiew need or need

 

 

 

1.97 - .
1.96 v r
19530 31 32 33 34 35 36

Bottom Layer Plate Conductivity (%IACS)

Figure 4.43 Peak value of the signal magnitude vs. bottom layer plate conductivity

84

Table 4.4 Peak value of the signal magnitude verses bottom layer plate conductivity

 

 

 

 

 

 

 

 

CONDUCTIVITY CONDUCTIVITY MAGNITUDE
(% IAC S) (SIEMENS/M) (T ESLA)
30 1.739167 2.03906-8
31 1.798067 2.03466-8
32 1.856067 2.02716-8
33 I 1.914067 2.02166-8
34 1.971067 2.01536-8
35 2.029067 2.00916-8
36 2.086967 1.98676-8

 

 

 

 

 

From the results presented above, it is seen that the peak value of the signal decreases
when bottom layer conductivity increases. When the bottom layer plate conductivity
increases, the induced eddy current becomes more concentrated around the bottom

surface and the observed signal on the top surface is reduced.

4.4.5 Conductivity of Fastener

The conductivity of fastener was varied as described in the parameter Table 4.1. The
operation frequency is chosen as 400 Hz and liftoff is chosen as 0.005 inch. The
different conductivity values for the fastener were chosen to be 1.0% IACS, 2.2%
IACS, and 3.1% IACS. The real, imaginary and magnitude values of magnetic flux

density at the GMR sensor for different conductivity values are plotted in Figure 4.44.

85

 

Figure 4.44 Magnitude, Real, and Imaginary components of magnetic flux density for
conductivity values of Ti fastener: (a) 1.0% IACS, (b) 2.2% IACS, and (c) 3.1% IACS

The variation of the peak values are summarized in Table 4.5. The fastener

conductivities considered in this study does not affect the signals with any

significance.

Table 4.5 Effect of fastemer conductivity on peak values of signal magnitude

 

 

 

 

 

FASTENER CONDUCTIVITY MAGNITUDE
(% IACS) (TESLA)
1.0 2.01856-008
2.2 2.0212e-OO8
3.1 2.01876-008

 

 

 

86

4.4.6 Fastener to Edge Distance

A commonly encountered problem in airframe geometry is the influence of surface
and subsurface edge on fastener and crack signals. The edge discontinuity behaves as
a large defect and generates its own signature that can affect the defect signal and
thereby lead to misinterpretation of GMR imaging data. This section describes a
systematic study of the effect of fastener-to-edge distance on the defect signal. The
edge is in the second layer and the nominal “fastener to edge” distance for test sample

is 0.6 inch. The range of values for parametric variations is 0.4 to 0.8 inch.

Simulations were performed at 400Hz excitation frequency and 0.0095 inch lift-off.

The normal component of the magnetic flux density 32 is calculated at the prescribed

liftoff and the resulting image data are plotted for each simulation. The real,
imaginary, magnitude and demodulated images are presented. The demodulated
image is derived using optimum detection angle (ODA) which is 70 degrees. Figures
4.45 to 4.49 show the real, imaginary, magnitude and demodulated images for varying

“fastener to edge” distance in the range 0.4 to 0.8 inch.

87

 

Figure 4.45 Real(a), Imaginary(b), Magnitude(c), Demodulated(d) magnetic flux density
(2 component) for fastener edge distance equals to 0.4 inch

 

Figure 4.46 Real(a), Imaginary(b), Magnitude(c), Demodulated(d) magnetic flux density
(2 component) for fastener edge distance equals to 0.5 inch

88

  

(c) (d)

Figure 4.47 Real(a), Imaginary(b), Magnitude(c), Demodulated(d) magnetic flux density
(2 component) for fastener edge distance equals to 0.6 inch

 

(C) ((0

Figure 4.48 Real(a), Imaginary(b), Magnitude(c), Demodulated(d) magnetic flux density
(2 component) for fastener edge distance equals to 0.7 inch

8

\O

 

(C) ((1)

Figure 4.49 Real(a), Imaginary(b), Magnitude(c), Demodulated(d) magnetic flux density
(2 component) for fastener edge distance equals to 0.8 inch

The line plots are extracted from all the image data and plotted. Figures 4.50 and 4.51
show the line scans across the center of the rivet image of real and imaginary parts.
Figure 4.52 shows the line scan of the demodulated signal using ODA. When the
“fastener edge distance” increases, the edge effect on the right side of the fastener is

decreased, so the magnitude of the signal on the right lobe will decrease slightly.

90

 

)(10'8 Fastener-edge-dist (Real. Bottom 03' crack)

1.5

     

 

 

 

 

—0.4inch
-1 _.. .. —0.5inch m
— 0.6inch
— 0.7inch
— 0.8inch

 

 

Signal (Tesla)

-1‘5”. ,. ,

 

 

 

 

 

 

 

33.02 -0015 -001 -0505 o 0.005 0.01 0.015 0.02
x distance (m)

 

Figure 4.50 Line scans across the center of the fastener for real part of magnetic flux
density (2 component)

 

x10'8 Fastener-edge-dist (lmag, Bottom 0.3' crack)
1 v Y Y ‘. '

    

 

 

Signal (Tesla)

 

 

 

"-3.02 .0015 .0101 -0305 0 0.005 0.01 0.015 0.02
x distance (m)

Figure 4.51 Line scans across the center of the fastener for imaginary part of magnetic
flux density (2 component)

91

x10’8 Fastener-edge—dist (ODA, Bottom 03" crack)

 

 

 

 

 

 

 

1.5 I
— 0.4 inch
— 0.5ind1
1 — 0.6inch ..
— 0.7inch
-\ — 0.3mm
0.5

.\. .............. (\

\\

 

 

 

Signal (Tesla)

 

 

 

 

 

 

 

 

 

 

 

 

 

112.02 -0.015 ~0.01 -0.005 0 01115 0.01 0.015 0.02
x distance (m)

Figure 4.52 Line scans across the center of the fastener for mixed signal using ODA

To further examine the effect of edge on defect signal amplitude, the signal from a
defective fastener for the largest fastener-to-edge distance is calculated and subtracted
from the defect signal, in each case, at the defect location. Assuming that there is
minimal effect on defect signal when the edge is the farthest, this value reflects the
effect of an edge in the proximity of the fastener. The defect signal amplitude

calculated in this manner is summarized in Table 4.6 and plotted in Figure 4.53.

92

x10"8 Fastener-edge-dist (Bo

ttom 0.3" crack)

 

 

Subtraction Signal (Tesla)

 

I I
o u

 

l'
I

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

i
045

i i
0 5 0.55

i
06 065 0.7

Fastener-edge—distance (inch)

Figure 4.53 Edge effect on defect signal amplitude: Real, Imaginary, and Mixed signals
vs. fastener-to-edge distance

Table 4.6 Peak to peak signal magnitude vs. fastener to edge distance

 

 

 

 

 

 

FASTENER TO REAL IMAGINARY MIXED AFTER
EDGE COMPONENT COMPONENT ODA
DISTANCE (TESLA) (TESLA) (TESLA)
(INCH)
0.4 -0.2593 6-7 0.3132 6-6 0.2430 6-6
0.5 0.1501 6-7 0.2369 6-6 0.2326 6-6
0.6 0.2031 6-7 0.1898 6-6 0.2352 e-6
0.7 0.1908 6-7 0.0810 6-6 0.0936 6-6

 

 

 

 

93

 

From Figure 4.53, it is Observed that the closer the fastener is to an edge the more
effect on defect signal. The signal due to the edge adds to the defect signal which

increases with proximity of fastener to the edge.

4.4.7 Crack Dimension

Four different trapezoidal crack dimensions were modeled including the case of a
crack-free fastener. The radial dimensions of the top of crack top were chosen to be
0.2, 0.22, 0.25 and 0.3 inches. The crack geometry was tapered towards the bottom
surface. The sensor liftoff was kept at 0.0095 inch and the frequency was chosen to be
500Hz. The normal component of the magnetic flux density Bz is calculated and
plotted. The real, imaginary and magnitude of the complex flux density are plotted in

Figures 4.54 through 4.58.

 

(a) (b) (C)

Figure 4.54 Magnetic flux density (2 component) for 0.2 inch crack (a) real, (b)
imaginary, and (c) magnitude

94

   
  
 

]
(b) (C)

Figure 4.55 Magnetic flux density (2 component) for 0.22 inch crack (a) real, (b)

imaginary, and (c) magnitude
3 J
(a) (B) (c)

Figure 4.56 Magnetic flux density (2 component) for 0.25 inch crack (a) real, (b)

(a)

     

imaginary, and (c) magnitude

  
 

 

I" . D
(a) (b) (c)

Figure 4.57 Magnetic flux density (2 component) for 0.3 inch crack (a) real, (b)
imaginary, and (c) magnitude

95

 

(a) (b) (c)

Figure 4.58 Magnetic flux density (2 component) for no crack (a) real, (b) imaginary, and
(c) magnitude

In order to see the variation in the signals with increasing crack dimensions, the line
scans across the center of the fastener data are plotted. The real part of the magnetic
flux density is shown in Figure 4.59 and Figure 4.60 shows the image of the
imaginary part. The demodulated signal using an optimum detection angle of 70

degree at 500 Hz is shown in Figure 4.61.

x10“8 Crack Dimension (Real)

 

 

 

 

 

 

 

 

'.s I l I
No defect
‘ t .20“ defect
.22” defect
.25” defect
0.5 ~ .30“ defect
0 ~ 4
-0.5 ~
.1 . a
4.5 . «
_2 1 I I I I i
-3 -2 -1 o 1 2 3 4
.4
x10

Figure 4.59 Line scans across the center of the fastener for real image data

96

x10’9 Crack Dimension (lmag)
8 . , , .

 

 

   
  
     
 

 

 

 

 

 

 

—l No defect

5 . — .20“ defect

.22” defect

4 > — .25“ defect

.30" defect

2 - ..

0 ~ .

.2 . _

.4 - -(

'6 ’ -4
.8 i i i I L
3 -2 -1 0 1 2

3
x10‘4
Figure 4.60 Line scans across the center of the fastener for imaginary image data

x10'9 Crack Dimension (ODA)
3 , , fl .

r

 

 

   

No defect
.20“ defect
2"“ defect
.25" defect
.30" defect

-2 -1 0 1 2

l

 

 

 

 

 

 

Figure 4.61 Mixed line scans using ODA of 70 degree

97

The defect signal amplitude is calculated using the peak-to-peak value of all three
signals, real, imaginary, and mixed. The defect signal values are also re-calculated

after subtracting the value obtained for the crack-free fastener. These values are

summarized in Table 4.7 and plotted in Figure 4.62.

Table 4.7 Peak to peak values of Magnetic flux density (2 component) for various crack
dimensions after substracting the crack-flee signal

 

 

 

 

 

 

 

 

 

CRACK REAL REAL REAL
WIDTH ON COMPONENT COMPONENT COMPONENT
TOP (TESLA) (TESLA) (TESLA)

(INCH)

0.2 9.896-010 -5.656—010 4.496-010
0.22 1.256-009 -8.586—010 7.746-010
0.25 1.5086-009 -1 . 1996-009 1.12596-009

0.3 2.6126-009 -2.l636—009 3.0055e-009

 

98

 

x10‘9 Peak to peak Values vs. Crack Dimension
4 . a . t
—0- Real
3 - —04 Imaginary .
—>— Mixed (ODA=70)

 

 

 

 

 

 

Peak to peak Signal (T)

 

l

0 0 02 0.04 0% 0.08 0.1

 

 

.3 1

Crack area (square inch)

Figure 4.62 Peak-to-peak values of Real, Imaginary, and Mixed MR signals vs. crack
area

4.4.8 Conclusions

The effect of parametric variations on the GMR measurements, presented in this
Chapter can be utilized in optimization of the GMR imaging system. The parameters
considered are frequency, sensor liftoff, conductivity of top and bottom layer,
conductivity of fastener, fastener to edge distance, and crack dimensions effect. Some
of these parameters are controllable, 6. g., fi'equency, liftoff, and other parameters are

uncontrollable (material conductivity, fastener tilt, etc.).

The study helps understand the effect of each parameter on the signal. As expected,
the model does predict that the measured signal decreases when the sensor liftoff

increases. Secondly, since the optimum detection angle is selected based on minimum

99

peak-tO-peak value in the image data Of the fastener, it is seen that in the mixed
signals the overall spread in the peak value with parametric variation is lower that the
variation in raw signal magnitude. This implies that the signal derived using optimum

detection angle is relatively insensitive to parametric variations.

Studies on conductivity of sample conductivity show that the peak value of the signal
increases when top layer conductivity increases. This is as expected since increase in
the top layer conductivity increases the induced eddy currents Of that layer thereby
increasing the amplitude of the measured signal. When the conductivity of bottom
layer increases, the induced eddy current becomes more concentrated in the bottom
surface and the observed signal on the top surface is reduced. The fastener
conductivities considered in this study does not affect the signals with any
significance. In the case of fastener to edge distance, it is observed that the closer the
fastener is to an edge the more effect on defect signal. The signal due to the edge adds
to the defect signal which increases with proximity of fastener to the edge. Results
from parametric study can be quantified in terms of skewness functions and

probability of detection curves.

100

\. Ii
'1:er --——-—‘ 'm

CHAPTER 5. INVERSE PROBLEMS IN

ELECTROMAGNETIC IMAGING

5.1 Introduction

The Chapter covers the various inversion techniques using in MO and GMR sensor
data. Sections 5 .2 to 5 .5 present the inverse problems for M0 imaging and Sections
5.6 to 5.7 present the corresponding results of inverse problems for GMR imaging.
The objective of inverse problems is to detect second and third layer cracks in

multilayer structures using MO and GMR sensor measurement.

5.2 Data Analysis for M0 Imaging

The analysis of MO imaging uses skewness functions defined for quantifying the
asymmetry of MO image. These fimctions are used in classifying a rivet as defect-free
or with defect. Classification algorithms were developed and optimized using
simulation data obtained via finite element models and the performance was
evaluated on experimentally measured MO images thereby substantiating the validity
of the skewness definitions and its usefulness in developing automatic MOI

classification system.

101

5.2.1 Automatic Rivet Detection

The overall approach for automated rivet classification is depicted in Figure 5.1. The
raw image obtained from the M01 image acquisition system is applied to the
motion-based filterer [32][33] to remove the background noise associated with
domain structures in the MO sensor. The filtered image is enhanced and thresholded
to obtain a binary image, which is used in the subsequent rivet detection and

classification modules.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

| Raw MO |
Images
1
Motion Based
Filter
.1
Rivet
Detection
Connected Circular . .
Rivet Image Split Rivet image
Edge Rivet Interior Rivet
v i i
Rivet Rivet Rivet
Classification Classification Classification
using 8; using 8. using 8;

 

 

 

 

 

 

 

 

Figure 5.1 The overall approach for automated rivet classification

102

The detection and classification of a rivet is based on the skewness value associated
with the binary image. The rivet center detection is performed using one of two
different approaches, namely Hough transformation and morphological image
processing [31]. All the classification results in this dissertation are obtained using the
morphological image processing operators of dilation and erosion with varying
structural elements [2]. Asymmetry quantification is first performed using skewness
feature of the data. The ideal rivet image is circular and hence its skewness must be
zero. If the computed skewness value is greater than some tolerance, the rivet is

classified as defective; otherwise the rivet site is deemed good.
5.2.2 Skewness Definition and automated classification

The skewness of the circular MOI rivet image and skewness quantification plays a
critical role in the automatic rivet classification. The skewness function should
capture the deviation from circularity of the rivet image, be robust under small
variations in the measured image and at the same time have discriminatory
information between good and defective rivets. This dissertation investigates several

definitions of skewness associated with an M0 image.

5.2.2.1 Skewness function S1

The first definition studied was designed to be independent of the frequency of

inspection and is expressed as:

103

Its:

Iii

Ei‘l

h

:0

f% R Z(Di_r)2
= 100 R+Br2+Z(Di—r)2

 

1

(5.1)
where f is the frequency parameter which is 50 kHz for the surface defects and 1.5, 3

or 5 kHz for subsurface defects. Here, B is the ‘width’ of the region outside the
rivet in the image, r is the average measured radius of rivet and {D1, i = 1,2,. . .n} are
the distances from the rivet center to the edge pixels. Finally, R is computed by
calculating the histogram of D,- and picked at its peak location. All these parameters

are illustrated in Figure 5.2.

 

Figure 5.2 Feature Parameters definition in skewness function S 1

5.2.2.2 Skewness function Sz

The second definition of skewness fimction that was studied is based on the

auto-covariance function, cov(D,-, Dj) of the series {Di }. The vector-valued skewness

function is derived from features in the auto-covariance plot and is hence defined in a

two or higher dimensional space for improved separability between classes. Consider
a set of n distances denoted by{ D1 D2 ,0 3 ,..., Dn }, as illustrated in Figure 5.3 in

which the circle represents the boundary of the binary MO image of a rivet.

Starting point for
”’7 {Di}

   
 
 
 
 

Starting point for Staging}: <})1nt
{Dj},where j > i D) ’’’’’’ ---------
/ Starting point
for {Du}

Figure 5.3 A set of boundary distances defined in skewness function Sz

The continuous form of normalized covariance definition is as follows:

(310(0) — r)(D(6 + 9') — r)d6
[02” D(6l)2 d0

 

cov(6l') =

(5.2)

105

The corresponding discrete form of normalized covariance definition is given by

213(6).) - r) - (0w.- + 0.) — 066.-
Cov(9k) = *‘=1

 

mnem-

i=1

(5.3)

where 19,, = 6k_1+§6k, 56?,- is the incremental angle between Di+1 and Di, and

the initial angle 00 is set to 0. Typical auto-covariance curves of the series {0,},

derived fi'om MO images are shown in Figure 5.4. Three skewness parameters are
derived fi'om the auto-covariance function, namely, Central Peak Ratio (CPR), Noise

to Signal Ratio (NSR) and Central Peak Width (CPW), shown in Figure 5.4. CPR is

the ratio between the central peak height Hcp and total height of the covariance curve
H T, defined as:

H
CPR = —CP
HT
(5.4)
CPW is the width of the central peak at its 50% peak value. NSR is the noise to signal

ratio defined fi'om continuous form of normalized covariance definition:

(D09) — r)(D((9) — r)dl9

[2%
NSR= 0 2
IO”D(9)2d6

 

(5.5)

106

which is also the first point of the covariance curve cov(0).

 

[ I V U V V Y ‘

 

~ NSR = cov(O)

 

 

 

 

 

 

 

 

——>CPW<—

 

 

Covariance plots for experimental images

 

 

 

   

 

 

 

\
.\
‘\
0
0
\
' ‘x
5‘ X 1 (S1=0.162) 2 (81:0.263)
‘\ i") t' o
\ g ‘ I a
\‘ " = {II
V" t.» 3 (31:0.317) 4 (31:03.59)

 

Figure 5.4 Parameters definition in skewness function Sz and corresponding results

Using these skewness parameters, a 2D feature space spanned by CPR and NSR or a
3D feature space spanned by CPR, CPW and NSR can be constructed, where the

coordinates in 2D space is (CPR, NSR) and in 3D space is (CPR, CPW, NSR). In

107

the 2D feature space, according to the definition and characteristic of the M01 binary
images, it is seen that for small defects the value of CPR is close to 1.0; also, the
value of NSR is proportional to the size (radial length) of the crack. The skewness

vector values in this case are defined based on the multi-dimensional coordinates,

2..

(CPR, NSR) for 21) Space
(CPR,NSR,CPW) for 3D Space

(5.6)

5.2.2.3 Skewness function S3

A third definition of skewness is based on the normalized central moments of
two-dimensional (2D) binary valued rivet images [2]. For a 2D digital image function

I (x, y), the discrete central moment can be expressed as:
(in. = 2 20¢.- - fem,- — wax, y)
x y

(5.7)
where (37c, 33c) is the center of the MO rivet image obtained using the morphological

image processing operators of dilation and erosion. For the binary valued rivet images

after enhancement and thresholding, I (x, y) is simplified to

x-, - is ' ' l
MAI My» we

0 otherwise

(5.8)

The normalized central moments, denoted Mpg, is then computed as

108

-111

M _
#700

P4

(5.9)

+4

 

where y=p +l,for p+q=2,3,....

The skewness function S3 was chosen from the various normalized central moments,

based on their discriminatory abilities, as the vector S3 in Equation 5.10:

S3 = {MzoaMzzaM4o}

(5.10)

5.2.2.4 Skewness Function Se

In the case of the distorted rivet image close to a structural discontinuity, we define

the skewness feature as:

 

 

lAr-A2|+le—le

S.=k<

)

(5.11)
where k is a scaling factor so that the skewness values for edge rivet images are in the

same range as that of interior rivet images. In the current classification procedure, we

set k = 0.5. A1 and A2 are the areas of the two largest lobes. Xm and Ym are the

dimensions of the largest lobe calculated from the center in the x and y direction

respectively, and X T is the total dimension of lobes in the x direction.

109

5.3 MO Image Classification

The rivet detection and classification algorithms were applied to experimental video
data obtained from scanning rivets on an aircraft panel using the M01 303 instrument
[27][28][29]. The video data was first segmented into single rivet images, which

were, then analyzed using the automated signal classification algorithm. The

classification results using the skewness functions 51 and 52 are shown in l, 2 and 3

dimensional skewness spaces in Figure 5.5. The skewness feature corresponding to

good rivets are represented by the circles and data points corresponding to cracked

rivet images by triangles. It is seen that overlap in the 1D space of S 1 is eliminated in

the higher order feature space. In the multi-dimensional skewness spaces the classes
are better separated and hence provide more accurate detection and classification
performance.

Figure 5.5 Classification results in multi-dimensional skewness space: (a) 1D Space (b)
2D Space (6) 3D Space

10 Skewness Space

 

 

 

 

 

 

 

 

 

 

 

""" ONo-CrackRivet(Good)"Ԥ' -------- i
.......... A Crack Rivet(Bad)
.......... ommAm
.. .......... Overlap ..........
———i— Classification Decision of S1 —

(a)

110

Figure 5.5 continued

 

2-D Skewness Space (CPR, NSR)

006

 

5 5 5 QANo-Ciack Rivet(Cvood) :
005- '''''' --------- --------- i--- A Crack Rivet(Bad) bi -------- i ---------

 

 

 

.0

O

b
l
I
i
:
1

Noise to signal ratio
2 8
to (i)
l i
i

a 2 a 2 E E %.
001—-------§ ......... ......... g--.fit“&-.é.----i ............. 9.-

~33

 

 

   

D i i . i i iQi
0 01 02 03 04 05 06 0

Central Peak Ratio
(b)
3-D Skewness Space (CPR, CPW, NSR)

 

09 1

.0.

| \
5
1 "~
0
I

‘\
I x
a l ‘
.-r 1‘
n‘. “ ' P‘
I s
“| I s
\I
‘ I

-
.....

 

.
n‘-
it"
If
.--
.-
40'
-9,
.--
______
-
I.-
C'-
-0."
--
.“
1"
--a
,-

0.06\

 

 

 

 

 

 

4“ Central Peak Width
(C)

111

2-D Normalized Central Moment Space (M20-M40)

 

 

 

 

 

 

 

 

 

 

 

 

2m ' .
o 1&1 -------1 ------- ------- ------- ------ -
£15m ______ i _______ i.-. O No-Creck Rlvet(Good)L _______ J _______ j ______ _
E E i I Crack Rlvei(Bad) . g
g 140 - ------- i ------ -
o E i E 5 E E E E E
2 120 -------- I~ ------- 1 ------- 1 ------- é ------- 5 ------- i ------- 1 ------- i ------- £---'---—
a = s 2 2 a a 1 i
‘5 1CD— ------ 1 ------- 1 ------- 1 ------- 1 ------- 1 ------- 1 ------- 1 ------- 1 ------- 1 ------ —
8 S E E E E E '5

1 . . . I- i
1: 3°" """ «i """"""" T
0 : : : : If :
1 1... ------ ~ ------- ~ ------- ;~ ------- 1 - , - . - ,.1 -- 1+1 ------- 1 ------- , ------- , ------ ~
g E ’ . E E
O 40",--.,--‘L-'.---i.-.-----r--.---.?.--....E...--.-E. .............................. -—1
Z :
20..----.-..------.-------.-------r-------.----.--§ .............................. _
9m 25) 240 200 l 300 350 :10 3150 300 400
Normalzed Central Moment M20
(a)
3-D Normalized Central Moment Space (M20-M40- M22)
"1' ~~~~~~~~
i = ~~~~~~~~~
. .L. ' I : ~~~~~~~~
' l ~~~~~~~ : -------
’7' ' J" i ......... I: "“~

21.;~ i - 1' - l : '~-.."~I- E E

, . v ...... , ,

0 No-Crack Rivet(Good).._‘1 l“ 1 g ........ 1

20” I Crack Rivet(Bad) E " v» i i
I g-T . ".7 : --“~-"‘

 

Normalized Central Moment M22

 

250
Normalzed Central Moment M40 50 nplormalized Central Moment M20
(b)

Figure 5.6 Classification results in normalized central moment space: (a) M20-M40
Space (b) M20-M40-M22 Space

The classification results using the normalized central moment vectors are shown in
Figure 5.6. The moments corresponding to good rivets are represented by the circles
and data points corresponding to cracked rivet images by squares. The advantage of

moment-based features is that it is intuitive and allows fast computation. These results

112

have also been extended to the calculation of Probability of Detection (POD) of
critical flaws along with the associated Probability of False alarm (PFA). The

definitions of POD and PFA for the binary image are defined as:

N
2a,.

P0D=i—=1—x100%
N

(5.11)
N *
26:,
N
(5.12)

where N and N * are the total number of data points corresponding to flawed

*
and good rivet images respectively and a,‘ , 051' are the corresponding

classification decisions of the automated rivet classification system. ai=1 implies

that the rivet image is classified as flawed where as (1,- =0 means it is a good rivet.

For real aircraft rivet images, statistical analysis was performed using the “POD
Version 3.0” software [36], which is a C-l-t program that interacts with a Microsoft

Excel 97 workbook to perform a POD analysis on the results with varying threshold

values of S1. When the value of $1 is greater than or equal to the classification

thresholds, the rivet site was classified as having a crack. POD plots with varying

classification threshold values are shown in Figure 5.7.

113

 

POD of S1

 

 

 

 

+ Tn: .150
Th= .210
+ Th= .230
—1— Th: .2675
—— Th= .310
0 0.05 0.1 0.15 0.2
Flaw Length (inches)

Figure 5.7 POD curves from response data of automated MOI inspections of surface layer
cracks around fasteners

POD curves are presented for four different threshold settings, which show how the
POD increases with reduction in threshold settings. However it should be kept in
mind that the PFA values also increases with increasing POD. The Receiver Operator

Curve (ROC) for different flaw sizes a is shown in Figure 5.8.

 

 

 

 

 

o 0.2 0.4 0.6 0.3 1
PFA

Figure 5.8 ROC curves for Automated MOI on surface cracks with different flaw size a

114

POD plots for hit/miss data from Automated MOI inspections of surface cracks are
presented in Figure 5.9 from the multi-flaw model fits. Note that both 2-parameter
and 4- parameter fits were performed. For S1 method, a threshold of 0.210 achieved a
0.9 POD value and a false call rate of about 1% for a flaw of length 0.094 inches. A

threshold of 0.150 achieved a 0.9 POD value and a false call rate of almost 21% for a

0.063 inches crack. The S2 method can achieve the same POD value and was shown

that it was effective in removing the false calls being made at edges. More image
processing and classification results are shown in Figure 5.10 to Figure 5.12. The

connected circular images are recognized as interior ones, while split images can be

recognized by the proportion of Al/Az (A1 3A2) as either rivet images close to the

edge when Al/Az is great than 0.4 or interior rivet images otherwise. Skewness

features S,- and Se are then applied respectively.

115

POD

— 31 call 2-parameter
— S1 call 4-parameter
—- 82 call 2-parameter

82 call 4—parameter

 

0 005 01 015 02 025 03 035 0.4

Flaw Length (inches)

Figure 5.9 POD fits for Automated MOI inspections of surface cracks. (POD for S 1 (dark,
red), for 52 (green, blue), both 2-pararneter and 4-parameter fits were performed)

 

Bad Good
(bl , (c)

I.

 

 

 

 

 

 

Figure 5.10 Three-step characterization for surface interior MO images (a) raw images (b)
enhanced images using motion-based filtering (6) detection and classification of images

116

 

 

Good
(3) ’ ~ ‘ (b)

 

 

 

 

 

 

Figure 5.11 Characterization for structural surface edge MO images (a) bad edge rivet
image detection and classification (b) good edge rivet image detection and classification

 

 

 

Figure 5.12 Characterization for subsurface interior MO images obtained under 5 kHz
frequency and various intensities of magnetic field: (a)-(b) are raw and classified images
(high intensity), (c)-(d) are raw and classified images (low intensity)

117

5.4 Inverse Problems of MO imaging Conclusion

This dissertation addresses some of the issues in automatic classification of MO
imaging data for aircraft inspections. Three definitions of skewness functions are
investigated to automatically quantify and classify aircraft defects. The moment based
method is computationally fast and the classification performance is as good as 81
method. The classification algorithms have been tested on experimentally measured
MOI aircraft rivet images and were shown to be robust. The performance of the
algorithm has been quantified in terms of the POD and PFA values.

The advantages of this method are that in addition to enhancing the MO images, the
technique also provides a quantitative basis for characterizing the images, which in
turn will yield consistent reject or accept decisions. The classification algorithm also
reduces the false calls due to the distortion of MO images near structural boundaries.

The overall automated classification system is implemented in real-time.

5.5 Inversion and Data analysis for GMR imaging

The inversion for GMR imaging data in the follow sections involves the image

classification, enhancement and characterization.

118

5.5.1 Optimization of Frequency Parameter in model predition - Skewness

Calculations

Although the defect indication is visible in most of the images, a simple method for
determining the optimum value of fi'equency is presented in this section via a
quantitative measure of the asymmetry in the two lobes of the image. A simple
skewness function for quantifying the asymmetry is calculated based on the peak

values of the two lobes of the fastener image, and is defined in Equation 5.13.

 

= max {I31 I,
a min{]Bl

le}
le}

 

 

 

9

(5.13)

The parameters of this fimction are shown pictorially in F igur65.13. The value of S0

is calculated for the fastener image obtained at each frequency and is plotted in Figure

5.14.

119

I31 — le

 

 

Figure 5.13 Surface plot of the image data showing the asymmetry in two lobes of
fastener image

Skewness vs. Frequency

Skewness S1

 

2 3 4
Frequency(kHz)

Figure 5.14 Skewness vs. Frequency for results with 0.3” crack

120

From Figure 5.14 it is seen that at lower frequencies the skewness values are much
higher than 1 and the crack is easily detected. At higher frequency, the crack is hard to
see due to the skin depth effect. The skewness plots show that the optimal operating
frequency is 0.4 kHz. Further, the modeling results indicate that as the frequency
increases, the images are sharper and the edge effect vanishes, but the crack detection

ability decreases because of the skin depth effect.

5.5.2 Optimization of frequency parameter in model prediction-SNR

Calculations

An alternate method for optimizing frequency of operation is to select a fiequency
that maximizes the POD or equivalently the SNR of the signal. In this study, SNR
values are calculated with simulation results obtained at frequencies: 100Hz, 400Hz,
SOOHz, 700Hz and lOOOHz. Both crack-free fastener and the case of a fastener with
0.3 inch subsurface crack were simulated. The fastener image at each frequency is
analyzed using the procedure described in section 5.5.1. The image data was
preprocessed using the optimum detection angle at each frequency and the line scans
across the center of the fastener was extracted. The MAE and MSE feature defined in
Equations 5.14 and 5.15 using experimental data from test sample are applied for
simulated signal.

1 N 1 M

FMAE = —Z LDIi]—— ZLND_j[i]
N i=1 M j=l

(5.14)

121

F 44:1 L [ii-Lg} [112
MSE N = D szl ND_]
(5.15)

F MAE is the feature calculated from mean absolute error and the F MSE is calculated
fiom mean square error. In Equations 5.14 and 5.15, N is the length of the mixed 1D
data vector, M is the total number of non-defective rivets. LD[i] is the ith pixel of

mixed data for different frequencies and LNDJ-[i] is set to zero vector as the baseline

signal. The line scans at different frequencies for the 0.3” subsurface crack signals are

shown in Figure 5.15.
x10'9 Parameter-Frequency (ODA. Subtraction)

 

  
  
   

 

 

5 ' 1 1
— 100 Hz f 5

5 """" 400 HZ """"""""""""""""""""""""""

4 ....... — 500 Hz ....... .................... .3 ................
— 700 HZ

3 ..... — 1000 Hz

 

  

 

 

 

 

Signal (Tesla)

-1

-2

 

 

 

.3 I 1
-0.01 -0.005 0 0.005 0.01
X Distance (m)

Figure 5.15 Simulated line scans at different frequencies for 0.3” subsurface cracks after
DA processing

122

The SNR definition used earlier is modified. The definition used in analyzing
experimental data is provided in Equation (5.16). The modified equation for analysis

of simulated data is given in (5.17), where the mean and variance are set to 0 and 1

respectively.

 

M 2
SNRM =\/Z(Fi ‘m0i)
i=1 0'01'

(5.16)
M 2
SNRM= 2E
i=1
(5.17)

where M = 2 in this case.

A scatter plot of 2D feature vector in the feature space is shown in Figure 5.16. A
classification rule can be devised based on partitioning the defect free fastener
(green dot) from the feature vector corresponding to the 0.3 inch crack data (red
triangle). The SNR value is proportional to Euclidean distance between (0, 0) and
the feature at each frequency. The SNR vs. frequency plots are shown in Figure

5.17 and it is seen that the SNR is maximized at 400Hz frequency.

123

 

x1 0'18 20 Feature Space for SIM data at diff. freq
I r i I

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

s «w 400 H2 11
5 "' """""" Is 500 H2
4 .. ................... _
12‘ s s 5
LL 3 ,,,,,,,,,,,,,,,,, p; 700 HZ ..... _
2 ----- n 1000 H2 1 -------------- -
1_ .................... i ................... 5‘ 100 HZ ...... i .................... _
03 015 i 115 2
FMAE x10'9

Figure 5.16 2D Feature Space for simulated data at different frequencies

124

x10‘9SNR

vs. Frequency for Simulated Data

 

 

 

 

 

2 ~ .
. ifs . :
i l E \ E s
. . \ . .
: I : : 1
1-3 r """"""""" [‘0‘ """""""""
6? E I t E 5
a 1' \s
616.. .............. .1, ............. . .............. 1,1, .............. .; ................
% ,.I §\
fl: 1 1 3 1
o .. ........... Li ............................... E ............. i ...............
,9 1'4 I : : ‘5 : '
2 I S E ‘k
g l E E E“
..- __ ........ In“: ______________________________ : _______________ 5.....fi
_ 1.2 , . , , ~~ .
1' ' : : :
.5; 1
1- ----- é ...... 1 --------------- 1 --------------- i --------------- 1 ................
0.8 i l 1 i
0 21:0 400 60) 8(1) 1010
Frequency(l-lz)

Figure 5. 17 SNR vs. Frequency for simulated signal for test sample with 0.3 inch

subsurface crack

From Figure 5.17, the model predicted optimum fiequency is 400 Hz, which is close

to the optimum frequency of 450 Hz obtained fiom the experimental signals for test

sample.

5.5.3 Magnitude Based Approach (MAG) data analysis for experimental data

A schematic of the test samples used in for collecting the GMR imaging data, are

shown in Figure 5.18. The wing and stringer parts are Aluminum, which are

fastened by flush head Steel (left two columns) or Titanium (right two columns)

rivets. The sample has 11 outside Titanium rivets and 10 inside Ti rivets with

125

various subsurface crack sizes as shown: 0.20 inch, 0.22 inch, 0.25 inch and 0.30
inch. The Steel side has the same layout of rivets and cracks.

The in-phase and quadrature components for test sample obtained at 400Hz are
shown in Figure 5.19. The segmented in-phase images at 400Hz for outside
Titanium rivets are shown in Figure 5.20 for different crack sizes, from
defect-free to the largest one, 0.30 inch subsurface crack. As we can see, it is
difficult to differentiate cracks in the raw images.

A simple intuitive Magnitude Based (MAG) approach involves deriving the
corresponding “zero-mean” signals (Figure 5.21) in which the progression of the

cracks is observed from the shape of the right lobes.

126

 

A 0.30 inch
0.25 inch

0.22 inch

0.20 inch

. O? O (P OCR? O OO 1
‘O Q 0 Q O C? O C) O O O

 

 

 

 

 

Figure 5.18 The test standard schematic

127

 

7.: .
C _-

Figure 5.19 Raw GMR output images for the test standard at 400Hz: (from left to right)
Steel in-phase, Steel quadrature, Titanium in-phase and Titanium quadrature components

1"” _ A __ . .
CD @0310 CD @Dl
(a) (b) (C) (d) (C)

Figure 5.20 Raw output images for the test standard at 400Hz: (a) defect-free, (b) 0.20”,
(c) 0.22”, (d) 0.25” and (e) 0.30” subsurface cracks

 

 

 

   

(a) (b) (C) (d) (C)

Figure 5.21 MAG images at 400Hz: (a) defect-free, (b) 0.20”, (c) 0.22”, (d) 0.25” and (e)
0.30” subsurface cracks

128

n. 1111113

 

(c) (d)

Figure 5.22 Definitions of features to quantify the MAG signals: (a) F; and F2, (b) F3, (6)
F4 and (d) F5

Several features are extracted from the line plots across the center of the images
in Figure 5.21. Equations 5.18 to 5.22 give the mathematical forms for those

features, which are illustrated in Figure 5.22.

W .
F1 = J‘- (5.18)
i=1 WLi
N
F2 = 2(WRi ‘ WLi) (5.19)

129

Area __ (Xs-X2)'(Y3 ‘Y2)

 

_-_-____ (5.20)
3 lSlopeI IY3 —Y4I/|X3 —X4I

F4 = HMHW (5.21)

FS=HTaiI=Y4—Y2 (5-22)

The two dimensional scatter plots for classifying the test sample Ti rivets using the
above features are shown in Figure 5 .23. Although the MAG signals for outside and
inside Ti rivets overlap, the features can separate the defective and defect-free data.
However, this approach fails in the case of Steel rivets data due to the strong signal

generated by the magnetic steel rivet.

Figure 5.23 MAG signals for the test standard Ti rivets and 2-D classification scatter
plots

 

2-D Feature Scatter for CT. S-2

 

 

 

 

 

 

 

0.12 . . .
0.115 - - ' ' ' ' ...... .......... i ......
Circle-Defect free g V
0"" Triangle-Defective --------- ~
o.1os~ . . . . ...... .......... ......... ..
(‘0 s s s a a s
(D 0.1 - ......... ..
" = = = = s s
a 0.095 .. ......... ---------- --------- 1
g a i s s s
0.09 .......... .......... .......... p ......... _
0185- --------- 3. ---------- >--1 ---------- é-p ----- «f ---------- i --------- -
emf. ---------- --------- —
0.075 1 i i i i
2

Feature 1
(a)

130

Figure 5.23 continued

 

2-D Feature Scatter for N. S-2

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

0.055
0.1 ~11 Circle-Defect free P
Triangle-Defective
<1- 0.045 -...
9 5
$3 0.01%,
u. .
0.035 »
0.03 --------- C.
0.025 i
0.2 0.25 0.3 0.35 0.4
Featurez

(1))
MAG Line Plots for OT. 32

Dash-Defect free
Solid-Defective

 

2.6 2.8 3 3.2 3.4 3.6 3.8 4
(C)

131

Figure 5.23 continued

MAG Line Plots for IN. S-2

Dash-Defect free
Solid-Defective

012

01

008

005

004

002

 

5.5.4 Optimum Detection Angle (ODA) Based Approach for experimental data

The detection angle based approach derives a weighted sum of the in-phase and

quadrature components according to Equation 5.23.

 

 

< (S, , SQ ), (Cos 6?, sin (9 >= [AB/1” cos((p,, ):| cos 19 + [A024, sin((o,, )] sin 61

(5.23)

where 0 is the DA. The operator <*,*> is the inner product.

 

(a) (b) (C) (d) (C)

Figure 5.24 GMR mixed images using ODA = 70 deg.: (a) defect-free, (b) 0.20”, (c)
0.22”, (d) 0.25” and (e) 0.30” subsurface cracks

132

The right hand of equation 5.23 is the “mixed” EC-GMR data. The optimum

detection angle (ODA) a is calculated as:

a = arg min(Spp (6))
0

(5.24)

where S pp ((9) is the peak to peak value in the mixed signal.

The mixed GMR images after applying the ODA a in Equation 5.23 are shown in
Figure 5.24. Figure 5.25 showing the central line signals across the images shows
significant separation between defective and defect-fies rivet data. To quantify the

difference between defect-free and defective mixed signals, two features, mean
absolute error feature (F MAE) and mean square error feature (FMSE) are extracted as

Equations 5.14 and 5.15 in the previous section. In Figure 5.26, the 2D scatter plots

clear show better separation than that obtained using MAG approach.

Figure 5.25 Central line plots of images in Figure 5.24 (a) raw image-inside rivet (b) raw
image-outside rivet, (6) mixed image-inside rivet ((1) mixed image-outside rivet

RAW Line plots for OT. S-2

Dash-Defect free
Solid-Defective

 

133

Figure 5.25 continued

 

RAW L'ne plots for lN. S-2

I i r

: Dash-Defect free
. Solid-Defective
0 as — ----- , _

 

 

 

  

 

 

 

 

  

 

 

 

 

 

 

 

 

 

 

 

 

 

(b
ODA Line plots for CT. S-2

Dash-Defect free
Solid-Defective

 

(C)
ODA Line plots for IN. 32
0.025

Dash-Defect free
Solid-Defective

 

134

 

x1 0‘4 2D ODA Feature Scatter for OT. 8+2
1.5 .' : :r i 4F"—

 

 

Circle-Defect free
Triangle-Defective

 

 

 

 

 

 

 

 

 

 

0.2 h
.0 .1?

 

 

 

 

 

 

O

 

0.004 0.” 0.008 0.01 0.012
FM°E
(a)
x104 2D ODA Feature Scatter for OT. S-2
1.5 1 . .
‘ P

 

 

Circle-Defect free
Triangle-Defective

 

 

 

 

 

 

 

 

 

 

 

 

 

0—0

0 0.032 0.004 0.0116 0.008 0.01 0.012
l:MAE
(b)
Figul‘e 5.26 2D classification scatter plots for ODA approach: (a) OT and (b) IN

135

5.6 Discussion and Conclusion for Inversion of GMR imaging

The performance of MAG and ODA approaches is compared quantitatively using

a signal to noise ratio SNR], metric defined as

0.12

0.115

0.11

0.105

0.1

0.095

Feature 3

0.09

0.085

0.08

0.075
0.2

 

L 2
SNRL = 2(Fi _m0i)
i=1 0'01

(5.25)

0'03 in F 3 direction

. 0'02 in F direction

 

0.25 0.3 0.35 0.4 0.45 0.5 0.55
Featue 2

Figure 5.27 Illustration of SNRL definition

where L is the dimensionality of signals, m0 and 0'02 are the mean and

variance of the defect-free data cluster as illustrated in Figure 5.27. The SNR

136

in)

‘11

comparison between conventional ECT and GMR approaches is shown in Table
5.1. Dramatic improvement has been observed by applying the automated data
analysis methods: MAG and ODA methods. This dissertation develops two
innovative data analysis approaches for enhancing the GMR imaging data by
increasing the SNR and hence maximizing the POD. Using these methods, GMR
sensors show potential for detecting subsurface cracks in multi-layered
structures. Automated crack detection is performed via optimum feature

estimation and simple clustering.

137

Table 5.1 Signal to Noise Ratio Comparison between conventional ECT and GMR

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

approaches
TEST SIGNAL TO NOISE RATIO (EQ. 5.25)
STANDARD
Thickness Row Methods 0.30 inch 0.25 inch 0.22 inch 0.20 inch
0.16” Outside ECT 3.00 1.50 1.50 1.50
MAG 45.0 18.0 14.0 9.60
ODA 317 49.8 16.4 17.1
0.16” Inside ECT 2.50 1.50 1.00 1.00
MAG 20.6 8.20 4.40 4.00
ODA 292 64.5 17.1 18.2
0.25” Outside ECT 2.00 l .00 1 .00 1 .00
MAG 30.1 11.4 7.70 5.50
ODA 114 13.3 4.60 11.4
0.25” Inside ECT 2.00 1.50 1.00 1.00
MAG 7.20 4.60 2.50 3.30
ODA 44.1 9.70 4.10 5.10

 

 

 

 

 

 

 

138

 

“Ir-I'll

. , tun”
p
.

CHAPTER 6. POSITRON EMISSION TOMOGRAPHY

IMAGING

6.1 Introduction

Positron Emission Tomography (PET) [l][3][7][8][15] imaging is a functional and
metabolic imaging technique that generates images showing the functional details in
patients by depicting the distributions of positron emitting radioisotopes.

. . 11 15 18 13 . .
Neutron-deficrent Isotopes, 6. g. C, O, F, N, that are incorporated into

metabolically relevant compounds can decay by emitting positron, i.e. fl+ decay.

Two 51 lkeV photons are emitted in anti-parallel directions after the positron has been

emitted fiom the ,B+ decay and subsequently annihilated with an electron. If both of

these annihilation photons are detected by the annihilation coincidence detection
(ACD) of the PET system, the projection data is sorted in the Radon space, i.e.
Sinogram, and reconstructed into a, transverse slice of the image using filtered
back-projection algorithm or the iterative methods [1][3][15]. Thus ACD sorts
detected photons in the Sinogram list-mode data, establishes the trajectories of the
detected photons and finally forms the PET images. Modern PET scanners are
multi-slice devices; permitting simultaneous 2D/3D acquisition of multiple slices over

156m to 206m of axial distance. PET imaging has proved to be effective and valuable

139

in clinical work and research particularly in cardiac and brain imaging, cancer
detection and delineation and small animal studies [1][3][15].

The schematic diagram of a PET imaging with its physics is show in Figure 6.1. More
details of the PET imaging system and a brief history of PET in medicine and

positron-emitting tracers are discussed in the following sections.

 

 

 

 

 

 

Figure 6.1 Schematic diagram of the PET imaging system with physics of PET
6.2 History of PET in Medicine

Positron emission from radioactive nuclei was discovered by Thibaud and Joliot in
1933 [10][11]. Subsequent research showed that alter positron emission, two photons
were emitted simultaneously at nearly opposite direction to each other [12]. The
potential importance of positron-emitting radionuclides in medicine was suggested in

1946 [13]. The first use of positron in medicine is in the study in 1945 by Tobias et a1.
using llC-CO in single photon mode [14] was to conduct an in-vivo studying of

carbon monoxide. The Massachusetts General Hospital (MGH) first in -1951

140

published an image of positron-emitting tracer study in man for localization of brain
tumor [15]. The coincidence detection techniques were also employed at the same
time.

In the late 19503 studies of metabolism in murine tumors with autoradiography, lung
ventilation oxygen metalbolism in human were performed at the Hammersmith
Hospital in London [16][17]. The inventor of the scintillation “Anger” camera, Hal
Anger discussed the practicalities of coincidence detection in 1959 [18]. A circular
“section scanner” was developed at Brookhaven National Laboratories also in the late
19508 [19] for localizing brain tumors, which appears to be the first device used to
measure regional cerebral blood flow with positron emitters [20].

Due to the limitation of computational techniques, tomographic images and
reconstruction algorithms were not developed until a decade later Kuhl 6t al. in
Philadelphia explored transaxial SPECT in the 19603 and early 19703 [21][22]. Their
approach was based primarily on reconstruction by back projection without
pre-filtering compensation for the inherent blurring effect. Chesler adopted a filtered
back projection approach to produce transaxial tomographic reconstructions in 1973

using data collected from the rotation MGH positron camera, presented convincing
transaxial distributions of l3N-NH3 in the dog heart [23]. The development in Single
Photon Emission Computed Tomography (SPECT): Anger cameras and scintillation
probes helped in purchasing a PET scanner later. In 1974, a dedicated single-plane

Positron Emissing Transaxial Tomography (PETT) was developed at the Mallinckrodt

Institute of Radiology, St Louis The name of PETT was shortened to PET later.

141

Ter-Pogossian, Phelps, Hoffman, and Mullani had re-addressed the physical criteria
that determined the accuracy of PET tomography [24] [25] and optimized the physical
design in order to minimize registration of scattered and random coincidences, and
also to ensure adequate spatial sampling. In 1978, Phelps and Hoffman collaborated
with Douglas and Williams and designed a commercial, single-plane PET scanner
known as the Emission Computerized Axial Tomography (ECAT) [64]. In the late
19703, PET made its biggest impact in realizing absolute quantifiable tomographic
physiological data by minimizing the recording of scattered coincidences as well as
random events. It is realized by using parallel delayed coincidence circuits [26]. In the
early 1980s, a number of new radioactive tracers were introduced and most research
was focused on 18F and 11C tracers. Till the year 1985, a number of different PET
scanners had been developed including whole body scanners [15]. In the late 19808,
when block detector design was used, septa were placed between the transaxial planes
in order to reduce the registration of scattered coincidences as well as the random
ones. The dead time loss was also overcome by the shielding technique. Nevertheless,
the counting statistics lost because of the small solid angle for coincidence detection
and the sensitivity of 3D PET was restricted. The “open” mode PET was implemented
by removing the septa and all possible coincidences were recorded. The basic
sensitivity was increased by a factor of five while it brought new challenge for
reconstructing the full volumetric 3D PET data. The value of 3D PET was evaluated
using the criteria of Noise Equivalent Count (NEC) rate, as had been derived by
Strother [63]. The use of 3D PET reduced the radiation dose required to achieve a

significant statistical result. In 19905, 3D PET/CT system was introduced and thought

142

 

to be a fundamental advance. As we see in the rich history of PET, a wide range of
clinical research studies have been on going extensively but there still remains many
challenges. Three different generations of PET scanners are shown in Figure 6.2 [80],

in which the last one is the PET/CT fusion scarmer developed by GE (Milwaukee,

WI).

Figure 6.2 Three generations of PET scanner: the last is the PET/CT fusion system

GE Advance PET scanner

 

143

Figure 6.2 continued

  

(C)

6.3 Physics and instrumentation of PET

PET and SPECT are both nuclear medicine imaging techniques. Similar to SPECT,
PET is used to measure the functional and metabolic details of human body instead of
anatomy. The main difi'erence between PET and SPECT is that the injected or
ingested radioisotopes used in PET emit positrons and generate two y rays in nearly
opposite direction after positron-electron annihilation. Two photons are detected in
PET instead of the single photon in SPECT. Consequently PET offers higher SNR and
spatial resolution than SPECT.

The positron emitting radioisotopes must be synthesized using a cyclotron, and are
structural analogs of a biologically active molecule, such as glucose, in which one or
more of the atoms has been replaced by a radioactive atom [1]. Because of the short
half-life of these nuclides, the cyclotron must be on-site to produce the radioisotopes,

which is seen as one of the disadvantages of PET system. The most used radioisotope

is 18F-fluorodeoxyglucose (FDG) due to its relative long half-life, which is 109.7
. 11 15 13 . .
minutes. Isotopes such as C, O, and N can also undergo radroactrve decay by
emitting positrons. Take 18F as an example, the radioactive decay is like Equation 2.1:
198F—>1880 + ,B+ + v
(6.1)

where fl+ is a positron and V is a neutrino. The radioactive decay is illustrated in

Figure 6.3.

 

Anti-neutrino

ll/

\Positron
fl 4-

 

 

 

 

Figure 6.3 Illustration of radioactive decay in which the positron is generated

The positron travels a few millimeters in tissue randomly before it annihilates with an
electron and results in the formation of two y rays (or simply say photons) in nearly
opposite direction with energy of SllkeV each. The annihilation is expressed in

Equation 6.2:

145

,6 + + e” —-) y + y

(6.2)
A PET system has a complete ring of scintillation crystals, usually Bismuth
Germanate (BGO), surrounding the patients emitting 51 lkev photon pairs, which is
illustrated in Figure 6.4. The dashed lines in Figure 6.4 represent the scintillation
crystals detection ring with ACD, which forms the signal localization in the system.
The green line is referred to as the line of response (LOR) along which the
annihilation must have occurred. Ideally, two simultaneously emitting photons are
detected at the crystals at the two ends of LOR within a certain time window (typical
time window for BGO detectors is 12 nsec [3]), which is referred to as a coincidence
event. Other issues, such as scatter, random noise, that introduce errors in this
process, are discussed in detail in the next chapter on modeling of PET imaging
system. The schematic of the PET detector structure with a detector element [80] is
shown in Figure 6.5. The Gamma ray (Gamma photon) hits the scintillation crystal
and is converted into optical photons, whose energy is proportional to the gamma ray
energy. The photons are collected at the 6nd of the crystal and passed into
photomultiplier tubes (PMTs), where the light is converted to an electrical signal and

amplified.

146

      
   

   

 

I ~ \
I”’ ~\ “
I, I ‘\ \
I ” SllkeV \ \
’0’ photon \\\

I I \ \
I ’ \ \
I I . . . \\
, l Annihilation ‘ ‘
1: 13+ 11
h ‘ Posrtron l:
11511k6 :1
\ ‘ photon I I
\\ , I
\\ I I
\\ I I
\\ I I
\ \ I I
x \ I
x \ a I

\ s ‘,
\"~~-_-" ’/

Figure 6.4 Illustrations of the overall photon emission and detection process of PET
scanner

In contrast to SPECT, which requires collimation of single photon, PET gives much
higher detection efficiency [1] than SPECT. Furthermore the energy of photons in
PET is 511 k6V, which is much higher than the 140 keV in conventional nuclear
medicine. This in turn implies less attenuation in tissue and hence results in higher
sensitivity of PET system. The resolution of a PET system depends on a number of

factors and the typical values of the spatial resolution are 4.5 to 5mm [3].

147

, Gamma Ray
~. (Photons)

O
\
I

I

‘0

N ,u

 

 

Scintillation Crystal

1

PMT

l

Pro-amplifier +
Electronics

i

Computers and
Image Processing

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Figure 6.5 Schematic of the PET detector structure

To understand the physics of PET better, Table 6.1 [3] compares PET and SPECT

systems in detail.

148

Table 6.1 Comparison between PET and SPECT

 

 

 

 

PET SPECT
Principle of data collection ACD Collimation
Image reconstruction Filtered back projection Filtered back projection
(FBP) or iterative methods (F BP) or iterative methods
Radionuclides Positron emitters only Any emitting x-rays,

gamma rays, or

annihilation photons

 

 

Spatial resolution

Relatively constant across
transaxial image, best at
center. Typically 4.5 to

5mm FWHM at center

Resolution in the radial
direction is relatively
uniform but the tangential
resolution is degraded
toward the center.
Typically 10 mm FWHM
at center for a 30 cm

diameter orbit and

 

 

Tc-99m.
¥ Attenuation More severe Less severe
Cost 1 Million to 2 Million Typically 0.5 Million for

 

 

double-head variable-angle

 

149

 

6.4 Challenges in PET imaging

PET is widely used due to its uniqueness and is becoming more and more important
in clinical and research work after its invention. However PET image acquisition
takes a long period of time (usually greater than 10 minutes) for each scan while the
patient is breathing. The respiratory motion that includes pulmonary and cardiac
motion introduces severe artifacts in the final reconstructed PET images, as described
further in the following sections. Hence for accurate interpretation of PET data, these
artifacts have to be eliminated. Extensive motion correction techniques are
investigated and can be found in the following literatures [82] [83] [84] [85][86]

[87][89][9l] [108][111]. We will discuss these issues in the next two Chapters.

150

CHAPTER 7. FORWARD MODELS FOR POSITRON

EMISSION TOMOGRAPHY IMAGING

7.1 Introduction

In this dissertation Monte Carlo (MC) methods are used for solving forward problem
in PET imaging. MC simulation methods are finding increasing applications in
nuclear medical imaging systems such as PET and SPECT [65] [66] [69]. The
computational PET model developed in this dissertation is based on MC methods.
The name “Monte Carlo” was coined by Von Neumann [66] and chosen during the
W'VVII Manhattan Project [70] because of the similarity to “random” properties of
casino games based on chance. They have been applied to simulate a process with
random behavior and in which the physical parameters are hard and even impossible
to measure and quantify. The availability of an accurate computational model using
MC technique helps immensely in understanding the undesirable effects of several
opel‘ational parameters such as blurring effects of patient movement due to respiratory

or other physiological motions in PET.

In a statistical sense, the Monte Carlo method can be explained simply as an

evaluation method for an integral:

I[f(x)] = 3%; [flow

(7.1)

151

Uniform random variables X1, X2,..., XN are generated in the interval (a, b) and

substituted into:
. 1 N
I[f(x)] = NZflXi)
i=1

(7.2)

Based on the Law of Large numbers (LLN), we get:

1
b-a

 

it/(xn z E[f(x)] = [:f(x)° dx = 1mm

(7.3)
In medical physics applications, the systematic development of MC methods did not
start until 1948 by Kahn and Harris [71]. Fermi, Metropolis and Ulam used MC to
estimate the eigenvalues of the Schrodinger equation during the same year. The
applications of MC in this field increased after the review paper by Raeside [72] and
several books [70] [73] [74]. The methods were applied to all aspects of nuclear
medicine imaging, including Planner Gamma Camera (PGC) imaging, SPECT, PET,
and multiple emission tomography (MET) [66]. However, MC methods did not find
widespread use until high performance computers became more common in recent
years.
In the PET system described in this dissertation, due to the complex detector
geometries, inhomogeneous attenuation within the body or phantom structures and
stochastic nature of nuclear medicine systems, analytical solutions are seldom
possible and an accurate, fast, and efficient numerical model is necessary. This

chapter first presents a 2D MC model, which simulates the radiation emission,

152

transport, and detection process of PET systems to generate data in the Sinogram
space. The reconstructed images are then obtained using filtered back projection
(FBP) method [1] [75] to validate the model. A 3D MC model using the SimSET [68]

is then introduced.

7.2 Monte Carlo Methods in PET

Monte Carlo (MC) methods are statistical methods that use random numbers to
simulate any specified situation and physical process. The major components of MC
method include

1. The probability density functions (pdf’s) of the physical system, which is PET

system;

2. Random number generators;

3. Sampling rules for the specific system pdf’s;

4. Tallying or scoring methods to accumulate the outcomes and

5. Error estimation as a function of the number of trials and other quantities.
Details of the general principles of MC method are well established and described in
several publications [70][72][73][74]. The following sections present a discussion of

the general principles.
7.2.1 System Model

In a PET scanner system, the relationship between the measurement xi}- and the

object f (7c) can be described by a linear model that represents the summation or

153

integration of spatially varying detector sensitivity function SOC") [76] representing
the object. The sensitivity function is different for scatter-free or with scatter and
random noise situations. The detector measurement can be represented
mathematically as follows:

N

22,. = [FOV f(x).s(r)dr
j=1

(7.4)

Where )7 is Rm , m = 2 or 3 for m dimensional object function, FOV is the field of
view. In PET image reconstruction, a number of approaches have been presented [77]
to estimate the object f (3?) by assuming that a discrete number of PET scanning
' measurements are obtained accurately. In this Chapter, we present a forward
probabilistic model that uses a known sensitivity function s(3c') to predict the

simulated measurements which are then used as input for image reconstruction.
7.2.2 Random Number Generator

The random number generator (RNG) is an important and fundamental part of MC
simulation models, in which the RNG is used to control decision making for physical
events having a number of possibilities. Generally, the sequence of random numbers
used in the MC model should have the following properties: uncorrelated, long
period, uniformity, reproducibility and speed [66].

The sequences of random numbers are uncorrelated and independent. Actually all

RNGs based upon mathematical algorithms are “pseudo-random” and are repeatable.

154

However the period of the RNG is long enough and has the appearance of
randomness. The repetition only occurs after a very large set of random numbers are
generated. Also the sequence should be uniform and unbiased. The reproducibility is
important because we need to repeat the simulation when debugging the programs or
transfer the programs to a different machine. Speed is of course a concern to make the
whole simulation efficient and fast.
A large number of RNGs are available for implementation [78]. The most commonly
used RNGs include Linear Congruential Random Number Generator (LCRNG),
which is also tested in the dissertation, and Lagged-Fibonacci Random Number
Generator (LFRNG). The LRRN G and LCRNG formulations are presented in
Equation 7.5 and 7.6.
RN.“ =(RN.-1_1 ®RN.-k-1)mod(c),I > k

(7.5)

where the operation ® may be one of the binary arithmetic operators, 1 and k are

the lags and c is a power of 2.

(Cl-RN" +b mod c
RIvn+1= ) ()

 

c
(7.6)
where a is the multiplier, b is the additive constant or addend and c is the

modulus. The initial random number RN 0 is the seed. If b is zero then the

LCRNG is called a multiplicative congruential random number generator

(MCRNG). The pseudo-code of MATLAB implementation for LCRNG is as follow:

155

IF n = 0

RN(n+1) = SEED (real number from [0, 1])
ELSE

RN(n+1) = LCRNG_iter(RN(n));
END
Function RN_n+1 = LCRNG_iter(RN_n)

A = 7000; B = 55000; C = 260000;

RN_n+I_temp = mod(floor(RN_n) *A +B, C);

RN_n+I = RN_n+I_temp./C
End
In the 2D model, a uniform RNG is used in deciding the rotation angle of LOR 9
and photon path length d before scatter is used. The uniform RNG outputs are also
used in deciding the annihilation position while sampling the cumulative distribution
function (cdt) calculated based on the input 2D phantom intensity. The uniform RNG

in MATLAB was also tested in the current model with good results.
7.2.3 Sampling Techniques

In MC simulation, the pdf’s of the PET system needs to be randomly sampled so that
the physical events can be described and evolution of the overall system can be

simulated. The illustration of the basic idea of sampling is shown in Figure 7.1.

156

 

 

 

 

 

 

 

 

 

 

PdeU)

Figure 7.1 Illustration of sampling the cumulative distribution function (cdf) for
simulating the annihilation event

The sampling methods include direct method, rejection method and mixed method
which combines the two above methods [70]. In this dissertation, the simplest direct
method is used since the inverse of the cumulative distribution function (cdf) for
determining the annihilation position is easily obtainable using the true image
intensity distribution. This cdf is then uniformly sampled to provide the pdf of the
region where the annihilation occurs (discretized pixel index in this 2D model). The
corresponding cdf for determining the annihilation position is calculated according to
the 2D phantom intensity and distribution. This is reasonable because the more the
radioactive uptake concentrated in the distribution; the higher the annihilation density.

The cdf of the photon pair emitting angle is simply a straight line because the

probability of emitting direction is distributed uniformly in[0° ,3 60°] .
7.3 Two dimensional (2D) model

In all the simulations described in Chapter 7, we assume that the patient is at rest, i.e.

by holding breath and the simulated projection data will be consistent with the output

157

obtained using Equation (7.13). The simulated projection data are called
“motion-free” Sinogram. The “motion-encoded” sinograrn [82] is generated to
simulate the patient data with thoracic motion and these details are discussed in the
following Chapter.

In the probabilistic model for motion-free case, the object f (3?) is a 2D phantom,
which is defined as a nonnegative intensity function p(R, 0) and discretized in a
bounded region 0. For the purpose of simplicity of mathematical calculations, all
functions and equations are derived in the Radon space (RS) 93(R, 6) coordinates.
The 2D phantom p(R, 0) is treated as a known system input with intensity p,- , i =
1, 2, ..., M at each pixel, where M is the total number of pixels in the bounded

regionQ. Any detector pair D J- , j =1, 2, ..., M is defined in the RS coordinates as a
vector (til-1,6712) , where 671-1 =(—Rj,9j) and 671-2 =(Rj,0j+fl') . The

geometry and the new RS coordinate introduced in this dissertation with the
discretized region and the detector pairs are illustrated in Figure 2(a). There are three
different cases considered in the probabilistic model, namely, i) scatter-free, ii) single

scattering and iii) random noise. Generally, in all three situations the relationship

between the measurement xlj of jth detector pair and the 2D intensity function

p(R, 6) is defined as
21,- = llp(R.-,6i)p,-,-de6
1' S)

(7.7)

158

where pij , the specific form of the detector sensitivity function 5(3?) in our model,

is the probability of a 511keV photon pair generated by annihilation at the ith pixel

and detected by the jth detector pair. The value of pi}- is calculated as the

exponential of a line integral of the attenuation coefficient of each pixel along the

photon emission path as follows:

1 ‘71'2
Pi,- = 76X“- _l #(Ri: 90611)
m: djl
(7.8)

where the attenuation coefficients ,u is also predefined in the bounded region Q ,

shown in Figure 7.2(a) and r is the radius of the PET scanner. However, the sensitivity
function in Equation 7.8 for different scenarios is calculated differently according to
the apparent lines of response (LORs) for different situations. By applying MC

technique, the integral in the Equation 7.7 is evaluated by generating a large number

of trials and then evaluating the summation. The Sinogram data is generated by ’1 J-

and then used to reconstruct the phantom using the FBP method. The illustrations of
scatter-free, single scatter and random noise models are shown in Figure 7.2(b) to
7.2(d). In the single scattering case, the photon path length d before the scattering is
calculated using MC technique according to Equation 7.9, where U is a uniform
random number in (O, l).

d = —-1—ln(l—U) = —-1—1nU
# .U

(7.9)

159

The scattering angle 05 and the photon energy after scattering E, are calculated
based on Klein-Nishina cross-section equation and the sampling algorithm originally
developed by Kahn [71]. Em and E s , the energy before and after the scattering and

the photon path I, afler the scattering are calculated through Equations 7.10 to 7.12.

 

Em = 511e_’7dkeV (7.10)
E, =511keV/(2—cost93) (7.11)
I, = M cos2 a, —d2 +R2 —dcos0, (7.12)

Figure 7 .2 (a) 2D MC Model Geometry, (b) scatter-free (c) single scatter and (d) random
noise

 

 

 

160

Figure 7 .2 continued

 

161

Figure 7.2 continued

 

 

 

 

 

(d)

7.4 Three dimensional (3D) model

A number of three dimensional MC models can perform accurate simulations of PET
systems, such as SimSET [68], Geant4 Application for Tomographic Emission
(GATE) [79]. The emission images were integrated into the simulation similar the 2D
case, which is generated by applying the motion models. The attenuation images are
also generated by assigning the appropriate attenuation coefficients in the 3D
phantom and incorporated in to the model. Finally the projection data, Sinogram or

list mode data can be output and stored. More details are provided in Chapter 8.

162

7 .5 System Data Acquisition

7.5.1 Sinogram

In order to understand motion correction in projection space, it is important to discuss
the projection or Sinogram data representation. A schematic of the overall procedure
for projection acquisition and sorting in Sinogram are depicted in Figure 7.3, in which

two groups of lines of response (LORs) at different projection angles are shown. A

PET imaging system is described by a linear model so that the measurement A, , the

photons detected at the jth detector pair, is represented by the integral along the
LOR of the product f (7c ) (radioactivity ) and SOC") (detector sensitivity). The
sensitivity function is determined by the scanner geometry and detector
characteristics. Mathematically, the detector measurement also called the ray sum is

represented simply in Equation 7.4.

163

 

Figure 7.3 Projection data acquisition and sorting in sinograrn: dashed lines and solid
lines represent the projections at an angle of -90 deg. and an arbitrary angle, where the
angle is between -90 deg. and 90 deg. Each projection reprents one complete view of the
object at one particular angle by lines of response (LORs)

The collection of II, with the same angle (9 is the projection P(R,-,t9), where
Ri is the trans-axial distance between each ray and the system geometric center.

Each point of the object f (7c) is included in only one ray of each projection, and
N

the sum over all rays for all the projections, 21?. J- is proportional to the total
i=1

radioactivity distribution of the 2D slice of the object. This entire collection of

164

projection data is defined as the Sinogram which is plotted in the 2D space spanned by

trans-axial distance R,- and projection angle0. The 2D projection is the 2D Radon

transform described as [81]:

P(R,0) = ] f(x, y)s(x, y)5(R — (xcosB + ysin 0))dxdy
FOV

(7.13)

where P(R,0) is a Radon transform (projection) of function p(x, y). R is the

coordinates of the LOR (ray) and 0 is the view angle or rotation angle. The

organization of P(R,6) , which is the record of coincidence events in PET by storing

R and 9, is called Sinogram.
7.5.2 List-Mode Acquisition

In the list-mode data acquisition, the pairs of X- and Y- position signals are stored in a
list, which is time dependent. The position signals can be obtained using position
circuits of the scintillation detectors from pulses in photomultiplier tubes (PMTs).
Gated signal indicators can also be included in the list-mode data. For example, if the
ECG signal is monitored in the gated cardiac PET imaging, the “trigger” signals are
also stored in the list. The acquisition of a gated cardiac image sequences is illustrated
in Figure 7.4. It should be noted, that gated acquisition is one of the three types of
frame-mode acquisition in nuclear medicine [3].

After the acquisition, the list mode can be reconstructed to obtain PET images for
display, or converted into projection data (Sinogram). The illustration of list-mode

acquisition is shown in Figure 7.4 with the gated timing signals also being stored.

165

The advantage of list-mode acquisition is that it is flexible for data manipulation and

image formation. The disadvantage is that it requires more disk storage space.

 

 

 

Memory

 

 

 

' 'Timing ' "
LORn LORn+liLORn+2J1 Marks LORn+3J LORn-I-4l \

Circuits-->

l"."""""

 

 

 

 

 

 

 

 

 

 

I
I
I
I
l
I
I
I
I
I
I
I
I
I
I
I
I
\\

 

 

 

LORn+5 3:125 LORD“; L0Rn+7L

 

 

Figure 7.4 Illustration of list-mode acquisition

7.6 PET Image Performance and Quality

There are several factors affecting the PET image quality besides the patient motion
and breathing noise. Figure 7.5 illustrates three difl‘erent coincidences in PET
imaging, which is also modeled in the 2-D geometry presented in this Chapter. Scatter
and random noise reduces the image contrast to noise ratio (CNR) as we can see in

the simulation results, are therefore one of main sources of PET imaging noise.

166

Figure 7.5 (a) True, (b) scatter and (c) random coincidences affecting the PET Image
Quality

 

167

Figure 7.5 continued

 

 

 

 

7.6.1 Scatter Noise

The scatter coincidence means one or both of the emitting photons from a single
annihilation are detected by ACD, which gives “misplaced” LOR The scatter noise

can be suppressed with energy resolution and collimation.

168

7.6.2 Random Noise

The random noise means emission photons from different annihilations are
simultaneously detected by ACD, which gives the “false” LOR. It can be suppressed

with small coincidence time and collimation.
7.7 Results and Discussion

7.7.1 2D results

We have developed a reliable 2D MC model in MATLAB environment for the study
of PET imaging with different noise levels [155]. In the 2D model the simulated

object is placed in a discretized region Q and surrounded by the detectors with

detector pair defined as (d 11,312) and detector ring radius r , detailed

configuration is shown in Figure 7.6.

169

 
 
    

Y axis (m)

..................

X axis (m)

Figure 7.6 2D Monte Carlo Model for patient motion study with discretization grids
illustrated

170

 

Y dimension (m)

 

X dimension (m)

(a)

Figure 7.7 (a) 2D phantom object and detectors, (b) discretized solution domain with
detection tube(red dotted line) and (c) attenuation coefficients map

 

Figure 7.8 Simulation results with one million trials: (a) true image (b) simulated
sonogram and (c) reconstructed image

The 2D phantom and the detectors are shown in Figure 7.7(a). The discretized region
Q has the dimension 600mm X 640mm and the radius of the PET scanner is 0.5m,

as illustrated in Figure 7.7(b). The simulated Sinogram and the reconstruction result

171

are shown in Figure 7.8(b) and Figure 7.8(c) respectively. Two types of noise, namely,
1) scatter noise and 2) random measurement noise were introduced in a controlled

manner. The noise level was defined by the parameter a as:

a : Nnoise/(N

noise

+Ntrue)

(7.14)

where N noise is the number of noise (scatter or measurement noise) events and NM“,

is the number of true annihilation events. The values of a = O, 0.2, 0.4 and 0.6 were
simulated and the results are shown in Figure 7.9(a) showing effect of noise due to
single scatter. Figure 7.10(a) shows the effect of measurement noise. A comparison of
the line plots in Figures 7.9(b) and 7.10(b) show the general similarity of true and
simulated signals. The difference for a = 0 case is due to the small number of

simulations used in this feasibility study.

172

 

Figure 7.9 Reconstructed images at different scatter noise level from (a) to (d):
reconstructed images at alpha = O, 0.2, 0.4 and 0.6

173

 

((1)

Figure 7.10 Reconstructed images at different random noise level from (a) to (d):
reconstructed images at alpha = 0, 0.2, 0.4 and 0.6

174

 

 

 

Modelling. , I - 7'

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

(MC Simulation) -
: 1
r
I

. ‘ ' '- ii'List-‘mode'fii‘53."
Sin ram ., , .2»::..-,:;
of . . ‘1 (Gating). : :jfgzrg :
Filtered Back Iterative EM
Projection ' ' 'Reeon‘ struction

 

 

 

L i
l

‘ Motion-free ? '
..lm‘age‘ i

 

 

 

 

Figure 7.11 Motion-free PET images modeling and image reconstructions procedure

The basic idea of PET modeling without motion is illustrated in Figure 7.11. The
“static” object is passed into the models and generated the projection data, either
Sinogram or list-mode data. Afier applying the proper reconstruction schemes, the
motion-free images are obtained. This approach is presented and validated in this

Chapter, which is very important before introducing any motion into this frame. An

175

updated schematic will be presented in Chapter 8, which incorporates the motion

detection and correction algorithms.

7.7.2 3D Results

Although only 2D sliced PET images and related Sinogram are analyzed in this
dissertation, a more realistic 3D human torso is constructed as shown in Figure 7.12.
This scenario has been modeled using the Simulation System for Emission
Tomography (SimSET) solver [68] developed by Division of Nuclear Medicine
University of Washington and originally published by R. L. Harrison [156].

For the 3D human torso model, the lungs are modeled as cylindrical objects with
inhalation and exhalation respiratory motion at a fi'equency of 12 cycles per minute;
the heart is modeled as spherical object with diastole and systole cardiac motion at a
frequency of 60 cycles per minute. The elliptical chest wall shape is also changing in
phase with the respiratory motion, but is set to a static cylindrical object here for the
sake of simplicity. Also this motion is seen to have little effect on the final Sinogram
data of internal organs. The spine is also simulated using a static cylindrical shape

without motion.

176

3D View

 

Figure 7. 12 3D view of the human torso model in the motion model, the chest wall,
lungs, heart and spine are simulated with 3D objects with varying shape at various
frequencies

7.8 Motion Models

The pulmonary and cardiac motion are intrinsically non-rigid deformable and difiicult
to be characterized using simple models. In the previous study for CT respiratory
motion, the non-linear synthetic motion functions were applied [142]. Because PET
images have less anatomical details than CT images, we apply a linear motion
function for approximating the changing anatomical shapes of internal organs, i.e.
lungs and heart. The motion patterns for the lungs are modeled linearly and a

quasi-linear function is used to model the motion of heart as illustrated in Figure 7.13.

177

Motion Signals

—e— Ling Motion

§ § +Cardiac Motion
9--------- - ----------:~ ---------- I ---------- -. -------------------- -

 

 

...t
O

 

 

 

 

 

 

Organ Radius(em): Lung(Cylinder), Heart(Sphere)
V

(a)

1
1O 1 5 20
Tirne(second)

D
01

Figure 7. 13 Motion signals for both pulmonary (circle & solid line) and cardiac (star &
solid line) motion

We generate twenty-one phases for each cycle of pulmonary and cardiac motion. In
the motion model, the lung radius R, varies fi'om 4cm (minima of the circle & solid

line) to 9cm (maxima of the circle & solid line) with a spatial step size of 0.5cm and

the temporal step size of 0.258. The heart outer wall radius RC varies from 3.6cm to

6.6cm with uneven spatial step size and 0.058 temporal step size (the star & solid line

in Figure 7.13). The cross section of the 3D human torso model is depicted in Figure

7.14, in which the chest wall size Ru and R0, remain the same throughout the

motion cycles.

178

 

Figure 7. 14 Cross section of the 3D human torso modeled as cylindrical and elliptical
objects with different sizes

In contrast to the simple x- and y- direction scaling method, a 2D algorithm
incorporating elastic deformation factors is introduced in this dissertation. The cardiac

motion, lung motion and chest motion are treated independently with different

contraction frequencies. RC , R1, (Ru , R0,) are the periodic radius of heart, lung
and human torso respectively and defined in equation (7.15) to (7.17):

R = Rcmx —RcMIN sin(27gfct +¢c) + Rem +RcMIN

 

 

 

 

c 2 2
(7.15)
R —R . R +R
R1: 1W2 Wsm(2nf,t+¢,)+ W2 [WV
(7.16)
R —R . R +R
Rt = W W sm(27?731+¢t)+ (MAX 2 tMN
(7.17)

179

where RM and RMHV are the maximum and minimum radius for the organs and

torso size. fc, f], f, are the frequencies for the cardiac motion, lung motion and
chest motion, in which the cardiac motion frequency is about 60 beats per minute,
and the latter two frequencies are assumed to be equal to the respiratory motion

frequency which is about 12 cycles per minute.

  
 

(Rrymx,
RtyMlN)
Cardiac Motion
Detection and .
Estimation '_' . . .
Rang!” ’43 '. l. .)
(R ) ' . (RcMAX,
“MIN ' 7 RcMIN)
let-'-
. 1 Cardiac Motion
j Detection and
Estimation

Figure 7.15 Illustration of motion encoded Sinogram with motion detection and
estimation

180

The motion-encoded sinograrn is illustrated in Figure 7.15, in which the motion
effects are exaggerated and shown as dotted region. The width of the motion-encoded

region is estimated and the motion information is extracted.

There are several approaches to evaluate motion in Sinogram space. In the R MAX and

RMHV estimation and organ boundary detection, the maximum and minimum radius

(amplitude of the motion) can be estimated according to the width of the motion

blurred region in the Sinogram.
7.9 Motion Effects on Sinogram

Mathematically, the Sinogram data for 2D projection are obtained using Equation

7.13. For the sake of simplicity, we assume the object function p(x, y) is a product

of f (x, y) and s(x, y). The following derivations are for 2D Sinogram only.

181

 

 

 

0%
“y
D

Figure 7.16 Mathematical parameters illustration for motion models

 

 

 

We start with one point at a spatial location (x, y) on the moving organ of interest
M a. After a period of time, it is moved to a new location whose coordinates are

represented by (x',y')=(x+Ax, y+Ay) where Ax and Ay are the motion in

L-R and A-P directions that need to be compensated. The corresponding sinograms

are represented as P(R,0) and P(R',0‘) at two different times with the
relationship as:
R = xcosB+ysin6l
(7.18)
R'= (x + Ax)cost9'+(y + Ay)sini9'= R + AR

(7.19)

182

By choosing 9 = 0' when sorting the Sinogram data in the data acquisition process, the
effect of motion reduces to simply scaling in trans-axial direction as:
AR 2 Axcos0+ Aysin0

(7.20)

For the points on the organs, the following relationship should also be satisfied:

 

AB=\/Kx .sz +K, ~Ay2
(7.21)

where K x and K y depend on the shapes of organs only. AR and AB are the

changes in spatial location and anatomical boundaries obtained using the fast CT
protocol. From Equations 7.20 and 7.21, we can calculate the motion in the two
dimensional space with the knowledge of motion-corrupted sinograrn and CT scans or

analytical Radon transform. The correction terms are expressed by Equation 7.22:
{Ax,Ay} = <D(AR,AB)
(7.22)

where AR, AR are estimated using simulated motion corrupted Sinogram and

motion models hypothesized in Chapter 8. (Dis the function incorporating the

mathematical operations presented in motion correction section.
7.10 Motion Synthesis Results

The torso phantom is discretized in 3D using brick elements in the solution domain

183

Q with a resolution of 2.5 X 2.5 X 4 mm in the lateral, anterior-posterior and
superior-inferior directions respectively.
Discretizing a complete respiratory cycle into twenty-one temporal phases, six phases

of the lung motion in the first half cycle are shown in Figure 7.17 (exhalation minima

to inhalation maxima). The radius of lungs R, varies from 40 mm to 90 mm in 10

mm steps. The radius of the heart outer wall RC is 40 mm with the vessel radius set

to 20 mm. The chest wall radius is set to be 240 mm (L-R) and 160 mm (A-P). The
spine cylinder is 20 mm in radius.

The material properties and attenuation coefficients are assigned appropriately in each
object and surrounding air. The amount of radioactivity in each object and air are also
defined so that the dose concentration in heart, bone, body, blood, lungs and air is in

descending order. The trans-axial and azimuthal directions in the sinogram are

grouped into 288 bins (720 mm FOV) and 360 bins (half circle FOV, 5° -bin‘l ). The

number of bins determines the resolution in the sinogram space and hence the

reconstruction resolution. The acceptance angle for the PET detector is set to be 5°

and minimum acceptance photon energy is 53 keV for a 511 keV emitted photon pair.

184

 

Figure 7. 17 Six respiratory phases with different lungs sizes for input phantom images:
exhalation minima (upper left) and inhalation maxima (lower right)

Figure 7.18 Sliced simulation results from 3D model: (a)-(b) sinogram and (c)
corresponding reconstructed images using FBP

 

Figure 7.18 continued

 

(C)
The reconstructed PET images using FBP are shown in Figure 7.18(b) with the

corresponding sinogram data in Figure 7.18(a). The radioactivity concentration in the

186

lungs is artificially set to be high for the purpose of enhancing the motion effects in
sinogram data. Attenuation correction has not been applied since only geometry
dependent motion is of interest. It is clearly seen that pulmonary motion changes the
pattern of sinogram intensity distribution derived in Section 7.9 in the trans-axial
direction (R-axis). The torso-air boundaries in the sinogram remain the same since the
chest wall was assumed to be static under the simplification hypothesis. The motion
corrupted image is obtained by a weighted sum of all the sinograms at different

phases for multiple cycles.

187

CHAPTER 8. INVERSE PROBLEMS IN POSITRON

EMISSION TOMOGRAPHY IMAGING

8.1 Introduction

Patient motion, including body movement, head movement and thoracic motion due
to respiratory and cardiac cycles, is a major source of artifacts in nuclear medicine
imaging [124][l25][126][129][131][134][135]. In the past few years extensive
research has been conducted for the purpose of motion correction for PET imaging
[83][84][85][87][88][89][91][101][102][104][106][lO7][108][110][1l6][127][133]
[138][139][l40][146][147].

We propose a new approach to correct for thoracic motion in PET images using
sinogram data in the projection space. This approach involves detecting motion
components in the projection space and correcting them before image reconstruction
is performed. This sino gram based method employs the entire projection data set, in
contrast to the multiple acquisition frames (MAP) [127] and gating [101][116]
methods, which involve time averaging for increasing the SNR. Once motion
compensation is performed in sinogram space, standard reconstruction algorithms
[15][143][l44] can be applied to obtain a reconstructed image flee from motion
effects.

The methodology presented in this dissertation extends the motion correction in

188

projection space to 3D where the elastic organ deformation and thoracic motion are
simulated using motion models presented in Chapter 7. The objective of this work is
to not only compensate for motion in PET imaging but also retain the high SNR,
which will benefit the subsequent diagnosis of lesions/tumors, treatment planning and
therapy. For simplicity, in the later sections, we use the term “motion” to stand for
“thoracic motion” unless we mention the specific motion type, such as respiratory,

cardiac, body, or head motion.
8.2 PET imaging with motion

Thoracic motion effects on PET imaging data can be traced back to two types of
sources: the cardiac cycle, namely, motion of heart at about 1 second per cycle and
respiratory cycle, namely, motion of the organs in the abdominal region including
lungs, stomach, liver and spleen. The period for the respiratory motion is about 5
seconds per cycle. Patients usually lie in the PET gantry for a few minutes, which is
much longer than the cycle period. Consequently, the acquired PET data usually
containing tens of millions of events are degraded by the thoracic motion. The
“motion corrupted” PET image data leads to errors in diagnosis, treatment planning,
dose estimation and subsequent therapy. In the worst case scenario a small tumor
could go undetected due to blurring in the motion-corrupted images. Failure of early
detection and diagnosis can in turn decrease the probability of therapy success. Also
motion effects can cause error in estimating the stage of a tumor since the size of

lesions and tumors are wrongly reconstructed without motion correction. Significant

189

volume changes of lung lesions were observed by Nehrnen et al. leading to
misdiagnosis [101]. In the fusion of PET and CT images investigated by Osman et al.,
mis-localizations of lesions and tumors were found to be due to motion [116]. The
impact of motion can be reduced in several ways as reported in previous studies. In
the case of respiratory motion there are several ways to diminish motion artifacts as

discussed in Section 8.3. Typical motion artifacts are shown in Figure 8.1 [89].

 

Figure 8.1 Motion artifacts in PET in comparison to CT images: (a) coronal view and (b)
sagittal view

190

8.3 Current motion correction approaches

In the case of respiratory motion several methods to diminish motion artifacts have
been reported in the previous studies with their associated advantages and limitations.

Multiple acquisition fi'ame (MAF) methods proposed by Picard and Thompson
involve regrouping the projections into multiple subsets after detecting motion [127].
However the MAF results can be noisy due to the small numbers of projections in
each frame, used for reconstruction. To overcome this problem, higher radiation dose
might be needed. Furthermore, normal clinical setup does not allow acquisition of
MAF data and motion correction in a retrospective way. One common approach used
in clinical work is to acquire data in each position, multiple times and discard the
apparent motion distorted images using visual examination by radiologists or
technicians.

Post-processing methods involve applying transformations and image registrations on
reconstructed images. These methods were typically proposed to improve the image
fusion between different modalities, such as PET/CT, PET/MRI and the latest
performance of these approaches for PET images can be found in [104][150]
published by Berger and McGough. Rigid and affine transformations results were
presented by Camara and Mattes [148][149]. More systematic introduction of nuclear
medicine image registration in post-processing stage is found in Maintz [151] and
Hill’s review articles [152]. Post-processing methods in the image space, however,
must be off-line and constrained by the time flames of reconstruction and registration

algorithms. This procedure is time consuming due to the nonlinear, elastic nature of

191

organ and thoracic cavity deformation. The intrinsically distorted projection data used
in subsequent reconstruction also render compensation techniques in post-processing
stage unreliable.

Breath-hold methods which involve holding breath during data acquisition with fast
imaging techniques [145][146] and new scanning protocols is an on-going research
area and the latest techniques were just published by Nehmeh [147]. However, hold
the breathing for long times can make patients uncomfortable; especially those
subjects having lung lesions and/or respiratory complications.

Gating methods performs a respiratory-gated acquisition employing external signals
dependent upon the state of respiration, e.g. abdominal belt, skin marker, etc.
proposed by Nehmen et al. [101] and Visvikis et al. [116]. For cardiac motion, an
external signal, e. g. the electrocardiagram (ECG) signal can be obtained
simultaneously to perform a gated acquisition. The PET data for a heart beat cycle is
stored in gated frames to eliminate the motion artifacts. Improvements in reducing the
blurring due to motion using gating methods have been demonstrated. Gating with
breathing cycle by tracking external source points situated over the patients and
respiration correlated dynamic PET have been developed recently. Apparently, gating
methods require a priori information on respiratory and cardiac cycles, which implies
need for additional hardware and complex implementation and synchronization. A
posteriori gating by estimating the respiratory frequency from acquired PET dynamic
data was published by Visvikis [117]. By using the low number of events detected by
ACD in a phase-of-interest in the entire cycle for image reconstruction, the image

quality suffers from the low statistics like MAF as well.

192

Model-based motion compensation methods that incorporate motion models directly
into the image reconstruction were studied by Jacobson and Fessler [84], Qiao [112],
and Li [153] recently. Progresses of these techniques are published by Reyes [107].
The main problem in these approaches is the need for a 4D CT scan for constructing
the motion model at each phase of the motion cycle. Further, these techniques are
associated with the reconstruction algorithms rather than motion compensation in the
projection space.

Recent research on the 5-dimensional (5D) models that incorporate two more
independent variables, i.e. tidal volume and breathing air flow besides the three
spatial dimensions are developed at Washington University in St. Louis School of
Medicine. This approach has been validated for computed tomography (CT) studies,
also referred to as 5-D CT [154] where motion correction for PET imaging research is

undergoing.

8.4 Motion Correction in Projection Space

The overview of the motion correction scheme is shown in Figure 8.2. The methods
proposed in this dissertation are depicted in the center box of this scheme figure. The
filtered back projection is implemented for sinogram projection data correction. The
iterative expectation maximization (EM) reconstruction methods will be implemented
for further list-mode projection data correction (See Chapter 9 Conclusions and

Future work).

193

Model-based Motion
Correction Approach

 
   

OI" I

Figure 8.2 Overview of motion correction scheme

To overcome the limitations of existing motion correction techniques, this dissertation
proposes an innovative projection space correction approach prior to reconstruction.
A schematic of the proposed approach contrasted with image space correction

approach is illustrated in Figure 8.3.

194

PET Scans

r:>

Patient wl
motion

  

 

Prior correction in
projection space:
SINOGRAM

 

 

 

 

/ \

Raw Data w! motion

Projection Space:
{SINOGRAM LIST-
MODE")

1

Gotten corrected\

Raw Data

 

 

/

 

Projection Space:
{SINOGRAM, LIST-

 

Before image reconstruction
(no intrinsic error)

 

MODE,...}
\ /

    
   
     
   

Reconstruction
Motion
"HQ Impaired
? Images

5

Posterior
correction in
Image space

E

Corrected
PET Images

 

 

 

 

‘

Figure 8.3 Schematic of projection and image space motion correction in PET imaging

Compensating for motion in the raw data/proj ection space (sinogram or list mode),

has several advantages:

1) Implementing the processing algorithms in PET scanner

hardware can make this approach fast and potentially real time in contrast to image

space registration and processing; 2) Offers higher accuracy since intrinsic errors

introduced in image reconstruction is eliminated; 3) In contrast to the gating

techniques, correction in projection space retains the full statistics and thus improves

the image resolution and target to background (T/B) ratio levels. Particularly, cancer

imaging, tumor detection or delineation applications will be enhanced using this

approach.

195

8.5 Motion Tracking, Estimation and Extraction

To understand the motion better, we introduce some fiducial points with 3D radius

equals to 5mm X 5mm X 4mm that are attached to the lung-blood interface. In order to

track the motion in both L-R and A-P directions, three fiducial points fa , fb , and
fc are assigned as illustrated in Figure 8.4, in which the distance from the phantom

center to the lung center is D, and the lung radius is R, .

 

 

 

 

 

 

Figure 8.4 Fiducial points allocation for motion tracking

Since in most practical reconstructions, 3D PET volumetric data also needs to be
approximately rebinned into a 2D format, it is reasonable to assume that no fiducial
point is needed in the S-I direction (i.e. for out-of-plane motion). The coordinates

(L-R, A-P, S-I) for the fiducial points used in this dissertation are obtained as follows:

196

f(1 = (D1 + R: ,0.0)

(8.1)
fb = (—D,,R,,O)
(8.2)
fc = (—D, — R, cos45°,—R, sin 45°,0)
(8.3)

The simulation result with fiducial points {fa, fb,fc} is shown in Figure 8.5, in
which the radioactivity at these points are set to be relatively high. The tracking of
lung motion in the sinogram shown in Figure 7.16 can be simplified to tracking of the
high-intensity pixels associated with the fiducial points. The motion pattern is

estimated and extracted such that AR, AF in Equation 7.19 can be predicted.

Reconstructed PETimage using FBP

 

Figure 8.5 Reconstructed image with fiducial points for motion tracking and estimation

197

Motion extraction in sinogram is based on prior information providing the starting
point for searching motion, as illustrated in Figure 8.6. The analytical Radon
transform is used to generate the sinogram of just fiducial points using Equation 7.13.
Fast CT acquisition protocol, namely, the breath-hold helical CT scans that offer rich
anatomical information which provides “landmarks” that serve the same purpose as
fiducial points used in these studies. Figure 8.6(a) is the simulated prediction for the
three fiducial points in one particular breathing phase that is obtained using fast
Radon projection method. The correlated PET sinogram data is shown in Figure
8.6(b), in which lines of strong intensity indicate the locations of the fiducial points. It
is clearly demonstrated that the prior information obtained using CT can serve as a
reasonable starting point for motion tracking. We adopted fast Radon transform
simulation in our study to mimic the fast CT acquisitions. By incorporating more
motion phases and information, boundaries of the organ Of interest can be estimated
using fast radon transform simulator. Figure 8.7 shows the effect of motion on the
sonogram; Figure 8.7(a) is obtained by considering just two phases, i.e. extreme
inhalation and exhalation time instances and Figure 8.7 (b) is obtained with
twenty-one phases characterizing the full respiratory cycle. The intensity distribution
of the sinogram not only reveals the boundaries of organs but also provides the
information for further compensation in the same projection space. The motion
blurring effect is clearly seen in Figure 8.7(b) at the outer boundaries of the bands,

which is due to the motion at the organ-blood interface.

198

 

(a) (b)

Figure 8.6 Motion extraction in sinogram based on prior information obtained using fast
Radon transform result in (a) and high intensity extraction in simulated sinogram (b)

199

Azimuthal Angle

 

Trans-axial distance

Figure 8.7 Motion effect estimation in sinogram for right lung PET image in noise-free
case

Once the organ of interest is defined, e.g. the right lung in Figure 8.7, the radon

transform simulator and MC simulator results are combined using maximum
correlation method to provide motion corrupted organ boundary M , reference organ

boundary S and the center for the ROI RC. To simplify the mathematical

A

equations,M , S , and Re are all integer numbers and defined as fraction of the
total length trans-axial distance vector (number of data points in the R direction of

sinogram projection space).

200

For different azimuthal angles in the sinogram, the estimated boundaries and center

coordinates are different. In the following derivations, the angle dependent term is

omitted but actually the estimated variables M , 5' , and RC are all ftmctions of
- -M95(o) i119 9:- --

azrmuthal angle, i.e. ( ), and c( ). For an angle , the mtensrty

plots in the motion corrupted and reference projection are illustrated by the

dashed-line and solid-line respectively in Figure 8.8.

Sinogram Intensity

 

 

 

h
. . A . R-direction Pixel Number
-M —s S M

- — — Motion-corrupted

 

Motion-free/reference

Figure 8.8 Illustration of the line plots of motion-corrupted and motion-free sinogram
intensity with boundaries estimation labeled

8.6 Motion Correction Implementation and Results

Let the reference sinogram (target), motion-corrupted sinogram and motion-corrected

201

sinogram be represented as P(R,9) , R(R,6) and R(R,6) respectively. The
resolutions in the trans-axial distance direction (R direction of sinogram projection
space) are denoted as AR and AR for reference and motion-corrupted sinogram.

The width of ROI in the above two cases are easily obtained when AB is estimated

by equations 8.4 and 8.5:

W=2AR§

(8.4)

NA

(8.5)

Since M and S are integers representing the motion-corrupted and reference ROI

boundaries, we can label each data point in projection space fi'om center to the
boundaries as j=l,2,...,M—l and k =1,2,...,§—l. Since we discretized
the range of azimuthal angle (projection angel) (- 90° to 90°) into N A bins, we
can represent any pixel intensity in the sinogram as P91. (R j) or P91, (Rk) , where
1': l,2,...,NA.

The motion correction algorithm can be expressed in a closed form function shown in

Equation 8.6 for each projection angle 6,- :

W

P9i(R$‘—k)=P9i(W_W 'Rj)+AP0,-,M-j

(8.6)

202

Substitute Equation 8.4 and 8.5 into Equation 8.7 and we get (derived in Section 8.8).

1 N—1~
F 12139, (RM-j+1)
=0

 

AR-S‘ ~

13.R~ =F . .
0,( 5-1,) 0<AR°°M~ARS

l

(8.7)

112—1

A

 

where j =

 

. kl and H is the rounding operator to get closest integers. The

first term involves a scaling operation and mapping from the corrupted sinogram to
1 N—l~

the target one. The averaging term X,- ngi (R M_ j + ,) is applied to make sure the

corrected sinogram intensity distribution is smooth. Additional details and pseudo

codes for the correction algorithm are attached in the Appendix.

Typical sinogram images before and after motion correction are shown in Figure

8.9(a) and 8.9(b), respectively. The sinogram intensity is smeared with wider ROI as

shown in figure 8.9(a) which is due to the cumulative motion effects in the respiratory

cycles. The corresponding reconstructed images are shown in Figure 8.10.

203

 

(a) (b)

Figure 8.9 Sinogram before and after the motion correction. The cumulative motion
effects with widen ROI can be seen in (a) due to five cycles of respiration simulated

 

Figure 8.10 Reconstructed PET images: (a) target image at one particular respiratory
phase, (b) motion blurred image and (c) corrected image using the scaling-mapping
algorithm

204

The reference/target image is shown in Figure 8.10(a). Figure 8.10(b) and 8.10(c) are
motion-corrupted image after five cycles of respiration and motion-corrected image
using the scaling-mapping algorithm. The image blurring effect is clearly
demonstrated by the model. It should also be mentioned that additional procedures
were conducted in the image space, such as image de-blurring, non-linear image
registration, etc for the purpose of better anatomical delineation and multi-modality
image fusion. The scaling-mapping method presented here not only performs
innovative motion tracking, estimation and extraction, but also overcomes difficulties
in the image space manipulation and preserves clear organ boundaries by reducing the

motion contamination in the raw projection space.
8.7 Discussion and Conclusion

There are several issues and practical constraints of the proposed algorithm that
deserve to be mentioned:

1) Standardized uptake value (SUV) Accuracy

The SUV accuracy is achieved by the conservation law, i.e. keeping the sorted photon
pairs number constant, for the radioactive event counts in the sinogram. The SUV
based on patient’s body weight is the most common parameter extracted fi'om nuclear
medicine images for quantitative evaluation, which is related to the injected dose per

body mass and defined in Equation 8.8.

C (0

SUV =
W InjectedDose/PatientsWeight

 

205

(8.8)
where C (t) is the tissue radioactivity concentration in kBq/ml at time t. The patient
motion introduces the artifacts that decrease the tissue radioactivity concentration
term which can be demonstrated by the spreading effect of the corrupted sinogram in
Figure 8.9(a). Using the prior-correction method and the conservation law by
maintaining the total radioactive event counts in one scanning session, the SUV is the
same or is very close to that of motion-free images taken with patient holding breath.
This is one advantage for the correction in projection space because it is difficult to
manipulate the image pixel intensity in image space to get the accurate estimate of
SUV.

2) Residual Artifacts

Some residual artifacts can still be seen in the projection space corrected image in
Figure 8.10(c) on the reconstructed organ boundaries and the low intensity region in
the image center. A possible reason for these artifacts might be due to the Filtered
Back Projection (FBP) method used. A more accurate but time consuming iterative
reconstruction algorithm such as Expectation Maximization (EM) iterative
reconstructions need to be investigated to reduce the residual artifacts.

3) 3D Volumetric Data

In real applications we encounter 3D volumetric data. The 3D PET acquisition with
the inter-slice septa removed greatly increases the scanner sensitivity and image
signal to noise ratio (SNR). The patient motion is inherently in three dimensions and

the 3D radioactive isotope distribution results in the 3D projection data. However,

206

since the computational load for 3D data processing is much greater than that of 2D,
this requires use of rebinning. The method for converting the projections from 3D
I data into 2D approximations is widely used and one sinogram per transverse slice is
obtained. The out-of-plane motion tracking and estimation is simplified in this
dissertation by in-plane motion using three fiducial points tracking. The fully 3D
iterative approaches are also available and the out-of-plane motion correction is

expected to be taken care of in the list-mode space.
8.8 Scaling-mapping method pseudo codes

The derivation of motion correction algorithm and its pseudo codes are attached here.
LOOP A: FOR 1' =1,2,...,NA
LOOP B: FOR j =1,2,...,M —1
Intensity difference at jth R direction pixel in motion-corrupted sinogram

(updating step calculated in Equation 8.9)

~ 1 N-1~
APBi M-j : P9: (Rf) _ X," £100: (RM-141)

(8.9)
LOOP c: for k =1,2,...,§ —1
Updating target sinogram by linear mapping pixels from

motion-corrupted space to reference space by equation (8.10):

207

 

 

 

1' = fl ' k
(8.10)
Corrected projection data is calculated as equation 8.11 as:
1’59. (R§_k)=f59. (TE—«jump, r4--
1 r W _ W 1, J
(8.11)

END LOOP C (Trans-axial correction)

 

CONTINUE (Forced to end Loop B)
END LOOP B
END LOOP A (Azimuthal correction)

Equation 8.11 can also be rewritten using all the known quantities in a closed-form

 

as:
P R —F ———W R AP
0i( S-k)— 0i(W—W. J)+ 95,M—j
N [512.3 N 1 N-1~
=Poi(AR.M_AR.§'Rj)+Po,-(Rj)—figpoi(RM-j+1)
(8.12)
where j = 5.53-:11. kl and H is the rounding operator to get closest integers.

 

208

CHAPTER 9. CONCLUSIONS AND FUTURE WORK

9.1 Accomplishments and Conclusions

(a) Formulated a general mathematical model-based approach for noninvasive
imaging in NDE and biomedical application

(b) Implemented qualitative and quantitative image analysis procedures

(c) Developed computational models for simulating M01 and GMR NDE techniques.
Used these models to conduct a parametric study for optimization of sensors and
systems

((1) Results of parametric study and data classification algorithms were used in
estimating Probability of Detection (POD) curves for the NDE techniques

(e) Developed computational models based on Monte Carlo (MC) technique for
simulating PET imaging

(1) Developed and incorporated motion models into the MC-PET model

(g) Developed motion compensation techniques for correction in sinogram space

9.2 Future Work

In this dissertation, a new patient motion compensation method in projection
space-the sinogram is proposed and initial results validate this approach using Monte
Carlo simulation that incorporates respiratory motion. Innovative motion tracking,

estimation, and extraction is presented with the knowledge of prior anatomical

209

information and forward Radon transform. The respiratory motion artifacts are
significantly reduced by correcting the motion of each 2D sliced image from the
entire 3D volumetric PET data after rebinning procedure. Future work is needed to
improve the reconstructed image quality by using the iterative reconstruction
algorithms with huge size list-mode data. More realistic motion models need to be
tested so the robustness of this approach can be evaluated.

Besides the organ boundary extraction and detection based correction, there are other
approaches need to be investiaged: a) Sinogram stratification and motion frequency
estimation. The sinogram is possibly stratified into different layers based on the
intensity. Due to the higher intensity in the sinogram data means more appearances of
the organ (dose distribution). The motion frequency is able to be estimated based on
the intensity distribution of the blurred regions; b) Gradient-based method for motion
vector estimation. In this method, the motion vectors also need to be estimated using
gradient-based method if the stratification of sinogram is correct and accurate; c)
List-mode based motion correction methods. Since projection data is stored as
list-mode, which is large in size, efficient methods need to be proposed to process and
correct the motion in list-mode. The list mode data can be simply treated as “time

tagged” projection data, so it can also be converted into sinogram and then corrected.

210

“Motion”

Breathing Signal
“Motion- (oracle/Lung Volume)

free”

I‘_’l‘ ’I‘—’I‘ ‘ “““

Rest Inhalation Rest Inhalation
& &

Expiration Expiration

Figure 9.1 List-mode correction using rest motion phase

When the list-mode data is acquired with some external hardware, which can output
the breathing signal, it can be corrected by only averaging the “motion—free” part of
data, as we can see in Figure 8.11. The rest duration is about one half of the whole
respiratory cycle, and the “motion-flee” data can be averaged after the start and end
time of this rest periods are detected. The inhalation and expiration phase data of the
breathing cycle needs to be discarded.

The other future work includes:

(a) BM imaging sensors and system optimization through POD studies.

(b) More quantitative feature extraction for EM imaging data.

(c) More accurate motion models development to describe the patient motion not

limited in the lungs and heart, but extended to upper abdominal regions.

211

(d) Quantitative comparison of motion correction methods between sinogram
approach and current correction results through commercial software and real
patient data.

(e) Model-based attenuation correction for PET/CT after motion correction.

212

REFERENCES

[1] A. Webb, Introduction to Biomedical Imaging, Wiley-Interscience, 2003.

[2] R. C. Gonzalez, and R. E. Woods, Digital Image Processing, 2nd Ed., Prentice Hall,
2002.

[3] J. T. Bushberg et al., The Essential Physics of Medical Imaging, 2nd Ed., Lippincott
Williams & Wilkins, 2002.

[4] M. Bertero et al., Introduction to Inverse Problems in imaging, Institute of Physics
Publishing, Bristol and Philadelphia, 1998.

[5] C. T. A. Johnk, Engineering Electromagnetic Fields and Waves, 2nd Ed., John Wiley
& Sons, 1976.

[6] N. N. Rao, Elements of Engineering Electromagnetics, Prentice Hall, 2000.
[7] http://www.j1ab.org
[8] http://en.wikipedia.org/wiki/Positron_emission_tomography

[9] J. B. West, Respiratory physiology. The Essentials, 6th Ed., Lippincott Williams &
Wilkins, 2000.

[10] J. Thibaud, L’ annihilation des positrons au contact de la matiere et la radiation qiu en
resulte, C R Acad Sci, 197(1933), 1629-1632.

[11] F. J oliot, Preuve experimentale de l’annihilation des electons positives, C R Acad Sci,
197(1933), 1622-1625.

[12] R. Beringer, et al., Phys. Rev. 61(1942), 222.

[13] J. Mitchell, Applications of recent advances in nuclear physics to medicine, Brit. J.
Radiol. , 19(1946), 481-497.

[14] C. Tobias et al., The elimination of carbon monoxide fi'om the body with refernce to
the possible conversion of CO to CO2, Amer: J. Physiol., 145(1945), 253-263.

[15] P. E. Valk, Positron Emission Tomography: Basic Science and Clinical Practice, P. E.
Valk, D. L. Bailey, et al. (eds.), Springer, 2003.

[16] N. A. Dyson et al., The preparation and use of Oxygen-15 with particular reference
to its value in the study of pulmonary malfunction, 1958 Peaceful Uses of Atomic Energy,
Geneva, 1958, 103-115.

[17] C. T. Dollery, Metabolism of oxygen—1 5, Nature, 189(1960), 1121.

213

[18] H. Anger, Scintillation camera and positron camera in: Medical Radioisotope
Scanning, Vienna: IAEA, (1959), 59-73.

[19] S. Rankowitz et al., Positron scanner for locating brain tumors, IRE International
Convention Record, 10(1962), 49.

[20] Y. L. Yamarnoto et al., Study of quantitative assessment of section microregional
blood flow in man by multipole positron detecting system using Krypton-79, BNL Med.
Dept. Circ. Brookhaven, Report no. 28, 1966.

[21] D. E. Kuhl et al., Image separation radioisotope scanning, Radiology, 80(1963),
653-662.

[22] D. E. Kuhl et al., Quantitative section scanning using orthogonal tangent correction,
J. Nucl. Med, 14(1973), 196-200.

[23] D. A. Chesler et al., Transverse section imaging of myocardium with 13NH4, J. Nucl.
Med. 14(1973), 623.

[24] M. M. Ter-Pogossian et al., A positron emissin transaxial tomograph for nuclear
imaging, Radiology, 114(1975), 89-98.

[25] M. E. Phelps et al., Application of annihilation coincidence detection to transaxial
reconstruction tomography, J. Nucl. Med, 16(1975), 210-224.

[26] N. A. Dyson, The annihilation coincidence method of localizing positron emitting
isotopes and a comparison with parallel counting, Phy. Med. Biol, 4(1960), 376-390.

[27] GL. Fitzpatrick, Flaw Imaging in Ferrous and Nonferrous Materials Using
Magneto-Optic Visualization, US. Patents 4,625,167, 4,755,752, 5,053,704 and
associated foreign patents.

[28] GL. Fitzpatrick et al., Magneto-optic/Eddy Current Imaging of Aging Aircraft: a
New NDI Technique, Material evaluation, 1993, 1402-1407.

[29] GL. Fitzpatrick et al., Magneto-Optic/Eddy Current Imaging of Subsurface
Corrosion and Fatigue Cracks in Aging Aircraft, Review of Progress in Quantitative
Nondestructive Evaluation, 15(1996).

[30] K. Zvezdin and V. A. Kotov, Modern magneto-optics and magneto-optical
materials, Institute of Physical Publishing, Bristol and Philadelphia, 1997.

[31] P. Rarnuhalli et al., Enhancement of magneto-optic images, The 9th International
Workshop on Electromagnetic Nondestructive Evaluation, Saclay, France, 2003, 199-205.

[32] U. Park, L. Udpa, G C. Stockmart, Motion-based Filtering of Magneto-Optic
Irnagers, Image and Vision Computing, 2004, 243-249.

214

[33] U. Park, L. Udpa, G C. Stockrnan, B. Shih, J. Fitzpatrick, Real-time Implementation
of Motion-based Filtering in Magneto-Optic Irnager, 30th Annual Review of Progress in
Quantitative Nondestructive Evaluation, Green bay, Wisconsin, 2003.

[34] L. G Shapiro, G C. Stockrnan, Computer Vision, Prentice Hall, 2001.

[3 5] Y. Fan et al., Automated rivet inspection for aging aircraft with magneto-optic
imager, Electromagnetic Nondestructive Evaluation (IX), IOS Press, 2005, 185-194.

[36] Y. Deng et al., Automatic classification of M01 for aircraft rivet inspection,
Proceedings of 12th International Symposium on Applied Electromagnetics and
Mechanics, Bad Gastein, Austria, 2005.

[37] S. Udpa, L. Udpa, “Eddy Current Nondestructive Evaluation,” Encyclopedia of
Electrical and Electronics Engineering, John G Webster, Editor, John Wiley, 1999

[3 8] S. Udpa, Nondestructive Testing Handbook, Electromagnetic Testing, American
Society for Nondestructive Testing, Technical Editor, Columbus, 2004

[39] J. Salon and J. Schneider, “A Comparison of Boundary Integral and Finite Element
Formulations of the Eddy Current Problems,” IEEE Transactions on Power Systems,
100(1981), 1473-1479

[40] N. Ida, Numerical Modeling for Electromagnetic Non-Destructive Evaluation
Chapman & Hall, McGraw Hill, 1984

[41] O. Riibenkonig, The Finite Difference Method (F DM) - An Introduction, Albert
Ludwigs University of Freiburg, 2006

[42] Y. K. Na and M. A. Franklin, “Detection of subsurface flaws in metals with GMR
sensors,” Review Of Progress in Quantitative Nondestructive Evaluation, 760(2005),
1600-1607.

[43] L. Xuan, “Finite Element and Meshless Methods in NDT Applications,” Ph. D
Dissertation, Iowa State University, 2002

[44] Paul Rutherford, Boeing Phantom Works, Seattle, WA
[45] J. Jin, The Finite Element Method in Electromagnetics, New Your: Vlfrley, 2002

[46] Y. Tian et al., Simulation of the world federation's first eddy current benchmark
problem," Review Of Progress in Quantitative Nondestructive Evaluation, 700(2004),
1560-1566

[47] Z. Zeng et al., “Probability of Detection Model for Gas Transmission Pipeline
Inspection”, Research in Nondestructive Evaluation 3(2004).

[48] L. Xuan et al., “Developments of a Meshless Finite Element Model for NDE

215

Applications,” Proceedings of the 29th Annual Review of Progress in Quantitative
Nondestructive Evaluation, Bellingharn, WA, 2002

[49] Z. Zeng et al., "Modeling microwave NDE using the element-free Galerkin method,"
Electromagnetic Nondestructive Evaluation (120, 108 Process, 2005, 41-48

[50] L. N. Trefethen and D. Bau, Numerical Linear Algebra, SIAM, 1997
[51] T. Dogaru and S. T. Smith, IEEE. Trans. Magn. 37(2001), 3831-3838.

[52] T. Dogaru et al., Deep crack detection around fastener holes in airplane
multi-layered structures using GMR-based eddy current probes,” in Review of Progress in
QNDE, NY 700(2004), 398-405.

[53] A. J ander, C. Smith and R. Schneide, Magnetoresistive Sensors for Nondestructive
Evaluation, in 10th SPIE International Symposium, 2005.

[54] N. Nair, V. Melapudi, H. Jimenez, X. Liu, Y. Deng, Z. Zeng, L. Udpa, T. Moran and
S. Udpa, IEEE. Trans. Magn. 42(2006), 3312-3314.

[55] L. Xuan et a1, “Finite Element Model for MOI Applications Using A-v
Formulation,” Review of Progress in Quantitative Nondestructive Evaluation, 557(2001),
385-391.

[56] G L. Fitzpatrick, "hnaging Near Surface Flaws in Ferromagnetic Materials Using
Magneto-Optic Detectors," Review of Progress in Quantitative Nondestructive
Evaluation, 1985.

[57] R. Albanese et al, A Comparative Study of Finite Element Models for Magneto Optic
Imaging Applications, 6th International Workshop on Electromagnetic Nondestructive
Evaluation, Budapest, Hungary, 2000.

[58] Z. Zeng et al., Application of Taguchi Methodology for Optimizing Test Parameters
in Magneto-Optic Imaging, Proceedings of the 29th Annual Review of Progress in
Quantitative Nondestructive Evaluation, 22(2003), 651-658.

[59] B. E. C. Koltenbah, “Demodulation Algorithm for the Boeing Magnetoresistive
Sensors Program”, Boeing Physical Sciences Research Center.

[60] T. Belytschko, Y. Krongauz, D. Organ, M. Fleming, and P. Krysl, “Meshless
methods: an overview and recent developments,” Comput. Methods Appl. Mech. and
Engn, 139(1996), 3-47.

[61] J. Dolbow and T. Belytschko, “An introduction to programming the meshless
element fiee Galerkin method,” Archives in Computational Mechanics, 5(1998), 207-241.

[62] L. Xuan, Z. Zeng, B. Shanker, and L. Udpa, Element-free Galerkin method for static
and quasi-static electromagnetic field computation, IEEE Trans Magn., 40(2004), 12-20.

216

[63] S. C. Strother, et al., Measuring PET scanner sensitivity: relating countrates to image
signal-to-noise ratios using noise equivalent counts, IEEE Trans Nucl. Sci, 37(1990),
783-788.

[64] CW. Williams, et al., Design and performance characteristics of a positron emission
computed axial tomography-ECAT-II, IEEE Trans Nucl. Sci, 1979, 619-627.

[65] L. Lupton and N. Keller, Performance Study of Single Slice Positron Emission
Tomography Scanners by Monte Carlo Techniques, IEEE Trans on Med. Imag. 2(1983),
154-167.

[66] H. Zaidi, Relevance of accurate Monte Carlo modeling in nuclear medical imaging,
Med. Phys. 26(1999), 574-608.

[67] http://www.bme.unc.edu/~wsegars/NCAT_instructions.htm

[68] http://depts.washington.edu/simset/html/user_guide/user_guide_index.htrnl

[69] C. J. Thompson, J. Moreno-Cantu and Y. Picard, PETSIM: Monte Carlo simulation
of all sensitivity and resolution parameters of cylindrical positron imaging systems, Phys.

Med. Biol, 37(1992), 73 l-749.

[70] M. Ljungberg, et al., Monte Carlo Calculations in Nuclear Medicine: Applications in
diagnostic imaging, Institute of Physics Publishing, Bristol and Philadelphia, 1998.

[71] H. Kahn, Use of different Monte Carlo sampling techniques, Monte Carlo Methods,
H. A. Meyer (edited), Wiley, New York, 1956.

[72] D. E. Raeside, Monte Carlo principles and applications, Phys. Med. Biol, 21(1976),
181-197.

[73] R. L. Morin, Monte Carlo Simulation in the radiological sciences, CRC, Boca Raton,
FL, 1988.

[74] I. Lux, et al., Monte Carlo Particle Transport Methods: Neutron and Photon
Calculations, Chemical Rubber Company, Boca Raton, FL, 1991.

[75] R. Lewitt and S. Matej, Overview of Methods for Image Reconstruction from
Projections in Emission Computed Tomography, Proceedings of the IEEE, 91(2003),
1588-1611.

[76] S. Vandenberghe, et. al., Reconstruction of 2D PET data with Monte Carlo generated
system matrix for generalized natural pixels, Phys. Med. Biol, 51(2006), 3105-3125.

[77] M. Defiise, et al., Image Reconstruction Algorithms in PET, Positron Emission
Tomography. Basic Science and Clinical Practice, P. E. Valk, D. L. Bailey, et al. (eds.),
Springer, 4(2003), 91-114.

217

[78] G Marsaglia and A. Zarnan, Some portable very-long-period random number
generators, Comput. Phys. 8(1994), 117—121.

[79] httpz/lopengatecollaboration.hea1thgrid.org/

[80] http://www.gehealthcare.com/

[81] S. Deans, The Radon Transforms and some of it applications, Wiley, New York,
1983.

[82] W. Lu et al., Tomographic motion detection and correction directly in sinogram
space, Phys. Med. Biol, 47(2002), 1267-1284.

[83] P. Bruyant, et al., Correction of the respiratory motion of the heart by tracking of the

center of mass of thresholded projections: a simulation study using the dynamic MCAT
phantom, IEEE Trans. Nucl. Sci, 49(2003), 2159-2166.

[84] M. Jacobson, et al., Joint estimation of image and deformation parameters in
motion-corrected PET, IEEE Nucl. Sci. Symp. Conf Rec. 5(2003), 3290-3294.

[85] A. Rahrnim, et al., Motion compensation in histogram-mode and list-mode
reconstructions: beyond the event-driven approach, IEEE Trans. Nucl. Sci. 51(2004),
2588—2596.

[86] F. Lamare, et al., Incorporation of elastic transformations in list-mode based

reconstruction for respiratory motion correction in PET, IEEE Nucl. Sci. Symp. Conf Rec.
3(2005), 1740-1744.

[87] L. Liveratos, et al., Rigid-body transformation of list-mode projection data for
respiratory motion correction in cardiac PET, Phys. Med. Biol. 50(2005), 3313-3322.

[88] B. Thomdyke, et al., Reducing respiratory motion artifacts in positron emission
tomography through retrospective stacking, Med. Phys. 33(2006), 2632-2641.

[89] M. Dawood, et al., Lung motion correction on respiratory gated 3-D PET/CT images,
IEEE Trans. Med. Imaging 25(2006), 476-485.

[90] T. Li, et al., Model-based image reconstruction for four-dimensional PET, Med. Phys.
33 (2006), 1288-1298.

[91] F. Lamare, et al., Respiratory motion correction for PET oncology applications using
affine transformation of list mode data, Phys. Med. Biol. 52(2007), 121-140.

[92] WP. Segars. Development of a new dynamic NURBS-based cardiac-torso (NCAT)
phantom, PhD dissertation, The University of North Carolina, May 2001.

[93] C. J. Groiselle, et al., Monte-Carlo Simulation of the Photodetection Systems
Prototype PET Scanner Using GATE: A Validation Study, 2004 IEEE Nuclear Science
Symposium Conference, 5(2004), 3130-3132.

218

[94] J. C. Yanch, et al., Physically realistic Monte Carlo simulation of source, collimator

and tomographic data acquisition for emisssion computed tomography, Phys. Med. Biol.
37(1992), 853-870.

[95] Y. K. Dewaraja, et al., 3-D Monte Carlo-Based Scatter Compensation in Quantitative
I-131 SPECT Reconstruction, IEEE Hans. Nucl. Sci, 53(2006), 181-188.

[96] G W. Goerres, et al., PET-CT image co-registration in the thorax: influence of
respiration, Eur: J. Nucl. Med. Mol. Imaging 29(2002), 351-360.

[97] J. Kybic, et al., Fast parametric elastic image registration, IEEE Trans. Image
Processing, 12(2003), 1427-1442.

[98] F. Lamare, et al., Validation of a Monte Carlo simulation of the Philips
Allegro/GEMINI PET systems using GATE, Phys. Med. Biol. 51(2006), 943-962.

[99] J. R. McClelland, et al., A continuous 4D motion model from multiple respiratory
cycles for use in lung radiotherapy, Med. Phys. 33(2006), 3348—3358.

[100] S. A. Nehmeh, et al., Four-dirnensional (4D) PET/CT imaging of the thorax, Med.
Phys. 31(2004), 3179—3186.

[101] S. A. Nehmeh et al, Effect of respiratory gating on quantifying PET images of lung
cancer, J. Nucl. Med, 43(2002), 876-881.

[102] S. A. Nehmeh et a1, Reduction of respiratory motion artifacts in PET imaging of
lung cancer by respiratory correlation dynamic PET: Methohology and comparison with
respiratory gated PET, J. Nucl. Med, 44(2003), 1644-1648.

[103] M. Osman et al, Clinically significant inaccurate localization of lesions with
PET/CT: frequency in 300 patients, J. Nucl. Med, 30(2003), 240-243.

[104] K. Berger, K. Khurshid, R. McGough, J. L. Seamans. CTAC Misregistration in
PET/CT Cardiac Imaging: Correction with a Dedicated CTAC/PET Alignment
Visualization Application, and Automated Registration Algorithm, and ECG- Gated
CTAC. 12th Annual American Society of Nuclear Cardiology Scientific Session, San
Diego, CA, 2007.

[105] T. Pan, et al., Attenuation correction of PET images with respiration-averaged CT
images in PET/CT, J. Nucl. Med. 46(2005), 1481-1487.

[106] J. Qi, et al., List Mode reconstruction for PET with motion compensation: a
simulation study, Proc. IEEE Int. Symp. Biomed. Imaging, 2002, 413—416.

[107] M. Reyes et al, Model-based respiratory motion compensation for emission
tomography image, reconstruction, Phys. Med. Biol, 52(2007), 3579-3600.

[108] M. Reyes et al, Respiratory motion correction in emission tomography image

219

reconstruction, MICCAI 2005, J. Duncan and G Gerig (Eds), 2005, 369-376.

[109] W. Lu et al, Motion-encoded dose calculation through fluence/sinogram
modification, Med. Phys. 31(1), 2005, 118-127.

[110] W. Lu et al, Reduction of motion blurring artifacts using respiratory gated CT in
sinogram space: A quantitative evaluation, Med. Phys., 32(11), 2005, 3295-3304.

[111] F. Qiao et al, Region of interest motion compensation for PET image reconstruction,
Phys. Med. Biol, 52(2007), 2675-2689.

[112] F. Qiao, et al., A motion-incorporated reconstruction method for gated PET studies,
Phys. Med. Biol. 51(2006), 3769—3783.

[113] A. Reader, et al., One-pass list-mode EM algorithm for high-resolution 3D PET
image reconstruction into large array, IEEE Trans. Nucl. Sci. 49(2002), 693—699.

[114] H. Schoder, et al., Clinical implications of different image reconstruction
parameters for interpretation of whole-body PET studies in cancer patients, J. Nucl. Med.
45(2004), 559—566.

[115] W. Segars, et al., Modelling respiratory mechanics in the MCAT and spline-based
MCAT phantoms, IEEE Trans. Nucl. Sci. 48(2001), 89—97.

[116]D. Visvikis et al, A posteriori respiratory motion gating of dynamic pet images,
IEEE Nucl. Sci. Symp. Conf. Rec., 5(2003), 3276-3280.

[117] D. Visvikis, et al., Evaluation of respiratory motion effects in comparison with
other parameters affecting PET image quality, IEEE Nucl. Sci. Symp. Conf Rec. 6(2004),
3668—3672.

[118] L. A. Shepp and Y. Vardi, Maximum likelihood reconstruction for emission
tomography, IEEE Trans. Med. Img., 2(1982), 113-122.

[119] P. J. Green, Bayesian reconstructions from emission tomography data using a
modified EM algorithm, IEEE Trans. Med. 1mg. 9(1), 1990, 84-93.

[120] H. M. Hudson et al, Accelerated image reconstruction using ordered subsets of
projection data, IEEE Trans. Med. 1mg, 13(4), 1994, 601-609.

[121] http://www.owlnet.rice.edu/~elec539/Projects97/cult/report.html

[122] D. Visvikis, et al., Influence of OSEM and segmented attenuation correction in the
calculation of standardised uptake values for [18F ] FDG PET, Eur. J. Nucl. Med.
28(2001), 1326—1335.

[123] J. M. Ollinger, Model-based scatter correction for fully 3D PET, Phys. Med. Biol.
41(1996), 153-176.

220

[124] A. Gottschalk, P. V. Harper, F. F. Jimenez, and J. P. Petasnick, “Quantification of
the respiratory scanning with the rectilinear focused collimator scanner and the gamma
scintillation camera,” J. Nucl. Med, 7(1966), 243—251.

[125] M. V. Green, J. Seidel, S. D. Stein, T. E. Tedder, K. M. Kernpner, C. Kertzrnan, and
T. A. Zeffiro, “Head movement in normal subjects during simulated PET brain imaging
with and without head restraint,” J. Nucl. Med, 35(1994), 1538—1546.

[126] M. Menke, M. S. Atkins, and K. R. Buckley, “Compensation methods for head
motion detected during PET imaging,” IEEE Trans. Nucl. Sci, 43(1996), 310—317.

[127] Y. Picard, and C. J. Thompson, “Motion Correction of PET Images Using Multiple
Acquisition Frames,” IEEE Trans. Medical Imaging, 16(1997), 137-144.

[128] P. E. Kinahan, D. W. Townsend, T. Beyer, and D. Sashin, “Attenuation correction
for a combined 3D PET/CT scanner’, Med. Phys., 25(1998), 2046-2053.

[129] B.M.W. Tsui, W.P. Segars, and D.S. Lalush, “Effects of Upward Creep and
Respiratory Motion in Myocardial SPECT,” IEEE Trans. Nucl. Sci, 47(2000),
1192-1195.

[130] P. P. Bruyant et al., “Correction of the Respiratory Motion of the Heart by Tracking
of the Center of Mass of Thresholded Projections: A Simulation Study Using the
Dynamic MCAT Phantom,” IEEE Trans. Nucl. Sci, 49(2002), 2159-2166.

[131] DR Gilland, BA Mair, J E Bowsher, and RJ J aszczak. Simultaneous reconstruction
and motion estimation for gated cardiac ECT. IEEE Trans. Nucl. Sci, 49(2002),
2344—2349.

[132] R. R. Fulton et al., “Correction for Head Movements in Positron Emission
Tomography Using an Optical Motion-Tracking System,” IEEE Trans. Nucl. Sci,
49(2002), 116-123.

[133] W. Lu et al., "Tomographic motion detection and correction directly in sinogram
space,” Phys. Med. Biol, 47(2002), 1267-1284.

[134] S. S. Vedarn et al., “Quantifying the predictability of diaphragm motion during
respiration with a noninvasive external marker,” Med. Phys., 30(2003), 505-513.

[135] M. Osman et al, “Clinically significant inaccurate localization of lesions with
PET/CT: frequency in 300 patients,” J. Nucl. Med, 30(2003), 240-243.

[136] Cao et al., “Three-dimensional motion estimation with image reconstruction for
gated cardiac ECT,” IEEE Trans. Nucl. Sci, vol. 50, pp. 384-388, 2003.

[137] Paul Bi'rhler, Uwe Just, Edmund Will, Jorg Kotzerke, and Jorg van den Hoff, “An
Accurate Method for Correction of Head Movement in PET,” IEEE Trans. Medical
Imaging, 23(2004), 1 176-1 185.

221

[138] E. Gravier and Y. Yang, “Motion-compensated reconstruction of tomographic
image sequences,” IEEE Tran. Nucl. Sci, 52(2005), 51-56.

[139] W. Lu et al., “Reduction of motion blurring artifacts using respiratory gated CT in
sinogram space: A quantitative evaluation,” Med. Phys., 32(2005), 3295-3304.

[140] R. Zeng et al., “Respiratory motion estimation from slowly rotating x-ray
projections: Theory and simulation,” Med. Phys., 32(2005), 984-991 .

[141] F. Lamare et al. “List-mode-based reconstruction for respiratory motion correction
in PET using non-rigid body transformations,” Phys. Med. Biol, 52(2007), 5187-5204.

[142] R. Zeng et al., “Estimating 3-D Respiratory Motion From Orbiting Vrews by
Tomographic Image Registration,” IEEE Trans. Medical Imaging, 26(2007), 153-161.

[143] R. L. Wahl et al., Principles and Practice of Positron Emission Tomography,
Lippincott Williams & Wilkins, 2002.

[144] D. L. Bailey et al., Positron Emission Tomography: Basic Sciences,
Springer-Verlag London Limited, 2005.

[145] D. Mah et al., “Technical aspects of the deep inspiration breath-hold technique in
the treatment of thoracic cancer,” Int J Radiat Oncol Biol Phys, 48(2000), 1175-1185.

[146] K. E. Sixel et al., “Deep inspiration breath hold to reduce irradiated heart volume in
breast cancer patients,” Int J Radiat Oncol Biol Phys, 49(2001), 199-204.

[147] S. A. Nehmeh et al., “Deep-Inspiration Breath-Hold PET/CT of the Thorax,” J.
Nucl. Med, 48(2007), 22-26.

[148] O. Camara et al., “Evaluation of a thoracic elastic registration method using
anatomical constraints in oncology”, 2nd. Joint Conf. of the IEEE Engineering in
Medicine and Biology Society, 2002.

[149] Mattes D, Haynor D, Vesselle H, Lewellen T and Eubank W, PET-CT image
registration in the chest using free-form deformations, IEEE Trans. Med. Imaging,
22(2003), 120-128.

[150] K. Khurshid et al., “Automated cardiac motion compensation in PET/CT for
accurate reconstruction of PET myocardial perfusion images,” Phys. Med. Biol, 53(2008),
5705-5718.

[151] J. B. A. Maintz et al., “A survey of medical imaging registration,” Medical Image
Analysis, 2(1998), 1-37.

[152] D. L. G Hill et al., “Medical image registration,” Phys. Med. Biol, 46(2001), 1-45.

[153] T. Li, et al., “Model-based image reconstruction for four-dimensional PET,” Med.

222

Phys., 33(2006), 1288-1298.

[154] D. A. Low et al., “Novel breathing motion model for radiotherapy,” Int. J. of Rad
Onc. Bio. Phys. , 63(2005), 921-929.

[155] Y. Deng et al., “A Monte Carlo model for positron emission tomography imaging,”
Int. J. of App]. Electrom, 28(2008), 79-85.

[156] R. L. Harrison et al., “A public-domain simulation system for emission tomography:
photon tracking through heterogeneous attenuation using importance sampling,” J. Nucl
Med, 34(1993), 6OP.

[157] H. Anger, Scintillation camera and positron camera in: Medical Radioisotope
Scanning, Vrenna: IAEA, (1959), 59-73.

[158] CW. Williams, et al., Design and performance characteristics of a positron emission
computed axial tomography-ECAT-II, IEEE Trans Nucl Sci, 1979, 619—627. °

[159] L. Lupton and N. Keller, Performance Study of Single —Slice Positron Emission
Tomography Scanners by Monte Carlo Techniques, IEEE Trans on Med. Imag. 2(1983),
154-167.

[160] R. Lewitt and S. Matej, Overview of Methods for Image Reconstruction from
Projections in Emission Computed Tomography, Proceedings of the IEEE, 91(2003),
1588-1611.

[161] S. Vandenberghe, et. al., Reconstruction of 2D PET data with Monte Carlo
generated system matrix for generalized natural pixels, Phys. Med. Biol, 51(2006),
3105-3125.

[162] M. Defrise, et al., 4-Irnage Reconstruction Algorithms in PET, Positron Emission
Tomography: Basic Science and Clinical Practice, P. E. Valk, D. L. Bailey, et al. (eds.),
Springer, 4(2003), 91-114.

[163] G Marsaglia and A. Zarnan, Some portable very-long-period random number
generators, Comput. Phys. 8(1994), 117—121.

223