Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Real-time imaging of elastic properties of soft tissue with ultrasound Zahiri-Azar, Reza 2005

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
831-ubc_2005-0720.pdf [ 17.97MB ]
Metadata
JSON: 831-1.0065559.json
JSON-LD: 831-1.0065559-ld.json
RDF/XML (Pretty): 831-1.0065559-rdf.xml
RDF/JSON: 831-1.0065559-rdf.json
Turtle: 831-1.0065559-turtle.txt
N-Triples: 831-1.0065559-rdf-ntriples.txt
Original Record: 831-1.0065559-source.json
Full Text
831-1.0065559-fulltext.txt
Citation
831-1.0065559.ris

Full Text

Real-Time Imaging of Elastic Properties Soft Tissue with Ultrasound by R e z a Z a h i r i - A z a r B . S c , Shar i f Un ive r s i ty of Technology, 2003 A T H E S I S S U B M I T T E D I N P A R T I A L F U L F I L L M E N T O F T H E R E Q U I R E M E N T S F O R T H E D E G R E E O F M A S T E R O F A P P L I E D S C I E N C E i n T h e Facu l ty of Gradua te Studies (Elec t r ica l and C o m p u t e r Engineering) T H E U N I V E R S I T Y O F B R I T I S H C O L U M B I A September 2005 © R e z a Z a h i r i - A z a r , 2005 Abstract Curren t imaging devices such as computed tomography ( C T ) , u l t rasound (US) and magnetic resonance imaging ( M R I ) are not d i rect ly capable of measuring the mechanical properties of soft tissue even though such measurement wou ld have a h igh c l in ica l demand. Elas tography w i t h the a id of ul t rasound has been wel l established i n the l i terature as a s t ra in imag ing technique. Unde r cer ta in condit ions, these s t ra in images can give a clear i l lus t ra t ion of the under ly ing tissue stiffness dis t r ibut ions which has been shown to provide useful c l in ica l informat ion. V i b r o Elas tography is another new imaging system that performs a transfer function analysis of the tissue mot ion . T h e shape of the transfer function can be analyzed further and the stiffness of tissue can be est imated from the magnitude of the transfer functions at low-frequencies. T h i s thesis introduces a fast and accurate mot ion t rack ing a lgor i thm wh ich is at the heart of b o t h s t rain imaging and stiffness imaging. T h e a lgor i thm achieves real-t ime performance (> 20 fps) wi thout any need for add i t iona l hardware and its overhead. T h e performance of the proposed method is evaluated quant i ta t ive ly according to its signal-to-noise ra t io , contrast-to-noise ra t io , dynamic range, resolut ion and sensi t ivi ty w i t h bo th s imula t ion da ta and phan tom data . A l s o , the computa t iona l efficiency of the a lgor i thm is compared w i t h current real- t ime mot ion t racking algori thms. T h e results show that i t is the most t ime efficient a lgor i thm to date. Fur thermore the performance of the proposed me thod is evaluated qual i ta t ively from the real-t ime images that are generated i n bo th tissue m i m i c k i n g phantoms and real tissues i n v ivo . B y using this method two real-t ime elastography packages have been implemented which can easily be c l in ica l ly applied. These implementat ions run at 35fps for s t ra in images and 2fps for transfer function images of 16,000 pixels on an Ul t r a son ix R P 5 0 0 ul t rasound machine. ii Contents Abstract ii Table of Contents iii List of Tables vii List of Figures viii Glossary xiv Acknowledgments xvi 1 Introduction 1 1.1 Background 1 1.2 Ul t r a sound Imaging 2 1.3 Elas tography M e t h o d s 4 1.4 A x i a l S t ra in 4 1.5 Previous Research on S t ra in Imaging 6 1.5.1 Gradient-based S t ra in E s t i m a t i o n 7 1.5.1.1 Displacement E s t i m a t i o n Techniques 7 1.5.1.2 S t ra in E s t i m a t i o n Techniques 10 1.5.2 Direc t S t ra in E s t i m a t i o n 10 1.5.3 M u l t i D imens iona l S t ra in Imaging 11 1.5.3.1 2D S t ra in E s t i m a t i o n 12 1.5.3.2 3D S t ra in E s t i m a t i o n 13 i i i CONTENTS 1.6 Compar i son M e t r i c s i n Elas tography 14 1.7 Performance O p t i m i z a t i o n i n Elas tography 17 1.8 Rea l - t ime S t ra in Imaging and Its C o m p u t a t i o n a l Issues 19 1.8.1 Current R e a l - T i m e M o t i o n T rack ing A l g o r i t h m s 20 1.9 V i b r o Elas tography 21 1.9.1 Transfer F u n c t i o n A n a l y s i s 22 1.10 Thesis Object ives 23 1.11 Thesis Ou t l ine 25 2 M o t i o n E s t i m a t i o n 27 2.1 S tandard T i m e D o m a i n Cross Cor re la t ion Technique 27 2.1.1 Pros and Cons of the S tandard M e t h o d 29 2.2 T i m e D o m a i n Cross Corre la t ion w i t h P r i o r Es t imates 32 2.2.1 T h e M a i n Structure of the A l g o r i t h m 34 2.2.2 P red i c t i on A l g o r i t h m 35 2.2.3 P ros and Cons of the A l g o r i t h m 37 2.3 Conc lus ion 37 3 C o m p u t a t i o n a l Cost 40 3.1 P r e l i m i n a r y Considerat ions and Compar i son Measures 40 3.2 C o m p u t a t i o n a l Cos t of the T D E a lgor i thm 41 3.3 C o m p u t a t i o n a l Cos t of the T D P E a lgor i thm 42 3.4 C o m p u t a t i o n a l Cos t of the P R S a lgor i thm 43 3.5 C o m p u t a t i o n a l Cos t of the C A M a lgor i thm 43 3.6 Compar i son of the C o m p u t a t i o n a l Costs 44 3.7 Conc lus ion 48 4 Stra in E s t i m a t i o n 49 4.1 Gradien t Opera tor 49 4.2 Least Squares S t r a in E s t i m a t i o n ( L S E ) 51 4.3 Higher Order N u m e r i c a l Differentiators 52 4.4 L S E w i t h Higher Order Po lynomia l s 53 CONTENTS v 4.5 Conc lus ion 54 5 S imulat ion Results 56 5.1 V a l i d a t i o n 56 5.2 I D S imula t ion of a 2D Image 57 5.2.1 Eva lua t i on of S t ra in Es t imators 58 5.2.2 S t ra in F i l t e r S tudy (SF) 62 5.2.2.1 Effect of L S E and M e d i a n F i l t e r on the S t ra in F i l t e r 64 5.2.3 Contrast- to-Noise R a t i o S tudy 65 5.2.3.1 Effect of L S E on Contras t - to-Noise R a t i o 68 5.2.4 A x i a l Reso lu t ion S tudy (Ra) 69 5.2.4.1 Effect of L S E and M e d i a n F i l t e r i n g on A x i a l Reso lu t ion 70 5.2.5 Qua l i t a t ive Eva lua t ion 72 5.3 Conc lus ion 79 6 Exper iments 80 6.1 E x p e r i m e n t a l Setup 80 6.2 P h a n t o m Mate r ia l s 82 6.3 S t ra in F i l t e r S tudy 83 6.4 Qua l i t a t ive Resul ts 84 6.5 Conc lus ion 88 7 Implementat ion 89 7.1 R e a l T i m e S t ra in Imaging 89 7.1.1 B-mode T h r e a d 90 7.1.2 S t ra in Imaging T h r e a d 90 7.1.2.1 S t ra in E s t i m a t i o n • • 92 7.1.2.2 Pos t Processing 94 7.1.2.3 D i sp l ay ing 97 7.1.3 Software Features 98 7.2 R e a l T i m e V i b r o Elas tography 100 7.2.1 Stiffness Imaging Th read 100 CONTENTS v i 7.2.1.1 Transfer Func t i on E s t i m a t i o n 102 7.2.1.2 Stiffness E s t i m a t i o n 103 7.2.1.3 Pos t Process ing and Di sp l ay ing 103 7.2.2 Software Features 105 7.3 Ver i f ica t ion 106 7.4 Conc lus ion I l l 8 S u m m a r y and Future W o r k 112 8.1 Summary 112 8.2 Future W o r k 114 Bibl iography 116 A R e a l T i m e M o t i o n T r a c k i n g A lgor i thms 126 A . 1 Phase R o o t Seeking ( P R S ) 126 A . 2 C o m b i n e d Autocor re l a t ion M e t h o d ( C A M ) 128 B Corre la t ion Based M o t i o n Track ing M e t h o d s 130 B . l Cor re la t ion E s t i m a t i o n M e t h o d s 130 B . 2 Subpixe l l ing 133 B . 3 Cor re la t ion Reca lcu la t ion 135 C C o m p u t a t i o n a l Issues 137 C l C o m p u t a t i o n a l Cos t of the L S E M e t h o d 137 List o f Tables 5.1 Compar i son of gradient operator w i t h numerica l differentiators. A l l values are av-eraged over 100 realizat ions. To compare the method the result from ideal s t ra in estimator is also inc luded. T h e table compares the a lgor i thms ' s tandard devia t ion and their resul t ing signal-to-noise ratios ( S N R ) 59 5.2 Compar i son of LSEi(N). A l l values are averaged over 100 realizations 60 5.3 Compar i son of LSE2(N). A l l values are averaged over 100 realizations 60 7.1 T h e current frame rate of the real t ime s train imaging system on U l t r a s o n i x R P 5 0 0 w i t h a single P e n t i u m 4, C P U 3 . 2 G H Z , 1 G B R A M for a window size of 1mm (52 samples) 100 v i i List of Figures 1.1 A n example of the ul t rasound B-mode image. T h e image is capture from the R P 500 ul t rasound machine (Ul t rasonix C o r p . , Burnaby , Canada) 3 1.2 I l lus t ra t ion of the static elastographic process. P r e - and post-compressed tissue states are also shown. T h e figure is reproduced, w i t h permission, from [38] 5 1.3 Af te r d i v i d i n g the R F A- l ines to a number of windows the mot ion at each window is est imated by using a mot ion t racking a lgor i thm. Taken from reference [38] 8 1.4 Simpl i f ied block d iagram of the elastographic process. T h e compression, together w i t h the in t r ins ic stiffness d i s t r ibu t ion (a) i n the tissue and the boundary condit ions, causes a s t ra in d i s t r ibu t ion to be set up i n the tissue (b). The s tat is t ical properties of the sonographic tissue mot ion es t imat ion process give rise to the S F . T h e s t ra in d i s t r ibu t ion filtered by the S F is the elastogram (c) 15 1.5 A sample of s t ra in filter is shown. For clar if icat ion dynamic range, sensi t ivi ty and sa turat ion have also been shown on the figure 16 1.6 Sketch showing the s implif icat ion of a l inear dynamic system into a static one under low frequency excitat ions [117] 22 1.7 Sketch showing the s implif icat ion of bo th s t ra in imaging and stiffness imaging process. 24 2.1 T h e T D E me thod splits the pre-compression R F A - l i n e into a number of over lapping windows and finds their corresponding windows i n the post-compression R F A - l i n e by cross correlat ion method. T h e impor tan t parameters are shown i n the figure. W is the length of each window and AW is the shift between windows. Some researchers use the overlap between windows instead of AW as a signal processing parameter. . 28 v i i i LIST OF FIGURES i x 2.2 P re - and post-compressed R F da ta (Left) , and the corresponding normal ized cross-correlations (Right ) . T h e peak i n the correla t ion function shows the displacement. . 29 2.3 Theore t ica l displacement image on a pair of s imulated R F frames on a homogeneous tissue (Top) and est imated displacement image on a pair of s imula ted R F frames when T D E w i t h cross correlat ion is used (Bo t tom) . T h e z-axis shows the ax ia l displacement, the y-axis shows the R F A- l ines and the x-axis shows the windows s tar t ing from the window prox imal to the probe to the d is ta l window deep i n the tissue. Large peaks are the result of selecting false peaks. T h e image shows how often does this happen when non-normal ized cross correlat ion is used 30 2.4 E s t i m a t e d displacement image on a pair of s imulated R F frames when T D E w i t h normal ized cross correlat ion is used (Top) fol lowing a 2D median f i l ter ing w i t h a kernel size of 3x3 (Bo t tom) . T h e remain ing false peaks are removed by using a median filter 31 2.5 In comput ing the displacement of the center window, T D P E exploi ts the fact tha t displacements of two of its neighbors (leading window to the left and upper above) have already been calculated 33 2.6 A sample of the correlat ion function used by T D P E is shown on the left. T h e correla t ion is not normal ized and i t has only been calculated for a smal l number of time-shifts. A sample of est imated displacement vectors generated by T D E and T D P E is shown on the right. E v e n though T D E uses normal ized correlat ion, i t may s t i l l f ind false peaks while T D P E can work even w i t h non-normalized correla t ion and estimates the mot ion correctly. 34 2.7 E s t i m a t e d displacement image on a pair of s imula ted R F frames when T D P E w i t h non-normal ized cross correlat ion is used (median filter is not applied) . T h e image shows that even though T D P E uses the non-normal ized correlat ion i t s t i l l finds the correct displacements 35 2.8 T h e search region of the T D E a lgor i thm needs to be greater than the largest possible mo t ion (a), whi le the search region of the T D P E a lgor i thm is fixed (b) 38 4.1 Sketch showing the process of es t imat ing the s t ra in from the displacement estimates. 50 LIST OF FIGURES x 4.2 L S E w i t h a kernel size of 5 is shown. A line is fit ted to the da ta points w i t h a least squares error and the value of the slope is then reported as the estimate of the gradient at the center of the kernel 51 5.1 E s t i m a t e d displacement at 1% st ra in on bo th s imulated homogeneous phan tom and phantom w i t h inclus ion when T D P E is used as a mot ion t rack ing a lgor i thm ( D = 40mm, W = l m m , overlap=75%, fa=5 M H z and fa=40 M H z ) 58 5.2 A sample of displacement estimates of a three layered phan tom at 1% compression w i t h its corresponding elastogram when different s t ra in estimators are used 61 5.3 T h e s t ra in filter of the T D P E a lgor i thm. Compress ion is done once from 0.001% to 100% i n logar i thmic steps (Top) and once from 0% to 15% i n 0.2% steps to increase the accuracy (Bo t tom) . 100 realizations are used to show the mean values (Top, B o t t o m ) and the devia t ion from the mean values (Bo t tom) . T h e gradient operator is used as a s t ra in est imator wi thout any fil tering ( D = 4 0 m m , W = l m m , overlap=75%, f0=5 M H z and / a = 4 0 M H z ) 63 5.4 Improvement i n s t ra in filters for T D P E when (a,b) L S E , (c,d) gradient operator fol lowing a 2D median filter and (e,f) 2 D median filter fol lowing a gradient operator is used as a s t ra in est imator. Compress ion is done from 0% to 15% i n 0.2% steps. 100 realizations are used to show the mean values. T h e corresponding improvements i n S N R e at 1% compression are also shown i n (b, d and f) (D=40 m m , W = l m m , overlap=75%, fQ=5 M H z and / s = 4 0 M H z ) 66 5.5 C N R e of T D P E plot ted as a function of contrast ra t io . C N R e is calcula ted from the S F and the strain and devia t ion values up to 6% stra in are used. U s i n g 0.2% s t ra in i n the inclus ion and 0.2% to 6% s t ra in i n the background ( D = 4 0 m m , W = l m m , overlap=75%, fQ=5 M H z and / 5 = 4 0 M H z ) 67 5.6 C N R e of T D P E using I D simulated da ta and LSE(2-12) as an s t ra in estimator. C N R e is calculated from S F and the s t ra in and devia t ion values up to 6% stra in are used, by using 0.2% s t ra in i n the inclus ion and 0.2% to 6% stra in i n the background (a). 100 realizations are used to show the mean values. T h e improvement i n C N R e at contrast rat io of 5 is also shown i n (b) (D=40 m m , W = l m m , overlap=75%, f0=5 M H z and fs =40 M H z ) 68 LIST OF FIGURES x i 5.7 To simulate a 3 layered phan tom, the stiffness of the inside layer was set to be 4 t imes stiffer than its background. T h e stiffness map is shown i n (a) and its corresponding estimated elastogram image at 1% compression is shown i n (b). L S E w i t h a kernel size of 6 is used as an s t ra in est imator and a 3 x 3 median fi l tering was appl ied after the s t ra in es t imat ion 70 5.8 Effect of fi l tering on the ax i a l resolution. For L S E the kernel size was varied from 2 to 30 (a) and for 2D M e d i a n fi l tering the kernel size was changed from 3 x 3 to 11 x 11 (b). 100 real izat ions are used to show the mean values. (D=40 m m , W = l m m , overlap=75%, / 0 = 5 M H z and / s = 4 0 M H z ) 71 5.9 T o simulate a phan tom w i t h a circle inclusion the stiffness of a circle inside a phan tom was set to be 4 t imes stiffer than its background. T h e stiffness map is shown i n (a) and its corresponding idea l s t ra in image at 1% compression is shown i n (b,c) 72 5.10 To demonstrate the dynamic range of the T D P E , the elastograms at different com-pression ratios are shown i n this figure. T h e inclus ion is s t i l l detectable at 10% compression (D=40 m m , W = l m m , overlap=75%, / Q = 5 M H z and / s = 4 0 M H z ) . . . 74 5.11 T h e improvement i n elastogram when the Least Squares S t r a in Es t ima to r ( L S E ) w i t h different kernel sizes are used instead of gradient estimator. A s the kernel size increases the elastograms improve rapidly. T h e smoother elastograms also show the loss of resolution. T h e boundary of the inclus ion fades as the kernel size becomes bigger (D=40 m m , W = l m m , overlap=75%, / D = 5 M H z and / s = 4 0 M H z ) 75 5.12 To demonstrate the effect of median fi l tering i n the elastograms' (SNRe) the or ig inal elastogram is compared w i t h median filtered elastograms w i t h different kernel sizes. A s the kernel size increases the elastograms improve rapidly. T h e filtered elastograms also show the loss of resolut ion since the boundary of the inc lus ion fades as the kernel size becomes bigger (D=40 m m , W = l m m , overlap=75%, / 0 = 5 M H z and / s = 4 0 M H z ) . 76 5.13 T h e improvement i n elastogram when the Least Squares S t ra in E s t i m a t o r ( L S E ) is used instead of gradient estimator. A s the kernel size increases the elastograms' (SNRe) improves r ap id ly whi le the resolution degrades (the boundary of the layeres fade) (D=40 m m , W = l m m , overlap=75%, fa=5 M H z and / s = 4 0 M H z ) 77 L I S T OF FIGURES x i i 5.14 T h e or ig ina l elastogram is compared w i t h M e d i a n F i l t e red elastograms w i t h different kernel sizes. A s the kernel size increases the elastograms' (SNRe) improve r ap id ly whi le the resolut ion degrades (the boundary of the layers fade) (D=40 m m , W = l m m , over lap=75%, / Q = 5 M H z and / s = 4 0 M H z ) 78 6.1 Overv iew of the exper imental set-up 81 6.2 T h e s t ra in filter of the T D P E a lgor i thm w i t h real da ta from experiments. Compres-sion is done from 0.2% to 5% i n steps of 0.2%. T h e result from 100 R F lines are used to show the mean values and the devia t ion from the mean values. T h e gradient operator is used as a s t rain estimator wi thout any fil tering 83 6.3 Image on right shows the phantom w i t h cy l ind r i ca l inclus ion and the image on left shows its corresponding B-mode image 84 6.4 Elas tograms of a phan tom w i t h 14 m m inc lus ion at 0.1%, 1% and 2%. L S E w i t h kernel size of 6 is used to estimate the s t rain. Elas tograms wi thout app ly ing a med ian filter (left) and median filtered w i t h kernel size of 3 x 3 (right) are shown ( D = 5 0 m m , W = l m m , overlap=75%, / o = 1 0 M H z and / s = 4 0 M H z ) 85 6.5 T h e image on the right shows the phan tom w i t h the cy l indr ica l inclus ion and the image on the left shows its corresponding B-mode image 86 6.6 Elas tograms of a phan tom w i t h 14mm inc lus ion at 0.5%, 1% and 2%. L S E w i t h kernel size of 6 is used to estimate the s t ra in . Elas tograms wi thout app ly ing a median filter (Left) and median filtered w i t h kernel size of 3 x 3 (right) are shown ( D = 5 0 m m , W = l m m , overlap=75%, / o = 1 0 M H z and / s = 4 0 M H z ) 87 7.1 F l o w C h a r t of the real t ime s t ra in imaging system 91 7.2 In the L S E method, neighboring kernels have large overlap w i t h eachother. Three overlapping kernels of size 5 is shown 93 7.3 Curren t interface of the real t ime s t ra in imaging system. T h e control box is on the left side of the screen and B-mode and s t ra in images are d rawn i n the right side of the screen 99 7.4 F l o w C h a r t of the real t ime v ibro elastography system 101 LIST OF FIGURES x i i i 7.5 Cur ren t interface of the real t ime s train imaging system. T h e control box is on the left side of the screen and B-mode and s t ra in images are drawn i n the right side of the screen 105 7.6 Screen shots of the interface and acquired images olbained w i t h two developed soft-ware. T w o types of tissue m i m i c k i n g phantoms are tested. T h e s t ra in image and stiffness image of a phantom w i t h inclus ion (top), and a three layered phan tom (bot tom) are shown 107 7.7 Pho tog raph of actuat ion system, w i t h the t ransrectal ul t rasound ( T R U S ) probe mounted i n i t (left). A sketch of the procedure is also shown (right) 108 7.8 Prosta te transverse B-mode images (left), transfer function images (middle) and s t ra in images (right) . T h e mot ion appl ied was 2 m m peak to peak, band-pass filtered (0.5-4.5Hz) whi te noise. T h e images show the prostate i n apex (top), mid-g land (center) and base (bottom) 109 7.9 Prosta te sagi t ta l B-mode images (left), transfer function images (middle) and s t ra in images (right) . T h e mot ion appl ied was 2 m m peak to peak, band-pass filtered (0.5-4.5Hz) whi te noise 110 A . 1 A M o d i f i e d N e w t o n search far the es t imat ion of the t ime delay. T h e root of the phase is i terat ively searched. In each step of the i terat ion, the phase is calcula ted and approximated by a linear function w i t h a slope equal to the transducers centroid frequency. T h e intercept of this l inear funct ion w i t h the abscissa is the new estimate for the t ime delay [90] 127 B . l E x a c t pos i t ion of the peak of cross correlat ion function is typ ica l ly es t imated by using in te rpola t ion techniques. T h e image is provided by Turgay E 134 Glossary I D One dimensional . 2 D T w o dimensional . 3 D Three dimensional . A C C Cost of conversion from real s ignal to analyt ic signal for a R F A - l i n e . A C C W Cost of conversion from real s ignal to analyt ic signal for a single R F window. B C C Cost of conversion from real s ignal to base-band signal for a R F A - l i n e . B C C W Cost of conversion from real s ignal to base-band signal for a single R F window. B-mode Brightness mode ul t rasound image. C C Cor re la t ion Cost . C A M C o m b i n e d Autocor re l a t ion M e t h o d . C N R e Contras t to Noise R a t i o i n Elastography. C S D Cross Spectra l Densi ty. C T C o m p u t e d tomography. D D e p t h of U l t r a sound Image. D R e D y n a m i c Range i n Elastography. E las tography Techniques for imaging mechanical properties of tissue. E l a s tog ram A n image obtained by elastography. L Leng th of an R F A - l i n e . L S E Least Square Es t ima to r . M R I Magne t i c resonance imaging. N C C Norma l i zed Cor re la t ion Cos t . x i v GLOSSARY x v P R S Phase R o o t Seeking. P S A Peak Search A l g o r i t h m . P S D Power Spectra l Densi ty. P S F Poin t Spread Func t ion . R F R a d i o Frequency D a t a . S N R e Signal to Noise R a t i o i n Elastography. S N R s Ul t r a sound Signal to Noise R a t i o . T D E T i m e D o m a i n M o t i o n E s t i m a t i o n . T D P E T i m e D o m a i n M o t i o n E s t i m a t i o n w i t h P r i o r Est imates . U S Ul t rasound . Z C T Zero Cross Track ing . Ra A x i a l Reso lu t ion i n Elastography. Ri La te r a l Reso lu t ion i n Elastography. c Speed of U l t r a sound i n tissue. AW W i n d o w Shift. W W i n d o w Leng th . n W N u m b e r of W i n d o w s i n a single R F line. fs Sampl ing Frequency of U l t r a sound Machine . fo Cent ro id Frequency of U l t r a sound Probe . X Ul t ra sound Pulse wavelength. E Young ' s modulus . di Displacement measurement of w indow i. Stiffness between windows i and i + H{{w) Transfer function from window i to window j. Acknowledgments • I would p r i m a r i l y l ike to thank my supervisor D r . T i m Salcudean for his guidance and support . W i t h o u t his guidance and inspi ra t ion throughout the process, none of this w o u l d have been possible. • I would also l ike to thank Trevor Hansen, K r i s D ick i e and Laurent Pelissier from U l t r a s o n i x M e d i c a l Corpora t ion , for their technical support and help. • M a n y thanks to a l l m y friends i n the Robo t i c s and C o n t r o l lab, J u l i a n for his help on the basics of U l t r a s o u n d imaging, E h s a n for l inear algebra, H a n i for signal processing, Hassan for continuous mechanics, Neerav for p rogramming and D a n n y for his cont inual ly helpful advices. • I owe m y loving thanks to m y wife Sara. W i t h o u t her support , encouragement and under-s tanding i t wou ld have been impossible for me to finish this work. • Las t ly , I a m sincerely grateful to m y parents who have always been the biggest mot iva t ion throughout m y studies and whose years of love and hard work created the oppor tun i ty for me to chase m y dreams. x v i Chapter 1 Introduction 1.1 B ackground T h e basis of medica l imaging is the measurement of a proper ty of tissue that varies w i t h t is-sue composi t ion. M e d i c a l images are formed by d isp lay ing these properties measured at mul t ip le locations i n the body. F r o m such images, a depic t ion of anatomy or pathology is gained. E a c h imaging moda l i ty i n common use, such as X - r a y , computed tomography ( C T ) , u l t rasound (US) and magnetic resonance imaging ( M R I ) , measures a different property of tissue. B u t none of the properties measured by those modali t ies describe d i rec t ly the mechanical properties of tissue. Mechan ica l properties are properties that describe the deformation and mechanical behavior of tissue under a stat ic or va ry ing force. E la s t i c i t y (Young's modulus) , compressibi l i ty (Poisson's rat io) , stiffness, viscosi ty and density are examples of mechanical properties. Crea t ing images f rom one or more mechanical properties is useful i n medica l practice, i f the var ia t ion of the measured property creates sufficient contrast i n the image to depict anatomy or pathology. Elas tography was first in t roduced i n [81] and i t aims to produce a new type of image that depicts the mechanical properties of tissue. T h e idea is to apply an external deformation to the tissue (such as ax ia l compression) or to use an available in ternal deformation (like the heart beat or deformation of ar ter ia l walls from osci l la t ing b lood pressure) and investigate its corresponding interactions inside the tissue. A n imaging device is employed to record these deformations and 1 1.2 Ultrasound Imaging 2 interactions. A mot ion est imat ion a lgor i thm is then appl ied to find the displacement of each region i n sequences of images. F i n a l l y different mechanical properties of each region are est imated from the recorded displacements. E las tography has been wel l established i n the l i terature as a s t ra in imaging technique [18,78-82]. T h e s t ra in dis t r ibut ions i n tissues i n response to an external deformation are closely related to the d i s t r ibu t ion of tissue elastici ty [108]. It has been shown that the s t ra in images are useful i n different c l in ica l appl icat ions such as tumor detection i n the breast [30,39,40] , b ra in [8] and prostate [69], mon i to r ing h igh intensity focused ul t rasound lesions [95], imaging the m y o c a r d i u m [54], s tudy ing deep venous thrombosis [1,32,100], s tudy ing renal pathology [33], moni to r ing the rmal changes and abla t ion [114,128], intravascular applicat ions such as f inding vulnerable plaque [25,26], mechanical behavior of sk in [134] and obta ining the mechanical s t ructura l properties of tissue [50]. 1.2 Ultrasound Imaging U l t r a s o u n d has become a widely used medica l scanning method due to its cost, ab i l i ty for real t ime performance, quick setup procedures and easy access to its d ig i ta l data . Therefore i t has been wide ly used i n mo t ion est imation, flow es t imat ion and elastography. In this thesis, the tissue's deformation is imaged by a conventional u l t rasound system. A n u l t rasound imaging system acquires da t a through the generation of an u l t rasound wave directed toward the area to be examined, followed by measurement of the echoes generated by the interact ion of the ul t rasound wave w i t h the tissue. T h e system consists of an u l t rasound transducer for generating the ul t rasound pulses and measurement of the echoes, and a computa t ion system to convert the echo ampli tudes into an image. T h e ul t rasound transducer sends out a short burst of u l t rasound, followed by a per iod of silence to l is ten for the re turning echoes. Fo r example, the u l t rasound wave t ransmit ted to the body strikes an interface and is pa r t i a l ly reflected back to the transducer. The ampli tude of the reflected echo depends on the ul t rasonic properties of the tissue at the interface. The t ime between the sent and the received echo is used to calculate the depth of the interface assuming the speed of sound is constant throughout the media . T h e form of the t ransmi t ted wave is an ampl i tude modula ted signal w i t h a fixed carrier frequency 1.2 Ultrasound Imaging 3 determined by the probe. T h e re turning echoes are rap id ly sampled dur ing the l is tening interval . T h e unprocessed samples are known as R F data . A set of R F da ta collected along one ax ia l l ine is cal led an R F A - l i n e (where A stands for ampl i tude) . The set of R F lines (collected on a co lumn- to-column basis) are referred to as RF- images throughout the thesis. E a c h R F da ta set is first amplif ied. T h e amplif icat ion increases as a function of t ime to compensate for the fact that echoes from deep locations are weaker and more at tenuated than echoes from shallow locations. T h i s operat ion is called t ime gain compensat ion ( T G C ) and is necessary to compensate for the a t tenuat ion of ul t rasound i n tissue. T h e next steps are envelope detection of the R F da ta and conversion to regular spat ia l coordinates (called scan conversion). A n image formed this way is cal led a B-mode (Brightness mode) image or B-scan (see F igure 1.1) because the image pixels have brightness propor t iona l to each ampl i tude [21,61,73]. F igu re 1.1: A n example of the ul t rasound B-mode image. T h e image is capture from the R P 500 ul t rasound machine (Ul t rasonix Corp . , Burnaby , Canada) . 1.3 Elastography Methods 4 1.3 Elastography Methods St ra in est imat ion methods based on ul t rasound fal l in to two m a i n groups depending on their external deformation method: 1) Me thods where a quasi-static compression is appl ied to the tissue and the result ing components of the s t ra in tensor are est imated; and 2) methods where a low-frequency v ib ra t ion is appl ied to the tissue and the result ing tissue behavior is inspected. T h e first method, s t ra in imaging w i t h a quasi-static compression [18,75,78-81], is performed by (i) obta in ing a set of ul t rasonic echo signals from a target, (ii) subject ing the target to a smal l ax ia l deformation and (iii) ob ta in ing a second set of echo signals. M o t i o n s along the direct ion of the appl ied load are est imated by performing piecewise m o t i o n es t imat ion on corresponding pairs of s ignal segments (Figure 1.2). S t ra in es t imat ion algori thms can then be appl ied to generate the s t ra in images named elastograms. It has been shown that the s t ra in dis t r ibut ions i n tissues are closely related to the d i s t r ibu t ion of tissue elast ici ty [108]. Fast acquis i t ion and /o r computa t ion are not an impor tant factor i n this k i n d of imaging, because the scatterer m o t i o n can be control led by the appl ied compression and a l l the computa t ion is done off l ine. S t ra in imaging can also be performed by app ly ing a dynamic compression (free hand deformation or low frequency v ib ra t ion to the tissue) [62,64,65,83,84,133,135] . U l t r a sound echo signals from the target are obtained simultaneously. T h e echo signals can be recorded for off line processing or they can be processed on the fly i n real t ime. Elas tograms are generated s imi la r ly to previous methods. Fast da ta acquis i t ion is an impor tan t factor i n th is k i n d of imaging and, i f real-t ime processing is preferred, fast computa t ion is also necessary. 1.4 Axial Strain A l t h o u g h images of ax ia l s t ra in , la teral s t ra in , Poisson's ra t io and Young ' s modulus can be generated [56,78, 82] i n this thesis we deal only w i t h the es t imat ion and imaging of ax ia l s t ra in wh ich is called an elastogram. 1.4 Axial Strain 5 Pre-compression Post-compression RFA-Iine RF A-iine Figure 1.2: I l lus t ra t ion of the static elastographic process. P re - and post-compressed tissue states are also shown. T h e figure is reproduced, w i t h permission, from [38]. W h e n an elastic med ium, such as tissue, is compressed by a constant uniaxia l stress, a l l points i n the med ium experience a result ing level of longi tud ina l s t ra in whose ma in components are along the axis of compression. If one or more of the tissue elements has a different stiffness parameter t han the others, the level of s t ra in i n that element w i l l generally be different. A stiffer tissue element w i l l generally experience less s t ra in that the softer one. T h e ax ia l s t ra in is est imated i n one dimension from the analysis of ultrasonic signals obta ined from standard medical u l t rasound. T h i s is accomplished by (i) acquir ing a set of digi t ized R F echo A- l ines from the tissue region of interest, (ii) compressing the tissue w i t h the ultrasonic transducer (or w i t h a combinat ion of transducer and compressor) along the ul trasonic radia t ion axis by a smal l amount (about 1%, or less, of the tissue depth) and (iii) acquir ing a second, post-compression set of echo lines from the 1.5 Previous Research on Strain Imaging 6 same region of interest. A s explained before this approach can also be used i n real t ime by acqui r ing continuous R F da ta from the tissue while an external deformation is appl ied to the tissue. Cor responding echo lines are then subd iv ided into smal l t empora l windows w h i c h are compared pair-wise by using mot ion est imat ion techniques, from wh ich the change i n a r r iva l t ime of the echoes before and after compression can be est imated (F ig . 1.3). D u e to the smal l magni tude of the appl ied compression i n static compression, or h igh da ta acquis i t ion rate i n dynamic compression, there are only smal l distort ions of the echo lines, and the changes i n arr ival t imes are smal l . T h e loca l long i tud ina l s t ra in is then estimated as E q u a t i o n (1.1): where t\a is the ar r iva l t ime of the pre-compression echo from the current window, tn, is the arr ival t ime of the pre-compression echo from the next window, t^a is the a r r iva l t ime of post-compression echo from the current window and t^h is the ar r iva l t ime of the post-compression echo from the next window [81]. Therefore E q u a t i o n (1.1) can be rewri t ten i n the form of time-shifts as E q u a t i o n (1.2) where t (i) is the time-shift between windows i n the indexed segment pair and A T is the spacing between windows. T h e windows are usual ly t ranslated i n overlapping steps a long the axis of the echo l ine, and the calcula t ion is repeated for a l l depths. It has been assumed that speckle m o t i o n adequately represents the under ly ing tissue m o t i o n for sma l l un iax ia l compressions [18]. 1.5 Previous Research on Strain Imaging en = (tlb — ha) — (t2b — tla) (1.1) At (i) - At (i - 1) A T (1.2) Fo l lowing the seminal paper by O p h i r et al [81], a number of methods have been developed and a number of methods have been borrowed from other research areas to evaluate the s t ra in d i s t r ibu t ion due to an applied external force. M o s t methods fall in to two ma in categories; gradient-1.5 Previous Research on Strain Imaging 7 based s t ra in es t imat ion and direct s t ra in es t imat ion [121]. E a c h of these categories is discussed i n detai l i n fol lowing sections. 1.5.1 Gradient-based Strain Estimation A x i a l s t ra in i n elastography is typ ica l ly est imated w i t h gradient-based methods from the gradi-ent of the displacement (time-shift or phase-shift) estimates. Therefore m o t i o n t r ack ing algori thms are the m a i n part of these methods. A l t h o u g h gradient operat ion can be swapped w i t h other s t ra in es t imat ion techniques, the te rm gradient-based methods are most ly used for this family of a lgori thms. C o m m o n displacement and s t ra in es t imat ion techniques are presented i n the following sections. 1.5.1.1 Displacement Estimation Techniques M o t i o n es t imat ion methods on ul t rasound images were first used for b lood flow es t imat ion [12, 115] where a window i n one R F A - l i n e is t racked i n i ts corresponding R F A - l i n e (Figure 1.3). T i m e domain cross correlat ion was among the first methods [18,37,43,47,78-81,110] i n wh ich the time-shift is est imated as the peak of the cross-correlation function between the pre- and post-compression R F A- l ines (see A p p e n d i x B . l ) . Al te rna t ive ly , cross correlat ion can be done on the envelope of the signals [59,119,120] instead of R F da ta to estimate the displacements. A s the evaluated cross correlat ion function is spa t ia l ly discrete, i t must be interpolated to detect the exact pos i t ion of the peak. In order to do subpixe l ing a number of in terpola t ion techniques have been in t roduced i n the l i terature (see A p p e n d i x B . 2 ) . T h e quadratic parabola fitting me thod was in t roduced i n [37] which is equivalent to finding the zero crossings of the derivat ive of the signal . T h e Sine in terpola t ion was in t roduced i n [16], and because of the sinusoidal nature of the u l t rasound signal , the cosine fitting method was used i n [16,24]. If needed, the correct value of the correlations can also be recalculated to measure the accuracy of the mot ion t rack ing process at each window (see A p p e n d i x B . 3 ) . In add i t ion to t ime domain methods, a number of phase domain methods have also been intro-duced. These involve the est imat ion of the phase-shift from the complex cross correlat ion of the 1.5 Previous Research on Strain Imaging 8 Pre-Compression Data W i n d o w "Post-Compression Data Window F igu re 1.3: After d iv id ing the R F A- l ines to a number of windows the m o t i o n at each window is es t imated by using a mot ion t rack ing a lgor i thm. Taken from reference [38]. pre- and post-compression R F A- l ines . Due to aliasing, phase d o m a i n methods fai l when a single-step displacement exceeds a quarter of the ultrasonic wavelength. T o avoid al iasing, phase unwrap methods have been in t roduced to extend the al iasing l imi t from A / 4 to A / 2 [75-77,131]. It has been suggested i n Phase R o o t Seeking ( P R S ) [90] to use Newton's iteration to unwrap the phase-shift (see A p p e n d i x A . 1 ) . T o avoid al iasing, the C o m b i n e d Au toco r r e l a t i on M e t h o d ( C A M ) [102] suggested the use of the envelope correlat ion coefficient (see A p p e n d i x A . 2 ) . In [113] al iasing is s imply ignored and the gradient operator is used di rect ly on the phase-shift wi thou t convert ing the phase shift to time-shift. Because phase informat ion is used directly, the author calls this method "direct s t rain es t imat ion" . A c c o r d i n g to the categorizat ion i n this thesis, the method falls into gradient-based methods since i t uses the gradient operator on the phase shift estimates to estimate the s t ra in . 1.5 Previous Research on Strain Imaging 9 A number of s t ra in imaging algori thms, instead of using cross correlat ion, use feature extrac-t i o n methods to f ind the displacement. The Zero Cross T rack ing a lgor i thm ( Z C T ) [109] identifies the zero-crossing instants of the pre- and post-compression R F A- l ines . T h e sign change at the occurrence of the zeros is used to ob ta in two groups of zero-crossing locations (corresponding to a posi t ive or negative sign change) for each R F A - l i n e . A tempora l t rack ing of such zero-crossings for successive steps or a cumulat ive accumulat ion of the displacements of the zero-crossings can be performed at each compression step. T h e distance between the zeros at each window is equal to the displacement estimate. T h e Peak Search A l g o r i t h m ( P S A ) [34] is another mo t ion t rack ing method that works s imi la r ly to Z C T but instead of f inding the zeros of the R F signals i t finds the peaks of R F signals. A peak assignment a lgor i thm is then appl ied to assign each peak i n the pre-compression R F da t a to i ts corresponding peak i n the post-compression R F data. T h e shift between the peak i n pre- and post-compression R F da ta is equal to the displacement estimates. A n i terative a lgor i thm has also been in t roduced [117] that uses bo th time-shifts and loca l s tretching factors to increase the accuracy of mot ion es t imat ion. It estimates the displacements in i t i a l l y by doing t ime domain cross correlat ion. T h e n , loca l s tretching factors are calcula ted according to the es t imat ion of the displacements. Fo l lowing that , s tretching factors are appl ied to the pre-compression R F da ta to increase the correlat ion between the pre- and post-compression signals. T h e new pre-compression R F da ta is used i n the next i te ra t ion to find the displacements again. Es t ima ted mot ions at each i tera t ion step are added to the previously est imated mot ion . T h e i tera t ion stops after a predefined number of i terations or when a certain s topping condi t ion is satisfied. T h e staggered s t ra in est imator as another i terative method was in t roduced i n [110]. It estimates the s t rain i n nonoverlapping windows at each step and shifts the windows by a smal l po r t ion for the next step and estimate the s t ra in again. T h e s t ra in estimates for a l l such shifts are staggered to produce the final elastogram. 1.5 Previous Research on Strain Imaging 10 1.5.1.2 S tra in E s t i m a t i o n Techniques After finding the time-shifts with any of the above methods, the strain is then estimated according to Equation (1.2). Since estimating the strain from the motion estimate is a slope estimation problem, all previously used slope estimators can be used as strain estimators. The gradient operator was used in the first strain estimation method [81]. It was used on the previously obtained time-shift estimates to estimate the strain. Since the gradient operator introduces additional amplification of the noise into the strain estimation process [60,121], filtering methods such as averaging or median filtering, can be used before and/or after applying the gradient operator [57,109,118]. Alternatively the strain can be estimated by a least squares strain estimator [49] or higher order numerical differentiator with smaller error instead of simple gradient operator. 1.5.2 Di rec t S t ra in Es t imat ion Direct methods such as the adaptive strain estimator [4] and spectral strain estimation methods [59,60,121] estimate the strain directly without involving the use of motion estimation algorithms and the gradient operator. The main idea behind the spectral strain estimators [59,60,121] is based on the Fourier scaling property, i.e., that a compression or expansion of the time-domain signal should lead to an expansion or compression of its power spectrum. Unlike gradient-based phase domain methods, spectral methods do not use the phase shift of the echo signals. They use the whole power spectrum. The adaptive strain estimator [4] is another direct method that uses the local stretching factors as an estimator of the local strains. It uses an iterative algorithm that adaptively maximizes the correlation between the pre- and post-compression echo signals by appropriately stretching the latter. The final stretching factors are then reported as local strains. This method is different from [107,117], where they use the stretching factor to improve the correlation between pre- and post-compression signals and then estimate the displacement. Here, the stretching factors themselves represent the local strain. 1.5 Previous Research on Strain Imaging 11 T h e correlat ion coefficient i tself can be used to estimate the s t ra in direct ly. T h e correlat ion between the pre- and post-compression echoes can decrease w i t h appl ied s t ra in . However, decor-re la t ion itself has been used to estimate delay and /o r time-shifts. Var ious researchers used the correla t ion coefficient to estimate tissue mot ion [28,116]. In [6], i t was proposed to use the decorre-l a t ion coefficient (1 minus the correlat ion coefficient) for the envelope signals for free hand elast ici ty imaging but i t was demonstrated i n [122] that the correlat ion coefficient has a poor precision as a s t ra in est imator. Because of its speed and s impl ic i ty , this estimator may be a valuable too l for free-hand real t ime elast ici ty imaging. However by using s imulated data, i t was also shown i n [2] that changes i n the center frequency and S N R s introduce unknown errors i n the corre la t ion coefficient. 1.5.3 M u l t i Dimensional St ra in Imaging T y p i c a l l y i n elastography, the ax ia l component of the s t rain is est imated by t a k i n g the gradient of the ax ia l (along the beam propagat ion axis) displacement occurr ing after a quasi-static tissue compression [18,78-81]. In general, however, the tissue mot ion that occurs du r ing compression is three d imensional (3D) . Because the la teral (perpendicular to the beam propagat ion axis and i n the scan plane) and elevational (perpendicular to the beam propagat ion axis and to the scan plane) mot ions are not measured, two major disadvantages may occur [56]. F i r s t , the ax ia l elas-togram takes into account only a smal l part of the mechanical tissue mot ion informat ion . Second, undesirable la teral and elevational motions are the p r imary causes of signal decorrelat ion [51,52]. T h e la tera l or elevational decorrelation can be prevented by appropriate confinement of the tissue under s tudy [51,52]. However, the confinement may not always be prac t ica l i n c l in ica l appl i -cations, especially when the tissues under s tudy are not easily accessible. Therefore, confinement of the tissue has not been considered i n this thesis. Ini t ia l ly , mul t id imens iona l mot ion es t imat ion started i n the field of b l ood veloci ty es t imat ion, to estimate veloci ty components other than the ax ia l one. Later , i n the field of elast ici ty imaging, some of the methods developed i n the flow measurement field were borrowed to f ind ax ia l , la teral and even elevational components of the s t rain. M a n y papers i n the li terature have dealt w i t h the problem of mot ion es t imat ion i n 2D [10,12, 1.5 Previous Research on Strain Imaging 12 13 ,19 ,29 ,56 ,71 ,93 ,115 ,135] , that t r y to find components of the m o t i o n i n the scan plane (axia l and lateral) , and i n 3D [13,42,55,58,63,68,74] , that t ry to find a l l the spat ia l components of the m o t i o n (axial , la teral and elevational). These two categories are discussed i n de ta i l i n the following sections. 1.5.3.1 2 D Stra in E s t i m a t i o n Simi la r to I D s t ra in imaging, t ime domain cross correlat ion can be used i n 2 D mot ion es t imat ion by t rack ing the mot ion once i n ax ia l d i rect ion and once i n lateral d i rec t ion independent ly [56] or i t can be extended to 2D block match ing instead of I D vector match ing [10,12,19,42,115]. 2 D block match ing methods detect the peak pos i t ion of the 2D cross-correlation funct ion evaluated w i t h respect to each loca l R F echo data . Because of the trade-off between the smoothness of mo t ion t r ack ing and spat ia l resolut ion depending on the size of the 2D block, i t was suggested i n [135] to use a hierarchical coarse-to-fine approach i n the selection of the sizes of ma tch ing blocks. A me thod s imilar to 2D block matching , called 2D companding [19], was developed for elas-tography to reduce decorrelation noise i n elastograms using appl ied ax ia l strains and correct ing for the la teral decorrelation (elevational mot ion was ignored). It uses speckle b lock-match ing tech-niques [12,115] to track the 2 D mot ion of the scatterers after static compression. T h i s me thod has recently been extended to 3D companding [20,46], i n order to compensate for bo th lateral and elevational displacement. In [70] the speckle-tracking technique was used i n lateral d i rect ion. T h e y showed that la teral displacement images obta ined by t r ad i t iona l 2D correlat ion t racking are noisy. T h e y also showed theoret ical ly that , when envelope da ta are used to track la teral mot ion , the variance of the la teral displacement es t imat ion is larger t han the variance of the ax ia l displacement estimate by orders of magni tude. Based on these results, they proposed the est imat ion of la tera l displacement by m a k i n g the assumption of isotropic incompress ibi l i ty i n the tissue and, thus, m a k i n g use of the precision of the ax i a l measurement. In [56,58] the lateral component of the s t rain was est imated by t rack ing the phase of single A- l ines i n the lateral direct ion th rough the use of a weighted in terpola t ion method . Since the resolut ion i n the lateral d i rect ion is poor, they used in terpola t ion techniques to reconstruct the 1.5 Previous Research on Strain Imaging 13 R F signal i n the la tera l d i rect ion [56,70]. Accura te es t imat ion of la teral s t ra in made it possible to estimate the Poisson's ra t io , wh ich is denned as the ra t io of the la teral s t ra in over ax ia l s t ra in , and is an impor tan t mechanical proper ty of tissue [56,97,98]. It has been shown that , i f there is a s t ra in contrast between an inc lus ion and the background, the Poisson elastogram is able to show whether this s t rain contrast is due to a Poisson's rat io contrast or an elastic modulus contrast. Fur thermore , by detecting the m o t i o n i n the la teral and /o r elevational di rect ion, correct ion algori thms were used to increase the accuracy of mot ion detection i n the ax ia l di rect ion [56,58]. In addi t ion to 2D t ime domain methods, a number of 2 D frequency domain methods have also been in t roduced that involve the es t imat ion of the phase-shift from the complex correlat ion of the pre- and post-compression R F images. A phase correlat ion method was used i n [131] to f ind the mot ion i n each block. T h e 2 D Fourier t ransform of a 2 D block i n pre-compression R F da t a was mul t ip l i ed by the conjugate of a 2D Fourier t ransform of a 2D block i n post-compression R F data , and the result was normal ized by the no rm of the mul t ip l i ca t ion . T h e inverse Fourier t ransform was then appl ied to the result to find the estimated mot ion i n a 2D window. A l t e rna t ive ly an i terative phase match ing method was in t roduced i n [112], wh ich uses the 2D phase characteristics of the pre-compression R F as the index to i terat ively search for the corresponding loca l da ta i n the post-compression R F data . 1.5.3.2 3 D Stra in E s t i m a t i o n 2D mot ion t rack ing methods are used for es t imat ing the mot ion only i n ax ia l and lateral d i rec t ion i n the plane of ul t rasound image. In order to find a l l three components of the mot ion , elevational mo t ion is also need to be est imated. T w o dist inct estimators for t rack ing elevational mo t ion can be developed using (i) a single plane or (ii) mul t ip le planes. A m o n g the methods that use a single plane to track the elevational mot ion , i n [7] the correlat ion coefficient of the signal envelope was used to generate lateral displacement images. A s a s imi lar approach when on ly one plane is available, mot ion i n the elevational d i rec t ion can be tracked by using the decorrelation informat ion [66]. B u t these approaches have several l imi ta t ions . To track the mot ion i n mul t ip le plane one needs to capture mul t ip le planes. There are three dist inct ways i n which mul t ip le planes can be generated i n the elevational d i rect ion. F i r s t ly , this 1.6 Comparison Metrics in Elastography 14 may be done by d isp lac ing a I D linear array i n the elevational direct ion and acqui r ing successive frames at different para l le l elevational planes. T h i s me thod w i l l most probably not be workable i n an exper imental set t ing and medica l applicat ions since i t is not always possible to keep the tissue steady for a long t ime. Secondly, a 1.5D array can be used to provide three to four elevat ional planes by fir ing different sets of elements i n the elevational di rect ion. T h i s me thod m a y y i e ld the elevational component . Las t ly , a 2D (or N x N) array may be used, p rovid ing as many planes i n the elevational d i rec t ion as there are beams i n the la teral di rect ion. In this case, the elevat ional mo t ion can be es t imated i n a plane perpendicular to a cer ta in scanning plane. A s s u m i n g that mul t ip le planes of ul t rasound images are available, previous m o t i o n t r ack ing methods can be extended to 3 D . A n interplanar t r ack ing method, s imi lar to the one used for i n -terbeam lateral t r ack ing i n [56] was applied i n [58] to interpolate the R F line between the or ig ina l beams to increase the resolut ion i n bo th lateral and elevational di rect ion. A x i a l , la tera l and eleva-t ional motions are found independently by using I D t ime domain cross correlat ion and correct ion a lgor i thm is used to recorrelate a l l components of the mot ion . It should be noted that , a l though i n -terpolat ing R F lines may increase accuracy, i t is very t ime consuming and i t is not always appl icable i f the computa t iona l cost is an impor tant factor. In [63,68], 2 D block match ing is used i n each plane and its neighboring planes to f ind and correct for bo th la tera l and elevational motions of the current plane. The P R S , wh ich is a phase domain a lgor i thm, is then used i n the ax ia l d i rect ion to find the ax ia l component of the mot ion . If the 2D phase array u l t rasound probe is used to acquire R F data , the C A M method can be used to find the mot ion i n 3 D [74]. 1.6 Comparison Metrics in Elastography T h e entire elastographic process can be described using three fundamental concepts [127]. T h e block d iagram is shown i n F igu re 1.4. T h e first concept is the contrast transfer efficiency ( C T E ) [48,91], that determines the d i s t r ibu t ion of the ideal loca l ax ia l strains i n the compressed tissue according to its in t r ins ic tissue parameters. C T E is independent of the imaging system and s t ra in est imat ion a lgor i thm and i t describes the fundamental relat ionship between modulus and s t ra in 1.6 Comparison Metrics in Elastography 15 contrasts. Input Contrast Strain Filter Contrast to Tissue Transfer • Noise Modulus Efficiency | i Ratio (a) (b) (c) Figure 1.4: Simpl i f ied block d iagram of the elastographic process. The compression, together w i t h the in t r ins ic stiffness d i s t r ibu t ion (a) i n the tissue and the boundary condit ions, causes a s t ra in d i s t r ibu t ion to be set up i n the tissue (b). T h e s ta t is t ical properties of the sonographic tissue mot ion es t imat ion process give rise to the S F . T h e s t ra in d i s t r ibu t ion filtered by the S F is the elastogram (c). T h e second concept of the elastographic process characterizes the formation of the elastogram v i a the use of ul t rasonic pulse-echo techniques. T h i s process is described i n terms of a general theoretical framework, namely the s t ra in filter (SF) , wh ich characterizes the noise properties of the elastographic system [51,106,123,125]. T h e S F is a s ta t is t ical upper bound of the transfer characteristic that describes the relat ionship between actual tissue strains and the corresponding s t rain estimates corrupted by noise effects. It describes the fil tering process i n the s t ra in domain , wh ich allows qual i ty elastographic depic t ion of on ly a l imi t ed range of strains i n tissue. The S F provides upper and prac t ica l performance bounds for s t ra in est imat ion and their dependence on the various system parameters. T h e S F is formed by p lo t t ing the upper bound of the elastographic signal-to-noise ra t io (SNRe) [106,123] for each est imated value of the tissue s t ra in at a given ax ia l and lateral resolut ion (RaXiai,Riaterai) [15,96,99]. T h e S F allows only a restr icted range of s t ra in values to be inc luded i n the elastogram and has a band pass characteristic. A sample o f the S F is shown i n F igure 1.5. The sensi t ivi ty (smallest s t ra in that can be est imated 1.6 Comparison Metrics in Elastography 16 correctly) of the system is l im i t ed on the low end by sonographic noise effects. T h e saturat ion, (largest s t ra in that can be est imated correctly) of the system is l i m i t e d on the h igh end by signal decorrelat ion, resul t ing i n a band pass filter i n the s t ra in domain w i t h a specific dynamic range ( D R e ) . T h e devia t ion of the S F from an idea l all-pass characteristic i n the s t ra in doma in is due to the u l t rasound system parameters, the finite value of the sonographic S N R ( S N R s ) and the effects of s ignal decorrelation. T h e S F can be modif ied and improved u p o n w i t h i n cer ta in l imi t s by proper selection of the ul trasonic system parameters and by algori thms used for s ignal processing (mot ion t rack ing and s t rain est imat ion) . T h e trade offs between achievable S N R e and the Raxiai i n elastography has been explored i n [111]. .40 35 30 25 a? i :g 20 CO 15 10 5. 0: . , . „ . , 10 10 10' 10" Applied Strain (%! Figure 1.5: A sample of s t rain filter is shown. For clarif icat ion dynamic range, sensi t ivi ty and sa turat ion have also been shown on the figure. F ina l l y , the t h i rd concept encompasses the contrast characteristics of the image produced, the elastogram, whose lesion detectabi l i ty is characterized by the contrast-to-noise ra t io ( C N R e ) parameter [11,126]. Detai ls of the fundamental concepts of elastography are given i n the fol lowing chapters. Based on the comparison metrics that were in t roduced i n the S F ( S N R e , D R e , sensi t ivi ty and n n r r Ill' • • 1.7 Performance Optimization in Elastography 17 saturat ion) , the performance of different s t ra in es t imat ion algori thms have been studied. Gradien t -based s t ra in estimators that use R F da ta generally have the advantage of be ing h ighly precise and sensitive and are suitable for es t imat ing smal l s t ra in w i t h h igh accuracy [60,123]. However the s t ra in filter (SF) s tudy has shown that they are not very robust i n the presence of even a moderate amount of deccorrelat ion between the pre- and post-compression signals that most ly happens at h igh compression ratios [52,120]. Therefore, these methods mos t ly do not have a large D R e and fai l to estimate h igh strains. In contrast, direct s t ra in es t imat ion methods such as adaptive strain estimation, spectral strain estimators [4,59,60,121], and gradient-based s t ra in estimators that use the envelope of the signals [59,119,120], are not sensitive to smal l strains. However the S F s tudy has shown that they are robust i n the presence of the noise and decorrelat ion and mos t ly have a larger D R e and are capable of es t imat ing very large strains. 1.7 Performance Optimization in Elastography Signa l decorrelation due to tissue compression is a significant source of error i n tissue displace-ment estimates. T h i s reduct ion i n the value of the correlat ion can have a significant impact on the performance of the s t ra in est imator [123,127]. S imi la r to I D companding w h i c h is a noise-reduction technique (that applies single compression at the t ransmit ter and complementary expansion at the receiver) tempora l s tretching is wide ly used i n s t ra in imaging w i t h a static deformation to improve the S N R e and D R e [3,78,109,124]. Af te r apply ing a known compression to the tissue (1%), the post-compression signal is stretching g lobal ly (1%) to remove the mechanical artifacts. Tempora l s tretching was extended i n 2D [19] and 3D [20] mot ion t racking to remove the decorelation noise due to la teral and elevational motions [51]. Af ter the correction process, ax i a l s t ra in es t imat ion is appl ied to the signal . Loga r i t hmic ampl i tude compression is another technique wh ich is used wide ly [18, 90] that achieves the desired reduct ion of the var iab i l i ty of the R F echo signals wi thout sacrificing the reso-lu t ion and reduces the uncertainty i n the t ime shift es t imat ion. In order to implement logar i thmic compression the whole R F A- l ines needs to be compressed which is computa t iona l ly intensive. Mul t i compress ion is another performance op t imiza t ion technique that suggests d iv id ing a b ig 1.7 Performance Optimization in Elastography 18 static compression to several smal l compressions [118]. Ave rag ing mul t ip le displacement or s t ra in estimates reduces the variance by a factor of N (where N is the number of estimates) and improves the S N R e by a factor of y/N [118]. Since the s t ra in between the R F da ta can be kept w i t h i n the D R e , i f the mul t icompress ion technique is used, the S F w i l l have the shape of a h igh pass filter instead of band pass [109,118], wh ich is close to idea l s t ra in estimator. A s imilar me thod was in t roduced i n [57] tha t selects the op t ima l s t ra in corresponding to each compression level according to its value i n the S F . H i g h s t ra in regions are selected from smal l compressions and low s t ra in regions are picked from h igh compression regions and a final s t ra in image is then generated. S i m i l a r l y to this idea i n real t ime elastography, tempora l f i l tering can be used to generate the elastogram from previously generated s t ra in images [88]. A s mentioned before, t yp ica l ly i n elastography, t ime-domain techniques are used that involve the computa t ion of the time-shift to estimate the displacement following an appl ied compression, and the es t imat ion of the s t ra in by apply ing gradient operations on the previously obta ined t ime-shifts. Since the gradient operat ion introduces addi t iona l ampl i f ica t ion of the noise into the s t ra in est imat ion process the Least Squares S t ra in Es t i ma t o r ( L S E ) was introduced i n [49]. T h e y proposed to use a piecewise approach i n order to reduce the variance of the measured s t ra in profile. T h e improvement i n the S N R e was also shown i n the same paper. Post processing of the s t ra in images makes i t possible to improve the qual i ty of elastograms even further. T h e median filter is normal ly used to reduce noise i n an image, somewhat l ike the mean filter [72,86]. However, i t often does a better j ob t han the mean filter of preserving useful detai l i n the image. M e d i a n filters w i t h I D and 2 D kernels are wide ly used to remove the noisy regions i n elastograms [109]. 1.8 Real-time Strain Imaging and Its Computational Issues 19 1.8 Real-time Strain Imaging and Its Computational Issues M o s t s t ra in es t imat ion methods are used for off-line processing and are t yp i ca l l y compared according to their signal-to-noise rat io ( S N R e ) , contrast-to-noise rat io ( C N R e ) , d y n a m i c range ( D R e ) , sensi t ivi ty and resolut ion of elastograms [96,99,111,123,126,127] . For implementa t ions that preserve the real t ime aspect of u l t rasound imaging, the computa t iona l efficiency of such algori thms is also impor tan t . However, very l i t t le quant i ta t ive informat ion is available on t i m i n g and computa t iona l efficiency issues i n elastography. A l t h o u g h the t e rm time efficient has been wide ly used recently [89,103,109], a formal definit ion has not been used i n the l i terature to compare different methods. A lot of work has been done to have a real-time elastography system that has a h igh c l in i ca l u t i l i t y and to preserve the real-time aspect of u l t rasound imaging. In order to do h igh frame rate (> 20 fps [88]) processing, the computa t iona l cost of the signal processing a lgor i thms becomes an impor tan t factor. Fur thermore the amounts of s t ra in between the sequences of R F da t a are very smal l (< 1%) w i t h negligible decorrelat ion noise [87] therefore a sensitive s t ra in est imator is preferred. It should be noted that smal l strains lead to low S N R e s t ra in images [123], however, mul t ip le successive frames of s t ra in images can also t emporar i ly be filtered for an improvement i n S N R e [87]. Di rec t s t ra in estimators are most ly computa t iona l ly intensive and are not sensitive to smal l strains. Therefore for now, they do not have a potent ia l for real-time applicat ions. T h e current implementa t ion of the adaptive s t rain est imator takes an hour to generate a single elastogram. Gradient-based s t ra in estimators are better suited for real-t ime purposes due to the i r sensi t ivi ty and speed. Since es t imat ing the s t ra in from the displacement estimates has a negligible cost (using a simple differentiator or L S E ) , the mot ion es t imat ion process becomes the most impor tan t step i n real-time elastography. In order to speed up the mot ion t racking process, a number of hardware architectures have been in t roduced i n the l i terature [44,53,92,101]. U s i n g just one digi t of the R F da ta instead of a l l digits to estimate the mot ion was another approach that was in t roduced i n [17]. T h e y showed that i f just one bit of the R F da ta is considered (sign function i n match) it is s t i l l possible to estimate the 1.8 Real-time Strain Imaging and Its Computational Issues 20 mot ion . T h i s way the summat ion and mul t ip l i ca t ion operator can be replaced w i t h s imple X O R and A N D gates i n hardware. Based on this theory, a number of real-time correla t ion estimators architecture have been in t roduced [136]. Fur ther research showed that one bi t resolut ion for m o t i o n es t imat ion is not enough for accurate results [78-80,82]. Therefore, these architectures are not used anymore for elastography. A number of researchers have c la imed methods that are capable of real- t ime m o t i o n t rack ing , such as the Zero Cross Track ing ( Z C T ) me thod [109], Peak Search A l g o r i t h m s ( P S A ) [34]. However, implementat ions have not been reported yet. 1.8.1 Cur ren t Rea l -T ime M o t i o n Tracking Algor i thms Research on accurate real t ime s t ra in imag ing and mot ion est imat ion a lgor i thms has focused m a i n l y on phase domain algori thms tha t have been shown to be capable of r a p i d de terminat ion of displacement. A s mentioned before due to al iasing such methods fai l when a single-step displace-ment exceeds a quarter of the ul trasonic wavelength. To avoid aliasing, phase unwrap methods have been in t roduced to extend the aliasing l i m i t f rom A / 4 to A / 2 [77]. These include Phase R o o t Seeking ( P R S ) [90] that uses Newton ' s i te ra t ion to unwrap the phase-shift and C o m b i n e d Au toco r re l a t i on M e t h o d ( C A M ) [102,103] that uses the envelope correlat ion coefficient to avoid al iasing. P R S was used by the L P - I T company (Lorenz and Pesavento Ingenieurburo fur Information-stechnik, B o c h u m , Germany) and was the first a lgor i thm to calculate s t ra in images of significant size i n real t ime [89]. For the world 's first real t ime implementa t ion of s t ra in imaging , wh ich allows s t ra in imaging i n a c l in ica l setting, they won the 2nd award of the Innovat ion R e w a r d Ruhrgebie t i n November 2000. Cur r en t l y they can generate up to 30 frames of s t ra in per second. C A M is another s t ra in imaging method that achieves real t ime performance. T h e C A M me thod is used by H i t a c h i (H i t ach i M e d i c a l Corpora t ion , C h i b a , Japan) and current ly they are able to generate up to 20 frame of s t ra in images per second [104,105]. 1.9 Vibro Elastography 21 1.9 Vibro Elastography Since s t ra in is not an in t r ins ic property of soft tissue and its value at each point and each t ime depends on a number of parameters such as mechanical properties of tissue and amount of compression, a number of different methods have been in t roduced to produce images of more stable properties of tissue. V i b r o Elas tography ( V E ) is a new imaging system that estimates the elastic properties of soft tissue [117]. T h e tissue is v ib ra ted external ly w i t h known frequency components and a h igh frame rate u l t rasound machine is used to image the tissue. A mot ion es t imat ion me thod is then appl ied to the captured R F da ta and est imated displacements are buffered. Af te r this step, the tissue mot ion is analyzed w i t h a Transfer Func t ion . A c c o r d i n g to this method, the tissue dynamics between two locations along the axis of m o t i o n is considered as a linear dynamic system. T h e transfer funct ion between the two locations is obta ined by spectral analysis w i t h the buffered est imated motions used as inputs and outputs . T h e shape of the transfer function can then be analyzed further. T h e stiffness of tissue can be est imated from the magnitudes of the transfer funct ion at low-frequencies. S imi la r to s t ra in imaging, V E uses the informat ion of the displacement estimates and mot ion t racking plays an impor tant rule i n V E , but i n contrast to s t ra in imaging, m o t i o n t racking is not the only computa t iona l ly intensive part i n V E . Transfer function analysis is another t ime consuming part i n the V E and because of these two parts, un t i l now, V E was on ly used off-line and processing was only done on the previously captured R F data. 1.9 Vibro Elastography 22 Pressure Field Hi x 2 H 3 1 x 3 H 4 1 H ' Hn / / > / / / / / / / / / / / / / / / Pressure Field UUUIUUiU IIUlllUIUIl -1= 7 IHf "3 1 "4 Jxn-1 j x n I H n / ) ) ) ) ) ? } ) J ? } } } /' J V > Figure 1.6: Sketch showing the s implif icat ion of a l inear dynamic system into a static one under low frequency excitat ions [117]. 1.9.1 Transfer Func t ion Analys i s In F igure 1.6, the sequence of tissue displacements at regions i and j are Xi(t) and Xj(t) respec-tively. These displacements are the inputs and outputs of a l inear dynamic system shown as Hj i n F igure 1.6. T h e transfer function between element i and element j is calculated w i t h s tandard techniques [67]. F i r s t , the power spectral density PSDXitXi(uj) of element Xi(t) is computed. T h e n the cross spectral density CSDXitXj(u) between elements i and j is computed [117]. Af ter that the complex transfer function is computed as follows [132]: CSDxx.(u) T h e value of the transfer function is then averaged i n the range of frequency that we are inter-1.10 Thesis Objectives 23 ested i n . T h e shape of the transfer function can then be analyzed further. Mechan ica l properties of the tissue such as stiffness, damping and mass can then be est imated from these transfer funct ion values. For example, the stiffness kj of tissue at region j can be est imated from the magni tudes of the transfer function at low-frequencies: 3 H{-H{+1 V ' D y n a m i c elastography methods have been wel l s tudied i n [117]. V E is a model-free method , which uses spectral analysis of the tissue mot ion . It investigates the behavior of the tissue at different frequencies and finds the frequency response of the tissue. B o t h static and dynamic parameters of the system can then be extracted by focusing on different ranges of frequencies. T h i s way V E eliminates the need for solving the inverse p rob lem to identify the parameters of the system. So lv ing the inverse problem is computa t iona l ly intensive and not applicable i n real t ime . Fur thermore, the implementa t ion of the transfer funct ion according to E q u a t i o n (1.3) makes the a lgor i thm very immune to the noise whi le most parameter identif icat ion methods are very sensitive to i t . 1.10 Thesis Objectives T h e required processing components of bo th s t ra in imaging and stiffness imaging is shown i n F igure 1.7. T h e goal of this work is the real t ime implementa t ion of these packages wh ich are used to measure the mechanical properties of soft tissue inc lud ing stat ic and dynamic properties based on ul trasonic measurements of tissue mot ion under free hand deformations or low frequency mechanical v ibrat ions . To have b o t h systems i n real t ime, a l l blocks need to be opt imized . Since the m o t i o n t rack ing block is the most essential part of bo th systems, a fast and accurate mot ion t racking a lgor i thm that can be used for a l l elastography methods that require displacement es t imat ion is provided. To validate the a lgor i thm from its time-efficiency point of view, the computa t iona l cost of the int roduced a lgor i thm needs to be compared w i t h current real t ime algori thms such as P R S [90] and 1.10 Thesis Objectives 24 Strain Imaging Process Mechanical Echo * Mechanical Echo Properties Signal j Properties Signal D a t a Acqu is i t i on Displacement Estimates Distance from US probe (mm) Mot ion T r a c k i n g St ra in Es t imat ion Transfer Functran Estimates 0 Frequency^Hz} Transfer Function Analysis Stiffness Imaging Process Strain Estimates Distance from US probe (mm) Stiffness Estimates 4 Distance from US probe (mm) St i f fness Es t ima t ion Figure 1.7: Sketch showing the s implif icat ion of b o t h s t ra in imaging and stiffness imaging process. C A M [102,103]. In addi t ion to the computa t iona l cost, the mot ion t racking a lgor i thm also needs to be val ida ted according to its accuracy of es t imat ing the mot ion to make sure i t is b o t h fast and accurate. Fo l lowing the mot ion est imat ion block i n the s t ra in est imat ion process and fol lowing the transfer function analysis i n stiffness imaging process, s t ra in and stiffness need to be es t imated (Figure 1.7). Therefore current estimators like the gradient operator [81] and the least squares est imator [49] need to be s tudied. S imi la r to the mot ion t racking block, s t ra in and stiffness es t imat ion methods need to be fast and accurate to be applicable i n real t ime. A s imula t ion environment and an exper imental set-up need to be provided to validate the proposed real t ime imaging systems based on impor tan t metrics, namely signal-to-noise rat io [51, 106,123,125], resolut ion [15,96,99] and contrast-to-noise rat io [11,126]. S imula t ion results w i l l help us to determine the best parameter settings for each block of the system. Fo l lowing that , experiments w i l l val idate the performance of the whole system, main ly its m o t i o n t rack ing block w i t h real da t a when tissue m i m i c k i n g materials w i t h known mechanical properties are external ly 1.11 Thesis Outline 25 deformed and imaged w i t h an u l t rasound system. F i n a l l y bo th s t ra in and stiffness imaging packages need to be implemented and tested for real t ime performance that can be used i n c l in ica l applicat ions. 1.11 Thesis Outline In this chapter, a review of the current elastography approaches and m o t i o n t r ack ing methods i n successive ul t rasound images was given. Chapte r 2 presents the s tandard Time Domain Cross Correlation ( T D E ) technique as a mo t ion es t imat ion method. Af te r discussing the short comings of T D E , a new technique is in t roduced that we ca l l Time Domain Cross Correlation with Prior Estimates ( T D P E ) since the m a i n structure of the a lgor i thm is s imi lar to that of the T D E method but it is faster and more accurate. Chap te r 3 discusses the computa t iona l cost of the mot ion t rack ing a lgor i thm. It starts w i t h an in t roduc t ion of two popular mo t ion t rack ing algori thms that achieve real t ime performance. A comparison key is then in t roduced and is used to compare the in t roduced me thod i n Chapte r 2 w i t h current available real-time methods. In chapter 4, a number of s t ra in es t imat ion methods are in t roduced. T h e chapter starts w i t h a s tandard gradient operator wh ich is the fastest and simplest s t ra in est imator . It then introduces other versions of the differentiator that have smaller error. F ina l l y , the Least Squares method is explained as another opt ion to estimate the s t rain. T h e chapter finishes by compar ing these s t ra in es t imat ion algori thms. Chapte r 5 presents the s imula t ion s tudy of the T D P E a lgor i thm. T h e chapter starts by quan-t i ta t ive evaluat ion of the T D P E method by means of signal-to-noise ra t io (SNRe), dynamic range (DRe), sensitivity, contrast-to-noise rat io (CNRe) and ax ia l resolut ion (Ra) us ing s imulated R F data . Fur thermore the performance of the median filter and least squares s t ra in es t imat ion and their effect on S N R e , C N R e and resolut ion w i t h the available trade offs are investigated. In Chapte r 6, two types of tissue m i m i c k i n g materials (phantoms) are tested to validate the T D P E a lgor i thm w i t h real R F signals from ul t rasound machine. T h e chapter starts w i t h a descrip-1.11 Thesis Outline 26 t ion of the experimental setup. T h e n , the construct ion of two types of tissue m i m i c k i n g materials is detai led. T h e exper imental results are presented next, and remarks on the results conclude the chapter. In Chap te r 7, the current version of bo th s t ra in imaging and stiffness imaging packages are explained. T h e flow chart of b o t h programs is shown and details about each box i n the flow chart are explained. T h e interface of the current version of bo th packages w i t h their features is demonstrated and the same tissue m i m i c k i n g phantoms that were used i n previous chapter to generate the elastograms are used to validate the performance of the software. T h e chapter finishes by showing the current images that axe generated i n real t ime from the tissue m i m i c k i n g materials and real prostate tissue i n vivo when external deformation is appl ied to the phantoms. F i n a l l y Chapte r 8 presents the conclusions of the research together w i t h suggestions for future work. Chapter 2 Motion Estimation In this chapter, first the s tandard Time Domain Cross Correlation Estimator ( T D E ) technique is presented as a mo t ion es t imat ion method. T D E was among the first algori thms used to perform speckle t rack ing and its properties have been wel l established using the s t ra in filter, S N R e , D R e , sensitivity, resolut ion and C N R e measures [96,99,111,123,126,127] . Af ter discussing the short-comings of T D E , a new technique called Time Domain Cross Correlation with Prior Estimates ( T D P E ) is in t roduced. T h e m a i n structure of the a lgor i thm is s imi lar to that of the T D E m e t h o d but T D P E performs much faster and is more accurate. 2.1 Standard Time Domain Cross Correlation Technique T h e s tandard t ime domain cross correlat ion ( T D E ) was the first a lgor i thm used to perform speckle t racking [81]. W i t h this technique smal l displacements between pairs of ul t rasonic i m -ages that are acquired under different ax ia l compressions are determined using a cross-correlation analysis of the corresponding A- l ines of an R F - d a t a set. A- l ines of R F da ta are spli t into a number of overlapping windows w i t h a specific length ( F i g . 2.1). T h e loca l ax ia l displacements are then computed by cross correlat ing two corresponding windows over a predefined search area (F ig . 2.2(b)). A peak i n the correlat ion occurs at the actual displacement. 27 2.1 Standard Time Domain Cross Correlation Technique 28 W Pre- compression RF A-line AW Overlap Post-compression RF A-line Figure 2.1: The TDE method splits the pre-compression RF A-line into a number of overlapping windows and finds their corresponding windows in the post-compression RF A-line by cross corre-lation method. The important parameters are shown in the figure. W is the length of each window and AW is the shift between windows. Some researchers use the overlap between windows instead of A VF as a signal processing parameter. The exact position of the peak is usually calculated by subpixelling techniques to avoid re-sampling [16]. After this step a gradient operator or least squares strain estimation (LSE) algorithm [49] is used to calculate the strain from the displacement estimates (Equation (1.2)). The size of the search region is one of the parameters that needs to be initialized in the TDE algorithm and it depends on the range of possible motion. For small compression ratios the search region can be set to be small since the motions are also small but for higher compression ratios the search region should be selected much larger since the TDE fails to find the motions that happen outside of its search area. The temporal stretching of the post-compression A-lines is typically done to undo the effects of the mechanical compression on the signal [123,124,127]. As explained in the previous chapter, this stretching improves the correlation between the pre- and post-compression A-lines and reduces the 2.1 Standard Time Domain Cross Correlation Technique 29 s t ra in noise which leads to higher S N R e and D R e . 0 10 20 30 40 50 60 -15 -10 -5 0 5 10 15 sample number time-shifts (a) Echo signals (b) Normalized correlation 10 20 30 40 50 60 -20 -10 0 10 20 sample number time-shifts (c) Echo signals (d) Normalized correlation Figu re 2.2: P re - and post-compressed R F da ta (Left) , and the corresponding normal ized cross-correlations (Righ t ) . T h e peak i n the correlat ion function shows the displacement. 2.1.1 P ros and Cons of the Standard M e t h o d T h e properties of T D E have been wel l established i n the l i terature [96,99, 111, 123,126,127]. T D E is easy to implement . Its S t ra in F i l t e r [123] shows that i t is sensitive to smal l s t ra in and has a h igh S N R e for low strains which is what we need i n real t ime applicat ions. For higher strains, because of de-correlation noise, the method often fails to f ind the displacement correctly. Therefore, i f es t imat ing the s t ra in at h igh compression ratios are required T D E w i l l not perform very wel l . 2.1 Standard Time Domain Cross Correlation Technique 80 RF Lines u 0 windows (Depth) (a) Ideal displacement image on a pair of simulated RF frames on a homogeneous tissue. RF Lines u u windows (Depth) (b) Estimated displacement image with standard T D E with cross corre-lation. F igure 2.3: Theore t ica l displacement image on a pair of s imulated R F frames on a homogeneous tissue (Top) and est imated displacement image on a pair of s imulated R F frames when T D E w i t h cross correlat ion is used (Bo t tom) . T h e z-axis shows the ax ia l displacement, the y-axis shows the R F A- l ines and the x-axis shows the windows s tar t ing from the window p rox ima l to the probe to the d is ta l window deep i n the tissue. Large peaks are the result of selecting false peaks. T h e image shows how often does this happen when non-normal ized cross correlat ion is used. 2.1 Standard Time Domain Cross Correlation Technique 31 RF Lines 0 0 windows (Depth) (a) Estimated displacement image with standard T D E when normalized cross correlation. Hh Lines - » windows (Depth) (b) Estimated displacement image with standard T D E with normalized cross correlation following a 2D median filter. F igure 2.4: E s t i m a t e d displacement image on a pair of s imulated R F frames when T D E w i t h normal ized cross correlat ion is used (Top) fol lowing a 2D median fil tering w i t h a kernel size of 3x3 (Bot tom) . T h e remaining false peaks are removed by using a median filter. 2.2 Time Domain Cross Correlation with Prior Estimates 32 T h e dynamic range of 0% to 2% is enough for the h igh frame rates typ ica l ly encountered i n real t ime s train imaging since the s t ra in between consecutive sets of R F da ta is usual ly smal l and as the frame rate goes up the required dynamic range becomes smaller. T h e only reason w h y T D E is not used i n real t ime applicat ions is i ts computa t iona l inefficiency. Indeed, since T D E cross-correlates two signals over a large area to f ind displacements, i t is slow. Searching over a large area also causes another problem that is referred to i n most papers as false peak [130] or false match [85] and that is s imi lar to aliasing i n the t ime domain . Since the R F da ta are essentially sinusoidal i n shape, their cross correlat ion is also s imi lar to a sinusoid. T h u s searching over a large area may cause another large peak to appear i n the correla t ion coefficients (F ig . 2.2(d)). Choos ing the wrong correlat ion peak leads to false match ing wh ich appears i n the displacement vector as a sharp peak and causes the displacement images to become noisy (F ig . 2.3(b)). T h e problem of f inding false peaks can be al leviated somewhat i f the normal ized correlat ion [9,127,130] is used rather t han the cross correlat ion itself (F ig . 2.4(a)). However, the normal ized cross correlat ion does not e l iminate a l l false peaks and makes the a lgor i thm slower. For further improvements i n displacement es t imat ion, usual ly the normal ized correlat ion is followed by a I D or 2D median fi l tering [23] operat ion i n order to remove the remain ing false matches and noisy regions [96,109]. T h i s improves the performance of the me thod but adds more overhead and makes the a lgor i thm even slower (F ig . 2.4(b)). 2.2 Time Domain Cross Correlation with Prior Estimates A l t h o u g h T D E is easy to implement and very accurate when the cross-correlation funct ion has a single well-determined peak, i t is not suitable for real t ime applicat ions since i t is t ime consuming. T h e inefficiency of the a lgor i thm and the problem of f inding false peaks are T D E problems that are related to each other. B o t h result because the cross correlat ion is performed over a large interval . Since correlat ion is just loca l informat ion, i t is not always sufficient to f ind displacements accurately. Therefore T D E is slow and error-prone. A new mot ion t racking a lgor i thm is in t roduced i n this chapter. W e ca l l i t Time Domain Cross Correlation with Prior Estimates or T D P E since the ma in structure of the a lgor i thm is s imi lar 2.2 Time Domain Cross Correlation with Prior Estimates 3:5 to that of the s tandard T D E method. T D E uses previous est imated da ta to predict the mo t ion at a current posi t ion. T h i s predic t ion process makes the a lgor i thm much faster, and since a l l the predict ions are based on boundary condit ions and internal constraints, the est imated motions are more accurate. T h e T D P E a lgor i thm exploits the fact that neighboring windows i n an R F frame correspond to smal l regions that are close to each other (F ig . 2.5) so their displacements should also be close to each other. Th i s property can be seen i n a l l displacement estimate figures (F ig . 2.3(a) 2.3(b) 2.4(a) 2.4(b)). F igu re 2.5: In comput ing the displacement of the center window, T D P E exploits the fact that displacements of two of its neighbors (leading window to the left and upper above) have already been calculated. T D P E is s imi lar to T D E , but i t collects information about the displacements of neighboring blocks (global information) and uses i t to predict and reduce search areas to smal l regions. In this way T D P E becomes much faster t han T D E because i t searches a very smal l region. It also becomes more accurate and wi thout false peaks because, i n the smal l region employed, the correlat ion coefficient has a single peak, wh ich is easy to detect (F ig . 2.6). 2.2 Time Domain Cross Correlation with Prior Estimates 34 time-shifts window (Depth) Figure 2.6: A sample of the correlat ion function used by T D P E is shown on the left. T h e corre la t ion is not normal ized and i t has on ly been calculated for a smal l number of time-shifts. A sample of est imated displacement vectors generated by T D E and T D P E is shown on the right. E v e n though T D E uses normal ized correlat ion, i t may s t i l l f ind false peaks whi le T D P E can work even w i t h non-normal ized correlat ion and estimates the mot ion correctly. 2.2.1 The M a i n Structure of the A l g o r i t h m Simi l a r l y to the T D E method, T D P E splits each R F line into a number of windows. E a c h window has four adjacent neighbors - a leading window to the left and a lagging window to the right on the same R F l ine, and windows at the same pos i t ion on the upper and lower R F lines ( F i g . 2.5). In elastography, windows have large overlap w i t h each other (50% to 75%) and the spaces between the windows are very smal l . Because of this phys ica l proximity , adjacent windows have s imilar displacements and strains. W e use the displacement estimates of the leading window (previous window) and windows at the same posi t ion on the upper R F line (upper window) wh ich already have been calculated to bracket the displacement of the current window (internal condit ions) . T h e cross-correlation is then used to determine displacements w i t h i n the bracket. T y p i c a l l y , on ly three lags inside the predic ted brackets are checked to detect the displacement qu ick ly and w i t h h igh accuracy. Fur thermore , since the region is smal l , normal iza t ion is not necessary. 2.2 Time Domain Cross Correlation with Prior Estimates 35 F igure 2.7: E s t i m a t e d displacement image on a pair of s imulated R F frames when T D P E w i t h non-normal ized cross correlat ion is used (median filter is not applied). T h e image shows that even though T D P E uses the non-normal ized correlat ion i t s t i l l finds the correct displacements. 2.2.2 P red ic t ion A l g o r i t h m T h e predic t ion a lgor i thm is based on two types of information: (i) a boundary condi t ion at the probe where the displacement is almost zero since the es t imat ion is always done w i t h respect to the probe itself and, (ii) in ternal condit ions that state that the displacement of adjacent neighbors are always close to each other due to their physical proximity . F i n d i n g the displacement is done from the first window (the closest window to the probe) to the last window at the end of the R F line. R F lines are processed s tar t ing from the top. T h u s i n F i g . 2.5 R F windows are processed from the top left to the bo t tom right. Cons ider ing the boundary condit ions, at the first window of each R F line the correlations for the lag -1, lag 0 and lag +1 are calculated and saved. A three-point parabol ic in terpola t ion [16] is then appl ied to f ind the exact pos i t ion of the correlat ion peak. T h i s value is then recorded as the displacement of the current window. T h e closest integer value to this displacement is then used to bracket the search area for the next window (internal condi t ion) . 2.2 Time Domain Cross Correlation with Prior Estimates 36 For the first R F line, since the informat ion on the upper window does not exist, the search area w i l l be: [Round(previousWindowDis.) — 1, Round(previousWindowDis.) + 1], (2-1) wh ich are just three lags and for the rest of R F lines the search area w i l l be: \min{Round{previousWindowDis.),Round{upperWindowDis.)) — 1, (2-2) max(Round(previousWindowDis.), Round(upperWindowDis.)) + 1], (2-3) where previousWindowDis. is the displacement estimate from the previous window (lagging window on the same R F A- l ine ) and upperWindowDis. is the displacement estimate from the upper window (window at the same posi t ion on the preceding R F A- l ine ) (Figure 2.5). Because of physical proximity , the previous-window and upper-window typ ica l ly point to the same posi t ion. T h i s means that i n most cases a bracket of 3 lags is enough for the cross-correlation search area. In the cases i n wh ich the previous-window and upper-window point to different loca-t ions, more than 3 lags are used to bracket the displacement and the one that points to the highest correlat ion coefficient w i l l be selected. A g a i n a three-point parabol ic in terpola t ion is then appl ied to find the exact posi t ion of the correla t ion peak. B y using the informat ion of the upper window we make the predic t ion more precise but we also r u n the risk of car ry ing through an erroneous estimate from one R F line to the next. Therefore, as another implementa t ion version, E q u a t i o n (2.1) can be used for a l l R F lines. T h i s means only the informat ion form previous window w i l l be considered i n the predic t ion a lgor i thm and the upper window w i l l be ignored. T h i s way we may lose track of the mot ion for a single R F line but i t does not affect the rest of the R F lines. Fur thermore a single corrupted l ine can be filtered later w i t h the a id of median filter whi le a corrupted region can not be removed easily. A s another extension to the T D P E method the m a x i m u m correlat ion coefficient at each search can be used to validate the accuracy of the present mo t ion es t imat ion at that window. If this correlat ion coefficient is smaller t han what we expect the search bracket can be enlarged. For 2.3 Conclusion 3 7 example 5 lags can be used instead of 3 lags. It should be noted that i t is because of physical p rox imi ty fact that median filter is used i n the s tandard T D E me thod to remove the remaining false peaks, since i t replaces the value at each p ixe l by the median value of i ts neighborhood. B u t this fact is most ly used as a post processing key and i t has never been used inside the process of m o t i o n es t imat ion itself. T D P E uses th is fact inside the search process to speed up the search and to make i t more accurate. 2.2.3 Pros and Cons of the A l g o r i t h m T D P E is very fast since i t on ly searches very smal l regions to find the displacements. T D P E ca l -culates the correlat ion values at 3 lags and since a l l subpixe l l ing methods need at least 3 correlat ion estimates to find the exact pos i t ion of the correlat ion peak (see A p p e n d i x B . 2 ) , i ts computa t iona l cost is close to m i n i m a l among the correlat ion based a lgor i thm. A n d because i ts m a i n structure is the same as the s tandard T D E method, i t keeps its sensi t iv i ty and h igh S N R e at low compression ratios and i t is also easy to implement wi thout add i t iona l hardware overhead. It is because of these properties that the a lgor i thm offers an at tract ive op t ion for real t ime applicat ions since these applicat ions need a fast and sensitive a lgor i thm which is also easy to implement . T h e current version of the T D P E a lgor i thm predicts the mo t ion at each window based on i ts neighbor's previously est imated mot ion . Therefore, i f for any reason (like R F da ta cor rupt ion or high compression ratios) mo t ion estimates from neighboring windows become imprecise the predict ion a lgor i thm w i l l fa i l . In other words, the method relies on o ld informat ion and i f this informat ion is incorrect the a lgor i thm can not predict and estimate the mot ion correctly. Therefore, corrupted regions start to appear i n displacement image. 2.3 Conclusion In this chapter the s tandard T i m e D o m a i n Cross Cor re l a t ion ( T D E ) was overviewed and its shortcomings were pointed out. Fo l lowing that , the T D P E , as a new mot ion es t imat ion method 2.3 Conclusion 38 *• * « >• Search R e g i o n Search R e g i o n (a) The TDE algorithm. Search-Region Search R e g i o n (b) The TDPE algorithm. Figure 2.8: T h e search region of the T D E a lgor i thm needs to be greater t han the largest possible mo t ion (a), whi le the search region of the T D P E a lgor i thm is fixed (b). 2.3 Conclusion 39 was in t roduced that works s imi la r ly to the T D E but performs much faster and is more accurate. These properties make T D P E an excellent choice for real t ime mot ion es t imat ion w i t h u l t rasound images. Chapter 3 Computational Cost T h e chapter starts w i t h the in t roduc t ion of the computa t iona l cost as a new s t ra in imaging factor. Fo l lowing that it introduces a comparison factor and uses that to compare the current available algori thms that achieve real-time performance such as Phase Root Seeking ( P R S ) [90] and Combined Auto Correlation Method C A M [102] w i t h the T D P E . The computa t iona l cost of s tandard T D E method is also computed as a reference and is compared w i t h these methods. 3.1 Preliminary Considerations and Comparison Measures P R S and C A M are the only mo t ion t racking algori thms that are used i n a real-t ime elastography system. P R S is used by the L P - I T company (Lorenz and Pesavento Ingenieurburo fur Informat ion-stechnik, B o c h u m , Germany) company and was the first a lgor i thm to calculate s t ra in images of significant size i n real t ime [89]. Cur ren t ly they can generate up to 30 frames of s t ra in per second. T h e C A M method is used by H i t a c h i (Hi tachi M e d i c a l Corpora t ion , C h i b a , Japan) and current ly they are able to generate up to 20 frames of s t rain images per second [104,105]. Therefore i n order to validate the in t roduced T D P E a lgor i thm as a real-t ime mot ion t racking a lgor i thm, i t is better to compare i t w i t h these algori thms. For the algori thms considered here and for the s tandard T D E method, the computa t ion t ime 40 3.2 Computational Cost of the TDE algorithm 41 depends on d ig i t a l s ignal processing parameters and is independent of the acoustic and mechanical properties of the tissue. It should be noted that th is statement is not t rue for a l l m o t i o n t r ack ing algori thms. A s an example, i n some i terative a lgor i thms that t r y to maximize the correlat ion, the number of i terat ions depends on the amount of loca l compression which is h igh for softer regions and low for harder regions. Therefore the number of i terations at each region may depend on the mechanical properties of that region. C o m p u t a t i o n a l performance is determined by the number of R F lines, number of windows i n each R F l ine and the computa t ion needed to f ind the displacement at a par t icu lar window. L e t the window-length be W , the shift between windows be AW (Figure 2.1) and the number of windows nW. AW can be calculated di rect ly from W, nW and the R F da ta length L (AW% = [L/(nW x W)] x 100). R F da ta length depends on the depth of the tissue D be ing examined, the sampl ing frequency / s and the speed of sound c (L — fs x 2 x D/c). T h e T D E , T D P E , P R S and C A M algori thms are compared based on the number of multiplica-tions and summations they require [multiplies, adds]. Since summat ion and mu l t i p l i ca t i on are the m a i n operations i n hardware, by count ing t h e m a very accurate comparison can be made. 3.2 Computational Cost of the T D E algorithm T h e number of i terations for the cross correlat ion ca lcula t ion (nl) is a predefined variable i n the T D E a lgor i thm. It shows how many correlations are needed at each window to f ind the displacement of the window. T h e exact value of nl depends on the largest possible mot ion . W i t h o u t stretching, the largest displacement usual ly happens at the end of the R F line since displacement accumulates and i t is p ropor t iona l to the amount of compression appl ied to the tissue. T h i s fact can be seen i n a l l displacement estimate images (F ig . 2.3(a) 2.3(b) 2.4(a) 2.4(b)). E s t i m a t i n g the number of iterations for each a lgor i thm makes the compar ison easier. For a depth of 50 m m , c = 1540 m / s , fs = 40 M H z , the length of the R F - d a t a wou ld be more t han 2500 samples. For 1% s t ra in , wi thout tempora l stretching, the worst case speckle m o t i o n is almost 25 samples (2500 x 0.01). Since i n real t ime the external deformation could be bo th compression and re laxat ion the d i rec t ion of the mot ion is also unknown therefore the search area should be at least 3.3 Computational Cost of the TDPE algorithm 42 2 x 25 = 50 samples around its center point for each window. I n off-line processing since the deformation is mos t ly i n compression the search area can be reduced to 25. Therefore 25 is selected as a default size for search area for 1% compression. For higher compression ratios the search area is t yp ica l ly much larger. It should be noted that for the same system, i f the compression is increased to 10%, the search area for m o t i o n es t imat ion should be at least 250 samples i n each di rect ion (2500 x 0.1) wh ich is so large that i t makes the a lgor i thm very slow and error prone. In pract ice the size of the search region for T D E is a predefined value and set to be a constant wh ich defines the m a x i m u m displacement that can be captured correct ly by the a lgor i thm and the a lgor i thm s imp ly fails for larger displacements. Cons ider ing this the computa t iona l cost CTDE of the T D E for ca lcula t ing the displacement vector for a single R F line is given by E q u a t i o n (3.1) below: where nW is the number of windows, nl is the number of i terations and NCC is the normal ized correlat ion cost. CTDE is equal to the number of windows times the computa t iona l cost at each window, wh ich itself is equal to the number of correlations that need to be est imated at each window. S imi la r to T D E , i n the T D E a lgor i thm A- l ines of R F da ta are split in to a number of overlapping windows w i t h a specific length. The loca l ax ia l displacements are then computed by cross correlat ing two corresponding windows but i n contrast w i t h T D E it typ ica l ly checks jus t 3 lags to detect the displacement. Therefore its computa t iona l cost CTDPE c a n be- described s imi l a r ly as: CTDE = nWx (nl « 25) x NCC, (3.1) 3.3 Computational Cost of the TDPE algorithm CTPDE = nW x (nl « 3) x CC, (3.2) 3.4 Computational Cost of the PRS algorithm 43 where CC is the correlat ion cost (without normal iz ing) . If T D P E w i t h normal ized correlat ion is needed its cost should be considered i n the amount of computa t ion at each window which means CC should be replaced w i t h NCC wh i ch is the normal ized correla t ion cost. 3.4 Computational Cost of the P R S algorithm T h e P R S a lgor i thm has addi t iona l overhead i n order to convert the R F signals to base-band signals. It estimates the t ime shifts from phase shifts of the base-band signal by using a modified N e w t o n i tera t ion (please see A p p e n d i x ( A - l ) ) . T h e a lgor i thm iterates at least twice [90]. P R S needs to consider the subsample time-shifts of base-band signals. Therefore, at each i tera t ion step i t needs to resample the complex da ta of the window i n post-compression R F da ta and calculate i ts complex correlat ion w i t h the window i n pre-compression R F da ta i n order to estimate the phase shift between them. A s a result the computa t iona l cost CPRS of the P R S a lgor i thm can be wr i t t en as: CPRS = BCC + nW x (nl » 2) x (RCcompiex + CC^^), (3.3) where BCC is the base-band conversion cost for a whole R F line, and RCcompiex and CCcompiex are resampling cost and cross correlat ion costs for a complex s ignal , respectively. 3.5 Computational Cost of the C A M algorithm Simi la r to P R S , the C A M algor i thm has an addi t ional overhead i n order to convert the R F signals to analyt ic signals. It estimates t ime shifts from the phase shifts of analy t ic signals that are est imated from the correlat ion of two complex signals. T o unwrap the phase shift, the C A M a lgor i thm estimates the phase shift at several locations (M + N + 1 locations according to E q u a -t ion ( A - 8 ) i n the A p p e n d i x ) w i t h the spacing of A / 2 where A is the wave length of the ul t rasound signal . In paral le l w i t h this process i t also calculates the envelope correla t ion at the same locations. 3.6 Comparison of the Computational Costs 44 T h e phase shift corresponding to the m a x i m u m envelope correla t ion is unwrapped. T h i s phase shift is then used to estimate the displacement [102]. ./V and M are predefined variables i n C A M and they show how many correlations need to be calculated w i t h respect to the current posi t ion to the left ( M ) and to the right (TV) to extend the al iasing l i m i t from A / 4 i n bo th directions. T h e exact value of AT and M depends on the largest possible mot ion . S imi la r to the T D E a lgor i thm, i t is preferred to estimate the number of i terations for the C A M a lgor i thm to make the comparison easier. For c = 1540 m / s , f0 = 5 M H z and A = 0.31 m m , i f we want to capture the m o t i o n due to 1% compression (0.5 m m for the 50 m m phan tom depth) , the a lgor i thm should be able to unwrap the al iasing up to 3 A / 2 (0.46 mm) i n bo th direct ions (compression and re laxat ion) . T h i s means M and TV i n the C A M algor i thm should be at least set to 3, since increasing M and TV increases the unwrapped region by A / 2 . T h i s means that the autocorrelat ion function needs to be run TW + TV + l = 3 + 3 + l = 7 t imes. Cons ider ing that the deformation is just compression just one of M or TV should be set to 3 and the other one can be set zero and the ca lcula t ion of autocorrelat ion function reduces to 4 t imes ( M + TV + 1 = 3 + 0 + 1 = 4 ) . T h u s M + TV + l = 4 i s a n est imated value for the number of i terat ions for the C A M a lgor i thm for 1% compression, whi le for higher compression ratios M and TV should be considered much larger. A s a result for the C A M a lgor i thm the computa t iona l cost CCAM can be wr i t t en as: CCAM = ACC + nWx (nl « 4) x (CCCMNPLEX + ENC), (3.4) where ACC is the cost of convert ing to analyt ic signal for a complete R F line and ENC is the envelope normal ized correlat ion cost. 3.6 Comparison of the Computational Costs In order to compare the computa t iona l costs of these methods we need to go into further deta i l . Since windows have overlap, a more efficient conversion to base-band and analyt ic signals is done 3.6 Comparison of the Computational Costs 45 once per l ine before processing as shown i n E q u a t i o n (3.3) and E q u a t i o n (3.4). In order to consider the cost of conversion for each window ind iv idua l l y rather t han the whole R F line we assume that windows have 50% overlap w i t h each other. So the cost of conversion is d iv ided between neighboring windows. Therefore E q u a t i o n (3.3) and E q u a t i o n (3.4) can be revised as follows: CPRS = nWx [BCCW/2 + (nl « 2) x (RCcomplex + CCcompiex)}, (3.5) CCAM = nWx [ACCW/2 + (nl « 4) x (CCcompiex + ENC)}, (3.6) where BCCW is the cost of base-band conversion of a single window, ACCW is the cost of analyt ic conversion of a single window. T h e cost of a correla t ion is that of inner product of vectors (Equa t ion (B-2)) . For the normal ized correlat ion one needs to calculate the autocorrelat ion of the two windows as wel l , thus m a k i n g i t necessary to calculate 3 inner products of vectors (Equa t ion (B-3)) . Fur thermore l inear resampl ing calculations consist of 2 complex mul t ip l ica t ions and one complex summat ion [89, 90] for each sample. Moreover complex mul t ip l i ca t ion needs 4 mul t ip l ica t ions and 2 summations and complex summat ion needs 2 summat ions for each element: (a,jb) x (c,jd) = (a x c — b x d, j (a x b + c x d)), (3.7) (a,jb) + (c,jd) = (a + c,j (b + d)), \(a,jb)\ = y/(axa + bxb), Therefore the cost of each operat ion according to the [multiplies, adds] on a vector w i t h length N can be wr i t t en as follow: • Element-wise M u l t i p l i c a t i o n of two real signals: [TV, 0] • Inner P r o d u c t of two real signals: [N, N — 1] • Element-wise M u l t i p l i c a t i o n of two complex signals: [4N, 2N] 3.6 Comparison of the Computational Costs 46 • Inner P r o d u c t of two complex signals: [47V, 47V — 2] • L inea r Resampl ing of two real signals: [27V, TV] • L inea r Resampl ing of two complex signals: [8TV, 6TV] • N o r m calcula t ion of a complex signal ( ignoring the cost of sqrt): [27V, TV] Furthermore , as shown i n A p p e n d i x (Equa t ion ( A - 3 ) ) , the conversion from ana ly t ic signals to base-band signals consists on ly of element-wise mul t ip l i ca t ion of two complex signals of the size of each window. B u t the conversion from real signals to analyt ic signals is ca lcula ted by adding an imaginary part to the signal, wh ich is equal to its H i lbe r t t ransform. There are several ways of comput ing the H i l b e r t transform; i n this paper we w i l l use the most pertinent 14-point F I R filter approach presented i n [89]. Ca l cu l a t i ng the 14 point F I R i n t ime domain for each point requires 14 mul t ip l ica t ions and 14 summations. Since b o t h P R S and C A M algori thms use the analy t ic signals this conversion overhead should be added to the overall cost of each a lgor i thm. T h e current implementa t ion of C A M uses a quadrature detector to capture real-time analyt ic signals [102] but since this feature is not available on a l l u l t rasound machines the analyt ic conversion cost is considered here for the C A M a lgor i thm. B u t i f a quadrature detector is already available on the ul t rasound machine the cost of conver t ing the real s ignal to the analy t ic signal should be omi t t ed from the overall cost for b o t h C A M and P R S algori thms. Cons ider ing these assumptions the computa t iona l requirement of each a lgor i thm, based on [multiplies, adds] reduces to: 3.6 Comparison ofthe Computational Costs 47 SearchRegion NormalizedCorrelation CTDE = nWx ( n / « 2 5 ) x 3 x [Wlength, Wlength - 1], (3.8) SearchRegion Correlation CTDPE = nW x (nl « 3) x [Wlength, Wlength - if, .Ana(!2/£ic-r-i?ase—bandConversion , * s Cpf ts = nW x [[14 + 4 ,14 + 2] x Wlength/2 +(nl « 2) x Resampling ComplexCorrelation , * v , * •> ([8,6] x Wlength + [4 x Wlength, 4 x Wlength - 2])], .47ia /y tocCora;ers ion / A V CCAM = nW X [[14,14] x Wlength/2 +(nl « 4 ) x ComplexCorrelation EnvelopeN or malizedC or relation ([4 x Wlength, 4 x Wlength - 2] + 3 x [2Wlength, Wlength] )], N o w i t is easy to compare the computa t iona l efficiency of these algori thms. Since nW (number of windows i n each line) and Wlength ( length of each window) are a l l variables common to a l l methods, their product can be replaced by a constant r nW x Wlength = r , (3.9) Therefore the simplest form for the number of operations [multiplies, adds] for each a lgor i thm reduces to: CTDE = [75 x T , 75 x T ] , (3.10) CTDPE = [3 x T, 3 x r ] , C P R S = [33 x T,28 x T] , C = [ 2 6 x r , 2 1 x r ] , CCAM = [47 X T , 3 5 X T ] , C ^ M = [ 4 0 x r , 2 8 x r ] 3.7 Conclusion 48 where CQAM and CpRS are the computa t iona l costs of the C A M and the P R S algori thms when a quadrature detector is available on the ul t rasound machine. F r o m the above, the computa t iona l cost of the T D P E a lgor i thm is almost 1/25 of that of T D E and 1/10 of those of C A M and P R S algori thms for sma l l compressions (< 1%). For h igh compression ratios, as ment ioned before, the search area for T D P E and P R S remains the same since the number of i terat ions for bo th of them is constant whi le the search area increases for T D E and C A M wh ich requires t hem to run more calculations. T h e result also shows that number of operations for C A M is almost close to that of T D E whi le i t was reported to be 6 t imes faster [102]. 3.7 Conclusion In th is chapter the computa t iona l cost was in t roduced as a new s t ra in imaging factor. T h e number of mul t ip l ica t ions and summations [multiplies, adds] was used as a comparison key. B y using this factor, the T D P E a lgor i thm (which was in t roduced i n Section 2.2) was compared w i t h current real-t ime algor i thms and i t was shown that i t estimates the mot ion w i t h 1/10 of the com-puta t iona l cost of b o t h Phase Root Seeking [90] and Combined Auto Correlation Method [102] and w i t h 1/25 of the computa t iona l cost of the s tandard T D E method. Chapter 4 Strain Estimation A wel l denned approach does not exist to validate the m o t i o n t rack ing a lgor i thm independent of the strain es t imat ion process i n elastography. However, the compar ison factors for s t ra in es t imat ion such as the strain filter [51,106,123,125], axial resolution and lateral resolution [15,96,99] and contrast-to-noise ratio [11,126] have been wel l specified i n the li teratures. Therefore, th is chapter discusses the es t imat ion of s t ra in from the displacement estimates and the T D P E a lgor i thm w i l l be val idated on the next chapter based on the s t ra in images that i t generates. A s explained i n the In t roduct ion , i n gradient-based s t ra in es t imat ion methods, s t ra in is ca lcu-la ted from the est imated displacement. The chapter starts w i t h the overview of two s tandard s t ra in est imat ion techniques namely, the gradient operator [81] and the least squares strain estimator [49]. Fo l lowing those, new extensions for bo th techniques w i l l be in t roduced, that can be used i n the s t ra in es t imat ion process. It should be noted that i n this chapter strain estimation is the name given to the process of es t imat ing the value of s t ra in from the displacement estimates wh ich have already been est imated, instead of the whole process of es t imat ing s t rain from the sequences of raw R F data . 4.1 Gradient Operator After es t imat ing the displacement at each window the s t ra in needs to be est imated at that 49 4.1 Gradient Operator 50 A. Displacement Estimates . • > Depth F igure 4.1: Sketch showing the process of es t imat ing the s t ra in from the displacement estimates. region. A c c o r d i n g to E q u a t i o n (1.2) loca l s t ra in is equal to the difference between the est imated displacement at the current window and its previous window, normal ized by the space between the windows. T h i s is shown i n E q u a t i o n (4.1): . t o - " " - * ' - ' * . ( « ) where d (i) is the displacement estimate i n the indexed segment and A T is the spacing between windows. To estimate the s t ra in we can s imply differentiate the displacement estimates. Therefore the s t rain es t imat ion problem is s imilar to the slope es t imat ion problem and a l l slope es t imat ion techniques can be appl ied. T y p i c a l l y i n elastography, the s train is est imated by app ly ing a gradient operat ion on the 4.2 Least Squares Strain Estimation (LSE) 51 previously obta ined time-shift estimates. S t r a in es t imat ion by a simple gradient operator is very fast. 4.2 Least Squares Strain Estimation (LSE) Least Squares Strain Estimator ( L S E ) was in t roduced i n [49] to improve the process of s t rain es t imat ion. T h e m a i n goal of L S E is the reduct ion of noise ampli f icat ion due to the gradient operat ion and this is achieved by a piecewise l inear curve fit to the est imated displacement field. L S E finds the best-fi t t ing line to a set of points inside a region around the point of interest (kernel), by m i n i m i z i n g the sum of the squares of the offsets of the points from the fi t ted l ine (Figure 4.2). T h e slope of the l ine is then reported as the estimate value of the gradient of the time-shifts at the center of the kernel . S imi l a r ly to the gradient approach, this value is then normal ized by the space between the windows to estimate the real value of the s t ra in at each window. A displacement vector d w i t h length M, a round a given kernel of N points , m a y be modeled by: Displacement Estimates Fitted Line » Depth F igure 4.2: L S E w i t h a kernel size of 5 is shown. A line is fitted to the da ta points w i t h a least squares error and the value of the slope is then reported as the estimate of the gradient at the center of the kernel. 4.3 Higher Order Numerical Differentiators 52 d(i) = a x z(i) + b, (4.2) where z is the tissue depth and a and b are constants to be est imated. E q u a t i o n (4.2) can be wr i t t en i n ma t r i x format as: (4.3) where A is an N x 2 m a t r i x for wh ich the first co lumn represents the depth pos i t ion d(i) and the second is a co lumn of ones. T h e classical least squares solut ion of E q u a t i o n (4.3) is given by [14]: = [ATA]-1ATd, (4.4) where d contains the noisy displacement data, a and b , are the L S E estimates of a and b and a is an estimate of the gradient at the middle of the kernel. In other words i n L S E method the differentiation of the displacement estimates is calculated by f i t t ing a straight line to a number of neighboring windows that has a least squares error. T h e slope of this line is then reported as loca l differentiation. 4.3 Higher Order Numerical Differentiators T h e gradient operat ion amplifies smal l differences i n displacements and has a l inear frequency response (JOJ). T h i s operat ion can be thought of as the crudest form of a high-pass filter. It has low gains at low frequencies and h igh gains at h igh frequencies and therefore amplifies high frequency noise [60,121]. N u m e r i c a l differentiation methods w i t h smaller error can also be used i n the s t ra in es t imat ion process. These techniques are wide ly used for slope es t imat ion but they have not been used i n elastography yet. N u m e r i c a l differentiation formulas can be derived by (i) construct ing the Lagrange interpolat ing 4.4 LSE with Higher Order Polynomials 53 po lynomia l th rough a number of points, (ii) differentiating the Lagrange po lynomia l , and (iii) f inally evaluat ing the differentiation at the desired point . A three points, five points and seven points rule for ca lcula t ing the differentiation for s t ra in es t imat ion can be wr i t t en as E q u a t i o n (4.5), E q u a t i o n (4.6) and E q u a t i o n (4.7): .w-" ( , + 1 ^"- 1 ) . («) _ - d ( t + 2) + 8d ( t + l ) - 8 d ( i - l ) + d ( t - 2 ) s W ~~ 1 2 A T ' . , d {i + 4) - 40d (i + 2) + 256d (i + 1) - 256d (i - 1) + 40d (i - 2) - d (i - 2) , . S W = 3 6 0 A T ' ( 4 ' ? ) where again d (i) is the displacement estimate i n the indexed segment and A T is the spacing between windows. B y using more displacement estimates from neighboring windows (more points) numer ica l dif-ferentiators can estimate the gradient more accurately. B u t since displacement estimates are not accurate the interpolated p o l y n o m i a l w i l l not be accurate either and therefore increasing the number of points may not necessarily improve the s t rain es t imat ion process. 4.4 LSE with Higher Order Polynomials T h e idea of L S E can be extended to higher order po lynomia l f i t t ing. T h i s means instead of best-f i t t ing line (first order po lynomia l ) we can fit a higher order po lynomia l to da ta points (LSEi(N) where i is the order of p o l y n o m i a l and N is the size of the kernel) to improve the b a n d w i d t h of the L S E . For example for a second order po lynomia l a displacement vector d, a round a given kernel of N points, can be modeled by: d(i) = a x z(i)2 + b x z(i) + c, (4.8) 4.5 Conclusion 54 where z is the tissue depth and a, b and c are constants to be estimated. Equation (4.8) can be written in matrix format as: d = A (4.9) where A is an N x 3 matrix for which the first column represents the squares of the depth d(i)2, the second column represents the depth position d(i) and the third is a column of ones. The classical least squares solution of Equation (4.9) is given by [14]: (4.10) where d contains the noisy displacement data, a, b and c , are the LSE2 estimates of a, b and c and 2 x a x [ t + J^ - 1 J + b is an estimate of the gradient at the middle of the kernel where |_*"2^J is the center of the kernel. Although going to higher order polynomials can help increase the bandwidth of the strain estimator and increase its resolution by providing faster transitions from one region to another, but since the displacement estimates themselves are not very accurate, fitting a higher order polynomials to these inaccurate data is not the best solution. The trade-off between the improvement in signal-to-noise ratio and loss of resolution needs to be studied in order to choose an optimum kernel size and polynomial. 4.5 Conclusion Standard strain estimation techniques have been overviewed in this chapter. Furthermore two slope estimation techniques have also been introduced, that can be used for strain estimation. Numerical differentiation methods have been introduced as an extension to the gradient operator and the least squares strain estimation method was extended to higher order polynomials. 4.5 Conclusion 55 It should be noted that est imating the stiffness i n v ib ro elastography according to the E q u a -t ion (1.4) is also a slope est imat ion process: (4.11) Therefore, a l l s t ra in est imat ion techniques can be used to estimate the stiffness f rom the mag-nitudes of the transfer function at low-frequencies. Chapter 5 Simulation Results T h i s chapter presents a s imula t ion s tudy of the T D P E a lgor i thm. T h e chapter starts by a quant i ta t ive evaluat ion of the T D P E method by means of signal-to-noise ra t io (SNRe), dynamic range (DRe), sensitivity, contrast-to-noise ra t io (CNRe) and ax ia l resolut ion (Ra) using s imulated R F data . Fur thermore , the performance of the least squares s t ra in es t imat ion ( L S E ) is investigated. T h e effect of the L S E on S N R e have been s tudied i n [49] but its effect on the ax ia l resolut ion and C N R e w i t h the available trade-offs have not been studied yet. In add i t ion to that , for the first t ime, the effect of the median filter on the S N R e and ax ia l resolut ion w i t h the available trade-offs is investigated. 5.1 Validation T h e performance of T D P E is val idated using I D s imula t ion i n M a t l a b ( M a t h W o r k s , Inc., N a t -ick, M A , U S A ) v i a three impor tan t metrics, namely (i) the s t ra in filter (SF) [51,123,125] wh ich represents the behavior of S N R e as a function of compression ratios and includes some informat ion about dynamic range and sensi t ivi ty of the a lgor i thm, (ii) the elastographic contrast-to-noise ra t io ( C N R e ) wh ich is propor t ional to the detectabi l i ty of the lesion i n the tissue [11,126,127] and (iii) elastographic ax ia l resolution (Ra) that shows the smallest lesion that can be visual ized i n the 56 5.2 ID Simulation of a 2D Image 57 elastogram [5,96]. 5.2 ID Simulation of a 2D Image I D simulations were preferred over 2D simulat ions for two reasons: 1. la teral mo t ion and beam effects have been accounted for as derat ing factors i n the S F [51,125] and, hence, the results are expected to be s imi lar to those detai led i n those papers, and 2. significantly smaller s imula t ion t ime is required when performing a s ta t is t ical analysis over several parameters, l ike the w i n d o w size, overlap, f i l tering etc. 2 D R F frames were generated w i t h I D simulations for each R F A - l i n e . S imi l a r l y to previous work [117], pre- and post-compression R F signals corresponding to a 40-mm (D) uni formly elastic target segment were generated using a convolut ion between the impulse-response of the system (point spread function, P S F ) and a no rma l d i s t r ibu t ion of scatterer ampli tudes [130]. T h e point scatterers were uni formly d i s t r ibu ted w i t h density of at least 10 scatterers per pulse w i d t h . T h e speed of sound i n tissue was assumed to be constant at 1540 m / s . T h e P S F was s imulated using a Gauss ian-modula ted cosine pulse w i t h a 5 M H z center frequency ( / „ ) . T h e sampl ing frequency (f3) was set to 40 M H z . A l t h o u g h the value of each of these parameters has a direct effect on the performance of the s t ra in est imator [123], these parameters are chosen to simulate the R P 500 ul t rasound system (Ul t rasonix , Burnaby , Canada) as closely as possible. To compute the post compression signal the phantom was modeled as a serial connection of springs [81]. A stiffness map was assigned to each l ine and the displacements of each point scatterers were calculated by compressing the model and ca lcula t ing the displacement at each region. For s imula t ing a homogeneous phan tom the value of the stiffness is set to be the same throughout the image which leads to uni form compression of the point scatterers [15]. For s imula t ing a phan tom w i t h a circle inclus ion, the stiffness of a circle inside a phantom was set to be 4 t imes stiffer t h a n its background and for s imula t ing a three layered phantom, the stiffness of the middle layer was set to be 4 times stiffer t han other two layers. T h e post-compression signal is then calculated by convolving the compressed point scatterers w i t h the or ig inal P S F . 5.2 ID Simulation of a 2D Image 58 150 RF Lines 0 0 Window (Depth) RF Lines 0 0 Window (Depth) (a) Homogeneous Phantom. (b) Phantom with Inclusion. Figure 5.1: E s t i m a t e d displacement at 1% strain on b o t h s imulated homogeneous phan tom and phantom w i t h inc lus ion when T D P E is used as a mot ion t racking a lgor i thm ( D = 40mm, W = l m m , overlap=75%, / c = 5 M H z and / s = 4 0 M H z ) . A sample of the displacement estimate from bo th s imulated homogeneous phantom and s imu-lated phantom w i t h inc lus ion is shown i n F i g . 5.1(a) and 5.1(b). T D P E was used to estimate the mot ion from s imulated R F images. T h e inclusion is almost detectable i n the displacement estimates itself before s t ra in es t imat ion process. 5.2.1 Eva lua t ion of S t ra in Est imators T h e choice of s t ra in est imator has a direct effect on a l l of the metrics. Therefore, before a quant i ta t ive evaluat ion of the T D P E method by means of SNRe, DRe, sensitivity, CNRe and Ra, it is impor tant to evaluate the s t ra in estimator that we are going to use. S tandard s t ra in es t imat ion techniques have been overviewed and new strain es t imat ion tech-niques have also been in t roduced i n the previous chapter. T h i s section evaluates these methods based on their s t ra in images' signal-to-noise rat io. T h e elastographic signal-to-noise rat io characterizes the noise at which a value of s t ra in is est imated and it is denned by [123]: 5.2 ID Simulation of a 2D Image 59 SNRe = — , (5.1) where ms denote the s tat is t ical mean s t ra in estimate and ers denotes the s tandard dev ia t ion of the s t ra in noise es t imated from the elastogram. In order to compare the s train es t imat ion algori thms, a three layered phan tom (4 x 4 c m 2 ) w i t h a four t imes stiffer layer inside was s imulated and two R F images were generated w i t h 100 R F A-l ines each (1% compression). Displacements were est imated at each of the lines by using the T D P E a lgor i thm ( W = l m m , overlap=75%). T h e displacement estimates of the first soft layer were then used to compare the s t ra in est imat ion methods (50 windows) . T h e performance of the numerical differentiators is shown i n Table 5.1. A s expected, the results from the numer ica l differentiators show a l i t t le improvement w i t h respect to the gradient me thod (smaller deviat ions and larger S N R e ) . Strain Estimator Standard Deviation Mean SNR - * T D M f.nn. Ideal 0 1.4134 O O Gradient 0.9283 1.4175 1.4959 3 Points Differentiator 0.7055 1.4239 1.9735 5 Points Differentiator 0.8455 1.4256 1.6504 7 Points Differentiator 0.8755 1.4271 1.5963 Table 5.1: C o m p a r i s o n of gradient operator w i t h numer ica l differentiators. A l l values are averaged over 100 realizat ions. T o compare the method the result from ideal s t ra in est imator is also inc luded . T h e table compares the algori thms ' s tandard devia t ion and their resul t ing signal-to-noise rat ios ( S N R ) . T h e performance of the LSE\ is shown i n Table 5.2. A s the size of the kernel increases the algori thm's performance improves (smaller deviat ions and larger SNRe): F i n a l l y the performance of the LSE2 is shown i n Table 5.3. S imi la r ly to the LSE\, as the size of the kernel increases the algori thm's performance improves (smaller deviat ions and larger SNRe). However, the improvement is less marked as i n the case of the LSE\ method . T h i s could be explained by the fact that i n LSE2 we are f i t t ing a second order po lynomia l (to increase the bandwid th of the estimator) wh ich is not a l ine anymore and have its own r ipple . 5.2 ID Simulation of a 2D Image 60 S t r a i n E s t i m a t o r S t a n d a r d D e v i a t i o n M e a n Gradient 0.9283 1.4175 1.4959 LSEXIA) 0.5489 1.4216 2.5293 L S £ i ( 6 ) 0.3204 1.4224 4.2678 L5£*i(8) 0.1950 1.4232 6.9533 LSEI{IQ) 0.1298 1.4223 10.4761 L S £ i ( 1 2 ) 0.0917 1.4215 14.7107 LSE^U) 0.0673 1.4200 19.9715 LSExO-6) " 0.0520 1.4188 26.0171 L 5 £ i ( 1 8 ) 0.0436 1.4164 31.3489 Table 5.2: Compar i son of LSE\(N). A l l values are averaged over 100 real izat ions. S t r a i n E s t i m a t o r S t a n d a r d D e v i a t i o n M e a n SNR - STD M pn.n Gradient 0.9283 1.4175 1.4959 LSE2{A) 1.2878 1.6617 1.2596 LSE2(6) 0.9718 1.5700 1.5828 LSE2{8) 0.7302 1.4658 1.9201 LSE2(W) 0.5465 1.4287 2.4496 LSE2(12) 0.4012 1.4233 3.2980 LSE2{14) 0.3073 1.4225 4.3012 LSE2(16) 0.2360 1.4198 5.4440 LSE2{18) 0.1860 1.4175 6.8615 Table 5.3: Compar i son of LSE2(N). A l l values are averaged over 100 real izat ions. Fo r the qual i ta t ive comparison of the s t ra in es t imat ion method, a result of m o t i o n es t imat ion of a s imulated three layered phan tom was used as i n i t i a l da ta i n the s t ra in es t imat ion process. T h e result of s t ra in es t imat ion w i t h the gradient operator is shown i n F igu re 5.2(b). It has been shown that wi thou t any filtering this operator results i n a noisy elastogram w i t h h igh resolut ion (sharp edges). T h e result of s t rain es t imat ion w i t h a numerica l differentiator is also shown i n F igure 5.2(c) and 5.2(d). B o t h methods show and improvement w i t h respect to the simple gradient me thod (smaller upper bound) but the difference between three points differentiator and seven points differentiator is not quite clear. F i n a l l y the performance of LSE\ m e t h o d is shown i n for two kernel sizes of 6 (F ig . 5.2(e)) and 12 (F ig . 5.2(f)). T h e improvements i n the elastograms are evident. A s the kernel size increases the r ipple i n the est imated s t rain decreases wh ich shows the improvement i n S N R e . 5.2 ID Simulation of a 2D Image GI Windows (Depth) Windows (Depth) (c) Elastogram with three points differentiator. (d) Elastogram with seven points differentiator. Windows (Depth) Windows (Depth) (e) Elastogram with LSEi with kernel size of 6. (f) Elastogram with LSE\ with kernel size of 12. Figure 5.2: A sample of displacement estimates of a three layered phantom at 1% compression with its corresponding elastogram when different strain estimators are used. 5.2 ID Simulation of a 2D Image 62 To be consistent w i t h other researchers, i n this thesis a s tandard L S E me thod (LSE{) and the gradient operator w i l l be used as s t ra in estimators. T h e trade-offs involved i n their use w i l l be presented i n the fol lowing sections. S tudy ing the trade-offs between the S N R e and resolut ion for higher degree p o l y n o m i a l fittings and LSE-2 w i l l be the topic of future work. 5.2.2 Strain Filter Study (SF) T h e S F [123] describes the impor tan t relat ionship among the resolution, dynamic range ( D R e ) , sensi t ivi ty and elastographic S N R e (defined as the ra t io of the mean s t ra in to the s tandard devia t ion of the s t ra in i n the elastogram), and may be plot ted as a graph of the upper bound of the S N R e vs. the s t ra in experienced by the tissue, for a given elastographic ax ia l resolut ion (as defined by the da ta window length [96]). T h e S F is a s ta t is t ical upper bound of the transfer characterist ic that describes the relat ionship between actual tissue strains and the corresponding s t ra in estimates depicted on the elastogram. It describes the filtering process i n the s t ra in domain , wh ich allows qua l i ty elastographic depic t ion of only a l im i t ed range of strains from tissue. T h i s l imi t ed range of strains is due to the l imi ta t ions of the ul t rasound system and of the s ignal processing parameters and algori thms. T h e S F is obtainable as the ra t io between the mean s t ra in estimate and the appropriate lower bound on its s tandard devia t ion. T h e low-stra in behavior of the S F is determined by the acoustic properties of the u l t rasound machine and the high-strain behavior of the S F is determined by the rate of decorrelat ion of a pair of corresponding signals due to tissue d is tor t ion that arises due to several reasons. These include nonuniform and /o r nonsl ip boundary condit ions, la teral and elevational mot ion , differing angular or ientat ion of the transducers relative to the direct ion of compression (result ing i n depth-dependent decorrelation), and stress-strain nonlinearit ies i n tissues and t i ssue-mimicking phantoms [111]. T h e general appearance of the S F i n three dimensions, demonstra t ing the trade offs between S N R e and s t ra in at a l l resolutions is available i n [79] which investigates the behavior of the S F for different signal processing parameters such as window size (W) and window overlap. In add i t ion to the selection of the signal processing parameters, the S F can be derated by other cor rupt ing processes such as frequency dependent at tenuation [125] and la teral and elevational signal decor-5.2 ID Simulation of a 2D Image 63 0 1 2 3 4 5 6 7 8 9 1 0 11 12 1 3 14 15 S t r a i n % Figure 5.3: T h e s t ra in filter of the T D P E a lgor i thm. Compress ion is done once from 0.001% to 100% i n logar i thmic steps (Top) and once from 0% to 15% i n 0.2% steps to increase the accuracy (Bot tom) . 100 realizations are used to show the mean values (Top, Bo t tom) and the devia t ion from the mean values (Bo t tom) . The gradient operator is used as a s t ra in est imator wi thout any filtering ( D = 4 0 m m , W = l m m , overlap=75%, fQ=5 M H z and / s = 4 0 M H z ) . 5.2 ID Simulation of a 2D Image 64 relat ion [51]. T o generate the S F of the T D P E method and compare it w i t h the S F of the s tandard T D E a lgor i thm, a homogeneous 40 m m phan tom was s imula ted using a I D s imula t ion . T h e phan tom was compressed from 0.001% to 100% i n logar i thmic scale. B o t h T D E and T D P E w i t h W = 1 m m and 75% window overlap (160 window i n 40 m m ) were used to estimate the m o t i o n at each step and the gradient operator was used to estimate the s t rain. T h e S N R e was est imated at each compression ra t io according to E q u a t i o n (5.1) and the S F s are shown i n (F ig . 5.3(a)). T h e mean values of the S N R e for 100 simulations are used as the S N R e at each compression. It can be seen from F igure 5.3(a) that for sma l l compression ratios the T D P E and the T D E performs s imi la r ly (sensitive to 0.01% compression). For higher compression ratios, the T D P E shows a larger dynamic range (DRe) compared to the T D E . T h e S N R e of the T D E a lgor t ihm drops for compressions higher than 3% while the T D P E a lgor i thm starts to saturate at compressions higher t han 10%. T o investigate the behavior of the T D P E i n the h igh S N R e region (0.2% to 15%) another s imula t ion was run . A homogenous 40 m m phan tom was simulated using a I D s imula t ion . T h e phantom was compressed from 0% to 15% i n steps of 0.2% s t ra in (76 frames). T D P E w i t h W = 1 m m and 75% window overlap (160 window i n 40 mm) was used to estimate the m o t i o n at each step and the gradient operator was used to estimate the s t rain. T h e S N R e was est imated at each compression ra t io (F ig . 5.3(b)). The M a t l a b boxplot function was used to show the mean value and the devia t ion from the mean value of the S N R e for 100 simulations. It can be seen from F igure 5.3(b) that T D P E is capable of es t imat ing strains of up to 8% compression. 5.2.2.1 Effect of L S E and M e d i a n F i l t er on the S tra in F i l ter A s discussed i n the previous chapter, the L S E can be appl ied instead of a s imple gradient operator to improve the S N R e wh ich leads to elastograms w i t h higher quali ty. T h e effect of the L S E on the S F of the s tandard T D E a lgor i thm has been presented i n [49]. To study the effect of the L S E on the proposed T D P E a lgor i thm, the improvement i n S F is shown i n F i g 5.4(a) and an improvement i n S N R e at a single compression ratio is shown i n F i g 5.4(b) when the L S E me thod w i t h different kernel sizes (2-12) are used instead of the gradient operator on the same da ta that was used i n 5.2 ID Simulation of a 2D Image 65 F igure 5.3(b). T h e result shows more than linear improvement i n S N R e when L S E is used w h i c h confirms the results f rom [49] that state that a kernel size of N improves the S N R e by TV3/2. T o show this , the improvement i n S N R e at a single compression ra t io (1%) is shown i n F i g 5.4(b). L S E w i t h a kernel size of 2 is the same as the gradient operator and as the kernel size increases, the S N R e improves. M e d i a n filtering is another popular f i l tering type i n elastography that helps to remove the noisy regions and improves the qual i ty of the elastograms. M e d i a n f i l tering can be appl ied either after m o t i o n es t imat ion to remove false peaks [130] or after s t rain es t imat ion to remove noisy regions [109]. A l t h o u g h median fil tering has been wide ly used i n elastography its effect on the S F and resolu-t ion has not been investigated i n the l i terature yet [111]. T o demonstrate the performance of median fi l tering, the s imulated R F da t a were used by T D P E to estimate the motions. 2D median filters w i t h variable kernel sizes were used once to filter the displacement estimates before apply ing the s t ra in estimator (F ig . 5.4(c), 5.4(d)) and once after the s t ra in es t imat ion process (F ig . 5.4(e), 5.4(f)). T h e gradient operator was used i n b o t h cases as a s t ra in est imator . A 2 D median filter was preferred over I D because i t has a larger effect and i t makes i t possible to remove corrupted lines while a I D filter selects the median value from the same corrupted l ine. A s we expected, an improvement i n b o t h cases is not l inear since the median is a nonlinear filter. T h e results also shows that the median filter performs better i f i t is appl ied on the s t ra in estimates instead of displacement estimates wh ich could be explained by the fact that i n a single neighborhood s t ra in estimates are much closer than displacement estimates since the displacement accumulates whi le the strains stay the same. Therefore, we w i l l apply it after the s t ra in es t imat ion process for the rest of simulations. 5.2.3 Contrast-to-Noise Ratio Study T h e contrast-to-noise rat io ( C N R e ) i n elastography is an impor tant quant i ty that is related 5.2 ID Simulation of a 2D Image 66 Strain % (a) Improvement in SF with LSE. (b) Improvement in SNRe at 1%. 2D Median Kernel Size 3 5 7 9 2D Median Kernel S ize (c) Median filter is applied to displacement estimates. (d) Improvement in SNRe at 1%. 2D Median Kernel Size 3 5 7 9 2D Median Kernel Size (e) Median filter is applied to strain estimates. (f) Improvement in SNRe at 1%. Figure 5.4: Improvement in strain filters for TDPE when (a,b) LSE, (c,d) gradient operator follow-ing a 2D median filter and (e,f) 2D median filter following a gradient operator is used as a strain estimator. Compression is done from 0% to 15% in 0.2% steps. 100 realizations are used to show the mean values. The corresponding improvements in SNRe at 1% compression are also shown in (b, d and f) (D=40 mm, W = l mm, overlap=75%, f0=5 MHz and / s=40 MHz). 5.2 ID Simulation of a 2D Image 67 to the detectabi l i ty of a lesion [79]. T h e properties of the ul t rasound imaging system and signal processing algori thms described by the S F can be combined w i t h the elastic contrast properties of tissues w i t h s imple geometry [48], enabl ing predic t ion of the elastographic contrast-to-noise ra t io ( C N R e ) parameter. T h i s combined theoret ical mode l enables predic t ion of the elastographic C N R e for s imple geometry such as layered ( I D models) or circular lesions (2D models) embedded i n a uni formly elastic background. T h i s quant i ty includes the contrast characteristics of the output elastogram and is given by [11,126]: CNRe = 2 X [ S t ' f , (5.2) °t - n where St and Sf, represent the mean values of the est imated s t ra in i n the target and the back-ground, respectively and o~t and Ob are the s tandard deviations i n the est imated strains. F igure 5.5: C N R e of T D P E plot ted as a funct ion of contrast ra t io . C N R e is calcula ted from the S F and the s t rain and devia t ion values up to 6% s t ra in are used. Us ing 0.2% s t ra in i n the inclus ion and 0.2% to 6% s t ra in i n the background ( D = 4 0 m m , W = l m m , overlap=75%, / c = 5 M H z and / s = 4 0 M H z ) . T h e same da ta from the homogeneous phan tom at different compression rat ios were used to calculate the C N R e of T D P E at different contrast ratios between the target and background (F ig . 5.2 ID Simulation of a 2D Image 68 5.5). A strain of 0.2% was used as the target strain (inclusion) and 0.2% — 6% strain was selected as the background strain. Similarly to SNRe, the Matlab boxplot function was used to show the mean value and the deviation from the mean value of CNRe for 100 simulations. 5.2.3.1 Effect of L S E on Contrast-to-Noise Ratio The effect of the LSE on CNRe has not been presented in the literature yet. The improvement in CNRe is shown in Figure 5.6(a) when LSE with different kernel sizes (2-12) is used instead of the gradient operator on the same data that was used in Figure 5.5. As we expected it shows more than quadratic improvement in CNRe which confirms the results from [111], that estimates that CNRe w SNRe2, and [49] concludes that for a kernel size of N, the SNRe improves with TV3/2. Therefore the CNRe should improve with TV3. To investigate the improvement of CNRe more precisely when LSE is used, CNRe at a single compression ratio (1%) is shown in Figure 5.6(b). It shows a large improvement in CNRe as the kernel size of LSE goes up which means the detectability of the lesion increases rapidly as the kernel (a) Improvement in CNRe with LSE. (b) Improvement in CNRe at contrast ratio o f 5. Figure 5.6: CNRe of TDPE using ID simulated data and LSE(2-12) as an strain estimator. CNRe is calculated from SF and the strain and deviation values up to 6% strain are used, by using 0.2% strain in the inclusion and 0.2% to 6% strain in the background (a). 100 realizations are used to show the mean values. The improvement in CNRe at contrast ratio of 5 is also shown in (b) (D=40 mm, W = l mm, overlap=75%, f0=5 MHz and / s=40 MHz). 5.2 ID Simulation of a 2D Image 69 size of the L S E increases. A s a result i t seems that i t is better to use L S E w i t h large kernel size du r ing the detection process. B u t as we w i l l discuss later, increasing the L S E kernel size degrades the resolution. Therefore the operator needs to change the kernel size i n real t ime to opt imize the elastogram from b o t h resolut ion and C N R e points of v iew. 5.2.4 A x i a l Resolut ion S tudy (Ra) A number of definitions for the es t imat ion of the ax ia l resolution exist i n the l i terature [5,15,96]. It has recently been demonstrated i n [96] that the upper bound on the elastographic ax ia l resolut ion is p ropor t iona l to the length of the point spread function ( P S F ) of the transducer. T h e constant of propor t iona l i ty depends on the exact defini t ion that is used for the ax ia l resolut ion, but is generally a number between 1 and 2 [79]. T h i s l i m i t is not always achieved due to noise considerations [5]. T h e prac t ica l ly at tainable resolut ion for a given system, however, is generally l i m i t e d by the choice of s ignal processing parameters such as window size and overlap, as long as the window size is larger than the P S F of the system and other signal processing parameters such as s t ra in es t imat ion a lgor i thm. S imi la r to [15], i n this thesis ax ia l resolut ion is defined as the length of the t ime segment over wh ich the s t ra in changes from the value i n the background to the value i n the target (90%-10%). T h i s defini t ion is close to the definit ion of the filter accuracy i n signal processing. S imi l a r ly to signal processing where filters w i t h sharp t rans i t ion from pass-band to s top-band are preferred, elastograms w i t h h igh resolutions are also desired. T o derive the ax ia l resolut ion of our s t ra in est imator a 3 layered phan tom was s imulated w i t h a four t imes stiffer homogeneous layer inside ( F i g . 5.7). S imi la r ly to previous sections, the R F da ta were s imula ted at different compression ratios and time-delays were est imated by using the T D P E a lgor i thm. S t ra in was then est imated as the gradient of the t ime-delay estimates. T h e resolut ion was est imated at each l ine and the results were averaged over 100 independent lines. For W = 1 m m , over lap= 75% (AW= 25%) the ax ia l resolut ion est imated to be Ra = 0.3 m m . T h i s result is consistent w i t h results i n the l i terature, since it has been predicted that for a gradient 5.2 ID Simulation of a 2D Image 70 Figure 5.7: To simulate a 3 layered phantom, the stiffness of the inside layer was set to be 4 t imes stiffer t han its background. T h e stiffness map is shown i n (a) and its corresponding est imated elastogram image at 1% compression is shown i n (b). L S E w i t h a kernel size of 6 is used as an s t ra in estimator and a 3 x 3 median fi l tering was appl ied after the s t rain est imat ion s t ra in estimator the ax ia l resolut ion is almost equal to the window shift (Ra « AW) [5, 111, 127]. 5.2.4.1 Effect of L S E and M e d i a n Fi l ter ing on A x i a l Resolut ion T h e trade-off between achievable S N R e ( C N R e ) and ax ia l resolut ion i n elastography has been wel l s tudied i n [111] w i t h respect to the window size and window overlap. T h e trade-off was demonstrated for two impor tan t s ignal processing parameters (W, AW) but other signal processing parameters have not been s tudied yet. In this section, for the first t ime the trade-off w i t h two other parameters, namely kernel size of L S E and median filter is presented. T h e use of filtering typ ica l ly improves the S N R e at the expense of spatial resolution. Therefore, we expect the commonly used least-squares s t ra in est imator to have a poorer spat ial resolut ion than those est imated using the gradient operat ion. However the trade-off between S N R e and Ra for L S E has not been s tudied i n l i terature [111]. In order to measure the effect of L S E on the ax ia l resolut ion of the elastogram, the same simulated data from the three layered phantom used, but instead of the gradient estimator, L S E w i t h different kernel sizes va ry ing form 2 to 30 were used to estimate the strain. T h e resolut ion 5.2 ID Simulation of a 2D Image 71 4 8 12 16 20 24 28 1 3 5 7 9 11 LSE Kernel Size 2D Median Filter Kernel Size (a) Resolution of LSE at different kernel sizes. (b) Axial Resolution when elastogram is median fil-tered. Figure 5.8: Effect of f i l tering on the ax ia l resolution. For L S E the kernel size was var ied from 2 to 30 (a) and for 2 D M e d i a n fi l tering the kernel size was changed from 3 x 3 to 11 x 11 (b). 100 realizations are used to show the mean values. (D=40 m m , W = l m m , overlap=75%, / 0 = 5 M H z and / s = 4 0 M H z ) . was then est imated as the length of the t ime segment over w h i c h the s t ra in changes from the value i n the background to the value i n the target (90%-10%) ( F i g . 5.8(a)). A s predicted by [111] the result shows a l inear degradat ion i n resolution. Therefore, for a known window-shift ( A W ) , and L S E kernel size (N), the ax ia l resolut ion can be est imated according to E q u a t i o n (5.3): Ra « AW x N, (5.3) T h e gradient s t ra in es t imat ion technique results i n a l inear t rans i t ion of s t ra in from the value i n the background to the value i n the stiff layer inside, however, i t amplifies the noise. O n the other hand L S E results i n a nonlinear t rans i t ion from soft region to hard inclus ion but i t reduces the noise significantly [111]. T h e performance of the L S E can be concluded according to the results i n Sections 5.2.3 and 5.2.2. For a kernel size of N, on can state that: • S N R e improves w i t h TV3/2 • C N R e improves w i t h A ^ 3 5.2 ID Simulation of a 2D Image 72 • Ra degrades l inear ly w i t h N To s tudy the effect of median filtering, the 3 layered elastogram generated by the gradient operator was median filtered w i t h a number of different kernel sizes varying from 3 x 3 to 1 1 x 1 1 . The median filter was appl ied after the s t ra in es t imat ion process. A g a i n the resolut ion was est imated as the length of the t ime segment over wh ich the s t ra in changes from the value i n the background to the value i n the target (90%-10%). T h e result is shown i n F igure 5.8(b). A s we expected due to the nonl inear i ty of the median filter the resolut ion degrades nonlinearly. It is also clear from Figure 5.8(b) that whi le the median filter improves the elastogram (SNRe) i t doesn't degrade the resolut ion very much. Therefore it is because of these property that median filter is being widely used i n elastography. 5.2.5 Qual i ta t ive Evalua t ion For a qual i ta t ive evaluat ion of the T D P E a lgor i thm a 2D phantom (40 m m x 40 m m ) w i t h a 4 t imes stiffer circular inclus ion was s imulated using a I D s imulat ion to generate the R F da ta (F ig . 5.9). These types of phantoms are widely used since they simulate the tumors inside the body. In order to see the dynamic range of the a lgor i thm, the phantom was compressed from 0.5% to 10% w i t h a 0.5% strain step (20 frames). T D P E w i t h W = 1 m m and 75% overlap was used to estimate the mot ion and the s tandard gradient operator was used to estimate the s t ra in at different 1.2-5 0.8 06 0.4 20 40 60 SO 100 Depth (a) Stiffness map. (b) Ideal elastogram. (c) Ideal strain at center. F igure 5.9: To simulate a phantom w i t h a circle inc lus ion the stiffness of a circle inside a phantom was set to be 4 t imes stiffer than its background. T h e stiffness map is shown i n (a) and its corresponding ideal s t ra in image at 1% compression is shown i n (b,c). 5.2 ID Simulation of a 2D Image 73 compression ratios. A s shown i n F igu re 5.10, the a lgor i thm is capable of es t imat ing strains at h igh compression ratios. T h e colorbar besides each elastogram shows the m a x i m u m available s t ra in wh ich is p ropor t iona l to the compression rat io . A s the compression goes up the elastogram starts to degrade due to the decorrelat ion but the inc lus ion is s t i l l detectable even at 10% compression ra t io (F ig . 5.10(f)). T D P E fails at very h igh compression ratios because the displacements of the neighbors are not accurate anymore. Since the T D P E me thod depends on this informat ion i t fails to predict the search area correct ly and cannot estimate the mot ion . A s shown i n F igu re 5.10 corrupted regions start to appear i n the elastograms at very h igh strains. T o show the effect of the Least Squares S t ra in Es t ima to r ( L S E ) , the time-shift estimates at 1% compression were selected and elastograms were obtained w i t h L S E of different kernel sizes ( F i g . 5.11). It can be seen that as the kernel size (TV) of the L S E increases, the elastogram improves rap id ly and looks more like the ideal elastogram (F ig . 5.9) which is due to the improvement i n S N R e (AT3/2). F i g . 5.11 also shows the trade-off between the S N R e and ax ia l resolut ion (Ra). A s the kernels size of L S E increase the boundary of the inc lus ion also fades. To clarify th is , the s t ra in at the center l ine corresponding to each elastogram is also shown below the elastograms. C o m p a r i n g i t w i t h an ideal s t ra in (F ig . 5.9(c)), the ripples i n the h igh s t rain and low s t ra in regions are inversely propor t iona l to the S N R e and the slope of the t ransi t ions between h igh s t ra in i n the background ^ and low s t ra in i n the inclus ion are also inversely propor t iona l to the resolut ion. Therefore, as the kernel size of the L S E increases, the r ipple i n b o t h regions decreases which results i n an improvement i n S N R e . However, as the kernel size increases, the slope of the t rans i t ion region between h igh s t rain and low s t ra in region decreases which means the resolut ion degrades. S imi la r to L S E , to show the effect of med ian fi l tering on the est imated elastograms, the t ime-shift estimates at 1% compression were selected and the gradient operator was used to create the s t ra in image. To generate the final elastogram, the s t rain image was median filtered w i t h different 2D kernel sizes (F ig . 5.12). It can be seen that as the kernel size (N x Af) of the 2 D median filter increases, the elastogram improves rapidly, wh ich is due to the improvement i n S N R e . F i g . 5.12 also shows the trade-off between the S N R e and Ra. S imi la r to L S E , as the kernel 5.2 ID Simulation of a 2D Image 74 (d) 7% compression. (e) 9.5% compression. (f) 10% compression. Figure 5.10: To demonstrate the dynamic range of the T D P E , the elastograms at different com-pression ratios are shown i n th is figure. T h e inclusion is s t i l l detectable at 10% compression (D=40 m m , W = l m m , overlap=75%, / 0 = 5 M H z and / 3 = 4 0 M H z ) . size of the median filter increases, the boundary of the inc lus ion fades. T h e s t rain at the center line of each elastogram shows this fact better. A s the kernel size of the median filter increases, the r ipple i n bo th regions decreases wh ich shows the improvement i n S N R e . However the slope of the t rans i t ion region decreases wh ich means the resolution degrades. B u t this t ime the resolution does not degrade w i t h the same rate as L S E . For further qual i ta t ive evaluat ion of the T D P E , another 2D phantom was s imulated which consists of three layers w i t h a harder layer inside (F ig . 5.7). These types of phantoms are also wide ly used since they s imulate different layers of tissue inside the body. T y p i c a l l y u l t rasound images from the body start w i t h a layer of fat, which is soft, and another layer of muscle, wh ich is hard , following another soft layer. I D s imula t ion was used to generate the R F da ta at 1% compression and T D P E was appl ied to estimate the time-shifts. T h e effect of the L S E on the elastogram and the trade-off between S N R e and Ra is shown i n F igure 5.13 and the performance of the median filter on the s t ra in image, est imated by the gradient operator, is shown i n F igure 5.14. 5.2 ID Simulation of a 2D Image 75 (a) Elastogram with gradient. (b) LSE with kernel size of 4. (c) LSE with kernel size of 8. (d) Center line of elastogram. (e) Center line of LSE (4). (f) Center line of LSE(8). so too (g) LSE with kernel size of 12. (h) LSE with kernel size of 16. (i) LSE with kernel size of 24. (j) Center line of LSE(12). (k) Center line of LSE(16). (1) Center line of LSE(24). Figure 5.11: The improvement in elastogram when the Least Squares Strain Estimator (LSE) with different kernel sizes are used instead of gradient estimator. As the kernel size increases the elas-tograms improve rapidly. The smoother elastograms also show the loss of resolution. The boundary of the inclusion fades as the kernel size becomes bigger (D=40 mm, W = l mm, overlap=75%, / 0=5 MHz and /5=40 MHz). 5.2 ID Simulation of a 2D Image 7G (a) Elastogram with gradient (b) Median with 3 x 3 kernel. (c) Median with 5 x 5 kernel. (d) Center line of elastogram. (e) Center line of 3 x 3 kernel. (f) Center Line of 5 x 5 kernel. 20 •10 i« m i J i - H S,oo £ 5 ioo 120 0.5 120] 140 MJ 20 40 60 80 100 R F U n M (g) Median wi th 7 x 7 kernel. (h) Median wi th 9 x 9 kernel. (i) Median wi th 1 1 x 1 1 kernel. i . . . (j) Center line of 7 x 7 kernel. (k) Center line of 9 x 9 kernel. (1) Center line of 11 x 11 kernel. Figure 5.12: To demonstrate the effect of median fi l tering i n the elastograms' (SNRe) the or ig ina l elastogram is compared w i t h median filtered elastograms w i t h different kernel sizes. A s the kernel size increases the elastograms improve rapidly. T h e filtered elastograms also show the loss of resolution since the boundary of the inclusion fades as the kernel size becomes bigger (D=40 m m , W = l m m , overlap=75%, / e = 5 M H z and / s = 4 0 M H z ) . 5.2 ID Simulation of a 2D Image 77 (a) Elastogram with gradient. (b) L S E with kernel size of 4. (c) LSE with kernel size of 8. (d) Center line of elastogram. (e) Center line of LSE 4. (f) Center line of LSE 8. j I ; H i (g) LSE with kernel size of 12. (h) L S E with kernel size of 16. (i) LSE with kernel size of 24. (j) Center line of L S E 12. (k) Center line of LSE 16. j u t (1) Center line of LSE 24. Figure 5.13: The improvement in elastogram when the Least Squares Strain Estimator (LSE) is used instead of gradient estimator. As the kernel size increases the elastograms' (SNRe) improves rapidly while the resolution degrades (the boundary of the layeres fade) (D=40 mm, W = l mm, overlap=75%, / 0=5 MHz and / s=40 MHz). 5.2 ID Simulation of a 2D Image 78 (a) Elastogram with gradient (b) Median with 3x3 kernel. (c) Median with 5x5 kernel. (1%). 1-4 (d) Center line of elastogram. (e) Center line of 3 x 3 kernel. (f) Center Line of 5 x 5 kernel. 20 40 60 SO 100 Median with 7x7 kernel. (h) Median with 9x9 kernel. (i) Median with 11 X 11 kernel. (j) Center line of 7 x 7 kernel. (k) Center line of 9 X 9 kernel. (1) Center line of 11x11 kernel. Figure 5.14: T h e or ig inal e lastogram is compared w i t h M e d i a n F i l t e r ed elastograms w i t h different kernel sizes. A s the kernel size increases the elastograms' (SNRe) improve rap id ly while the reso-lu t ion degrades (the boundary of the layers fade) (D=40 m m , W = l m m , overlap=75%, / „ = 5 M H z and fs=A0 M H z ) . 5.3 Conclusion 79 Simi la r to previous simulations, F i g . 5.13 and 5.14 show the trade-offs between the elastographic signal-to-noise ra t io ( S N R e ) and ax ia l resolut ion (Ra). A s the kernel size increases for b o t h L S E and M e d i a n F i l t e r the elastogram improves r ap id ly due to the improvement i n S N R e whi le the boundary between different layers start to fade due to the degradation i n ax ia l resolut ion. 5.3 Conclusion T h e T D P E a lgor i thm was val idated i n this chapter. W h i l e the computa t iona l cost of T D P E was shown to be at least 1/25 of T D E and 1/10 of C A M and P R S , i t maintains the h igh S N R e , C N R e and Ra. I D s imulat ions of R F da ta showed that T D P E is even sensitive to 0.01% compression and is capable of es t imat ing s t ra in up to 10% whi le the T D E fails. S tudy of the ax i a l resolut ion of the T D P E showed that i t is capable of detecting the inclusions as smal l as 0.3 m m ( W = l m m , 75% window overlap and fa = 5 M H z ) . H a v i n g these properties makes the T D P E a suitable me thod for real t ime m o t i o n t r ack ing and s t rain imaging. It was also shown that the S N R e of T D P E can be improved i f the Least Squares E s t i m a t o r is used instead of the gradient operator but a trade-off exists, between the improvement i n S N R e and degradation i n ax ia l resolut ion. W e showed that for a kernel size of N the S N R e of T D P E improves w i t h N3/2 whi le the resolut ion degrades l inearly w i t h N. T h e effect of median fi l tering was also presented i n this chapter. It was shown that i t causes bo th a nonlinear improvement i n S N R e and a nonlinear degradation i n ax ia l resolut ion. It was shown that the effect of median fil tering on the ax ia l resolut ion was negligible whi le i t p rovided a good improvement i n S N R e . Therefore i t seems to be a good filter for elastogram. Chapter 6 Experiments In this chapter, two types of tissue m i m i c k i n g mater ia ls (phantoms) are tested to val idate the T D P E a lgor i thm w i t h real R F signals from the U l t r a s o n i x Rp500 ul t rasound machine (Ul t rason ix M e d i c a l Corpora t ion , Burnaby , B C ) . T h e chapter starts w i t h a descript ion of the exper imenta l setup. T h e n , the construct ion of two types of tissue m i m i c k i n g materials is detai led. T h e experi-menta l results are presented next, and remarks on the results conclude the chapter. 6.1 Experimental Setup In order to val idate the a lgor i thm qual i ta t ive ly w i t h real R F da ta rather t han s imula ted R F data, a number of experiments were run . A n overview of the set-up is i l lus t ra ted i n F igu re 6.1. T h e exper imental set-up is composed of an u l t rasound machine, an accurate compression device, a platform for ana lyz ing tissue m i m i c k i n g materials and a P C for off line processing. A 5-12 M H z linear probe connected to a P C - b a s e d ul t rasound machine (Ul t r a son ix R P 5 0 0 , U l t r a son ix M e d i c a l Corpora t ion , Burnaby , B C ) was used for a l l scans and R F da t a cap tur ing . T h e probe was placed on one side of the rectangular tissue m i m i c k i n g mater ia l , and the compressor was placed on the opposite side. N o lateral confinement was used i n the whole experiment . T h e tissue m i m i c k i n g mater ia l was lifted from the bo t t om surface dur ing the experiments to prevent f r ic t ion 80 6.1 Experimental Setup 81 Figure 6.1: Overv iew of the experimental set-up between the tissue m i m i c k i n g mate r ia l and the test platform. T h e tissue m i m i c k i n g mater ia l was squeezed from bo th sides between the probe and the compressor, to reduce the la teral mot ion on the contact surfaces. T h e phan tom was placed i n a pre-compression mode and pre-compression R F da ta were captured by ul t rasound machine. T o capture the post-compression R F da ta at a certain compression rat io, by considering the size of the phantom, the amount of compression was converted to an ax ia l displacement. Fo r instance, i f the phan tom is 6 c m i n its i n i t i a l state, i n order to apply a 1% ax ia l compression the phantom needs to be compressed 0.6 m m i n the ax ia l direct ion. B y using the mot ion vector the phantom was compressed i n the ax ia l d i rec t ion un t i l the required compression was appl ied to the phantom. After this step the second set of R F da ta were captured as a post-compression data . These two sets of R F da ta were used later i n a mot ion t racking a lgor i thm to generate the stat ic elastogram on another PC i n an off line process by using M a t l a b code. 6.2 Phantom Materials 82 T h e same software that had been used for s imulated R F da ta was used for off line processing. T h e software reads the two sets of R F da ta (one from pre- and one from post-compression) a n d runs a T D P E mot ion t rack ing a lgor i thm (Section 2.2) to f ind the displacement estimates. L S E was used as a s t ra in es t imat ion method. T h e generated elastograms were then color coded a n d displayed on the screen. Because of the s l ippery surface of the agar phantom, and to reduce the exper imental error, each experiment was done separately, w i t h separate pre- and post-compression R F data, instead of cap tur ing one pre-compression da ta and several post-compression data . 6.2 Phantom Materials A phan tom is a tissue m i m i c k i n g mater ia l . In the thesis, bo th acoustic and mechanical properties of tissue are considered i n phan tom construct ion. T h i s includes the texture of the ul t rasound image, the speed of sound, a t tenuat ion and stiffness [18,78-82]. T h e phantoms used i n this work are made from agar. A g a r is wide ly used i n the field of u l t rasound elastography because i t is easy to produce, non-toxic and its acoustic and mechanical properties is close to the human tissue [41]. A g a r is a water-based phantom. In order to produce an agar phantom the water is heated up to 90° C . T h e required percentage of agar powder is then added to the water. T h i s step needs to be done very slowly to make sure a l l the powder is disolved i n the water. T h e stiffness is control led by adjusting the p ropor t ion of agar powder i n the aqueous solut ion [41]. T h e mix ture is then cooled down to 80° C . T h e n a realist ic speckle pat tern of the B-mode image is produced for agar phantoms by adding cellulose scat ter ing particles (2% cellulose i n these experiments) into the l i q u i d mix tu re solut ion. T h e n the final mix tu re is cooled i n a m o u l d u n t i l i t congeals (around 26 ° C [41]). To be consistent w i t h the simulat ions, two types of phantoms were made, one w i t h a harder inclus ion and one layered phantom. 6.3 Strain Filter Study 83 6.3 Strain Filter Study T o generate the S F of the T D P E a lgor i thm a homogeneous phan tom (6 x 6 x 6 c m 3 ) was made. T h e phantom was prepared using a 2% agar solut ion w i t h 2% cellulose as scatterers. T h e phan tom was imaged to a dep th of 40 m m w i t h a linear array of 100 elements w i t h a 1 0 - M H z centroid frequency, 80% fract ional bandwid th transducer and digi t ized at 40 M H z . T h e phan tom was compressed from 0.2% to 5% i n steps of 0.2% s t ra in (25 experiments) and i n each compression rat io i t was imaged once for pre-compression R F da ta and once for post-compression R F data . T D P E w i t h W = 1 m m and 75% window overlap (160 window i n 40 mm) was used to estimate the mot ion at each step and the gradient operator was used to estimate the s t rain. T h e S N R e was estimated at each compression ra t io according to E q u a t i o n (5.1). 1.8h 1 2 3 4 5 Strain % Figure 6.2: T h e s t ra in filter of the T D P E a lgor i thm w i t h real da ta from experiments. Compress ion is done from 0.2% to 5% i n steps of 0.2%. T h e result from 100 R F lines are used to show the mean values and the devia t ion from the mean values. T h e gradient operator is used as a s t ra in est imator wi thout any fi l tering. T h e M a t l a b boxplot function was used to show the mean value and the devia t ion from the mean value of the S N R e for 100 R F A- l ines . It can be seen from F igure 6.2 that T D P E is capable of es t imat ing strains of up to 4% compression. 6.4 Qualitative Results 81 A s we expected the result from s imula t ion data over-estimates the performance of the T D P E a lgor i thm. T h i s could be explained by the fact that the tissue mot ion that occurs dur ing compres-sion is three dimensional (3D) whi le i t was assumed to be I D i n s imulat ion data. Fur thermore agar phantoms have s l ippery boundaries and they t rend to s l ip mos t ly at high compression ratios. T h i s change i n the probe's field of view w i l l cause the mot ion t rack ing to fai l . 6.4 Qualitative Results A n experiment was r u n for a phantom w i t h a harder inc lus ion. The experiment was performed on a 6 x 6 x 6 c m 3 uni formly elastic phantom w i t h a cy l ind r i ca l inclusion. The center of the cy l ind r i ca l inclus ion (radius of 7 m m ) was placed 15 m m below the surface. T h e phantom was prepared using a 2% agar solut ion w i t h 2% cellulose as scatterers for the background and 4% agar and 2% cellulose for the hard inclusion. F igure 6.3: Image on right shows the phantom w i t h cy l ind r i ca l inclusion and the image on left shows its corresponding B-mode image. T h e phan tom was imaged to a depth of 50 m m w i t h a linear array of 100 elements w i t h a 1 0 - M H z centroid frequency, 80% fractional bandwid th transducer and digi t ized at 40 M H z . A n image of the phan tom w i t h a regular u l t rasound B-mode image is shown i n F igure 6.3. Since the 6.4 Qualitative Results 85 RF Lines RF lines (a) Elastogram at 0.1% compression. (b) Median filtered elastogram at 0.1%. (c) Elastogram at 1% compression. (d) Median filtered elastogram at 1%. (e) Elastogram at 2% compression. (f) Median filtered elastogram at 2%. Figure 6.4: Elastograms of a phantom with 14 mm inclusion at 0.1%, 1% and 2%. LSE with kernel size of 6 is used to estimate the strain. Elastograms without applying a median filter (left) and median filtered with kernel size of 3 x 3 (right) are shown (D=50 mm, W = l mm, overlap=75%, / o=10 MHz and / s=40 MHz). 6.4 Qualitative Results 86 amount of cellulose (scatterers) i n b o t h the background and the inclus ion is the same, the inclus ion is not visible i n the regular B-mode . T h e phan tom was imaged once i n rest pos i t ion (pre-compression data) . It was then compressed uni formly and it was imaged again to capture the post-compression R F data . T h e experiment was repeated once for 0.1%, once for 1% and once for 2% compression. Af te r captur ing R F da ta t ime shifts were est imated by app ly ing the T D P E a lgor i thm. S t ra in was then est imated by the L S E method w i t h a kernel size of 6. T h e elastograms at strains of 0.1%, 1% and 2% are shown i n F igure 6.4. D u e to the stiffness contrast between the background and inclus ion, the inclusion is quite vis ible i n a l l elastograms. F igure 6.5: T h e image on the right shows the phantom w i t h the cy l indr ica l inclus ion and the image on the left shows its corresponding B-mode image. A n o t h e r experiment was r u n for a three layered phantom. T h e experiment was performed o n a 6 x 6 x 6 c m 3 uniformly elastic phan tom w i t h a harder layer inside. T h e center of the harder layer was placed 15 m m below the surface. T h e phantom was prepared using a 2% agar solut ion w i t h 2% cellulose as scatterers for upper and lower layers and 4% agar and 2% cellulose for the hard layer inside. S imi l a r ly to previous experiments, the phantom was imaged to a depth of 50 m m w i t h a linear array of 100 elements w i t h a 1 0 - M H z centroid frequency, 80% fractional bandwid th transducer and digi t ized at 40 M H z . A n image of the phantom w i t h a conventional u l t rasound B-mode image is 6.4 Qualitative Results 87 60 I RF Lines 100 120 (a) Elastogram at 0.5% compression. (b) Med ian filtered elastogram at 0.5%. (c) Elastogram at 1% compression. (d) Median filtered elastogram at 1%. (e) Elastogram at 2% compression. (f) Median filtered elastogram at 2%. Figure 6.6: Elastograms of a phantom with 14mm inclusion at 0.5%, 1% and 2%. LSE with kernel size of 6 is used to estimate the strain. Elastograms without applying a median filter (Left) and median filtered with kernel size of 3 x 3 (right) are shown (D=50mm, W=lmm, overlap=75%, /o=10 MHz and / a=40 MHz). 6.5 Conclusion 88 shown i n F igu re 6.5. A s i n the previous experiment, since the amount of cellulose (scatterers) i n bo th background and inc lus ion is the same the inc lus ion is not visible i n the regular B - m o d e image. T h e phan tom was imaged once i n rest pos i t ion (pre-compression data) and once after the compression (post-compression data). T h e phan tom was compressed uni formly once 0.5%, once 1% and once 2%. Af te r cap tur ing R F da ta t ime shifts were est imated by app ly ing the T D P E a lgor i thm. S t ra in was then est imated by the L S E method w i t h a kernel size of 6. T h e elastograms at strains of 0.5%, 1% and 2% are shown i n F igure 6.6. D u e to the stiffness contrast between the inside layer and background, the hard layer is detectable i n a l l elastograms. T h e elastograms also show the improvement when median filter is appl ied to the s t rain images. 6.5 Conclusion In this chapter, two types of tissue m i m i c k i n g materials (phantoms) were tested to val idate the T D P E a lgor i thm w i t h real R F signals from ul t rasound machine. T h e results show that the T D P E a lgor i thm is capable of es t imat ing the mot ion from the sequences of R F data, captured from an ul t rasound machine. Cons ider ing its computa t iona l cost, S N R e , C N R e , Ra and sensi t ivi ty T D P E seems to be an at t ract ive op t ion for real t ime appl ica t ion . Chapter 7 Implementation In this chapter, the current version of b o t h s t ra in imaging and stiffness imag ing packages are explained. T h e flow chart of bo th programs is shown and the details i n each box i n the flow chart are explained. T h e interface to the current version of bo th programs w i t h their features is demonstrated and the same tissue m i m i c k i n g phantoms that were used i n the previous chapter to generate the elastograms are used to val idate the performance of the software. T h e chapter finishes w i t h showing the current images tha t are generated i n real t ime from the tissue m i m i c k i n g materials when external deformations are appl ied to the phantoms. 7.1 Real Time Strain Imaging A real t ime s train imaging system has been implemented on the U l t r a son ix 5 0 0 R P ul t rasound machine (Ul t rason ix M e d i c a l Corpora t ion , Burnaby , B C ) , which is a P C - b a s e d machine w i t h a single 3.2 G H z C P U that provides R F da t a i n real t ime. 30% of the C P U is used by the m a i n program (Ul t rason ix service pack) itself and the remaining 70% of the C P U was used for R F da ta acquisi t ion, mo t ion t rack ing w i t h T D P E a lgor i thm, s t ra in est imation w i t h L S E , scal ing, color-coding and superimposed display on the B-mode images, or on the screen beside B - m o d e images. T h e flow chart of the software is shown i n F igure 7.1. T h e ma in server wh ich is the or ig ina l soft-89 7.1 Real Time Strain Imaging 90 ware on ul t rasound machine should be started before the s t ra in imaging software since i t provides b o t h B-mode and raw R F da ta for elastography programs. In i t i a l ly i n the software a connection establishes w i t h the m a i n server. Fo l lowing the in i t i a l i za t ion two independent threads are started (F ig . 7.1). One of the threads is used to generate B-mode and the other one is used for s t ra in imaging . T h e details about each thread are explained i n fol lowing sections. 7.1.1 B-mode Thread T h e B-mode thread captures the B-mode da ta from the m a i n server. S i m i l a r l y to a l l conven-t iona l u l t rasound machines, B-mode is then color coded i n gray scale i n R G B format. In order to display the image correctly i n b o t h B-mode and strain image, the software should know the informat ion about the current u l t rasound probe such as curvature, number of elements, p i t ch and so on . T h i s informat ion is then used to display the image on the screen correctly. Ideally the informat ion about the probe should be retrieved from the server, but i n the current version of the software this informat ion is set inside the software for different u l t rasound probes and the operator needs to select the type of probe manually. 7.1.2 Strain Imaging Thread St ra in is est imated from the displacement estimate. To estimate the displacement two sets of R F da ta are needed. A s shown i n F i g 7.1, the thread is an infinite loop. It always uses the current R F frame and the R F frame from the previous cycle to do the mot ion t racking . Before performing mot ion est imat ion, the thread checks whether two R F frames have the same number of lines and number of samples at each l ine. T h i s s imi la r i ty check makes i t possible to change the parameters of the m a i n server such as depth, line density or even changing the probe wi thout any need to close the s t ra in imaging software. Therefore a l l parameter settings can be done i n real t ime. After the da ta acquis i t ion process, R F frames are sent to the mot ion es t imat ion block. T h e fl 1-1 CD a CU STRAIN'IMAGINGiPRdCESS MOTION STRAIN ; | ^ ^ P L A C » f e < T ; ; V :::::: STRAINv::: •ESTIMATION POST ' S C A J . E j -:j C1W. ::;V: SUPER IMPOSING »™:-: : : aspu«MG ON SCREEN S COLOR CODING: B-MODE IMAGING PROCESS :coraNcs;:::;:: OtSSLAVINSONSCRSgN: fl 'Sb 03 fl '3 CS 03 . f l o o i—i fl bo 7.1 Real Time Strain Imaging 92 T D P E was used as a m o t i o n estimator due to i ts fast performance. Quadra t i c parabola f i t t ing (please see A p p e n d i x B.2) was used as a subpixe l l ing method to increase the accuracy of m o t i o n t racking and the three point Lagrang ian in terpola t ion me thod was used to recalculate the value of the correlat ion at the exact posi t ion of the correlat ion peak (please see A p p e n d i x B . 3 ) . I D median filter was then appl ied on the displacement estimates to improve the process of m o t i o n t racking [130]. W h i l e i t is not necessary to calculate the correla t ion image, i t was added to the software since it shows the accuracy of the displacement es t imat ion process. A l t h o u g h the cross correlat ion was set as a default correlat ion estimator, other correlat ion est imator methods (see A p p e n d i x B . l ) are implemented i n the software and i t can be set for the T D P E to use each of them. In add i t ion to the selection of the correlat ion estimator, the number of windows, the length of each window and the size of the search region can also be set for the a lgor i thm. Therefore the mot ion es t imat ion block generates bo th an est imated displacement image and its corresponding correlat ion image. 7.1.2.1 Strain Estimation L S E finds the best-fi t t ing l ine to a set of points inside each kernel. T h e slope of the line is then reported as the estimate value of the gradient of the time-shifts at the center of the kernel [49]. A s shown i f F igu re 7.2, these kernels have a large overlap w i t h each other. Therefore, computa t ion at each kernel may be useful for i ts ovelapping kernels. Therefore, i n order to investigate this idea, L S E was implemented. A c c o r d i n g to E q u a t i o n (4.3) we can wri te A and d wh i ch are N x 2 and N x 1 matrices as follows: i 1 i + 1 1 i + N-1 1 d(i) d(i + 1 ) d(i + N- 1) (7.1) then [A-^A] 1 reduces to a 2 x 2 ma t r i x and ATd to a 2 x 1 mat r ix . N o w i f we define 7.1 Real Time Strain Imaging 93 F igure 7.2: In the L S E method, neighboring kernels have large overlap w i t h eachother. Three overlapping kernels of size 5 is shown. ATA = u V , ATd = y W X z (7.2) the Least Squares E s t i m a t i o n (Equa t ion (4.4)) reduces to & u V -1 y b w X z ux — VW X —V (7.3) Since we are look ing for the value of a that is the slope of the crossing l ine, i ts value can be calculated as a = ^ ^ i . (7.4) ux — VW T h e value of a l l u,v,w,x,y and z can be calculated from E q u a t i o n (7.1) and E q u a t i o n (7.2) wh ich w i l l reduce to: 7.1 Real Time Strain Imaging 94 u = i 2 + (i + l ) 2 + . . . + (t + W - l ) 2 , (7.5) v = w = i + (i + l) + ... + (i + N-l), x = l + l + ... + l = N, y = d(i)2 + d(i + l ) 2 + . . . + d(i + N - l ) 2 , z = d(i) + d{i + 1) + . . . + d(i + N - 1), All of the above equations can be calculated in a single iteration loop which iterates up to the size of the kernel. It should be noted that the iteration comes from a matrix multiplication in which each row should be multiplied by its corresponding column. As can be seen from Equation(7.5), the values of u, v, w, x, y for estimating the slope of the current position can be used to calculate the values of u, v, w, x, y for the next position. Therefore after first calculation of u, v, w, x, y, the rest of them can be calculated as follows: = V>0ld = Wold %new = x o l d Vnew = Void znew — z o l d (i + N)2, (7.6) This implementation of L S E makes it possible to change the kernel size of the least squares method in real time to improve the quality of the elastogram's or its resolution. The computational cost of this implementation of the L S E method has been presented in Appendix (C.1). 7.1.2.2 Post Processing Post processing of the strain images is becoming a new field of research in elastography [72,86]. Therefore, this block is considered in the current implementation of the software and a number of 7.1 Real Time Strain Imaging 95 linear and nonlinear filters are implemented i n this block. Fo l lowing the s t ra in es t imat ion block a control swi tch is used to select between the available images from correlat ion, displacement and strain. T h e same control swi tch is then used to scale the da ta [0 255]. In this way the software is capable of showing s t ra in , displacement and correlat ion images i n real t ime. 1. M e d i a n F i l t e r ing : A 2 D median filter has been implemented i n the software and is being used as a post processing filter. T h e kernel size can be set i n the software and the operator has an opt ion to apply i t to the s t ra in image or not. T h e M e d i a n F i l t e r is normal ly used to reduce noise i n elastography [23,109]. It often does a better job than the M e a n F i l t e r of preserving useful detai l i n the elastogram. It considers each pixel i n the image i n t u r n and looks at its nearby neighbors to decide whether or not it is representative of i ts surroundings [36]. Instead of s imp ly replacing the pixel value w i t h the mean of neighboring p ixe l values, it replaces i t w i t h the median of those values. T h e median is calculated by first sor t ing a l l the p ixe l values from the surrounding neighborhood into numerical order and then replacing the pixel being considered w i t h the middle p ixe l value. T h e most t ime-consuming part of the median fi l tering is the sor t ing process such that for a median filter w i t h a kernel size of n x n one needs to sort n2 numbers and select i ts median for each p ixe l . For smal l kernel sizes the selection of the sort a lgor i thm has not much effect on the performance of the filter but as the kernel size of the increases the correct selection of the sort a lgor i thm becomes impor tan t . If fast sort ing algori thms such as Heap sort, Quick sort and Merge sort are used the order of complexi ty of median filter wou ld be n2Log2(n2) for each p ixe l [22] which makes i t too intensive and not applicable for real-t ime applicat ions. A fast 2 D median fi l tering a lgor i thm was introduced i n [45] that reduces the complexi ty of 2 D median fi l tering to 2n + 10 wh ich is much faster than the or ig inal version. T h e y used the previously sorted da ta of the neighbor p ixe l to speed up the sort ing process for the current p ixe l . In elastography software since the required kernel size was smal l (3 x 3) a simple implemen-ta t ion of the M e d i a n F i l t e r was used by app ly ing insertion sort [22]. 7.1 Real Time Strain Imaging 96 2. Tempora l F i l t e r i ng : S t ra in is not an intr insic proper ty of the tissue and the amount of s t ra in at each region depends on the amount of compression that is appl ied at that moment [88]. A s a result, i n h igh frame rate s t ra in imaging systems, the s t ra in images are not stable and they change r ap id ly on the screen. In order to have a more stable elastogram on the screen, several approaches have been implemented and embedded i n the post processing block. In order to make the s t ra in images more stable, s t ra in images can be averaged and the average s t rain can be shown instead of showing each image independently. To avoid the ex t ra buffer for this process, averaging was implemented according to the E q u a t i o n (7.7) and the user can set the number of frames to be averaged: o K X Saverage,K—1 ~f" Scurrent ^average,K — K -\-\ ' where i f is a counter that resets to one, whenever i t reaches its predefined l im i t , Saverage,K-i is the averaged s t ra in at step K — 1 and S c u r r e n t is the latest est imated s t rain. In this way the software starts to average two images then three and so on un t i l i t reaches the m a x i m u m averaging. T h e n i t starts again and averages the first two images. In the current implementa t ion , averaging does not register the elastograms which leads to losing the spat ia l resolut ion for large compressions, but since smal l compressions are mos t ly used i n elastography this is not a significant problem. Exponen t i a l averaging is another method to generate more stable s t ra in images [88]. It is simple to implement and does not need ex t ra buffering. It works according to E q u a t i o n (7.8). B y changing the value of the a f rom 0 to 1 different elastograms can be generated. Fo r values of a close to 0 the elastogram w i l l most ly show the most current s t ra in image wh ich is close to unfiltered elastograms and for the values of a close to 1 the elastogram w i l l mos t ly show the history of the s t ra in which is weighted average and is much more stable. T h i s post processing function was implemented i n the software and the amount of filtering can be selected by sett ing the value of a i n real-t ime. Sexponential,k = (1 Oi)Scurrent ~t~ OLS'exponential,k—\i (T•ty In addi t ion to these linear filters, nonlinear filters can also be appl ied to generate more stable 7.1 Real Time Strain Imaging 97 elastograms. For example, instead of just average, m a x i m u m or median value of the s t ra in among a number of images can be shown and the result can be called median strain imaging and maximum strain imaging. B u t since these kinds of imaging have not been s tudied yet they are not implemented i n the current version of the software. 7.1.2.3 Displaying Fol lowing the post processing procedure the final image is color coded. F i v e s tandard colormaps are implemented i n the software ( N O R M A L , H O T , H S V , J E T and G R A Y ) and the operator can easily swi tch between these colorcodes i n real-t ime. C o l o r coding provides much more dynamic range t han the gray scale color mapping . For example i n the H O T colormapping the input da ta should have a range of [0 768], i n H S V i t should have a range of [0 1280] whi le i n the gray scale the da ta needs to be i n the range of [0 255]. Ini t ia l ly , to superimpose the s t ra in image on top of the B-mode image, the B - m o d e image was drawn on the screen first. T h e strain image was then converted to transparent and the transparent s t ra in image was d rawn on top of the B-mode image at the same locat ion on the screen. T h i s implementa t ion made the software very slow and inappropr ia te for real t ime performance. To fix this problem the superimposing was s imp ly implemented by averaging B - m o d e and s t ra in images that were already converted to R G B format. To keep the flow chart s imple , the B - m o d e and s t ra in threads are shown to be quite independent. B u t i n order to do the averaging and super-imposing, the two threads need to communicate w i t h each other. In the current implementa t ion , the s t ra in thread provides the B-mode thread w i t h i ts s t ra in image i n R G B format. In the B - m o d e thread, after d rawing the or ig inal B-mode image on one side of the screen, to draw the super im-posed image on the other side of the screen, the B - m o d e thread shrinks the same da t a to make the size of bo th B-mode and s t ra in images equal. It then averages bo th images w i t h adjustable opacity. It then draws the final image on the other side of the screen beside the or ig inal B -mode . T o have an adjustable opaci ty the averaging was s imply replaced w i t h weighted averaging. Therefore to see more s t ra in image the operator needs to increase the weight of the s t rain image and to see more B-mode the operator needs to increase the weight of the B-mode i n the averaging process. 7.1 Real Time Strain Imaging 98 Af te r color coding the da ta i t is superimposed on the B-mode w i t h adjustable opacity. It is then d rawn on the screen correctly depending on the probe's type and specifications. F i n a l l y at the end of the cycle the previous R F frame is dropped and the current R F frame is recorded i n its place since it should be used as a previous frame for the next m o t i o n t rack ing cycle. 7.1.3 Software Features Cur ren t l y the system is capable of comput ing h igh frame rate s t ra in images wh ich are color coded and then superimposed on the B-mode image w i t h adjustable opacity. T h e current frame rate of the s t rain imaging is provided i n Table 7.1. T h e present interface of the software is also shown i n F igu re 7.3. Impor tant settings are gathered i n the control console that can be located on the left side of the screen. T h e control box is d iv ided into smaller boxes. T h e first box on top (Configuration Options) provides the operator to switch from simple averaging (Equa t ion (7.7)) to exponent ia l averaging (Equa t ion (7.8)). T h e first slider on top makes i t possible to change the value of K or a depending whether the simple or exponent ial averaging is desired. T h e second slider changes the kernel size of the L S E method which can be appl ied i n real t ime. T h e operator can also select to activate the super imposing or deactivate i t . If super imposing is active the t h i r d slider adjusts the opaci ty of the super imposing process. T h e second box on the control box (Color Map) provides the user to select the appropriate colormap for the s t ra in images. B y c l ick ing on each color map its corresponding colorbar appears i n a box below. Cur ren t ly the range of [0% to 2%] is colorcoded and higher strains are colorcoded s imi la r ly to 2% strain. T h i s range can be set i n the software. T h e t h i r d box on the control box (settings) provides a number of options for the user. T h e first check box selects the processing of the full image or just a smaller rectangle at the center of the image. T h i s helps the viewer to focus on the s t ra in d i s t r ibu t ion at the center of the image. T h e second check box switches between cross correlat ion and normal ized cross correla t ion (see A p p e n d i x B ) i n the T D P E method. T h e t h i r d check box activates or deactivates the med ian filter as a post processing filter. T h e fourth check box is considered for prostate imaging to swi tch the drawing 7.1 Real Time Strain Imaging 99 rawing |i,eastSt)Hdift£!iifiMtert|r sftrt*;?*^ r i m <5 Horn Hp $? V«rt«(np ttMbtfOllrtR 127 Mwrtw. Ol WneJo** 80 W«Sc*.L*nglh I X mm £ tetteraphy&rw* £l*Wogn£*» Group Figure 7.3: Current interface of the real time strain imaging system. The control box is on the left side of the screen and B-mode and strain images are drawn in the right side of the screen. from a rectangular view (which is required for a linear probe in sagittal view), to triangular view (which is necessary for a sector probe in transverse view). The last two check boxes flip both B-mode and strain images in horizontal and vertical direction. This feature helps the doctor to understand the anatomy by aligning physical probe orientation with the image displayed. There is another window (Information) which does not provide any setting and it only shows pertinent information in real time, such as frame rate, number of RF lines and number of windows which are being processed, the length of each window, window overlap in percentage and the sampling frequency of the ultrasound machine. Finally, at the bottom of the control box (Developing Group) the information about the devel-opers is included in front of their institute's logo. 7.2 Real Time Vibro Elastography 100 R F Lines x windows Pixels Processing time Algorithm frame rate Sys. frame rate 256 x 256 65,000 55 ms 18 fps 14 fps 192 x 192 37,000 34 ms 29 fps 20 fps 128 x 128 16,000 20 ms 50 fps 30 fps 100 x 100 10,000 13 ms 77 fps 40 fps 64 x 64 4,000 5 ms 200 fps 60 fps Table 7.1: T h e current frame rate of the real t ime s t ra in imaging system on Ul t r a son ix R P 5 0 0 w i t h a single P e n t i u m 4, C P U 3 . 2 G H Z , 1 G B R A M for a w indow size of 1mm (52 samples). 7.2 Real Time Vibro Elastography Simi l a r l y to real t ime s t ra in imaging, a real t ime stiffness imaging system based on v ib ro elas-tography [117] has been implemented. T h e flow chart of the software is shown i n F igure 7.4 and the process at each box w i t h its sett ing parameters are explained below. L i k e s t ra in imaging, the m a i n server should be started before the stiffness imaging software since i t provides bo th B - m o d e and raw R F data. Ini t ia l ly , a connection is established w i t h the m a i n server. Fo l lowing the in i t i a l i za t ion two independent threads are started. One of the threads is used to generate B-mode and the other one is used for stiffness imaging . T h e B-mode thread is the same thread as we used for s t ra in imaging (see Sect ion 7.1.1) except super imposing was not considered for stiffness images due to i ts low frame rate. T h e stiffness thread is discussed i n fol lowing section. 7.2.1 Stiffness Imaging Thread Stiffness imaging thread buffers the displacement estimates and runs a transfer analysis on the es t imat ion displacement. T h e stiffness is then est imated by using the value of the transfer funct ion at a specific range of frequency. Therefore, l ike s t ra in imaging , two sets of R F da ta are used at each buffering cycle to estimate the mot ion . Due to its real t ime nature, T D P E is used as a m o t i o n est imator and Quadra t i c P a r a b o l a f i t t ing B . 2 was used as a subpixel l ing method to increase the accuracy of mo t ion t racking. S imi l a r ly to s t ra in imaging, the cross correlat ion was set as a default correlat ion est imator while a l l of the correlat ion est imator methods that in t roduced i n B . l have o VIBRO ELASTOGRAPHY IMAGING PROCESS" CAPTURE RF I is;:-. OAWii't ;i P R E V I O U S R F l , Kx.-:;:rMTA¥-::?;r SAVE CURRENT :i;:;:;RFOATA::';:j' TRANSFER >::'?ANM.VBS"'':::: STIFFNESS : /•M*ion :;'. ..'Mqtjbo:::. juntos: 1 •:*i:i-$£:;.r EXEAN PREVIOUS RF.OATA CURRENT R F . O A T A * BUFFER. DISPLAVINQ OW SCttBBN* c p t o R c o o i N S ; MODE IMAGING PROCESS :; CAPTURE K aS'MOdei:?: C C C O R C C O N S ' DISPIAVWGON S C R E E N ; cS 6X5 o C O J33 "CD O SH <P a CP CD a3 o CD M cS CD at, 7.2 Real Time Vibro Elastography 102 been implemented in the software and can be used. The motion estimation images are recorded in a large buffer. The software fills the buffer until it is full. Therefore, in addition to the selection of the correlation estimator, number of windows, length of each window and size of the search region, the size of the buffer, length of F F T window and percentage of overlap for F F T windows, can also be set for the algorithm. Currently the software buffers 16 frames of displacement estimates and F F T windows are 8 points with 50% overlap (3 windows). Similarly to strain imaging, the software always checks two sets of R F data to make sure they are similar according to their number of lines and length of the lines. This similarity check makes it possible to change the parameters of the main server such as depth, line density or even changing the probe in real time without any need to close the software. 7.2.1.1 Transfer Function Estimation As can be seen from Equation (1.3) £,(£) is repeated three times in the calculation of the PSD and CSD therefore if we combine two functions with a single Transfer Function Estimation, the number of required F F T calculations reduces to two series instead of four series since we just need to calculate the F F T s once for Xi(t) and once for Xj(t). Therefore the transfer function estimate at each window can be rewritten as: The transfer function is then estimated by averaging this partial estimate of the transfer function at each window [132]. The performance of the transfer function estimate can be improved even further by noticing the fact that the transfer function is always estimated with respect to the reference input (Hrefi(w)). This means that during the whole process, the input of all transfer functions is the same therefore we just need to calculate its F F T s once and use it to estimate the transfer functions for other blocks. This implementation reduces the number of F F T calculations to one series for each block. FFT{xj(t)) x conj(FFT(Xj(t))) FFT(xi(t)) x conj(FFT{xi{t))y (7.9) 7.2 Real Time Vibro Elastography 103 7.2.1.2 Stiffness Estimation Afte r es t imat ing the value of the transfer funct ion at each block w i t h respect to the reference block, the relative stiffness value at each block is est imated according to [117]: k2 Hf - Hf 0 Hf Hf Hf-Hf ki(l-Hf) 0 (7.10) T h i s approach to estimate the relative stiffness has a lot of computa t iona l cost and is not applicable for real t ime purposes. Therefore i n order to have this implemented i n real t ime: k2 ' M l - f f ? ) / ( H f - H f ) k3 k^l-Hf)/(Hf-Hf) = k,(l-Hf)/(Hf-Hf) . h(l-Hf)/(H?) T h i s implementa t ion for the stiffness es t imat ion was used i n real t ime appl icat ions. For further improvement i n the result of stiffness estimates, s imi lar to s t rain es t imat ion (as expla ined i n Section 4.5), the difference between the values of the transfer function was est imated by using the L S E method. T h i s method results i n a much smoother stiffness estimate. T h e tradeoff for using least squares me thod to estimate the stiffness has not been studied yet. 7.2.1.3 Post Processing and Displaying Fol lowing the stiffness es t imat ion block is the post processing block wh ich is much simpler compared to the s t ra in imaging software. Because of the buffering process, the frame rate of the stiffness imaging is slow (1 fps for 128 lines of R F da ta and 2 fps for 64 lines of R F da t a w i t h 16 frames 7.2 Real Time Vibro Elastography 104 buffering). Therefore the averaging was not implemented and the post processing block for stiffness imaging only consists of a 2D median filtering. T h e kernel size of the filter can be set i n the software. T h e frame rate of the stiffness imaging improves i f a r i ng buffer is used instead of simple buffer wh ich means that each t ime, instead of removing a l l the da ta from the buffer, on ly a number of displacement estimates is removed and replaced w i t h new displacement estimates. T h i s imple-menta t ion may help to increase the frame rate and w i l l be considered i n the next version of the software. S imi la r to s t ra in imag ing software, after post processing the stiffness estimate, the da ta are scaled and colorcoded. T o increase the contrast another ampl i f ica t ion step was added to this process so that the user can scale the result i n real t ime by changing the pos i t ion of a slider i n the control box. F i v e s tandard colormap was implemented i n the software ( N O R M A L , H O T , H S V , J E T and G R A Y ) and the operator can easily switch between these colorcodes. T h e same disp laying process as i n Section 7.1.1 was used for the stiffness image. In the current implementa t ion , the transfer function es t imat ion and m o t i o n t racking r u n i n the same thread (sequential) and the transfer function runs i f the m o t i o n t racking buffer is fu l l . Therefore, it is impor tan t to note that transfer function analysis is computa t iona l ly intensive and takes some t ime. Therefore the current R F frame is not va l id anymore for the next mot ion detect ion cycle and the software should drop b o t h current and previous R F frames at the end of transfer function est imat ion cycle. T h e thread w i l l start to capture R F da ta i n its next cycle and w i l l s tart to f i l l the buffer from the next stiffness es t imat ion. 7.2 Real Time Vibro Elastography 105 7.2.2 Software Features Cur ren t l y the system is capable of comput ing real t ime stiffness images. T h e most recent interface of the software is shown i n F igure 7.5. T h e control box on the left is very s imilar to the s t ra in imaging software except for the C o n -figuration B o x . In the v ib ro elastography software, besides the kernel size of the least squares and amount of ampl i f ica t ion, the user can also choose the start frequency and end frequency of w h i c h to average the value of the transfer function wi th . T h e user can also change the length of the F F T window used i n the transfer function estimate. 7.3 Verification 106 T h e real value of the start frequency and end frequency depends on the buffering frequency of the software. For example i f the software is buffering w i t h 16 fps, according to N y q u i s t rate the m a x i m u m frequency component that we can capture is 8 H z . Therefore for a F F T window length of 8 (4 posi t ions are used for real parts and 4 posit ions are used for imaginary parts) . T h e pos i t ion of 0 w i l l be 0 H z , pos i t ion of 1 w i l l be 2 H z , pos i t ion of 2 w i l l be 4 H z , pos i t ion of 3 w i l l be 6 H z and pos i t ion of 4 w i l l be 8 H z . Increasing the F F T length may help increase the resolut ion of each frequency components but it also requires more da ta to be buffered wh ich makes the a lgor i thm even slower. 7.3 Verification T o show the current performance of bo th packages, they were tested on the same phantoms that were used i n Section 6.4. One sample of bo th real t ime s train image and stiffness image from the p h a n t o m w i t h inc lus ion and the three layered phan tom is shown i n F igure 7.6. T o generate the s t ra in image bo th phantoms were deformed w i t h an external v ib ra to r (single frequency s inusoid w i t h ampl i tude of 1.5 m m and frequency of 0.5 H z ) . A s i t can be seen from the control box, the kernel size of L S E was set to be 6, exponent ial averaging ( a = 0.79) was used to stabil ize the images and a 2 D median filter (wi th a kernel size of 3 x 3) was appl ied to the s t ra in images. To generate the stiffness image b o t h phantoms were deformed w i t h an external v ibra tor (am-plified whi te Gauss ian noise bandpassed from 4 H z to 6Hz) . T h e value of the transfer funct ion was averaged from 4 H z to 6Hz to include the exci ta t ion frequency component. A s i t can be seen from the control box, the kernel size of L S E was set to be 6 and a 2D median filter (wi th a kernel size of 3 x 3) was appl ied to the stiffness images. T h e results from F igure 7.6 confirm the correct implementa t ion of b o t h s t ra in imag ing and stiffness imag ing programs. It should be noted that the same colorcoding for b o t h s t ra in image and stiffness image result i n reverse colormap. T h i s means dark region i n s t ra in image appears as a bright region i n stiffness image. T h i s phenomenon is due to the fact that soft region experience more s t ra in and hard region experience less s t ra in. Therefore hard inc lus ion w i l l appear dark i n 7.3 Verification 107 (a) Sample of real time strain image of a phantom (b) Sample of real time stiffness image of a phantom with inclusion. with inclusion. (c) Sample of real time strain image of a three lay- (d) Sample of real time stiffness image of a three ered phantom. layered phantom. F igure 7.6: Screen shots of the interface and acquired images olbained w i t h two developed software. T w o types of tissue m i m i c k i n g phantoms are tested. T h e s t ra in image and stiffness image of a phantom w i t h inc lus ion (top), and a three layered phan tom (bottom) are shown. the s t ra in image and w i l l appear bright i n stiffness image. For further verif ication the real t ime system was taken to the operat ing room and i t was used on a real tissue i n v ivo . T h e idea was to check i f the prostate w i l l show up i n the s t rain and transfer function images or not. Genera l ly prostate tissue is stiffer than its surrounding tissue. The ac tuat ion system is shown in F igure 7.7. T h e dua l imaging probe w i t h two sets of crystals make i t possible to image i n bo th sagit tal d i rect ion (wi th linear array crystals) and i n transverse direct ion (wi th convex array crystals) . T h e probe was mounted i n housing and was v ibra ted along 7.3 Verification 108 its sagi t ta l plane by two motor -dr iven four-bar linkages, along a direct ion set by the ro ta t ional stepper, w i t h its depth adjusted by the t ransla t ional stepper. A S imul ink code was used as a controller to vibrate the probe w i t h adjustable ampl i tude and range of frequencies. F igure 7.7: Pho tograph of ac tua t ion system, w i t h the t ransrectal ul t rasound ( T R U S ) probe mounted i n i t (left). A sketch of the procedure is also shown (right). A t each step, the probe was posi t ioned w i t h the doctor into the patient 's rec tum and the prostate was deformed w i t h an external v ibra tor (amplified whi te Gauss ian noise bandpassed from 0.5Hz to 5Hz) w i t h m a x i m u m 2 m m peak to peak ampli tude. To calculate the transfer function, the length of F F T windows was set to be 8 and 50% overlap between windows were considered. T h e value of the transfer function was then averaged from 1Hz to 6Hz to include the exc i ta t ion frequency components. T h e difference between the value of the transfer function was est imated w i t h the L S E ( I O ) and this difference was then shown as a transfer function image. T h e s t rain images were est imated w i t h L S E ( 8 ) and exponent ia l ly filtered ( a = 0.9) to generate stable images. The transfer function ( T F ) and s t rain images from the convex crys ta l array i n the transverse direct ion are shown i n F igure 7.8 besides their corresponding B-mode images. A s shown i n F i g . 7.8, bo th elastograms (transfer function images and s t rain images) show the boundary of the prostate clearly mos t ly i n apex and mid-gland. 7.3 Verification 109 (a) Apex B-mode image. (b) Apex T F image. (c) Apex strain image. (d) Mid-gland B-mode image. (e) Mid-gland T F image. (f) Mid-gland strain image. (g) Base B-mode image. (h) Base T F image. (i) Base strain image. F igu re 7.8: Pros ta te transverse B-mode images (left), transfer function images (middle) and strain images (right) . T h e mot ion appl ied was 2 m m peak to peak, band-pass filtered (0.5-4.5Hz) whi te noise. T h e images show the prostate i n apex (top), mid-gland (center) and base (bot tom). Fo r further evaluat ion of the software, the images from sagit tal planes are also presented. T h e elastograms are shown i n F igure 7.9 w i t h their corresponding B-mode images. A s the probe rotates (counter clockwise), the prostate starts to move out of the probe's filed of v iew. T h i s fact is quite vis ible i n bo th transfer function and s t ra in images while it is hard to see the boundary of the prostate i n the B-mode images. A s shown i n F igure 7.9, bo th elastograms show much more contrast compare to the or iginal B-mode images. 7.3 Verification 110 (a) Sagittal plane B-mode. (d) Sagittal plane -10° B-mode. (g) Sagittal plane -20° B-mode. (j) Sagittal plane -30° B-mode. (b) Sagittal plane T F . (e) Sagittal plane -10° T F . (h) Sagittal plane -20° TF . (k) Sagittal plane -30° TF . (c) Sagittal plane strain. (f) Sagittal plane -10°strain. (i) Sagittal plane -20°strain. (1) Sagittal plane -30° strain. F igure 7.9: Pros ta te sagit tal B-mode images (left), transfer function images (middle) and s t ra in images (right) . T h e mot ion applied was 2 m m peak to peak, band-pass filtered (0.5-4.5Hz) whi te noise. 7.4 Conclusion 111 7.4 Conclusion In this chapter, the current version of bo th s t ra in imaging and stiffness imag ing programs were ex-pla ined i n deta i l according to their flow chart . T h e interface of the current version of b o t h programs w i t h their features was demonstrated. T h e result from tissue m i m i c k i n g phantoms demonstrates the correct implementa t ion of b o t h programs. T h e software can perform i n real t ime and i n c l in ica l applicat ions wi thout ex t ra hardware overhead. Chapter 8 Summary and Future Work 8.1 Summary T w o real t ime elastography packages (strain imaging and stiffness imaging) have been presented. B o t h programs include a new m o t i o n t rack ing a lgor i thm. T h e s t ra in imag ing a lgor i thm estimates the s t ra in di rect ly from the displacement estimates. T h e stiffness imag ing buffers the displacement estimates and analyze the tissue m o t i o n w i t h a Transfer Func t ion . T h e stiffness of tissue is then est imated from the magnitudes of the transfer function. T h e specific contr ibut ions of this thesis can be summarized as follows: • Fast Motion Tracking Algorithm: T h i s thesis in t roduced a new method for real t ime mot ion t racking i n u l t rasound images which is very fast. T h e a lgor i thm was named Time Domain Cross Correlation with Prior Estimates since i t has a s imi lar structure to that of s tandard Time Domain Cross Correlation. T h e s tandard method is very accurate for smal l s t ra in but i t is computa t iona l ly intensive and not applicable for real t ime performance. • Comparison of Computational Cost: Computational cost is in t roduced as a new compar-ison factor for different m o t i o n t rack ing algori thms and the numbers of mul t ipl ies and adds was used as a comparison factor. W e compared the computa t iona l cost of the in t roduced method w i t h the P R S [68] and C A M [102] methods that are the state-of-the-art methods i n real t ime mot ion t racking algori thms [86,105]. W e showed that the proposed T D P E method is the most efficient a lgor i thm among the current real t ime methods and is at least 25 times 112 8.1 Summary 113 faster than the s tandard T D E method . Therefore T D P E offers an at tract ive op t ion for real t ime mot ion t racking . • Stra in E s t i m a t i o n Algor i thms: In addi t ion to the gradient operator and L S E m e t h o d as s tandard s t ra in es t imat ion algori thms, numerical differentiation methods have been intro-duced as an extension to the gradient operator and the least squares s t ra in es t imat ion me thod was extended to higher order polynomials . • S tudy ing the Stra in F i l ter of T D P E : T h e performance of the T D P E was s tudied based on three impor tan t metrics, namely the s t ra in filter, contrast-to-noise ra t io and resolut ion w i t h s imula t ion data . T h e s tudy of the s t ra in filter showed that the proposed a lgor i thm is very sensitive and that it can track motions due to very smal l compressions ( « 0.01%) a n d has a large dynamic range so i t can track mot ions due to compressions as large as 10%. Fur thermore by s tudying the ax ia l resolution of the a lgor i thm we showed that the T D P E a lgor i thm is capable of detecting an inclus ion as smal l as 0.3 m m (for W = 1 m m , 75% window overlap and f0 = 5 M h z ) . In add i t ion to the s imula t ion data, results from tissue m i m i c k i n g phantoms showed tha t the proposed a lgor i thm is capable of accurately es t imat ing the mot ion w i t h smal l amount of computa t ion . These properties make the proposed method an at tract ive op t ion for the applicat ions that need real t ime mot ion t racking. • S tudy ing the M e d i a n Fi l ter and the Least Squares Est imator: For the first t ime the performance of the median filter has been presented i n elastography. Quant i ta t ive da ta showed that the median filter improves the s t ra in filter whi le i t has a smal l negligible effect on the ax ia l resolution. Qual i ta t ive da ta was then used to confirm these results. Moreover , the effect of the L S E on the contrast-to-noise and ax i a l resolut ion has been studied. • R e a l T i m e Elastography System: T h e proposed mot ion t rack ing a lgor i thm was used i n two real t ime elastography systems on the U l t r a son ix 500 R P machine that require fast and accurate displacement estimator. T h e first system is real t ime s t ra in imaging wh ich shows the d i s t r ibu t ion of the s t ra in due to an external deformation i n real t ime. T h e s t ra in is est imated from the displacement estimates. Therefore the most impor tan t part i n real t ime s t ra in imaging is its real t ime mot ion es t imat ion which needs to be as fast as possible. U s i n g T D P E 8.2 Future Work 114 as a m o t i o n est imator made i t possible to have a h igh frame rate s t ra in imag ing system. T h e second system is real t ime v ib ro elastography that extracts not on ly stat ic but also dy-namic properties of the tissue. A computer-control led vibra tor induces m o t i o n over a range of frequencies and simultaneously the resul t ing displacements are est imated and recorded at mul t ip le locat ions and t ime instants. T h e m o t i o n da t a is then used to extract the mechan-ica l properties of tissue. Therefore, s imi la r ly to s t ra in imaging the bot t le neck of the v ib ro elastography is i ts mot ion t racking part . T o have i t i n real t ime, the m o t i o n es t imat ion part needs to be as fast as possible. S t ra in imag ing software consists of da ta acquis i t ion, mot ion t racking, s t ra in es t imat ion, color coding, super imposing and displaying the images on the screen. In add i t ion to us ing T D P E as a m o t i o n est imator, the Least Squares S t r a in E s t i m a t o r was also implemented according to Sect ion 7.1.2.1 to speed up the s t ra in es t imat ion part . V i b r o elastography software consists of da ta acquis i t ion, mot ion t racking, buffering, transfer function analysis, stiffness es t imat ion and color coding, superimposing and d isp lay ing the images on the screen. S imi l a r ly to s t ra in imaging , T D P E was used as a m o t i o n est imator. Cur ren t ly b o t h elastography programs can perform i n real t ime and are ready to use i n c l in ica l applicat ions. 8.2 Future Work T h e idea used i n this paper to speed up the computa t ion of displacement using cross-correlation can also be used for phase domain algori thms to do the phase unwrapping. M o s t of the current phase domain a lgor i thms must perform addi t iona l calculat ions to avoid aliasing and to unwrap the phase. A s i n our method, the phase shift of neighbors can be used to unwrap the phase shift of the current window wi thou t addi t iona l calculat ions since the phase shifts of neighboring windows are expected to be close to each other. T y p i c a l l y i n elastography, the ax ia l component of the s t ra in tensor is est imated by t ak ing the gradient of the ax ia l (along the beam propagat ion axis) displacement occurr ing after a quasi-static 8.2 Future Work 115 tissue compression [18,78-81] which requires a I D m o t i o n t racking a lgor i thm as proposed i n this thesis. In general, however, the tissue mot ion that occurs dur ing compression is 3 D . Because the la teral (perpendicular to the beam propagat ion axis and i n the scan plane) and elevational (perpendicular to the beam propagat ion axis and to the scan plane) mot ions are not measured, two major drawbacks are encountered [56]. F i r s t , the axia l elastogram takes into account on ly a sma l l part of the mechanical tissue m o t i o n informat ion. Second, undesirable la tera l and elevational mot ions are the p r imary causes of s ignal decorrelat ion [51,52]. Therefore 2 D and 3 D mot ion t rack ing algori thms are preferred and need to be considered i n the next version of the elastography programs. T h e s imula t ion environment w i l l be upgraded to 2 D and 3D model for generating the R F da ta wh ich are more realist ic compared to I D R F data . 2 D and 3 D R F da ta w i l l be used to investigate the performance of the T D P E a lgor i thm. T h e same da ta w i l l be used to investigate 2 D and 3 D mot ion t rack ing algori thms. Cur r en t l y the frame rate of the elastography software is l imi t ed by the u l t rasound machine that can not provide R F data at a faster rate. T h i s p rob lem w i l l be solved soon by the Ul t r a son ix Company . F ina l l y , c l in ica l applicat ions for in vitro and in vivo experiments w i l l be invest igated s tar t ing w i t h prostate imaging for cancer detection, and R F abla t ion for cancer treatment. Bibliography [1] S R A g l y a m o v , A R Skovoroda, J M R u b i n , M O ' D o n n e l l , and Eme l i anov S Y . Young ' s modulus reconstruction i n D V T elast ici ty imaging. In In: Linzer M, ed. 27th International Symposium on Ultrasonic Imaging and Tissue Characterization, page 174, A r l i n g t o n , V A , 2002. U l t r a son Imaging. S . K . A l a m and J . O p h i r . O n the use of baseband and R F signal decorrelat ion as tissue s t ra in estimators. Ultrasound in Medicine and Biology, 23:1427—1434, 1997. S . K . A l a m and J . O p h i r . Reduc t ion of signal decorrelat ion from mechanical compression of tissues by tempora l stretching: applicat ions to elastography. Ultrasound in Medicine and Biology, 23:95-105, January 1997. S . K . A l a m , J . O p h i r , and E . E . Konofagou . A n adaptive s t ra in est imator for elastography. IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control, 45:461-72, M a r c h 1998. S . K . A l a m , J . Oph i r , and T . Varghese. Elas tographic ax ia l resolut ion cr i ter ia : an exper imental s tudy. IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control, 47:304-309, January 2000. J C B a m b e r and N L B u s h . Freehand elast ici ty imaging using speckle decorrelat ion rate. In IEEE Symposium on Acoustical Imaging, pages 285-292 vo l . 22, Florence, Italy, 1996. I E E E . J C B a m b e r and N L B u s h . Freehand elast ici ty imaging using speckle decorrelat ion rate. In Acoustical Imaging, pages 285-292 Vol .22 , P l e n u m Press, N e w Y o r k , 1996. J . B a n g , T . Selbekk, and G . Unsgard . S t ra in Processing of Intraoperat ive Ul t r a sound Images of B r a i n Tumours . In Proceedings of the Third International Conference on the Ultrasonic Measurement and Imaging of Tissue Elasticity, page 56, 2004. J .S . Bendat and A . G . P ie rso l , editors. Random Data. Analysis and Measurement Procedures, 2nd Edition. J o h n W i l e y and Sons, 1986. M Ber t r and , M Meunie r , M Doucet , and G . Fer land. Ul t rason ic b iomechanical s t ra in gauge based on speckle t racking. In IEEE Ultrasonics Symposium, page 859864 vol .2. I E E E , 1989. M . B i l g e n and Insana M F . P r e d i c t i n g target detectabil i ty i n acoustic elastography. In IEEE Ultrasonics Symposium, pages 1427-1430 vol .2. I E E E , 1997. 116 BIBLIOGRAPHY 117 L . N . Bohs and G . E . Trahey. A novel me thod for angle independent ul trasonic imaging of b lood flow and tissue mot ion . IEEE Transactions on Biomedical Engineering, 38:280-6, M a r c h 1991. 0 . Bonnefous. Measurement of the complete (3D) veloci ty vector of b lood flows. In IEEE Ultrasonics Symposium, page 795799. I E E E , 1988. W . L . B rogan , editor. Modern Control Theory, Third Edition, chapter 6. Prent ice H a l l , 1991. 1. Cespedes. Elastography. Imaging of biological tissue elasticity. P h D thesis, Un ive r s i ty of B r i t i s h C o l u m b i a , 1993. I. Cespedes, Y . H u a n g , J . O p h i r , and S. Sprat t . M e t h o d s for the es t imat ion of subsample time-delays of d ig i t ized echo signals. Ultrasonic Imaging, 17:142-171, 1995. I. Cespedes and J . Oph i r . Reduc t ion of image noise i n elastography. Ultrasonic Imaging, 15:89-102, 1993. I. Cespedes, J . O p h i r , H . Ponnekant i , and N . F . M a k l a d . Elas tography: elastici ty imag ing using ul t rasound w i t h appl ica t ion to muscle and breast i n v ivo . Ultrasonic Imaging, 15:73-88, 1993. P . Cha tu rved i , M . Insana, and T . H a l l . 2D companding for noise reduct ion i n s t ra in imaging . IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control, 45:179191, 1998. P . Cha tu rved i , M . Insana, T . H a l l , and M . B i l g e n . 3D companding using l inear arrays for improved s t ra in imaging. In IEEE Ultrasonics Symposium, page 14351438, Toronto, C a n a d a , October 1997. I E E E . A . D . Chris tensen, editor. Ultrasonic Bioinstrumentation, chapter 6. J o h n W i l e y and Sons, 1988. T . Cormen , C . Leiserson, R . Rivers t , and C . Ste in , editors. Introduction to Algorithms. Second Edition, chapter 3. M I T Press, 2001. E . Davies , editor. Machine Vision: Theory, Algorithms and Practicalities, chapter 3. A c a d -emic Press, 1990. P . de Jong, T . A r t s , A . Hoks , and R . Reneman. De te rmina t ion of Tissue M o t i o n Ve loc i ty by Cor re la t ion Interpolat ion of Pu l sed Ul t rasonic E c h o Signals. Ultrasonic Imaging, 12:84-98, 1990. C . D e K o r t e , E . Cespedes, A . V a n Der Steen, and C . Lancee. Intravascular elast ici ty imaging using ul t rasound: Feas ib i l i ty studies i n phantoms. Ultrasound in Medicine and Biology, 23:735-746, 1997. C . L . D e K o r t e , A . F . W . van der Steen, E . I . Cspedes, G . Pas terkamp, S . G . Car l ie r , F . M a s t i k , A . H . Schoneveld, P . W . Serruys, and N B o m . Charac ter i sa t ion of plaque components and vulnerabi l i ty w i t h Intravascular U l t r a sound Elas tography. Physics in Medicine and Biology, 45:1465-1475, 2000. BIBLIOGRAPHY 118 P deJong, T . A r t s , A . Hoeks, and R . Reneman . E x p e r i m e n t a l evaluat ion of the correla t ion in terpola t ion technique to measure regional tissue veloci ty . Ultrasonic Imaging, 13:145-161, A p r i l 1991. R J D i c k i n s o n and H i l l C R . Measurement of soft tissue mot ion using correla t ion between A-scans. Ultrasound in Medicine and Biology, 8:263-271, 1982. D D o t t i , R L o m b a r d i , and P . P i a z z i . Vec to r i a l measurement of b lood veloci ty by means of ul t rasound. Medical and Biological Engineering and Computing, 30:219225, M a r c h 1992. M . Doyley, J . Bamber , F . Fuechsel, and N . B u s h . A freehand elastographic imaging approach for c l in ica l breast imaging: System development and performance evaluat ion. Ultrasound in Medicine and Biology, 27:1347-1357, 2001. V . N . Dvornychenko. Bounds on (deterministic) correlat ion functions w i t h app l ica t ion to registrat ion. IEEE Transactions on Pattern Analysis and Machine Intelligence, P A M I - 5 : 2 0 6 -13, M a r c h 1983. S Y . Emel i anov , X . Chen , M . O ' D o n n e l l , B . K n i p p , D . Myer s , T . Wakefield, and J M . R u b i n . T r ip l ex ul t rasound: E l a s t i c i t y imaging to age deep venous thrombosis. Ultrasound in Medicine and Biology, 28:757-767, 2002. S Y Emel ianov , M A L u b i n s k i , W F Wei tze l , and et a l . E la s t i c i t y imaging for early detect ion of renal pathology. Ultrasound in Medicine and Biology, 21:871-883, 1995. H . Eskanda r i and T . Salcudean. Tissue S t ra in Imaging U s i n g A Wavelet -Based Peak Search A l g o r i t h m . In Proceedings of the Third International Conference on the Ultrasonic Measure-ment and Imaging of Tissue Elasticity, page 31, 2004. A . Fertner and A . S jo lund. Compar i son of various t ime delay est imat ion methods by computer s imula t ion . In IEEE Transactions on Acoustics, Speech, and Signal Processing, pages 1329-1330 vol.34. I E E E , 1986. D . Fo r sy th and J . Ponce, editors. Computer Vision: A Modern Approach., chapter 3. Prent ice H a l l , 2003. S . G . Foster, P . M . Embree , and W . D . O ' B r i e n . F l o w veloci ty profile v i a t ime-domain cor-relat ion: error analysis and computer s imula t ion . Ultrasonics, Ferroelectrics and Frequency Control, 37:164-175, M a y 1990. D . F u . Spectral Strain Estimation Techniques for Tissue Elasticity Imaging. P h D thesis, D r e x e l Univers i ty , 2005. B S G a r r a , E I Cespedes, J O p h i r , and et a l . E las tography of breast lesions: In i t i a l c l in ica l results. Radiology, 202:79-86, January 1997. J T H a l l , Z h u Y , C S Spalding, and L T Cook . In-vivo results of real-time freehand elast ici ty imaging. In IEEE Ultrasonics Symposium, page 1649:1653. I E E E , 2001. T . J . H a l l , M . B i l g e n , M . F . Insana, and T . A . K r o u s k o p . P h a n t o m materials for elastogra-phy. IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control, 44:1355-65, November 1997. BIBLIOGRAPHY 119 [42] I A . H e i n . Mu l t i - d i r ec t i ona l ul trasonic b lood flow measurement w i t h a t r ip le beam lens. In IEEE Ultrasonics Symposium, page 10651069. I E E E , 1993. [43] I . A . H e i n , V . Suorsa, J . Zachary, R . F i s h , J . C h e n , W . K . Jenkins, and W . D . O ' B r i e n . Accura t e and precise measurement of b lood flow using u l t rasound t ime domain correlat ion. In IEEE Ultrasonics Symposium, pages 881-886 vol .2 . I E E E , Oc t 1989. [44] H . H o n m a , M . O h t a , and T . N i s h i t a n i . A simplif ied method of mot ion vector es t imat ion for V L S I implementa t ion . In ACM/SIGDA ninth international symposium on Field program-mable gate arrays, pages 216-219. I E E E Internat ional Conference on Systems Engineer ing , September 1992. [45] T . Huang , G . Y a n g , and G . Tang . A fast two-dimensional median f i l tering a lgor i thm. IEEE Acoustics, Speech, and Signal Processing, 27:13-18, Feb. 1979. [46] M . Insana, P . Cha tu rved i , T . H a l l , and M . B i l g e n . 3D companding using 1.5D arrays for improved s t ra in imaging. In IEEE Ultrasonics Symposium, page 14271430, Toronto , Canada , Oc tober 1997. I E E E . [47] J . A . Jensen and I .R. Lacasa . E s t i m a t i o n of b lood veloci ty vectors using transverse u l t rasound beam focusing and cross-correlation. In IEEE Ultrasonics Symposium, pages 1493-1497 vol .2 . I E E E , Oc t 1999. [48] F . K a l l e l , M . Ber t r and , and J . Oph i r . Fundamenta l l imi ta t ions on the contrast-transfer efficiency i n elastography: an analyt ic s tudy. Ultrasound in Medicine and Biology, 22:463-470, 1996. [49] F . K a l l e l and J . O p h i r . A least-squares s t ra in est imator for elastography. Ultrasonic Imaging, 19:195-208, 1997. [50] F . K a l l e l , J . O p h i r , K . Magee, and T . K r o u s k o p . Elas tographic imaging of low-contrast elastic modulus dis t r ibut ions i n tissue. Ultrasound in Medicine and Biology, 24:409-425, 1998. [51] F . K a l l e l , T . Varghese, J . Oph i r , and M . B i l g e n . T h e nonstat ionary s t ra in filter i n elastog-raphy, P a r t II . La t e r a l and elevational decorrelat ion. Ultrasound in Medicine and Biology, 23:1357-69, 1997. [52] O p h i r J . K a l l e l F . Three-dimensional tissue m o t i o n and its effect on image noise i n elastogra-phy. IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control, 44:12861296, 1997. [53] S. K h a w a m , T . A r s l a n , and F . Wes ta l i . E m b e d d e d reconfigurable array target ing m o t i o n es t imat ion applicat ions. In International Symposium on Circuits and Systems, pages 11-760-11-763 vol .2. I E E E , M a y 2003. [54] E . Konofagou, J . D'hooge, and J . O p h i r . M y o c a r d i a l elastography - A feasibil i ty s tudy i n v ivo . Ultrasound in Medicine and Biology, 28:475-482, October 2002. [55] E . Konofagou, F . K a l l e l , and J . Oph i r . Three-dimensional mot ion es t imat ion i n elastography. 1998 IEEE Ultrasonics Symposium. Proceedings, pages 1745-8 vol .2, October 1998. BIBLIOGRAPHY 120 [56] E . Konofagou and J . O p h i r . A new elastographic method for es t imat ion and imaging of la tera l displacements, la teral strains, corrected ax ia l strains and Poisson 's ratios i n tissues. Ultrasound in Medicine and Biology, 24:1183-99, October 1998. [57] E . Konofagou, J . O p h i r , F . K a l l e l , and T . Varghese. Elas tographic d y n a m i c range expansion using variable appl ied strains. Ultrasonic Imaging, 19:145-166, 1997. [58] E . E . Konofagou and J . O p h i r . P rec i s ion E s t i m a t i o n and Imaging of N o r m a l and Shear C o m p o -nents of the 3D S t ra in Tensor i n elastography. Physics in Medicine and Biology, 45:1553-1563, 2000. [59] E . E . Konofagou, T . Varghese, and J . O p h i r . Spectra l Es t imators i n Elas tography . Ultrasonics, 38:412-416, 2000. [60] E . E . Konofagou, T . Varghese, J . O p h i r , and S . K . A l a m . Power spectral s t ra in estimators i n elastography. Ultrasound in Medicine and Biology, 25:1115-29, September 1999. [61] F . K r e m k a u , editor. Diagnostic Ultrasound: Principles and Instruments, ^th Edition, chap-ter 1. W . B . Saunders Company, 1993. [62] T A K r o u s k o p , S V i n s o n , B Goode , and D . Dougherty. A pulsed Dopp le r ul t rasonic system for m a k i n g noninvasive measurements of the mechanical properties of soft tissue. Journal of Rehabilitation Research and Development, 24:1-8, 1987. [63] M . Krueger , A . Pesavento, H . E r m e r t , K . M . Hi l tawsky , L . Heuser, H . Rosenthal , and A . Jensen. Ul t rasonic s t rain imag ing of the female breast using phase root seeking and three-dimensional opt ica l flow. In IEEE Ultrasonics Symposium, pages 1757-1760 vol .2 . I E E E , Oc tober 1998. [64] R . M . Lerner , S .R. H u a n g , and K . J . Parker . 'Sonoelast ici ty ' images derived from ul t rasound signals i n mechanical ly v ib ra ted tissues. Ultrasound in Medicine and Biology, 16(3):231-239, 1990. [65] R M Lerner and K J . Parker . Sono-elast ici ty i n ul trasonic tissue character izat ion and echo-graphic imaging. In Proceedings ofthe 7th European Comm Workshop, Ni jmegen, T h e Nether-lands: Office of Pub l ica t ions of the E u r o p e a n Communi ty , 1987. [66] W L i , F M a s t i k , I Cespedes, S Car l i e r , and A Steen. Echo decorrelat ion est imated from signal powers Ul t rasound . Ultrasound in Medicine and Biology, 254:405409, 1999. [67] L . L j u n g . System Identification, Theory For The User. Prent ice H a l l P T R , N J , 1999. [68] A . Lorenz , A . Pesavento, M . Pesavento, and H . Ermer t . Three-dimensional s t ra in imaging and related s t rain artifacts us ing an ultrasonic 3D abdomina l probe. In IEEE Ultrasonics Symposium, pages 1657-1660 vol .2 . I E E E , October 1999. [69] A . Lorenz , H J Sommerfeld, M G Schurmann, and et a l . A new system for the acquis i t ion of ul t rasonic mult icompression s t ra in images of the human prostate i n v ivo . IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control, 46:969-1147-1154, 1999. BIBLIOGRAPHY 121 [70] A M L u b i n s k i , S Y Emel i anov , K R Raghavan, A E Yagle , A R Skovoroda, and M . O D o n n e l l . La t e r a l displacement es t imat ion using tissue incompressibi l i ty . IEEE Transactions on Ultra-sonics, Ferroelectrics and Frequency Control, 43:247255, 1996. [71] G E M a i l l o u x , F Langlo i s , P Y Simard , and M . Be r t r and . Res tora t ion of the veloci ty field of the heart from two-dimensional echocardiograms. IEEE Transactions on Medical Imaging, 8(2):143153, June 1989. [72] T . M a t s u m u r a , S. Tamano , R . Shinomura , T . M i t a k e , M . Yamakawa , T . Shi ina , A . I toh, and E . Ueno . P r e l i m i n a r y E v a l u a t i o n of Breast Disease Diagnosis Based on R e a l - T i m e E l a s t i c i t y Imaging. In Proceedings of the Third International Conference on the Ultrasonic Measurement and Imaging of Tissue Elasticity, page 58, 2004. [73] W . N . M c D i c k e n , editor. Diagnostic Ultrasonics: Principles and use of Instruments, 3rd Edition, chapter 1. C h u r c h i l l Liv ings tone , 1991. [74] N . N i t t a and T . Shi ina . M y o c a r d i a l s t ra in imaging based on three dimensional mot ion track-ing . In IEEE Ultrasonics Symposium, pages 1258-1261 V o l . 2 . I E E E , October 2003. [75] M . O ' D o n n e l , A . R . Scovoroda, B . M . Shapo, and S . Y . Emel ianov . Internal Displacement and S t ra in Imaging U s i n g Ul t rason ic Speckle Track ing . IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, 41:314-25, M a y 1994. [76] M . O ' D o n n e l l , S . Y . Emel ianov , A . R . Skovoroda, M . A . L u b i n s k i , W . F . Wei tze l , and R . C Wigg ins . Quant i ta t ive elast ici ty imaging. In IEEE Ultrasonics Symposium Proceedings, pages 893-903 vol .2 . I E E E , 1993. [77] M . O ' D o n n e l l , A . R . Skovorada, and B . M . Shapo. Measurement of ar ter ia l wa l l mot ion using Fourier based speckle t rack ing algori thms. In IEEE Ultrasonics Symposium Proceedings, pages 1101-1104 vol .2 . I E E E , 1991. [78] J . Oph i r , S . K . A l a m , B . G a r r a , F . K a l l e l , E . Konofagou , T . K r o u s k o p , and T . Varghese. Elas tography: ul t rasonic es t imat ion and imaging of the elastic properties of tissues. Ultrasonic Imaging, 213:203-33, 1999. [79] J . Oph i r , S . K . A l a m , B . S . G a r r a , F . K a l l e l , E . E . Konofagou , T . K r o u s k o p , C . R . B . M e r r i t t , R . R ighe t t i , R . Souchon, S. Srinivasan, and T . Varghese. Elas tography: Imaging the E las t i c Proper t ies of Soft Tissues w i t h Ul t r a sound . Medical Ultrasound, 29:155-171, 2002. [80] J . O p h i r , I. Cespedes, B . G a r r a , H . Ponnekant i , Y . H u a n g , and N . M a k l a d . Elas tography: u l -trasonic imaging of tissue s t ra in and elastic modulus i n v ivo . European Journal of Ultrasound, 3:49, January 1996. [81] J . Oph i r , I. Cespedes, H . Ponnekant i , Y . Y a z d i , and X . L i . Elas tography: a quant i ta t ive method for imaging the elast ici ty of biological tissues. Ultrasonic Imaging, 13:111-134, A p r i l 1991. [82] J . Oph i r , B . G a r r a , F . K a l l e l , E . E . Konofagou, T . K r o u s k o p , R . R ighe t t i , and T . Vargh -ese. Elas tographic Imaging. Ultrasound In Medicine and Biology, The Future of Biomedical Ultrasound, 26:S23-S29, 2000. BIBLIOGRAPHY 122 [83] J . Pa rk , S. K w o n , and M K . Jeong. A S tudy of V i b r a t i o n Character is t ics for Ul t r a son ic E l a s t i c i t y Imaging. In Proceedings of the Third International Conference on the Ultrasonic Measurement and Imaging of Tissue Elasticity, page 106, 2004. [84] K . J . Parker , S .R. H u a n g , R . A . M u s u l i n , and R . M . Lerner . Tissue response to mechanical v ibra t ions for 'sonoelastici ty imaging ' . Ultrasound in Medicine and Biology, 16(3):241-246, 1990. [85] F . Pascal . A para l le l stereo a lgor i thm that produces dense depth maps and preserves image features. In Machine Vision and Applications, pages 35-49 vol .6, 1993. [86] A . Pesavento and A . Lorenz . R e a l T i m e S t ra in Imaging a new Ul t rasonic M e t h o d for Cancer Detec t ion: F i r s t S tudy Resul ts . In IEEE Ultrasonics Symposium, pages 1647-1652. I E E E , 2001. [87] A . Pesavento, A . Lorenz , H . E r m e r t , H . Sommerfeld, M . Garc ia -Schurmann, T h . Senge, and S. P h i l i p p o u . Frame-to-frame statistics of real- t ime s t ra in images. In IEEE Ultrasonics Symposium, pages 1809-1812 vol .2. I E E E , October 2000. [88] A . Pesavento, A . Lorenz , Siebers S., and H . E r m e r t . F rame T o Frame F i l t e r i n g O f R e a l T i m e S t ra in Images. In IEEE Ultrasonics Symposium. I E E E , 2002. [89] A . Pesavento, A . Lorenz , S. Siebers, and H . E r m e r t . N e w real-time s t ra in imaging concepts using diagnostic u l t rasound. Physics in Medicine and Biology, 45:1423-1435, 2000. [90] A . Pesavento, C . Perrey, M . Krueger , and H . E r m e r t . A t ime efficient and accurate s t ra in es t imat ion concept for ul trasonic elastography using i terative phase zero es t imat ion. IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, 46:1057-1067, 1999. [91] H Ponnekant i , J O p h i r , Y H u a n g , and I Cespedes. Fundamenta l mechanical l imi ta t ions on the v isua l iza t ion of elast ici ty contrast i n elastography. Ultrasound in Medicine and Biology, 21:533-543, 1995. [92] S. Ramachandran and S. Sr inivasan. F P G A implementa t ion of a novel, fast mo t ion es t imat ion a lgor i thm for real-t ime video compression. In ACM/SIGDA ninth international symposium on Field programmable gate arrays, pages 213-219. Internat ional S y m p o s i u m on F i e l d P r o -grammable G a t e A r r a y s , 2001. [93] B S R a m a m u r t h y and G E . Trahey. Po ten t ia l and l imi ta t ions of angle-independent flow de-tect ion a lgor i thms using radio-frequency and detected echo signals. Ultrasonic Imaging, 13:252268, 1991. [94] J . Reve l l , M . M i r m e h d i , and D . M c N a l l y . Technica l report on Ul t r a sound Speckle T rack ing For S t ra in E s t i m a t i o n , December 2003. [95] R R ighe t t i , F K a l l e l , R J Stafford, and et a l . E las tographic characterizat ion of H I F U - i n d u c e d lesions i n canine l ivers. Ultrasound in Medicine and Biology, 25:1099-1113, 1999. [96] R . R ighe t t i , J . O p h i r , and P . K tonas . A x i a l Reso lu t ion i n Elas tography. Ultrasound in Medicine and Biology, 28:101-113, January 2002. BIBLIOGRAPHY 123 [97] R . R i g h e t t i , J . O p h i r , S. Srinivasan, and T . Krouskop . T h e Feasibi l i ty of U s i n g Elas tography for Imaging the Poissons R a t i o i n Porous M e d i a . Ultrasound in Medicine and Biology, 30:215-228, 2004. [98] R . R i g h e t t i , J . O p h i r , and K r o u s k o p T . Poisson's R a t i o Imaging and Poroelas tography i n B i o l o g i c a l Tissues: A Feasibi l i ty S tudy . In Proceedings of the Third International Conference on the Ultrasonic Measurement and Imaging of Tissue Elasticity, page 94, 2004. [99] R . R i g h e t t i , S. Srinivasan, and J . O p h i r . L a t e r a l Reso lu t ion i n Elas tography. Ultrasound in Medicine and Biology, 29:695-704, January 2003. [100] J M R u b i n , S R A g l y a m o v , T . Wakefield, M . O ' D o n n e l l , and S Y Emel i anov . C l i n i c a l A p p l i c a -t i o n of Sonographic E las t i c i ty Imaging for A g i n g of Deep Venous Thrombos i s : P r e l i m i n a r y F ind ings . Ultrasound Medicine, 22:443-448, 2003. [101] A . R y s z k o and K . W i a t r . A n assessment of F P G A sui tabi l i ty for implementa t ion of real-t ime m o t i o n est imat ion. In Euromicro Symposium on Digital Systems and Design, pages 364-367. I E E E , September 2001. [102] T . Sh i ina , M . Doyley, and J . Bamer . S t ra in imaging using combined R F and envelope auto-correla t ion processing. In IEEE Ultrasonics Symposium, pages 1331-1336. I E E E , 1996. [103] T . Sh i ina , N . N i t t a , E . Ueno, and J . C . Bamber . R e a l T i m e Tissue E l a s t i c i t y Imaging using C o m b i n e d Autocor re la t ion M e t h o d . Medical Ultrasonics, 26:57-66, 1999. [104] T . Sh i ina , M . Yamakawa , N . N i t t a , E . Ueno , T . M a t s u m u r a , S. Tamano , and T . M i t a k e . C l i n i c a l assessment of real-time, freehand elast ici ty imaging system based on the combined autocorre la t ion method. In IEEE Symposium on Ultrasonics, pages 664-667 V o l . 1 . I E E E , Oc tober 2003. [105] T . Sh i ina , M . Yamakawa , N . N i t t a , E . Ueno , T . M a t s u m u r a , S. Tamano , R . Shinomura , and T . M i t a k e . Breast and Prostate Cancer Diagnosis by R e a l - T i m e Tissue E l a s t i c i t y Imaging Sys tem - A Feasibi l i ty Study. In 2nd International Conference on the Ultrasonic Measurement and Imaging of Tissue Elasticity, page 60, 2003. [106] S. Sr inivasan, F . K a l l e l , and J . Oph i r . E s t i m a t i n g the elastographic signal-to-noise ra t io using corre la t ion coefficients. Ultrasound in Medicine and Biology, 28:359-368, M a r c h 2002. [107] S. Sr inivasan, F . K a l l e l , R . Souchon, and J . Oph i r . A n a l y s i s of an A d a p t i v e S t r a in E s t i m a t i o n Technique i n Elas tography. Ultrasonic Imaging, 24:109-118, A p r i l 2002. [108] S. Sr inivasan, T . Krouskop , and J . O p h i r . C o m p a r i n g Elas tographic S t r a in Images w i t h M o d u l u s Images Obta ined U s i n g Nano-Indenta t ion: P r e l i m i n a r y Resul ts U s i n g Phan toms and Tissue Samples. Ultrasound in Medicine and Biology, 30:329-3434, 2004. [109] S. Sr inivasan and J . Oph i r . A zero-crossing s t ra in est imator i n elastography. Ultrasound in Medicine and Biology, 29:227-238, 2003. [110] S. Sr inivasan, J . O p h i r , and S . K . A l a m . Elas tographic Imaging U s i n g Staggered S t ra in E s t i -mates. Ultrasonic Imaging, 25:229-245, M a y 2002. BIBLIOGRAPHY 124 [111] S. Sr inivasan, R . R ighe t t i , and J . Oph i r . Tradeoffs between the ax i a l resolut ion and the signal-to-noise rat io i n elastography. Ultrasound in Medicine and Biology, 29:847-866, 2003. [112] C . S u m i . F i n e E la s t i c i ty Imaging U t i l i z i n g the Iterative R F - e c h o Phase M a t c h i n g M e t h o d . IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control, 46:158166, January 1999. [113] C . S u m i . Di rec t S t ra in Measurement Based on an Autocor re la t ion M e t h o d . In Proceedings of the Third International Conference on the Ultrasonic Measurement and Imaging of Tissue Elasticity, page 64, 2004. [114] C . S u m i . T h e r m a l Propet ies Recons t ruc t ion Technique Based on Tempera ture Measurement . In Proceedings of the Third International Conference on the Ultrasonic Measurement and Imaging of Tissue Elasticity, page 73, 2004. [115] G . E . Trahey, J . W . A l l i s o n , and O . T . V o n R a m m . A n g l e independent ul t rasonic detection of b lood flow. IEEE Transactions on Biomedical Engineering, B M E - 3 4 : 9 6 5 - 7 , December 1987. [116] M T r i s t a m , D C Barbosa , D O Cosgrove, D K Nass i r i , J C Bamber , and H i l l C R . Ul t rason ic s tudy of i n v ivo kinet ic characteristics of human tissues. Ultrasound in Medicine and Biology, 12:927-937, December 1986. [117] E . Turgay. Imaging V i s c o - E l a s t i c Proper t ies of Soft Tissue w i t h U l t r a sound . Mas te r ' s thesis, Un ive r s i ty of B r i t i s h C o l u m b i a , 2004. [118] T . Varghese, , and J . O p h i r . Performance op t imiza t ion i n elastography: Mul t i compress ion w i t h t empora l stretching. Ultrasonic Imaging, 18:193-214, 1996. [119] T Varghese, M Bi lgen , and J . Oph i r . Mu l t i r e so lu t i on imaging i n elastography. IEEE Trans-actions on Ultrasonics, Ferroelectrics and Frequency Control, 45:6575, 1998. [120] T . Varghese and O p h i r J . Charac te r iza t ion of elastographic noise using the envelope of echo signals. Ultrasound in Medicine and Biology, 24:543-555, 1998. [121] T . Varghese, E . E . Konofagou, J . O p h i r , S . K . A l a m , and M . B i lgen . D i rec t s t ra in es t imat ion i n elastography using spectral cross-correlation. Ultrasound in Medicine and Biology, 26:1525-1537, November 2000. [122] T . Varghese and J . Oph i r . E s t i m a t i n g tissue s t ra in from signal decorrelat ion using the cor-re la t ion coefficient. Ultrasound in Medicine and Biology, 22:1249-1254, 1996. [123] T . Varghese and J . Oph i r . A theoret ical framework for performance character izat ion of elas-tography: the s t ra in filter. IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control, 44:164-72, January 1997. [124] T . Varghese and J . Oph i r . Enhancement of echo-signal correlat ion i n elastography using t empora l stretching. IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control, 44:173-180, January 1997. [125] T . Varghese and J . Oph i r . T h e nonstat ionary s t ra in filter i n elastography, Pa r t I Frequency dependent at tenuation. Ultrasound in Medicine and Biology, 23:1343-56, 1997. BIBLIOGRAPHY 125 [126] T . Varghese and J . O p h i r . A n analysis of elastographic contrast-to-noise rat io performance. Ultrasound in Medicine and Biology, 24:915-924, 1998. [127] T . Varghese, J . O p h i r , E . Konofagou, F . K a l l e l , and R . R i g h e t t i . Tradeoffs In Elas tographic Imaging, Ul t rason ic Imaging. Ultrasonic Imaging, 23:216-248, October 2001. [128] T . Varghese, J . Zagzebski , and J . Jee. Elas tography imag ing of the rmal lesion i n the l iver i n v ivo fol lowing radio frequency ablat ion: P r e l i m i n a r y results. Ultrasound in Medicine and Biology, 28:1467-1473, 2002. [129] F . V i o l a and W . F . Walke r . Compar i son of t ime delay estimators i n medica l ul t rasound. In IEEE Ultrasonics Symposium, pages 1485-1488 vol .2. I E E E , October 2001. [130] F . Walker and E . Trahey. A fundamental l im i t on delay es t imat ion using par t ia l ly correlated speckle signals. IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control, 42:301-308, 1995. [131] T . E . Weber and R . L . Ho l l i s . A vis ion based correlator to act ively damp vibra t ions of a coarse-fine manipula tor . In IEEE International Conference on Robotics and Automation Proceedings, pages 818-825 vol .2 . I E E E , 1989. [132] P . Welch . T h e use of fast Fourier t ransform for the es t imat ion of power spectra: A method based on t ime averaging over short, modified periodograms. IEEE Transactions on Audio and Electroacoustics, 15:70-73, Junuary 1967. [133] Y . Yamakosh i , J . Sato, and T . Sato. Ul t rasonic imaging of in ternal v ib ra t i on of soft tissue under forced v ib ra t ion . IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control, 37(2):45-53, M a r c h 1990. [134] M . Yassine, F . Ossant, C . Imberdis, and F . Pa ta t . N e w Elas tographic A l g o r i t h m for i n v ivo S tudy of the Mecha n i c a l Behavior of Sk in . In Proceedings of the Third International Conference on the Ultrasonic Measurement and Imaging of Tissue Elasticity, page 59, 2004. [135] F Yeung , S F Lev inson , and K J . Parker . M u l t i l e v e l and mot ion model-based ul trasonic speckle t racking a lgor i thm. Ultrasound in Medicine and Biology, 24:427441, 1998. [136] K o Y o u n g - K i , K H y u n - G y u , L . Jong, K Y o u n g - R o , O h Hyeong-Cheol , and K o ; Sung-Jea. N e w mot ion es t imat ion a lgor i thm based on bit-plane ma tch ing and i ts V L S I implementa t ion . In IEEE Region 10 Conference, pages 848-851 vol .2 . I E E E , September 1999. Appendix A Real Time Motion Tracking Algorithms A.1 Phase Root Seeking (PRS) T h e P R S a lgor i thm is i l lus t ra ted i n F i g . A . 1 . It uses a modif ied N e w t o n i tera t ion me thod at each window to estimate the time-shift from the phase-shift of the signals: tn+1 =tn- T 7 7 T ( A - l J . _ <P(tn) Simi l a r ly to the T D E method, t ime shifts are est imated by the P R S method using a discrete number of windows at each R F line [90]. T h e t ime shift of the k - th window of two A- l ines centered around t}- = kAT is est imated by the following i terat ive formula 126 A.1 Phase Root Seeking (PRS) 127 FirstGuess TimeShiftOfPreviousWindow '^kfi = (A-2) EstimatedPhaseShift CurrentGuess PreviousGuess 1 'Tk,l = + arg Sloped fPhaseShift ftk+Tw/2 Jtk-Twj1 where a0(t) and b0(t) denote the corresponding pre-compression and post-compression base-band echo data . I is the i tera t ion index and L is the number of i terations that is mos t ly between 2 to 6. wo denotes the transducers nomina l center frequency, Tw denotes the w indow length and the base-band signals are calculated from analyt ic signals by [90]: ab (t) = a+ (t) e~ju>ot = [a (t) + jHilbert (a (t))} e -*" 0*, (A-3) F igure A . 1 : A Modi f i ed N e w t o n search far the es t imat ion of the t ime delay. T h e root of the phase is i terat ively searched. In each step of the i tera t ion, the phase is calculated and approximated by a linear funct ion w i t h a slope equal to the transducers centroid frequency. T h e intercept of this l inear function w i t h the abscissa is the new estimate for the t ime delay [90]. Fur thermore the post-compression signal has to be re-sampled by a sub-sample t ime shift [90]. A m o n g a l l available re-sampling methods, linear interpolation is used for the purpose of real t ime applicat ions since i t is more accurate for base-band da ta [89]. A.2 Combined Autocorrelation Method (CAM) 128 A . 2 Combined Autocorrelation Method (CAM) T h e me thod has the merits of phase domain processing but wi thout al iasing, since i t combines the result of phase correlat ion w i t h that of envelope correlat ion, b o t h of w h i c h are calculated d i rec t ly from the R F signal using complex autocorrela t ion processing [102,103]. Phase correlat ion is used for displacement measurement and envelope correlat ion is used for phase unwrapping . In order to perform the complex autocorrelat ion, phase informat ion is needed. S i m i l a r l y to the P R S method , the H i l b e r t t ransform can be used to convert the t ime-domain signals to analyt ic signals. Af te r this step the complex cross-correlation function is used to detect the phase shift: /•to/2 Rab (t- n)= a+(t + v) b*+ (t + nT/2 + v) dv = Ru (t; r - n T / 2 ) e-^o(r-r«r/2) ^ ( A_ 4) J-to/2 where T is per iod of ultrasonic wave, UJQ is carrier angular frequency, r is the t ime shift and Ru(t;T) is the autocorrelat ion of the envelope. Fo r the special case of n = 0, the displacement d can be obtained for the phase shift ip = LJQT. If the displacement is less than A / 4 i t can be obtained di rect ly and wi thout ambigui ty : < A - 5 > Since large displacements may be necessary for elast ici ty imaging, C A M uses another phase unwrapp ing step by using the envelope normal ized correlat ion coefficient defined by ^ u.^ \Rab(t;n)\ Cu{t'n)-\a(t)\\b(t + nT/2)\> (A"6) A c c o r d i n g to E q u a t i o n (A-4) and E q u a t i o n (A-6) , for each t ime t, two sets of Cu and (p = UQT may be obta ined as c u (t) — {Cu M , . . . , Cu 1 , C ° , C „ , . . . C^ , } v(t) = {<p-M, (A-7) (A-8) A.2 Combined Autocorrelation Method (CAM) 129 where C™ = C™ (t; n) and ipn is the phase of Rab (t; n). If M and N (size of the search region) are selected to be sufficiently large, one component of {</? (t)}, ipk, among the n components is not wrapped because it is obta ined from two sequences i n which the displacement at t is less than A / 4 . A t the same t ime C™ becomes m a x i m u m at n = k (the smaller the displacement the higher the correla t ion coefficient). Therefore the first step is to determine the value of n (among N+M+l values) wh ich maximizes the C™ , after wh ich the unwrapped phase shift may be obtained as cp(t) = arctan(Rab(t;n)) and the displacement can be calcula ted from E q u a t i o n (A-5) . In other words, the C A M a lgor i thm estimates the phase-shift at a number of locat ions w i t h a A / 2 spacing between them. T h e phase shift that corresponds to the m a x i m u m envelope correlat ion at that pos i t ion is then reported as the unwrapped phase shift and is used to estimate the displacement. Appendix B Correlation Based Mot ion Tracking Methods B . l Correlation Estimation Methods M o t i o n es t imat ion is performed on a pair of signals wh ich can be mathemat ica l ly expressed as E q u a t i o n ( B -2 ) a{t) = r 1 ( t ) + m ( i ) ( B - l ) b(t) = r2(i)+n2(i), where a (t) and b (t) are pre-compression and post-compression d ig i t ized R F data . In this expression, r\ (i) and r2 (i) are the echo signals received by the transducer, wh ich may have been decorrelated by physical processes, whereas n\ (i) and n2 (i) represent addi t ive noise in t roduced from electronic sources. For a smal l window inside each R F da ta the s ta t ionary property of the process is va l id . Unde r these condit ions, a correlat ion es t imat ion a lgor i thm can be used to estimate the time-shift inside smal l windows i n a (t) and b(t). C o m m o n pat tern match ing methods can be categorized i n below. • Corre lat ion: T h e correlat ion (Corr) between the signal a(t) and s ignal b(t) is defined as: 130 B.l Correlation Estimation Methods 131 t=i+W Corr (T) = a (*) x 6 (* + T) , (B-2) t=* where i is the first element of the current w i n d o w and W is the length of the window. T h i s method is wide ly used i n signal processing appl icat ions [35,129]. • Normalized Correlation: The normal ized correlat ion ( N C o r r ) between the signal a(t) and signal b(t) is defined as: NCr (r) - , a " " " " " 1 " ' (B-3) T h e normal ized Cor re la t ion scales the correla t ion coefficients between 0 and 1 and makes it possible to judge how accurate the pa t te rn match ing is accomplished, according to the coefncinet value (closer to 1 means more accurate) [127]. • Normalized Covariance: T h e normal ized covariance ( N C o v ) between the signal a(t) and signal b(t) is defined as : NCO, (T) . ar ( .M-» ) »(»('•••*)-1) , V E K + " ' <« ( 0 - S) 2 TS*W (»(« + r) - 5 ) 2 where a and b are the means of the current pre- and post-compression windows. T h e normal -ized covariance is most ly used i n image processing applicat ions [85] but i t has also been used for u l t rasound signals as wel l [37,129]. • Hybrid-Log Method: T h e hybr id- log method ( H L M ) between the signal a(t) and signal b{t) is defined as: t=i+W HLM(T)= [(a(t)-b(t + T))-ln(l + exp(2(a(t)-b(t + T))))], (B-5) t=i T h e Log compression transforms the mul t ip l ica t ive noise to addit ive [94]. • Polarity-Coincidence Correlation: B.l Correlation Estimation Methods 132 T h e polar i ty-coincidence Cor re la t ion ( P C C ) between the signal a(t) and s ignal b(t) was first ly used i n [17] for elastography and is denned as: t=i+W PCC(r)= ^2 sign (a (t)) x sign (b [t + T)) , (B-6) t=i • Hybrid-Sign Correlation: T h e hybr id-s ign correlat ion ( H S C ) between the signal a(t) and signal b(t) is defined as: t=i+W HSC (r) = a (*) x s i 9 n (b (* + T)). (B"7) t=i • Meyr-Spies Method: T h e meyr-spies method correlat ion between the signal a(t) and signal b(t) is defined as: t=i+W MSM ( r ) = ] T [o (t - 2) - a (*)] x b (t - 1 + r ) , (B-8) t=i • Sum of Absolute Differences: T h e s u m of absolute differences ( S A D ) between the signal a(t) and signal b(t) is defined as: t=i+W SAD(T)= Y \a{t)-b(t + T)\, (B-9) In contrast to previous methods, to find the pa t te rn w i t h S A D we search for the pos i t ion that minimizes the S A D while i n previous methods we search for the posi t ion tha t maximizes the correlat ion [19,35,129]. T h e advantage of S A D as a s imi la r i ty measure over normal ized cross correlat ion a lgor i thm is its s impl ic i ty i n implementa t ion . L i k e correlat ion, S A D a lgor i thm requires one operat ion per p ixe l whereas the normal ized cross correlat ion a lgor i thm requires five and normal ized covariance method requires eight operations per p ixe l . In [12] i t was shown that S A D is as accurate as correlat ion for bo th R F da ta and B - m o d e da t a i n detecting bo th ax ia l and lateral mot ion . • Sum of Squared Differences: T h e sum of squared differences (SSD) between the signal a(t) and signal b(t) is defined as: B.2 Subpixelling 133 t=i+W SSD(T)= {a{t)-b{t + T))2, (B-10) t=i T h e difference between S S D and S A D is that i t amplifies the large differences larger t han S A D . S imi la r to S A D we search for the posi t ion that minimizes the S S D [35,129]. • Normalized Sum of Squared Differences: T h e normal ized s u m of squared differences ( N S S D ) between the signal a(t) and signal b(t) is defined as: NSSD (r) = Sgr((<.m- S)-(M< + T ) - S ) f («(<) - a) 2 T . V (Ht + r)-^ • Logarithmic Compression Method: Logar i t hmic compression helps to reduce the decorrelation noise [17] and can be used prior to a l l previous methods to improve the performance of pat tern match ing . It compresses b o t h pre- and post-compression signals according to E q u a t i o n (B-12) . T h e compressed signals are then used for pat tern match ing . T h e logar i thmic compression of s ignal a(t) is defined as: a-Log (t) = sign (a (i)) x Logw (1 + \a (t)\), (B-12) A m o n g a l l available methods i t was shown that correlation and normalized correlation are among the most accurate methods for s t ra in es t imat ion [35,85,129] and because of th is proper ty they have been typ ica l ly used i n s train es t imat ion methods. B.2 Subpixelling Since t ime delays are generally not integral mult iples of the sampl ing per iod , the loca t ion of the largest sample of the cross correlat ion function is an inexact est imator of the loca t ion of the peak. Therefore, in terpola t ion techniques must be used between the samples of the cross correlat ion to improve the es t imat ion precision. A number of sub-sampling methods have been in t roduced i n literatures that use different curve f i t t ing methods to do the sub-sampling. Curve-fitting in terpola t ion can yie ld biased t ime-delay estimates. T h e ar t i factual effect of these bias errors B.2 Subpixelling 134 on elast ici ty imaging by elastography is discussed i n [16]. In addi t ion to curve-f i t t ing methods, there are a number of reconstructive interpolation techniques that can be used to find the exact pos i t ion of the peak i n the cross correlat ion function [16] but since these methods are computa t iona l ly very intensive, they are not considered here. T h e most popular sub-sampling methods that use curve f i t t ing methods can be categorized as follow: • Quadrat ic P a r a b o l a fitting: Three values of the correlat ion coefficients are considered: the highest value and its two neighbors. Quadra t i c in terpola t ion is appl ied to the three points , and the loca t ion that gives the peak value of the curve is chosen. T h e peak of the curve is determined by f inding the zero crossing of the first derivative [31,37]. A c c o r d i n g to F igure B . l , i f A is the distance of the loca t ion of the m a x i m u m coefficient to the peak point , then A is given i n [31] as: Coefficient value _b_ Coefficient index n+1 Figure B . l : E x a c t pos i t ion of the peak of cross correlat ion function is t yp ica l ly est imated by using in terpola t ion techniques. T h e image is provided by Turgay E . where b is the value of the m a x i m u m correlat ion coefficient, and a and c are the values of the left and right samples respectively. • Cosine fitting: A n o t h e r curve that has been used i n [16,27] to interpolate the peak i n the cross correlat ion function is the cosine function. T h i s model is reasonable for some signals w i t h rectangular and Gauss ian spectra and it has been used i n b o t h b lood flow es t imat ion B.3 Correlation Recalculation 135 and elastography. S imi l a r to previous method, the peak of a cosinusoid fi t ted to the three largest samples of the d ig i ta l cross correlat ion values (a, 6, c) was derived to be: A = - - f , (B-14) where UJQ is the angular frequency of the cosinusoid given by: wo = cos'1 {^^- ] , (B-15) and 0 is the phase of the cosinusoid given by: 6 = tan'1 ( — ^ — ] , (B-16) These two curve-f i t t ing methods were evaluated and compared i n [16]. T h e y showed that cosine in terpola t ion performs better than parabol ic in terpola t ion i n an a rb i t ra r i ly smal l t ime interval near the peak and the goodness of the cosine-fit deteriorates w i t h increased length of this interval . B .3 Correlation Recalculation After f inding the exact pos i t ion of the peak w i t h one of the subpixel l ing methods the value of the correlat ion at that point also needs to be recalculated. A l t h o u g h the value of the corre la t ion has no effect on the result of mo t ion t racking but i t can be used later to validate the es t imated mot ion . In this thesis three points Lagrangian in terpola t ion method is used to recalculate the value of the correlat ion. S imi l a r to subpixel l ing, the same three values of the correlat ion coefficients are considered: the highest value and its two neighbors (F ig . B . l ) . Lagrang ian in terpola t ion is appl ied to the three points, and the value of the correlat ion at the peak pos i t ion is calculated according to E q u a t i o n (B-17) . Correlation = a x A x (A - l ) / 2 - b x (A + 1) x (A - 1) + c x (A - 1) x A / 2 , (B-17) B.3 Correlation Recalculation 136 T h e accuracy of the recalculat ion can be improved by increasing the number of points and using higher order in terpola t ion techniques wh ich requires es t imat ion of the correla t ion at more points . Appendix C Computational Issues C l Computational Cost of the LSE Method If the ca lcula t ion is r u n wi thout any s impl i f icat ion and according to the E q u a t i o n (4.4) from left to r ight , the complexi ty of the a lgor i thm for a kernel size of TV w i l l be: 1. AT: N x 2 M a t r i x Transposing, (cost is ignored). 2. ATA: N x 2 by N x 2 M a t r i x M u l t i p l i c a t i o n , [4iV,4(JV - 1)]. 3. [ . A ^ - A ] - 1 : 2 x 2 M a t r i x Inversion, (cost is ignored). 4. [ATA]-1AT: 2 x 2 by 2 x N M a t r i x M u l t i p l i c a t i o n , [4N,2N]. 5. [ATA]-1ATd: 2 x N by N x 2 M a t r i x M u l t i p l i c a t i o n , [4N,4(N - 1)]. wh ich is 12N mul t ip l i ca t ion and l O i V summations for each window ([12iV, 10A 7]) for the or ig inal implementa t ion of the L S E . T h e implementa t ion that is used i n this thesis (Section 7.1.2.1) reduces the complex i ty of L S E to 2 mul t ip l ica t ions and 8 summations (Equation(7.6)) for each window where N is the kernel size. Therefore the overall cost of the second fast version w i l l be [2, 8] for each window which is independent from N. Cons ider ing the number of mul t ip l ica t ions and summat ion, the current implementa t ion of the L S E performs almost 6N t imes faster t han s tandard implementat ion. T h i s means for the kernel 137 Cl Computational Cost of the LSE Method 138 size of 10 the second approach is almost 60 times faster than the or ig ina l implementa t ion . U s i n g this approach the ca lcula t ion of the least squares es t imat ion becomes possible for even very large kernel size (N > 50) i n the real t ime appl icat ions whi le the size of the kernel can also be adjusted i n real t ime. 

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0065559/manifest

Comment

Related Items